postgis: scaling / optimization of shortest distance between two geometries (seeking recommendations)

Question

I have two tables in postgis:

• one with 1B rows containing points (POINT, 4326)::geometry
• the other, with 250k rows containing linestrings also in 4326
• Both tables have a GIST index.

I am trying to construct a query that looks at the entire combined geometry of all points and linestrings and tells me which points are nearest the linestrings (I don't need exact distance, closeness is fine).

My question is how to restructure / optimize the approach so that if I expand the size of geography I'm evaluating, the query will also scale efficiently (a return in minutes, few hours are fine...currently it takes days to run a full US query).

If I do the below query with a very small bounding box (ST_MakeEnvelope), I get results very quickly (less than 1sec).

SELECT id, lat, long, to_timestamp(unix_date) 

FROM point_data JOIN linestring_data on 

    ST_Intersects(linestring_data.wkb_geometry, ST_Buffer(point_data.geom, 0.0003)::geometry) 

WHERE geom && ST_MakeEnvelope(-75.188867,40.983303,-75.180389,40.989025, 4326);

Here is the query plan excluding WHERE (i.e. will query all points -- what I would like to do):

    Gather  (cost=1000.28..10429611075.47 rows=13454214512 width=64)
       Workers Planned: 2
       ->  Nested Loop  (cost=0.28..9084188624.27 rows=5605922713 width=64)
             ->  Parallel Seq Scan on point_data  (cost=0.00..54331353.71 rows=715436671 width=96)
             ->  Index Scan using test_wkb_geometry_geom_idx on linestring_data  (cost=0.28..12.61 rows=1 width=425)
                   Index Cond: (wkb_geometry && st_buffer(point_data.geom, '0.0003'::double precision))
                   Filter: _st_intersects(wkb_geometry, st_buffer(point_data.geom, '0.0003'::double precision))

I've tried below as a way to limit the number of points in the form of a buffer around the entire linestring geometry,

SELECT  distinct on (id) id, lat, long, unix_date, 
    ST_Distance(point_data.geom, linestring_data.wkb_geometry) 

FROM point_data JOIN linestring_data on 
    ST_Intersects(linestring_data.wkb_geometry, ST_Buffer(point_data.geom, 
        0.0003)::geometry) 

WHERE geom && 
    (SELECT ST_Buffer(ST_Collect(wkb_geometry), 0.0003) 
         FROM linestring_data);

Here is the query plan from explain:

    Unique  (cost=3107502.58..3114229.70 rows=141065 width=72)
    InitPlan 1 (returns $0)
    ->  Aggregate  (cost=1554.84..1554.85 rows=1 width=32)
       ->  Seq Scan on linestring_data linestring_data_1  (cost=0.00..1496.07 rows=23507 width=425)
    ->  Sort  (cost=3105947.73..3109311.29 rows=1345423 width=72)
     Sort Key: point_data.ad_id
     ->  Nested Loop  (cost=5811.55..2913801.67 rows=1345423 width=72)
           ->  Bitmap Heap Scan on point_data  (cost=5811.27..662541.58 rows=171705 width=96)
                 Recheck Cond: (geom && $0)
                 ->  Bitmap Index Scan on geom_idx  (cost=0.00..5768.34 rows=171705 width=0)
                       Index Cond: (geom && $0)
           ->  Index Scan using test_wkb_geometry_geom_idx on linestring_data  (cost=0.28..12.61 rows=1 width=425)
                 Index Cond: (wkb_geometry && st_buffer(point_data.geom, '0.0003'::double precision))
                 Filter: _st_intersects(wkb_geometry, st_buffer(point_data.geom, '0.0003'::double precision))

Questions:

• Is there a more efficient way to construct indexes so these queries run more efficiently?
• Is there a better way to construct the queries?

Update:

I decided to simplify by using bounding boxes rather than buffer. Here is the updated query:

    spatial=# explain SELECT  distinct on (id) id, lat, long, unix_date,
    spatial-#     ST_Distance(point_data.geom, linestring_data.wkb_geometry)
    spatial-#
    spatial-# FROM point_data JOIN linestring_data on
    spatial-#     ST_Intersects(linestring_data.wkb_geometry, ST_Expand(point_data.geom,
    spatial(#         0.0003))
    spatial-#
    spatial-# WHERE geom &&
    spatial-#     (SELECT ST_Envelope(wkb_geometry)
    spatial(#          FROM linestring_data)
    spatial-# ;
                                                        QUERY PLAN
    -------------------------------------------------------------------------------------------------------------------
     Unique  (cost=3107502.57..3114229.68 rows=141065 width=72)
       InitPlan 1 (returns $0)
         ->  Seq Scan on linestring_data linestring_data_1  (cost=0.00..1554.84 rows=23507 width=32)
       ->  Sort  (cost=3105947.73..3109311.29 rows=1345423 width=72)
             Sort Key: point_data.ad_id
             ->  Nested Loop  (cost=5811.55..2913801.67 rows=1345423 width=72)
                   ->  Bitmap Heap Scan on point_data  (cost=5811.27..662541.58 rows=171705 width=96)
                         Recheck Cond: (geom && $0)
                         ->  Bitmap Index Scan on geom_idx  (cost=0.00..5768.34 rows=171705 width=0)
                               Index Cond: (geom && $0)
                   ->  Index Scan using test_wkb_geometry_geom_idx on linestring_data  (cost=0.28..12.61 rows=1 width=425)
                         Index Cond: (wkb_geometry && st_expand(point_data.geom, '0.0003'::double precision))
                         Filter: _st_intersects(wkb_geometry, st_expand(point_data.geom, '0.0003'::double precision))
    (13 rows)

Seems like it may have made an improvement?

Answer 1

Some initial notes:

• using geographic CRS is a tricky business for anything distance related; a degree will not represent the same surface distance at different latitudes. You´d want to either use a projected CRS or cast to geography type for the implicit distance measurements
• in terms of performance, using buffers is generally a bad idea for proximity search (I think there´s even a hint at using ST_DWithin instead on the doc page for ST_Buffer - yes)

Two different approaches:

You can get a list of all POINTs that are within a threshold around a LINESTRING, as you tried in your queries by using a buffer; this will produce duplicates for POINTs found in thresholds around multiple LINESTRINGs (use DISTINCT ON (pt.id) as you did above, but the resulting POINT then is somewhat random) and exclude POINTs not in any of the thresholds (I included the ln.id to have a reference to the LINESTRING each POINT was found around):

SELECT pt.id AS point_id,
       ln.id AS line_id,
       pt.lat,
       pt.long,
       pt.to_timestamp(unix_date) AS date
FROM point_data AS pt
JOIN linestring_data AS ln
  ON ST_DWithin(ln.wkb_geometry, pt.geom, <distance_in_CRS_units>);

Or you can assign the nearest LINESTRING to each POINT; this effectively divides all POINTs into 'clusters' of their nearest LINESTRING, but will be slower:

SELECT pt.id AS point_id,
       ln.id AS line_id,
       pt.lat,
       pt.long,
       pt.to_timestamp(unix_date) AS date
FROM point_data AS pt
JOIN LATERAL (
    SELECT id
    FROM linestring_data
    ORDER BY pt.geom <#> linestring_data.wkb_geometry
    LIMIT 1
) AS ln
ON true;

Both queries are optimized to utilize the index: in the first query ST_DWithin does the index supported proximity search, in the second query I used the BBox KNN operator <#> to sort by distance. Here, if you are confident that no line is further apart than a certain distance, you could include a WHERE ST_DWithin(pt.geom, ln.geom, <distance>) in the subquery to filter the points before sorting by distance).
If true ground distance is of concern to you, do a geom::geography cast and provide the threshold distance in meter instead of degree. If you do this, consider ST_DWithin(pt.geom, ln.geom, <distance_in_meter>, false) to let measures be made on a sphere instead of the spheroid; this should increase speed for the otherwise a little more time consuming geography distance measuerement.