title: How to crop a polygon from an underlying grid with the R sf package?
body: <p>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:</p>
  
  <pre><code>library(sf)
  
  ## create a multipolygon
  p1 &lt;- rbind(c(7, 48), c(8, 52), c(11, 49), c(9, 48), c(7, 48))
  p2 &lt;- rbind(c(8, 50), c(8, 51), c(9, 50), c(8, 50))
  mpol &lt;- st_sfc(st_multipolygon(list(list(p1, p2))), crs=st_crs(4326))
  
  ## create a half degree grid covering the multipolygon
  grid &lt;- st_make_grid(mpol, cellsize=0.5)
  
  ## Create a sf object with the required attributes
  ## package 'lwgeom' reqiured
  area &lt;- st_area(grid)
  grid &lt;- 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)
  </code></pre>
  
  <p>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.</p>
  
  <pre><code>## create a rectangular bounding box around the polygon
  bbox &lt;- st_bbox(mpol)
  bpol &lt;- 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 &lt;- st_difference(bpol, mpol)
  
  ## substract the inverted polygon from the grid
  gridinpoly &lt;- 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'])
  </code></pre>
  
  <p>And it can get even more complicated if a gridcell is covered twice, then cropping the inverted polygon from the grid breaks</p>
  
  <pre><code>p3 &lt;- rbind(c(9.75, 50.5), c(10, 51), c(10.5, 50), c(9.75, 50.5))
  mpol &lt;- st_sfc(st_multipolygon(list(list(p1, p2, p3))), crs=st_crs(4326))
  ipol &lt;- st_difference(bpol, mpol)
  gridinpoly &lt;- st_difference(grid, ipol)
  </code></pre>
  
  <p>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."</p>
  
  <p>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?</p>
  
  <p>Thanks &amp; kind regards
  Jörg</p>