# Cropping sf object in R? [duplicate]

This question already has an answer here:

Newish to R Very new to GIS

I've grabbed some vector spatial census data using the cancensus package. I now have a spatial data frame with variables and geometry I plot it and decide to focus on a part of the whole province

My thinking here is to create a rectangle with the coordinates that I am interested in using extent(xmin, xmax, ymin, ymax)

Assign it the same CRS as the spatial data frame and then use st_intersection. Is this the right approach?

``````extent(box) <- c(-81, -74, 41.68, 45)
st_set_crs(box, "+proj=longlat +datum=WGS84 +no_defs") # Fails
ont_crop <- st_intersection(census_data, box)
``````

I'm not sure how to show you a reproducible example since census data is large. Please let me know what I am doing wrong.

## marked as duplicate by dbaston, Evil Genius, xunilk, Andre Silva, VinceJun 19 '18 at 2:43

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

• Welcome to GIS.se! Your approach broadly sounds OK - so what happened when you tried it? What went wrong? Edit and post the error message or explain why it fails! – Simbamangu May 11 '18 at 5:30
• I believe that in the development version you can just use sf::crop – Jeffrey Evans May 11 '18 at 7:03

## 3 Answers

One possibility is to define a polygon with sf, set the crs and intersect the features

``````library(sf)
# 5 points to close the geometry
mat <- list(matrix(c(-81, 41.68,
-74, 41.68,
-74, 45,
-81, 45,
-81, 41.68),
ncol = 2,
byrow = TRUE))
box <- st_polygon(mat)
box <- st_geometry(box)
box <- st_set_crs(box, "+proj=longlat +datum=WGS84 +no_defs")
st_intersection(census_data, box)
``````
• I am using the CRAN version of sf so st_crop is not an option. When I try this solution I get the warning, "although coordinates are longitude/latitude, st_intersection assumes that they are planar." So I used st_transform to move both into a projected CRS. Still, I was left with the warning, "attribute variables are assumed to be spatially constant throughout all geometries". I'm not sure what that means. – ixodid May 12 '18 at 1:51
• See r-spatial.github.io/sf/articles/… for an explanation of that warning, @ixodid. – obrl_soil May 12 '18 at 23:17

I was able to get some census data from the cancensus package and demonstrate how you might crop it here:

get data as sf object:

``````census <- get_census(dataset='CA16', regions=list(CMA="59933"), vectors=c("v_CA16_408","v_CA16_409","v_CA16_410"), level='CSD', geo_format = "sf")
``````

Create a box first by converting the extent object to SpatialPolygons then to sf object and set crs to 4326:

``````extent_ras <- extent(-122.9, -122.5, 49.1, 49.3)
extent_sf <- st_set_crs(st_as_sf(as(extent_ras, "SpatialPolygons")), 4326)
``````

Finally use st_intersection:

``````cropped <- st_intersection(census, extent_sf)
``````

Here is the cropped census data: Keep it simple using `st_crop` (credits to @JeffreyEvans's comment above):

``````library(cancensus)
library(sf)

## example from ?get_census
options(cancensus.api_key = "XXX") # insert own api token
census_data <- get_census(dataset='CA16', regions=list(CMA="59933"),
vectors=c("v_CA16_408","v_CA16_409","v_CA16_410"),
level='CSD', geo_format = "sf", labels="short")

ont_crop <- st_crop(census_data
, xmin = -122.9, xmax = -122.5, ymin = 49.1, ymax = 49.3)
`````` 