Real (Great Circle) Abstand in PostGIS mit lat / long SRID?
Frage
Ich bin mit einem lat / long SRID in meiner PostGIS-Datenbank (-4326). Ich möchte die nächsten Punkte an einem bestimmten Punkt in einer effizienten Weise zu finden. Ich versuchte dabei ein
ORDER BY ST_Distance(point, ST_GeomFromText(?,-4326))
, das gibt mir ok Ergebnisse in den unteren 48 Staaten, aber bis in Alaska gibt es mir Müll. Gibt es eine Möglichkeit reale Entfernungsberechnungen in PostGIS zu tun, oder werde ich habe einen recht großen Puffer zu geben und dann die großen Kreis Entfernungen berechnen und die Ergebnisse in dem Code danach sortieren?
Lösung
Sie suchen ST_distance_sphere (Punkt, Punkt) oder st_distance_spheroid (Punkt, Punkt).
Siehe auch:
http://postgis.refractions.net/documentation/manual- 1.3 / ch06.html # distance_sphere http://postgis.refractions.net/documentation/manual-1.3/ch06 .html # distance_spheroid
Dies ist in der Regel auf einem geodätischen oder geodätischen Abstand bezeichnet ... während die beiden Begriffe etwas andere Bedeutung haben, neigen sie dazu, synonym verwendet werden.
Alternativ können Sie die Daten projizieren und die Standard-ST_Distance Funktion ... Dies ist nur praktisch über kurze Distanzen (UTM oder Zustandsebene verwendet wird) oder wenn alle Entfernungen sind relativ zu einer ein oder zwei Punkte (gleich weit Projektionen).
Andere Tipps
PostGIS 1.5 Griffe wahr Globus Entfernungen mit lat longs und Meter. Es ist bekannt, dass lat / long Winkel ist in der Natur und hat eine 360-Grad-Linie
Dies ist von SQL Server, und ich verwende Haversine für einen lächerlich schnell Abstand, der von Ihrem Alaska Problem leiden kann (kann durch eine Meile aus sein):
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 ist langsam, aber genau innerhalb von 1 mm (und ich fand nur einen JavaScript-imp davon):
/*
* 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;
}