issues with st_crop on s2_data_countries()

[MatthieuStigler]

When using st_crop() on st_as_sf(s2_data_countries()), I get some strange results. I faced the same problem with rnaturalearth (see issue ropensci/rnaturalearth#78), and assumed the problem was that rnaturalearth contained invalid geometries (which s2::s2_data_countries doesn't seem to contain), but it seems even with valid geometries, there is an issue?

Note I am posting here and might need to post instead on sf?

Thanks!

library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.0; sf_use_s2() is TRUE
library(s2)

co <- st_as_sf(s2_data_countries())

pol <- st_sfc(st_polygon(list(cbind(c(-120, -120, 30, 30, -120),
                                      c(30, 55, 55, 30, 30)))))

plot(co)
plot(pol, add=TRUE, border ="blue")
co |> 
  st_crop(st_bbox(pol)) |>
  plot(border="red", lwd=2, add=TRUE)

^{Created on 2023-05-25 with reprex v2.0.2}

[paleolimbot]

I'm pretty sure what's happening is that the box is getting interpreted to have a geodesic edge. This is almost certainly not what anybody doing st_crop() is asking for...sf could tessellate the box edges before sending to s2 just use GEOS?

[edzer]

Yes, something along these lines should work:

library(sf)                                                                           
#> Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.0; sf_use_s2() is TRUE                
library(s2)                                                                           
                                                                                      
co <- st_as_sf(s2_data_countries())                                                   
                                                                                      
pol <- st_sfc(st_polygon(list(cbind(c(-120, -120, 30, 30, -120),                      
                                      c(30, 55, 55, 30, 30)))))                       
                                                                                      
pol1 <- st_segmentize(pol, 1)                                                         
plot(co)                                                                              
plot(pol, add=TRUE, border ="blue")                                                   
st_crs(pol1) = st_crs(co)                                                             
co |>                                                                                 
  st_intersection(pol1) |>                                                            
  plot(border="red", lwd=2, add=TRUE)   

[edzer]

(It works because pol has an NA crs when passed to st_segmentize)