3

我想知道它是否比下面使用 rstats sf 的示例更简单，以计算底层网格中多边形 shapefile 的（分数）面积覆盖率。到目前为止，我还没有找到可以为我完成这项工作的功能。这是我的例子：

library(sf)

## create a multipolygon
p1 <- rbind(c(7, 48), c(8, 52), c(11, 49), c(9, 48), c(7, 48))
p2 <- rbind(c(8, 50), c(8, 51), c(9, 50), c(8, 50))
mpol <- st_sfc(st_multipolygon(list(list(p1, p2))), crs=st_crs(4326))

## create a half degree grid covering the multipolygon
grid <- st_make_grid(mpol, cellsize=0.5)

## Create a sf object with the required attributes
## package 'lwgeom' reqiured
area <- st_area(grid)
grid <- st_as_sf(data.frame(ID=c(1:length(area)), area=area, geometry=grid))

## check 'map' elements
plot(grid['area'])
plot(mpol, col="#FF0000AA", add=TRUE)

现在开始计算，我首先通过从边界框（'bpol'）中减去'mpol'来反转多边形（'ipol'）并计算新的覆盖区域。

## create a rectangular bounding box around the polygon
bbox <- st_bbox(mpol)
bpol <- st_sfc(st_polygon(list(rbind(c(bbox['xmin'], bbox['ymin']),
                                     c(bbox['xmax'], bbox['ymin']),
                                     c(bbox['xmax'], bbox['ymax']),
                                     c(bbox['xmin'], bbox['ymax']),
                                     c(bbox['xmin'], bbox['ymin'])))), crs=st_crs(4326))

## invert the multipolygon, by substracting it from the bounding box
ipol <- st_difference(bpol, mpol)

## substract the inverted polygon from the grid
gridinpoly <- st_difference(grid, ipol)

## calculate new 'cropped' area
gridinpoly$croppedArea = st_area(gridinpoly)
plot(gridinpoly['area'])
plot(gridinpoly['croppedArea'])

## Finally the fractional coverage is calculated:
gridinpoly$frac = gridinpoly$croppedArea / gridinpoly$area
plot(gridinpoly['frac'])

如果一个网格单元被覆盖两次，然后从网格中断处裁剪倒多边形，它会变得更加复杂

p3 <- rbind(c(9.75, 50.5), c(10, 51), c(10.5, 50), c(9.75, 50.5))
mpol <- st_sfc(st_multipolygon(list(list(p1, p2, p3))), crs=st_crs(4326))
ipol <- st_difference(bpol, mpol)
gridinpoly <- st_difference(grid, ipol)

错误消息：“CPL_geos_op2(op, st_geometry(x), st_geometry(y)) 中的错误：评估错误：TopologyException：输入几何 1 无效：在 9.75 50.5 处或附近的嵌套壳在 9.75 50.5。

现在我的问题是，是否已经有一个功能可以完成这项工作？我错过了文档中的某些内容吗？当网格单元被裁剪多次时，如何解决错误？

谢谢和亲切的问候约尔格

4

1 回答 1

3

对我来说似乎更自然/更容易使用st_intersection。

复制和示例几乎与您的相同。唯一的区别是没有指定 crs。它应该适用于任何 crs 但不适用于经度纬度（见下文）

library(sf)
## create a multipolygon
p1 <- rbind(c(7, 48), c(8, 52), c(11, 49), c(9, 48), c(7, 48))
p2 <- rbind(c(8, 50), c(8, 51), c(9, 50), c(8, 50))
mpol <- st_sfc(st_multipolygon(list(list(p1, p2))))

## create a half unit grid covering the multipolygon
grid <- st_make_grid(mpol, cellsize=0.5)
area <- st_area(grid)
grid <- st_as_sf(data.frame(ID=c(1:length(area)), area=area, geometry=grid))

计算多边形内的网格表面的分数：

tmp <- st_intersection(grid, mpol)
tmp$area <- st_area(tmp)
tmp$frac <- tmp$area/unique(area)
plot(tmp['frac'])

# If needed, you can remove the non polygons :
tmp[tmp$area > 0,]

由于两个原因，这不适用于经度纬度：

  1. 当st_area用于经度纬度时，sf调用另一个geosphere适用于sp对象的包。然后，您的sf对象会即时转换为sp对象。
  2. 在您的情况下，网格和多边形相交的结果不是多多边形，而是由点、线和多边形（类型 = 几何）组成的复杂对象。并且从这种类型的sf对象转换为sp不起作用

但是如果您只想要表面的百分比，则可能不建议并且不需要测量经度/纬度上的区域。

至于您使用 p3 的错误消息，我不确定您要模拟什么。如果您想要一个由 p1 和 p3 组成的多面体，在 p1 中有孔 (p2)，那么要创建的语法mpol必须稍有不同（并且区域计算工作正常）：

p3 <- rbind(c(9.75, 50.5), c(10, 51), c(10.5, 50), c(9.75, 50.5))
mpol <- st_sfc(st_multipolygon(list(list(p1, p2), list(p3))))
plot(mpol)

tmp <- st_intersection(grid, mpol)
tmp$area <- st_area(tmp)
tmp$frac <- tmp$area/unique(area)
plot(tmp['frac'])
于 2018-01-25T09:51:13.813 回答