3

I want to clip a polygon shapefile using sf::st_intersection, however the territory is not cut straight.

Here is a RE. I use a Canadian province shapefile that can be downloaded easily:

library(sf)
library(dplyr)
    
the_zip <- tempfile(fileext=".zip")
download.file("http://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/files-fichiers/2016/lpr_000b16a_e.zip",the_zip)
unzip(the_zip, exdir = dirname(the_zip))

can <- st_read(paste0(dirname(the_zip), "/lpr_000b16a_e.shp")) %>% 
  filter(PRNAME %in% c("Quebec / Québec", "Ontario")) %>% 
  st_simplify %>% 
  st_transform(st_crs(4326)) 

clip <- st_as_sfc(st_bbox(c(xmin = -100, xmax = -60, ymax = 50, ymin = 40), crs = st_crs(4326)))

To help visualize the problem here is the shape I want to cut and the clipping polygon in blue:

dev.off()
plot(can[,1], reset=F, axes=T)
plot(clip, add=T, border="blue", lwd=2)

enter image description here

Now, I can clip the polygon:

can_clipped <- st_intersection(can, clip)

However, when I plot it, I get:

dev.off()
plot(can_clipped[,1], reset=F, axes=T)
plot(clip, add=T, border="blue", lwd=2)

enter image description here

As you can see, the clip do not fit well with the desired output. Also, both polygons in the shapefile seems bent.

What causes this problem and can it be fixed?

Could it be linked with the fact I'm clipping on shapes in EPSG:4326?

I do want to cut my shape on a strait line at the 50th parallel, so I don't really want to play around another projection.

Or maybe it's related to the warning I get when doing the st_intersection?

Warning message:
attribute variables are assumed to be spatially constant throughout all geometries

1 Answer 1

9

This is a consequence of sf now working with spherical geometry by default. The real straight line between the two top corner points of your rectangle are great circles, which will be bulging north.

Switch that off and sf will assume lat-long is planar and you get the result you expected:

> sf_use_s2(FALSE)
Spherical geometry (s2) switched off
> can_clipped2 <- st_intersection(can, clip)
although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning message:
attribute variables are assumed to be spatially constant throughout all geometries 

enter image description here

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.