st_union breaks when some edges are crossing, but only when S2 is on

[courtiol]

I have noticed that st_union() no longer works with some of the polygons I have when s2 is on.

This seems to be because the polygons have some genuine glitches, which make the internal call to st_as_s2.sfc() crash.

Using st_as_s2() with the argument rebuild = TRUE seems to avoid the problem (see below), but st_union() won't forward the parameters till that call.

They probably is a better way altogether to fix polygons (ideally a warning or the doc could suggest what to do), but I thought it could be useful to point to somewhat unexcepted behaviours caused by activating S2.

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

mask_future <- st_read("mask_future.shp")$geometry
#> Reading layer `mask_future' from data source 
#>   `mask_future.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 6 features and 22 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -5.185555 ymin: 17.47994 xmax: 56.09433 ymax: 63.4352
#> Geodetic CRS:  WGS 84

sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off
union_sf <- st_union(mask_future)
#> although coordinates are longitude/latitude, st_union assumes that they are planar
union_sf
#> Geometry set for 1 feature 
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -5.185555 ymin: 17.47994 xmax: 56.09433 ymax: 63.4352
#> Geodetic CRS:  WGS 84
#> MULTIPOLYGON (((-4.707389 55.62114, -4.698309 5...

sf_use_s2(TRUE)
#> Spherical geometry (s2) switched on
st_union(mask_future)
#> Error in s2_geography_from_wkb(x, oriented = oriented, check = check): Evaluation error: Found 2 features with invalid spherical geometry.
#> [3] Loop 0 is not valid: Edge 1604 crosses edge 1606
#> [5] Loop 0 is not valid: Edge 18946 crosses edge 18948.

traceback()
#> 15: stop(structure(list(message = "Evaluation error: Found 2 features with invalid spherical geometry.\n[3] Loop 0 is not valid: #> Edge 1604 crosses edge 1606\n[5] Loop 0 is not valid: Edge 18946 crosses edge 18948.", 
#>         call = s2_geography_from_wkb(x, oriented = oriented, check = check), 
#>         cppstack = NULL), class = c("Rcpp::eval_error", "C++Error", 
#>      "error", "condition")))
#> 14: s2_geography_from_wkb(x, oriented = oriented, check = check)
#> 13: new_s2_xptr(s2_geography_from_wkb(x, oriented = oriented, check = check), 
#>         "s2_geography")
#> 12: as_s2_geography.WKB(st_as_binary(x), ..., oriented = oriented)
#> 11: s2::as_s2_geography(st_as_binary(x), ..., oriented = oriented)
#> 10: st_as_s2.sfc(x, ..., oriented = oriented)
#> 9: as_s2_geography.sfc(x)
#> 8: as_s2_geography(x)
#> 7: s2_union(x, options = options)
#> 6: cpp_s2_union_agg(s2_union(x, options = options), options, na.rm)
#> 5: new_s2_xptr(cpp_s2_union_agg(s2_union(x, options = options), 
#>        options, na.rm), "s2_geography")
#> 4: s2::s2_union_agg(x, ...)
#> 3: st_as_sfc(s2::s2_union_agg(x, ...), crs = st_crs(x))
#> 2: st_union.sfc(mask_future)
#> 1: st_union(mask_future)

union_s2 <- st_as_sf(s2::s2_union_agg(st_as_s2(mask_future, rebuild = TRUE)))
union_s2
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -5.185555 ymin: 17.47994 xmax: 56.09433 ymax: 63.4352
#> Geodetic CRS:  WGS 84
#>                         geometry
#> 1 MULTIPOLYGON (((17.10703 57...

sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off
library(ggplot2)
ggplot() +
  geom_sf(data = st_difference(union_sf, union_s2), colour = "red", size = 5) +
  geom_sf(data = st_difference(union_s2), colour = "blue", alpha = 0.5) +
  coord_sf(ylim = c(48, 65))
#> although coordinates are longitude/latitude, st_difference assumes that they are planar
#> although coordinates are longitude/latitude, st_difference assumes that they are planar

I am attaching the files used here: mask_future.zip

[courtiol]

In passing, I noticed that the content of union_sf and union_s2 are equivalent (beyond the differences shown in red in the plot), but interestingly the polygons are not sorted in the same order.

I have no idea if that matters for some people or not... (it does not for me).

