技术标签: 给定经纬度计算距离
通过经纬度坐标计算距离的方法(经纬度距离计算)
最近在网上搜索“通过经纬度坐标计算距离的方法”,发现网上大部分都是如下的代码:
#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是什么意思”。但从来没有见到有人回答。
提问的链接:
我刚开始接触这个问题,在中文世界中也只能搜到这些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 doubleBWD)
{
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 doublenewLat)
{
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 doubleminLongitude)
{
double temp = 0.0;
GetlatLon(centorlatitude, centorLogitude, distance, 0, out temp, outmaxLatitude);
GetlatLon(centorlatitude, centorLogitude, distance, 180, out temp, outminLatitude);
GetlatLon(centorlatitude, centorLogitude, distance, 90, out maxLongitude, outtemp);
GetlatLon(centorlatitude, centorLogitude, distance, 270, out minLongitude , outtemp);
}
这里得到的maxLatitude, minLatitude, maxLongitude, minLongitude就组成矩形的四个顶点。
网上另外的方法,
这里的“Destination point given distance and bearing from start point”一节介绍了。我直接把代码贴上来:
这里GetRectRange 这个函数 也是以正北方(due north)为起始角度,顺时针方向,求得maxLatitude, minLatitude, maxLongitude, minLongitude 这样一个MBR。2种方法的误差很小。我觉得都是可以用的公式。
///
///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
///
private static void GetNewLatLon(double lat, double lon, double d, double bearing, out double lat2, out doublelon2)
{
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
完整代码:
usingSystem;
usingSystem.Collections.Generic;
usingSystem.Linq;
usingSystem.Text;
namespaceGetJWD
{
public classGetMBR
{
public doubleMaxLatitude;
public doubleMinLatitude;
public doubleMaxLongitude;
public doubleMinLongitude;
public doubleMaxLatitude2;
public doubleMinLatitude2;
public doubleMaxLongitude2;
public doubleMinLongitude2;
public GetMBR(double centorlatitude, double centorLogitude, doubledistance)
{
GetRectRange(centorlatitude, centorLogitude, distance, out MaxLatitude, out MinLatitude, out MaxLongitude, outMinLongitude);
GetRectRange2(centorlatitude, centorLogitude, distance, out MaxLatitude2, out MinLatitude2, out MaxLongitude2, outMinLongitude2);
}
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 doublenewLat)
{
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, doubledistance,
out double maxLatitude, out double minLatitude, out double maxLongitude, out doubleminLongitude)
{
double temp = 0.0;
GetlatLon(centorlatitude, centorLogitude, distance, 0, out temp, outmaxLatitude);
GetlatLon(centorlatitude, centorLogitude, distance, 180, out temp, outminLatitude);
GetlatLon(centorlatitude, centorLogitude, distance, 90, out maxLongitude, outtemp);
GetlatLon(centorlatitude, centorLogitude, distance, 270, out minLongitude, outtemp);
}
public static void GetRectRange2(double centorlatitude, double centorLogitude, doubledistance,
out double maxLatitude, out double minLatitude, out double maxLongitude, out doubleminLongitude)
{
double temp = 0.0;
GetNewLatLon(centorlatitude, centorLogitude, distance, 0, out maxLatitude, outtemp);
GetNewLatLon(centorlatitude, centorLogitude, distance, 180, out minLatitude, outtemp);
GetNewLatLon(centorlatitude, centorLogitude, distance, 90, out temp, outmaxLongitude);
GetNewLatLon(centorlatitude, centorLogitude, distance, 270, out temp, outminLongitude);
}
///
///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
///
private static void GetNewLatLon(double lat, double lon, double d, double bearing, out double lat2, out doublelon2)
{
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随处可见的时代,希望此文能引导后人,少走我的弯路。如果有更好的方案,也欢迎留言。值此庆祝五一佳节!
文章浏览阅读1.6k次。安装配置gi、安装数据库软件、dbca建库见下:http://blog.csdn.net/kadwf123/article/details/784299611、检查集群节点及状态:[root@rac2 ~]# olsnodes -srac1 Activerac2 Activerac3 Activerac4 Active[root@rac2 ~]_12c查看crs状态
文章浏览阅读1.3w次,点赞45次,收藏99次。我个人用的是anaconda3的一个python集成环境,自带jupyter notebook,但在我打开jupyter notebook界面后,却找不到对应的虚拟环境,原来是jupyter notebook只是通用于下载anaconda时自带的环境,其他环境要想使用必须手动下载一些库:1.首先进入到自己创建的虚拟环境(pytorch是虚拟环境的名字)activate pytorch2.在该环境下下载这个库conda install ipykernelconda install nb__jupyter没有pytorch环境
文章浏览阅读5.2k次,点赞19次,收藏28次。选择scoop纯属意外,也是无奈,因为电脑用户被锁了管理员权限,所有exe安装程序都无法安装,只可以用绿色软件,最后被我发现scoop,省去了到处下载XXX绿色版的烦恼,当然scoop里需要管理员权限的软件也跟我无缘了(譬如everything)。推荐添加dorado这个bucket镜像,里面很多中文软件,但是部分国外的软件下载地址在github,可能无法下载。以上两个是官方bucket的国内镜像,所有软件建议优先从这里下载。上面可以看到很多bucket以及软件数。如果官网登陆不了可以试一下以下方式。_scoop-cn
文章浏览阅读4.5k次,点赞2次,收藏3次。首先要有一个color-picker组件 <el-color-picker v-model="headcolor"></el-color-picker>在data里面data() { return {headcolor: ’ #278add ’ //这里可以选择一个默认的颜色} }然后在你想要改变颜色的地方用v-bind绑定就好了,例如:这里的:sty..._vue el-color-picker
文章浏览阅读640次。基于芯片日益增长的问题,所以内核开发者们引入了新的方法,就是在内核中只保留函数,而数据则不包含,由用户(应用程序员)自己把数据按照规定的格式编写,并放在约定的地方,为了不占用过多的内存,还要求数据以根精简的方式编写。boot启动时,传参给内核,告诉内核设备树文件和kernel的位置,内核启动时根据地址去找到设备树文件,再利用专用的编译器去反编译dtb文件,将dtb还原成数据结构,以供驱动的函数去调用。firmware是三星的一个固件的设备信息,因为找不到固件,所以内核启动不成功。_exynos 4412 刷机
文章浏览阅读2w次,点赞24次,收藏42次。Linux系统配置jdkLinux学习教程,Linux入门教程(超详细)_linux配置jdk
文章浏览阅读3.3k次,点赞5次,收藏19次。xlabel('\delta');ylabel('AUC');具体符号的对照表参照下图:_matlab微米怎么输入
文章浏览阅读119次。顺序读写指的是按照文件中数据的顺序进行读取或写入。对于文本文件,可以使用fgets、fputs、fscanf、fprintf等函数进行顺序读写。在C语言中,对文件的操作通常涉及文件的打开、读写以及关闭。文件的打开使用fopen函数,而关闭则使用fclose函数。在C语言中,可以使用fread和fwrite函数进行二进制读写。 Biaoge 于2024-03-09 23:51发布 阅读量:7 ️文章类型:【 C语言程序设计 】在C语言中,用于打开文件的函数是____,用于关闭文件的函数是____。
文章浏览阅读3.4k次,点赞2次,收藏13次。跟随鼠标移动的粒子以grid(SOP)为partical(SOP)的资源模板,调整后连接【Geo组合+point spirit(MAT)】,在连接【feedback组合】适当调整。影响粒子动态的节点【metaball(SOP)+force(SOP)】添加mouse in(CHOP)鼠标位置到metaball的坐标,实现鼠标影响。..._touchdesigner怎么让一个模型跟着鼠标移动
文章浏览阅读178次。项目运行环境配置:Jdk1.8 + Tomcat7.0 + Mysql + HBuilderX(Webstorm也行)+ Eclispe(IntelliJ IDEA,Eclispe,MyEclispe,Sts都支持)。项目技术:Springboot + mybatis + Maven +mysql5.7或8.0+html+css+js等等组成,B/S模式 + Maven管理等等。环境需要1.运行环境:最好是java jdk 1.8,我们在这个平台上运行的。其他版本理论上也可以。_基于java技术的停车场管理系统实现与设计
文章浏览阅读3.5k次。前言对于MediaPlayer播放器的源码分析内容相对来说比较多,会从Java-&amp;gt;Jni-&amp;gt;C/C++慢慢分析,后面会慢慢更新。另外,博客只作为自己学习记录的一种方式,对于其他的不过多的评论。MediaPlayerDemopublic class MainActivity extends AppCompatActivity implements SurfaceHolder.Cal..._android多媒体播放源码分析 时序图
文章浏览阅读2.4k次,点赞41次,收藏13次。java 数据结构与算法 ——快速排序法_快速排序法