I have a national dataset of address points (37 million) and a polygon dataset of flood outlines (2 million) of type MultiPolygonZ, some of the polygons are very complex, max ST_NPoints is around 200,000. I am trying to identify using PostGIS (2.18) which address points are in a flood polygon and write these to a new table with address id and flood risk details. I have tried from an address perspective (ST_Within) but then swapped this starting from the flood area perspective (ST_Contains), rationale being that there are large areas with no flood risk at all. Both datasets have been reprojected to 4326 and both tables have a spatial index. My query below has been running for 3 days now and shows no signs of finishing any time soon!
select a.id, f.risk_factor_1, f.risk_factor_2, f.risk_factor_3
into gb.addresses_with_flood_risk
from gb.flood_risk_areas f, gb.addresses a
where ST_Contains(f.the_geom, a.the_geom);
Is there a more optimal way to run this? Also, for long running queries of this type what is the best way of monitoring progress other than looking at resource utilisation and pg_stat_activity?
My original query finished OK albeit for 3 days and I got sidetracked with other work so I never got to dedicate the time to trying out the solution. However I have just re-visited this and working through the recommendations, so far so good. I have used the following:
- Created a 50km grid over the UK using the ST_FishNet solution suggested here
- Set the SRID of the generated grid to British National Grid and built a spatial index on it
- Clipped my flood data (MultiPolygon) using ST_Intersection and ST_Intersects (only gotcha here was I had to use ST_Force_2D on the geom as shape2pgsql added a Z index
- Clipped my point data using the same grid
- Created indexes on the row and col and spatial index on each of the tables
I am ready to run my script now, will iterate over the rows and columns populating results into a new table until I have covered the whole country. But just checked my flood data and some of the very largest polygons appear to have been lost in translation! This is my query:
SELECT g.row, g.col, f.gid, f.objectid, f.prob_4band, ST_Intersection(ST_Force_2D(f.geom), g.geom) AS geom
INTO rofrse.tmp_flood_risk_grid
FROM rofrse.raw_flood_risk f, rofrse.gb_grid g
WHERE (ST_Intersects(ST_Force_2D(f.geom), g.geom));
My original data looks like this:
However post clipping it looks like this:
This is an example of a "missing" polygon: