Cropping sf object in R?

Question

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.

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 05:30
I believe that in the development version you can just use sf::crop — Jeffrey Evans, May 11 '18 at 07:03

score 5 · Answer 1 · answered May 11 '18 at 09:51

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)

score 3 · Answer 2 · answered May 11 '18 at 06:30

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:

score 2 · Accepted Answer · answered May 11 '18 at 07:25

2

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)

answered May 11 '18 at 07:25

nebi

611
5
16

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 01:51
See https://r-spatial.github.io/sf/articles/sf1.html#how-attributes-relate-to-geometries for an explanation of that warning, @ixodid. – obrl_soil May 12 '18 at 23:17

Cropping sf object in R?

3 Answers3