经纬度计算

开发经常会遇到经纬度计算的相关场景。这次对相关知识做了下整理。

首先回顾一下科普知识:

1，经度： 英文 longitude 缩写 lng；纬度：英文 latitude  缩写 lat。

2，经度 是地球上一个地点离一根被称为本初子午线的南北方向走线以东或以西的度数。本初子午线的经度是0°。

    纬度 是指某点与地球球心的连线和地球赤道面所成的线面角。

3，角度的是60进制的，在程序中为了计算方便都存成了10进制的形式。

4，程序中的经纬度范围: 经度[-180,180] 纬度[-90,90]

5，在坐标系中怎么根据经纬度标柱一个点。

如图，建立一个空间直角坐标系，在赤道平面上，X轴指向本初子午线，Y轴在赤道平面垂直X轴，Z轴指向北极。

A[θ1,θ2]，经度θ2 ，纬度θ1。假设半径为单位1，A在坐标系的坐标（x,y,z）=(cosθ1*cosθ2, cosθ1*sinθ2, sinθ1)

6，正弦sin ，在坐标系中的曲线图

7，余弦cos,  在坐标系中的曲线图

9，空间两向量的夹角 

10，圆的周长：d=2∏r

    弧度: 弧长等于半径的弧，其所对的圆心角为1弧度。一个周长即2∏ 弧度(单位缩写是rad)。

    根据角度计算弧度: θ/360 *2∏= θ*∏/180 

    根据弧度计算角度: Ø*180 /∏ (Ø为弧度)

    根据弧度计算弧长: 弧度*半径

11，地球的赤道半径:  6377.830 = 6356.9088 + 20.9212  千米

      极地半径:6356.9088 千米

 

根据两点经纬度 A[lat1,lng1],B[lat2,lng2],计算两点间的球面距离。

方法：先计算出两点的弧度，然后再算球面距离。弧度可过通过两点对球心的夹角的余弦求出。

先假设地球半径为1

A点的坐标表示为(cos lat1*cos lng1, cos lat1*sin lng1, sin lat1)

B点的坐标表示为(cos lat2*cos lng2, cos lat2*sin lng2, sin lat2)

cos A^B= cos lat1*cos lng1 *cos lat2*cos lng2 +cos lat1*sin lng1*cos lat2*sin lng2+sin lat1*sin lat2

          = cos lat1*cos lat2(cos lng1*cos lng2+sin lng1*sin lng2) +sin lat1*sin lat2

          = cos lat1*cos lat2*cos(lng1-lng2) +sin lat1*sin lat2

球面距离= R*arccos(cos lat1*cos lat2*cos(lng1-lng2) +sin lat1*sin lat2)

注：lat1,lng1,lat2,lng2 为弧度，数学计算都是用弧度单位来计算的。

计算方法一：

代码如下：

/**
  *   方法一：
  */
private const double EARTH_RADIUS = 6377.830;//地球半径(千米)
        
        /// <summary>
        /// 角度换算成弧度
        /// </summary>
        /// <param name="d">角度</param>
        /// <returns>弧度</returns>
        private static double rad(double d)
        {
            return d * Math.PI / 180.0;
        }
         /// <summary>
        /// 计算距离
        /// </summary>
        /// <param name="lat1">纬度1</param>
        /// <param name="lng1">经度1</param>
        /// <param name="lat2">纬度2</param>
        /// <param name="lng2">经度2</param>
        /// <returns>距离（千米）</returns>
        public static double GetDistance(double lat1, double lng1, double lat2, double lng2)
        {
            double radLat1 = rad(lat1);
            double radLat2 = rad(lat2);
            double radLng1 = rad(lng1);
            double radLng2 = rad(lng2);
            double s = Math.Acos(Math.Cos(radLat1) * Math.Cos(radLat2) * Math.Cos(radLng1 - radLng2) + Math.Sin(radLat1) * Math.Sin(radLat2));
            s = s * EARTH_RADIUS;
            s = Math.Round(s * 10000) / 10000;
            return s;
        }


