3

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?

7
  • Please Edit the question to specify the exact versions of PostgreSQL and PostGIS in use. – Vince Mar 11 '19 at 11:04
  • Done! thank you! – alasarr Mar 11 '19 at 11:58
  • You should try to make it easier to help you. You know that something happens close to point 2.373905, 48.54429 so you could try to clip some data around that point and test if you get the same error. – user30184 Mar 11 '19 at 12:05
  • That would be good for a "one shot" run, but it's not the case, it's dynamic algorithm where intersection layer changes, so the point could change – alasarr Mar 11 '19 at 13:04
  • 1
    try st_makevalid() – ziggy Apr 17 '20 at 17:09
1

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... enter image description here

5
  • Not sure how you were looking at the data but what you posted is 1:1 the geometry of the first feature in the SQL file. Maybe your editor failed at displaying things correctly? I verified it with the less tool. – bugmenot123 Aug 23 '20 at 7:41
  • @bugmenot123, I published the first line of corrected geometry... – Cyril Mikhalchenko Aug 23 '20 at 8:46
  • Yes but that's exactly what is in the file. There is nothing "torn apart" as you say. – bugmenot123 Aug 23 '20 at 9:19
  • Hm, sorry for the noise then. But it definitely wasn't OPs issue :) Cheers! – bugmenot123 Aug 23 '20 at 20:03
  • @bugmenot123, my answer clarified the problem of initial OP geodata? – Cyril Mikhalchenko Aug 29 '20 at 14:10
0

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.

0

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!)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.