Sunday 18 December 2016

postgresql - How get intersection of polygons in the same table by PostGIS?



I have a single table of over 1000 polygons in postgreSQL (with PostGIS), and I have a field "evaluation_type" that only receives two values, either 'good'(green) or 'bad'(red), as you can see in the picture below.


enter image description here Once some 'bad' and 'good' polygons overlap, I'd like to use a PostGIS function where I can get the following:


3 feat


And in the attribute table I'd have - 'good' - 'bad' - 'fair': For the intersections


I'd like the result layer (the one containing the intersection) in a "2D" layer, in other words, I want no overlap between the intersects. It's because, in the end, the thing I'm looking for is to count, for each small little polygon, how many 'good' and 'bad' polygons it represents, see below the hypothetical attribute table I want to have.


id| evaluation | count_good | count_bad


1 |  good    |   223    |   0 
2 | fair | 108 | 288
3 | bad | 0 | 278
4 | fair | 193 | 132

5 | fair | 183 | 198
6 | ...


  • If the small polygon only represents a intersection of 'good' polygons it will receive an value of 'good' and will have 0 in the field 'count_bad'

  • If the small polygon only represents a intersection of 'bad' polygons it will receive an value of 'bad' and will have 0 in the field 'count_good'

  • If the small polygon represents a intersection of 'bad' and 'good' polygons it will receive an value of 'fair' and will have 0 counts > 0 in both fields 'count_good' and 'count_bad'


As far as I remember in ArcGIS Union tool, it would break the intersections. I realized that ST_Union behaves a bit different. I used the following code and it returned a dissolved polygon.


    CREATE TABLE out_table AS  

SELECT ST_Union(geom) as geom
FROM areas_demo;

Do you have a hint, how could I get such result I described?



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 a self contained example of a problem similar to your one:



WITH geomtable AS (
SELECT 1 id, 1 val, ST_GeomFromText('POLYGON((0 1, 3 2, 3 0, 0 1), (1.5 1.333, 2 1.333, 2 0.666, 1.5 0.666, 1.5 1.333))') geom
UNION ALL
SELECT 2 id, 0 val, ST_GeomFromText('POLYGON((1 1, 3.8 2, 4 0, 1 1))') geom
UNION ALL
SELECT 3 id, 1 val, ST_GeomFromText('POLYGON((2 1, 4.6 2, 5 0, 2 1))') geom
UNION ALL
SELECT 4 id, 0 val, ST_GeomFromText('POLYGON((3 1, 5.4 2, 6 0, 3 1))') geom
UNION ALL
SELECT 5 id, 1 val, ST_GeomFromText('POLYGON((3 1, 5.4 2, 6 0, 3 1))') geom

UNION ALL
SELECT 6 id, 0 val, ST_GeomFromText('POLYGON((1.75 1, 1 2, 2 2, 1.75 1))') geom
), parts AS (
SELECT a.id, a.val val, unnest(ST_SplitAgg(a.geom, b.geom, 0.00001)) geom
FROM geomtable a,
geomtable 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, ST_AsEWKB(a.geom), a.val
)
SELECT avg(val*1.0) val, ST_Union(geom) geom
FROM parts
GROUP BY ST_Centroid(geom)

"good" has been replaced with 1 and "bad" with 0. Each polygon part gets the mean of the overlapping parts so that every polygon having val between 0 and 1 is actually "fair".


Should work with thousands of polygons and when there are more than two overlaps.


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...