通过经纬度坐标计算距离的方法(实为通过一个经纬度和距离角度求另一个经纬度)

article/2025/8/27 16:45:19

转自:https://www.cnblogs.com/softfair/p/lat_lon_distance_bearing_new_lat_lon.html

通过经纬度坐标计算距离的方法(经纬度距离计算)

通过经纬度坐标计算距离的方法(经纬度距离计算)

最近在网上搜索“通过经纬度坐标计算距离的方法”,发现网上大部分都是如下的代码:

#define PI 3.14159265

static double Rc = 6378137;  // 赤道半径

static double Rj = 6356725;  // 极半径

class JWD

{

public:

double m_Longitude, m_Latitude;

double m_RadLo, m_RadLa;

double Ec;

double Ed;

 

public:

JWD(double longitude, double latitude)

{

  m_Longitude = longitude;

  m_Latitude = latitude;

      m_RadLo = longitude * PI/180.;

  m_RadLa = latitude * PI/180.;

  Ec = Rj + (Rc - Rj) * (90.-m_Latitude) / 90.;

  Ed = Ec * cos(m_RadLa);

}

~JWD() {};

};

 

static JWD GetJWDB(JWD A, double x,double y)

{

double dx=x;

double dy=y;

double BJD = (dx/A.Ed + A.m_RadLo) * 180./PI;

double BWD = (dy/A.Ec + A.m_RadLa) * 180./PI;

JWD B(BJD, BWD);

return B;

}

void main()

{

double referla=30.0;

double referlo=60.0;

double dx=500.0;

double dy=60.0;

    JWD A(referla,referlo),B(0.0,0.0);

    B=GetJWDB(A,dx,dy);

cout < < "  LA = " < < B.m_Latitude < < "  LO= " < < B.m_Longitude < < endl;

}

 

上面这段与之类似的代码是最容易通过搜索引擎找到的。大部分都是相互的转载,从来没有说明和注释。给初学者造成了很大的困扰。特别是:

Ec = Rj + (Rc - Rj) * (90.-m_Latitude) / 90.;

Ed = Ec * cos(m_RadLa);

Ec、Ed这2个参数,有人还在CSDN上发帖询问“函数中Ec,Ed是什么意思”。但从来没有见到有人回答。

提问的链接:

http://www.gisforum.net/bbs/TopicOther.asp?t=5&BoardID=33&id=155609

http://bbs.csdn.net/topics/320024634

我刚开始接触这个问题,在中文世界中也只能搜到这些Ctrl+C 到Ctrl+V的东西。全球最大的中文互联网上也没有任何解答。已经明白的人知道很简单,但初学者肯定很疑惑。更何况现在LBS应用已经非常多了。肯定有非常多的人已经找到了更好的解答方式。 但对于我来说,用最常用的关键词,最容易的找到了这些复制的答案,但疑点重重。本着好奇的心,经过一番了解,我把结果给大家分析一下。下面是C#的代码:

        

复制代码

        public const double Ea = 6378137;     //   赤道半径  public const double Eb = 6356725;     //   极半径 private static void GetJWDB(double GLAT, double GLON, double distance, double angle, out double BJD, out double BWD){double dx = distance * 1000 * Math.Sin(angle * Math.PI / 180.0);double dy = distance * 1000 * Math.Cos(angle * Math.PI / 180.0);//double ec = 6356725 + 21412 * (90.0 - GLAT) / 90.0;// 21412 是赤道半径与极半径的差double ec = Eb + (Ea-Eb) * (90.0 - GLAT) / 90.0;double ed = ec * Math.Cos(GLAT * Math.PI / 180);BJD = (dx / ed + GLON * Math.PI / 180.0) * 180.0 / Math.PI;BWD = (dy / ec + GLAT * Math.PI / 180.0) * 180.0 / Math.PI;}

复制代码

 

上面这个函数一看就是懂中文的人搞的,名字叫GetJWDB,取得经纬度。他根据输入的经度、纬度、距离和一个角度,得到另外一个经纬度坐标,输出参数为BJD、BWD。

