Real (Great Circle) distance in PostGIS with lat / long SRID?

Question

I am using lat / long SRID in my PostGIS database (-4326). I would like to find the most suitable points for this point. I tried to do

ORDER BY ST_Distance(point, ST_GeomFromText(?,-4326)) 

which gives me good results in the lower 48 states, but in Alaska it gives me garbage. Is there a way to do real distance calculations in PostGIS, or will I need to provide a buffer of a reasonable size and then calculate the large circle distances and sort the results in the code after that?

Answer 1

You are looking for ST_distance_sphere (dot, dot) or st_distance_spheroid (dot, dot).

Cm:

http://postgis.refractions.net/documentation/manual-1.3/ch06.html#distance_sphere http://postgis.refractions.net/documentation/manual-1.3/ch06.html#distance_spheroid

This usually refers to the geodetic or geodetic distance ... while the two terms have slightly different meanings, they are usually used interchangeably.

Alternatively, you can project data and use the standard function st_distance ... it is practical only for short distances (using the UTM or state plane) or if all distances are relative to one or two points (equidistant projections).

Answer 2

PostGIS 1.5 processes the true distances of the globe using long armor and meters. He knows lat / long is angular in nature and has a 360 degree line

Answer 3

This is from SQL Server, and I'm using Haversine for a ridiculously fast distance that could suffer from your problem in Alaska (maybe shut down by a mile):

 ALTER function [dbo].[getCoordinateDistance] ( @Latitude1 decimal(16,12), @Longitude1 decimal(16,12), @Latitude2 decimal(16,12), @Longitude2 decimal(16,12) ) returns decimal(16,12) as /* fUNCTION: getCoordinateDistance Computes the Great Circle distance in kilometers between two points on the Earth using the Haversine formula distance calculation. Input Parameters: @Longitude1 - Longitude in degrees of point 1 @Latitude1 - Latitude in degrees of point 1 @Longitude2 - Longitude in degrees of point 2 @Latitude2 - Latitude in degrees of point 2 */ begin declare @radius decimal(16,12) declare @lon1 decimal(16,12) declare @lon2 decimal(16,12) declare @lat1 decimal(16,12) declare @lat2 decimal(16,12) declare @a decimal(16,12) declare @distance decimal(16,12) -- Sets average radius of Earth in Kilometers set @radius = 6366.70701949371 -- Convert degrees to radians set @lon1 = radians( @Longitude1 ) set @lon2 = radians( @Longitude2 ) set @lat1 = radians( @Latitude1 ) set @lat2 = radians( @Latitude2 ) set @a = sqrt(square(sin((@ lat2-@lat1 )/2.0E)) + (cos(@lat1) * cos(@lat2) * square(sin((@ lon2-@lon1 )/2.0E))) ) set @distance = @radius * ( 2.0E *asin(case when 1.0E < @a then 1.0E else @a end ) ) return @distance end 

Vicenty is slow but accurate to 1 mm (and I only found javascript imp):

 /* * Calculate geodesic distance (in m) between two points specified by latitude/longitude (in numeric degrees) * using Vincenty inverse formula for ellipsoids */ function distVincenty(lat1, lon1, lat2, lon2) { var a = 6378137, b = 6356752.3142, f = 1/298.257223563; // WGS-84 ellipsiod var L = (lon2-lon1).toRad(); var U1 = Math.atan((1-f) * Math.tan(lat1.toRad())); var U2 = Math.atan((1-f) * Math.tan(lat2.toRad())); var sinU1 = Math.sin(U1), cosU1 = Math.cos(U1); var sinU2 = Math.sin(U2), cosU2 = Math.cos(U2); var lambda = L, lambdaP = 2*Math.PI; var iterLimit = 20; while (Math.abs(lambda-lambdaP) > 1e-12 && --iterLimit>0) { var sinLambda = Math.sin(lambda), cosLambda = Math.cos(lambda); var sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) + (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda)); if (sinSigma==0) return 0; // co-incident points var cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda; var sigma = Math.atan2(sinSigma, cosSigma); var sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma; var cosSqAlpha = 1 - sinAlpha*sinAlpha; var cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha; if (isNaN(cos2SigmaM)) cos2SigmaM = 0; // equatorial line: cosSqAlpha=0 (§6) var C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha)); lambdaP = lambda; lambda = L + (1-C) * f * sinAlpha * (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM))); } if (iterLimit==0) return NaN // formula failed to converge var uSq = cosSqAlpha * (a*a - b*b) / (b*b); var A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq))); var B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq))); var deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)- B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM))); var s = b*A*(sigma-deltaSigma); s = s.toFixed(3); // round to 1mm precision return s; }