[postgis-users] distance query is not returning expected result

[Mark Cave-Ayland]

Jimmy Forrester wrote:

> Hey all, I'm new to postgis and implementing a few geographic features 
> in to my site. I'm having a bit of difficulty getting my distance 
> calculations working:
> 
> select ST_Distance(    ST_Transform(ST_SetSRID(ST_Point(53.7996374, 
> -1.54911),     4326), 2163),
>             ST_Transform(ST_SetSRID(ST_Point(53.6451179, -1.784876),     
> 4326), 2163)
>         ) / 1000
> 
> 
> In the query above I'm testing the distance between Leeds & 
> Huddersfield. I'm trying to convert the 4326 srid geom to  srid 2163 
> which is meters (i think!). For an unknown reason this is returning 64 
> km which i know is incorrect.
> 
> I've plotted the coords on google maps and they are definitely correct, 
> google maps says the distance is about 24 km. see:
> 
> http://www.google.com/maps?f=d&source=s_d&saddr=53.7996374,+-1.54911&daddr=53.6451179,+-1.784876&hl=en&geocode=FdXqNAMdylzo_w%3BFT6PMgMd1MPk_w&mra=ls&dirflg=w&doflg=ptk&sll=53.722107,-1.645889&sspn=0.161297,0.418854&ie=UTF8&ll=53.72231,-1.667175&spn=0.161296,0.418854&t=h&z=12 
> <http://www.google.com/maps?f=d&source=s_d&saddr=53.7996374,+-1.54911&daddr=53.6451179,+-1.784876&hl=en&geocode=FdXqNAMdylzo_w%3BFT6PMgMd1MPk_w&mra=ls&dirflg=w&doflg=ptk&sll=53.722107,-1.645889&sspn=0.161297,0.418854&ie=UTF8&ll=53.72231,-1.667175&spn=0.161296,0.418854&t=h&z=12>
> 
> Any help on this would be great,
> 
> Thanks,
> Jimmy

Hi Jimmy,

 From your above query you appear to have the latitude/longitude the 
wrong way around - PostGIS follows the normal computing convention of X 
then Y (or longitude then latitude). Also the normal SRID for usage 
within the UK is 27700 which matches the standard Ordnance Survey grid 
so you'll want to use this instead of 2163.

Note that the spatial_ref_sys table supplied with PostGIS contains an 
error for the 27700 transformation, and so first of all you'll need to 
do this on your PostGIS database:


postgis=# UPDATE spatial_ref_sys SET proj4text = '+proj=tmerc +lat_0=49 
+lon_0=-2 +k=0.999601 +x_0=400000 +y_0=-100000 +ellps=airy +units=m 
+no_defs +datum=OSGB36' WHERE srid = 27700;
UPDATE 1


And now your query looks like this:


select ST_Distance(
	ST_Transform(
		ST_SetSRID(ST_Point(-1.54911, 53.7996374), 4326),
		27700),
	ST_Transform(
		ST_SetSRID(ST_Point(-1.784876, 53.6451179), 4326),
		27700)
	) / 1000;

which returns:

      ?column?
------------------
  23.1860213462966
(1 row)


HTH,

Mark.

-- 
Mark Cave-Ayland
Sirius Corporation - The Open Source Experts
http://www.siriusit.co.uk
T: +44 870 608 0063