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

2024-07-03 04:15| 来源: 网络整理| 查看: 265

已知两点经纬度，求两点之间绝对距离（Java）本文给出两种距离算法（Haversine formula；Vincenty formula）Haversine formulaVincenty formula坐标转换（参考，原文地址没有找到）

本文给出两种距离算法（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) 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公式计算椭球面两点间绝对距离 * 理论计算结果精度可达到1米内 * * 参考文献： * 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 137.8347) { return true; } else if (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

【本文地址】

公司简介

联系我们

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

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