I have two tables that I join in my existing workflow. One of these includes flood events (was a shapefile that has been imported to PostgreSQL), the other one is a grid that holds the amount of precipitation for each day and each grid. It is a grid spawned across the world, each pattern one degree high and one degree wide (64800 in total). In my workflow, I combine the flood event with the summed up precipitation in the affected area. And this is how I do it: At first I calculate the centroid of each pattern of the grid (st_centroid).

Then I check if the centroid is within an event of the flood-polygon, using the st_within function:

    SELECT  cell_id, (id of the grid)
            geometry AS cell_geom,
            gid, (id of the flood event)
            geom AS dfo_geom
    FROM    precipitation_grid,
          ST_Within(precipitation_grid.centroid, flood_event_lib.geom);

CREATE INDEX within_gix ON within USING GIST (dfo_geom);

In the next step, I combine the two libraries, based on the cell_id that I used in the first step.

    SELECT cell_id,
           gid AS dfo_id,
    FROM within, precipitation_lib
    WHERE within.cell_id = precipitation_lib.cell AND
(precipitation_lib.date >= within.started_at AND precipitation_lib.date <= within.finished_at) 
    AND precipitation > 0 

ORDER BY gid, cell_id, date;

In the final two steps, I calculate first the average of the daily precipitation

 CREATE TABLE libjoin_avg AS
   SELECT *,
   AVG(precipitation) OVER (PARTITION BY date ORDER BY dfo_id)
   FROM libjoin

And then I sum it up by adding the averages of each day within the flood event

 CREATE TABLE libjoin_sum AS
   SELECT *, 
   SUM(avg) OVER (PARTITION BY started_at ORDER BY dfo_id)
   FROM libjoin_avg

This workflow is working but I would like to optimise it. My grid is not very precise and, as I am using the centroid of a pattern, I miss pretty many flood events because the centroid of a pattern is not always within the flood event. Another problem is, that once the centroid is within the flood event, the pattern is used for the calculation, no matter how much of its area is actually covered. You can see that in the following image, especially on the upper right side. These patterns are used for calculation because the centroid is within the flood event geometry.

This image shows my different layers, imported in QGIS. The green background is the actual grid and the highlighted blue ones are those with the centroid within the flood event. So they are used for the calculation of the summarised precipitation of the flood event shown in light blue on top. enter image description here

Do you have any ideas of how I could optimise my workflow? I would like to weight the patterns by the amount of area that is covered by the flood event polygon and use that weighted amount for the calculation of the overall precipitation for each flood event but I don't know how to do that.

I am open for other solutions, not based on the centroid of the pattern, as well.

  • Hello! Welcome to GIS.StackExchange. Can you reference explicit table names in your paragraphs? I'm having difficulty following which is which. E.g.: What's the difference between precipitation_lib and precipitation_grid?
    – raphael
    Nov 24, 2015 at 21:28
  • Hello, thank you for the welcoming words @raphael! I am sorry, that was some bad explaining: I got one table with the flood events, I think that one is pretty clear. And then there are two more tables, one is the precipitaion_grid which holds the unique cell_id for all 64800 cells aswell as geometry and centroid of the geometry.
    – marius
    Nov 26, 2015 at 7:44
  • @raphael And then there is precipitaion_lib which holds the cell_id aswell and the actual precipitaion value for each cell_id and every day over a 19 years period. In my first step I collect the cells within the flood geometry of the precipitation_grid table and in the second I select those of my precipitation_lib where the cell_id matches and the date is within the time period of the flood event. I hope this makes it clear to you now.
    – marius
    Nov 26, 2015 at 7:50

2 Answers 2


It sounds like you want to intersect the flood shapefiles with the grid-cells, and then weight the precipitation by the proportion of the flood area that intersects each grid cell.
Step 1: Store the grid area.

ALTER TABLE precipitation_grid
UPDATE TABLE precipitation_grid
SET grid_area = ST_AREA(geometry);

Step 2: Intersect, using ST_Intersects, ST_Intersection, and SELECT INTO

    ST_AREA(ST_INTERSECTION(g.geom, f.geom)) as x_area,
    ... other columns
    INTO grid_flood_xsect
FROM precipitation_grid g, 
    flood_event_lib f
WHERE ST_INTERSECTS(g.geom, f.geom)

And then you can perform your calculations by weighting observations based on the area overlap, e.g.:

SELECT gid, cell_id, x_area / grid_area * some_value
FROM grid_flood_xsect
  • Thank you very much for that idea, I think that could work out! But when I try to run the second query, I get the following error: ERROR: Error performing intersection: TopologyException: Input geom 1 is invalid: Self-intersection at or near point -1.7085519303179475 55.098976701523 at -1.7085519303179475 55.098976701523 I will try to fix it
    – marius
    Nov 26, 2015 at 8:44
  • 1
    Okay, it worked, using ST_MakeValid: ST_AREA(ST_INTERSECTION(ST_MakeValid(g.geom), ST_MakeValid(f.geom))) as x_area made it work! Thank you so much, I can implement this in my workflow now!!
    – marius
    Nov 26, 2015 at 11:33

similar answer. select the polygons which intersect more area.

a.cell_id, d.gid, st_area(st_intersection(a.geom, d.geom)) as x_area,
FROM precipitation_grid a, flood_event_lib d
WHERE st_intersects(a.geom, d.geom)
AND (st_area(st_intersection(a.geom, d.geom))/st_area(a.geom)) > .5

another approach

SELECT a.cell_id, d.gid, DISTINCT ON (d.gid) st_area(st_intersection(a.geom, d.geom))::float8
FROM precipitation_grid a, flood_event_lib d
WHERE st_intersects(a.geom, d.geom)
--check sliver intersection
AND st_area(st_intersection(a.geom, d.geom))::float8 > 0.000000000000000001
ORDER BY a.geoid, st_area(st_intersection(a.geom, d.geom))::float8 DESC

Your Answer

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

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