企业网站管理系统免费,东莞营销外包公司,有什么做动画的网站,做网站营销发布文章经纬度与ECEF直角坐标的基本换算
我们目前最常用的全球坐标系是WGS-84坐标系#xff0c;各种手机、地图基本用经纬度来标记位置。然而#xff0c;经纬度对于空间的计算是很复杂的#xff0c;需要很多三角函数操作。平面直角坐标系利用向量的运算#xff0c;可以非常方便的…经纬度与ECEF直角坐标的基本换算
我们目前最常用的全球坐标系是WGS-84坐标系各种手机、地图基本用经纬度来标记位置。然而经纬度对于空间的计算是很复杂的需要很多三角函数操作。平面直角坐标系利用向量的运算可以非常方便的计算角度、距离等参数在实际应用中往往作为中间计算的工具。
目前用到的很多GIS、遥感与测绘工具里都有这种功能比如利用 libproj、RTK等工具直接进行转换。 文章目录 经纬度与ECEF直角坐标的基本换算1. WGS-84 坐标系2 从经纬度转换到ECEF3 从 ECEF 到经纬度(1) 求取经度(2) 纬度、海拔计算 4 代码与工程 1. WGS-84 坐标系
WGS-84 坐标系是目前使用最为广泛的全球坐标系统。其示意图如下 WGS-84坐标系里地球等效为一个椭圆直角坐标系原点O位于地球的中心。
坐标轴 Z ⃗ O N ⃗ \vec{Z} \vec{ON} Z ON 指向北极点。坐标轴 X ⃗ O P ⃗ \vec{X} \vec{OP} X OP 指向格林威治0度经线与赤道面的交点P。坐标轴 Y ⃗ X ⃗ ⊗ Z ⃗ \vec Y\vec X\otimes \vec Z Y X ⊗Z 呈右手螺旋关系。
当前位置 T 的经纬度定义
从高度 H 开始定义
T是一个距离参考海平面椭球面高度为 H ⃗ t T ⃗ \vec{H} \vec{tT} H tT 的点。t是点T在椭球面上的正射投影矢量 t T ⃗ \vec{tT} tT 是椭球面上t点切面的法向量。需要注意的是和正球面不同法向量 t T ⃗ \vec{tT} tT 的反向延长线与Z轴的交点是 O’和原点O并不重合。
对于经度 θ \theta θ 的定义是这样的
0度经度大圆 N O P ⌢ \stackrel\frown{NOP} NOP⌢是由北极点 N, 地心 O, 格林威治经线与赤道面的交点 P 组成的平面。当前经度大圆 N O ′ T ⌢ \stackrel\frown{NOT} NO′T⌢是由北极点 N, 交点O’, 当前位置T组成的平面。O、 O’ 、t、 T 都在当前经度面上直线 O’T 与赤道面的交点为u矢量 O u ⃗ \vec{Ou} Ou 与赤道的交点是Q。O’、u、t、T四点共线。经度 θ ∡ P O Q \theta\measuredangle {POQ} θ∡POQ
对于纬度 φ \varphi φ 的定义是这样的
椭球的切平面的法线与Z的交点 O’ 和原点并不重合。纬度定义为法向量 u T ⃗ \vec{uT} uT 与赤道面的夹角。在北半球为正南半球为负纬度 φ ∡ T u Q ∡ O u O ′ \varphi\measuredangle {TuQ}\measuredangle {OuO} φ∡TuQ∡OuO′
2 从经纬度转换到ECEF
ECEF 坐标系有利于利用矢量计算简化角度、距离的求取。
已知:
赤道半径 a 6378137 米对应上图 O P ⃗ \vec{OP} OP 的长度;椭球扁率 f ( a − b ) / a f (a-b)/a f(a−b)/a 1.0/298.257223563。扁率就是衡量椭球是不是很“扁”正球的扁率为0.椭球偏心率 e a 2 − b 2 / a e{\sqrt {a^2-b^2} }/{a} ea2−b2 /a 0.0818191908426e和f的关系为 e 2 f ( 2 − f ) e^2f(2-f) e2f(2−f)极地半径 b 6356752.3142 米对应上图 O N ⃗ \vec{ON} ON 的长度;
设经度 θ \theta θ 纬度 φ \varphi φ 海拔 H则有 r ∣ R ⃗ ∣ ∣ O ′ t ⃗ ∣ a 1 − e 2 sin 2 φ a 1 − f ( 2 − f ) sin 2 φ r \left | \vec R\right |\left | \vec{Ot} \right | \frac{ a}{\sqrt{1-e^2 \sin^2\varphi }} \frac{ a}{\sqrt{1-f(2-f) \sin^2\varphi }} r R O′t 1−e2sin2φ a1−f(2−f)sin2φ a
经纬度极坐标可以利用矢量 R ⃗ \vec R R 构造一个Z轴平移了的正球进行极坐标运算。这个正球的球心是 O’。由于x,y没有平移故而两个坐标系的x、y是重合的。Z坐标z’与z的关系是平移平移量为OO’的长度 ∣ O O ′ ⃗ ∣ r e 2 sin φ \left | \vec {OO} \right | re^2\sin \varphi OO′ re2sinφ
下图是把经度大圆 N O Q ⌢ \stackrel\frown {NOQ} NOQ⌢切割出来观察 在正球模型下根据极坐标的基本定义可以直接计算x,y: x ( r H ) ⋅ cos φ ⋅ c o s θ x\left ( r H \right ) \cdot \cos \varphi \cdot cos \theta x(rH)⋅cosφ⋅cosθ y ( r H ) ⋅ cos φ ⋅ s i n θ y\left ( r H \right ) \cdot \cos \varphi \cdot sin \theta y(rH)⋅cosφ⋅sinθ
Z的坐标要先计算后再平移即可得到: z ′ ( r H ) ⋅ s i n φ z\left ( r H \right ) \cdot sin \varphi z′(rH)⋅sinφ z z ′ − r e 2 sin φ zz - re^2\sin \varphi zz′−re2sinφ
带入后可写成如下等效形式 z ( r ( 1 − e 2 ) H ) sin φ z(r(1-e^2)H) \sin \varphi z(r(1−e2)H)sinφ z ( r ( 1 − f ( 2 − f ) ) H ) sin φ ( r ( 1 − f ) 2 H ) sin φ z(r(1-f(2-f))H) \sin \varphi (r(1-f)^2H) \sin \varphi z(r(1−f(2−f))H)sinφ(r(1−f)2H)sinφ
含有扁率的化简是基于扁率 f ( a − b ) / a 1 − b / a f (a-b)/a 1 - b/a f(a−b)/a1−b/a 计算得到的。
通过上面的处理极坐标的 T ( θ , φ , H ) T(\theta, \varphi,H) T(θ,φ,H) 便转换为平面直角坐标 T ( x , y , z ) T(x,y,z) T(x,y,z) 。
3 从 ECEF 到经纬度
对上述计算过程而言逆向运算可以立刻得到经度却不能解析得到纬度。纬度需要进行迭代。
主要原因是 ∣ O O ′ ⃗ ∣ r e 2 sin φ \left | \vec {OO} \right | re^2\sin \varphi OO′ re2sinφ
的具体取值不知道导致的。
(1) 求取经度
由于 x / y tan θ x/y\tan \theta x/ytanθ
故而 θ tan − 1 ( x / y ) \theta\tan^{-1} (x/y) θtan−1(x/y)
注意处理好象限问题即可得到经度。
(2) 纬度、海拔计算
对于纬度则比较复杂。求解三角方程 ( r H ) sin φ z ∣ O O ′ ⃗ ∣ (r H) \sin \varphi z \left | \vec {OO} \right | (rH)sinφz OO′
由于 H 不知道r 也不知道加入第二条件相切H向量与r共线后得到的是二元四次三角方程很难解析求取。
可以首先假设 ∣ O O ′ ⃗ ∣ 0 \left | \vec {OO} \right | 0 OO′ 0, 求取一个粗略的纬度并得到 d利用d反过来更新其他数值进行迭代。这种迭代收敛的前提是因为每次计算出的纬度一定小于真实的纬度这是由三角关系决定的。
迭代的效果是不断地把原点向真实的原点移动直到z达到一个很小的误差。迭代方法
设起始偏移 d ∣ O O ′ ⃗ ∣ 0 d \left | \vec {OO} \right | 0 d OO′ 0迭代开始 R ′ r H x 2 y 2 ( z d ) 2 RrH \sqrt{x^2y^2(zd)^2} R′rHx2y2(zd)2 φ ′ sin − 1 z d R ′ \varphi \sin^{-1} {\frac {zd}{R}} φ′sin−1R′zd r ′ a 1 − e 2 sin 2 φ ′ r \frac{ a}{\sqrt{1-e^2 \sin^2\varphi }} r′1−e2sin2φ′ a d ′ r ′ e 2 sin φ ′ d re^2\sin \varphi d′r′e2sinφ′ z ′ R ′ s i n φ ′ − d ′ zRsin \varphi-d z′R′sinφ′−d′ d d ′ d d dd′
当 e r r ∣ z − z ′ ∣ ε err \left | z-z \right | \varepsilon err∣z−z′∣ε 时停止迭代。此时 φ sin − 1 z d x 2 y 2 ( z d ) 2 \varphi \sin^{-1} {\frac {zd}{\sqrt{x^2y^2(zd)^2}}} φsin−1x2y2(zd)2 zd H R ′ − r ′ H R-r HR′−r′
4 代码与工程
代码与工程参考
https://gitcode.net/coloreaglestdio/geocalc/-/blob/master/geocalc.h
核心转换逻辑如下
/*!* GEOCalc 实现了基本的WGS-84坐标系的计算。* by CEStdio* 1997-2023*
*/
#include cassert
#include cmath
namespace CES_GEOCALC {inline const double a 6378137;
inline const double pi 3.14159265358979323846;
inline const double f 1.0 / 298.257223563;
inline const double e sqrt(f * (2 - f));
inline const double es f * (2 - f);
inline const double deg2rad pi / 180.0;
inline const double rad2deg 180.0 / pi;/*!* \brief lla2ecef 经纬度坐标到ECEF* \param lla 纬经高(默认)/经纬高, 量纲是度(默认)/弧度、米* \param ecef xyz量纲是米* \param rad 角度量纲开关false 是度true 是弧度* \param latfirst 经纬度顺序false 是经度\纬度\高度true 是纬度\经度\高度
*/
inline void lla2ecef(const double lla[/*3*/],double ecef[/*3*/],const bool rad false,const bool latfirst true)
{const double lat latfirst ? (rad ? (lla[0]) : (lla[0] * deg2rad)): (rad ? (lla[1]) : (lla[1] * deg2rad));const double lon latfirst ? (rad ? (lla[1]) : (lla[1] * deg2rad)): (rad ? (lla[0]) : (lla[0] * deg2rad));const double H lla[2];const double r a / sqrt(1 - es * sin(lat) * sin(lat));const double d r * es * sin(lat);ecef[0] (r H) * cos(lat) * cos(lon);ecef[1] (r H) * cos(lat) * sin(lon);ecef[2] (r H) * sin(lat) - d;
}/*!* \brief ecef2lla ECEF到经纬度坐标* \param ecef xyz量纲是米* \param lla 纬经高(默认)/经纬高, 量纲是度(默认)/弧度、米* \param maxiter 最大迭代次数* \param piter 迭代次数输出可以为null* \param eps Z误差门限* \param rad 角度量纲开关false 是度true 是弧度* \param latfirst 经纬度顺序false 是经度\纬度\高度true 是纬度\经度\高度* \return 迭代收敛标志
*/
inline bool ecef2lla(const double ecef[/*3*/],double lla[/*3*/],const int maxiter 32,int * piter nullptr,const double eps 1e-10,const bool rad false,const bool latfirst true)
{const double x ecef[0], y ecef[1], z ecef[2] 0 ? -ecef[2] : ecef[2];const int south ecef[2] 0 ? -1 : 1;double d 0;double err 0;int iter 0;double Rp 0, sinphi 0, rp 0, dp 0, zp 0;do {Rp sqrt(x * x y * y (z d) * (z d));sinphi (z d) / Rp;rp a / sqrt(1 - es * sinphi * sinphi);dp rp * es * sinphi;zp Rp * sinphi - dp;d dp;assert(z zp);err z - zp; //z always zpiter;} while (err eps iter maxiter);const double lat asin(sinphi) * south;const double lon atan2(y, x);const double H Rp - rp;lla[latfirst?0:1] lat * (rad?1:rad2deg);lla[latfirst?1:0] lon * (rad?1:rad2deg);lla[2] H;if (piter)*piter iter;return err eps;
}}; // namespace CES_GEOCALC
对各种边界和迭代次数进行测试, 以确定迭代的收敛性:
using namespace CES_GEOCALC;for (double lon -180; lon 180; lon60){for (double lat -90; lat 90; lat 15){double LLA[] {lat,lon,10000};double ECEF[] {0, 0, 0};double LLA2[] {0, 0, 0};lla2ecef(LLA, ECEF);int iter 0;ecef2lla(ECEF,LLA2,32,iter);printf(LAT%12.7lf,LON%12.7lf,dLAT%10.7lf,dLON%10.7lf,dALT%5.3lf, iter %d\n,lat,lon,LLA[0]-LLA2[0], LLA[1]-LLA2[1], LLA[2]-LLA2[2],iter);}}输出
LAT -90.0000000,LON-180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 2
LAT -75.0000000,LON-180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT -60.0000000,LON-180.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT -45.0000000,LON-180.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 7
LAT -30.0000000,LON-180.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT -15.0000000,LON-180.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 0.0000000,LON-180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 1
LAT 15.0000000,LON-180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 30.0000000,LON-180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 45.0000000,LON-180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 7
LAT 60.0000000,LON-180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT 75.0000000,LON-180.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT 90.0000000,LON-180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 2
LAT -90.0000000,LON-120.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 2
LAT -75.0000000,LON-120.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 6
LAT -60.0000000,LON-120.0000000,dLAT-0.0000000,dLON-0.0000000,dALT0.000, iter 6
LAT -45.0000000,LON-120.0000000,dLAT-0.0000000,dLON-0.0000000,dALT0.000, iter 7
LAT -30.0000000,LON-120.0000000,dLAT-0.0000000,dLON-0.0000000,dALT0.000, iter 8
LAT -15.0000000,LON-120.0000000,dLAT-0.0000000,dLON-0.0000000,dALT0.000, iter 8
LAT 0.0000000,LON-120.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 1
LAT 15.0000000,LON-120.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 8
LAT 30.0000000,LON-120.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 8
LAT 45.0000000,LON-120.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 7
LAT 60.0000000,LON-120.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 6
LAT 75.0000000,LON-120.0000000,dLAT-0.0000000,dLON-0.0000000,dALT0.000, iter 6
LAT 90.0000000,LON-120.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 2
LAT -90.0000000,LON -60.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 2
LAT -75.0000000,LON -60.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 6
LAT -60.0000000,LON -60.0000000,dLAT-0.0000000,dLON-0.0000000,dALT0.000, iter 6
LAT -45.0000000,LON -60.0000000,dLAT-0.0000000,dLON-0.0000000,dALT0.000, iter 7
LAT -30.0000000,LON -60.0000000,dLAT-0.0000000,dLON-0.0000000,dALT0.000, iter 8
LAT -15.0000000,LON -60.0000000,dLAT-0.0000000,dLON-0.0000000,dALT0.000, iter 8
LAT 0.0000000,LON -60.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 1
LAT 15.0000000,LON -60.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 8
LAT 30.0000000,LON -60.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 8
LAT 45.0000000,LON -60.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 7
LAT 60.0000000,LON -60.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 6
LAT 75.0000000,LON -60.0000000,dLAT-0.0000000,dLON-0.0000000,dALT0.000, iter 6
LAT 90.0000000,LON -60.0000000,dLAT 0.0000000,dLON-0.0000000,dALT0.000, iter 2
LAT -90.0000000,LON 0.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 2
LAT -75.0000000,LON 0.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT -60.0000000,LON 0.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT -45.0000000,LON 0.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 7
LAT -30.0000000,LON 0.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT -15.0000000,LON 0.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 0.0000000,LON 0.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 1
LAT 15.0000000,LON 0.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 30.0000000,LON 0.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 45.0000000,LON 0.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 7
LAT 60.0000000,LON 0.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT 75.0000000,LON 0.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT 90.0000000,LON 0.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 2
LAT -90.0000000,LON 60.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 2
LAT -75.0000000,LON 60.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT -60.0000000,LON 60.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT -45.0000000,LON 60.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 7
LAT -30.0000000,LON 60.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT -15.0000000,LON 60.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 0.0000000,LON 60.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 1
LAT 15.0000000,LON 60.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 30.0000000,LON 60.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 45.0000000,LON 60.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 7
LAT 60.0000000,LON 60.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT 75.0000000,LON 60.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT 90.0000000,LON 60.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 2
LAT -90.0000000,LON 120.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 2
LAT -75.0000000,LON 120.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT -60.0000000,LON 120.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT -45.0000000,LON 120.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 7
LAT -30.0000000,LON 120.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT -15.0000000,LON 120.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 0.0000000,LON 120.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 1
LAT 15.0000000,LON 120.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 30.0000000,LON 120.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 45.0000000,LON 120.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 7
LAT 60.0000000,LON 120.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT 75.0000000,LON 120.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT 90.0000000,LON 120.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 2
LAT -90.0000000,LON 180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 2
LAT -75.0000000,LON 180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT -60.0000000,LON 180.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT -45.0000000,LON 180.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 7
LAT -30.0000000,LON 180.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT -15.0000000,LON 180.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 0.0000000,LON 180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 1
LAT 15.0000000,LON 180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 30.0000000,LON 180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 8
LAT 45.0000000,LON 180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 7
LAT 60.0000000,LON 180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT 75.0000000,LON 180.0000000,dLAT-0.0000000,dLON 0.0000000,dALT0.000, iter 6
LAT 90.0000000,LON 180.0000000,dLAT 0.0000000,dLON 0.0000000,dALT0.000, iter 2
