20

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.

17

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.

  • Thanks so much @mdsumner, it worked like a charm. I spent hours on st_intersection but couldn't solve it myself. – Kazuhito Mar 6 '17 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 '18 at 10:36
  • 1
    sf now includes st_crop, see my answer for details. – AF7 Apr 23 '18 at 7:34
22

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.

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)
  • 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 '17 at 7:33
  • That's interesting - the st_intersection approach takes around 80s for me. – pbaylis Apr 3 '17 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. – Jeffrey Evans Apr 14 '17 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 '17 at 10:16
  • 1
    sf now includes st_crop, see my answer for details. – AF7 Apr 23 '18 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))
}
  • sf now has an st_crop function which is probably worth checking out. – cmc Jul 23 at 16:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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