根据地球上任意两点的经纬度计算两点间的距离

地球是一个近乎标准的椭球体,它的赤道半径为6378.140千米,极半径为6356.755千米,平均半径6371.004千米。如果我们假设地球是一个完美的球体,那么它的半径就是地球的平均半径,记为R。如果以0度经线为基准,那么根据地球表面任意两点的经纬度就可以计算出这两点间的地表距离(这里忽略地球表面地形对计算带来的误差,仅仅是理论上的估算值)。设第一点A的经纬度为(LonA, LatA),第二点B的经纬度为(LonB, LatB),按照0度经线的基准,东经取经度的正值(Longitude),西经取经度负值(-Longitude),北纬取90-纬度值(90-Latitude),南纬取90+纬度值(90+Latitude),则经过上述处理过后的两点被计为(MLonA, MLatA)和(MLonB, MLatB)。那么根据三角推导,可以得到计算两点距离的如下公式:

C = sin(MLatA)*sin(MLatB)*cos(MLonA-MLonB) + cos(MLatA)*cos(MLatB)

Distance = R*Arccos(C)*Pi/180

这里,R和Distance单位是相同,如果是采用6371.004千米作为半径,那么Distance就是千米为单位,如果要使用其他单位,比如mile,还需要做单位换算,1千米=0.621371192mile 如果仅对经度作正负的处理,而不对纬度作90-Latitude(假设都是北半球,南半球只有澳洲

具有应用意义)的处理,那么公式将是:

C = sin(LatA)*sin(LatB) + cos(LatA)*cos(LatB)*cos(MLonA-MLonB)

Distance = R*Arccos(C)*Pi/180

以上通过简单的三角变换就可以推出。

如果三角函数的输入和输出都采用弧度值,那么公式还可以写作:

C = sin(LatA*Pi/180)*sin(LatB*Pi/180)

cos(LatA*Pi/180)*cos(LatB*Pi/180)*cos((MLonA-MLonB)*Pi/180)

Distance = R*Arccos(C)*Pi/180 +

也就是:

C = sin(LatA/57.2958)*sin(LatB/57.2958) + cos(LatA/57.2958)*cos(LatB/57.2958)*cos((MLonA-MLonB)/57.2958)

Distance = R*Arccos(C) = 6371.004*Arccos(C) kilometer = 0.621371192*6371.004*Arccos(C) mile = 3958.758349716768*Arccos(C) mile

在实际应用当中,一般是通过一个个体的邮政编码来查找该邮政编码对应的地区中心的经纬度,然后再根据这些经纬度来计算彼此的距离,从而估算出某些群体之间的大致距离范围(比如酒店旅客的分布范围-各个旅客的邮政编码对应的经纬度和酒店的经纬度所计算的距离范围-等等),所以,通过邮政编码查询经纬度这样一个数据库是一个很有用的资源。 根据经纬度计算该位置点与附近某一经纬度位置点的位置偏移角

如何计算?

设某一城市有一个城市地标,可以视为城市的中心点 ,其经纬度已知。

分散在城市中心分别

(1)由多层环形交通公路(环线),环城市地标为基准,在不同层的环线上分布着重要的设施。

(2)分布着多个卫星城、主要居民区和学校、基础设施。

现控制无人机分别自城市地标出发,依次探访和航拍环城公路上的重要设施。

无人机配置有GPS导航和电子罗盘。

自城市中心地标处控制无人机起飞。

1 计算任意经纬度点间的距离

用GPS测出两个点的经纬度后,如何计算这两个点之间的距离呢。用的是便携式GPS,只有经纬度显示,无坐标显示和输出。

设两点A、B的经、纬度分别为(jA,wA)(jB,wB),则半径为R的球面上两点间的最短距离(大圆弧)为:

弧AB=R*arccos[sin(wA)sin(wB)+cos(wA)cos(wB)*cos(jA-jB)]

google maps的脚本里代码。

private const double EARTH_RADIUS = 6378.137;

private static double rad(double d)

{

return d * Math.PI / 180.0;

}

public static double GetDistance(double lat1, double lng1, double lat2, double lng2)

{

double radLat1 = rad(lat1);

double radLat2 = rad(lat2);

double a = radLat1 - radLat2;

double b = rad(lng1) - rad(lng2);

double s = 2 * Math.Asin(Math.Sqrt(Math.Pow(Math.Sin(a/2),2) +

Math.Cos(radLat1)*Math.Cos(radLat2)*Math.Pow(Math.Sin(b/2),2)));

s = s * EARTH_RADIUS;

s = Math.Round(s * 10000) / 10000;

return s;

}

;该函数用IDL语言编写,利用GIS中根据两点经纬度计算距离公式

;计算经纬度两点间的距离

Function Length2LatLon,pts0,pts1

radius=6378.14d;地球平均半径,单位公里

; print,pts0,pts1

; length=radius * acos(sin(pts0[1]*!PI/180.0)*sin(pts1[1]*!PI/180.0)+cos(pts0[1]*!PI/180.0)*cos(pts1[1]*!PI/180.0)*cos(pts0[0]*!PI/180.0-pts1[0]*!PI/180.0))

