我在我的PostGIS數據庫中使用緯度/長度SRID(-4326)。我想以有效的方式找到給定點的最近點。我試圖做一個使用lat/long SRID的PostGIS中的真實(大圓圈)距離?
ORDER BY ST_Distance(point, ST_GeomFromText(?,-4326))
這給我好結果在較低的48個州,但在阿拉斯加,它給了我垃圾。有沒有辦法在PostGIS中進行實際的距離計算,或者我需要給出一個合理大小的緩衝區,然後計算大圓距離並在代碼中對結果進行排序?
我在我的PostGIS數據庫中使用緯度/長度SRID(-4326)。我想以有效的方式找到給定點的最近點。我試圖做一個使用lat/long SRID的PostGIS中的真實(大圓圈)距離?
ORDER BY ST_Distance(point, ST_GeomFromText(?,-4326))
這給我好結果在較低的48個州,但在阿拉斯加,它給了我垃圾。有沒有辦法在PostGIS中進行實際的距離計算,或者我需要給出一個合理大小的緩衝區,然後計算大圓距離並在代碼中對結果進行排序?
您正在尋找ST_distance_sphere(點,點)或st_distance_spheroid(點,點)。
參見:
http://postgis.refractions.net/documentation/manual-1.3/ch06.html#distance_sphere http://postgis.refractions.net/documentation/manual-1.3/ch06.html#distance_spheroid
這通常被稱爲一個短程線或測地距離......而兩個術語具有稍微不同的含義,它們傾向於被互換使用。
或者,您可以投影數據並使用標準的st_distance函數......這隻適用於短距離(使用UTM或狀態平面)或所有距離都相對於一個或兩個點(等距投影) 。
這是從SQL Server,我用半正矢一個可笑的快速距離可以從阿拉斯加的問題受苦(可通過一英里關閉):
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((@[email protected])/2.0E)) +
(cos(@lat1) * cos(@lat2) * square(sin((@[email protected])/2.0E))))
set @distance =
@radius * (2.0E *asin(case when 1.0E < @a then 1.0E else @a end))
return @distance
end
Vicenty是緩慢的,但精確到內1毫米(我只找到一個javascript小鬼吧):
/*
* 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;
}
感謝您的迴應,但我仍然希望有人在PostGIS中有一個解決方案,或者至少使用它的同一個庫。 – 2008-09-23 19:40:24
我的錯誤 - 我查閱了PostGIS,發現它最初運行在PostgreSQL上。我以爲你可以利用這些邏輯來編寫你需要的東西。 – nathaniel 2008-09-24 16:22:32
PostGIS 1.5使用經緯度和緯度來處理真實的地球距離。它知道lat/long本質上是有角度的,並且有360度的線條
不幸的是,它不在Debian Stable(1.1.6-2)附帶的PostGIS版本中。我猜,是時候開始尋找backports了。 – 2008-09-23 22:58:17