PostGIS distance query using a dynamic radius

Question

(Background: I'm more of a MySQL guy, I've only been using PostGIS for a few months to solve this particular problem set.)

I have a table that contains a point and a radius:

id      | location                                            | radius
1469989 | 0101000020E6100000091B9E5E292651C0DCF4673F529E5140  | 100
1948448 | 0101000020E6100000BB0A293F299863C03B8908FF22D35140  | 10

I want to do a search that basically says "get me all points within (20 km + radius) from me". In other words, some points need to "reach" farther than others and need to be returned from farther away, whereas others have a smaller radius.

My understanding is that I can't use spatial indexes for this since the DB doesn't know in advance what that radius will be.

My other option is to use ST_Buffer to create a buffer with the radius. But this seems to be pretty slow since there are hundreds of thousands of points, each of which get turned into a polygon with X points (default is 36). This is powering a real-time feature, so it needs to be really fast.

Any tips on what the best way to do this is? Or is this just not going to work and I need to look into more aggressive caching / in-memory stores?

Update

The current index looks like this:

CREATE INDEX IF NOT EXISTS "places_location" ON "places" USING gist("location");

and here's a sample query (20 is a dynamic value that IRL would be parameterized):

explain analyze SELECT id from places
WHERE ST_DWithin(ST_MakePoint('-79.387', '43.642')::geography, places.location, 20 *radius*1000);

Results:

Gather  (cost=1000.00..7076917.51 rows=68 width=8) (actual time=387.909..5444.545 rows=458668 loops=1)
  Workers Planned: 2
  Workers Launched: 2
  ->  Parallel Seq Scan on places  (cost=0.00..7075910.71 rows=28 width=8) (actual time=1017.367..5100.683 rows=152889 loops=3)
"        Filter: st_dwithin('0101000020E610000054E3A59BC4D853C01904560E2DD24540'::geography, location, (('20'::double precision * radius) * '1000'::double precision), true)"
        Rows Removed by Filter: 72442
Planning Time: 0.886 ms
JIT:
  Functions: 12
"  Options: Inlining true, Optimization true, Expressions true, Deforming true"
"  Timing: Generation 29.416 ms, Inlining 593.395 ms, Optimization 1217.528 ms, Emission 941.965 ms, Total 2782.304 ms"
Execution Time: 5598.359 ms

This is what it looks like without dynamic radius:

explain analyze SELECT id from places
WHERE ST_DWithin(ST_MakePoint('-79.387', '43.642')::geography, places.location, 20*1000);

Bitmap Heap Scan on places  (cost=19.15..10423.16 rows=68 width=8) (actual time=38.297..242.233 rows=8849 loops=1)
"  Filter: st_dwithin('0101000020E610000054E3A59BC4D853C01904560E2DD24540'::geography, location, '20000'::double precision, true)"
  Rows Removed by Filter: 2958
  Heap Blocks: exact=6553
  ->  Bitmap Index Scan on places_location  (cost=0.00..19.13 rows=363 width=0) (actual time=31.454..31.456 rows=11807 loops=1)
"        Index Cond: (location && _st_expand('0101000020E610000054E3A59BC4D853C01904560E2DD24540'::geography, '20000'::double precision))"
Planning Time: 8.443 ms
Execution Time: 247.098 ms

Answer 1

It seems like an odd detail of the engine that an ad-hoc lookup with varying filter parameters is not sargable - technically, the index wouldn't care. But an RDBS is a system with persisted data in mind, and this limitation is required to support functionality, and boost performance, elsewhere.

That being said, what you really need in order to gain access to index lookups is to persist those varying radii into known and types and sargable arguments - and about the only suitable way in this scenario would be to use a functional index on the expanded POINT geometries:

-- for GEOMETRY; note the implications of CRS and their units for <radius>!!!
CREATE INDEX ON <table> USING GIST ( ST_Expand(location, radius + <extra>) );

-- for GEOGRAPHY; unit of <radius> is meter!!!
CREATE INDEX ON <table> USING GIST ( _ST_Expand(location, radius + <extra>) );

and

SELECT
  *
FROM
  <table> AS t
WHERE
  -- either for GEOMETRY or GEOGRAPHY; note that the expression needs to match that of the index definition EXACTLY
  [_]ST_Expand(t.geom, t.radius + <extra>) && <INPUT_POINT>
  AND
  ST_DWithin(
    <INPUT_POINT>,
    t.geom,
    t.radius + <extra>
  )
;

It seems equally odd that, given the functional index is in place, ST_DWithin cannot use it by itself even though internally it used to do the exact same index invocation.

Update:

If you need extra to be dynamic at query time, but essentially static per query, you'd simply want to change the above to this:

-- for GEOMETRY; note the implications of CRS and their units for <radius>!!!
CREATE INDEX ON <table> USING GIST ( ST_Expand(location, radius) );

-- for GEOGRAPHY; unit of <radius> is meter!!!
CREATE INDEX ON <table> USING GIST ( _ST_Expand(location, radius) );

and

SELECT
  *
FROM
  <table> AS t
WHERE
  -- either for GEOMETRY or GEOGRAPHY; note that the expression needs to match that of the index definition EXACTLY
  [_]ST_Expand(t.geom, t.radius) && [_]ST_Expand(<INPUT_POINT>, <extra>)
  AND
  ST_DWithin(
    <INPUT_POINT>,
    t.geom,
    t.radius + <extra>
  )
;

---

However:

Point to point distance is cheap. Always, even for spherical/spheroidal calculations. Index lookups are only cheap in relation, e.g. when their initial load up and random row access penalty are neglectable compared to their impact on e.g. large Cartesian joins. Their efficiency degenerates fast with increasing matches.

With that in mind, for a single input POINT against a few million stored POINTs it will likely be faster to see 2 or more parallel worker nodes spun up and do sequential scans for higher match counts!

Answer 2

If the radius values are fairly similar, select first using the greatest radius (always, as a constant, so the index can be used) and refine the results using the variable width radius

...
where a.geom && st_expand(g.geom,25000) 
and st_dwithin(a.geom,b.geom, a.radius)

Answer 3

• First is your data stored as geography type or geometry type. I would suggest changing it to geography especially if you are a newbie. Otherwise you'd need to understand spatial reference systems, which is another topic altogether.

• I'm also assuming your radius is stored in kilometers thus the need to multiply by 1000, since geography always measures in meters

• Make sure your points are in order of longitude,latitude NOT longitude,latitude, which is a common mistake of people especially mappers (since many mapping platforms expect latitude,longitude)

Here are some self contained examples using your example points

    CREATE TABLE pois(id integer, geog geography(POINT,4326), radius float8, CONSTRAINT pk_pois PRIMARY KEY(id));
    
    INSERT INTO pois(id, geog, radius)
    VALUES (1469989,'0101000020E6100000091B9E5E292651C0DCF4673F529E5140'::geography(POINT,4326), 100 ), 
    (1948448,'0101000020E6100000BB0A293F299863C03B8908FF22D35140'::geography(POINT,4326), 10 );
    
    CREATE INDEX ix_pois_geog ON pois USING gist(geog);
    
    -- checking a table of points
    SELECT p.*, pois.id
    FROM points_to_check AS p INNER JOIN pois ON ST_DWithin(pois.geog, pop.geog, pois.radius*1000 + 20*1000);
    
    -- checking a single point NOTE: coordinates are longitude,latitude
    SELECT pois.*
    FROM pois 
    WHERE ST_DWithin(pois.geog, ST_Point(70.10,-68)::geography(POINT,4326), pois.radius*1000 + 20*1000);