ST_Union Topology exception despite of valid geometries

Question

I have an isochrones layer and I want to calculate the population lays on the overlapping regions.

My first step is to create an intersection between all the isochrones. After that, I try to create a union with the intersections.

Unfortunately, ST_Union throws GEOSUnaryUnion: TopologyException.

I've checked that all the geometries used for the union are valid.

Union query:

select st_union(the_geom) from isochrones_intersection;

ERROR: GEOSUnaryUnion: TopologyException: Input geom 0 is invalid: Self-intersection at or near point 2.3739051999999998 48.544292458332933 at 2.3739051999999998 48.544292458332933

Check geometries:

select count(*) from isochrones_intersection where st_isvalid(the_geom)=false;

count
-------
     0
(1 row)

Here a pg_dump of isochrones_intersectionlayer

I've tried several workarounds without success:

1. Reproject the layer to 3857.
2. ST_SnapToGrid to reduce precision.

Something similar appears at this post, but it looks they implemented a workaround which it's not useful for me.

I agree it's not an easy union.

EDIT

Adding versions:

PostgreSQL 10.1 on x86_64-pc-linux-gnu, compiled by gcc (GCC) 4.9.2 20150212 (Red Hat 4.9.2-6), 64-bit


POSTGIS="2.4.4 r16869" PGSQL="100" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.9.1, 04 March 2015" GDAL="GDAL 2.2.2, released 2017/09/15" LIBXML="2.9.4" LIBJSON="0.11" TOPOLOGY RASTER

I've also created a smaller version of the layer with the issue

Any ideas?

Answer 1

This kind of issue is usually caused by numerical precision errors affecting the robustness of the overlay code in GEOS (used by PostGIS).

The latest version of PostGIS (3.1) uses GEOS 3.9, which contains a much-improved overlay engine (OverlayNG). This will likely fix this issue. (It should, since I used your dataset for testing during the development of the OverlayNG code!)

Answer 2

Here is my opinion, for 'iso_intersection_small' not sure what the problem is, but your geometry in some places look like this ... think about how to arrange it, i.e. the problem is most likely the length of coordinate values, for example ...

(4 row)
2.24052434999272 48.633728,
2.2494507 48.633728,
2.2515106 48.6344147,
2.25202558437494 48.633728,
(6 row) 
2.3648071 48.5733032,
2.36618038333414 48.5733032
(7 row)
2.42197043750091 48.6472892624991, 
...etc.

As a result, your data is not identical,

as evidenced by the artifacts in the figure below, resulting from the difference...

Answer 3

I've tried several things and all of them worked:

1. Transform geometries to a local CRS.
2. Clean the geometries in GRASS by creating a topology.
3. Simplify geometries using ST_SimplifyPreserveTopology with a tolerance of 0.001.

However, the way I was trying to accomplish my algorithm was too hard. What I need is to count population in overlapping areas, so what I've done is to just intersect isochrones against the grid and count grids which appear more than once, if this happens it's a grid which lay on an overlapping area.

Another way to solve this kind of problems could be a raster strategy. It will be much faster.