/**
  *  方法二：
  */
/// <summary>
        /// 计算距离
        /// </summary>
        /// <param name="lat1">纬度1</param>
        /// <param name="lng1">经度1</param>
        /// <param name="lat2">纬度2</param>
        /// <param name="lng2">经度2</param>
        /// <returns>距离（千米）</returns>
        public static double GetDistance2(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;
        }

根据某点A[lat,lng]，从地标库中查找周边N公里的地标。

方法：根据点A算出经线上N公里的角度偏差，和纬线上N公里的角度偏差。然后根据最大最小经纬度，画出矩形范围。再遍历矩形范围中的坐标计算出周边N公里的地标。

设地球半径为Re，则经线上的半径为R，纬线上的半径为 R*cos(lat) 。

经线上1公里对应角度=360/2∏R=180/∏R，

纬线上1公里对应角度=360/[2∏R*cos(lat)]=180/[∏R*cos(lat)]

代码如下：

/// <summary>
        /// 根据距离计算某点角度偏差
        /// </summary>
        /// <param name="lat">纬度</param>
        /// <param name="distance">距离（千米）</param>
        /// <param name="latDeviation">纬度偏差</param>
        /// <param name="lngDeviation">经度偏差</param>
        public static void GetMaxDeviation(double lat, double distance, out double latDeviation, out double lngDeviation)
        {
            double radLat1 = rad(lat);
            //另一种根据纬度计算地球半径的写法，因为极地半径和赤道半径有偏差。
            //EARTH_RADIUS = 6356.9088 + 20.9212 * (90.0 - lat) / 90.0; 
            double latRatio = 180 / (Math.PI * EARTH_RADIUS); //经线上1公里对应纬度偏差
            double lngRatio = latRatio / Math.Cos(radLat1);//纬线上1公里对应经度偏差
            latDeviation = distance * latRatio;
            lngDeviation = distance * lngRatio;
        }

 

另外博客上有人分享了3中计算经纬度距离的方法

计算方法二：

反余弦方法

// 分析 https://blog.csdn.net/jk940438163/article/details/83147557#commentsedit
package com.zhiliyuchi.web.rest.util;
 
/**
 * @创建人：Young
 * @时 间： 2019/3/13
 * @描 述: TODO
 */
public class Test {
    private static final double EARTH_RADIUS = 6371393; // 平均半径,单位：m；不是赤道半径。赤道为6378左右
 
    /**
     * @描述 反余弦进行计算
     * @参数  [lat1, lng1, lat2, lng2]
     * @返回值  double
     * @创建人  Young
     * @创建时间  2019/3/13 20:31
     **/
    public static double getDistance(Double lat1,Double lng1,Double lat2,Double lng2) {
        // 经纬度（角度）转弧度。弧度用作参数，以调用Math.cos和Math.sin
        double radiansAX = Math.toRadians(lng1); // A经弧度
        double radiansAY = Math.toRadians(lat1); // A纬弧度
        double radiansBX = Math.toRadians(lng2); // B经弧度
        double radiansBY = Math.toRadians(lat2); // B纬弧度
 
        // 公式中“cosβ1cosβ2cos（α1-α2）+sinβ1sinβ2”的部分，得到∠AOB的cos值
        double cos = Math.cos(radiansAY) * Math.cos(radiansBY) * Math.cos(radiansAX - radiansBX)
            + Math.sin(radiansAY) * Math.sin(radiansBY);
//        System.out.println("cos = " + cos); // 值域[-1,1]
        double acos = Math.acos(cos); // 反余弦值
//        System.out.println("acos = " + acos); // 值域[0,π]
//        System.out.println("∠AOB = " + Math.toDegrees(acos)); // 球心角 值域[0,180]
        return EARTH_RADIUS * acos; // 最终结果
    }
    public static void main(String[] args) {
        //121.717594,31.12055    121.817629,31.090867
        double distance = getDistance(31.12055, 121.717594,
            31.090867, 121.817629);
        System.out.println("距离" + distance  + "米");
    }
}

