PostGIS Spatial Relationships query performance

Question

I have a very large dataset contains over 700 million points and a polygon dataset as a buffer zone.

My task is to extract all points inside the buffer zone and create a new table.

Below is my code. I test it with a small point dataset and it works fine.

create table schema1.result as

select point.* from

schema1.site as point, schema2.buffer as poly

Where ST_Intersects(point.geo_loc,poly.wkb_geometry);

Unfortunately, the query lasted for 1 day and showed no signs to finish.

Is there any advice to optimise my code to speed up the query?

Update: The output of Explain

"Nested Loop (cost=0.41..17773703.88 rows=6789472 width=208)"

" -> Seq Scan on buffer poly (cost=0.00..18.50 rows=850 width=32)"

" -> Index Scan using idx_site on site point (cost=0.41..20902.23 rows=799 width=208)"

" Index Cond: (geo_loc && poly.wkb_geometry)"

" Filter: st_intersects(geo_loc, poly.wkb_geometry)"

"JIT:"

" Functions: 6"

" Options: Inlining true, Optimization true, Expressions true, Deforming true"

Answer 1

I think what's going on here is that when you write:

from schema1.site as point, schema2.buffer as poly

PostgreSQL is doing a CROSS JOIN between the two tables. When multiple tables are listed in the FROM clause postgres uses a CROSS JOIN source that results in a table with a number of rows equal to the Cartesian product of the two tables source.

to avoid that you can use WHERE EXISTS as:

SELECT
    *
FROM schema1.site AS point_table
WHERE EXISTS(
    SELECT 1 
    FROM schema2.buffer AS poly_table
    where ST_Intersects(point_table.geom, poly_table.geom)
);

This allows PostgreSQL to calculate only the spatial join and not have to build a temporary table of the full cross join which takes forever.

Answer 2

It might be faster to use the polygons as the driving table, using the point table index to filter down the large number of records. This allows PostGIS to optimize the ST_Intersects spatial predicate by preparing each polygon.

This can be forced using LATERAL:

SELECT pt.* FROM schema2.buffer AS poly
JOIN LATERAL (SELECT * FROM schema1.site) AS pt 
  ON ST_Intersects(poly.wkb_geometry, pt.geo_loc);

If the polygons have a large number of vertices, it can be faster to use ST_Subdivide to fragment them before doing the point query:

WITH poly AS (
  SELECT ST_Subdivide(wkb_geometry) AS geom FROM schema2.buffer 
)
SELECT pt.* FROM poly
JOIN LATERAL (SELECT * FROM schema1.site) AS pt 
  ON ST_Intersects(poly.geom, pt.geo_loc);