I have a table with 20000 records and with those columns Shop, longitude, latitude, geom

I would like to find a way to find the pair of ths farthest shops between a selection of shops

For example we have outletnames : A,B,C,D,E,F,G,H,I,J,K,L,M,N,O...

I would like to know what is the bigest distance between all those shops A B C D and E A-B = 4 km A-C- = 2 km A-D = 1 km A-E = 1,3 km B-C = 20 km C-D = 50 km ....

I would like to have an optimised query that would give me as a result 50 km

Thank you

1 Answer 1


To measure distances in Km you need to transform your geometry's coordinates from long/lat (spherical) to x/y (planar). In my example I am using UTM zone 35N, EPSG 32635 - you will have to pick something appropriate based on your data's location.

The query below can be optimized by avoiding the st_transform call (e.g. by creating an intermediate table with geometry in EPSG 32635).

The other thing that can be improved is that distances are computed as a query over a cartesian product (the table is joined with itself) so you unnecessarily compute both the distance from A to B and the distance from B to A.

select a.outletname, b.outletname, 
st_distance(st_transform(a.geom, 32635), st_transform(b.geom, 32635)) 
table_name a,
table_name b,
where a.outletname != b.outletname
and a.outletname in ('A', 'B', 'C', 'D')
and b.outletname in ('A', 'B', 'C', 'D')
and st_distance(st_transform(a.geom, 32635), st_transform(b.geom, 32635)) > 50000 
order by st_distance(st_transform(a.geom, 32635), st_transform(b.geom, 32635)) desc;

whuber's idea with the convex hull makes sense when the data set is large enough to prevent computing distances on the cartesian product.

I've checked on my dataset and, as expected, the greatest distance is equal when computed by the two different methods. For reference, the distances between the convex hull points query:

select st_distance(a.geom, b.geom) from
(select (st_dumppoints(st_convexhull(st_collect(st_transform(way, 32635))))).geom from planet_osm_point a where a.amenity='pharmacy' and st_contains((select way from planet_osm_polygon where name='Sector 1'), a.way)) a,
(select (st_dumppoints(st_convexhull(st_collect(st_transform(way, 32635))))).geom from planet_osm_point a where a.amenity='pharmacy' and st_contains((select way from planet_osm_polygon where name='Sector 1'), a.way)) b
order by st_distance(a.geom, b.geom) desc;

There's an additional problem with the convex hull approach - it works on aggregated geometry hence it anonymizes which point belongs to which feature, making it more difficult to compute outlet names.

The solution would probably be to compute the convex hull and then identify (by spatial queries), which feature geometries belong to the convex hull and run the distance query on those features only. E.g.:

The geometry of the first point of the convex hull:

osm=# select (st_dumppoints(st_convexhull(st_collect(way)))).geom from planet_osm_point a where a.amenity='pharmacy' limit 1;
(1 row)

The corresponding feature from the table:

osm=# select name, osm_id, way from planet_osm_point a  where a.way = '010100002031BF0D00D120685EC2064241D3F3D3FDA9AE5541';
   name    |   osm_id   |                        way                         
 Dentafarm | 1395848741 | 010100002031BF0D00D120685EC2064241D3F3D3FDA9AE5541
(1 row)
  • Out of curiosity, how long would this query take, given that it appears to compute about 20,000,000 distances? There are many ways to avoid computing all these distances: a quick and dirty one, if it's available, is to extract the vertices of the convex hull of all the points and then query the greatest distance among just those vertices. This would likely be two to three orders of magnitude faster.
    – whuber
    May 9, 2012 at 14:53
  • For point geometries that do not require re-projection st_distance will be somewhere between 5ms and 20ms per call. However, filtering based on outlet name has a chance to greatly reduce the dataset, hence the number of st_distance calls.
    – diciu
    May 9, 2012 at 15:00
  • For the convex hull idea, there's built in support - there's a st_convexhull function that takes in an aggregated geometry. The geometry returned by st_convexhull can then be parsed, one vertex at a time using st_npoints.
    – diciu
    May 9, 2012 at 15:06

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.