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

Question

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.

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

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 1

you can use ST_INTERSECTION using the source table twice. Make sure to avoid comparing the a polygon to itself. Also avoid comparing a pair of polygon twice.

SELECT ST_INTERSECTION(a.geom, b.geom), 'fair'
FROM mytable a, mytable b
WHERE a.ID < b.ID
AND ST_INTERSECTS(a.geom, b.geom);

Answer 2

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)

In the query above I've replaced "good" 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.