基于WGS-84坐标系的GPS坐标下两点间距离计算 VB.net实现

 

精度更高的一种计算方式.

 

WGS－84坐标系是一种国际上采用的地心坐标系。坐标原点为地球质心，其地心空间直角坐标系的Z轴指向国际时间局（BIH）1984.0定义的协议地极（CTP）方向，X轴指向BIH1984.0的协议子午面和CTP赤道的交点，Y轴与Z轴、X轴垂直构成右手坐标系，称为1984年世界大地坐标系。这是一个国际协议地球参考系统（ITRS），是目前国际上统一采用的大地坐标系。

 

 

''' <summary> ''' 通过经纬度计算两点间距离 ''' </summary> ''' <param name="lat1">出发点纬度</param> ''' <param name="lon1">出发点经度</param> ''' <param name="lat2">目标点纬度</param> ''' <param name="lon2">目标点经度</param> ''' <returns>距离 单位()</returns> ''' <remarks>采用WGS-84 椭球模型</remarks> Public Function GPSDistance(ByVal lat1 As Double, ByVal lon1 As Double, ByVal lat2 As Double, ByVal lon2 As Double) As Double Dim a As Double = 6378137.0, b As Double = 6356752.3142, f As Double = 1 / 298.257223563 Dim rlat1 As Double = GetRad(lat1) Dim rlon1 As Double = GetRad(lon1) Dim rlat2 As Double = GetRad(lat2) Dim rlon2 As Double = GetRad(lon2) Dim L As Double = GetRad(lon2 - lon1) 'lonDiff Dim U1 As Double = Math.Atan((1 - f) * Math.Tan(rlat1)) Dim U2 As Double = Math.Atan((1 - f) * Math.Tan(rlat2)) Dim sinU1 As Double = Math.Sin(U1) Dim cosU1 As Double = Math.Cos(U1) Dim sinU2 As Double = Math.Sin(U2) Dim cosU2 As Double = Math.Cos(U2) Dim lambda As Double = L, lambdap As Double = 0.0 Dim iterLimit As Integer = 100 Dim sinLambda As Double Dim cosLambda As Double Dim sinSigma As Double Dim cosSigma As Double Dim Sigma As Double Dim sinAlpha As Double Dim cosSqAlpha As Double Dim cos2SigmaM As Double Dim C As Double Do sinLambda = Math.Sin(lambda) cosLambda = Math.Cos(lambda) sinSigma = Math.Sqrt((cosU2 * sinLambda) ^ 2 + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) ^ 2) If sinSigma = 0.0 Then Return 0.0 End If cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda Sigma = Math.Atan2(sinSigma, cosSigma) sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma cosSqAlpha = 1.0 - sinAlpha ^ 2 cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha If (Double.IsNaN(cos2SigmaM)) Then cos2SigmaM = 0 C = f / 16.0 * cosSqAlpha * (4.0 + f * (4.0 - 3.0 * cosSqAlpha)) lambdap = lambda lambda = L + (1 - C) * f * sinAlpha * (Sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2.0 * cos2SigmaM ^ 2))) iterLimit = iterLimit - 1 Loop While ((Math.Abs(lambda - lambdap) > 0.000000000001) And (iterLimit > 0)) If iterLimit = 0 Then Return Double.NaN Dim uSq As Double = cosSqAlpha * (a ^ 2 - b ^ 2) / (b ^ 2) Dim aA As Double = 1 + uSq / 16384.0 * (4096.0 + uSq * (-768.0 + uSq * (320.0 - 175.0 * uSq))) Dim bB As Double = uSq / 1024.0 * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq))) Dim deltaSigma As Double = bB * sinSigma * (cos2SigmaM + bB / 4.0 * (cosSigma * (-1 + 2 * cos2SigmaM ^ 2) - _ bB / 6.0 * cos2SigmaM * (-3 + 4 * sinSigma ^ 2) * (-3 + 4 * cos2SigmaM ^ 2))) Dim s As Double = b * aA * (Sigma - deltaSigma) Return s End Function ''' <summary> ''' 通过经纬度计算两点间距离 ''' </summary> ''' <param name="p0">起始点坐标</param> ''' <param name="p1">目的点坐标</param> ''' <returns>两点间距离,单位(米)</returns> ''' <remarks>采用WGS-84 椭球模型</remarks> Public Function GPSDistance(ByVal p0 As StationPoint, ByVal p1 As StationPoint) As Double Return GPSDistance(p0.X, p0.Y, p1.X, p1.Y) End Function ''' <summary> ''' 获取角度的弧度制表示 ''' </summary> ''' <param name="DEG">角度制表示(单位度)</param> ''' <returns></returns> ''' <remarks></remarks> Function GetRad(ByVal DEG As Double) As Double Return (Math.PI * DEG / 180) End Function