1

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:

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.

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

sventechie
  • 113
  • 7
  • why are you using a lateral in the first query to subdivide it? why are you not just using st_subdivide on the geometry in a normal select statement? – ziggy Mar 06 '18 at 19:43
  • its not clear what you are trying to do – ziggy Mar 06 '18 at 19:56
  • So, the various cleanup steps didn't work in some way? – John Powell Mar 06 '18 at 20:58
  • I'm trying to divide a multipolygon and have only valid, non-overlapping sub-parts (database rows). – sventechie Mar 06 '18 at 21:30
  • 1
    is simply dumping the multipolygon (via ST_Dump) not an option? – geozelot Mar 06 '18 at 22:08
  • That wasn't the question. You list the 3 steps above to cleanup errors, but fail to mention what happened and how it didn't work. Also, please don't start questions with I have a huge this, enormous that -- it's kind of meaningless in the days of cloud computing and/or parallel queries. For some 1 million rows is big, for others it is a billion, for a few it is trillions. It isn't useful -- the algorithm is more important than the supposed size. – John Powell Mar 06 '18 at 22:19
  • @ThingumaBob Added more detail, some polygons are way too big so ST_Dump is not sufficient. I had tried SELECT ST_GeometryN(region.geom, series.gid)... (SELECT generate_series(1, ST_NumGeometries(region.geom)) AS gid but it was super slow and large polygons were still unmanageable. – sventechie Mar 06 '18 at 23:03
  • 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 07 '18 at 07: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 – geozelot Mar 07 '18 at 09:29
  • @ThingumaBob Guess I could do the fishnet with ST_SnapToGrid, but how would I union over each grid cell -- do you mean just union everything, or join by grid cell location somehow? – sventechie Mar 08 '18 at 03:58
  • Can you post the actual geometries. It is a bit hard to tell what is really happening from the above. – John Powell Mar 11 '18 at 15:59
  • I'm sorry, the data is private, but I would consider ideas for problem solving in either PostGIS or ArcGIS as a solution. I need a way to either decrease precision on import or to flatten the geometries after splitting them. – sventechie Mar 12 '18 at 18:39
  • 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. – geozelot Mar 13 '18 at 08: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
  • @JohnPowellakaBarça Could you elaborate on the ST_Snap() approach? It seems to be snapping to another shape, which I don't have. Were you suggesting snap to grid? – sventechie Mar 14 '18 at 19:40
  • Yes, ST_SnapToGrid, apologies. Given that 1 degree at the equator is approximately 100km (declining with a factor equivalent to the cosine of latitude), so you get down to a meter with the 5th decimal place. I don't know what precision you are working with, but using ST_SnapToGrid with a precision of 7 or 8 will give you cm or mm precision. The nano meters of shp files are probably overdoing it. – John Powell Mar 15 '18 at 07:17

0 Answers0