民航局MRV系统中的大圆距离算法 Vincenty算法 （java）

本文给出两种距离算法（Haversine formula；Vincenty formula）

Haversine formula

本计算式选取地球模型为球模型，以赤道半径为基准，故计算时纬度越高误差会越大，但胜在快速，具体推导公式可自行Google或参考美团技术文档（引：https://tech.meituan.com/lucene-distance.html （地址已失效））

package com.segment.position.calculate;


import com.segment.position.entity.Circle;
import com.segment.position.entity.Point;

import java.text.DecimalFormat;

/**
 * Haversine计算式，速度快误差较小，球模型
 */
public class DistanceCalculate
{
    private static final double EARTH_RADIUS = 6371e3;
    //private static final Logger logger = LoggerFactory.getLogger(DistanceCalculate.class);

    /**
     * 根据提供的两点经纬度获取两点间距离
     *
     * @param lat_one
     * @param lon_one
     * @param lat_two
     * @param lon_two
     * @return 两点距离
     */
    public static double getDistance(Double lat_one, Double lon_one, Double lat_two, Double lon_two)
    {
        double latOne = RadiansAngleCalculate.toRadians(lat_one);
        double latTwo = RadiansAngleCalculate.toRadians(lat_two);
        double latDiff = RadiansAngleCalculate.toRadians(lat_two - lat_one);
        double lonDiff = RadiansAngleCalculate.toRadians(lon_two - lon_one);

        double a = Math.sin(latDiff / 2) * Math.sin(latDiff / 2) + Math.cos(latOne) * Math.cos(latTwo) * Math.sin(lonDiff / 2) * Math.sin(lonDiff / 2);
        double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

        DecimalFormat format = new DecimalFormat("0.00");
        //logger.info("两点间距离：{" + format.format(EARTH_RADIUS * c) + "米}");
        return Double.parseDouble(format.format(EARTH_RADIUS * c));
    }

    /**
     * 计算两点间距离
     *
     * @param point_one
     * @param point_two
     * @return
     */
    public static double getDistance(Point point_one, Point point_two)
    {
        //获取点内信息，计算两点距离
        return getDistance(point_one.getLatitude(), point_one.getLongitude(), point_two.getLatitude(), point_two.getLongitude());
    }

    public static double getDistance(Circle circle, Point point)
    {
        //计算点到圆心的距离
        double point_to_center = getDistance(circle.getCenter(), point);
        //返回点距圆环的距离（圆内则为负，圆外则为正）
        return point_to_center - circle.getRadius();
    }

    /**
     * 根据提供的圆和点信息判断点是否在圆内
     *
     * @param circle
     * @param point
     * @return
     */
    public static boolean isPointInCircle(Circle circle, Point point)
    {
        return getDistance(circle, point) < 0;
    }

}


package com.segment.position.calculate;

public class RadiansAngleCalculate
{

    /**
     * 根据提供的角度值，将其转化为弧度
     *
     * @param angle 角度值
     * @return 结果
     */
    public static double toRadians(Double angle)
    {
        double result = 0L;
        if (angle != null) {
            result = angle * Math.PI / 180;
        }
        return result;
    }

    /**
     * 根据提供的经纬度，转化为弧度
     *
     * @param latitude  纬度
     * @param longitude 经度
     * @return 结果集
     */
    public static double[] latitudeLongitudeToRadians(Double latitude, Double longitude)
    {
        double[] result = new double[2];
        if (latitude != null && longitude != null) {
            result[0] = latitude * Math.PI / 180;
            result[1] = longitude * Math.PI / 180;
        }

        return result;
    }

    /**
     * 将给定的弧度转为角度
     *
     * @param radians 弧度值
     * @return 转换结果
     */
    public static double toAngle(Double radians)
    {
        double result = 0L;
        if (radians != null) {
            result = radians * Math.PI / 180;
        }
        return result;
    }

    /**
     * 将指定的经纬度弧度值转为角度
     *
     * @param latitude  纬度
     * @param longitude 经度
     * @return 角度
     */
    public static double[] latitudeLongitudeToAngle(Double latitude, Double longitude)
    {
        double[] result = new double[2];
        if (latitude != null && longitude != null) {
            result[0] = latitude * Math.PI / 180;
            result[1] = longitude * Math.PI / 180;
        }
        return result;
    }
}

Vincenty formula

