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:
- I have a set of points.
- For each point there's a radius, somewhere from 1 to 100 kms.
- I have other given point.
- 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.
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.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.