1)angle * Math.PI / 180.0 作用是将角度转换为弧度,经纬度坐标是角度值,计算时需要换为弧度。这里所有的计算都是用弧度。

2)函数以正北方(due north) 也就是指南针的方向为0度,顺时针方向增加。如下图,Distance距离如果是d的话,dx就是x轴方向的长度,即longitude经度方向的长度;dy就是y轴方向的长度,即latitude纬度方向的长度。

dx、dy的计算方式也可以是以正东(due east)方向为0度。那:dx=distance* Cos(θ),dy=distance*Sin(θ)。区别是Cos 与Sin互换。(图上未标)

图1

 

3)Ea 表示赤道半径,Eb表示极半径,地球是一个近似球体,Ea与Eb稍微有点差距。ec的作用就是修正因为纬度不断变化的球半径长度。如果在GLAT=0,即在赤道上的时候,ec=Eb+(Ea-Eb)*(90-0)/90=Ea,那ec就刚好是赤道半径Ea;如果在极点GLAT=90,ec=Eb+(Ea-Eb)*(90-90)/90=Eb,那ec 就刚好是极半径Eb。

4)Ed是GLAT所在纬度的纬度圈的半径,如下图:

 

图2

截面过球心,此时截面的面积最大,此圆叫球的大圆(Great Cycle),沿着经线进行截面,得到的都是大圆(Great Cycle)。球面被不经过球心. 的截面所截得的圆. 叫做小圆。纬度圈所在的圆是一个小圆。地球半径R,平均值R=6371.0Km, Ed如果用地球半径R表示,那就是Ed=R*Cos(θ),可以参看