本计算式理论模型为椭球模型，计算时会多次迭代，理论计算精度较高（为测地学里常用计算式推导结果）详细推导过程见wiki（引：https://en.wikipedia.org/wiki/Vincenty%27s_formulae）

package com.hna.utils;

import static java.lang.Float.NaN;
import java.text.DecimalFormat;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

/**
 * 
 * Vincenty公式计算椭球面两点间绝对距离<br/>
   * 理论计算结果精度可达到1米内<br/>
 * 
   * 参考文献：
 * https://www.movable-type.co.uk/scripts/latlong.html
 * https://www.movable-type.co.uk/scripts/latlong-vincenty.html
 * https://blog.csdn.net/qweasdzxc01233210/article/details/103405798
 * 
 * @author Dong
 *
 */
public class DistanceCalculateOfVincentyUtil {
	private static final Logger logger = LoggerFactory.getLogger(DistanceCalculateOfVincentyUtil.class);

	/** 长半径a=6378137 */
	private static final double a = 6378137;
	/** 短半径b=6356752.314245 */
	private static final double b = 6356752.314245;
	/** 扁率f=1/298.257223563 */
	private static final double f = 1 / 298.257223563;
	
	/**
     * 根据提供的经纬度计算两点间距离
     *
     * @param lat_one 坐标1纬度
     * @param lon_one 坐标1经度
     * @param lat_two 坐标2纬度
     * @param lon_two 坐标2经度
     * @return 两点间距离
     */
    public static double getDistance(Double lat_one, Double lon_one, Double lat_two, Double lon_two)
    {
    	logger.debug("坐标1经纬度("+lon_one+", "+lat_one+")");
    	logger.debug("坐标2经纬度("+lon_two+", "+lat_two+")");
        double L = DistanceCalculateOfVincentyUtil.toRadians(lon_one - lon_two);
        double U1 = Math.atan((1 - f) * Math.tan(DistanceCalculateOfVincentyUtil.toRadians(lat_one)));
        double U2 = Math.atan((1 - f) * Math.tan(DistanceCalculateOfVincentyUtil.toRadians(lat_two)));
        double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1),
                sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
        double lambda = L, lambdaP = Math.PI;
        double cosSqAlpha = 0L, sinSigma = 0L, cos2SigmaM = 0L, cosSigma = 0L, sigma = 0L;
        int circleCount = 40;
        //迭代循环
        while (Math.abs(lambda - lambdaP) > 1e-12 && --circleCount > 0) {
            double sinLambda = Math.sin(lambda), cosLambda = Math.cos(lambda);
            sinSigma = Math.sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda) +
                    (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda));
            if (sinSigma == 0) {
                return 0;  // co-incident points
            }
            cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
            sigma = Math.atan2(sinSigma, cosSigma);
            double alpha = Math.asin(cosU1 * cosU2 * sinLambda / sinSigma);
            cosSqAlpha = Math.cos(alpha) * Math.cos(alpha);
            cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
            double C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
            lambdaP = lambda;
            lambda = L + (1 - C) * f * Math.sin(alpha) *
                    (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
        }
        if (circleCount == 0) {
            return NaN;  // formula failed to converge
        }
        double uSq = cosSqAlpha * (a * a - b * b) / (b * b);
        double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
        double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
        double deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
                B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)));

        double result = b * A * (sigma - deltaSigma);
        DecimalFormat format = new DecimalFormat("0.000");
        double distance = Double.parseDouble(format.format(result));
        logger.debug("两点间距离：{" + distance + "米}");
        return distance;
    }

    /**
     * 根据提供的角度值，将其转化为弧度
     *
     * @param angle 角度值
     * @return 结果
     */
    public static double toRadians(Double angle)
    {
        double result = 0L;
        if (angle != null) {
            result = angle * Math.PI / 180;
        }
        return result;
    }
    
    
	public static void main(String[] args) {
		// TODO Auto-generated method stub
		//坐标1经度
		double lon_one = -5.71475;
		//坐标1纬度
		double lat_one = 50.06632;
		//坐标2经度
		double lon_two = -3.07009;
		//坐标2纬度
		double lat_two = 58.64402;
		DistanceCalculateOfVincentyUtil.getDistance(lat_one, lon_one, lat_two, lon_two);
	}

}

坐标转换（参考，原文地址没有找到）

本例提供火星坐标，百度坐标，WGS84坐标间互相转换方法

package com.segment.position.calculate;


import com.segment.position.entity.Point;

