3

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, Vince Jun 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.

  • 1
    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
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)
  • 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
3

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: enter image description here

3

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)

st_crop

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