I am working with Postgres 9.6 and PostGIS 2.2.2. These are my tables:

                                     Table "public.allparcels"
  Column  |            Type             |                        Modifiers
 gid      | integer                     | not null default nextval('allparcels_gid_seq'::regclass)
 the_geog | geometry(MultiPolygon,4326) |
 county   | character varying           |

                                       Table "public.parcels"
    Column    |          Type           |                            Modifiers
 ogc_fid      | integer                 | not null default nextval('parcels_ogc_fid_seq'::regclass)
 wkb_geometry | geometry(Polygon,4326)  |
 county       | character varying       |

When I run this (via psycopg2):

INSERT INTO allparcels(county, the_geog) 
SELECT 'Leeds', ST_MakeValid(ST_Union(wkb_geometry)) FROM parcels WHERE county='Leeds'

I get this error:

Traceback (most recent call last):
  File "combine.py", line 21, in <module>
    cursor.execute(q, (county, county))
psycopg2.InternalError: GEOSUnaryUnion: TopologyException: Input geom 1 is invalid: Self-intersection at or near point -1.394465365015132 53.741246412906918 at -1.394465365015132 53.741246412906918

I thought I must have invalid underlying geometries, but I can't see any:

SELECT * from parcels WHERE NOT ST_isvalid(wkb_geometry);

produces no results. As the query shows, I'm already running ST_MakeValid to make sure the unioned polygon is valid before inserting it.

What am I doing wrong?


This is the latest query I have tried. Debugging the query indicates that it's the ST_Union that fails. Decreasing the precision (even to 0.01) does not help:

INSERT INTO allparcels(county, the_geog) 
SELECT 'Leeds', ST_MakeValid(ST_Union(ST_SnapToGrid(wkb_geometry, .000001))) 
FROM parcels WHERE county='Leeds'
  • If you're not mapping CPU die masks in WGS84 your could try using ST_SnapToGrid() to impose a reasonable scale factor to the coordinates (0.0000001 is ~1cm). You might find more errors that way.
    – Vince
    Feb 22 '17 at 1:28
  • @Vince thanks. 1m accuracy should be fine. Just so I understand the syntax, would you do ST_Union(ST_SnapToGrid(wkb_geometry, .00001))?
    – Richard
    Feb 22 '17 at 8:44
  • 1
    Could you possibly upload this table somewhere, or enough of a subset to produce the error. I am reasonably confident that the answer below will work, but if it doesn't, there is another possible approach, but it would be easier to experiment locally, than speculate. Mar 1 '17 at 9:37
  • 2
    @RoVo. ST_Buffer(geom, 0) works for fixing self-intersections in input geometries (often coming from shp files), but I don't think it will help you when the input geometries are all reported as valid, and the self-intersection is coming from the intersection itself. Mar 1 '17 at 9:52
  • 1
    Richard, I'm very keen to help get to the bottom of this. If the below doesn't work, I'd be very happy to work on your data this weekend, if you could upload somewhere. Mar 3 '17 at 9:22

This happens often with ST_Intersection, irrespective of whether you use ST_SnapToGrid (which is more useful for ensuring a certain precision than for fixing geometry errors) and ST_MakeValid. The problem is to do with the fact that when you intersect polygons, often they will meet at a single point or along a line, as well as producing a (Multi)Polygon intersection(s). That is, when you intersect any two polygons, they can have multiple intersections, not all of which will be polygonal. Consequently, the data type for that particular intersection will be a GeometryCollection. As you are then attempting to insert this into a MultiPolygon or Polygon field, you will see errors about non-noded intersection or self-intersections. The simplest way I have found to fix this issue is using ST_CollectionExtract which allows you to extract only Points, Linestrings or Polygon types from the GeometryCollection. In your case, using the parameter 3 (for polygons) and dropping ST_MakeValid ought to fix it:

INSERT INTO allparcels(county, the_geog) 
SELECT 'Leeds', ST_CollectionExtract(ST_Union(wkb_geometry), 3) 
 FROM parcels 
WHERE county='Leeds';

The error you are seeing is consistent with your polygons all being valid in the first place, as you state. I generally think it is better to run ST_IsValid to check and in the case of errors, ST_MakeValid, on geometries before doing any intersections -- as you have done. ST_MakeValid and ST_SnapToGrid have other uses, but are generally not the right tool for fixing geometry collection issues within a query.

EDIT. After a lot of (unsuccessful) fiddling in an attempt to narrow down a specific subset of the polygons that cause this error, we had an insight from the comments. The input data, from the UK Ordnance Survey, originally in CRS 27700 had been reprojected to 4326. As the algorithm for this conversion is based on convergence, rather than being one step, it introduces arbitrary rounding errors into the reprojected geometries. From reading various GEOS mailing lists about the cause of the OP's error, relating to precision and rouding, we came to this realization. The OP has since re-run this job using the original 27700 data with no error.

  • as far as I understood, the problem happens with the use of ST_Union. If you put whatever around ST_Union() that cannot change anything ?!
    – pLumo
    Mar 1 '17 at 9:01
  • @RoVo. Have you tried it? I have seen this error (and related GEOS intersection errors) a million times and this fixed it in my case. Admittedly, more for the non-noded intersection one, but I believe the cause here to be the same. Mar 1 '17 at 9:29
  • 1
    HI Richard, you say that the source of the data is the OS, so my guess is that the data have been reprojected. I tried projecting back to 27700 and using ST_SnapToGrid, which worked for the gid > 140000 and gid < 170000 range, but failed on the whole table. I suspect the original 27700 to 4326 reprojection may have introduced some spurious precision, which is causing the various non-noded intersection type errors. Do you by any chance have the original data in 27700, just to rule this out? Mar 20 '17 at 14:25
  • 1
    Great spot - yes, I just tried re-importing the original data to PostGIS without reprojecting first and it worked like a charm (and very quickly, as well). THANK YOU SO MUCH!
    – Richard
    Mar 22 '17 at 1:28
  • 1
    @Richard. Great. It was funny, OS data reprojected to 4326 makes me feel a bit nauseous, as the buildings look slightly "wrong", like you are sliding off the page somehow. It has all been a bit of a wild goose chase, and my original answer is totally wrong in this context, but an interesting insight. Mar 22 '17 at 18:55

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.