Obtaining Each Unique Area of Overlapping Polygons in Postgres 9.6/Postgis 2.3

Question

I have recently finished compiling Postgres 9.6/GDAL 2.2 and Postgis 2.3 and am looking to build a generic, efficient solution for obtaining each unique area of polygon overlap within a given dataset using a feature balanced grid where said dataset can have many overlapping polygons -some of which will have holes. I have seen a great many guides that attempt this but none that actually work on the kinds of rugged geometry that one encounters in the wilderness of public data.

However, the process of intersecting the data with each grid cell and the polygonizing of the linework seems to result in a dumpster fire of a result set where portions of the grid remain embedded in the final geometry and polygons are inexplicably missing. The following sections will describe the steps I took in full detail (with images).The entire process takes about 30 seconds to run -which would not be terrible if the results were actually correct. I do hope to take advantage of paralellization once I actually figure out how to get the process working at a baseline.

Postgis Full Version

"POSTGIS="2.3.0dev r15032" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 2.2.0dev, released 2016/99/99" LIBXML="2.9.1" LIBJSON="0.11.99" TOPOLOGY RASTER"

Dataset

The Austin Community Registry (you can view and download here) shows boundaries for various community organizations in Austin, TX; has a great deal of overlapping polygons and polygons with holes. At ~ 1.2 MB it is not gigantic enough to spend hours of life waiting for processing to complete, but also large and complex enough to make sure the algorithm functions correctly.

Step 1: Import and Transform

Shapefile is imported using the awesome, life-changing ogr_fdw. st_curvetoline() and st_multi() are used to ensure that every geometry is a proper multipolygon that can function in Postgis and st_setsrid() is employed to set the correct SRID -since the underlying OGR/GDAL driver fails abysmally at reading from .prj files and cannot be trusted under any circumstances.

DROP SERVER IF EXISTS shpfile CASCADE;
CREATE SERVER shpfile
  FOREIGN DATA WRAPPER ogr_fdw
  OPTIONS (
    datasource '/tmp/community_registry',
    format 'ESRI Shapefile'
  );
  DROP SCHEMA IF EXISTS shpschema CASCADE;
  CREATE SCHEMA shpschema;

IMPORT FOREIGN SCHEMA community_registry FROM SERVER shpfile INTO shpschema;

DROP SCHEMA IF EXISTS austin CASCADE;
CREATE SCHEMA austin;
CREATE TABLE austin.community_registry(
  id uuid PRIMARY KEY DEFAULT uuid_generate_v4(),
  identifier text,
  geom geometry(MultiPolygon,4326)
);

INSERT INTO austin.community_registry(
  identifier,
  geom
)


SELECT
  community1 identifier,
  st_transform(st_curvetoline(st_multi(st_setsrid(geom,2277))),4326) geom

FROM

  shpschema.community_registry;

CREATE INDEX ON austin.community_registry(identifier);
CREATE INDEX ON austin.community_registry USING GIST (geom);

Step 2: Subdivide Geometry and Dump Inner/Outer Rings

In order to reduce number of veriticies of large polygons, I subdivide the geometry for efficient geoprocessing. Both inner and outer rings are dumped as one must account for holes. Finally, ring geometries are cleaned and made valid and then dumped again -to deal with cases of multipolygons or geometry collections.

DROP TABLE IF EXISTS grid;
CREATE TABLE grid(
  id bigserial PRIMARY KEY,
  isect_count bigint,
  geom geometry(Polygon,4326)
);

DROP TABLE IF EXISTS subdivided_rings;
CREATE TABLE subdivided_rings(
  id bigserial PRIMARY KEY,
  geom geometry(Polygon,4326)
);

INSERT INTO subdivided_rings(
  geom
)

SELECT

  (st_dump(st_makevalid(st_subdivide(
    (st_dumprings(
      (st_dump(geom)).geom
    )).geom
  )))).geom

FROM

(
SELECT DISTINCT ON (geom)

  geom
FROM

  austin.community_registry

) subdivided_rings;

