R – st_intersection from sf Package Not Clipping Territory Correctly

intersectionrsf

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

Best Answer

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

Related Question