《根据2个经纬度点,计算这2个经纬度点之间的距离(通过经度纬度得到距离)》这里提到的“A、D所在纬度圆圈的半径(AO`)”。放到上面函数里,因为它不断修正地球半径ec,那就是ed = ec * Math.Cos(GLAT * Math.PI / 180);

 

5)按照地球经纬度的划分,如下图:

每等分的纬度之间,经度线段的长度是一定的。 A段,B段长度是一样的。

每一等分的经度之间,纬度线段的长度是从赤道向2极点减小的。C端,D段的长度是不一样的。

图3

 

参看上图,那dx / ed 就相当于是在GLAT这个纬度上dx长度与总长度的占比,算出来应该是个经度跨度。如果这个经度跨度加上起始给定的经度就是最终的经度。

同理 dy/R就是在GLON这个经度上的dy长度与地球平均半径R的占比,算出来应该是一个纬度跨度。如果这个纬度跨度加上起始给定的纬度就是最终的纬度。这里使用了R,取地球平均半径。

dy/ec 就是用不断修正的ec半径替代了平均半径R。

(dy / ec + GLAT * Math.PI / 180.0) 就是起始纬度加上distance距离的最终纬度,同时需要将该结果转换为角度。 转换角度方法是:弧度* 180.0 / Math.PI。

BWD = (dy / ec + GLAT * Math.PI / 180.0) * 180.0 / Math.PI;

(dx / ed + GLON * Math.PI / 180.0)就是起始经度加上distance距离的最终经度,同时需要将该结果转换为角度。

BJD = (dx / ed + GLON * Math.PI / 180.0) * 180.0 / Math.PI;

 

这个根据一个经纬度坐标、距离然后求另外一个经纬度坐标的作用,主要就是确定一个最小外包矩形(Minimum bounding rectangle,简称MBR)。例如,我要找一个坐标点(lat,lon)的5公里范围内的所有商户信息、景点信息等。这个MBR就是一个最大的范围,这个矩形是包含5公里范围内所有这些有效信息的一个最小矩形。利用公式,求出四个方向0度、90度、180度、270度方向上的四个坐标点就可以得到这个MBR。

 

       

复制代码

        public const double Ea = 6378137;     //   赤道半径  public const double Eb = 6356725;     //   极半径 private static void GetlatLon(double LAT, double LON, double distance, double angle, out double newLon, out double newLat){double dx = distance * 1000 * Math.Sin(angle * Math.PI / 180.0);double dy = distance * 1000 * Math.Cos(angle * Math.PI / 180.0);double ec = Eb + (Ea-Eb) * (90.0 - LAT) / 90.0;double ed = ec * Math.Cos(LAT * Math.PI / 180);newLon = (dx / ed + LON * Math.PI / 180.0) * 180.0 / Math.PI;newLat = (dy / ec + LAT * Math.PI / 180.0) * 180.0 / Math.PI;}public static void GetRectRange(double centorlatitude, double centorLogitude, double distance, out double maxLatitude, out double minLatitude, out double maxLongitude, out double minLongitude){double temp = 0.0;GetlatLon(centorlatitude, centorLogitude, distance, 0, out temp, out maxLatitude);GetlatLon(centorlatitude, centorLogitude, distance, 180, out temp, out minLatitude);GetlatLon(centorlatitude, centorLogitude, distance, 90, out maxLongitude, out temp);GetlatLon(centorlatitude, centorLogitude, distance, 270, out minLongitude , out temp);}

复制代码

 

这里得到的maxLatitude, minLatitude, maxLongitude, minLongitude就组成矩形的四个顶点。

 

网上另外的方法,

http://www.movable-type.co.uk/scripts/latlong.html

这里的“Destination point given distance and bearing from start point”一节介绍了。我直接把代码贴上来:

这里GetRectRange 这个函数 也是以正北方(due north)为起始角度,顺时针方向,求得maxLatitude, minLatitude, maxLongitude, minLongitude 这样一个MBR。2种方法的误差很小。我觉得都是可以用的公式。

复制代码

         /// <summary>/// where    φ is latitude, λ is longitude, θ is the bearing (clockwise from north),/// δ is the angular distance d/R; d being the distance travelled, R the earth’s radius/// bearing 方位 0,90,180,270/// </summary>private static void GetNewLatLon(double lat, double lon, double d, double bearing, out double lat2, out double lon2){lat2 = 0.0;lon2 = 0.0;double R = 6378.137;var φ1 = ConvertDegreesToRadians(lat);var λ1 = ConvertDegreesToRadians(lon);var θ = ConvertDegreesToRadians(bearing);var φ2 = Math.Asin(Math.Sin(φ1) * Math.Cos(d / R) +Math.Cos(φ1) * Math.Sin(d / R) * Math.Cos(θ));var λ2 = λ1 + Math.Atan2(Math.Sin(θ) * Math.Sin(d / R) * Math.Cos(φ1),Math.Cos(d / R) - Math.Sin(φ1) * Math.Sin(φ2));λ2 = (λ2 + 3 * Math.PI) % (2 * Math.PI) - Math.PI; // normalise to -180..+180°lat2 = ConvertRadiansToDegrees(φ2);lon2 = ConvertRadiansToDegrees(λ2);}

复制代码

 

如果有一个应用,表里存有100万的数据,数据包含一个lat、lon的经纬度信息。就可以先根据输入的经纬度和距离得到一个MBR,然后通过类似

SELECT Id

FROM IdInfoTable

WHERE latitude >= minLat AND latitude < maxLat

AND longitude >= minLon AND longitude < maxLon

的方式过滤掉大部分的数据,然后再通过《根据2个经纬度点,计算这2个经纬度点之间的距离(通过经度纬度得到距离)》这里提到的方法进行精细过滤。

 

完整代码:

复制代码

using System;using System.Collections.Generic;using System.Linq;using System.Text;namespace GetJWD{public class GetMBR{public double MaxLatitude;public double MinLatitude;public double MaxLongitude;public double MinLongitude;public double MaxLatitude2;public double MinLatitude2;public double MaxLongitude2;public double MinLongitude2;public GetMBR(double centorlatitude, double centorLogitude, double distance){GetRectRange(centorlatitude, centorLogitude, distance, out MaxLatitude, out MinLatitude, out MaxLongitude, out MinLongitude);GetRectRange2(centorlatitude, centorLogitude, distance, out MaxLatitude2, out MinLatitude2, out MaxLongitude2, out MinLongitude2);}public const double Ea = 6378137;     //   赤道半径  public const double Eb = 6356725;     //   极半径 private static void GetlatLon(double LAT, double LON, double distance, double angle, out double newLon, out double newLat){double dx = distance * 1000 * Math.Sin(angle * Math.PI / 180.0);double dy = distance * 1000 * Math.Cos(angle * Math.PI / 180.0);double ec = Eb + (Ea - Eb) * (90.0 - LAT) / 90.0;double ed = ec * Math.Cos(LAT * Math.PI / 180);newLon = (dx / ed + LON * Math.PI / 180.0) * 180.0 / Math.PI;newLat = (dy / ec + LAT * Math.PI / 180.0) * 180.0 / Math.PI;}public static void GetRectRange(double centorlatitude, double centorLogitude, double distance,out double maxLatitude, out double minLatitude, out double maxLongitude, out double minLongitude){double temp = 0.0;GetlatLon(centorlatitude, centorLogitude, distance, 0, out temp, out maxLatitude);GetlatLon(centorlatitude, centorLogitude, distance, 180, out temp, out minLatitude);GetlatLon(centorlatitude, centorLogitude, distance, 90, out maxLongitude, out temp);GetlatLon(centorlatitude, centorLogitude, distance, 270, out minLongitude, out temp);}public static void GetRectRange2(double centorlatitude, double centorLogitude, double distance,out double maxLatitude, out double minLatitude, out double maxLongitude, out double minLongitude){double temp = 0.0;GetNewLatLon(centorlatitude, centorLogitude, distance, 0, out maxLatitude, out temp);GetNewLatLon(centorlatitude, centorLogitude, distance, 180, out minLatitude, out temp);GetNewLatLon(centorlatitude, centorLogitude, distance, 90, out temp, out maxLongitude);GetNewLatLon(centorlatitude, centorLogitude, distance, 270, out temp, out minLongitude);}/// <summary>/// where    φ is latitude, λ is longitude, θ is the bearing (clockwise from north),/// δ is the angular distance d/R; d being the distance travelled, R the earth’s radius/// bearing 方位 0,90,180,270/// </summary>private static void GetNewLatLon(double lat, double lon, double d, double bearing, out double lat2, out double lon2){lat2 = 0.0;lon2 = 0.0;double R = 6378.137;var φ1 = ConvertDegreesToRadians(lat);var λ1 = ConvertDegreesToRadians(lon);var θ = ConvertDegreesToRadians(bearing);var φ2 = Math.Asin(Math.Sin(φ1) * Math.Cos(d / R) +Math.Cos(φ1) * Math.Sin(d / R) * Math.Cos(θ));var λ2 = λ1 + Math.Atan2(Math.Sin(θ) * Math.Sin(d / R) * Math.Cos(φ1),Math.Cos(d / R) - Math.Sin(φ1) * Math.Sin(φ2));λ2 = (λ2 + 3 * Math.PI) % (2 * Math.PI) - Math.PI; // normalise to -180..+180°lat2 = ConvertRadiansToDegrees(φ2);lon2 = ConvertRadiansToDegrees(λ2);}public static double ConvertDegreesToRadians(double degrees){return degrees * Math.PI / 180;}public static double ConvertRadiansToDegrees(double radian){return radian * 180.0 / Math.PI;}}class Test{static void Main(string[] args){double latorg = 22.54587746, lonorg = 114.12873077;var gpsdis = new GetMBR(latorg, lonorg, 5);Console.WriteLine("maxlat:{0}, minlat:{1}, maxlon:{2}, minlon:{3}",gpsdis.MaxLatitude, gpsdis.MinLatitude, gpsdis.MaxLongitude, gpsdis.MinLongitude);Console.WriteLine("maxlat:{0}, minlat:{1}, maxlon:{2}, minlon:{3}",gpsdis.MaxLatitude2, gpsdis.MinLatitude2, gpsdis.MaxLongitude2, gpsdis.MinLongitude2);}}}

复制代码

 

 

谨以此文纪念那篇CSDN上因为 “本帖子已过去太久远了,不再提供回复功能。”而永远至今晚为止都还没有答案的帖子!

如今LBS应用泛滥,JavaScript到处都能看到源码,gitHub上sourceCode随处可见的时代,希望此文能引导后人,少走我的弯路。如果有更好的方案,也欢迎留言。值此庆祝五一佳节!


http://chatgpt.dhexx.cn/article/EkYXMKwU.shtml

相关文章

经纬度计算距离

通过经纬度坐标计算距离的方法&#xff08;经纬度距离计算&#xff09; 最近在网上搜索“通过经纬度坐标计算距离的方法”&#xff0c;发现网上大部分都是如下的代码&#xff1a; #define PI 3.14159265 static double Rc 6378137; // 赤道半径 static double Rj 6356725; …

Mysql 根据经纬度计算距离

方式1&#xff1a;st_distance_sphere st_distance_sphere函数返回以米为单位 SELECT*,st_distance_sphere (point ( longitudes, latitudes ),point ( 113.264435, 23.129163 )) AS juli FROMzxh_distance ORDER BYjuli ASC;方式2&#xff1a;st_distance st_distance返回的…

计算经纬度点之间的距离

球面上任意两点之间的距离计算公式可以参考维基百科上的下述文章。 Great-circle distance Haversine formula 值得一提的是&#xff0c;维基百科推荐使用Haversine公式&#xff0c;理由是Great-circle distance公式用到了大量余弦函数&#xff0c; 而两点间距离很短时&#x…

计算两个经纬度的距离

一 点睛 计算两个经纬度的距离的方法这里提供两种。 1 基于 googleMap 中的算法得到两经纬度之间的距离&#xff0c;计算精度与谷歌地图的距离精度差不多。 2 计算中心经纬度与目标经纬度的距离(米) 二 代码 public class LongitudeLatitude {/*** 默认地球半径*/private …

根据两点经纬度计算距离

根据两点经纬度计算距离 这些经纬线是怎样定出来的呢&#xff1f;地球是在不停地绕地轴旋转&#xff08;地轴是一根通过地球南北两极和地球中心的 假想线&#xff09;&#xff0c;在地球中腰画一个与地轴垂直的大圆圈&#xff0c;使圈上的每一点都和南北两极的距离相等&#xf…

通过经纬度坐标计算距离的方法(经纬度距离计算)

通过经纬度坐标计算距离的方法&#xff08;经纬度距离计算&#xff09; 最近在网上搜索“通过经纬度坐标计算距离的方法”&#xff0c;发现网上大部分都是如下的代码&#xff1a; #define PI 3.14159265 static double Rc 6378137; // 赤道半径 static double Rj 635672…

根据经纬度计算两点之间的距离

前言 在我们平时使用美团&#xff0c;饿了么等app进行订餐&#xff0c;或者使用猫眼进行订电影票的时候&#xff0c;都有一个距离的排序&#xff0c;表明该家店距离我们当前的位置&#xff0c;这种基于地理位置的服务&#xff0c;统一被称为LBS&#xff08;Location Based Ser…

logistic模型原理与推导过程分析(1)

从线性分类器谈起 给定一些数据集合&#xff0c;他们分别属于两个不同的类别。例如对于广告数据来说&#xff0c;是典型的二分类问题&#xff0c;一般将被点击的数据称为正样本&#xff0c;没被点击的数据称为负样本。现在我们要找到一个线性分类器&#xff0c;将这些数据分为两…

阻滞增长模型--Logistic模型

文章目录 写在前面一、思路二、原始数据三、源代码四、结果 写在前面 做了一个作业&#xff0c;觉得挺有意思&#xff0c;有空就发上来了。 阻滞增长模型–Logistic模型的部分推导、思路&#xff0c;并给出参考代码。 参考书籍&#xff1a;《数学模型》第四版&#xff0c;姜启…

logistic模型原理与推导过程分析(2)

二项逻辑回归模型 既然logistic回归把结果压缩到连续的区间(0,1)&#xff0c;而不是离散的0或者1&#xff0c;然后我们可以取定一个阈值&#xff0c;通常以0.5为阈值&#xff0c;如果计算出来的概率大于0.5&#xff0c;则将结果归为一类&#xff08;1&#xff09;&#xff0c;…

从logit变换到logistic模型

从logit变换到logistic模型 logit变换几率logistic模型 前面我们知道对数函数和对数函数的一些基本性质&#xff0c;也许你会问&#xff0c;为什么要引入对数函数&#xff1f;而且还是一个基本初等函数&#xff1f;这就要从logit变换说起。 logit变换 我们在研究某一结果&…

Gompertz模型绘图 matlab,Logistic模型matlab求解

Logistic模型求解怎么用matlab求解啊&#xff1f; 悬赏分&#xff1a;100 - 解决时间&#xff1a;2008-11-17 23:09 已知 x0:1:12 y[43.65 109.86 187.21 312.67 496.58 707.65 960.25 1238.75 1560.00 1824.29 2199.00 2438.89 2737.71] yL/(1a*exp(-k*x)) 利用线性回归模型所…

数学建模-Logistic模型

文章目录 Malthus模型模型假设建模与求解模型评价 Logistic模型模型假设建模与求解模型检验 为了更好地理解Logistic模型&#xff0c;我们先看看Malthus模型 Malthus模型 这是英国神父Malthus通过对一百多年人口统计资料的分析之后提出的人口模型假设 模型假设 设x(t)表示t时…

数学建模-Logistic模型附Matlab代码

目录 一、Logistic模型介绍 二、Logistic模型实例 三、Logistic模型原理 3.1 Logistic 方程定义 3.2 Yule算法 3.2 Rhodes算法 3.3 Nair算法 4、Logistic模型Matlab部分代码 4.1 Yule算法 4.2 Rhodes算法 4.3 Nair算法 一、Logistic模型介绍 logistic回归又称logi…

logistic回归模型—基于R

logistic回归模型—基于R 数据理解和准备一. 对缺失值的处理二.虚拟变量的赋值三.箱线图四.相关性分析 训练集与测试集的划分模型构建与评价一.logistic回归模型二.检查模型在训练数据集和测试数据集上的表现使用交叉验证的logistic回归 logistic回归又称logistic回归分析&…

数学模型——Logistic回归模型(含Matlab代码)

写在前面 Logistic回归模型是一种非常常见的统计回归模型&#xff0c;在处理大量数据&#xff0c;揭示各自变量如何作用于因变量&#xff08;描述X与Y之间的关系&#xff09;时有着十分重要的作用。笔者在写Logit回归模型前参加了一次市场调研比赛&#xff0c;在这次比赛中学到…

菜鸟的数学建模之路(五):Logistic模型

matlab实现 Logistic回归跟多元线性回归差不多&#xff0c;但是有区别&#xff1a; &#xff08;1&#xff09; 线性回归&#xff1a;y是一个定量的变量&#xff0c;这时y对于不同的自变量来说有相应的值。 &#xff08;2&#xff09; Logistic回归&#xff1a;y是一个定性的变…

自学鸿蒙应用开发(25)- 基本的CommonDialog

动作演示 对话框是应用程序的主要输入手段之一&#xff0c;但是遗憾的是目前鸿蒙的开发网站上只有ToastDialog用法&#xff0c;其他类型的对话框则只能参考为数不多的英文文档。 以下是作者经过不断尝试&#xff0c;终于鼓捣出来的CommonDialog对话框。 CommonDialog1 代码实…

Builder设计模式构建通用型Dialog

目录 写在前面 一、什么是Builder模式 二、AlertDialog源码分析 2.1、源码阅读 2.2、Builder模式工作流程 三、代码实战——Builder模式构建通用型Dialog 3.1、基本框架搭建 3.2、完善Builder 3.3、完善真正的构建器 3.4、自定义参数配置 四、使用Dialog 写在前面 …

自学鸿蒙应用开发(26)- 自定义CommonDialog

执行效果 上一篇文章中说过&#xff0c;直接使用鸿蒙系统中的CommonDialog大致是下面的效果&#xff1a; 这个效果实在是无法用于实际的应用开发。本文介绍如何定制自己的CommonDialog。还是先看演示视频&#xff1a; CustomizeCommonDialog 准备布局 定制CommonDialog的第一…