Finding all circles point lies in using PostGIS

Question

Suppose you have a task to find all points within a certain radius. That's where you use ST_DWithin which is indexable and it will deal with it.

My task is the opposite from that:

1. I have a set of points.
2. For each point there's a radius, somewhere from 1 to 100 kms.
3. I have other given point.
4. I want to select every point from p.1 where my given point lies within a radius from p.2 of that point.

On a figure below, green point is given, and red points are those fitting the condition.

The table setup code would be:

create table test_coordinates(
    coordinates geography,
    -- in meters
    radius integer not null default random() * 20000 + 1000
);
create index on test_coordinates using gist(coordinates);

I filled the coordinates from another table, so can't post them there. Some random points would do well, I have a couple hundred of them covering single country.

My first thought was to use ST_DWithin, like that:

select * from test_coordinates where st_dwithin(coordinates, '0101000020E6100000DA6A20A7F9F44D4004221F4F26812540', test_coordinates.radius);

But it won't be indexed anyhow because it will still need to do sequential scan. I tested it using set enable_seqscan = off and it still resorted to sequential scan, so I suppose it means that this is the only way.

So the next thought was about polygons, because they're indexable and I can do point-in-polygon using index (supposedly).

-- 1. geodesic buffer function taken from https://stackoverflow.com/a/13872887/5748383
CREATE OR REPLACE FUNCTION geodesic_buffer(geom geometry, dist double precision,
                                           num_seg_quarter_circle integer)
  RETURNS geometry AS $$
  SELECT ST_Transform(
    ST_Buffer(ST_Point(0, 0), $2, $3),
      ('+proj=aeqd +x_0=0 +y_0=0 +lat_0='
       || ST_Y(ST_Centroid($1))::text || ' +lon_0=' || ST_X(ST_Centroid($1))::text),
      ST_SRID($1))
  $$ LANGUAGE sql IMMUTABLE STRICT COST 100;
CREATE OR REPLACE FUNCTION geodesic_buffer(geog geography, dist double precision,
                                           num_seg_quarter_circle integer)
  RETURNS geography AS 'SELECT geodesic_buffer($1::geometry, $2, $3)::geography'
  LANGUAGE sql IMMUTABLE STRICT COST 100;
-- 2. index
create index on test_coordinates using gist(geodesic_buffer(test_coordinates.coordinates, test_coordinates.radius, 32));

However, there are 2 problems with that.

First, none of these queries use index by their query plans.

set enable_seqscan = off;
explain analyse select * from test_coordinates where st_intersects('0101000020E6100000DA6A20A7F9F44D4004221F4F26812540'::geography, geodesic_buffer(test_coordinates.coordinates, test_coordinates.radius, 32));
explain analyse select * from test_coordinates where st_coveredby('0101000020E6100000DA6A20A7F9F44D4004221F4F26812540'::geography, geodesic_buffer(test_coordinates.coordinates, test_coordinates.radius, 32));
explain analyse select * from test_coordinates where st_covers(geodesic_buffer(test_coordinates.coordinates, test_coordinates.radius, 32), '0101000020E6100000DA6A20A7F9F44D4004221F4F26812540');

The resulting query plan for all of them is as follows:

Seq Scan on test_coordinates  (cost=10000000000.00..10000000115.57 rows=67 width=36) (actual time=3.193..62.728 rows=132 loops=1)
  Filter: st_covers((st_transform(_st_buffer('010100000000000000000000000000000000000000'::geometry, (radius)::double precision, 'quad_segs=32'::cstring), ((('+proj=aeqd +x_0=0 +y_0=0 +lat_0='::text || (st_y(st_centroid((coordinates)::geometry)))::text) || ' +lon_0='::text) || (st_x(st_centroid((coordinates)::geometry)))::text), st_srid((coordinates)::geometry)))::geography, '0101000020E6100000DA6A20A7F9F44D4004221F4F26812540'::geography)
  Rows Removed by Filter: 69
Planning time: 2.965 ms
Execution time: 64.014 ms

What prevents it from using index, despite I index by the same result of geodesic_buffer call as in select statement?

The second problem is that I can't get ideally circular polygon. Using 32 segments per quarter circle, I still have 60 meters of error on 200 km radius. At the same time, a single circle already has a size of 2 kilobytes. This approach seems to bloat as soon as I have tens of thousands of records.

How do I solve this task in a performant (using indexes) and precise (ideally having true circle from radius, not an approximation) way?

I'm working with geography type because the points are on the Earth's surface.

Answer 1

I see 2 main ways of resolving this:

First, you can add like you said a geography field with the circle corresponding to your radius, and index it. If you're concerned about precision and size, you should add a buffer around your radius (for exemple +100m, a perfect one should take into account the radius and the theorical imprecision to be sure) and you can use less than 32 points if you play with this buffer. This way, you can prefilter easily your points by using the index on this polygon, and then filter again with more precision by using st_distance. This first filter will greatly help in performance, and because you added a buffer you can be sure not to loose any point.

The other way that I see is to add a field "radius_category" to your p1. For exemple, every point that have a radius between 10 and 100m will be marked "100". Then you can index this field, and you basically do a request for each of these categories (with a st_within 100m for all the points that are in the 100 category for exemple, and then a st_distance for the precision test). You will have to do a for loop probably, or you can manage with a good cte before, but it's better (and more efficient) than to do a for loop for all the points. But to be honest I cannot know how the planner will react if you filter with category then with a st_within with your geom, so I'm not sure how efficient this would be.

My favorite way of doing this would definitly be the first way. And if you're really concerned about size, you can just do a bounding box around the point instead of a st_buffer (it's just a squared circle after all ;) ), and it will work the same, you will just have a first filter less efficient.