Monday 25 January 2016

postgis - ST_Union fails with TopologyException despite valid polygons and using ST_SnapToGrid


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
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'

Answer



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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...