4

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.

Green: the given point; red: the points and circles which fit 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.

2
  • it is required to scan the table sequentially, you are effectively looping over it ('for each point in the table, see if it is within radius of the given geometry'); an index would only make sense if you compared against a table of geometries! PG's planner is correct here, there is no performance gain, and I suggest to not work with polygons at all.
    – geozelot
    Sep 28, 2019 at 9:03
  • since you didn't ask ;-): if you'd pass in that static geometry as CTE or subquery table, a JOIN ON ST_DWithin(...) would enable the use of the index; however, if the planner decides that most of your points are affected, it would still seq scan.
    – geozelot
    Sep 28, 2019 at 9:29

1 Answer 1

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.

7
  • I like the idea with buffer around radius. I'd prefer this approach, however, as stated in the question, I haven't managed to make the query use index. Do you have any idea how to fix that? And can you please show how to make a bounding box for my radius? I haven't found examples with geography and in the docs the functions take only geometry type.
    – 1valdis
    Sep 23, 2019 at 13:22
  • @1valdis In your example, you just created an index using a buffer. This kind of index is not really easy to trigger and maybe the planner didn't understood? Also be sure to analyze your table after the index creation. In any case if you create the field, and use a classical index on it (and analyze) it's sure that it will be used. Sep 23, 2019 at 14:30
  • @1valdis For the bounding box of the circle, you can do box2d(st_buffer(point,radius))::geometry Sep 23, 2019 at 14:34
  • Are you sure it'll take into account the Earth's curvature? Your functions calls don't seem to make any use of geography; although I see that I may use custom geodesic_buffer function to solve this. But still, box2d applies to geometries only... Also, I ran analyze and my index still isn't used in queries above. By "classical index" you mean GiST one?
    – 1valdis
    Sep 23, 2019 at 14:56
  • The earth curvature is taken into account in the st_buffer function, and box2D only take the max x and y (for a geom in 4326 it's longitude and latitude) Geodesic_buffer is only for really precise calculus but I don't see real practical use for this kind of precision (I mean a classical GPS have 5m precision, and the geodesic itself is not perfect eather). Plus you take a buffer around so no need for perfect precision. Just be sure to cast ::geography whenever you can: (box2d(st_buffer(point::geography,radius)::geometry)::geometry)::geography as geog and for your index: USING GIST (geog); Sep 23, 2019 at 15:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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