st_crop and st_intersection do not cut decimal degrees assuming planar cartesian coordinates

[MikkoVihtakari]

Apologies for opening a new issue. Please close if already solved.

I thought this was a feature of sf and therefore made an SO question: https://stackoverflow.com/questions/69027136/how-to-crop-non-projected-epgs4326-polygons-along-straight-lines

After studying the behavior more carefully, perhaps this is, at least, an unwanted feature?

Run using the CRAN version of sf 1.0-2. and R version 4.1.1 (2021-08-10); Platform: x86_64-apple-darwin17.0 (64-bit); Running under: macOS Big Sur 10.16

[edzer]

On the globe, there are no straight lines. Since sf 1.0-0, "lines" between points for ellipsoidal (geographic) coordinates are great circle distances, not straight lines in Plate Carree. You can set sf_use_s2(FALSE) to get the pre-1.0 behavior (use GEOS rather than s2geometry). See also https://r-spatial.org/r/2020/06/17/s2.html

[MikkoVihtakari]

Thanks for the clarification. Flat Earth was indeed practical when trying to break things into regions but an old dog can learn new tricks. sf_use_s2(FALSE) seems to work with st_intersection but not with st_crop. Is/will there (be) a way to add sf_use_s2(FALSE) as an argument inside the st_intersection and st_crop functions? Round earth seems like a great idea but right now the cropping is confusing and works suboptimally. Turning the feature off only during cropping would be nice.

[jlacko]

This is an interesting problem; I feel that there may be an issue with the plot part of the linked question.

Consider this code: the red polygon is taken from the SO question, and when plotted looks square (i.e. Plate Caree). Of the two cities - Being and Hanoi, sorry for not hardcoding the coordinates - one appears inside and the other outside of the polygon.

But when I run sf::st_contains() on the cities both are inside of the polygon. This makes perfect sense in context of the s2 behavior of the polygon, where the top of the polygon is actually a section of a great circle connecting the two points, but then the plot feels wrong.

library(sf)
library(dplyr)
library(spData)

red_polygon <- matrix(c(60, 0, 180, 0, 180, 30, 60, 30, 60, 0), byrow = TRUE, ncol = 2) %>%
  list() %>% 
  st_polygon() %>% 
  st_sfc() %>% 
  st_set_crs(4326) 

cities <- tidygeocoder::geo(c("Bejing, China", "Hanoi, Vietnam")) %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326)

plot(st_geometry(world), xlim = c(60,180), ylim = c(0, 60))
plot(red_polygon, add = T, border = "red")
plot(st_geometry(cities), add = T, col = "blue", pch = 4)

st_contains(red_polygon, cities, sparse = F)
#     [,1] [,2]
# [1,] TRUE TRUE

[edzer]

library(sf)
library(dplyr)
library(spData)

red_polygon <- matrix(c(60, 0, 179.9, 0, 179.9, 30, 60, 30, 60, 0), byrow = TRUE, ncol = 2) %>%
  list() %>% 
  st_polygon() %>% 
  st_sfc() %>% 
  st_set_crs(4326) %>%
  st_segmentize(units::set_units(1, degree))

cities <- tidygeocoder::geo(c("Bejing, China", "Hanoi, Vietnam")) %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326)

plot(st_geometry(world), xlim = c(60,180), ylim = c(0, 60))
plot(red_polygon, add = T, border = "red")
plot(st_geometry(cities), add = T, col = "blue", pch = 4)

[jlacko]

Exactly, and the question kind of morphs from an unexpectedly looking intersection to "how does a polygon specified by four points on a sphere look on a map" (having the bottom section sat on the equator spoils half of the fun).

[MikkoVihtakari]

I feel that st_crop philosophy has not been brought up to the round-Earth thinking. How about using the bbox limits something like in A of this figure?

Caption for the figure can be found from here.

[edzer]

For plotting, orthographic projections seem to work better for smaller areas, but also not really for an area this size (although the red one includes the points):

library(sf)
library(dplyr)
library(spData)

red_polygon <- matrix(c(60, 0, 179.9, 0, 179.9, 30, 60, 30, 60, 0), byrow = TRUE, ncol = 2) %>%
  list() %>% 
  st_polygon() %>% 
  st_sfc() %>% 
  st_set_crs(4326) 

cities <- tidygeocoder::geo(c("Bejing, China", "Hanoi, Vietnam")) %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326)

png("x.png")
crs = st_crs("+proj=ortho +lon_0=120 +lat_0=30")
plot(st_geometry(st_transform(world, crs)))
plot(st_transform(red_polygon, crs), add = T, border = "red")
plot(st_geometry(st_transform(cities, crs)), add = T, col = "blue", pch = 4)
plot(st_transform(st_segmentize(red_polygon, units::set_units(1, degree)), crs), add = T, border = 'green')

[tim-salabim]

Yet another reason to have a cesium port for R

[MikkoVihtakari]

Should this issue be closed or left open? I have a few suggestions on how to address it:

1. Add an attribute to st_crop and st_intersection which allows turning spherical geometry off during the operation.
2. Modify the xmin, ymin, xmax, ymax thinking of bbox arguments within st_crop
3. Add "arguments passed on to s2_options" to ... of st_crop like in st_intersection