在PostGIS中使用纬度/经度SRID计算真实（大圆）距离？

Question

我正在使用PostGIS数据库中的经纬度SRID（-4326）。我希望以高效的方式找到给定点最近的点。我尝试过

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

在美国的48个州，我的代码可以给出较好的结果，但是在阿拉斯加州，我的代码就会出现错误。请问是否有一种方法可以在PostGIS中进行真正的距离计算，还是我需要给出合理大小的缓冲区，然后在代码中计算大圆距离，并在之后对结果进行排序？

- Paul Tomblin

Answer 1

你正在寻找ST_distance_sphere（point，point）或st_distance_spheroid（point，point）。

请参阅：

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或州平面）或如果所有距离都相对于一个或两个点（等距投影）。

- James Schek

---

很不幸，Debian Stable 自带的 PostGIS 版本（1.1.6-2）中没有这个功能。看来是时候开始寻找后向移植了。 - Paul Tomblin

Answer 2

PostGIS 1.5使用经纬度和米来处理真实的地球距离。它知道经纬度是角度性质的，并且有一个360度的线。

- TheSteve0

---

怎么回事？对我来说行不通。 - Bartek Banachewicz

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毫米以内的精确度（我只找到了一个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;
}

- nathaniel

---

谢谢回复，但我仍然希望有人能够在PostGIS中找到解决方案，或者至少使用与其相同的库。 - Paul Tomblin

我的错误——我查了一下PostGIS，发现它最初是在PostgreSQL上运行的。我以为你可以利用其中的一些逻辑来编写你需要的内容。 - nathaniel