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;

(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.


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?

  • 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
  • I've just added an small version of the layer with the issue – alasarr Mar 11 '19 at 13:31

Here is my opinion:

I don’t know where you took this data, but your geometry is torn apart, this is how the first line table "iso_intersection_small" should look like, think how to put it in order, i.e. what you think corresponds to the value of the "id" field is actually the tails of the previous geometry ... 0106000020E6100000010000000103000000010000002C00000053FE5A0300FC0240FDB50600D84B4840CDCF68064CFB0240AF8DC23FB64B4840747918FFE3F90240DE4724C09F4B4840671C34057CF8024042F79B3F5C4B48400EC6E3FD13F7024071B1FDBF454B48400169FF03ACF502406A204C40024B4840A812AFFC43F40240041BD7BFEB4A484021E4BCFF8FF302404BB26900CA4A484021E4BCFF8FF30240A8D02F00434A48409BB5CA02DCF202405AA8EB3F214A4840425F7AFB73F1024089624DC00A4A4840BC3088FEBFF002403B3A0900E9494840BC3088FEBFF0024005EFF5FFBB494840350296010CF00240B7C6B13F9A494840DCAB45FAA3EE0240E68013C083494840CF4E61003CED0240DFEF614040494840F0C91EFC1FEB0240620DBCFF344948408A16EAFA4FE80240620DBCFF34494840F78A130434E60240DFEF6140404948402463B5F97FE502409858CFFF61494840F78A130434E60240E68013C08349484004E8F7FD9BE70240B7C6B13F9A4948408A16EAFA4FE8024005EFF5FFBB4948404D2BB2FDD7E70240A094FF7FD24948408A16EAFA4FE802403B3A0900E949484004E8F7FD9BE7024089624DC00A4A484091D7DE0264E302405AA8EB3F214A48400BA9EC05B0E20240A8D02F00434A48400BA9EC05B0E20240DF1B4300704A484091D7DE0264E302409884B0BF914A48405D3E480504E90240041BD7BFEB4A4840CF4E61003CED02406A204C40024B4840567D53FDEFED02402389B9FF234B4840567D53FDEFED02405AD4CCFF504B4840DCAB45FAA3EE0240A8FC10C0724B4840350296010CF002407842AF3F894B4840425F7AFB73F10240149337C0CC4B48409BB5CA02DCF20240E5D8D53FE34B4840A812AFFC43F4024081295EC0264C48400EC6E3FD13F7024023B59ABF534C4840ED4A260230F90240A09740005F4C4840CDCF68064CFB024023B59ABF534C484053FE5A0300FC0240694C2D00324C484053FE5A0300FC0240FDB50600D84B4840


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.

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.