|
|
|
|

大地坐标转换直角坐标的C++代码(高斯克吕格投影)

代码如下://CPP文件#i nclude "CoorTrans.h"////////////////////////////////////////////// Common functions////////////////////////////////////////////double Dms2Rad(double Dms){double Degree, Miniute;double Second;int S...

作者:未知来源:未知|2006年11月25日
代码如下:
//CPP文件

#i nclude   "CoorTrans.h"
////////////////////////////////////////////
//   Common   functions
////////////////////////////////////////////
double   Dms2Rad(double   Dms)
{
double   Degree,   Miniute;
double   Second;
int   Sign;
double   Rad;
if(Dms   >=   0)
Sign   =   1;
else
Sign   =   -1;
Dms   =   fabs(Dms);
Degree   =   floor(Dms);
Miniute   =   floor(fmod(Dms   *   100.0,   100.0));
Second   =   fmod(Dms   *   10000.0,   100.0);
Rad   =   Sign   *   (Degree   +   Miniute   /   60.0   +   Second   /   3600.0)   *   PI   /   180.0;
return   Rad;
}
double   Rad2Dms(double   Rad)
{
double   Degree,   Miniute;
double   Second;
int   Sign;
double   Dms;
if(Rad   >=   0)
Sign   =   1;
else
Sign   =   -1;
Rad   =   fabs(Rad   *   180.0   /   PI);
Degree   =   floor(Rad);
Miniute   =   floor(fmod(Rad   *   60.0,   60.0));
Second   =   fmod(Rad   *   3600.0,   60.0);
Dms   =   Sign   *   (Degree   +   Miniute   /   100.0   +   Second   /   10000.0);
return   Dms;
}
///////////////////////////////////////////////////
//   Definition   of   PrjPoint
///////////////////////////////////////////////////
BOOL   PrjPoint::BL2xy()
{
double   X,   N,   t,   t2,   m,   m2,   ng2;
double   sinB,   cosB;
X   =   A1   *   B   *   180.0   /   PI   +   A2   *   sin(2   *   B)   +   A3   *   sin(4   *   B)   +   A4   *   sin(6   *  
B);
sinB   =   sin(B);
cosB   =   cos(B);
t   =   tan(B);
t2   =   t   *   t;
N   =   a   /   sqrt(1   -   e2   *   sinB   *   sinB);
m   =   cosB   *   (L   -   L0);
m2   =   m   *   m;
ng2   =   cosB   *   cosB   *   e2   /   (1   -   e2);
x   =   X   +   N   *   t   *   ((0.5   +   ((5   -   t2   +   9   *   ng2   +   4   *   ng2   *   ng2)   /   24.0   +   (61   -  
58   *   t2   +   t2   *   t2)   *   m2   /   720.0)   *   m2)   *   m2);
y   =   N   *   m   *   (   1   +   m2   *   (   (1   -   t2   +   ng2)   /   6.0   +   m2   *   (   5   -   18   *   t2   +   t2   *   t
2   +   14   *   ng2   -   58   *   ng2   *   t2   )   /   120.0));
y   +=   500000;
return   TRUE;
}
BOOL   PrjPoint::xy2BL()
{
double   sinB,   cosB,   t,   t2,   N   ,ng2,   V,   yN;
double   preB0,   B0;
double   eta;
y   -=   500000;
B0   =   x   /   A1;
do
{
preB0   =   B0;
B0   =   B0   *   PI   /   180.0;
B0   =   (x   -   (A2   *   sin(2   *   B0)   +   A3   *   sin(4   *   B0)   +   A4   *   sin(6   *   B0)))   /   A1;
eta   =   fabs(B0   -   preB0);
}while(eta   >   0.000000001);
B0   =   B0   *   PI   /   180.0;
B   =   Rad2Dms(B0);
sinB   =   sin(B0);
cosB   =   cos(B0);
t   =   tan(B0);
t2   =   t   *   t;
N   =   a   /   sqrt(1   -   e2   *   sinB   *   sinB);
ng2   =   cosB   *   cosB   *   e2   /   (1   -   e2);
V   =   sqrt(1   +   ng2);
yN   =   y   /   N;
B   =   B0   -   (yN   *   yN   -   (5   +   3   *   t2   +   ng2   -   9   *   ng2   *   t2)   *   yN   *   yN   *   yN   *   yN   /
12.0   +   (61   +   90   *   t2   +   45   *   t2   *   t2)   *   yN   *   yN   *   yN   *   yN   *   yN   *   yN   /   360.0)
*   V   *   V   *   t   /   2;
L   =   L0   +   (yN   -   (1   +   2   *   t2   +   ng2)   *   yN   *   yN   *   yN   /   6.0   +   (5   +   28   *   t2   +   24  
*   t2   *   t2   +   6   *   ng2   +   8   *   ng2   *   t2)   *   yN   *   yN   *   yN   *   yN   *   yN   /   120.0)   /   cosB
;
return   TRUE;
}
BOOL   PrjPoint::SetL0(double   dL0)
{
L0   =   Dms2Rad(dL0);
return   TRUE;
}
BOOL   PrjPoint::SetBL(double   dB,   double   dL)
{
B   =   Dms2Rad(dB);
L   =   Dms2Rad(dL);
B   =   dB;
L   =   dL;
BL2xy();
return   TRUE;
}
BOOL   PrjPoint::GetBL(double   *dB,   double   *dL)
{
*dB   =   Rad2Dms(B);
*dL   =   Rad2Dms(L);
return   TRUE;
}
BOOL   PrjPoint::Setxy(double   dx,   double   dy)
{
x   =   dx;
y   =   dy;
xy2BL();
return   TRUE;
}
BOOL   PrjPoint::Getxy(double   *dx,   double   *dy)
{
*dx   =   x;
*dy   =   y;
return   TRUE;
}
///////////////////////////////////////////////////
//   Definition   of   PrjPoint_IUGG1975
///////////////////////////////////////////////////
PrjPoint_IUGG1975::PrjPoint_IUGG1975()
{
a   =   6378140;
f   =   298.257;
e2   =   1   -   ((f   -   1)   /   f)   *   ((f   -   1)   /   f);
e12   =   (f   /   (f   -   1))   *   (f   /   (f   -   1))   -   1;
A1   =   111133.0047;
A2   =   -16038.5282;
A3   =   16.8326;
A4   =   -0.0220;
}
PrjPoint_IUGG1975::~PrjPoint_IUGG1975()
{
}
///////////////////////////////////////////////////
//   Definition   of   PrjPoint_Krasovsky
///////////////////////////////////////////////////
PrjPoint_Krasovsky::PrjPoint_Krasovsky()
{
a   =   6378245;
f   =   298.3;
e2   =   1   -   ((f   -   1)   /   f)   *   ((f   -   1)   /   f);
e12   =   (f   /   (f   -   1))   *   (f   /   (f   -   1))   -   1;
A1   =   111134.8611;
A2   =   -16036.4803;
A3   =   16.8281;
A4   =   -0.0220;
}
PrjPoint_Krasovsky::~PrjPoint_Krasovsky()
{
}