sin/cosin

package com.zhiliyuchi.web.rest.util;
 
/**
 * @创建人：Young
 * @时 间： 2019/3/13
 * @描 述: 高德地图对应经纬度计算距离
 */
public class LocationUtils {
    // 地球赤道半径
    private static double EARTH_RADIUS = 6378.137;
 
    private static double rad(double d) {
        return d * Math.PI / 180.0;
    }
 
    /**
     * @描述 经纬度获取距离，单位为米
     * @参数  [lat1, lng1, lat2, lng2]
     * @返回值  double
     * @创建人  Young
     * @创建时间  2019/3/13 20:33
     **/
    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 * 10000d) / 10000d;
        s = s * 1000;
        return s;
    }
 
    public static void main(String[] args) {
        double distance = getDistance(31.12055, 131.717594,
            21.090867, 111.817629);
        System.out.println("距离" + distance + "米");
    }
 
}

 

第三方jar包

// maven 依赖
        <dependency>
            　　<groupId>org.gavaghan</groupId>
            　　<artifactId>geodesy</artifactId>
            　　<version>1.1.3</version>
        </dependency>

// ==================================================================

package com.zhiliyuchi.web.rest.util;
 
import org.gavaghan.geodesy.Ellipsoid;
import org.gavaghan.geodesy.GeodeticCalculator;
import org.gavaghan.geodesy.GeodeticCurve;
import org.gavaghan.geodesy.GlobalCoordinates;
 
/**
 * @创建人：Young
 * @时 间： 2019/3/13
 * @描 述: TODO
 */
public class test0 {
    public static void main(String[] args)
    {
        // //121.717594,31.12055    121.817629,31.090867
        GlobalCoordinates source = new GlobalCoordinates(31.12055, 121.717594);
        GlobalCoordinates target = new GlobalCoordinates(31.090867, 121.817629);
 
        double meter1 = getDistanceMeter(source, target, Ellipsoid.Sphere);
        double meter2 = getDistanceMeter(source, target, Ellipsoid.WGS84);
 
        System.out.println("Sphere坐标系计算结果："+meter1 + "米");
        System.out.println("WGS84坐标系计算结果："+meter2 + "米");
    }
 
    public static double getDistanceMeter(GlobalCoordinates gpsFrom, GlobalCoordinates gpsTo, Ellipsoid ellipsoid)
    {
        //创建GeodeticCalculator，调用计算方法，传入坐标系、经纬度用于计算距离
        GeodeticCurve geoCurve = new GeodeticCalculator().calculateGeodeticCurve(ellipsoid, gpsFrom, gpsTo);
 
        return geoCurve.getEllipsoidalDistance();
    }
}

进入jar中的源码

