TopologyException in multipolygons (intersection function)

[carlosmolinaguerra]

Hi,
I have a large sfc (MULTIPOLYGON) objects which throw several TopologyException errors when I try to use the intersection function on them. I have a couple of questions so, If anyone has a moment, I would really appreciate it if you could take a look at this. I've attached the files in RDS format (all). Below is some example code.

all<-readr::read_rds("all.rds")

None of these lines seems to work:

all %>% st_intersection
all %>% st_set_precision(1000000) %>% lwgeom::st_make_valid %>% st_intersection
all %>% st_set_precision(1000) %>% lwgeom::st_make_valid %>% st_intersection

1. sf::st_is_valid returns TRUE for all objects. I've tried using mapview, lwgeom::st_make_valid and modifying the precision of the calculations. Only when I set precision to 100 or lower, the "rounded coordinates" become valid geometries to execute the st_intersection function. This issue is related to discussions in #600 and #603. The problem for me is that rounding the coordinates (two decimal places) does not seem appropriate for the work I'm doing. I tried several strategies including rounding (and using st_make_valid) the coordinates just for next rows (sfg):
   all %>% st_set_precision(100) %>% lwgeom::st_make_valid %>% sf::st_is_valid
   However, the issue remains.

2. A subset of the data that generates one of the TopologyException errors is:

subset<-all %>% subset(id %in% c(47,98,158))
st_intersection(subset)

I imagine I am off-track here but It seems tricky to me that the st_intersection of these three multipolygons do not work but pairwise-intersections of these rows do not generate any error. Could you please let me know what I am missing?

int12<-st_intersection(subset[c(1,2),])
int13<-st_intersection(subset[c(1,3),])
int23<-st_intersection(subset[c(2,3),])

Apologies if I've posted this in the wrong place. I hope I am not completely off-track. In case, feel free to delete this.

[dbaston]

If you can output the offending geometries in WKB format, you could submit it as a bug report to GEOS.

[edzer]

Read the documentation section starting with "when called with a missing ‘y’,", and consider calling

out <- st_intersection(all, all)

to avoid the behavior described there (and get pairwise intersections, which you seem to want).

[carlosmolinaguerra]

Thanks for the reply and your time.
@edzer: Apologies if was not clear. I do not want pairwise intersections. Each polygon in input "all" represents a spoken language. But, since several languages may be spoken in a location, the polygons overlap. My ideal output is all non-empty intersections of the geometries and thus, each location would only have one polygon that indicates which language(s) is(are) spoken, which as far as I understand it is made by "st_intersection" with a missing ‘y’ (and the field "origins" indicating the indexes of the languages).

I have been working in this a bit and below what I find. I'll assume just for a moment that the input has 3 features as here:

all<-readr::read_rds("all.rds")
x<-all %>% subset(id %in% c(47,98,158),select=geometry)

In this case, I create all non-empty intersections of these 3 geometries "manually" as:

i123<-st_intersection(x[1,],x[2,]) %>% st_intersection(.,x[3,])
i12<-st_intersection(x[1,],x[2,]) %>% st_difference(.,i123 %>% st_combine %>% st_union)
i13<-st_intersection(x[1,],x[3,]) %>% st_difference(.,i123 %>% st_combine %>% st_union)
i23<-st_intersection(x[2,],x[3,]) %>% st_difference(.,i123 %>% st_combine %>% st_union)
i1<-st_difference(x[1,],x[2:3,] %>% st_combine %>% st_union)
i2<-st_difference(x[2,],x[c(1,3),] %>% st_combine %>% st_union)
i3<-st_difference(x[3,],x[1:2,] %>% st_combine %>% st_union)

So that the final output is:
manual<-rbind(i123,i12,i13,i23,i1,i2,i3)

Which should be equivalent to:
st_intersection(x)

However, note that unlike the "manual" program, the line below throws a TopologyException. I wonder if "st_intersection" should be throwing that TopologyException. Any answer here would be super appreciated.

In additional, all non-empty intersections of the main database "all" were created using the first method ("manual").

[rsbivand]

Please try in rgeos to see whether the precision model affects this, after converting all to sp classes. Further, try lwgeom::st_make_valid(all) to see whether there is something in all that is a topological issue.

[rsbivand]

With:

library(sp)
library(sf)
all_sp <- as(all, "Spatial")
x_sp <- all_sp[all_sp$id %in% c("47", "98", "158"),]
library(rgeos)
xt1 <- gIntersection(x_sp, x_sp, byid=TRUE, drop_lower_td=TRUE)

apart from generating geometry collections (point and line intersections needing dropping as having no area) with no apparent problems. Is id the variable identifying languages? However:

> i1 <- gIntersects(all_sp, byid=TRUE, returnDense=FALSE, checkValidity=TRUE)
> table(sapply(i1, length))

   1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
 299  968 1035 1113 1075  865  615  415  306  201  129   81   91   55   37   18 
  17   18   19   20   21   22   23   24   25   26   27   28   29   30   31   32 
  22   29   28    9   11    9   15    9    8    7    3    6    3    6    2    4 
  34   36   37   38   39   41   42   44   46   47   49   56   57   62   65   73 
   2    2    5    3    1    2    3    1    1    1    1    1    1    1    1    2 
  88   91   92   95  127  141  205  288 
   1    1    1    1    1    1    1    1 

so one of your polygons intersects with 280 others. How precise do you need your output to be?

[rsbivand]

Further, this completes:

t1 <- gIntersection(all_sp, all_sp, byid=TRUE, drop_lower_td=TRUE)

