Sf: st_intersection() produces error: TopologyException: found non-noded intersection

[markyxe]

Hello,

First, thank you for making sf, it's really great! I'm running into a problem, however. I'm trying to intersect multiple overlapping polygons using st_intersection() and getting the following error:

Error in CPL_nary_intersection(x) :
Evaluation error: TopologyException: found non-noded intersection between LINESTRING (-105.084 51.9804, -105.084 51.9804) and LINESTRING (-105.084 51.9804, -105.084 51.9804) at -105.08405553932499 51.980378885448744.

Following advice on github and stackoverflow, I've tried zero buffering the polygons, and also using st_st_precision() and st_make_valid() before the st_intersection() call, but neither of those approaches worked consistently (the latter approach works for some sets of polygons, but not others).

I can intersect the polygons using the raster package, and also in ArcGIS and QGIS, so the features seem okay, but they aren't playing nicely with sf. The example posted below is for a set of three particular polygons, but I'm experiencing this with other sets of polygons too.

Here is a simplified, reproducible example. I'm sorry I can't seem to post a csv here but I saved the input data as an xlsx.

Thank you for any help or advice. Best wishes,

Mark

points <-read.csv("dat.csv") # please see attached dat.xlsx
polys<- points %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326) %>%
  group_by(SitePondGpsRep) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON")
polys<-tidyr::separate(polys,SitePondGpsRep,"-",into=c("Site","Pond","Gps","Rep"),remove=F)
polys$Gps<-as.numeric(polys$Gps)
polys$SitePond<-paste(polys$Site,polys$Pond)
polys<-filter(polys,Gps==1|Gps==2|Gps==3) #select the three problematic polygons
polys_intersect<-st_intersection(polys)

points.xlsx

[edzer]

Thanks. Try setting a useful value for precision, e.g.

polys %>% st_set_precision(1e5) %>% st_intersection

[markyxe]

Hi Edzer, thanks for this suggestion. I did try setting precision before the intersection, and the value you suggested does work for the specific set of polygons I posted, but it doesn't work for other sets of polygons...so setting precision doesn't seem like a general solution. (My code performs 1,375 intersections of ~500 polygons so I'd prefer to hold precision constant across all of these operations.) Any other ideas? Many thanks for your help.

[edzer]

No other ideas here.

[markyxe]

Hi Edzer, here are examples of when setting precision works and doesn't work. However, the error I'm getting in the second case (SitePond=BURR 105) is different, now it refers to a ring self-intersection (whereas for SitePond = BURR 10 the error is due a non-noded intersection). Is there some way to clean up the geometry to avoid these types of problems, before calling st_intersection()? Thank you.

library(tidyverse)
library(dplyr)
points <-read.csv("points.csv") #please use uploaded file points.xlsx
polys<- points %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326) %>%
  group_by(SitePondGpsRep) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON")
polys<-tidyr::separate(polys,SitePondGpsRep,"-",into=c("Site","Pond","Gps","Rep"),remove=F) #separate id into multiple columns
polys$Gps<-as.numeric(polys$Gps)
polys$SitePond<-paste(polys$Site,polys$Pond)

#example where setting precision works
polys %>% filter(SitePond=="BURR 10", Gps==1|Gps==2|Gps==3) %>% st_intersection #doesn't work
polys %>% filter(SitePond=="BURR 10", Gps==1|Gps==2|Gps==3) %>% st_set_precision(1e5) %>% st_intersection #works!

#example where setting precision DOESN'T work
polys %>% filter(SitePond=="BURR 105", Gps==5|Gps==6) %>% st_intersection #doesn't work
polys %>% filter(SitePond=="BURR 105", Gps==5|Gps==6) %>% st_set_precision(1e5) %>% st_intersection #still doesn't work

points.xlsx

[markyxe]

HI Edzer. Further to this, I followed your advice here:
https://www.r-spatial.org/r/2017/03/19/invalid.html

and detected many invalid geometries which I fixed with st_make_valid. However, st_intersection still does not work on some of these features, even though they are now valid. In some cases, altering precision will make it work, but not in all cases. Do you have any advice for how to proceed? Thank you.

[edzer]

No. Maybe reconsider what you're trying to do.

[mdsumner]

Perhaps reduce by intersecting pairs sequentially?

points <-readxl::read_excel("points.xlsx") # please see attached dat.xlsx
polys<- points %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326) %>%
  group_by(SitePondGpsRep) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON")
g <- purrr::reduce(st_geometry(polys), st_intersection)


plot(st_geometry(polys))

plot(g, border = "red", add = TRUE)

(g is now purely a sfg, not a vector of those or sf data frame).

[edzer]

Any progress here? Not sure sf can fix this, or whether this points to an issue with sf, rather than with your data.

[jebyrnes]

BTW: I tried the polygon cast and then reduce idea with another data set, and it yields

GEOMETRYCOLLECTION EMPTY

Try it below:

test_set <- devtools::source_gist("c95701bd444cda3e342838fd9a090fb3",
                                  filename = "test_set.R")$value

polys <- test_set  %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON")

plot(polys)

g <- purrr::reduce(st_geometry(polys), st_intersection)

g

[dvictori]

I'm also struggling with this issue. Tried setting the precision; using st_make_valid; filter only valid geometries but in the end, I could not get st_intersection it to work with the very (very) bad shapefiles that I've got.

What I ended up doing was first run st_intersects() to get a list of the polygons that will intersect with another. Then I will run st_intersection() for each item of the list. This allows me to identify the group of features that are giving me a headache. The funny thing is that I separated a group of geometries that st_intersection will not work and they are all valid. I'm attaching a shapefile with the troublesome geometries, as an example.

