11

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)?

2 Answers 2

10
+50

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
8
  • I have looked at your addons git repo before. Inspirstionsl stuff. Nov 6, 2018 at 9:13
  • Wow great solution. You did a pretty amazing job for creating this addons, as well. Before I click to award this answer, I just one to make sure about one thing that's bugging me. Running the code you provide, the 'polygon 5' (of the second figure of the question) doesn't seem to recognize the overlap with another polygon ('ct_overlap_cat2=1'; 'ct_overlap_cat2=0'). Therefore, this polygon ends up classified as 'cat1' and the 'sum=2', instead of 'cat3' and the 'sum=7'. I'm facing a bit of difficulty to debug this little issue. Could you help me?
    – Matt_Geo
    Nov 6, 2018 at 15:40
  • 1
    The only problem with this solution is that the case statements are hard coded. In principle, it should probably be able to handle an arbitrary number of categories. Nov 6, 2018 at 16:20
  • @Matt_Geo I get 6 polygons in the resulting table. The triangle polygon get splitted in three. One with sum=2, one with sum=7 and one with sum=8 as in your wished figure. Nov 6, 2018 at 21:14
  • 1
    What if you replace ST_Centroid(geom) with ST_Area(geom)? Nov 7, 2018 at 22:15
1

I suppose if you use polygon geometry type instead of MultiPolygon everything will fall into place:

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

INSERT INTO tbl_foo (id, geom, att_category, att_value) VALUES 
    (1, ST_SetSRID('POLYGON ((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('POLYGON ((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('POLYGON ((4 4, 3 8, 4 12, 7 14,10 12, 11 8, 10 4, 4 4))'::geometry,4326) , 'cat2', 5 );

The result is 9 entries that correspond to the different intersection options in the example given by you.

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.