Unexpected TopologyException during ST_Union aggregate function with valid polygons

[blaggner]

When aggregating geometries of a large table (53.5 Mio. polygons) with ST_Union grouped by an ID column (resulting rows at the end: 2958073), ST_Union threw an error: ERROR: GEOSUnaryUnion: TopologyException: found non-noded intersection between LINESTRING (3.3987e+06 5.33572e+06, 3.3987e+06 5.33572e+06) and LINESTRING (3.3987e+06 5.33572e+06, 3.3987e+06 5.33572e+06) at 3398704.354833886 5335719.4814453134 I tested the geometries of the table with ST_IsValid and there were no invalid geometries. With the help of a wrapper function, I was able to identify the rows that were causing the problem (see attached data). However, when trying to reproduce the error on the selected geometries, ST_Union was able to aggregate the geometries without a problem.

[blaggner]

CSV-output of geometries possibly causing an TopologyException with aggregate function ST_Union

[robe]

This might be a geos issue. Please provide output of

SELECT postgis_full_version(), version();

[blaggner]

Replying to robe:

This might be a geos issue. Please provide output of

SELECT postgis_full_version(), version();

Output is:

"POSTGIS="2.3.3 r15473" GEOS="3.5.1-CAPI-1.9.1 r4246" SFCGAL="1.2.2" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 2015/09/16" LIBXML="2.9.3" LIBJSON="0.11.99" TOPOLOGY RASTER";"PostgreSQL 9.5.14 on x86_64-pc-linux-gnu, compiled by gcc (Ubuntu 5.4.0-6ubuntu1~16.04.10) 5.4.0 20160609, 64-bit"

[pramsey]

It's cool you could find rows that, when removed, caused the problem to stop... it's not surprising that they on their own do not cause the problem, because the problem is an interaction effect, probably part-way up the tree that union is traversing as it unions objects, and then unions those unions, etc, etc. It *might* be possible to make a small set that replicates the error, by adding into your extraction the features in the general area of the problematic features you found. Maybe. It's also possible the error happens far enough up the tree that a much larger chunk of data is needed to replicate it. I'm moving this over to PostGIS GEOS milestone and hoping that with GEOS 3.8 we'll have robust overlay and all these go away.

[blaggner]

This was a great idea, Paul, I added all polygons intersecting with the bounding box of my previous extraction (see new attachment). This time, ST_Union ran into the topology exception again, so this might be an usable test case.

[blaggner]

CSV-output of valid geometries causing an TopologyException with aggregate function ST_Union (new test data, including previous test data plus polygons intersecting with their bounding box)

[mdavis]

Thanks for the test case. It does indeed trigger the error, in GEOS and in JTS.

The reason for the error is that there are a LOT of "slivers" in the dataset - i.e. polygons whose boundaries almost line up, but not quite (the vertices are very slightly different to the adjacent polygon).