CREATE INDEX ON subdivided_rings USING GIST (geom);

Step 3: Create Feature Balanced Grid

Because the linework of all polygons is going to be sent to st_polygonize() it is necessary to create a grid which can be imposed over the dataset in order to limit the amount of linework being fed to the operation at one time and improve performance. To do this, I create a feature balanced grid (roughly following the same approach take here) by generating a recursive CTE that splits the extent into 4 smaller pieces and continues splitting until the number of features that intersect each grid cell is less than or equal to 50. The filter (where not st_contains(box2d(dataset.geom),box2d(gridsplit.geom))) is necessary because otherwise the function may never complete due to the high degree of overlapping polygons. The WHERE st_area(gridsplit.geom) > 0 is necessary because otherwise linestrings can be generated(!) instead of polygons due to what I imagine is due to a precision issue in GEOS. In any event, there do not appear to be any issues with the resulting grid as all cells can be st_union()ed back into the extent geometry.

WITH RECURSIVE grid as (
SELECT
  st_setsrid(st_extent(geom),4326) geom,
  count(*) isect_count
FROM
  subdivided_rings

UNION ALL

SELECT
  gridsplit.geom geom,
  count(*) filter (where not st_contains(box2d(dataset.geom),box2d(gridsplit.geom)))

FROM

(
SELECT
  st_setsrid(unnest(array[
    st_makebox2d(st_makepoint(st_xmin(geom),(st_ymin(geom) + st_ymax(geom))/2.0),st_makepoint((st_xmax(geom) + st_xmin(geom))/2.0,st_ymax(geom))),
    st_makebox2d(st_makepoint(st_xmin(geom),st_ymin(geom)),st_makepoint((st_xmax(geom) + st_xmin(geom))/2.0,(st_ymax(geom) + st_ymin(geom))/2.0)),
    st_makebox2d(st_makepoint((st_xmax(geom) + st_xmin(geom))/2.0,(st_ymax(geom) + st_ymin(geom))/2.0),st_makepoint(st_xmax(geom),st_ymax(geom))),
    st_makebox2d(st_makepoint((st_xmax(geom) + st_xmin(geom))/2.0,st_ymin(geom)),st_makepoint(st_xmax(geom),(st_ymax(geom) + st_ymin(geom))/2.0))
  ]),4326)::geometry geom

FROM

  grid

WHERE

  grid.isect_count > 50

) gridsplit
JOIN subdivided_rings dataset ON (gridsplit.geom && dataset.geom)

WHERE

  st_area(gridsplit.geom) > 0

GROUP BY

  gridsplit.geom
)

INSERT INTO grid(
  isect_count,
  geom
)

SELECT
  isect_count,
  geom

FROM

  grid

WHERE

  isect_count <= 50;

CREATE INDEX ON grid USING GIST (geom);

Step 4: Intersect Subdivided Ring Polygons With Grid Cells and Polygonize Each Cell

This is where things seem to go way off the rails -like Trump after Hillary received a larger than expected convention bounce. Basically I use st_intersection() to clip (tried using ST_ClipByBox2D() but this function breaks the result even further) each geometry to its respective cell. I then use st_makevalid(),st_buffer(geom,0) and st_collectionextract(geom,3) to ensure as much as possible that the results of the intersection are not complete garbage useless garbage. I then group the clipped features by grid cell id and feed the exterior rings to st_polygonize() in order to create each discrete area. Polygons are then dumped and collection extracted. As you can see from the image below, several polygons are missing in the resultant dataset, and holes exist where they should not.

DROP TABLE IF EXISTS raw_portions;
CREATE TABLE raw_portions(
  id bigserial PRIMARY KEY,
  geom geometry(Polygon,4326)
);

INSERT INTO raw_portions(
  geom
)

SELECT
  st_collectionextract((st_dump(st_makevalid(st_polygonize(st_makevalid(st_exteriorring(geom)))))).geom,3) geom

FROM

