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.

Kazuhito
  • 30,746
  • 5
  • 69
  • 149

4 Answers4

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.

AF7
  • 703
  • 1
  • 8
  • 13
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.

mdsumner
  • 8,196
  • 1
  • 27
  • 30
  • Thanks so much @mdsumner, it worked like a charm. I spent hours on st_intersection but couldn't solve it myself. – Kazuhito Mar 06 '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 05 '18 at 10:36
  • 1
    sf now includes st_crop, see my answer for details. – AF7 Apr 23 '18 at 07: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)
pbaylis
  • 210
  • 1
  • 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 01 '17 at 07:33
  • That's interesting - the st_intersection approach takes around 80s for me. – pbaylis Apr 03 '17 at 03: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 https://cran.r-project.org/web/packages/sf/news.html so is probably now faster than sp/rgeos. – mdsumner Jul 03 '17 at 10:16
  • 1
    sf now includes st_crop, see my answer for details. – AF7 Apr 23 '18 at 07: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))
}
cmc
  • 131
  • 4