r - 如何使用 R sf 包从底层网格裁剪多边形？

Question

我想知道它是否比下面使用 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:输入 geom 1 无效:9.75 50.5 处或附近的嵌套壳。”

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

谢谢和亲切的问候 约尔格

Answer 1

对我来说，使用 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,]

由于 2 个原因的组合，这不适用于经度纬度:

1. 当 st_area 用于经度纬度时，sf 调用另一个适用于 sp 对象的包 geosphere。然后，您的 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'])