/**
 * 本工具类提供百度坐标（加密后的火星坐标BD09），高德坐标（火星坐标，GCJ02），WGS84坐标间互相转换的方法
 */
public class LocationConvert
{
    //private static final Logger logger = LoggerFactory.getLogger(LocationConvert.class);

    private static final double X_PI = 3.14159265358979324 * 3000.0 / 180.0;
    // π,Math库中的经度不够，转换误差较大
    private static final double PI = 3.1415926535897932384626;
    // 长半轴
    private static final double A = 6378245.0;
    // 扁率
    private static final double EE = 0.00669342162296594323;

    /**
     * 百度坐标系(BD-09)转WGS坐标
     *
     * @param pointBD 百度坐标点
     * @return WGS84坐标点
     */
    public static Point BDToWGS(Point pointBD)
    {
        Point point = BDToGCJ(pointBD);
        return GCJToWGS(point);
    }

    /**
     * WGS坐标转百度坐标系(BD-09)
     *
     * @param pointWGS WGS84坐标系点
     * @return 百度坐标数组
     */
    public static Point WGSToBD(Point pointWGS)
    {
        Point point = WGSToGCJ(pointWGS);
        return GCJToBD(point);
    }

    /**
     * 火星坐标系(GCJ-02)转百度坐标系(BD-09)
     * 谷歌、高德——>百度
     *
     * @param point 火星坐标点
     * @return 百度坐标数组
     */
    public static Point GCJToBD(Point point)
    {
        //logger.info("转换前坐标，GCJ-02：{" + point.getLongitude() + "," + point.getLatitude() + "}");
        double z = Math.sqrt(point.getLongitude() * point.getLongitude() +
                point.getLatitude() * point.getLatitude()) + 0.00002 * Math.sin(point.getLatitude() * X_PI);
        double theta = Math.atan2(point.getLatitude(), point.getLongitude()) + 0.000003 * Math.cos(point.getLongitude() * X_PI);
        double bd_lng = z * Math.cos(theta) + 0.0065;
        double bd_lat = z * Math.sin(theta) + 0.006;
        //logger.info("转换后坐标，BD-09：{" + bd_lng + "," + bd_lat + "}");
        Point result = new Point();
        result.setLatitude(bd_lat);
        result.setLongitude(bd_lng);
        return result;
    }

    /**
     * 百度坐标系(BD-09)转火星坐标系(GCJ-02)
     * 百度——>谷歌、高德
     *
     * @param bdPoint 百度坐标点
     * @return 火星坐标点
     */
    public static Point BDToGCJ(Point bdPoint)
    {
        //logger.info("转换前坐标，BD-09：{" + bdPoint.getLongitude() + "," + bdPoint.getLatitude() + "}");
        double x = bdPoint.getLongitude() - 0.0065;
        double y = bdPoint.getLatitude() - 0.006;
        double z = Math.sqrt(x * x + y * y) - 0.00002 * Math.sin(y * X_PI);
        double theta = Math.atan2(y, x) - 0.000003 * Math.cos(x * X_PI);
        double gg_lng = z * Math.cos(theta);
        double gg_lat = z * Math.sin(theta);
        //logger.info("转换后坐标，GCJ-02：{" + gg_lng + "," + gg_lat + "}");
        Point result = new Point();
        result.setLatitude(gg_lat);
        result.setLongitude(gg_lng);
        return result;
    }

    /**
     * WGS84转GCJ02(火星坐标系)
     *
     * @param pointWGS WGS84坐标系点
     * @return 火星坐标数组
     */
    public static Point WGSToGCJ(Point pointWGS)
    {
        //logger.info("转换前坐标，WGS-84：{" + pointWGS.getLongitude() + "," + pointWGS.getLatitude() + "}");
        if (outOfChina(pointWGS.getLongitude(), pointWGS.getLatitude())) {
            return new Point(pointWGS.getLongitude(), pointWGS.getLatitude());
        }
        double dlat = transformLat(pointWGS.getLongitude() - 105.0, pointWGS.getLatitude() - 35.0);
        double dlng = transformLng(pointWGS.getLongitude() - 105.0, pointWGS.getLatitude() - 35.0);
        double radlat = pointWGS.getLatitude() / 180.0 * PI;
        double magic = Math.sin(radlat);
        magic = 1 - EE * magic * magic;
        double sqrtmagic = Math.sqrt(magic);
        dlat = (dlat * 180.0) / ((A * (1 - EE)) / (magic * sqrtmagic) * PI);
        dlng = (dlng * 180.0) / (A / sqrtmagic * Math.cos(radlat) * PI);
        double mglat = pointWGS.getLatitude() + dlat;
        double mglng = pointWGS.getLongitude() + dlng;
        //logger.info("转换后坐标，GCJ-02：{" + mglng + "," + mglat + "}");
        Point result = new Point();
        result.setLongitude(mglng);
        result.setLatitude(mglat);
        return result;
    }

