PostGIS: Convert meters to arbitrary spatial reference units?

Question

I'm writing a query that needs to be able to work with data of any SRID. The basic problem is very simple - I get two pieces of information:

1. A point (SRID = 900913)
2. A distance (unit = meters)

I [would like to] transform these into the native SRID, and then run a query to find features within the given distance of the point:

proj_point = ST_Transform(point, (SELECT ST_SRID(geom) FROM table LIMIT 1));
proj_dist  = Something...  But what?

SELECT * FROM table
    WHERE ST_Distance(geom, proj_point) < proj_dist
    ORDER BY geom <-> proj_point
    LIMIT 10;

My problem is the conversion of the units. I can look up the SRID of the native data and transform the point no problem. But how do I do the same with the distance?

I've seen answers to related questions, but they have all involved reprojecting the queried data. Holy overkill! This should be doable with a single lookup and multiplication...

Pointers, anyone? Thanks very much in advance!

Answer 1

Rather than calling st_distance on all your objects, which most likely alread prevents your index from being used, why not make your point into a "disc" and then find all enclosed features.

Geometries are usually very bad at handling measures like meters. However, geographies are good with them. So what we can do is transform your point to a geography, buffer it a few metres and then make it a geometry again. Then use that polygon to do a quick indexed search and finally do an accurate search to omit false positives.

As you already successfully transformed your point to the correct srid, I will use that.

proj_disc = geometry(st_buffer(geography(proj_point), distance))

SELECT * FROM table
WHERE geom && proj_disc
  AND st_within(geom, proj_disc)
ORDER BY st_distance(geography(geom), geography(proj_point)) DESC
LIMIT 10;

Note that I'm ordering by the distance after translating the points to geographies, to get a distance in meters as not to get skewed results because of the used projection.

Also, I'm using st_within. This will work fine for points, though if you're using polygons or lines, you might want to check for intersection or something else, since this only shows results that are fully within the distance of the point.

The first part of the where clause is what will make it all be quick, since that will trigger the spatial index. The slower operations (st_within and st_distance) are then only done for the ones that are likely to match.

-update-

According to Paul Ramsey there is a simpler way of achieving the desired result, which should still make use of a spatial index:

SELECT * FROM table
WHERE ST_DWithin(geom::geography, proj_point, distance)
ORDER BY ST_Distance(geom::geography, proj_point) DESC
LIMIT 10

(You may have to cast your point to a geography as well, haven't tried.)

Answer 2

My answer is that I dont know how to transform distance but this one works.

SELECT ST_DISTANCE(ST_TRANSFORM(a.geom , 4326), ST_TRANSFORM(b.geom, 4326))     
FROM tablea AS a. tableb AS b 
WHERE a.id =1 AND b.id = 2 

It just converts geoms to wanted SRID and measures distance, which is always measured in using SRID. So if you use SRID 3067, then you use meters; if you use 4236 then you use degree (I think it uses degrees) as unit.

Your code changes to :

SELECT a.* FROM tablea as a , tableb as b
WHERE ST_Distance(ST_TransForm(a.geom, 3067), ST_Transform(b.geom, 3067) < 10 
AND b.id = 1 -- you compare all data from table a to this 
ORDER BY ST_Distance(ST_TransForm(a.geom, 3067), ST_Transform(b.geom, 3067) ASC 
LIMIT 10;

Code Translates everything to 3067 (tm35fin) and returns everything from a table with distance to b.geom in b table with id 1. This assumes that you have the correct SRID in every geom. If not you need to set it with SetSRID (st_?). If you want the answer in degree decimal then use 4326 (I'm not 100% sure about that, but I recall it uses degree units).