st_intersection: found non-noded intersection between LINESTRING

[jepol77]

There seems to be an issue with st_intersection that I can not solve with the usual methods.

download polygons from here: https://github.com/jepol77/test/blob/main/invalid_polygons.rds
bbx <- readRDS(".../invalid_polygons.rds")
bbx <- bbx %>% st_set_precision(1e5) %>% st_make_valid()
i <- st_intersection(bbx)

I get this error: Fejl i CPL_nary_intersection(x) :
Evaluation error: TopologyException: found non-noded intersection between LINESTRING (497314 6.31343e+06, 497314 6.31343e+06) and LINESTRING (497314 6.31343e+06, 497314 6.31343e+06) at 497313.60609830072 6313426.1116788359.

I know there has been some recent work on this in PostGIS; any news when that is gonna trickle down to R?

[rsbivand]

Please see #794 and the extensive discussion there. You are attempting n-ary intersection, which recurses, intersecting intersection products. The input geometries are valid, but some intersection product ceases to be valid. Using GEOS-devel with OverlayNG does not help, although different invalid intermediate outputs are found. Further, your input objects, while being valid polygons, are actually:

[jepol77]

Could you put it in a for loop somehow, or would that make it very slow?

[rsbivand]

Not only very slow, but also using very much memory, apparently. N-ary intersection is not provided by JTS or GEOS; sf added it as a speculation, but because it grosses up imprecision as you reach greater depths of overlaps, it is not as such covered by the available computational geometry models. rgeos does not permit it. The depth seems to be up to 14 iterations, but I maxed out 16GB on the fifth with recusing binary st_intersection()s. Slicing out the categories as separate objects might help, but the number of possible combinations of n=34 is still very large (points in the surface intersecting 0, 1, combinations of 2, ... ) and so on. Maybe sampling with a very dense grid?

[rsbivand]

This gives 215 distinct areas:

bbx <- readRDS("all.rds")
library(sf)
bbxo <- st_union(bbx)
grd <- st_sample(bbxo, size=100000, type="regular")
grd_sp <- as(grd, "Spatial")
library(sp)
gridded(grd_sp) <- TRUE
grd_SP <- as(grd_sp, "SpatialPolygons")
grd_sf_p <- st_as_sf(grd_SP)
ovs <- st_intersects(grd_sf_p, st_geometry(bbx), sparse=FALSE)
codes <- apply(ovs, 1, function(x) {paste(sapply(x, ifelse, "1", "0"), collapse="")})
grdd_sf_p <- st_as_sf(grd_sf_p)
grdd_sf_p$codes <- codes
bbx_i <- aggregate(grdd_sf_p, list(grdd_sf_p$codes), head, 1)
bbx_i$fcodes <- factor(bbx_i$codes)
levels(bbx_i$fcodes)
dim(bbx_i) # [1] 215   4
i <- st_intersection(bbx, bbx_i[,"fcodes"])
dim(i) # [1] 417   3

Since the grid cells are about 1m^2, the intersection output less than 1 could perhaps be discarded.This is with about 100000 points, 1 million might also be feasible, but I feel that this is as good as you'll get.

[edzer]

This branch gives us

library(sf)
# Linking to GEOS 3.8.1, GDAL 3.1.3, PROJ 7.1.1
bbx <- readRDS("invalid_polygons.rds")
bbx <- bbx %>% st_set_precision(1e5) %>% st_make_valid()
(i <- st_intersection(bbx))
# Simple feature collection with 212 features and 3 fields
# geometry type:  GEOMETRY
# dimension:      XY
# bbox:           xmin: 495762 ymin: 6311742 xmax: 497599.9 ymax: 6314854
# projected CRS:  ETRS89 / UTM zone 32N
# First 10 features:
#     no_pnts n.overlaps origins                       geometry
# 1        82          1       1 POLYGON ((496433.2 6312980,...
# 2        82          1       2 MULTIPOLYGON (((497346 6314...
# 3        74          1       3 MULTIPOLYGON (((497260.1 63...
# 2.1      82          2    2, 4 POLYGON ((497297.5 6313184,...
# 4        85          1       4 MULTIPOLYGON (((497325.1 63...
# 2.2      82          3 2, 4, 5 POLYGON ((497295.3 6313180,...
# 4.1      85          2    4, 5 POLYGON ((497299.3 6313214,...
# 2.3      82          2    2, 5 MULTIPOLYGON (((497308 6313...
# 5        76          1       5 MULTIPOLYGON (((497281.9 63...
# 4.2      85          2    4, 6 MULTIPOLYGON (((497314.4 63...

basically by switching off GEOS error reporting and ignoring geometries causing trouble.

[edzer]

I'm still trying to count, and report, the amount of troubling geometries...

[rsbivand]

The nary branch and the grid give very similar output, with some missing combinations both ways (6.512942 [m^2] in nary intersection of 109888.1 [m^2], 28.57107 [m^2] in gridded binary intersection of 93167.42 [m^2]). I needed to aggregate my binary intersection as well, to combine separated geometries.

[edzer]

I merged the nary branch: better to return something than nothing with an error message.