(
SELECT
  grid.id,
  (st_dump(
    st_collectionextract(
      st_buffer(st_makevalid(st_intersection(subdivided_rings.geom,grid.geom)),0)
      ,3
      )
    )
  ).geom geom
FROM
  grid
  JOIN subdivided_rings ON (st_intersects(grid.geom,subdivided_rings.geom))
) clipped_ring_polys

WHERE

  st_geometrytype(geom) = 'ST_Polygon'

GROUP BY

  id;

CREATE INDEX ON raw_portions USING GIST (geom);
--CREATE INDEX ON raw_portions USING GIST (st_pointonsurface(geom));

Step 5: Intersect Discrete Area Polygons Surface Points with Original Polygons and Use st_union() to Erase Grid.

This step basically just takes raw distinct area and intersects their surface points using st_pointonsurface(). Polygons that are actually holes are simply ignored, and do not become a part of the final dataset. An array of feature ids is generated by grouping on the raw_portions id and this array is then used to union together all features that are part of the same overlap in order to erase the grid. However, like in the previous step, this fails to produce the expected result and leaves a significant amount of gridlines still intact as you can see below.

DROP TABLE IF EXISTS reg_pieces;
CREATE TABLE reg_pieces(
  id bigserial PRIMARY KEY,
  id_array uuid[],
  geom geometry(MultiPolygon,4326)
);

INSERT INTO reg_pieces(
  id_array,
  geom
)


SELECT
  id_array,
  st_multi(st_collectionextract(st_union(geom),3)) geom

FROM

(
SELECT DISTINCT ON (raw_portions.id,id_array)
  raw_portions.id id,
  array(SELECT DISTINCT ON (ids.id) ids.id FROM unnest(array_agg(orig_dataset.id)) ids(id))  id_array,
  raw_portions.geom

FROM

  raw_portions
  JOIN austin.community_registry orig_dataset ON (raw_portions.geom && orig_dataset.geom and st_intersects(st_pointonsurface(raw_portions.geom),orig_dataset.geom))

GROUP BY

  raw_portions.id
) valid_parts

GROUP BY

  id_array;

CREATE INDEX ON reg_pieces USING GIST (geom);

Update: Attempt Using PostGIS Topology

I took a stab using Postgis Topology by creating an empty Topology and then just adding each polygon from the original dataset to it in order to auto-polygonize the data. However after attempting various tolerances and trying without any tolerance, I was unable to get the process to complete without throwing an error:

ERROR:  Side-location conflict: new edge starts in face 5857 and ends in face 7476

It's catch 22 that you can't add anything to the topology unless it works correctly, when you are adding the data for the explicit purposes of FIXING the invalid topological relationships! Postgis is really unforgiving of imperfection in geometries.

SELECT CreateTopology('comm_reg_topo',find_srid('austin','community_registry','geom'));

SELECT

  TopoGeo_AddPolygon('comm_reg_topo',(st_dump(geom)).geom)

FROM

  austin.community_registry;

Answer 1

(This is more of a comment, since I don't have time to fully reproduce this, but too long for the comment box.)

This is a challenging problem, because ST_Polygonize requires perfectly noded inputs. Since you're processing your data by grid cell, ST_Polygonize can succeed only if the results of overlay operations in adjacent grid cells produce identical floating point numbers for generated coordinates. If you copy-paste the WKB for some of your inputs to ST_Polygonize into JTS TestBuilder, you'll no doubt see that this is not the case. Wrapping the inputs to ST_Polygonize with a ST_UnaryUnion or ST_Node could help, but only some of the time, and will decrease performance. Without support for GEOS precision models in PostGIS (and possibly even with them), I would not be optimistic that this limitation can be overcome.

Beyond that, I'd want to look into your usage of ST_Buffer(geom, 0) and ST_MakeValid. These functions will gladly discard parts of their inputs in order to make the result valid, which may explain some of your missing polygons.

If you're open to different approaches, I'll mention that I've had excellent results for problems like this using CGAL's 2D Snap Rounding package. This functionality cannot be accessed through PostGIS at the moment (and I'm not aware of any plans to, though I think it would be a fantastic addition), so you'd have to write your own (short) C++ program.