المسافة الحقيقية (الدائرة الكبرى) في PostGIS مع خطوط العرض/الطول SRID؟

Question

أنا أستخدم SRID خطوط العرض/الطول في قاعدة بيانات PostGIS الخاصة بي (-4326).أرغب في العثور على أقرب النقاط إلى نقطة معينة بطريقة فعالة.حاولت القيام ب

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

مما يعطيني نتائج جيدة في الولايات الـ 48 السفلى، ولكن في ألاسكا يعطيني نتائج سيئة.هل هناك طريقة لإجراء حسابات المسافة الحقيقية في PostGIS، أم سأضطر إلى توفير مخزن مؤقت بحجم معقول ثم حساب مسافات الدائرة الكبرى وفرز النتائج في الكود بعد ذلك؟

Answer 1

أنت تبحث عن ST_distance_sphere(point,point) أو st_distance_spheroid(point,point).

يرى:

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

ويشار عادة إلى المسافة الجيوديسية أو الجيوديسية...في حين أن المصطلحين لهما معاني مختلفة قليلاً، إلا أنهما يميلان إلى استخدامهما بالتبادل.

وبدلاً من ذلك، يمكنك عرض البيانات واستخدام الدالة st_distance القياسية...هذا عملي فقط على المسافات القصيرة (باستخدام UTM أو مستوى الحالة) أو إذا كانت جميع المسافات مرتبطة بنقطة واحدة أو نقطتين (إسقاطات متساوية البعد).

Answer 2

يتعامل PostGIS 1.5 مع المسافات الحقيقية للكرة الأرضية باستخدام خطوط الطول والأمتار.إنه يدرك أن خطوط العرض/الطول زاويّة بطبيعتها ولها خط بزاوية 360 درجة

Answer 3

هذا من SQL Server، وأنا أستخدم Haversine لمسافة سريعة يبعث على السخرية قد تعاني من مشكلة ألاسكا الخاصة بك (يمكن أن تكون بعيدة عن ميل):

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 بطيء، ولكنه دقيق في حدود 1 مم (ولقد وجدت فقط نسخة جافا سكريبت منه):

/*
 * 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;
}