The good news is that this should be handled by the forthcoming overlay code. The bad news is that the slivers will very much show up in the output (typically as tiny holes in the middle of a large unioin polygon. This may or may not be an issue depending on your intended use of the data.

[blaggner]

That is interesting to know! Even though the tiny holes will surely be unwanted for further processing, I immediately have ideas how to possibly deal with them as long as the ST_Union is able to run without an error. But, weren't there improvements in planning that would enable ST_Union and other functions to run with a user-defined reduced accuracy? That would surely solve (or reduce) the sliver problems, don't you think?

[fengel]

I have similar TopologyExceptions on valid geometries quite often and I am looking forward to the approvements that are on the way in JTS and GEOS. In my experience large polygon tables in the real world (data from institutions, geo portals etc.) are very often subject to this kind of problem. Since the exception results in aborting the process it can become quite problematic if you have to come up with a solution. I was wondering if there are any recommendations for workarounds to this issue. I currently use the functions st_safe_difference and st_safe_intersection from this extension at the moment: ​https://pgxn.org/dist/lostgis/1.0.2/ plus a self created safe_union function. The main mechanism is exception handling which allows you to ignore exceptions that cannot be solved. Using these functions I was able to work with any kind of data so far and I didn't notice any major drawbacks. But I am not sure if these solutions are reasonable and it would be very interesting to see workaround strategies from other people.

[mdavis]

Of course one way to deal with the "tiny holes" is to simply remove holes of less than a given area. However, it can happen (although not necessarily in your dataset) that the spurious holes are actually long slivers, which can have significant area but a very small cross-section. A while ago I did some work on code to identify this sort of geometry, and it worked pretty well. (One way is to buffer the sliver ring with a small negative distance and see if it disappears)

There is also the potential for gores (inward thin spikes). These require a different approach to deal with.

Which brings us to the question about whether the explicit precision model in the new overlay code can help with this. Unfortunately, the answer is "not entirely". Many times it does cause slivers to collapse away. But there can be cases where the sliver is actually widened by the grid snapping, in which case of course it still shows up as a hole.

The technical reason for this is that the snap-rounding is space-driven, rather than data-driven. So if two nearby vertices happen to lie in adjacent grid cells rather than the same one, they are snapped apart rather than together.

(It would be a very interesting project to try and implement a data-driven snapping facility. This would be relatively easy to do up to the 95% level, but getting to 100% in a safe way is a challenge. Other systems such as Arc do provide something like this, so it is possible...)

[mdavis]

That LostGIS API is quite interesting - thanks for posting it.

At first glance it seems like that's probably about the best that can be done to make the current PostGIS overlay functions error-free. I'm guessing that there must be situations where they don't produce a "reasonable" result however? Have you observed anything like this?

[mdavis]

I will also note that the LostGIS approach doesn't help with aggregate ST_Union as it stands. But it might be worth looking at their approach to see if it can be incorporated. Although the new overlay code should solve this problem anyway, so perhaps that will arrive first.

[mdavis]

Looking at the LostGIS overlay functions in more detail, it appears that they drop back to using MakeValid, and failing that a very small buffer. This is certainly one way to deal with issues, and it's nice to have it wrapped up in a convenient function.

Actually using buffer( [A,B], 0 ) is an effective substitute for union(A, B). Do you use that inside your safe_union?

And that makes me think it might be a useful fix to add to ST_Union aggregate, to handle the failures which occur there.... Could be a quick short-term fix to help older PG versions.

[fengel]

I attached my safe_union. It's probably pretty far from optimal, because I am more a user than a developer. I was glad that my first custom aggregate function worked at all. Maybe you can give me some hints to make it better. I am not sure about the best methods to attemt to fix exceptions. In LostGIS the main approaches seem to be translating, buffering by a small value (but not buffering back with a negative value?) and force_rhr. In my function I tried snapping to different precisions. I have the impression that at least the snapping doesn't succeed to fix things very often. The main benefit for me is that exceptions are skipped. I plan to write the point coordinates and geometries of the unfixed exceptions to a (temporary) table so that I can check visually in my GIS-software where exeptions occured and if the output has significant gaps compared to the original polygons. A result with minor gaps would be a "reasonable result". But of course it is a subjective decision whether a result is still acceptable or not. As you said, it also depends on the intended use. So far using LostGIS and my functions the results are ok on that account. I often check total area before and after and decide if I can live with the result. I usually use the result for aggregating area sums by factors. So if the result is slighly wrong because there is 1% of the original area missing I can usually live with that.

Also I want to point out that I am probably not using the latest versions of PostGIS and GEOS. I should probably upgrade but my operating sytem and Postgres needs upgrade, too, and my time is limited at the moment.

select postgis_full_version()

POSTGIS="2.4.2 r16113" PGSQL="95" GEOS="3.5.1-CAPI-1.9.1 r4246" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 2.2.2, released 2017/09/15" LIBXML="2.9.3" LIBJSON="0.11.99" LIBPROTOBUF="1.2.1" RASTER

select version()

PostgreSQL 9.5.19 on x86_64-pc-linux-gnu, compiled by gcc (Ubuntu 5.4.0-6ubuntu1~16.04.10) 5.4.0 20160609, 64-bit

[komzpa]

Hey there,

(I'm original author of LostGIS).

I didn't do ST_Safe_Union because I didn't know back then how to build it as an aggregate that wraps aggregate.

The reason the ST_Safe_* buffers geometries and not buffers them back was that there was a separate "deoverlap" build phase, and it was more convenient to get slightly overlapping geometries. Result lives at ​https://eu.wargaming.net/globalmap/.

The problem of large slivers is a real thing. I had a lot of it when was trying to union and subdivide the GADM polygons.

Here's an archive of functions before it got packaged as extension: ​https://github.com/wgnet/globalmap/tree/master/code/postgis_wrappers - ST_FilterSmallRings was used to get rid of slivers.

What I've used for slivers was along the lines of ST_Safe_Repair(ST_Buffer(ST_Safe_Repair(ST_Collect(ST_Safe_Repair(geom))),0)).

A solution to the madness was ​https://github.com/postgis/postgis/pull/228/files - I believe Juno still uses a build with that thing enabled, but it was reverted upstream.