4

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

11

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

2
  • how do you return features that intersect a polygon without actually clipping them?
    – jjunju
    Jun 16, 2023 at 9:45
  • 1
    st_intersects returns a matrix or list of intersections, you can then extract the intersecting features by subsetting. See the help for details, search for more, or ask a complete new question if you get stuck.
    – Spacedman
    Jun 16, 2023 at 10:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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