I am working with this dataset. I have imported it into a Postgres 9.6 database with shp2pgsql as follows:
shp2pgsql -I -s 5070:5070 us_govt/PADUS1_4Fee.shp usgovt | psql -d mydb
I have made all the polygons valid:
UPDATE usgovt SET geom=ST_Multi(ST_CollectionExtract(ST_MakeValid(geom), 3)) WHERE NOT ST_IsValid(geom);
Then simplified them, and re-validated them:
UPDATE usgovt SET geom_simpler=ST_Multi(ST_SimplifyPreserveTopology(geom,0.5));
To confirm, all the simplified polygons are valid (the following produces zero):
SELECT COUNT(*) FROM usgovt WHERE NOT ST_IsValid(geom_simpler);
However, trying an
ST_Union on a group of the polygons fails:
# SELECT ST_Union(geom_simpler) AS geom from usgovt WHERE gap_sts='2'; ERROR: GEOSUnaryUnion: TopologyException: found non-noded intersection between LINESTRING (-3.36702e+06 4.94515e+06, -3.36699e+06 4.94512e+06) and LINESTRING (-3.36699e+06 4.94512e+06, -3.36702e+06 4.94515e+06) at -3367018.7763304636 4945150.2479274161
Why could this be failing, and how can I debug it?
UPDATE: Tried buffering with
SELECT ST_Union(ST_Buffer(geom_simpler, 0.000001)) AS geom from usgovt WHERE gap_sts='2';, but that fails with:
ERROR: GEOSUnaryUnion: TopologyException: Input geom 0 is invalid: Self-intersection at or near point 3216270.6962009948 -12370.16660008162 at 3216270.6962009948 -12370.16660008162.
UPDATE 2: I've tried creating my own version of
ST_Union with better exception handling, as suggested in this answer to a related question:
CREATE OR REPLACE FUNCTION safe_union(geom_a geometry, geom_b geometry) RETURNS geometry AS $$ BEGIN RETURN ST_Union(geom_a, geom_b); EXCEPTION WHEN OTHERS THEN BEGIN RETURN ST_Union(ST_Buffer(geom_a, 0.00001), ST_Buffer(geom_b, 0.00001)); EXCEPTION WHEN OTHERS THEN RETURN ST_GeomFromText('POLYGON EMPTY'); END; END $$ LANGUAGE 'plpgsql' STABLE STRICT;
CREATE FUNCTION but when I actually try to use it:
SELECT safe_union(geom_simpler) AS geom from usgovt WHERE gap_sts='2';
HINT: No function matches the given name and argument types. You might need to add explicit type casts..
Also, I'm not sure how safe this method will be, if it's unioning tens of thousands of polygons, because I don't understand the internals of how
Is it possible that it could successfully union 5000 polygons, then get to the 5001st, hit an error, and silently replace the union of the previous 5000 with a null polygon, ultimately producing an extremely inaccurate result?