【发布时间】：2010-11-28 01:10:37
【问题描述】：
我有一组位置的纬度和经度。
- 如何找到从集合中的一个位置到另一个位置的距离？
- 有公式吗？
【问题讨论】：
【发布时间】：2010-11-28 01:10:37
【问题描述】：
你在找吗
harsine 公式是一个方程 在导航中很重要，给予 两个之间的大圆距离 从他们的球体上的点 经度和纬度。它是一个 更一般公式的特例 在球面三角学中， hasrsines, 关联双方和 球面“三角形”的角度。
【讨论】：
看看这个.. 还有一个 javascript 示例。
【讨论】：
-
仅链接的答案会随着时间的推移而失效，因此用处不大。
【讨论】：
-
仅链接的答案会随着时间的推移而失效，因此用处不大。
这是一个：http://www.movable-type.co.uk/scripts/latlong.html
使用Haversine公式：
R = earth’s radius (mean radius = 6,371km)
Δlat = lat2− lat1
Δlong = long2− long1
a = sin²(Δlat/2) + cos(lat1).cos(lat2).sin²(Δlong/2)
c = 2.atan2(√a, √(1−a))
d = R.c
【讨论】：
Haversine 公式假设地球是球形的。然而，耳朵的形状更复杂。扁球体模型会得到更好的结果。
如果需要这样的精度，最好使用Vincenty 逆公式。 有关详细信息，请参阅http://en.wikipedia.org/wiki/Vincenty's_formulae。使用它，您可以获得 0.5mm 的球体模型精度。
没有完美的公式，因为地球的真实形状太复杂，无法用公式来表达。此外，地球的形状会因气候事件而发生变化（见http://www.nasa.gov/centers/goddard/earthandsun/earthshape.html），也会因地球自转而随时间变化。
您还应该注意，上述方法没有考虑海拔高度，并假设海平面为扁球体。
2010 年 7 月 10 日编辑： 我发现在极少数情况下 Vincenty 逆公式不会收敛到声明的准确性。更好的办法是使用 GeographicLib（参见 http://sourceforge.net/projects/geographiclib/），它也更准确。
【讨论】：
-
+1 这让很多人在以前的雇主那里措手不及。
-
确实如此。当值不能超过几米时，这个问题就会变得复杂得多。
-
I've provided C# code 实现了扁球体（我不确定它是否是 VIncenty）、Haversine 和两个经常提到的长达 1 公里左右的高性能近似值。考虑到扁球体是“正确的”，在瑞典的一个地区（纬度 ~58 度）显示了其他公式的测量误差。
应用Haversine 公式求距离。请参阅下面的 C# 代码以查找 2 个坐标之间的距离。如果您想说要查找某个范围内的商店列表，则更好的是，您可以在 SQL 中应用 WHERE 子句或在 C# 中应用 LINQ 过滤器。
这里的公式是以公里为单位的，您必须更改相关数字，它才能以英里为单位。
例如：将 6371.392896 转换为英里。
DECLARE @radiusInKm AS FLOAT
DECLARE @lat2Compare AS FLOAT
DECLARE @long2Compare AS FLOAT
SET @radiusInKm = 5.000
SET @lat2Compare = insert_your_lat_to_compare_here
SET @long2Compare = insert_you_long_to_compare_here
SELECT * FROM insert_your_table_here WITH(NOLOCK)
WHERE (6371.392896*2*ATN2(SQRT((sin((radians(GeoLatitude - @lat2Compare)) / 2) * sin((radians(GeoLatitude - @lat2Compare)) / 2)) + (cos(radians(GeoLatitude)) * cos(radians(@lat2Compare)) * sin(radians(GeoLongitude - @long2Compare)/2) * sin(radians(GeoLongitude - @long2Compare)/2)))
, SQRT(1-((sin((radians(GeoLatitude - @lat2Compare)) / 2) * sin((radians(GeoLatitude - @lat2Compare)) / 2)) + (cos(radians(GeoLatitude)) * cos(radians(@lat2Compare)) * sin(radians(GeoLongitude - @long2Compare)/2) * sin(radians(GeoLongitude - @long2Compare)/2)))
))) <= @radiusInKm
如果您想在 C# 中执行 Haversine 公式，
double resultDistance = 0.0;
double avgRadiusOfEarth = 6371.392896; //Radius of the earth differ, I'm taking the average.
//Haversine formula
//distance = R * 2 * aTan2 ( square root of A, square root of 1 - A )
// where A = sinus squared (difference in latitude / 2) + (cosine of latitude 1 * cosine of latitude 2 * sinus squared (difference in longitude / 2))
// and R = the circumference of the earth
double differenceInLat = DegreeToRadian(currentLatitude - latitudeToCompare);
double differenceInLong = DegreeToRadian(currentLongitude - longtitudeToCompare);
double aInnerFormula = Math.Cos(DegreeToRadian(currentLatitude)) * Math.Cos(DegreeToRadian(latitudeToCompare)) * Math.Sin(differenceInLong / 2) * Math.Sin(differenceInLong / 2);
double aFormula = (Math.Sin((differenceInLat) / 2) * Math.Sin((differenceInLat) / 2)) + (aInnerFormula);
resultDistance = avgRadiusOfEarth * 2 * Math.Atan2(Math.Sqrt(aFormula), Math.Sqrt(1 - aFormula));
DegreesToRadian是我自定义创建的一个函数，它是"Math.PI * angle / 180.0的简单1行
【讨论】：
只需使用距离公式Sqrt( (x2-x1)^2 + (y2-y1)^2 )
【讨论】：
-
不要这样做。它仅适用于赤道附近的小距离。远离赤道，即使是小距离也是错误的；在高纬度地区，经度 1 度的变化仅
* cos(latitude)与纬度 1 度变化的距离一样大。
这是一个通过给定 IP 查找位置/附近位置到 long/lat：
http://jsfiddle.net/bassta/zrgd9qc3/2/
这是我用来计算直线距离的函数：
function distance(lat1, lng1, lat2, lng2) {
var radlat1 = Math.PI * lat1 / 180;
var radlat2 = Math.PI * lat2 / 180;
var radlon1 = Math.PI * lng1 / 180;
var radlon2 = Math.PI * lng2 / 180;
var theta = lng1 - lng2;
var radtheta = Math.PI * theta / 180;
var dist = Math.sin(radlat1) * Math.sin(radlat2) + Math.cos(radlat1) * Math.cos(radlat2) * Math.cos(radtheta);
dist = Math.acos(dist);
dist = dist * 180 / Math.PI;
dist = dist * 60 * 1.1515;
//Get in in kilometers
dist = dist * 1.609344;
return dist;
}
它以公里为单位返回距离
【讨论】：
-
请说出这是基于什么技术。哈弗辛？文森蒂？其他一些近似？并提供指向原始来源的链接，以证明您的代码合理。
-
您好，代码是 正交距离公式 的 JavaScript 实现，用于查找完美球面上两点之间的距离。我知道地球不是一个完美的球体，但就我而言，我不想要更高的精度。
-
谢谢，这有助于我了解您所拥有的。我看到several different formulas for great-circle distance。您是否知道您的代码是源自 Haversine、Vicenty 或弦长公式，还是球体的其他数学表示？
-
这个答案没有说明正在使用什么距离公式，这会对准确性产生重要影响。
在此页面上，您可以查看 Android Location 类中如何计算位置距离的完整代码和公式
android/location/Location.java
编辑：根据@Richard 的提示，我将链接函数的代码放入我的答案中，以避免链接无效：
private static void computeDistanceAndBearing(double lat1, double lon1,
double lat2, double lon2, BearingDistanceCache results) {
// Based on http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
// using the "Inverse Formula" (section 4)
int MAXITERS = 20;
// Convert lat/long to radians
lat1 *= Math.PI / 180.0;
lat2 *= Math.PI / 180.0;
lon1 *= Math.PI / 180.0;
lon2 *= Math.PI / 180.0;
double a = 6378137.0; // WGS84 major axis
double b = 6356752.3142; // WGS84 semi-major axis
double f = (a - b) / a;
double aSqMinusBSqOverBSq = (a * a - b * b) / (b * b);
double L = lon2 - lon1;
double A = 0.0;
double U1 = Math.atan((1.0 - f) * Math.tan(lat1));
double U2 = Math.atan((1.0 - f) * Math.tan(lat2));
double cosU1 = Math.cos(U1);
double cosU2 = Math.cos(U2);
double sinU1 = Math.sin(U1);
double sinU2 = Math.sin(U2);
double cosU1cosU2 = cosU1 * cosU2;
double sinU1sinU2 = sinU1 * sinU2;
double sigma = 0.0;
double deltaSigma = 0.0;
double cosSqAlpha = 0.0;
double cos2SM = 0.0;
double cosSigma = 0.0;
double sinSigma = 0.0;
double cosLambda = 0.0;
double sinLambda = 0.0;
double lambda = L; // initial guess
for (int iter = 0; iter < MAXITERS; iter++) {
double lambdaOrig = lambda;
cosLambda = Math.cos(lambda);
sinLambda = Math.sin(lambda);
double t1 = cosU2 * sinLambda;
double t2 = cosU1 * sinU2 - sinU1 * cosU2 * cosLambda;
double sinSqSigma = t1 * t1 + t2 * t2; // (14)
sinSigma = Math.sqrt(sinSqSigma);
cosSigma = sinU1sinU2 + cosU1cosU2 * cosLambda; // (15)
sigma = Math.atan2(sinSigma, cosSigma); // (16)
double sinAlpha = (sinSigma == 0) ? 0.0 :
cosU1cosU2 * sinLambda / sinSigma; // (17)
cosSqAlpha = 1.0 - sinAlpha * sinAlpha;
cos2SM = (cosSqAlpha == 0) ? 0.0 :
cosSigma - 2.0 * sinU1sinU2 / cosSqAlpha; // (18)
double uSquared = cosSqAlpha * aSqMinusBSqOverBSq; // defn
A = 1 + (uSquared / 16384.0) * // (3)
(4096.0 + uSquared *
(-768 + uSquared * (320.0 - 175.0 * uSquared)));
double B = (uSquared / 1024.0) * // (4)
(256.0 + uSquared *
(-128.0 + uSquared * (74.0 - 47.0 * uSquared)));
double C = (f / 16.0) *
cosSqAlpha *
(4.0 + f * (4.0 - 3.0 * cosSqAlpha)); // (10)
double cos2SMSq = cos2SM * cos2SM;
deltaSigma = B * sinSigma * // (6)
(cos2SM + (B / 4.0) *
(cosSigma * (-1.0 + 2.0 * cos2SMSq) -
(B / 6.0) * cos2SM *
(-3.0 + 4.0 * sinSigma * sinSigma) *
(-3.0 + 4.0 * cos2SMSq)));
lambda = L +
(1.0 - C) * f * sinAlpha *
(sigma + C * sinSigma *
(cos2SM + C * cosSigma *
(-1.0 + 2.0 * cos2SM * cos2SM))); // (11)
double delta = (lambda - lambdaOrig) / lambda;
if (Math.abs(delta) < 1.0e-12) {
break;
}
}
float distance = (float) (b * A * (sigma - deltaSigma));
results.mDistance = distance;
float initialBearing = (float) Math.atan2(cosU2 * sinLambda,
cosU1 * sinU2 - sinU1 * cosU2 * cosLambda);
initialBearing *= 180.0 / Math.PI;
results.mInitialBearing = initialBearing;
float finalBearing = (float) Math.atan2(cosU1 * sinLambda,
-sinU1 * cosU2 + cosU1 * sinU2 * cosLambda);
finalBearing *= 180.0 / Math.PI;
results.mFinalBearing = finalBearing;
results.mLat1 = lat1;
results.mLat2 = lat2;
results.mLon1 = lon1;
results.mLon2 = lon2;
}
【讨论】：
-
仅链接的答案会随着时间的推移而失效，因此用处不大。
-
你说得对@Richard 我更新了链接并将链接功能代码添加到我的答案中。
-
谢谢！对文档进行更完整的引用会很有用，这样人们就可以通过其他方式找到它。
以下是包含前面答案中讨论的三个公式的模块（以 f90 编码）。您可以将此模块放在程序的顶部 （在 PROGRAM MAIN 之前）或单独编译并在编译期间包含模块目录。以下模块包含三个公式。前两个是基于地球是球形假设的大圆距离。
module spherical_dists
contains
subroutine great_circle_distance(lon1,lat1,lon2,lat2,dist)
!https://en.wikipedia.org/wiki/Great-circle_distance
! It takes lon, lats of two points on an assumed spherical earth and
! calculates the distance between them along the great circle connecting the two points
implicit none
real,intent(in)::lon1,lon2,lat1,lat2
real,intent(out)::dist
real,parameter::pi=3.141592,mean_earth_radius=6371.0088
real::lonr1,lonr2,latr1,latr2
real::delangl,dellon
lonr1=lon1*(pi/180.);lonr2=lon2*(pi/180.)
latr1=lat1*(pi/180.);latr2=lat2*(pi/180.)
dellon=lonr2-lonr1
delangl=acos(sin(latr1)*sin(latr2)+cos(latr1)*cos(latr2)*cos(dellon))
dist=delangl*mean_earth_radius
end subroutine
subroutine haversine_formula(lon1,lat1,lon2,lat2,dist)
! https://en.wikipedia.org/wiki/Haversine_formula
! This is similar above but numerically better conditioned for small distances
implicit none
real,intent(in)::lon1,lon2,lat1,lat2
!lon, lats of two points
real,intent(out)::dist
real,parameter::pi=3.141592,mean_earth_radius=6371.0088
real::lonr1,lonr2,latr1,latr2
real::delangl,dellon,dellat,a
! degrees are converted to radians
lonr1=lon1*(pi/180.);lonr2=lon2*(pi/180.)
latr1=lat1*(pi/180.);latr2=lat2*(pi/180.)
dellon=lonr2-lonr1 ! These dels simplify the haversine formula
dellat=latr2-latr1
! The actual haversine formula
a=(sin(dellat/2))**2+cos(latr1)*cos(latr2)*(sin(dellon/2))**2
delangl=2*asin(sqrt(a)) !2*asin(sqrt(a))
dist=delangl*mean_earth_radius
end subroutine
subroutine vincenty_formula(lon1,lat1,lon2,lat2,dist)
!https://en.wikipedia.org/wiki/Vincenty%27s_formulae
!It's a better approximation over previous two, since it considers earth to in oblate spheroid, which better approximates the shape of the earth
implicit none
real,intent(in)::lon1,lon2,lat1,lat2
real,intent(out)::dist
real,parameter::pi=3.141592,mean_earth_radius=6371.0088
real::lonr1,lonr2,latr1,latr2
real::delangl,dellon,nom,denom
lonr1=lon1*(pi/180.);lonr2=lon2*(pi/180.)
latr1=lat1*(pi/180.);latr2=lat2*(pi/180.)
dellon=lonr2-lonr1
nom=sqrt((cos(latr2)*sin(dellon))**2. + (cos(latr1)*sin(latr2)-sin(latr1)*cos(latr2)*cos(dellon))**2.)
denom=sin(latr1)*sin(latr2)+cos(latr1)*cos(latr2)*cos(dellon)
delangl=atan2(nom,denom)
dist=delangl*mean_earth_radius
end subroutine
end module
【讨论】：
-
此代码缺少格式。更重要的是，它缺乏必要的 cmets 来解释它的作用和使用的公式。
-
我编辑了答案以添加更多详细信息并格式化代码。
我已完成使用 SQL 查询
select *, (acos(sin(input_lat* 0.01745329)*sin(lattitude *0.01745329) + cos(input_lat *0.01745329)*cos(lattitude *0.01745329)*cos((input_long -longitude)*0.01745329))* 57.29577951 )* 69.16 As D from table_name
【讨论】：
如果您测量的距离小于（可能）1 度的纬度/经度变化，正在寻找性能非常高的近似值，并且愿意接受更多不准确除了Haversine 公式，请考虑以下两种选择：
(1) 来自Computing Distances的“极坐标平地公式”：
a = pi/2 - lat1
b = pi/2 - lat2
c = sqrt( a^2 + b^2 - 2 * a * b * cos(lon2 - lon1) )
d = R * c
(2) 勾股定理根据纬度调整，如in Ewan Todd's SO post:
d_ew = (long1 - long0) * cos(average(lat0, lat1))
d_ns = (lat1 - lat0)
d = sqrt(d_ew * d_ew + d_ns * d_ns)
注意事项：
与 Ewan 的帖子相比，我在 cos( lat0 ) 中用 lat0 替换了 average(lat0, lat1)。
#2 对值是度数、弧度还是公里数不明确；您还需要一些转换代码。在这篇文章的底部查看我的完整代码。
#1 的设计即使在两极附近也能正常工作，但如果您测量的距离其端点位于两极的“相反”两侧（经度相差超过 90 度？），建议使用 Haversine，即使适合小距离。
我还没有彻底测量这些方法的误差，所以你应该为你的应用取代表点，并将结果与一些高质量的库进行比较，以确定准确性是否可以接受。对于不到几公里的距离，我的直觉是这些距离正确测量的误差在 1% 以内。
另一种获得高性能的方法（如适用）：
如果您有大量 静态 点，在经度/纬度一到两度范围内，那么您将计算与少量 动态 点的距离（移动）点，请考虑将您的 静态 点转换为包含 UTM 区域（或任何其他本地笛卡尔坐标系），然后在该笛卡尔坐标系中进行所有数学运算。
笛卡尔=平坦地球=勾股定理适用，所以distance = sqrt(dx^2 + dy^2)。
那么准确将少数移动点转换为 UTM 的成本很容易负担。
#1（极地）的警告：对于小于 0.1 (?) 米的距离可能是非常错误的。即使使用双精度数学，以下坐标（其真实距离约为 0.005 米）在我的 Polar 算法实现中被指定为“零”：
输入：
lon1Xdeg 16.6564465477996 double
lat1Ydeg 57.7760262271983 double
lon2Xdeg 16.6564466358281 double
lat2Ydeg 57.776026248554 double
结果：
Oblate spheroid formula:
0.00575254911118364 double
Haversine:
0.00573422966122257 double
Polar:
0
这是由于 u 和 v 两个因素完全相互抵消：
u 0.632619944868587 double
v -0.632619944868587 double
在另一种情况下，当扁球体答案为0.002887 m 时，它给出的距离为0.067129 m。问题是cos(lon2 - lon1) 与1 太接近，所以cos 函数返回的正是1。
除了测量亚米距离之外，我发现的最大误差（与扁球体公式相比）是迄今为止我输入的有限小距离数据：
maxHaversineErrorRatio 0.00350976281908381 double
maxPolarErrorRatio 0.0510789996931342 double
其中“1”表示答案中有 100% 的错误；例如当它返回“0”时，这是一个错误“1”（从上面的“maxPolar”中排除）。所以“0.01”将是“100 分之一”或 1% 的错误。
在小于 2000 米的距离内比较 Polar 误差与 Haversine 误差，看看这个更简单的公式有多糟糕。到目前为止，我见过的最差的是 Polar 的每 1000 份 51 份，而 Haversine 的每 1000 份中的 4 份。大约在纬度 58 度。
现在实现了“带纬度调整的毕达哥拉斯”。
对于 我最初认为 Polar 问题仅在 但是下面显示的结果非常令人不安。
随着距离接近于零，毕达哥拉斯/纬度接近半正弦。 例如这个测量值 ~ 217 米：
lon1Xdeg 16.6531667510102 double
lat1Ydeg 57.7751705615804 double
lon2Xdeg 16.6564468739869 double
lat2Ydeg 57.7760263007586 double
oblate 217.201200413731
haversine 216.518428601051
polar 226.128616011973
pythag-cos 216.518428631907
havErrRatio 0.00314349925958048
polErrRatio 0.041102054598393
pycErrRatio 0.00314349911751603
Polar 对这些输入有更严重的错误；要么我的代码有错误，要么我运行的 Cos 函数有错误，或者我不得不建议不要使用 Polar，即使大多数 Polar 测量值都比这更接近。
OTOH，毕达哥拉斯，即使使用* cos(latitude) 调整，误差也会比距离增加得更快（max_error/distance 的比率随着距离的增加而增加），因此您需要仔细考虑要测量的最大距离，和可接受的错误。此外，不建议使用毕达哥拉斯比较两个几乎相等的距离，以确定哪个更短，因为不同方向的误差不同（证据未显示）。
最坏情况测量，errorRatio = Abs(error) / distance（瑞典；高达 2000 m）：
t_maxHaversineErrorRatio 0.00351012021578681 double
t_maxPolarErrorRatio 66.0825360597085 double
t_maxPythagoreanErrorRatio 0.00350976281416454 double
如前所述，极端极坐标误差适用于亚米距离，它可能报告零而不是 6 cm，或者报告超过 0.5 m 的距离为 1 cm（因此“66 x”最坏情况显示在t_maxPolarErrorRatio），但在较大距离处也有一些较差的结果。 [需要使用已知高度准确的余弦函数再次测试。]
在 Moto E4 上运行的 Xamarin.Android 中的 C# 代码中进行的测量。
C#代码：
// x=longitude, y= latitude. oblate spheroid formula. TODO: From where?
public static double calculateDistanceDD_AED( double lon1Xdeg, double lat1Ydeg, double lon2Xdeg, double lat2Ydeg )
{
double c_dblEarthRadius = 6378.135; // km
double c_dblFlattening = 1.0 / 298.257223563; // WGS84 inverse
// flattening
// Q: Why "-" for longitudes??
double p1x = -degreesToRadians( lon1Xdeg );
double p1y = degreesToRadians( lat1Ydeg );
double p2x = -degreesToRadians( lon2Xdeg );
double p2y = degreesToRadians( lat2Ydeg );
double F = (p1y + p2y) / 2;
double G = (p1y - p2y) / 2;
double L = (p1x - p2x) / 2;
double sing = Math.Sin( G );
double cosl = Math.Cos( L );
double cosf = Math.Cos( F );
double sinl = Math.Sin( L );
double sinf = Math.Sin( F );
double cosg = Math.Cos( G );
double S = sing * sing * cosl * cosl + cosf * cosf * sinl * sinl;
double C = cosg * cosg * cosl * cosl + sinf * sinf * sinl * sinl;
double W = Math.Atan2( Math.Sqrt( S ), Math.Sqrt( C ) );
if (W == 0.0)
return 0.0;
double R = Math.Sqrt( (S * C) ) / W;
double H1 = (3 * R - 1.0) / (2.0 * C);
double H2 = (3 * R + 1.0) / (2.0 * S);
double D = 2 * W * c_dblEarthRadius;
// Apply flattening factor
D = D * (1.0 + c_dblFlattening * H1 * sinf * sinf * cosg * cosg - c_dblFlattening * H2 * cosf * cosf * sing * sing);
// Transform to meters
D = D * 1000.0;
// tmstest
if (true)
{
// Compare Haversine.
double haversine = HaversineApproxDistanceGeo( lon1Xdeg, lat1Ydeg, lon2Xdeg, lat2Ydeg );
double error = haversine - D;
double absError = Math.Abs( error );
double errorRatio = absError / D;
if (errorRatio > t_maxHaversineErrorRatio)
{
if (errorRatio > t_maxHaversineErrorRatio * 1.1)
Helper.test();
t_maxHaversineErrorRatio = errorRatio;
}
// Compare Polar Coordinate Flat Earth.
double polarDistanceGeo = ApproxDistanceGeo_Polar( lon1Xdeg, lat1Ydeg, lon2Xdeg, lat2Ydeg, D );
double error2 = polarDistanceGeo - D;
double absError2 = Math.Abs( error2 );
double errorRatio2 = absError2 / D;
if (errorRatio2 > t_maxPolarErrorRatio)
{
if (polarDistanceGeo > 0)
{
if (errorRatio2 > t_maxPolarErrorRatio * 1.1)
Helper.test();
t_maxPolarErrorRatio = errorRatio2;
}
else
Helper.dubious();
}
// Compare Pythagorean Theorem with Latitude Adjustment.
double pythagoreanDistanceGeo = ApproxDistanceGeo_PythagoreanCosLatitude( lon1Xdeg, lat1Ydeg, lon2Xdeg, lat2Ydeg, D );
double error3 = pythagoreanDistanceGeo - D;
double absError3 = Math.Abs( error3 );
double errorRatio3 = absError3 / D;
if (errorRatio3 > t_maxPythagoreanErrorRatio)
{
if (D < 2000)
{
if (errorRatio3 > t_maxPythagoreanErrorRatio * 1.05)
Helper.test();
t_maxPythagoreanErrorRatio = errorRatio3;
}
}
}
return D;
}
// As a fraction of the distance.
private static double t_maxHaversineErrorRatio, t_maxPolarErrorRatio, t_maxPythagoreanErrorRatio;
// Average of equatorial and polar radii (meters).
public const double EarthAvgRadius = 6371000;
public const double EarthAvgCircumference = EarthAvgRadius * 2 * PI;
// CAUTION: This is an average of great circles; won't be the actual distance of any longitude or latitude degree.
public const double EarthAvgMeterPerGreatCircleDegree = EarthAvgCircumference / 360;
// Haversine formula (assumes Earth is sphere).
// "deg" = degrees.
// Perhaps based on Haversine Formula in https://cs.nyu.edu/visual/home/proj/tiger/gisfaq.html
public static double HaversineApproxDistanceGeo(double lon1Xdeg, double lat1Ydeg, double lon2Xdeg, double lat2Ydeg)
{
double lon1 = degreesToRadians( lon1Xdeg );
double lat1 = degreesToRadians( lat1Ydeg );
double lon2 = degreesToRadians( lon2Xdeg );
double lat2 = degreesToRadians( lat2Ydeg );
double dlon = lon2 - lon1;
double dlat = lat2 - lat1;
double sinDLat2 = Sin( dlat / 2 );
double sinDLon2 = Sin( dlon / 2 );
double a = sinDLat2 * sinDLat2 + Cos( lat1 ) * Cos( lat2 ) * sinDLon2 * sinDLon2;
double c = 2 * Atan2( Sqrt( a ), Sqrt( 1 - a ) );
double d = EarthAvgRadius * c;
return d;
}
// From https://stackoverflow.com/a/19772119/199364
// Based on Polar Coordinate Flat Earth in https://cs.nyu.edu/visual/home/proj/tiger/gisfaq.html
public static double ApproxDistanceGeo_Polar( double lon1deg, double lat1deg, double lon2deg, double lat2deg, double D = 0 )
{
double approxUnitDistSq = ApproxUnitDistSq_Polar(lon1deg, lat1deg, lon2deg, lat2deg, D);
double c = Sqrt( approxUnitDistSq );
return EarthAvgRadius * c;
}
// Might be useful to avoid taking Sqrt, when comparing to some threshold.
// Threshold would have to be adjusted to match: Power(threshold / EarthAvgRadius, 2)
private static double ApproxUnitDistSq_Polar(double lon1deg, double lat1deg, double lon2deg, double lat2deg, double D = 0 )
{
const double HalfPi = PI / 2; //1.5707963267949;
double lon1 = degreesToRadians(lon1deg);
double lat1 = degreesToRadians(lat1deg);
double lon2 = degreesToRadians(lon2deg);
double lat2 = degreesToRadians(lat2deg);
double a = HalfPi - lat1;
double b = HalfPi - lat2;
double u = a * a + b * b;
double dlon21 = lon2 - lon1;
double cosDeltaLon = Cos( dlon21 );
double v = -2 * a * b * cosDeltaLon;
// TBD: Is "Abs" necessary? That is, is "u + v" ever negative?
// (I think not; "v" looks like a secondary term. Though might be round-off issue near zero when a~=b.)
double approxUnitDistSq = Abs(u + v);
//if (approxUnitDistSq.nearlyEquals(0, 1E-16))
// Helper.dubious();
//else if (D > 0)
//{
// double dba = b - a;
// double unitD = D / EarthAvgRadius;
// double unitDSq = unitD * unitD;
// if (approxUnitDistSq > 2 * unitDSq)
// Helper.dubious();
// else if (approxUnitDistSq * 2 < unitDSq)
// Helper.dubious();
//}
return approxUnitDistSq;
}
// Pythagorean Theorem with Latitude Adjustment - from Ewan Todd - https://stackoverflow.com/a/1664836/199364
// Refined by ToolmakerSteve - https://stackoverflow.com/a/53468745/199364
public static double ApproxDistanceGeo_PythagoreanCosLatitude( double lon1deg, double lat1deg, double lon2deg, double lat2deg, double D = 0 )
{
double approxDegreesSq = ApproxDegreesSq_PythagoreanCosLatitude( lon1deg, lat1deg, lon2deg, lat2deg );
// approximate degrees on the great circle between the points.
double d_degrees = Sqrt( approxDegreesSq );
return d_degrees * EarthAvgMeterPerGreatCircleDegree;
}
public static double ApproxDegreesSq_PythagoreanCosLatitude( double lon1deg, double lat1deg, double lon2deg, double lat2deg )
{
double avgLatDeg = average( lat1deg , lat2deg );
double avgLat = degreesToRadians( avgLatDeg );
double d_ew = (lon2deg - lon1deg) * Cos( avgLat );
double d_ns = (lat2deg - lat1deg);
double approxDegreesSq = d_ew * d_ew + d_ns * d_ns;
return approxDegreesSq;
}
【讨论】：