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"
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_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_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;