GPS坐标转笛卡尔坐标

提要

我们知道GPS坐标是由经度，纬度，海拔组成，精度和纬度都是角度，海报是高度。

在进行基于地理的搜索的时候，常用到KDTree，在构建在KDTree的时候，不能直接用GPS的坐标，要将GPS坐标转换成笛卡尔坐标才能用于构建KDTree。下面就是相关的转换算法。

注：GPS信息由几种标准，这里的采用的是google map的经纬度信息。

js实现

geodecy.js

/*             geodesy routines in JavaScriptJames R. Clynch NPS / 2003Done for support of web education pages== must convert inputs to numbers for safety ==== if string comes in - sometimes works, sometimes not !==
*/<!--// =======================================================================function geodGBL()//             test and ensure geodesy globals loaded{var  tstglobaltstglobal = typeof EARTH_A;if ( tstglobal == "undefined" )  wgs84() }// =======================================================================function earthcon(ai,bi)/*        Sets Earth Constants as globals--  input a,b--  Leaves Globals EARTH_A      EARTH_B   EARTH_F  EARTH_Ecc
*/{var  f,ecc, eccsq, a,ba        =  Number(ai);b        =  Number(bi);f        =  1-b/a;eccsq    =  1 - b*b/(a*a);ecc      =  Math.sqrt(eccsq);EARTH_A  =  a;EARTH_B  =  b;EARTH_F  =  f;EARTH_Ecc=  ecc;EARTH_Esq=  eccsq;}// =======================================================================function wgs84()/*        WGS84 Earth Constants--  returns a,b,f,e  ----  Leaves Globals EARTH_A      EARTH_B   EARTH_F  EARTH_Ecc*/{var  wgs84a, wgs84b, wgs84fwgs84a         =  6378.137;wgs84f         =  1.0/298.257223563;wgs84b         =  wgs84a * ( 1.0 - wgs84f );earthcon (wgs84a, wgs84b );}          // =======================================================================function  radcur(lati)/*compute the radii at the geodetic latitude lat (in degrees)input:lat       geodetic latitude in degreesoutput:   rrnrm     an array 3 longr,  rn,  rm   in km*/{var rrnrm = new Array(3)var dtr   = Math.PI/180.0var  a,b,latvar  asq,bsq,eccsq,ecc,clat,slatvar  dsq,d,rn,rm,rho,rsq,r,z//        -------------------------------------geodGBL();a     = EARTH_A;b     = EARTH_B;asq   = a*a;bsq   = b*b;eccsq  =  1 - bsq/asq;ecc = Math.sqrt(eccsq);lat   =  Number(lati);clat  =  Math.cos(dtr*lat);slat  =  Math.sin(dtr*lat);dsq   =  1.0 - eccsq * slat * slat;d     =  Math.sqrt(dsq);rn    =  a/d;rm    =  rn * (1.0 - eccsq ) / dsq;rho   =  rn * clat;z     =  (1.0 - eccsq ) * rn * slat;rsq   =  rho*rho + z*z;r     =  Math.sqrt( rsq );rrnrm[0]  =  r;rrnrm[1]  =  rn;rrnrm[2]  =  rm;return ( rrnrm );}// =======================================================================//        physical radius of earth from geodetic latitudefunction  rearth (lati){var    rrnrm, r,latlat   =  Number(lati);rrnrm =  radcur ( lat );r     =  rrnrm[0];return ( r );}// =======================================================================function  gc2gd (flatgci, altkmi )/*        geocentric latitude to geodetic latitudeInput:flatgc    geocentric latitude deg.altkm     altitide in kmouput:flatgd    geodetic latitude in deg*/{var dtr   = Math.PI/180.0;var rtd   = 1/dtr;var  flatgd,flatgc,altkmvar  rrnrm = new Array(3)var  re,rn,ecc, esq;var  slat,clat,tlatvar  altnow,ratiogeodGBL();flatgc=  Number(flatgci);altkm =  Number(altkmi);ecc   =  EARTH_Ecc;esq   =  ecc*ecc;//             approximation by stages
//             1st use gc-lat as if is gd, then correct alt dependencealtnow  =  altkm;rrnrm   =  radcur (flatgc);rn      =  rrnrm[1];ratio   = 1 - esq*rn/(rn+altnow);tlat    = Math.tan(dtr*flatgc) / ratio;flatgd  = rtd * Math.atan(tlat);//        now use this approximation for gd-lat to get rn etc.rrnrm   =  radcur ( flatgd );rn      =  rrnrm[1];ratio   =  1  - esq*rn/(rn+altnow)tlat    =  Math.tan(dtr*flatgc)/ratio;flatgd  =  rtd * Math.atan(tlat);return  flatgd}// =======================================================================function  gd2gc (flatgdi, altkmi )/*        geodetic latitude to geocentric latitudeInput:flatgd    geodetic latitude deg.altkm     altitide in kmouput:flatgc    geocentric latitude in deg*/{var dtr   = Math.PI/180.0;var rtd   = 1/dtr;var  flatgc,flatgd,altkmvar  rrnrm = new Array(3)var  re,rn,ecc, esq;var  slat,clat,tlatvar  altnow,ratiogeodGBL();flatgd=  Number(flatgdi);altkm =  Number(altkmi);ecc   =  EARTH_Ecc;esq   =  ecc*ecc;altnow  =  altkm;rrnrm   =  radcur (flatgd);rn      =  rrnrm[1];ratio   = 1 - esq*rn/(rn+altnow);tlat    = Math.tan(dtr*flatgd) * ratio;flatgc  = rtd * Math.atan(tlat);return  flatgc;}// =======================================================================function  llenu ( flati,floni)/*        latitude longitude to east,north,up unit vectorsinput:flat      latitude in degees N[ gc -> gc enu,  gd usual enu ]flon      longitude in degrees Eoutput:enu[3[3]]  packed 3-unit vectors / each a 3 vector*/{var  flat,flon;var  dtr,clat,slat,clon,slon  ;var  ee = new Array(3);var  en = new Array(3);var  eu = new Array(3);var  enu = new Array(3);var dtr   = Math.PI/180.0//             --------------------------------flat = Number(flati);flon = Number(floni);clat = Math.cos(dtr*flat);slat = Math.sin(dtr*flat);clon = Math.cos(dtr*flon);slon = Math.sin(dtr*flon);ee[0]  =  -slon;ee[1]  =   clon;ee[2]  =   0.0;en[0]  =  -clon * slat;en[1]  =  -slon * slat;en[2]  =          clat;eu[0]  =   clon * clat;eu[1]  =   slon * clat;eu[2]  =          slat;enu[0] =  ee;enu[1] =  en;enu[2] =  eu;return  enu ;}// =======================================================================function llhxyz (flati,floni, altkmi )/*        lat,lon,height to xyz vectorinput:flat      geodetic latitude in degflon      longitude in degaltkm     altitude in kmoutput:returns vector x 3 long ECEF in km*/{var  dtr =  Math.PI/180.0;var  flat,flon,altkm;var  clat,clon,slat,slon;var  rrnrm = new Array(3);var  rn,esq;var  x,y,z;var  xvec = new Array(3);geodGBL();flat  = Number(flati);flon  = Number(floni);altkm = Number(altkmi);clat = Math.cos(dtr*flat);slat = Math.sin(dtr*flat);clon = Math.cos(dtr*flon);slon = Math.sin(dtr*flon);rrnrm  = radcur (flat);rn     = rrnrm[1];re     = rrnrm[0];ecc    = EARTH_Ecc;esq    = ecc*eccx      =  (rn + altkm) * clat * clon;y      =  (rn + altkm) * clat * slon;z      =  ( (1-esq)*rn + altkm ) * slat;xvec[0]  = x;xvec[1]  = y;xvec[2]  = z;return  xvec ;}// =======================================================================function xyzllh ( xvec )/*        xyz vector  to  lat,lon,heightinput:xvec[3]   xyz ECEF locationoutput:llhvec[3] with componentsflat      geodetic latitude in degflon      longitude in degaltkm     altitude in km*/{var  dtr =  Math.PI/180.0;var  flatgc,flatn,dlat;var  rnow,rp;var  x,y,z,p;var  tangc,tangd;var  testval,kount;var  rn,esq;var  clat,slat;var  rrnrm = new Array(3);var  flat,flon,altkm;var  llhvec = new Array(3);geodGBL();esq    =  EARTH_Esq;x      = xvec[0];y      = xvec[1];z      = xvec[2];x      = Number(x);y      = Number(y);z      = Number(z);rp     = Math.sqrt ( x*x + y*y + z*z );flatgc = Math.asin ( z / rp )/dtr;testval= Math.abs(x) + Math.abs(y);if ( testval < 1.0e-10){flon = 0.0 }else{flon = Math.atan2 ( y,x )/dtr } if (flon < 0.0 )  { flon = flon + 360.0 }p      =  Math.sqrt( x*x + y*y );//             on pole special caseif ( p < 1.0e-10 ){  flat = 90.0if ( z < 0.0 ) { flat = -90.0 }altkm = rp - rearth(flat);llhvec[0]  = flat;llhvec[1]  = flon;llhvec[2]  = altkm;return  llhvec;}//        first iteration, use flatgc to get altitude
//        and alt needed to convert gc to gd lat.rnow  =  rearth(flatgc);altkm =  rp - rnow;flat  =  gc2gd (flatgc,altkm);rrnrm =  radcur(flat);rn    =  rrnrm[1];for ( var kount = 0; kount< 5 ; kount++ ){slat  =  Math.sin(dtr*flat);tangd =  ( z + rn*esq*slat ) / p;flatn =  Math.atan(tangd)/dtr;dlat  =  flatn - flat;flat  =  flatn;clat  =  Math.cos( dtr*flat );rrnrm =  radcur(flat);rn    =  rrnrm[1];altkm =  (p/clat) - rn;if ( Math.abs(dlat) < 1.0e-12 ) { break }}llhvec[0]  = flat;llhvec[1]  = flon;llhvec[2]  = altkm;return  llhvec ;}// =======================================================================//-->

