Ways to speed up ST_DWithin between millions of segments and thousands of polygons?

Question

I have 3 PostgreSQL tables with the following data:

tab1 => It contains about 5.7 millions of segments (straight lines)

tab2 => It contains 6673 polygons

tab3 => It contains 19620 polygons

Projected Coordinate System (British National Grid)

All 3 tables have a spatial index

I am trying to create a new table with all segments from tab1 that are less than 250 metres away from either tab2 or tab3. In order to do so I am running the following query:

CREATE TABLE public.final_def_sgm AS
SELECT tab1.segment_id, tab1.gb_def_id, tab1.def_id, tab1.geom
FROM public.flood_def_sgm_50m tab1, public.sop_polyg tab2, public.coastal_1000_polyg tab3
WHERE ST_DWithin(tab1.geom, tab2.geom, 250) OR ST_DWithin(tab1.geom, tab3.geom, 250);

I launched this query about 48 hours ago and it is still running. The problem is that I have no idea whether it will be running for two more days or for twenty more days (who knows).

I previously tried to do the same thing by creating 250m buffers around the polygons and tried to select those segments that intersected with the polygons (all in postgres). Very slow too.

Is it normal that it´s taking days to run the query? Is there any better way to do what I am trying to do?

"explain" throws the following information:

Gather  (cost=3884.82..447434055344.37 rows=824566434 width=60)

  Workers Planned: 

  ->  Nested Loop  (cost=2884.82..447351597700.97 rows=485039079 width=60)

        ->  Nested Loop  (cost=0.00..6414604.08 rows=127287635 width=6444)

              ->  Parallel Seq Scan on sop_polyg tab2  (cost=0.00..1864.88 rows=6488 width=1435)

              ->  Seq Scan on coastal_1000_polyg tab3  (cost=0.00..792.20 rows=19620 width=5009)

        ->  Bitmap Heap Scan on flood_def_sgm_50m tab1  (cost=2884.82..3514.43 rows=1 width=60)

              Recheck Cond: ((geom && st_expand(tab2.geom, '250'::double precision)) OR (geom && st_expand(tab3.geom, '250'::double precision)))

              Filter: (((geom && st_expand(tab2.geom, '250'::double precision)) AND (tab2.geom && st_expand(geom, '250'::double precision)) AND _st_dwithin(geom, tab2.geom, '250'::double precision)) OR ((geom && st_expand(tab3.geom, '250'::double precision)) AND (tab3.geom && st_expand(geom, '250'::double precision)) AND _st_dwithin(geom, tab3.geom, '250'::double precision)))

              ->  BitmapOr  (cost=2884.82..2884.82 rows=1143 width=0)

                    ->  Bitmap Index Scan on flood_def_sgm_50m_geom_idx  (cost=0.00..14.70 rows=572 width=0)

                          Index Cond: (geom && st_expand(tab2.geom, '250'::double precision))

                    ->  Bitmap Index Scan on flood_def_sgm_50m_geom_idx  (cost=0.00..12.19 rows=572 width=0)

                          Index Cond: (geom && st_expand(tab3.geom, '250'::double precision))

---

I have also tried to obtain the same information using a JOIN in the following way:

CREATE TABLE final_def_sgm AS
SELECT t1.segment_id , t2.coastal_1000_id , t3.sop_id
FROM public.flood_def_sgm_50m t1
LEFT JOIN public.coastal_1000_polyg t2
ON ST_DWithin(t1.geom, t2.geom, 250)
LEFT JOIN public.sop_polyg t3
ON ST_DWithin(t1.geom, t3.geom, 250)

Same thing. Query has been running for more than 48 hours now.

Answer 1

this is my attempt without having the data and under the assumption that you do not care about any attributes at all from tables 2 or 3.

1. merge the geometries of t2 and t3
2. create cluster geometries where there are overlapping geometries to reduce complexity
3. create an index on this layer
4. run a basic st_dwithin on the lines to clustered polygons made of table 2 and 3

code

drop table if exists polys;
create table polys as
with a 
  as(select t2.geom geom from sop_polyg t2
    union
     select t3.geom geom from coastal_1000_polyg t3
    ),
clusters
  as(select st_clusterdbscan(geom,0,2) over() cluster_id,geom from a
  )
select st_union(geom) geom from clusters
where cluster_id is not null
group by cluster_id
union
select geom from clusters where cluster_id is null;
CREATE INDEX i ON polys USING gist (geom);

drop table if exists final_def_sgm;
CREATE TABLE public.final_def_sgm AS
SELECT tab1.segment_id, tab1.gb_def_id, tab1.def_id, tab1.geom
FROM public.flood_def_sgm_50m tab1 join polys on ST_DWithin(tab1.geom, polys.geom, 250);