0

I am having a performance issue with ST_Intersects. I know similar questions have been asked several times already, but I cannot seem to make use of the past answers in my case. I am trying to perform a very straightforward update to an addresses table using St_Intersects. I have 500,000 addresses (points) and I want to capture the ID of the flood zones which each address intersects. There are 5,000 flood zones (polygons).

This is my query:

UPDATE nfm.nfm_addresses a
SET fld_zone_id = f.ogc_fid
FROM nfm.s_fld_haz_ar_36103 f
WHERE ST_Intersects(a.geom, f.geom);

I cannot index the flood zone geometry due to the size of some polygons exceeding index-able limits. I try anyway and let the query run overnight to no avail. Here is the EXPLAIN:

"Update on nfm_addresses a  (cost=0.00..872514959.88 rows=1097483 width=184)"
"  ->  Nested Loop  (cost=0.00..872514959.88 rows=1097483 width=184)"
"        Join Filter: ((a.geom && f.geom) AND _st_intersects(a.geom, f.geom))"
"        ->  Seq Scan on nfm_addresses a  (cost=0.00..14264.08 rows=587308 width=174)"
"        ->  Materialize  (cost=0.00..1818.09 rows=5606 width=42)"
"              ->  Seq Scan on s_fld_haz_ar_36103 f  (cost=0.00..1790.06 rows=5606 width=42)"

I then used St_Subdivide to get the flood zones small enough to be indexed. This resulted in 50,000 flood zone polygons instead of the original 5,000. I am able to index the geometry in this case and attempt the same query. Here is the EXPLAIN:

"Update on nfm_addresses a  (cost=0.00..7579370470.20 rows=9533770 width=184)"
"  ->  Nested Loop  (cost=0.00..7579370470.20 rows=9533770 width=184)"
"        Join Filter: ((a.geom && f.geom) AND _st_intersects(a.geom, f.geom))"
"        ->  Seq Scan on nfm_addresses a  (cost=0.00..14264.08 rows=587308 width=174)"
"        ->  Materialize  (cost=0.00..8570.49 rows=48699 width=42)"
"              ->  Seq Scan on s_fld_haz_ar_36103 f  (cost=0.00..8326.99 rows=48699 width=42)"

No geometry index appears to be used. I've tried rearranging my query in several ways using CTEs and/or subqueries to no avail. What am I missing?

I'm using PostgreSQL 11.3, compiled by Visual C++ build 1914, 64-bit. I'm using POSTGIS="2.5.2 r17328" [EXTENSION] PGSQL="110" GEOS="3.8.0-CAPI-1.13.1 " PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.4.4, released 2020/01/08" LIBXML="2.7.8" LIBJSON="0.12" LIBPROTOBUF="1.2.1" RASTER.

My machine is a top of the line desktop with an i9, 64GB RAM, etc.

7
  • 1
    If your addresses are points I'd use contains instead of intersects
    – Ian Turton
    Aug 22, 2020 at 16:09
  • 2
    Update the table stats (VACUUM ANALYZE <table>).
    – geozelot
    Aug 22, 2020 at 16:32
  • Contour-following polygons are problematic for many reasons. One is the sheer volume of vertices, which slows O(N*M) point-in-poly. Another is their large extent, vice the area used within the extent, which produces false positives in the spatial index.
    – Vince
    Aug 23, 2020 at 12:02
  • @Vince Do you have any recommendations on alternative tools? PostGIS is my go-to spatial toolset, but I think this particular problem exceeds its capabilities...
    – Josh
    Aug 23, 2020 at 15:00
  • 1
    No, don't grid the points, grid the polygons. Gridded polygons generally save two of three orders of magnitude in query runtime. I once had a point-in-poly query that ran 22 minutes. When I reversed the order it took 506 milliseconds.
    – Vince
    Aug 23, 2020 at 19:10

1 Answer 1

2

Having no more information than provided, and assuming all tables have an ogc_fid column, I would do the following:

  1. Intersect nfm.s_fld_haz_ar_36103 with a regular fishnet polygon table of 200-1000 tiles to produce nfm.tmp_fld_haz_grid

  2. Make sure the nfm.nfm_addresses table is indexed USING gist(geom)

  3. Use the gridded table to drive a query against the points:

UPDATE  nfm.nfm_addresses a
SET     fld_zone_id = vt.f_ogc_fid
FROM (
    SELECT  f.ogc_fid as f_ogc_fid, a.ogc_fid as a_ogc_fid
    FROM    nfm.tmp_fld_haz_grid f
    JOIN    nfm.nfm_addresses a ON ST_Intersects(f.geom, a.geom)
    ORDER   BY a.ogc_fid 
) vt
WHERE   a.ogc_fid = vt.a_ogc_fid;

Instead of generating 500K point-in-complex-poly queries, this will generate 50k simple-poly-over-point queries, which should return much more quickly.

Note that the poly-point query may generate multiple point occurrences on grid boundaries, and will not populate NULL fld_zone_id values in points that don't overlap any zone. You can either calculate all values to NULL before this UPDATE or nest again to OUTER JOIN the addresses back to themselves by ID.

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.