Calculate total coverage of circles on map with PostGIS

Question

I am trying to calculate the total area covered by c.a. 400 Circles of diameter 15Km. The centers of these circles are stored in a column of type geometry:4326 in my database.

Some circles overlap, and I should not add to the total area the overlapping amount twice.

Until now I tried the follwing approach: 1) Use buffer to get circles: select st_buffer(get_location, 15000) from d_station; 2) Use st_union to join the circles into a single multipoint geometry 3) use st_area to calculate area.

Unfortunately this approach doesn't work, I think because I am not using st_union in the correct way?

Then I found an interesting tutorial here, which explains something very similar to what I need. Unfortunately I am not able to get a useful result using this approach either.

How can I obtain a list of circles and merged circles polygons, and then calculate the sum of their areas?

Update: I am trying Pierre suggestion. Here the steps I performed:

CREATE temporary TABLE pts AS
select ds.id, geo_location as geom from pa.d_station ds

CREATE temporary TABLE circles AS
SELECT id, ST_Buffer(geom, 15000) AS geom FROM pts;

create temporary table circles2 as 
SELECT a.id, pa.ST_DifferenceAgg(a.geogr, b.geogr) geogr
FROM circles a, 
     circles b
WHERE ST_Equals(a.geogr, b.geogr) OR 
   ((ST_Contains(a.geogr, b.geogr) OR 
    ST_Contains(b.geogr, a.geogr) OR 
    ST_Overlaps(a.geogr, b.geogr)) AND 
   (ST_Area(a.geogr) < ST_Area(b.geogr) OR 
    (ST_Area(a.geogr) = ST_Area(b.geogr) AND 
     ST_AsText(a.geogr) < ST_AsText(b.geogr))))
GROUP BY a.id;

select sum(st_area(c2.geogr)) from circles2 c2;

But I get a result which is way too low from the expected. Is my approach wrong?

Answer 1

Install the PostGIS Addons and look at the ST_DifferenceAgg() example to remove overlaps:

SELECT a.id, ST_DifferenceAgg(a.geom, b.geom) geom
FROM overlappingtable a, 
     overlappingtable 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)) AND 
       (ST_Area(a.geom) < ST_Area(b.geom) OR 
        (ST_Area(a.geom) = ST_Area(b.geom) AND 
         ST_AsText(a.geom) < ST_AsText(b.geom))))
GROUP BY a.id;

Then you should be able to easily compute the area.

The method described in the tutorial you mention works for a limited number of polygons.