Friday, 10 February 2017

postgis - How to recursively loop through parent polygons intersections to get smallest (child) polygons with no overlaps?


I'm struggling with a problem for a couple of days and realized many people also get stuck when the topic is intersections in PostGIS (v2.5). That's why I decided to ask for a more detailed and generic, common question.


I have the following table:


DROP TABLE IF EXISTS tbl_foo;
CREATE TABLE tbl_foo (
id bigint NOT NULL,
geom public.geometry(MultiPolygon, 4326),
att_category character varying(15),
att_value integer

);
INSERT INTO tbl_foo (id, geom, att_category, att_value) VALUES
(1, ST_SetSRID('MULTIPOLYGON (((0 6, 0 12, 8 9, 0 6)))'::geometry,4326) , 'cat1', 2 );
INSERT INTO tbl_foo (id, geom, att_category, att_value) VALUES
(2, ST_SetSRID('MULTIPOLYGON (((5 0, 5 12, 9 12, 9 0, 5 0)))'::geometry,4326), 'cat1', 1 );
INSERT INTO tbl_foo (id, geom, att_category, att_value) VALUES
(3, ST_SetSRID('MULTIPOLYGON (((4 4, 3 8, 4 12, 7 14,10 12, 11 8, 10 4, 4 4)))'::geometry,4326) , 'cat2', 5 );

It looks like this:


start



I want to get all the child polygons based on the intersection of the parent polygons. For the outcome, it'd be expected:



  • The child polygons with no overlap between them.

  • An column containing the sum of the value of their parent polygons,

  • An column containing the count of parent polygons of one category

  • An column containing the count of another category

  • An column containing the category of the child polygon, based on the following rule: -If ALL the parent polygons are from one class, the child polygon also have this class. Else, the category of the child polygon is a third category.


So it'd looks like it:


output



So, in the end, the output table generated (for this example) will have 7 rows (all the 7, non-overlapping, child polygons), containing columns of category, sum_value, ct_overlap_cat1, ct_overlap_cat2


The following code I started, gives me the individual intersections, comparing one parent with another.


SELECT
(ST_Dump(
ST_SymDifference(a.geom, b.geom)
)).geom
FROM tbl_foo a, tbl_foo b
WHERE a.ID < b.ID AND ST_INTERSECTS(a.geom, b.geom)
UNION ALL
SELECT

ST_Intersection(a.geom, b.geom) as geom
FROM tbl_foo a, tbl_foo b
WHERE a.ID < b.ID AND ST_INTERSECTS(a.geom, b.geom);

How do I recursively loop through the result of this mentioned code, that, independent of the number of overlap polygons I always get its 'smallest' (child) polygons (Fig. 2)?



Answer



Try this:


Download the PostGIS Addons from this link: https://github.com/pedrogit/postgisaddons


Install by running the postgis_addons.sql file to get the ST_SplitAgg() function.


Test by running the postgis_addons_test.sql file.



Here is your query:


WITH  result_table AS (
WITH parts AS (
SELECT a.att_value val,
CASE WHEN a.att_category = 'cat1' THEN 1 ELSE 0 END cat1,
CASE WHEN a.att_category = 'cat2' THEN 1 ELSE 0 END cat2,
unnest(ST_SplitAgg(a.geom, b.geom, 0.00001)) geom
FROM tbl_foo a,
tbl_foo b
WHERE ST_Equals(a.geom, b.geom) OR

ST_Contains(a.geom, b.geom) OR
ST_Contains(b.geom, a.geom) OR
ST_Overlaps(a.geom, b.geom)
GROUP BY a.id, a.att_category , ST_AsEWKB(a.geom), val
)
SELECT CASE WHEN sum(cat2) = 0 THEN 'cat1'
WHEN sum(cat1) = 0 THEN 'cat2'
ELSE 'cat3'
END category,
sum(val*1.0) sum_value,

sum(cat1) ct_overlap_cat1,
sum(cat2) ct_overlap_cat2,
ST_Union(geom) geom
FROM parts
GROUP BY ST_Area(geom)
)
SELECT category, sum_value, ct_overlap_cat1, ct_overlap_cat2,
(ST_Dump(result_table.geom)).geom as geom
FROM result_table

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...