28

Is there a function to crop sf map object, similar to maptools::pruneMap(lines, xlim= c(4, 10), ylim= c(10, 15)) used for SpatialPolygon or SpatialLine?

I am considering st_intersection() but there can be proper way.

4 Answers 4

35

Since today, there is a st_crop function in the github version of sf (devtools::install_github("r-spatial/sf"), probably on CRAN in the near future too).

Just issue:

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

The vector must be named with xmin xmax ymin ymax (in whichever order).

You can also use any object that can be read by st_bbox as cropping limits, which is very handy.

0
22

st_intersection is probably the best way. Find whatever way works best to get an sf object to intersect with your input. Here's a way using the convenience of raster::extent and a mix of old and new. nc is created by example(st_read):

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

I don't think you can coax st_intersection to not need an exact matching CRS, so either set both of them to NA or make sure they are the same. There's no easy tools for bbox/extent afaik, so using raster is good way to make things easy-ish.

3
  • Thanks so much @mdsumner, it worked like a charm. I spent hours on st_intersection but couldn't solve it myself.
    – Kazuhito
    Mar 6, 2017 at 11:48
  • You can now use spex::spex to replace the st_as_sf(as(...)) call. Also, tmaptools::crop_shape() can do this.
    – AF7
    Mar 5, 2018 at 10:36
  • 1
    sf now includes st_crop, see my answer for details.
    – AF7
    Apr 23, 2018 at 7:34
5

Another workaround, for me it was faster for larger shapefiles:

library(sf)
library(raster)
library(rgeos)
library(ggplot2)

# Load National Forest shapefile
# https://data.fs.usda.gov/geodata/edw/edw_resources/shp/S_USA.AdministrativeForest.zip
nf.poly <- st_read("S_USA.AdministrativeForest"), "S_USA.AdministrativeForest")

crop_custom <- function(poly.sf) {
  poly.sp <- as(poly.sf, "Spatial")
  poly.sp.crop <- crop(poly.sp, extent(c(-82, -80, 35, 36)))
  st_as_sf(poly.sp.crop)
}

cropped <- crop_custom(nf.poly)
5
  • Thanks. It is interesting workflow, combination of raster::crop() and st_as_sf()...+1 from me. I wish we can have this kind of readily accessible function like crop() in future versions of sf. As to the speed, a quick run of system.time with your function reported user: 5.42, system: 0.09, elapsed 5.52, while st_intersection() approach was user: 1.18, system: 0.05, elapsed 1.23 on your dataset. (Probably my environment is different with yours... not sure.)
    – Kazuhito
    Apr 1, 2017 at 7:33
  • That's interesting - the st_intersection approach takes around 80s for me.
    – pbaylis
    Apr 3, 2017 at 3:36
  • Keep in mind that the raster::crop function, when applied to sp geometry objects, acts a wrapper for rgeos functions. Albeit, a very convenient wrapper. The GEOS API operates on WKT objects so, will invariably will be a standard for sf overlay operations. Apr 14, 2017 at 15:46
  • 1
    And it changes over time, sf now has "spatial indexing" built-in in 0.5-1 cran.r-project.org/web/packages/sf/news.html so is probably now faster than sp/rgeos.
    – mdsumner
    Jul 3, 2017 at 10:16
  • 1
    sf now includes st_crop, see my answer for details.
    – AF7
    Apr 23, 2018 at 7:34
1

@mdsumner's solution as a function. Works if rasta is a RasterBrick, extent, bbox, etc.

# Crop a Simple Features Data Frame to the extent of a raster
crop.sf = function(sfdf, rasta) {
  st_intersection(sfdf, st_set_crs(st_as_sf(as(extent(rasta), "SpatialPolygons")), st_crs(sfdf)))
}

It throws away the crs information of the raster because I don't know how to convert a raster crs() into an st_crs()

On my machine and for my data sample this has equivalent performance to raster::crop with a SpatialLinesDataFrame version of the data.

@pbaylis's solution is about 2.5 times slower:

# Slower option.
crop.sf2 = function(sfdf, rasta) {
  st_as_sf(crop(as(sfdf, "Spatial"), rasta))
}

Edit: Somebodies comment suggests spex, which produces SpatialPolygons with the crs from the rasta, if it has a crs.

This code uses the same method as spex:

# Crop a Simple Features Data Frame to the extent of a raster
crop.sf3 <- function(sfdf, rasta) {
  # Get extent and crs
  ext.sp <- as(extent(rasta), "SpatialPolygons")
  crs(ext.sp) <- crs(rasta)

  # crop
  st_intersection(sfdf, st_as_sf(ext.sp))
}
1
  • sf now has an st_crop function which is probably worth checking out.
    – cmc
    Jul 23, 2019 at 16:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.