C++实现

就是根据将JS代码转成C艹.

GPS 转笛卡尔

void GPSTransform::gpsPoint2DescartesPoint(const double latitude, const double longitude, const double altitude, double &x, double &y, double &z)
{//wgs84 WGS84 Earth Constantsdouble wgs84a = 6378.137;double wgs84f = 1.0 / 298.257223563;double wgs84b = wgs84a * (1.0 - wgs84f);//earthcondouble f = 1 - wgs84b / wgs84a;double eccsq = 1 - (wgs84b* wgs84b) / (wgs84a * wgs84a);double ecc = sqrt(eccsq);double esq = ecc * ecc;//llhxyzdouble dtr = M_PI / 180.0;//qDebug() << dtr << gpsPoint.latitude << endl;double clat = cos(dtr * latitude);double slat = sin(dtr * latitude);double clon = cos(dtr * longitude);double slon = sin(dtr * longitude);//qDebug() << clat << slon << endl;//radcur compute the radii at the geodetic latitude lat (in degrees)double dsq = 1.0 - eccsq * slat *slat;double d = sqrt(dsq);//qDebug() << d;double rn = wgs84a / d;double rm = rn * (1.0 - eccsq) / dsq;double rho = rn * clat;double zz = (1.0 - eccsq) * rn *slat;double rsq = rho * rho + zz*zz;double r = sqrt(rsq);x = (rn + altitude) * clat * clon;y = (rn + altitude) * clat * slon;z = ((1 - esq)*rn + altitude) * slat;
}