    /**
     * GCJ02(火星坐标系)转GPS84
     *
     * @param pointGCJ 火星坐标系点
     * @return WGS84坐标点
     */
    public static Point GCJToWGS(Point pointGCJ)
    {
        //logger.info("转换前坐标，GCJ-02：{" + pointGCJ.getLongitude() + "," + pointGCJ.getLatitude() + "}");
        if (outOfChina(pointGCJ.getLongitude(), pointGCJ.getLatitude())) {
            return new Point(pointGCJ.getLongitude(), pointGCJ.getLatitude());
        }
        double dlat = transformLat(pointGCJ.getLongitude() - 105.0, pointGCJ.getLatitude() - 35.0);
        double dlng = transformLng(pointGCJ.getLongitude() - 105.0, pointGCJ.getLatitude() - 35.0);
        double radlat = pointGCJ.getLatitude() / 180.0 * PI;
        double magic = Math.sin(radlat);
        magic = 1 - EE * magic * magic;
        double sqrtmagic = Math.sqrt(magic);
        dlat = (dlat * 180.0) / ((A * (1 - EE)) / (magic * sqrtmagic) * PI);
        dlng = (dlng * 180.0) / (A / sqrtmagic * Math.cos(radlat) * PI);
        double mglat = pointGCJ.getLatitude() + dlat;
        double mglng = pointGCJ.getLongitude() + dlng;
        //logger.info("转换后坐标，WGS-84：{" + (pointGCJ.getLongitude() * 2 - mglng) + "," + (pointGCJ.getLatitude() * 2 - mglat) + "}");
        Point result = new Point();
        result.setLatitude(mglat);
        result.setLongitude(mglng);
        return result;
    }

    /**
     * 纬度转换
     *
     * @param lng
     * @param lat
     * @return
     */
    private static double transformLat(double lng, double lat)
    {
        double ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + 0.1 * lng * lat + 0.2 * Math.sqrt(Math.abs(lng));
        ret += (20.0 * Math.sin(6.0 * lng * PI) + 20.0 * Math.sin(2.0 * lng * PI)) * 2.0 / 3.0;
        ret += (20.0 * Math.sin(lat * PI) + 40.0 * Math.sin(lat / 3.0 * PI)) * 2.0 / 3.0;
        ret += (160.0 * Math.sin(lat / 12.0 * PI) + 320 * Math.sin(lat * PI / 30.0)) * 2.0 / 3.0;
        return ret;
    }

    /**
     * 经度转换
     *
     * @param lng
     * @param lat
     * @return
     */
    private static double transformLng(double lng, double lat)
    {
        double ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + 0.1 * lng * lat + 0.1 * Math.sqrt(Math.abs(lng));
        ret += (20.0 * Math.sin(6.0 * lng * PI) + 20.0 * Math.sin(2.0 * lng * PI)) * 2.0 / 3.0;
        ret += (20.0 * Math.sin(lng * PI) + 40.0 * Math.sin(lng / 3.0 * PI)) * 2.0 / 3.0;
        ret += (150.0 * Math.sin(lng / 12.0 * PI) + 300.0 * Math.sin(lng / 30.0 * PI)) * 2.0 / 3.0;
        return ret;
    }

    /**
     * 判断是否在国内，不在国内不做偏移
     *
     * @param lng
     * @param lat
     * @return
     */
    private static boolean outOfChina(double lng, double lat)
    {
        if (lng < 72.004 || lng > 137.8347) {
            return true;
        } else if (lat < 0.8293 || lat > 55.8271) {
            return true;
        }
        return false;
    }
}

两种计算方法可酌情选取，不要盲目选择，具体数据分析可参考美团技术文档！

参考文档：
1.https://tech.meituan.com/lucene-distance.html
2.https://www.movable-type.co.uk/scripts/latlong.html
3.https://www.movable-type.co.uk/scripts/latlong-vincenty.html