sf icon indicating copy to clipboard operation
sf copied to clipboard

st_difference broken with S2 on?

Open courtiol opened this issue 4 years ago • 2 comments

Updating to sf 1.0 or 1.0-1 breaks some of my code when (and only when S2 is on). One issue is that sf::st_difference now behave differently:

library(sf)
# Linking to GEOS 3.9.0, GDAL 3.2.2, PROJ 7.2.1

m1
# Geometry set for 1 feature 
# Geometry type: MULTIPOLYGON
# Dimension:     XY
# Bounding box:  xmin: -5.720722 ymin: 37.65219 xmax: 50.37008 ymax: 57.58809
# CRS:           +proj=longlat +datum=WGS84 +no_defs 
# MULTIPOLYGON (((-4.915853 52.03233, -4.91627 52...

m2
# Geometry set for 1 feature 
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 90
# CRS:           +proj=longlat +datum=WGS84 +no_defs 
# POLYGON ((-180 90, 180 90, 180 -90, -180 -90, -...

sf_use_s2(FALSE)
m3 <- st_difference(m2, m1)

sf_use_s2(TRUE)
m4 <- st_difference(m2, m1)

par(mfrow = c(2, 2))
plot(m1, col = 1, main = "Mask 1")
plot(m2, col = 1, main = "Mask 2")
plot(m3, col = 1, main = "Mask 2 - Mask 1 without S2")
plot(m4, col = 1, main = "Mask 2 - Mask 1 with S2")

image

PS: please do let me know if you need data to reproduce this (I failed to save them but should be able to recreate them).

courtiol avatar Jun 27 '21 08:06 courtiol

@courtiol , I had the same issue. For what it is worth, using st_make_valid() on both objects resolved the error in my case.

k6adams avatar Jun 28 '21 17:06 k6adams

It's still more clumsy with s2 "on", and needs some work to make this run smoother:

# using sf with s2 backend:                                                           
library(rnaturalearth)
library(sf)
library(s2)
co = ne_countries(returnclass = "sf") %>% 
        st_make_valid()
bb = st_as_sfc(st_bbox(c(xmin = -179.9, xmax = 179.9, ymin = -89.9, ymax = 89.9),
    crs = st_crs(co))) %>% st_segmentize(units::set_units(100, km))
d = s2_difference(as_s2_geography(bb, oriented = TRUE),
    s2_union_agg(as_s2_geography(co)))
plot(st_as_sfc(d), col = 'blue')

x

  • the bounding box polygon is by default (oriented = FALSE) reversed: the assumption is that polygons larger than half the globe have their ring direction wrong; we'd need some mechanism to override this default so that one can call st_difference under the assunmption that oriented = TRUE is used;
  • the bounding box needs to be slightly smaller than the max range, so that transforming back to R2 (only for plotting!) works.

edzer avatar Jun 28 '21 19:06 edzer