运行结果

验证网站

Geodetic to Cartesian Converter - http://www.apsalin.com/convert-geodetic-to-cartesian.aspx

atitude,Longitude,Height to/from ECEF (X,Y,Z) - http://www.oc.nps.edu/oc2902w/coord/llhxyz.htm

GPS坐标转笛卡尔坐标相关推荐

四旋翼利用mavros进行GPS坐标指点飞行
先介绍一般px4飞控的xyz坐标指点飞行: 利用mavros的 /mavros/setpoint_raw/local 话题可以发送东北天(ENU)坐标给px4飞控进行指点飞行.ENU坐标原点在起飞点, ...
gps两点距离 php,PHP应用：PHP计算百度地图两个GPS坐标之间距离的方法
<PHP应用:PHP计算百度地图两个GPS坐标之间距离的方法>要点: 本文介绍了PHP应用:PHP计算百度地图两个GPS坐标之间距离的方法,希望对您有用.如果有疑问,可以联系我们. 本文实 ...
[python] 从GPS坐标获取国家名
目标比较明确,就是从GPS坐标得到它所在的国家. 网上可以找的比较典型的解决方案是利用一些网站(例如Google)的webservice来完成这个任务,但是这些解决方案有一个比较大的弱点,就是这些we ...
PHP+百度地图API+JAVASCRIPT实现GPS坐标与百度坐标转换的实例
PHP+百度地图API+JAVASCRIPT实现GPS坐标与百度坐标转换的实例原文:PHP+百度地图API+JAVASCRIPT实现GPS坐标与百度坐标转换的实例 <!--小幅的坐标转换点位程 ...
php 地图距离,PHP计算百度地图两个GPS坐标之间距离的方法
这篇文章主要介绍了PHP计算百度地图两个GPS坐标之间距离的方法,是针对百度地图接口开发的典型应用,需要的朋友可以参考下本文实例讲述了PHP计算百度地图两个GPS坐标之间距离的方法.分享给大家供大家 ...
百度地图批量转换 GPS坐标转百度地图坐标问题
百度地图的官方网址官方批量转换的demo 花了几天时间了解了一下百度地图,之前是后端的一个小伙伴在负责,他跟我吐槽这是前端的东西,让我来写(之前他们老大交给他了,我也以为是后端的任务(๑′ᴗ‵๑). ...
计算坐标点的距离计算机公式,计算两个GPS坐标点的距离
原标题:计算两个GPS坐标点的距离在日常开发中,我们难免要计算两个左边之间的距离,但是地图软件api的接口普遍要求我们必须要先将坐标点传递到他们服务器,然后计算出一个距离返还给我们,使用起来太不方便 ...
java gps 距离计算_Java教程之地图中计算两个GPS坐标点的距离
原标题:Java教程之地图中计算两个GPS坐标点的距离在日常开发中,我们难免要计算两个左边之间的距离,但是地图软件api的接口普遍要求我们必须要先将坐标点传递到他们服务器,然后计算出一个距离返还给我 ...
GPS坐标对应地图坐标偏移问题
各个地图的坐标和GPS坐标计算方式是不同的,因此有时候需要进行坐标转化再进行定位,坐标转化算法如下: //定义一些常量 var x_PI = 3.14159265358979324 * 3000.0 ...

GPS坐标转笛卡尔坐标

提要

js实现

C++实现

验证网站

GPS坐标转笛卡尔坐标相关推荐

最新文章

热门文章