search

r-spatial / s2

Spherical Geometry Operators Using the S2 Geometry Library
https://r-spatial.github.io/s2/
70 stars 10 forks source link

Given edges do not form loops (indegree != outdegree) #83

Closed edzer closed 3 years ago

edzer commented 3 years ago

The example comes from here: https://github.com/r-spatial/sf/issues/1510

library(sf)
# Linking to GEOS 3.8.1, GDAL 3.1.3, PROJ 7.1.1
library(s2)
download.file("https://www.dropbox.com/s/izol20rxf90apye/Areas.gpkg?dl=1",
              destfile = "Areas.gpkg" , mode='wb')
# trying URL 'https://www.dropbox.com/s/izol20rxf90apye/Areas.gpkg?dl=1'
# Content type 'application/binary' length 11182080 bytes (10.7 MB)
# ==================================================
# downloaded 10.7 MB

Boundary1 <- st_transform(st_read("Areas.gpkg"),crs=4326)
# Reading layer `Areas' from data source `/tmp/Areas.gpkg' using driver `GPKG'
# Simple feature collection with 1000 features and 1 field
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -116.1928 ymin: 68218.7 xmax: 640329.7 ymax: 868949.7
# projected CRS:  OSGB 1936 / British National Grid

download.file("https://www.dropbox.com/s/z8pham6w9ldesfy/Local_Authority_Districts__May_2020__Boundaries_UK_BFE-shp.zip?dl=1",
              destfile = "Local_Authority_Districts__May_2020__Boundaries_UK_BFE-shp.zip" , mode='wb')
# trying URL 'https://www.dropbox.com/s/z8pham6w9ldesfy/Local_Authority_Districts__May_2020__Boundaries_UK_BFE-shp.zip?dl=1'
# Content type 'application/binary' length 34702508 bytes (33.1 MB)
# ==================================================
# downloaded 33.1 MB

unzip("Local_Authority_Districts__May_2020__Boundaries_UK_BFE-shp.zip", exdir = ".")
Boundary2 <- st_transform(st_read("Local_Authority_Districts__May_2020__Boundaries_UK_BFE.shp"),crs=4326)
# Reading layer `Local_Authority_Districts__May_2020__Boundaries_UK_BFE' from data source `/tmp/Local_Authority_Districts__May_2020__Boundaries_UK_BFE.shp' using driver `ESRI Shapefile'
# Simple feature collection with 379 features and 10 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -116.1928 ymin: 5333.81 xmax: 655989 ymax: 1220302
# projected CRS:  OSGB 1936 / British National Grid

b1 = st_as_s2(Boundary1)
b2 = st_as_s2(Boundary2)
i = lapply(b1, function(x) s2_intersection(x, b2))
ii = lapply(b1, function(x) s2_intersection(x, b2, s2::s2_options(model = "closed")))
# Error in cpp_s2_intersection(recycled[[1]], recycled[[2]], options) : 
#   Given edges do not form loops (indegree != outdegree)
# Calls: lapply ... s2_intersection -> new_s2_xptr -> cpp_s2_intersection
# Execution halted

the model = "closed" is what sf does by default.

paleolimbot commented 3 years ago

Slightly simplified:

library(sf)
#> Warning: package 'sf' was built under R version 4.0.2
#> Linking to GEOS 3.8.1, GDAL 3.1.1, PROJ 6.3.1
library(s2)

download.file(
  "https://www.dropbox.com/s/izol20rxf90apye/Areas.gpkg?dl=1",
  destfile = "Areas.gpkg",
  mode='wb'
)

Boundary1 <- st_transform(st_read("Areas.gpkg"), crs = 4326)
#> Reading layer `Areas' from data source `/private/var/folders/bq/2rcjstv90nx1_wrt8d3gqw6m0000gn/T/Rtmp5BSc0e/reprex8a8b8dc3a26/Areas.gpkg' using driver `GPKG'
#> Simple feature collection with 1000 features and 1 field
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -116.1928 ymin: 68218.7 xmax: 640329.7 ymax: 868949.7
#> projected CRS:  OSGB 1936 / British National Grid

download.file(
  "https://www.dropbox.com/s/z8pham6w9ldesfy/Local_Authority_Districts__May_2020__Boundaries_UK_BFE-shp.zip?dl=1",
  destfile = "Local_Authority_Districts__May_2020__Boundaries_UK_BFE-shp.zip" , 
  mode='wb'
)

unzip("Local_Authority_Districts__May_2020__Boundaries_UK_BFE-shp.zip", exdir = ".")
Boundary2 <- st_transform(st_read("Local_Authority_Districts__May_2020__Boundaries_UK_BFE.shp"),crs=4326)
#> Reading layer `Local_Authority_Districts__May_2020__Boundaries_UK_BFE' from data source `/private/var/folders/bq/2rcjstv90nx1_wrt8d3gqw6m0000gn/T/Rtmp5BSc0e/reprex8a8b8dc3a26/Local_Authority_Districts__May_2020__Boundaries_UK_BFE.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 379 features and 10 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -116.1928 ymin: 5333.81 xmax: 655989 ymax: 1220302
#> projected CRS:  OSGB 1936 / British National Grid

plot(c(Boundary1$geom[22], Boundary2$geometry[345]))


b1 = st_as_s2(Boundary1)
b2 = st_as_s2(Boundary2)

# ok
st_intersection(Boundary1$geom[22], Boundary2$geometry[345])
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Geometry set for 1 feature 
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -2.456499 ymin: 56.88294 xmax: -2.130036 ymax: 57.06899
#> geographic CRS: WGS 84
#> MULTIPOLYGON (((-2.254012 56.89502, -2.26544 56...

# ok
s2_intersection(b1[22], b2[345], options = s2_options(model = "semi-open"))
#> <s2_geography[1]>
#> [1] <GEOMETRYCOLLECTION (MULTILINESTRING ((-2.19669 56.9451, -2.19669 56.9451, -2.19669 56.9451), (-2.19445 56.959, -2.19445 56.959...>

# fails
s2_intersection(b1[22], b2[345], options = s2_options(model = "closed"))
#> Error in cpp_s2_intersection(recycled[[1]], recycled[[2]], options): Given edges do not form loops (indegree != outdegree)

# fails
s2_intersection(
  b1[22] %>% s2_snap_to_grid(1e-9), 
  b2[345] %>% s2_snap_to_grid(1e-9), 
  options = s2_options(model = "closed")
)
#> Error in cpp_s2_intersection(recycled[[1]], recycled[[2]], options): Given edges do not form loops (indegree != outdegree)

# ok
s2_intersection(
  b1[22] %>% s2_snap_to_grid(1e-8), 
  b2[345] %>% s2_snap_to_grid(1e-8), 
  options = s2_options(model = "closed")
)
#> <s2_geography[1]>
#> [1] <MULTIPOLYGON (((-2.20008 56.9588, -2.20017 56.9588, -2.2003 56.9588, -2.20043 56.9588, -2.20057 56.9587...>

Created on 2020-10-15 by the reprex package (v0.3.0)

Setting precision from the s2 end seems to work (I didn't try to see if rounding using st_set_precision() did the trick as well). It looks like the intersections that fail occur when there is a common edge (that in reality is slightly different).