so the underlying issue may be the precision setting.

> sum(sapply(i1, length))
[1] 44004
> length(t1)
[1] 15930

suggest that care may be needed, but we are dropping line and point intersections. I needed about 15GB on a 16GB machine. There are also 7510 self-intersections, I think.

[carlosmolinaguerra]

Thanks @rsbivand (and @dbaston for submitting the report)!

• Yes, "id" is the variable identifying languages.
  -There is a polygon that intersects with 280 others (correct me if I am wrong but I think this should not be a limitation to do the intersection, plus this does not means that 280 polygons relies in the same location). I agree that precision may be import here. In fact, in some cases, two polygons barely intersects (which would create sliver polygons if an intersection is carried out).
  -lwgeom::st_make_valid(all) does not make changes to the features.
  -Thanks for your explanation using gIntersection (if I understand correctly, it also shows that pairwise intersection between features should not throw a TopologyException).

[rsbivand]

The original file is on http://spatial.nhh.no/misc/i794/, together with a WKT version (zip-compressed, otherwise 13.7GB), and a GPKG (difficult0.gpkg). This was read into GRASS with snap=1e-09, and written out including unattached geometries -c as difficult_out.gpkg Then:

library(sf)
diff_all <- st_read("difficult_out.gpkg")
# drop geometries not assgned to GRASS cat
diff_all0 <- diff_all[!is.na(diff_all$cat),]
# create 4-character strings to sort properly
diff_all0$id4 <- formatC(as.integer(diff_all0$id), flag="0", width=4, format="d")
library(units)
lim <- 1
units(lim) <- as_units("km2")
units(lim) <- "m2"
diff_areas <- st_area(diff_all0)
# drop all geometries with  an area of less than 1 km2
diff_all0a <- diff_all0[diff_areas > lim,]
# aggregate GRASS output geometries by language code
diff_all1 <- aggregate(diff_all0a, list(id=diff_all0a$id4), head, n=1) 
# impose GEOS-validity on geometries (about 5 were invalid)
diff_all2 <- st_make_valid(diff_all1)
# carry out the unary intersection
diff_all_int <- st_intersection(diff_all2)
# create the language code look-up table
lut <- unique(sort(diff_all_int$id4))
# convert list column `origins` to a string column
origins_str <- sapply(diff_all_int$origins, function(x) paste(lut[x], collapse=":"))
diff_all_int$origins_str <- origins_str
# drop output areas less than 1 km2
out_areas <- st_area(diff_all_int)
diff_all_int1 <- diff_all_int[out_areas > lim,]
# extract polygons from GEOMETRYCOLLECTION output
gts <- st_geometry_type(diff_all_int1)
oo <- st_geometry(diff_all_int1)
oo_GC <- oo[gts == "GEOMETRYCOLLECTION"]
for (i in seq(along=oo_GC)) { oo_GC[[i]] <- st_sfc(st_union(st_collection_extract(oo_GC[[i]], "POLYGON")))[[1]] }
oo[gts == "GEOMETRYCOLLECTION"] <- oo_GC
# and restore to object
st_geometry(diff_all_int1) <- oo
st_write(diff_all_int1, "diff_all_int1.gpkg", delete_layer=TRUE)

and

diff_all_int1 <- st_read("diff_all_int1.gpkg")
library(mapview)
mapviewOptions(fgb = FALSE)
mapview(diff_all_int1)

It is not clear that the assignement of the language codes from diff_all_int$origins is correct, and one would have to have the actual language strings to know what they may mean.

The input polygons had multiple problems, not just that they overlapped. Comments welcome: @dr-jts @edzer @dbaston

[edzer]

I can run this, but cannot follow 100% of what you are doing. What is question?

[rsbivand]

See https://lists.osgeo.org/pipermail/geos-devel/2020-October/009766.html and subsequent. The issue 794 case is in part used in JTS and GEOS, and I followed up checking whether OverlayNG resolved all of 794. Roger Bivand Falsensvei 32 5063 Bergen tir. 20. okt. 2020, 18:25 skrev Edzer Pebesma <notifications@github.com>:

…

I can run this, but cannot follow 100% of what you are doing. What is question? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#794 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACNZ3BCD5A3NKLFAD5CTU43SLW2Y5ANCNFSM4FJZICDQ> .

[dr-jts]

The original file is on http://spatial.nhh.no/misc/i794/, together with a WKT version (zip-compressed, otherwise 13.7GB), and a GPKG (difficult0.gpkg).

Note that it may be the case that the WKT version of the data has been pointwise precision truncation. This means that it won't reproduce this issue correctly for two reasons:

• several of the polygons are now invalid, due to precision truncation causing vertices "jumping" across nearby line segments
• The change in precision tends to eliminate robustness occurs that might otherwise occur

Testing should be done with a full-precision version of the data.

[dr-jts]

I can run this, but cannot follow 100% of what you are doing. What is question?

@edzer the goal is to determine if the original TopologyException errors are eliminated by the use of OverlayNG. This appears to be the case when using JTS and GEOS directly (for example, this GEOS unit test is an extract of two of the polygons from the data in this issue).

My understanding is that @rsbivand reran this case using GEOS OverlayNG and still encountered TopoEx errors. If this is correct I would like to find out why that is. To do that I will need a WKB dump of the two geometries causing the error in a GEOS overlay function. (And if there are multiple errors, ideally a dump of all the cases). This might need to be done by adding some debug tracing to the rsf C code?