PostGIS: ST_Union fails with valid polygons?

Question

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;

It returns CREATE FUNCTION but when I actually try to use it:

SELECT safe_union(geom_simpler) AS geom from usgovt WHERE gap_sts='2';

I get 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 ST_Union works.

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?

Answer 1

Same issue with valid geometries, the solution for me was to run :

SELECT ST_Union(ST_SnapToGrid(geom, 0.00001))

So in your case :

SELECT ST_Union(ST_SnapToGrid(geom_simpler, 0.00001)) AS geom from usgovt WHERE gap_sts='2';

If it also failed, something else is wrong in topology.