7

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?

UPDATE:

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. – John Powell Mar 1 '17 at 9:37
  • 1
    @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. – John Powell 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. – John Powell Mar 3 '17 at 9:22
10
+50

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. – John Powell Mar 1 '17 at 9:29
  • @JohnBarça Thank you very much for this answer! Sadly however, it gives me another error: ERROR: GEOSUnaryUnion: TopologyException: Input geom 1 is invalid: Self-intersection at or near point -1.394465365015132 53.741246412906918 at -1.394465365015132 53.741246412906918. – Richard Mar 12 '17 at 22:24
  • I've exported the data to a shapefile using pgsql2shp -h localhost -u dbuser -P dbpass dbname "SELECT ogc_fid, wkb_geometry, county FROM parcels where county='Leeds'". It's here: dropbox.com/s/9cy5qqbae0s97bt/… (62MB) - thank you again for offering to debug, much appreciated! – Richard Mar 12 '17 at 23:36
  • 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

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.