Sf: Add st_crop()

[adrfantini]

Currently I crop sf objects this way:

example(st_read)
st_intersection(nc, st_set_crs(st_as_sf(as(raster::extent(-82, -80, 35, 36), "SpatialPolygons")), st_crs(nc)))

or, more briefly, using spex:

st_intersection(nc, st_set_crs(st_as_sf(spex::spex(c(-82, -80, 35, 36))), st_crs(nc)))

or, most briefly, using tmaptools:

tmaptools::crop_shape(nc, raster::extent(-82, -80, 35, 36))
#or
tmaptools::crop_shape(nc, matrix(c(-82, -80, 35, 36), nrow=2))

I think anyone using spatial data and sf comes around doing this sooner or later, would this warrant an st_crop() function, even though there are indeed ways to do this?

[edzer]

sf has a native CRS-aware bbox constructor:

> st_bbox(c(xmin=1, ymin=1, xmax=2, ymax=2), crs = 4326)
xmin ymin xmax ymax 
   1    1    2    2 

What is the benefit of having an st_crop if st_intersection already does the cropping?

[adrfantini]

So, just to complete the above, this also works:

st_intersection(nc, st_as_sfc(st_bbox(c(xmin=-82, ymin=35, xmax=-80, ymax=36), crs=st_crs(nc))))

The proposed approach, e.g.:

st_crop(nc, c(-82, -80, 35, 36))

Is more intuitive, cleaner, easier to understand and to remember. Also, it is much more easily discoverable looking through the manual and the functions provided by sf. Up to you, of course, feel free to close.

[edzer]

I see your point. Would you agree that requiring a named vector, as in

st_crop(nc, c(xmin = -82, xmax,= -80, ymin = 35, ymax = 36))

is better, in order to avoid the mess of bbox parameters order (across packages)?

[adrfantini]

Personally I would prefer simplicity (it's a lot of typing!) and assume the same order as raster::extent:

> extent(c(-82, -80, 35, 36))
#or
> extent(-82, -80, 35, 36)
class       : Extent 
xmin        : -82 
xmax        : -80 
ymin        : 35 
ymax        : 36

but I understand you might want to keep in line with the syntax of st_bbox.

[edzer]

Requiring names does no longer require any order.

[adrfantini]

I totally see why that would be more convenient given that different packages assume different orders in creating boxes and extents, I'm fine either way.

[adrfantini]

Thanks!

[adrfantini]

Getting:

> st_crop(nc, c(xmin = -82, xmax= -80, ymin = 35, ymax = 36))
although coordinates are longitude/latitude, st_intersection assumes that they are planar
Error: nrow(x) == length(value) is not TRUE

[lbusett]

Hi,

on this, I think it could be useful also to pass a generic "spatial" object as the second argument, from which the "cropping bbox" would be automatically computed. Automatic reprojection of the cropping bbox to the CRS of the first argument could also be useful, IMO. For example, I have a very similar functionality implemented here:

https://github.com/lbusett/sprawl/blob/master/R/crop_vect.R

What do you think? I could give it a go and do a PR, if you agree.

[edzer]

@lbusett the current implementation tries to make a bbox from the second argument, if it not already is one. Extending this would then work by adding st_bbox methods; it now e.g. already works with a raster::Extent object as second argument; we have

> methods(st_bbox)
 [1] st_bbox.CIRCULARSTRING*     st_bbox.COMPOUNDCURVE*     
 [3] st_bbox.CURVEPOLYGON*       st_bbox.Extent*            
 [5] st_bbox.GEOMETRYCOLLECTION* st_bbox.LINESTRING*        
 [7] st_bbox.MULTICURVE*         st_bbox.MULTILINESTRING*   
 [9] st_bbox.MULTIPOINT*         st_bbox.MULTIPOLYGON*      
[11] st_bbox.MULTISURFACE*       st_bbox.numeric*           
[13] st_bbox.POINT*              st_bbox.POLYGON*           
[15] st_bbox.POLYHEDRALSURFACE*  st_bbox.Raster*            
[17] st_bbox.sf*                 st_bbox.sfc*               
[19] st_bbox.Spatial*            st_bbox.TIN*               
[21] st_bbox.TRIANGLE*          

I'm reluctant to adopt automatic reprojection; if here, there are many other places where we'd have to do this.

[edzer]

... and yes, these warnings are annoying, but setting agr here to a value that suppresses the warning does not mean that your cropped polygon will have a correct value for spatially extensive variables, such as population count.

[lbusett]

@edzer

Oh, right: didn't notice the st_bbox call. If it already works with all those "spatial objects types" it's already good as is! Thanks!

Concerning automatic reprojection: I understand. Maybe it could be worth Aborting in case of mismatching CRSs ? Otherwise, I think that due to:

y = st_bbox(y, crs = st_crs(x))

passing a spatial object as the second argument would risk to give wrong results (though that would not work in the case of passing a "named array", since it is not needed to provide a CRS and it is assumed that the CRS is that of the object to be cropped)). Though maybe the crs = st_crs(x) argument is ignored in case "y" already has a CRS - in that case st_intersection would already take care of aborting and this can be ignored. Sorry - was not able to test.

Concerning the warnings: you are absolutely right. That lines of code should be removed. Did not even remember they were there.

[edzer]

Thanks; this should only inherit it from x when not already set by the st_bbox call.