blog
地理坐标系与地固坐标系转换(自编)
由geocentric与geodetic之间的转换是大地测量中的基本工具，后者转为前者是显示表达式，反正则是非线性方程，原则上要求迭代处理。 在sofa中已经提供了两者的转换，可以将其编译为常用小工具。 这里用C++自行编写软件进行转换. 这里算法参考见"satellite geodesy"第二版2.1.4节。 注意：因浏览器显示问题，实际代码有出入。 公式如下 \[{e^2} = 2f - {f^2}\] \[{e^2} = \frac{{{a^2} - {b^2}}}{{{a^2}}}\] \[\left[ \begin{gathered} X \hfill \\ Y \hfill \\ Z \hfill \\ \end{gathered} \right] = \left[ \begin{gathered} \left( {\bar N + h} \right)\cos \varphi \cos \lambda \hfill \\ \left( {\bar N + h} \right)\cos \varphi \sin \lambda \hfill \\ \left( {\left( {1 - {e^2}} \right)\bar N + h} \right)\sin \varphi \hfill \\ \end{gathered} \right]\] \[\bar N = \frac{a}{{\sqrt {1 - {e^2}{{\sin }^2}\varphi } }} = \frac{a}{{\sqrt {1 - f\left( {2 - f} \right){{\sin }^2}\varphi } }}\] 反之 纬度初值为 \[\varphi = \arctan \left[ {\frac{Z}{{\sqrt {{X^2} + {Y^2}} }}{{\left( {1 - {e^2}} \right)}^{ - 1}}} \right]\] \[\bar N = \frac{a}{{\sqrt {1 - {e^2}{{\sin }^2}\varphi } }} = \frac{a}{{\sqrt {1 - f\left( {2 - f} \right){{\sin }^2}\varphi } }}\] \[h = \frac{{\sqrt {{X^2} + {Y^2}} }}{{\cos \varphi }} - \bar N\] \[\varphi = \arctan \left[ {\frac{Z}{{\sqrt {{X^2} + {Y^2}} }}{{\left( {1 - {e^2}\frac{{\bar N}}{{\bar N + h}}} \right)}^{ - 1}}} \right]\] \[\lambda = \arctan \frac{Y}{X}\] 其中大地坐标并不是直角坐标的显函数，因此需要迭代处理，详见代码。 注意在收敛过程中，有些书上只用前后phi值作为收敛准则，建议同时判断lamda。只有当两者都收敛时 停止迭代，防止某一变量巧合收敛，而另一变量并未收敛情况。
makefile文件为
#------------------------------------------------------ # Versions and Changes : # v1.0----2015-7-18 # makefile #-------------------------------------------------------- # 如果有很多文件，则在obj变量后加入 # 如果有多个 makefile文件 # 可以用 gmake -f makefilename #-------------------------------------------------------- # Copyrigt(C) : Shanghai Astronomical Observatory, CAS # (All rights reserved) 2015 #------------------------------------------------------- cc= g++ #cflag= -g -c obj = gc2gd.cpp coordinates.cpp target=gc2gdnew.exe $(target) : $(obj) $(cc) -g -o $(target) $(obj) .PHONY : clean clean : rm -f $(obj) $(target)
基本函数为coordinates.cpp . 其头文件 coordinates.h为
#ifndef coord_h
#define coord_h
#include
#include
#include
using namespace std;
void ellipsoidal2Cartesian(int n ,double phi,double lamda, double height,vector &XYZ);
void cartesian2Ellipsoidal(int n, const vector &XYZ, double &phi,double &lamda,double &height);
void eform ( int n, double &a, double &f );
#endif
#include "coordinates.h"
/*------------------------------------------------------
Created : Song Yezhi 2023-4-1 2:47
transformation equation from the geographic ellipsoidal coordinates
to Cartesian coordinates xyz
Ref : chapter 2.1.4
--------------------------------------------------------
Input Parameters :
phi (arc)
lamda (arc)
height (meter)
Output Parameters :
XYZ (meter)
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2023
-------------------------------------------------------*/
void ellipsoidal2Cartesian(int n ,double phi,double lamda,
double height,vector &XYZ)
{
double a,f ;
eform(n,a ,f);
double e_square = 2.0*f - f*f ;
double N_bar = a/sqrt(1.0- e_square*sin(phi)*sin(phi)) ;
XYZ[0] = (N_bar+height)*cos(phi)*cos(lamda);
XYZ[1] = (N_bar+height)*cos(phi)*sin(lamda);
XYZ[2] = ((1-e_square)*N_bar +height)*sin(phi);
}
/*------------------------------------------------------
Created : Song Yezhi 2023-4-1 2:47
transformation equation from Cartesian coordinates xyz to the geographic ellipsoidal coordinates
Ref : chapter 2.1.4
GNSS Data Processing gAge slides
--------------------------------------------------------
Input Parameters :
XYZ (meter)
Output Parameters :
phi (arc)
lamda (arc)
height (meter)
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2023
-------------------------------------------------------*/
void cartesian2Ellipsoidal(int n, const vector &XYZ,
double &phi,double &lamda,double &height)
{
const double tol = 1e-5 ;
double a,f ;
eform(n,a ,f);
double e_square ;
e_square = 2.0*f - f*f ;
double N_bar ;
double tmp1,tmp2 ;
double tmp3,tmp4;
double height0,phi0 ;
int i ;
double p ;
p = sqrt(XYZ[0]*XYZ[0]+XYZ[1]*XYZ[1]);
phi = atan(XYZ[2]/p/(1-e_square));
//the initial value of phi
for(i =0 ;i<50 ;i++ )
{
N_bar = a/sqrt(1.0- e_square*sin(phi)*sin(phi)) ;
height = p/cos(phi) - N_bar;
tmp1 = XYZ[2]/p ;
tmp2 = 1.0- e_square*N_bar/(N_bar+ height);
phi = atan(tmp1/tmp2);
tmp3 = abs(height -height0);
tmp4 = abs((phi-phi0)*a );
// 同时对经纬度判断，防止某一个变量巧合收敛
// 因浏览器显示原因, 使用该程序时需把这里的"小于"改为 C语言中的小于号.
if ((tmp3 小于 tol )&&(tmp4 小于 tol))
{
break ;
}
phi0 = phi ;
height0 = height;
}
lamda = atan2(XYZ[1],XYZ[0]);
}
/*------------------------------------------------------
Created : Song Yezhi 2023-4-1 2:50
** 1) The identifier n is a number that specifies the choice of
** reference ellipsoid. The following are supported:
** n ellipsoid
**
** 1 WGS84
** 2 GRS80
** 3 WGS72
the parameter reference is IAU sofa
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2023
-------------------------------------------------------*/
void eform ( int n, double &a, double &f )
{
/* Look up a and f for the specified reference ellipsoid. */
switch ( n )
{
case 1:
a = 6378137.0;
f = 1.0 / 298.257223563;
break;
case 2:
a = 6378137.0;
f = 1.0 / 298.257222101;
break;
case 3:
a = 6378135.0;
f = 1.0 / 298.26;
break;
default:
/* Invalid identifier. */
a = 0.0;
f = 0.0;
}
}
/*------------------------------------------------------
Versions and Changes :
Transform geodetic coordinates to geocentric using WGS84
reference ellipsoid.
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Author : 宋叶志
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2023
-------------------------------------------------------*/
#include
#include
#include "coordinates.h"
using namespace std;
int main(int argc, char **argv)
{
if (argc == 1)
{
printf("%s","\n");
printf("%s","Transform geodetic coordinates to geocentric \n");
printf("%s","the reference ellipsoid is WGS84 \n");
printf("%s","usage : \n");
printf("%s","gc2gd longitude latitude height (in degree and meters) \n");
return 0;
}
//double XYZ[3];
vector XYZ(3);
double elong,phi,height;
int itmp;
const double pi = 3.141592653589793238462643 ;
sscanf(argv[1],"%lf",&elong);
sscanf(argv[2],"%lf",&phi);
sscanf(argv[3],"%lf",&height);
printf("%s","\n");
printf("%s","longitude latitude and height = (in degree and meters) \n");
printf("%18lf %18lf %18lf \n \n", elong,phi,height);
ellipsoidal2Cartesian ( 1, elong*pi/180.0, phi*pi/180.0, height, XYZ);
printf("%s","longitude latitude and height = (in arc and meters) \n");
printf("%18lf %18lf %18lf \n \n", elong*pi/180.0 ,phi*pi/180.0 ,height);
printf("%s","XYZ = (in meters) \n");
printf("%18.4lf %18.4lf %18.4lf \n \n", XYZ[0],XYZ[1],XYZ[2]);
return 0;
}
/*------------------------------------------------------
Versions and Changes :
Transform geocentric coordinates to geodetic using the WGS84 ellipsoid
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Author : 宋叶志
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2023
-------------------------------------------------------*/
#include
#include
#include "coordinates.h"
using namespace std;
int main(int argc, char **argv)
{
if (argc == 1)
{
printf("%s","\n");
printf("%s","Transform geocentric coordinates to geodetic \n");
printf("%s","the reference ellipsoid is WGS84 \n");
printf("%s","usage : \n");
printf("%s","gc2gd x y z (in meters) \n");
return 0;
}
//double XYZ[3];
vector XYZ(3);
sscanf(argv[1],"%lf",&XYZ[0]);
sscanf(argv[2],"%lf",&XYZ[1]);
sscanf(argv[3],"%lf",&XYZ[2]);
printf("%s","\n");
printf("%s","XYZ = (in meters) \n");
printf("%18lf %18lf %18lf \n \n", XYZ[0],XYZ[1],XYZ[2]);
double elong,phi,height;
//cartesian2Ellipsoidal( 1, XYZ, &elong, &phi, &height );
cartesian2Ellipsoidal(1, XYZ, phi, elong,height) ;
const double pi = 3.141592653589793238462643 ;
printf("%s","longitude latitude and height = (in arc and meters) \n");
printf("%18.10lf %18.10lf %18.4lf \n \n", elong,phi,height);
printf("%s","longitude latitude and height = (in degree and meters) \n");
printf("%18.10lf %18.10lf %18.4lf \n \n", elong*180.0/pi ,phi*180.0/pi ,height);
return 0;
}