ST_DWithin optimization for multiple SRIDs

Question

In my PostgreSQL 9.3 database I have a table named location that includes a coordinate column used to store POINT geometry. These points are created using SRID 4326. Under certain circumstances we ST_Transform these coordinates to SRID 900913 to filter them using distances in meters.

For example, to use ST_WDithin to find all locations with coordinates within 10000 meters of a given coordinate the query looks like this:

SELECT *
FROM LOCATION 
WHERE ST_DWithin(ST_Transform(location.coordinate, 900913), ST_Transform(ST_GeomFromText('POINT(-74.005941 40.712784)', 4326), 900913), 10000)
ORDER BY ST_Distance_Sphere(location.coordinate, ST_GeomFromText('POINT(-74.005941 40.712784)', 4326))
LIMIT 100

The query plan for this statement looks like this:

Limit  (cost=19207.98..19208.23 rows=100 width=36)
  ->  Sort  (cost=19207.98..19209.81 rows=729 width=36)
        Sort Key: (_st_distance(geography(coordinate), '0101000020E6100000282D5C56618052C0588E90813C5B4440'::geography, 0::double precision, false))
        ->  Seq Scan on location  (cost=0.00..19180.12 rows=729 width=36)
              Filter: ((st_transform(coordinate, 900913) && '010300002031BF0D000100000005000000D42FBDEAFB765FC187C3AD4ED1EB5241D42FBDEAFB765FC187C3AD4E59FF5241D42FBDEA73635FC187C3AD4E59FF5241D42FBDEA73635FC187C3AD4ED1EB5241D42FBDEAFB765FC187C3AD4ED1EB5241'::geometry) AND ('010100002031BF0D00D42FBDEA376D5FC187C3AD4E95F55241'::geometry && st_expand(st_transform(coordinate, 900913), 10000::double precision)) AND _st_dwithin(st_transform(coordinate, 900913), '010100002031BF0D00D42FBDEA376D5FC187C3AD4E95F55241'::geometry, 10000::double precision))

This works, but is very slow. All coordinates involved are transformed to SRID 900913 to be able to work in meters. Through testing I have found that removing ST_Transform speeds this query up dramatically.

I also experimented with creating a circular POLYGON using ST_Buffer and then testing to see if location.coordinate intersects with this polygon. To do this I ST_Transform the input coordinate to SRID 900913, use ST_Buffer to draw a circle with a radius in meters, and then ST_Transform that polygon to SRID 4326 for comparison against the coordinates in my location table. The query looks like this:

SELECT *
FROM LOCATION 
WHERE location.coordinate && ST_Transform(ST_Buffer(ST_Transform(ST_GeomFromText('POINT(-74.005941 40.712784)', 4326), 900913), 10000), 4326))
ORDER BY ST_Distance_Sphere(location.coordinate, ST_GeomFromText('POINT(-74.005941 40.712784)', 4326))
LIMIT 100

In my tests this second query runs much, much faster that the first. It even runs a bit faster that my version of the query using ST_DWithin minus the ST_Transform. From all I have read it, seems ST_DWithin should be the fastest way to perform such a search. The answer to this question indicates I should be able to create an index of my coordinates transformed to a different SRID: PosgtreSQL Optimize Query with st_transform, st_makepoint, and st_contains

I attempted this by running:

CREATE INDEX idx_location_coordinate_900913
 ON LOCATION
 USING gist
 (ST_Transform(coordinate, 900913))
 WHERE coordinate IS NOT NULL;

After creating this index I saw no speed improvements when running my original query. I found it odd that this command completed successfully very quickly and rebuilding this index happened very quickly. There are tens of thousands of rows in the location table, so I imagined creating this index would be a time consuming process. Did I create it incorrectly?

Is there something I can do to speed up ST_DWithin when transforming my points? Is there a major flaw to my approach here?

EDIT: I'm adding the execution plan for the initial query above.

Answer 1

I would suggest creating a function in PostgreSQL to handle the conversion from meters to decimal degrees. By doing this, you can avoid transforming every single row, and instead the ST_DWithin function is working in the native projection.

You'll want to double check the math and I would suggest comparing the results to using ST_Transform, but I'm pretty sure this is close. The function I use typically uses miles instead of meters, but I took a quick stab at adding the conversion from meters.

In the function, the value 3960 is an estimation of Earth's diameter, and 1609.34 handles changing from miles to meters. I make no guarantee that this is as precise or accurate as ST_Transform, but it should be much better for performance since it doesn't have to convert every row.

CREATE FUNCTION meters_to_decimal_degrees(meters double precision)
RETURNS double precision AS
$BODY$
    SELECT (($1 * 180 * 1609.34) / ( 3960 * pi() ) ) AS decimal_degrees
$BODY$
LANGUAGE sql IMMUTABLE SECURITY DEFINER
;


ALTER FUNCTION public.meters_to_decimal_degrees(double precision) SET search_path=public, pg_temp;

With this, your query could be altered to be:

SELECT *
    FROM LOCATION 
    WHERE ST_DWithin(location.coordinate, ST_GeomFromText('POINT(-74.005941 40.712784)', 4326), meters_to_decimal_degrees(10000))
    ORDER BY ST_Distance_Sphere(location.coordinate, ST_GeomFromText('POINT(-74.005941 40.712784)', 4326))
    LIMIT 100

Hopefully that helps.