length=radius * acos(sin(pts0[1]*!DTOR)*sin(pts1[1]*!DTOR)+cos(pts0[1]*!DTOR)*cos(pts1[1]*!DTOR)*cos(pts0[0]*!DTOR-pts1[0]*!DTOR))

return,length

end

pro test

pts0=[117.548802d,38.625238d]

pts1=[117.548625d,38.625139d]

print,length2latlon(pts0,pts1)

end

另外在IDL中提供了计算两点距离的方法map_2points可以方便的计算不同单位(Miles,Meters)的距离,还可以通过设置不同的参数来得到不同的距离、角度值,

Syntax:

Result = MAP_2POINTS( lon0, lat0, lon1, lat1 [, DPATH=value | , /METERS | , /MILES | , NPATH=integer{2 or greater} | , /PARAMETERS | , RADIANS=value] [, /RADIUS] [, /RHUMB] ) 具体可以参考IDL 的Help里面讲解

如果要算的距离是椭球面的距离,就难些。

因为大地线的方程是一个微分方程,所以大地线长度是对一个微分方程的开方的积分,无法写成简单的公式.通常是直接用差分代替微分,用求和代替积分算出,这是很严格的.

下面是近似地由大地经纬度求椭球面距离的逼近方法:

如两点在同一纬线上,距离为N*cos(phi)*(lambda2-lambda1)。

如两点在同一经线上,距离为M*dphi从phi1到phi2的积分。

其中N=a/sqrt(d), d=1-e*e*sin(phi)*sin(phi),

M=a*(1-e*e)/[d*sqrt(d)],

其中长短轴a,b与扁率偏心率的关系是:a/b=1-f=sqrt(1-e*e).

把积分区域分成等分n小段,积分可以写成求和,既可算出。这也是严格的。

如两点不在同一经纬线上,且两点距离很近,则小区域的椭球面可以视为平面,于是两点的

距离近似地是上面算出的距离的平方和的开方。如两点距离较远,把两点间的经纬差分成等分n小段,每小段的距离的和就是所求的距离。

google maps的脚本里代码。

private const double EARTH_RADIUS = 6378.137;

private static double rad(double d)

{

return d * Math.PI / 180.0;

}

public static double GetDistance(double lat1, double lng1, double lat2, double lng2)

{

double radLat1 = rad(lat1);

double radLat2 = rad(lat2);

double a = radLat1 - radLat2;

double b = rad(lng1) - rad(lng2);

double s = 2 * Math.Asin(Math.Sqrt(Math.Pow(Math.Sin(a/2),2) +

Math.Cos(radLat1)*Math.Cos(radLat2)*Math.Pow(Math.Sin(b/2),2)));

s = s * EARTH_RADIUS;

s = Math.Round(s * 10000) / 10000;

return s;

}

;该函数用IDL语言编写,利用GIS中根据两点经纬度计算距离公式

;计算经纬度两点间的距离

Function Length2LatLon,pts0,pts1

radius=6378.14d;地球平均半径,单位公里

; print,pts0,pts1

; length=radius * acos(sin(pts0[1]*!PI/180.0)*sin(pts1[1]*!PI/180.0)+cos(pts0[1]*!PI/180.0)*cos(pts1[1]*!PI/180.0)*cos(pts0[0]*!PI/180.0-pts1[0]*!PI/180.0))

length=radius * acos(sin(pts0[1]*!DTOR)*sin(pts1[1]*!DTOR)+cos(pts0[1]*!DTOR)*cos(pts1[1]*!DTOR)*cos(pts0[0]*!DTOR-pts1[0]*!DTOR))

return,length

end

pro test

pts0=[117.548802d,38.625238d]

pts1=[117.548625d,38.625139d]

print,length2latlon(pts0,pts1)

end

另外在IDL中提供了计算两点距离的方法map_2points可以方便的计算不同单位(Miles,Meters)的距离,还可以通过设置不同的参数来得到不同的距离、角度值,

Syntax:

Result = MAP_2POINTS( lon0, lat0, lon1, lat1 [, DPATH=value | , /METERS | , /MILES | , NPATH=integer{2 or greater} | , /PARAMETERS | , RADIANS=value] [, /RADIUS] [, /RHUMB] ) 具体可以参考IDL 的Help里面讲解

如果要算的距离是椭球面的距离,就难些。

因为大地线的方程是一个微分方程,所以大地线长度是对一个微分方程的开方的积分,无法写成简单的公式.通常是直接用差分代替微分,用求和代替积分算出,这是很严格的.

下面是近似地由大地经纬度求椭球面距离的逼近方法:

如两点在同一纬线上,距离为N*cos(phi)*(lambda2-lambda1)。

如两点在同一经线上,距离为M*dphi从phi1到phi2的积分。

其中N=a/sqrt(d), d=1-e*e*sin(phi)*sin(phi),

M=a*(1-e*e)/[d*sqrt(d)],

其中长短轴a,b与扁率偏心率的关系是:a/b=1-f=sqrt(1-e*e).

把积分区域分成等分n小段,积分可以写成求和,既可算出。这也是严格的。

如两点不在同一经纬线上,且两点距离很近,则小区域的椭球面可以视为平面,于是两点的距离近似地是上面算出的距离的平方和的开方。如两点距离较远,把两点间的经纬差分成等分n小段,每小段的距离的和就是所求的距离。