trouble_geom.zip

[rsbivand]

This is most likely a further case of this thread: https://lists.osgeo.org/pipermail/geos-devel/2019-January/008794.html. See the explanations by Paul Ramsey and Martin Davis. This is a direct result of the use of Simple Features rather than an ArcNode representation, further complicated by double precision not being absolute precision on the real line. As Paul and Martin relate, operations may introduce points that "miss" the lines they should be part of in a world with absolute precision. Your troubled polygons are:

I don't see precision/scale issues in rgeos, but I don't know what you are trying to do.

[dvictori]

Thanks for the input @rsbivand. It looks like I'll need to rethink my strategy.

The polygons I sent are part of a very messy database, with lots and lots of overlaps. Right now I'm interested in estimating how much overlap there is and, if possible, identifying those trouble areas.

For the overlap amount, I can calculate the difference between the sum(st_area(trouble_geom) and st_area(st_union(trouble_geom)). For the polygons in question, the overlap is about 6.2%.

But, in order to visually identify the trouble areas, I though about using sf::st_intersection, which brought me to this issue. For what it's worth, I'm able to get the overlap areas using the following SQL Query in QGis DB Manager:

select st_intersection(a.geometry, b.geometry)
from trouble_geom a, trouble_geom b
where a.FID != b.FID

The image bellow shows the results. Now, one thing that made me scratch my head. Geometry operations in SQL also rely on Simple Features right? So this issue with sf::st_intersection should also affect st_intersection in SQL.

[mdsumner]

I see, decompose to triangles in a single mesh with all edges and you can summarize n_intersections with the space bucket:

  library(sf)
#> Warning: package 'sf' was built under R version 3.5.2
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
## remotes::install_github("mdsumner/spacebucket")
  library(spacebucket)
  geom <- sf::read_sf("C:\\temp\\trouble\\trouble_geom.shp")
  bucket <- spacebucket(geom)
  par(mfrow = c(3, 2))
  for (i in seq(nrow(geom), 1)) {
    i_overlap <- n_intersections(bucket, i)
   plot(st_geometry(geom), reset = FALSE, main = sprintf("%i overlaps\n %f", i, sum(st_area(i_overlap))))

   plot(i_overlap[1], add = TRUE, reset = FALSE)

  }

  ## total area including overlaps
  sum(st_area(geom))
#> 1352693856 [m^2]

  ## unioned area 
  st_area(st_union(geom))
#> 1273933403 [m^2]
  sum(st_area(n_intersections(bucket, 1)))
#> [1] 1273933403

^{Created on 2019-02-06 by the reprex package (v0.2.1)}

Super raw and unpolished and only on github, but the "bucket" keeps a track of all available elements (triangles) and which objects they are part of.

[rsbivand]

With reference to sf and topological operations, the common element isn't so much the Simple Features definition, it is the use of the GEOS C++ library (based on JTS/Java) in sf, in PostGIS, in GDAL, and most likely in QGIS (and in parts of GRASS). So you actually have to try rather different software to find a different implementation of planar topological operations.

[lbusett]

I recently found out that it is somehow possible to solve/mitigate the problem by snapping the object on itself before the intersection:

suppressPackageStartupMessages(library(sf))
geoms <- sf::st_read("/home/lb/Downloads/trouble_geom.shp", quiet = TRUE)
geoms <- st_set_precision(geoms, 1000000)
geom_int <- st_snap(geoms, geoms, tolerance = 5) %>% 
  st_intersection() %>% 
  st_collection_extract("POLYGON") %>%
  dplyr::group_by(n.overlaps) %>% 
  dplyr::summarise() %>% 
  dplyr::mutate(area = sf::st_area(.))

geom_int
#> Simple feature collection with 5 features and 2 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 5630707 ymin: 8219579 xmax: 5708479 ymax: 8251878
#> epsg (SRID):    NA
#> proj4string:    +proj=poly +lat_0=0 +lon_0=-54 +x_0=5000000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs
#> # A tibble: 5 x 3
#>   n.overlaps area                                                  geometry
#>        <int> <S3: units>                                 <MULTIPOLYGON [m]>
#> 1          1 1223314090… (((5694044 8232024, 5694055 8232064, 5694052 8232…
#> 2          2   27035986… (((5693895 8232747, 5693893 8232746, 5693895 8232…
#> 3          3   20608066… (((5694129 8235259, 5694129 8235259, 5694129 8235…
#> 4          4    1407664… (((5693848 8232728, 5693848 8232728, 5693849 8232…
#> 5          5    1578034… (((5694126 8235825, 5694126 8235825, 5694126 8235…
sum(geom_int$area)
#> 1273943841 [m^2]
plot(geom_int["n.overlaps"])

I am not really sure WHY this works. In my head, I figured it as a kind of "snap_to_grid" with a dynamic grid "defined" by the vertexes of the object itself, but I do not know if that's the truth and I did not give it further thought....

Results in terms of intersection areas are very similar to @mdsumner ones. Note however that here we do not have the 6-fold intersections (which are removed by the snapping, I think - this may be an added advantage because we get an auto-removal of small slivers? ). Also, this is still a bit "brittle": for example, it works with most tolerances between 5 and 20, while it fails with tolerance < 5, or 6 with the uisual topology exceptions (I was lucky enough to randomly set tolerance to 10 on the first run....) . Resulting polys are also dependent on tolerance.

HTH!

^{Created on 2019-02-07 by the reprex package (v0.2.1)}