How to crop a polygon from an underlying grid with the R sf package?

Question

I wonder if it is simpler than in the below example with rstats sf, to calculate the (fractional) areal coverage of a polygon shapefile in an underlying grid. So far I haven't found a function, which does that job for me. Here is my example:

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)

Now comes the calculation, where I first invert the polygon ('ipol') by substracting 'mpol' from the bounding box ('bpol') and calculate the new covered area.

## 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'])

And it can get even more complicated if a gridcell is covered twice, then cropping the inverted polygon from the grid breaks

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)

ERROR message: "Error in CPL_geos_op2(op, st_geometry(x), st_geometry(y)) : Evaluation error: TopologyException: Input geom 1 is invalid: Nested shells at or near point 9.75 50.5 at 9.75 50.5."

Now my question is, is there already a function, which does that job? Did I miss something in the documentation? And how to solve the error, when a gridcell is cropped several times?

Thanks & kind regards Jörg

Answer 1

It seems more natural/easy to me to use st_intersection.

Reproduce and example almost identical to yours. The only difference is that the crs is not specified. It should work with any crs but not with longitude latitude (see below)

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))

Compute the fraction of the surface of the grid that is inside the polygon :

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,]

This would not work with longitude latitude for a combination of 2 reasons :

1. When st_area is used on longitude latitude, sf calls another package geosphere that works on sp objects. Your sf object is then transformed on the fly into a sp object.
2. In your case the result of the intersection of your grid and polygons is not a multipolygon but a complex object composed of points, lines and polygons (type = Geometry). And the conversion from this type of sf object to sp does not work

But measuring areas on longitude/latitude is probably not recommended and not needed if you only want a % of surface.

As for your error message with p3, I'm not certain about what you want to simulate. If you want a multipolygon composed of p1 and p3 with hole (p2) in p1, then the syntax to create mpol must be slightly different (and the areas computations work correctly) :

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'])