{
        double a = ellipsoid.getSemiMajorAxis();
        double b = ellipsoid.getSemiMinorAxis();
        double f = ellipsoid.getFlattening();
        double phi1 = Angle.toRadians(start.getLatitude());
        double lambda1 = Angle.toRadians(start.getLongitude());
        double phi2 = Angle.toRadians(end.getLatitude());
        double lambda2 = Angle.toRadians(end.getLongitude());
        double a2 = a * a;
        double b2 = b * b;
        double a2b2b2 = (a2 - b2) / b2;
        double omega = lambda2 - lambda1;
        double tanphi1 = Math.tan(phi1);
        double tanU1 = (1.0D - f) * tanphi1;
        double U1 = Math.atan(tanU1);
        double sinU1 = Math.sin(U1);
        double cosU1 = Math.cos(U1);
        double tanphi2 = Math.tan(phi2);
        double tanU2 = (1.0D - f) * tanphi2;
        double U2 = Math.atan(tanU2);
        double sinU2 = Math.sin(U2);
        double cosU2 = Math.cos(U2);
        double sinU1sinU2 = sinU1 * sinU2;
        double cosU1sinU2 = cosU1 * sinU2;
        double sinU1cosU2 = sinU1 * cosU2;
        double cosU1cosU2 = cosU1 * cosU2;
        double lambda = omega;
        double A = 0.0D;
        double B = 0.0D;
        double sigma = 0.0D;
        double deltasigma = 0.0D;
        boolean converged = false;
 
        for(int s = 0; s < 20; ++s) {
            double lambda0 = lambda;
            double sinlambda = Math.sin(lambda);
            double coslambda = Math.cos(lambda);
            double sin2sigma = cosU2 * sinlambda * cosU2 * sinlambda + (cosU1sinU2 - sinU1cosU2 * coslambda) * (cosU1sinU2 - sinU1cosU2 * coslambda);
            double sinsigma = Math.sqrt(sin2sigma);
            double cossigma = sinU1sinU2 + cosU1cosU2 * coslambda;
            sigma = Math.atan2(sinsigma, cossigma);
            double sinalpha = sin2sigma == 0.0D?0.0D:cosU1cosU2 * sinlambda / sinsigma;
            double alpha = Math.asin(sinalpha);
            double cosalpha = Math.cos(alpha);
            double cos2alpha = cosalpha * cosalpha;
            double cos2sigmam = cos2alpha == 0.0D?0.0D:cossigma - 2.0D * sinU1sinU2 / cos2alpha;
            double u2 = cos2alpha * a2b2b2;
            double cos2sigmam2 = cos2sigmam * cos2sigmam;
            A = 1.0D + u2 / 16384.0D * (4096.0D + u2 * (-768.0D + u2 * (320.0D - 175.0D * u2)));
            B = u2 / 1024.0D * (256.0D + u2 * (-128.0D + u2 * (74.0D - 47.0D * u2)));
            deltasigma = B * sinsigma * (cos2sigmam + B / 4.0D * (cossigma * (-1.0D + 2.0D * cos2sigmam2) - B / 6.0D * cos2sigmam * (-3.0D + 4.0D * sin2sigma) * (-3.0D + 4.0D * cos2sigmam2)));
            double C = f / 16.0D * cos2alpha * (4.0D + f * (4.0D - 3.0D * cos2alpha));
            lambda = omega + (1.0D - C) * f * sinalpha * (sigma + C * sinsigma * (cos2sigmam + C * cossigma * (-1.0D + 2.0D * cos2sigmam2)));
            double change = Math.abs((lambda - lambda0) / lambda);
            if(s > 1 && change < 1.0E-13D) {
                converged = true;
                break;
            }
        }
 
        double var96 = b * A * (sigma - deltasigma);
        double alpha1;
        double alpha2;
        if(!converged) {
            if(phi1 > phi2) {
                alpha1 = 180.0D;
                alpha2 = 0.0D;
            } else if(phi1 < phi2) {
                alpha1 = 0.0D;
                alpha2 = 180.0D;
            } else {
                alpha1 = 0.0D / 0.0;
                alpha2 = 0.0D / 0.0;
            }
        } else {
            double radians = Math.atan2(cosU2 * Math.sin(lambda), cosU1sinU2 - sinU1cosU2 * Math.cos(lambda));
            if(radians < 0.0D) {
                radians += 6.283185307179586D;
            }
 
            alpha1 = Angle.toDegrees(radians);
            radians = Math.atan2(cosU1 * Math.sin(lambda), -sinU1cosU2 + cosU1sinU2 * Math.cos(lambda)) + 3.141592653589793D;
            if(radians < 0.0D) {
                radians += 6.283185307179586D;
            }
 
            alpha2 = Angle.toDegrees(radians);
        }
 
        if(alpha1 >= 360.0D) {
            alpha1 -= 360.0D;
        }
 
        if(alpha2 >= 360.0D) {
            alpha2 -= 360.0D;
        }
 
        return new GeodeticCurve(var96, alpha1, alpha2);
    }