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 (
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:
CREATE MATERIALIZED VIEW IF NOT EXISTS public.inverse_coverage AS 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.