我有一组地点的纬度和经度。
- 如何找到从该组中一个位置到另一个位置的距离?
- 是否有公式?
Haversine公式假设地球是一个球体,但是地球的形状更加复杂。一个扁球体模型会给出更好的结果。
如果需要这样的精度,你应该使用Vincenty反算法。详情请参见http://en.wikipedia.org/wiki/Vincenty's_formulae。使用该算法,对于扁球体模型可以得到0.5mm的精度。
没有完美的公式,因为地球的真实形状太过复杂而无法用公式表达。此外,由于气候事件(参见http://www.nasa.gov/centers/goddard/earthandsun/earthshape.html)和地球自转,地球的形状也会随着时间发生变化。
你还应该注意,上述方法不考虑高度,并假设扁球体为海平面。
编辑 10-Jul-2010:我发现有少数情况下Vincenty反算法无法收敛到声明的精度。更好的想法是使用GeographicLib(请参见http://sourceforge.net/projects/geographiclib/),它也更准确。
这里有一个链接: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
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
。
您是否正在寻找
Haversine公式是一种在导航中非常重要的方程,可以根据两点的经纬度计算它们在球体上的大圆距离。它是球面三角学中更通用公式——haversines定理的特例,该定理涉及球面“三角形”的边和角。
这里有一个示例,可以根据给定的IP查找位置/附近位置的经纬度:
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;
}
它返回距离的公里数
如果你要测量小于(也许)1度纬度/经度的距离,需要一个非常高效的近似值,并且愿意接受比Haversine公式更多的不准确性,请考虑以下两个替代方案:
(1) 计算距离中的“极坐标平面公式”:
a = pi/2 - lat1
b = pi/2 - lat2
c = sqrt( a^2 + b^2 - 2 * a * b * cos(lon2 - lon1) )
d = R * c
(2) 勾股定理根据纬度调整,如Ewan Todd的SO帖子所示:
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)
内部用average(lat0,lat1)
替换了lat0
。
#2没有明确数值是以度、弧度还是公里为单位,您需要编写一些转换代码。请参见本帖子底部的完整代码。
#1旨在即使在极地附近也能正常工作,但是如果您正在测量两个端点位于“相对”极点的距离(经度相差超过90度?),则建议使用Haversine,即使对于小距离也是如此。
我没有彻底测量这些方法的误差,因此您应该为您的应用程序选择代表性点,并将结果与某些高质量库进行比较,以确定准确性是否可接受。对于少于几千米的距离,我的直觉是这些距离在正确测量的情况下的误差不到1%。
一种获得高性能的替代方法(适用于以下情况):
如果您有一个大型的静态点集,其经度/纬度在一到两度之间,并且您将从少量动态(移动)点计算距离,请考虑将您的静态点转换为包含UTM区域(或任何其他本地笛卡尔坐标系),然后在该笛卡尔坐标系中进行所有计算。
笛卡尔坐标系=平面地球=勾股定理适用,因此distance = sqrt(dx^2 + dy^2)
。
然后,准确转换少数移动点的成本是可以轻松承受的。
警告 #1(Polar):对于小于0.1(?)米的距离可能非常不准确。即使使用双精度数学,以下坐标的真实距离约为0.005米,但我的Polar算法实现将其给出为“零”:
输入:
lon1Xdeg 16.6564465477996 double
lat1Ydeg 57.7760262271983 double
lon2Xdeg 16.6564466358281 double
lat2Ydeg 57.776026248554 double
results:
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部分”或1%的误差。
比较极坐标误差与哈弗赛恩误差在距离小于2000米的范围内,以查看这个更简单的公式有多糟糕。到目前为止,我见过的最糟糕的情况是极坐标误差每1000份中有51份,而哈弗赛恩误差每1000份中只有4份。在大约58度纬度。
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
这些输入参数的误差在极坐标中较大;可能是我的代码有误,也可能是我正在运行的Cos函数有误,或者我不得不建议不使用极坐标,即使大多数极坐标测量结果比这个更接近。
另一方面,即使进行了* cos(latitude)
的调整,毕达哥拉斯定理的误差也比距离增长更快(最大误差/距离的比率随着距离的增加而增加),因此您需要仔细考虑您将要测量的最大距离和可接受的误差。此外,不建议使用毕达哥拉斯定理来比较两个几乎相等的距离以确定哪个更短,因为误差在不同的方向上是不同的(未显示证据)。
最坏情况下的测量结果为:errorRatio = Abs(error) / distance
(瑞典;最高达2000米):
t_maxHaversineErrorRatio 0.00351012021578681 double
t_maxPolarErrorRatio 66.0825360597085 double
t_maxPythagoreanErrorRatio 0.00350976281416454 double
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://dev59.com/BHVD5IYBdhLWcg3wTJrF#19772119
// 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://dev59.com/jkrSa4cB1Zd3GeqPThpg#1664836
// Refined by ToolmakerSteve - https://dev59.com/cnM_5IYBdhLWcg3wWRkF#53468745
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;
}
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
我已经完成使用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