2

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?

  • in other words: you want to find the nearest linestring to each point? – ThingumaBob Mar 15 '18 at 15:30
  • vice versa - I want the result to be all the points that are nearest to the linestrings. – Matt Johnson Mar 15 '18 at 15:54
  • okay, and by a threshold then I take it (assuming from your use of buffers)? or do you want something like the top ten? because, if you want to uniquely reference every point to one distinct linestring, you´d want to go the other way – ThingumaBob Mar 15 '18 at 15:56
  • After thinking about your question on buffers, I realized that for a given linestring, buffers is going to create thousands of points that have to be evalauted. So I decided to simplify to an envelope. See updated question above. – Matt Johnson Mar 15 '18 at 16:16
  • Yes, I am going to try your first approach. I'm currently converting geometry to geography and adding indexes. The linestring table took 1 second to convert, the point data is probably going to take a few hours. I'll let you know how it goes once I have everything set up. Thanks again for sharing your recommendations. – Matt Johnson Mar 16 '18 at 2:26
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.

  • as an addition to the second approach: the query needs around 6 seconds for 150k points and 350k lines on my mid-tech machine, and 4 with the WHERE ... filter; choose a distance so that you get each point returned in the result (fewer means the distance is too short). if you plan on inserting more points periodically, consider using a trigger to update those new points with their nearest line ids instead of running the query over and over. for updating lines however, you'd need to run it again. – ThingumaBob Mar 15 '18 at 17:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.