//H文件


#ifndef   _COORTRANS_H_INCLUDED
#define   _COORTRANS_H_INCLUDED
#i nclude  
const   double   PI   =   3.14159265353846;
class   PrjPoint
{
public:
double   L0;   //   中央子午线经度
double   B,   L;   //   大地坐标
double   x,   y;   //   高斯投影平面坐标
public:
BOOL   BL2xy();
BOOL   xy2BL();
protected:
double   a,   f,   e2,   e12;   //   基本椭球参数
double   A1,   A2,   A3,   A4;   //   用于计算X的椭球参数
public:
BOOL   SetL0(double   dL0);
BOOL   SetBL(double   dB,   double   dL);
BOOL   GetBL(double   *dB,   double   *dL);
BOOL   Setxy(double   dx,   double   dy);
BOOL   Getxy(double   *dx,   double   *dy);
};
class   PrjPoint_Krasovsky   :   virtual   public   PrjPoint
{
public:
PrjPoint_Krasovsky();
~PrjPoint_Krasovsky();
};
class   PrjPoint_IUGG1975   :   virtual   public   PrjPoint
{
public:
PrjPoint_IUGG1975();
~PrjPoint_IUGG1975();
};
double   Dms2Rad(double   Dms);
double   Rad2Dms(double   Rad);
#endif   /*   ndef   _COORTRANS_H_INCLUDED   */

上一篇:最短路径A*算法实现(Javascript)

下一篇:地图着色算法C语言实现