根据另一个 sf 对象计算 sf 对象的面积 [英] Calculate area of sf object based on another sf object

问题描述

我有一个研究区(左)，还有一条道路(中间).我想计算由两种模式(右)相交产生的面积——这将产生 5 个子区域，这些子区域的总和为研究区域对象的总面积.

I have a study area (left), and roads dissecting it (middle). I would like to calculate the areas that result from the intersection of both patterns (right) -- this would result in 5 subareas which sum to the total area of the study area object.

两个对象都是 sf 并且具有相同的投影(见下文).

Both objects are sf and have the same projection (see below).

如何使用 sf 库计算子区域?

How can I calculate the subareas using the sf library?

> str(myarea)
Classes ‘sf’ and 'data.frame':  1 obs. of  4 variables:
 $ Id      : num 0
 $ Habitat : num 1
 $ Area    : num 202
 $ geometry:sfc_POLYGON of length 1; first list element: List of 1
  ..$ : num [1:135, 1:2] 858895 859084 859358 859865 860105 ...
  ..- attr(*, "class")= chr  "XY" "POLYGON" "sfg"
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA
  ..- attr(*, "names")= chr  "Id" "Habitat" "Area"


> str(roads)
Classes ‘sf’ and 'data.frame':  12 obs. of  23 variables:
 $ LINKNO    : Factor w/ 13673 levels "A1-1-1-1","A1-1-1-2",..: 7483 7327 7326 7325 2433 2436 174 2438 2439 2434 ...
 $ ROADNO    : Factor w/ 2981 levels "A1","A104","A104A",..: 1444 1415 1415 1415 246 246 2 247 247 246 ...
 $ START_KM  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ END_KM    : num  14 8 4 1 5 7 16 10 7 6 ...
 $ LENGTHKM  : num  14.07 7.95 4.25 1.32 4.92 ...
 $ STRATDESC : Factor w/ 2908 levels "1005","1699A",..: 61 1778 1112 1112 628 1115 92 1112 1112 628 ...
 $ ENDDESC   : Factor w/ 2990 levels "1","A.P.LINE",..: 58 1667 1668 1668 1025 72 655 56 56 1025 ...
 $ ROADCLASS : Factor w/ 12 levels "A","B","C","D",..: 5 5 5 5 4 4 1 4 4 4 ...
 $ CLASS     : Factor w/ 3 levels "Primary (Trunk)",..: 3 3 3 3 2 2 1 2 2 2 ...
 $ REGION    : Factor w/ 9 levels "0","Central",..: 8 8 8 8 8 8 8 8 8 8 ...
 $ WIDTH     : num  6 5.5 5.5 5.5 5.5 5.5 7 5.5 5.5 6 ...
 $ LANES     : int  2 2 2 2 2 2 2 2 2 2 ...
 $ TYPE      : Factor w/ 10 levels "230","Concrete (Jt-Plain)",..: 7 4 4 4 4 6 7 6 4 10 ...
 $ SURFTYPE  : Factor w/ 6 levels "Gravel","Paved",..: 2 6 6 6 6 6 2 6 6 2 ...
 $ CONDITION : Factor w/ 6 levels "Excellent","Fair",..: 4 6 2 4 2 3 4 4 4 3 ...
 $ AADT      : num  0 0 69 69 50 ...
 $ iso3      : Factor w/ 1 level "KEN": 1 1 1 1 1 1 1 1 1 1 ...
 $ AICD_REG  : int  NA NA NA NA NA NA 1 NA NA NA ...
 $ TAH       : int  NA NA NA NA NA NA 1 NA NA NA ...
 $ WB_Project: int  NA NA NA NA NA NA NA NA NA NA ...
 $ Kamp_Momba: int  NA NA NA NA NA NA 1 NA NA NA ...
 $ Nairo_Addi: int  NA NA NA NA NA NA NA NA NA NA ...
 $ geometry  :sfc_GEOMETRY of length 12; first list element:  'XY' num [1:581, 1:2] 858210 858211 858215 858219 858223 ...
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
  ..- attr(*, "names")= chr  "LINKNO" "ROADNO" "START_KM" "END_KM" ...

推荐答案

您可以使用 lwgeom 库中的 st_split()，然后使用 st_area 计算面积.我只是编了一些数据来演示.请注意，这会将您的初始多边形剪裁成多个多边形.

You can use the st_split() from the lwgeom library and then calculate area with st_area. I just made up some data to demonstrate. Note that this will clip your initial polygon into multiple polygons.

library(sf)
library(lwgeom)

# Some sample data
poly <- st_bbox(c(xmin = 0, xmax = 15, ymin = 0, ymax=15)) %>% st_as_sfc()
lines <- data.frame(x = seq(-1, 18), y = seq(1,20)) %>%      st_as_sf(coords=c("x","y")) %>% st_union() %>% st_cast('LINESTRING')

初始

# split and add area to ssf
spltply <- lwgeom::st_split(poly, lines) %>% st_collection_extract(c("POLYGON")) %>% st_sf()
spltply['area'] <- st_area(spltply) 

拆分后

这篇关于根据另一个 sf 对象计算 sf 对象的面积的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！