> str(union_s2)
Classes ‘sf’ and 'data.frame':  1 obs. of  1 variable:
 $ geometry:sfc_MULTIPOLYGON of length 1; first list element: List of 13
  ..$ :List of 1
  .. ..$ : num [1:880, 1:2] 17.1 17.1 17.1 17.1 17.1 ...
  ..$ :List of 1
  .. ..$ : num [1:126, 1:2] 56.1 56.1 56.1 56.1 56 ...
  ..$ :List of 1
  .. ..$ : num [1:6244, 1:2] -2.8 -2.81 -2.81 -2.81 -2.82 ...
  ..$ :List of 1
  .. ..$ : num [1:4, 1:2] 0.108 0.108 0.112 0.108 53.632 ...
  ..$ :List of 1
  .. ..$ : num [1:4, 1:2] -4.55 -4.55 -4.55 -4.55 52.79 ...
  ..$ :List of 1
  .. ..$ : num [1:1670, 1:2] 18.8 18.8 18.8 18.8 18.8 ...
  ..$ :List of 1
  .. ..$ : num [1:20803, 1:2] 11.8 11.8 11.8 11.8 11.8 ...
  ..$ :List of 1
  .. ..$ : num [1:4, 1:2] 14.3 14.3 14.3 14.3 56 ...
  ..$ :List of 1
  .. ..$ : num [1:6649, 1:2] 21.7 21.7 21.7 21.7 21.7 ...
  ..$ :List of 1
  .. ..$ : num [1:7, 1:2] 21.8 21.8 21.8 21.8 21.8 ...
  ..$ :List of 1
  .. ..$ : num [1:7, 1:2] 11.8 11.8 11.8 11.8 11.8 ...
  ..$ :List of 1
  .. ..$ : num [1:7, 1:2] 17.6 17.6 17.6 17.6 17.6 ...
  ..$ :List of 1
  .. ..$ : num [1:3927, 1:2] 20.9 24.6 29.7 33.8 48.3 ...
  ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: 
  ..- attr(*, "names")= chr(0) 

> str(union_sf)
sfc_MULTIPOLYGON of length 1; first list element: List of 12
 $ :List of 1
  ..$ : num [1:9, 1:2] -4.71 -4.7 -4.7 -4.7 -4.7 ...
 $ :List of 1
  ..$ : num [1:6246, 1:2] -2.79 -2.79 -2.79 -2.78 -2.78 ...
 $ :List of 1
  ..$ : num [1:7, 1:2] 11.8 11.8 11.8 11.8 11.8 ...
 $ :List of 1
  ..$ : num [1:20804, 1:2] 17.6 17.6 17.6 17.7 17.7 ...
 $ :List of 1
  ..$ : num [1:131, 1:2] 12.5 12.5 12.5 12.5 12.5 ...
 $ :List of 1
  ..$ : num [1:880, 1:2] 17.1 17.1 17.1 17.1 17.1 ...
 $ :List of 1
  ..$ : num [1:7, 1:2] 17.6 17.6 17.6 17.6 17.6 ...
 $ :List of 1
  ..$ : num [1:1670, 1:2] 18.8 18.8 18.8 18.8 18.8 ...
 $ :List of 1
  ..$ : num [1:3927, 1:2] 20 20 20 19.9 19.9 ...
 $ :List of 1
  ..$ : num [1:7, 1:2] 21.8 21.8 21.8 21.8 21.8 ...
 $ :List of 1
  ..$ : num [1:6649, 1:2] 21.8 21.8 21.8 21.8 21.8 ...
 $ :List of 1
  ..$ : num [1:126, 1:2] 56.1 56.1 56.1 56.1 56.1 ...
 - attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"

[edzer]

Yes, being valid on R2 doesn't mean something is valid on S2. This might help:

> st_is_valid(st_make_valid(st_set_precision(mask_future, 1e6)))
[1] TRUE TRUE TRUE TRUE TRUE TRUE

[edzer]

close?

[courtiol]

I have not check it yet, but yes let's close this for now. Many thanks!

[eastclintw00d]

@courtiol Did you check it?

I have a similar issue. st_is_valid(st_make_valid(st_set_precision(mydata, 1e6))) does not solve it. It worked fine on 0.9-8.

[courtiol]

I did and on the data above it does work:

unzip("~/Downloads/mask_future.zip", exdir = "~/Downloads/")
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.3.1, PROJ 8.1.0
mask_future <- st_read("~/Downloads/mask_future.shp")$geometry
#> Reading layer `mask_future' from data source 
#>   `/home/alex/Downloads/mask_future.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 6 features and 22 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -5.185555 ymin: 17.47994 xmax: 56.09433 ymax: 63.4352
#> Geodetic CRS:  WGS 84
st_union(st_make_valid(st_set_precision(mask_future, 1e6)))
#> Geometry set for 1 feature 
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -5.185555 ymin: 17.47994 xmax: 56.09433 ymax: 63.4352
#> Geodetic CRS:  WGS 84
#> MULTIPOLYGON (((17.10703 57.33286, 17.09697 57....

^{Created on 2021-08-16 by the reprex package (v2.0.0)}