s2_intersection creates invalid geometries

[edzer]

See r-spatial/sf#1810 ; reprex through sf:

library(sf)
library(giscoR)
library(dplyr)

lander <- gisco_get_nuts(country = "DE", 
                         nuts_level = 1) %>% 
  mutate(value = runif(n(), min = 100, max = 500)) # some data to interpolate ...

grid <- st_make_grid(st_bbox(lander),
                     n = c(10, 10),
                     square = F) %>% 
  st_as_sf() %>% 
  mutate(id = 1:n()) # polygon ids
all(st_is_valid(st_intersection(lander, grid))) 

[paleolimbot]

I think this is because when the S2Point (xyz on the sphere) is converted to longitude/latitude, some sort of precision loss occurs and the polygon becomes invalid. I can trigger this without sf as well.

The workaround is to set the precision for the intersection operation (I don't know if there's a way to do this in sf). The value can be very small and it will apparently still generate valid geometries.

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.2.3, PROJ 7.2.1
library(giscoR)
library(dplyr, warn.conflicts = FALSE)
library(s2)

lander <- gisco_get_nuts(country = "DE", 
                         nuts_level = 1) %>% 
  mutate(value = runif(n(), min = 100, max = 500)) # some data to interpolate ...

all(s2::s2_is_valid(lander$geometry))
#> [1] TRUE

grid <- st_make_grid(st_bbox(lander),
                     n = c(10, 10),
                     square = F) %>% 
  st_as_sf() %>% 
  mutate(id = 1:n()) # polygon ids


all(s2::s2_is_valid(grid$x))
#> [1] TRUE

geoms <- expand.grid(
  lander_id = seq_len(nrow(lander)),
  grid_id = seq_len(nrow(grid))
) %>% 
  mutate(
    lander_geom_s2 = as_s2_geography(lander$geometry)[lander_id],
    grid_geom_s2 = as_s2_geography(grid$x)[grid_id]
  )

all(s2::s2_is_valid(geoms$lander_geom_s2))
#> [1] TRUE
all(s2::s2_is_valid(geoms$grid_geom_s2))
#> [1] TRUE

geoms$intersection <- s2::s2_intersection(
  geoms$lander_geom_s2,
  geoms$grid_geom_s2
)

all(s2::s2_is_valid(geoms$intersection))
#> [1] TRUE

geoms$geom <- wk::as_wkb(geoms$intersection)
all(s2::s2_is_valid(geoms$geom))
#> [1] FALSE

geoms$intersection <- s2::s2_intersection(
  geoms$lander_geom_s2,
  geoms$grid_geom_s2, options = s2_options(snap = s2_snap_precision(1e-16))
)

geoms$geom <- wk::as_wkb(geoms$intersection)
all(s2::s2_is_valid(geoms$geom))
#> [1] TRUE

^{Created on 2022-05-21 by the reprex package (v2.0.1)}

[paleolimbot]

In #185 I added the ability to export/import points in S2's unit vector space. It's not 100% fixing this issue, but it at least allows a path forward to "store" s2 geometries outside the external pointer structure (e.g., in an sf object). This should at least make it possible to circumvent precision loss by roundtripping to longitude/latitude.

library(s2)
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE

wk::wk_handle(
  s2_data_countries("Fiji"),
  wk::sfc_writer(),
  s2_projection = NULL
)
#> Geometry set for 1 feature 
#> Geometry type: MULTIPOLYGON
#> Dimension:     XYZ
#> Bounding box:  xmin: -0.9611549 ymin: -0.003467136 xmax: -0.9488733 ymax: 0.04511876
#> z_range:       zmin: -0.3137934 zmax: -0.2759877
#> CRS:           NA
#> MULTIPOLYGON Z (((-0.9541688 0.02709235 -0.2980...

# can also transform using a wk_trans
wk::wk_transform(
  wk::wkt("POLYGON Z ((1 0 0, 0 1 0, 0 0 1, 1 0 0))"),
  s2_trans_lnglat()
)
#> <wk_wkt[1]>
#> [1] POLYGON ((0 0, 90 0, 0 90, 0 0))

wk::wk_transform(
  sf::st_as_sfc(s2_data_countries("Fiji")),
  s2_trans_point()
)
#> Geometry set for 1 feature 
#> Geometry type: MULTIPOLYGON
#> Dimension:     XYZ
#> Bounding box:  xmin: -0.9611549 ymin: -0.003467136 xmax: -0.9488733 ymax: 0.04511876
#> z_range:       zmin: -0.3137934 zmax: -0.2759877
#> CRS:           NA
#> MULTIPOLYGON Z (((-0.9541688 0.02709235 -0.2980...

^{Created on 2022-06-23 by the reprex package (v2.0.1)}