I have a single MULTIPOLYGON from ESRI ArcGIS with 10,000+ sub-polygons in a Shapefile, some having > 1.4 million vertices. I'm importing the data into PostGIS (ogr2ogr or shp2pgsql, same results) so I can query it for intersection/difference with a LINESTRING. A giant polygon can't be indexed and takes hours to query directly as a single row/column, so dividing it (using default max 255 vertices per polygon) -- lateral join is for unique row IDs:

CREATE TABLE region_subdiv AS 
SELECT row_number() over() AS gid, label, segment.geom AS geom 
FROM region,
LATERAL (SELECT ST_Multi(ST_SubDivide(geom), 255) AS geom) AS segment;
CREATE INDEX region_subdiv_geom_idx ON region_subdiv USING GIST(geom);

When I export the polygons to a Shapefile to view, the cutouts (small negative spaces inside larger polygons) end up as filled polygons stacked on the larger polygons instead. Also, some small shapes become large rectangles, probably because of errors.

I want to calculate the ST_Difference of a bounding box to find parts of a LineString through the region which are not covered, but the subdivided region multi-polygons have spurious overlapping polygons and cause topology exceptions (GEOSDifference: TopologyException: Input geom 1 is invalid: Self-intersection at or near point ...). I've tried a number of methods to simplify the incoming data, including ST_GeomFromTWKB(ST_AsTWKB(geom, <xy_precision>)) with a few different precision levels. It did not appear to change the behavior.

I've tried the main data cleanup approaches:

UPDATE region_subdiv SET geom=ST_Multi(ST_MakeValid(geom));
UPDATE region_subdiv SET geom=ST_Buffer(geom, 0.0);
UPDATE region_subdiv SET geom=ST_SimplifyPreserveTopology(geom, 0.0001) WHERE NOT ST_IsValid(geom);

These cleanup attempts result in geometries for which ST_IsValid is true, but when I visualize my output data, there are still all kinds of overlapping polygons, up to 3 layers deep, and calculating intersections results in topology exceptions. I use the following materialized view to make an inverse of coverage to find the distance and times traveled outside the polygon regions:

    SELECT 1 AS gid, region.name AS region_name,
    ST_Difference(ST_MakeEnvelope(124.41,-16.6717, 125.05, -17.1609, 4326), coverage.geom) as geom
    FROM path, LATERAL (SELECT cover.name, ST_Union(geom) AS geom FROM region_subdiv cover) AS coverage;

I'm using the latest PostgreSQL 10.2 and PostGIS 2.4 on macOS High Sierra.

For more context, here is a question. about the path intersection query I'm using.

Here is a related question about the combination part I'm attempting.

  • 1
    is simply dumping the multipolygon (via ST_Dump) not an option? – ThingumaBob Mar 6 '18 at 22:08
  • 2
    You might need to use something like ST_Snap. Geometries coming from ArcGIS are notorious for having many spurious decimal places, which can cause issues with intersections. – John Powell Mar 7 '18 at 7:00
  • 1
    mmh...you possibly need to do a serious rework of those geometries. large multis always tend to be subject to weird behaviour and performance issues I guess. I would dump the large polygon prior to any further steps into a new table. you could even split them then by a fishnet first and do a union over each gridcell. this way, the polygons are easier to handle. with an index and bbox filter you could query that table for splitting your linestring – ThingumaBob Mar 7 '18 at 9:29
  • 1
    just more suggestions: [1] to get rid of unnecessary precision, as @JohnPowellakaBarça mentioned, you could update your geometries with ST_GeomFromTWKB(ST_AsTWKB(geom, <xy_precision>)). [2] drop union into a new table as the linked post suggested? if you need back attributes, you could create a separate point table with centroids or pointsonsurface and join those with the original table. should be faster that poly/poly intersections. [3] don´t do the LATERAL part, add row number later or from a subquery. cross joining that table with itself (or that generate_series) creates huge overhead. – ThingumaBob Mar 13 '18 at 8:59
  • 1
    Fair enough on the private data front. I really think you should try the ST_Snap approach. I have seen the issue you refer to many times and it has nearly always come from shp files with spurious precision, which is why the validation approaches fail. Please reply with @whoever, so that we know you have replied. – John Powell Mar 13 '18 at 14:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.