這段代碼是偶找來(lái)的，里面存在一些錯(cuò)誤的代碼偶沒(méi)有改正，哈哈，希望大家自己去找出來(lái)，因?yàn)榇a不是偶自己寫(xiě)的，如何調(diào)用也不要問(wèn)偶，如果需要的話(huà)偶可以直接提供dll（偶不想寫(xiě)注釋?zhuān)ＷC絕對(duì)正確，計(jì)算出來(lái)的數(shù)值和國(guó)家正在使用的軟件系統(tǒng)中的轉(zhuǎn)換保持一致。偶不是搞地理信息系統(tǒng)的，其他問(wèn)題不要找我。呵呵。  
  代碼如下：  
  //CPP文件  
   
  #include   "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  
  #include    
  const   double   PI   =   3.14159265353846;  
  class   PrjPoint  
  {  
  public:  
  double   L0;   //   中央子午線經(jīng)度  
  double   B,   L;   //   大地坐標(biāo)  
  double   x,   y;   //   高斯投影平面坐標(biāo)  
  public:  
  BOOL   BL2xy();  
  BOOL   xy2BL();  
  protected:  
  double   a,   f,   e2,   e12;   //   基本橢球參數(shù)  
  double   A1,   A2,   A3,   A4;   //   用于計(jì)算X的橢球參數(shù)  
  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   */ 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  //高斯投影正、反算
//////6度帶寬 54年北京坐標(biāo)系
//高斯投影由經(jīng)緯度(Unit:DD)反算大地坐標(biāo)(含帶號(hào)，Unit:Metres)
void GaussProjCal(double longitude, double latitude, double *X, double *Y)
{
int ProjNo=0; int ZoneWide; ////帶寬
double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;
double a,f, e2,ee, NN, T,C,A, M, iPI;
iPI = 0.0174532925199433; ////3.1415926535898/180.0;
ZoneWide = 6; ////6度帶寬
a=6378245.0; f=1.0/298.3; //54年北京坐標(biāo)系參數(shù)
////a=6378140.0; f=1/298.257; //80年西安坐標(biāo)系參數(shù)
ProjNo = (int)(longitude / ZoneWide) ;
longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
longitude0 = longitude0 * iPI ;
latitude0=0;
longitude1 = longitude * iPI ; //經(jīng)度轉(zhuǎn)換為弧度
latitude1 = latitude * iPI ; //緯度轉(zhuǎn)換為弧度
e2=2*f-f*f;
ee=e2*(1.0-e2);
NN=a/sqrt(1.0-e2*sin(latitude1)*sin(latitude1));
T=tan(latitude1)*tan(latitude1);
C=ee*cos(latitude1)*cos(latitude1);
A=(longitude1-longitude0)*cos(latitude1);
M=a*((1-e2/4-3*e2*e2/64-5*e2*e2*e2/256)*latitude1-(3*e2/8+3*e2*e2/32+45*e2*e2
*e2/1024)*sin(2*latitude1)
+(15*e2*e2/256+45*e2*e2*e2/1024)*sin(4*latitude1)-(35*e2*e2*e2/3072)*sin(6*l
atitude1));
xval = NN*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*ee)*A*A*A*A*A/120);
yval = M+NN*tan(latitude1)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
+(61-58*T+T*T+600*C-330*ee)*A*A*A*A*A*A/720);
X0 = 1000000L*(ProjNo+1)+500000L;
Y0 = 0;
xval = xval+X0; yval = yval+Y0;
*X = xval;
*Y = yval;
}


//高斯投影由大地坐標(biāo)(Unit:Metres)反算經(jīng)緯度(Unit:DD)
void GaussProjInvCal(double X, double Y, double *longitude, double *latitude)

{
int ProjNo; int ZoneWide; ////帶寬
double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;
double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;
iPI = 0.0174532925199433; ////3.1415926535898/180.0;
a = 6378245.0; f = 1.0/298.3; //54年北京坐標(biāo)系參數(shù)
////a=6378140.0; f=1/298.257; //80年西安坐標(biāo)系參數(shù)
ZoneWide = 6; ////6度帶寬
ProjNo = (int)(X/1000000L) ; //查找?guī)?hào)
longitude0 = (ProjNo-1) * ZoneWide + ZoneWide / 2;
longitude0 = longitude0 * iPI ; //中央經(jīng)線
X0 = ProjNo*1000000L+500000L;
Y0 = 0;
xval = X-X0; yval = Y-Y0; //帶內(nèi)大地坐標(biāo)
e2 = 2*f-f*f;
e1 = (1.0-sqrt(1-e2))/(1.0+sqrt(1-e2));
ee = e2/(1-e2);
M = yval;
u = M/(a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256));
fai = u+(3*e1/2-27*e1*e1*e1/32)*sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(
4*u)
+(151*e1*e1*e1/96)*sin(6*u)+(1097*e1*e1*e1*e1/512)*sin(8*u);
C = ee*cos(fai)*cos(fai);
T = tan(fai)*tan(fai);
NN = a/sqrt(1.0-e2*sin(fai)*sin(fai));
R = a*(1-e2)/sqrt((1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai))*(1-e2*sin
(fai)*sin(fai)));
D = xval/NN;
//計(jì)算經(jīng)度(Longitude) 緯度(Latitude)
longitude1 = longitude0+(D-(1+2*T+C)*D*D*D/6+(5-2*C+28*T-3*C*C+8*ee+24*T*T)*D
*D*D*D*D/120)/cos(fai);
latitude1 = fai -(NN*tan(fai)/R)*(D*D/2-(5+3*T+10*C-4*C*C-9*ee)*D*D*D*D/24
+(61+90*T+298*C+45*T*T-256*ee-3*C*C)*D*D*D*D*D*D/720);
//轉(zhuǎn)換為度 DD
*longitude = longitude1 / iPI;
*latitude = latitude1 / iPI;
}




































對(duì)利用EXCEL電子表格進(jìn)行高斯投影換算的方法進(jìn)行了較詳細(xì)的介紹，對(duì)如何進(jìn)行GPS坐標(biāo)系轉(zhuǎn)換進(jìn)行了分析，提出了一種簡(jiǎn)單實(shí)用的坐標(biāo)改正轉(zhuǎn)換方法，介紹了用EXCEL完成轉(zhuǎn)換的思路。

[關(guān)鍵字] 電子表格；GPS；坐標(biāo)轉(zhuǎn)換

作為尖端技術(shù)GPS，能方便快捷性地測(cè)定出點(diǎn)位坐標(biāo)，無(wú)論是操作上還是精度上，比全站儀等其他常規(guī)測(cè)量設(shè)備有明顯的優(yōu)越性。隨著我國(guó)各地GPS差分臺(tái)站的不斷建立以及美國(guó)SA政策的取消，使得單機(jī)定位的精度大大提高，有的已經(jīng)達(dá)到了亞米級(jí)精度，能夠滿(mǎn)足國(guó)土資源調(diào)查、土地利用更新、遙感監(jiān)測(cè)、海域使用權(quán)清查等工作的應(yīng)用。在一般情況下，我們使用的是1954年北京坐標(biāo)系或1980年西安坐標(biāo)系（以下分別簡(jiǎn)稱(chēng)54系和80系），而GPS測(cè)定的坐標(biāo)是WGS- 84坐標(biāo)系坐標(biāo)，需要進(jìn)行坐標(biāo)系轉(zhuǎn)換。對(duì)于非測(cè)量專(zhuān)業(yè)的工作人員來(lái)說(shuō)，雖然GPS定位操作非常容易，但坐標(biāo)轉(zhuǎn)換則難以掌握，EXCEL是比較普及的電子表格軟件，能夠處理較復(fù)雜的數(shù)學(xué)運(yùn)算，用它來(lái)進(jìn)行GPS坐標(biāo)轉(zhuǎn)換、面積計(jì)算會(huì)非常輕松自如。要進(jìn)行坐標(biāo)系轉(zhuǎn)換，離不開(kāi)高斯投影換算，下面分別介紹用 EXCEL進(jìn)行換算的方法和GPS坐標(biāo)轉(zhuǎn)換方法。

一、用EXCEL進(jìn)行高斯投影換算

從經(jīng)緯度BL換算到高斯平面直角坐標(biāo)XY（高斯投影正算），或從XY換算成BL（高斯投影反算），一般需要專(zhuān)用計(jì)算機(jī)軟件完成，在目前流行的換算軟件中，存在一個(gè)共同的不足之處，就是靈活性較差，大都需要一個(gè)點(diǎn)一個(gè)點(diǎn)地進(jìn)行，不能成批量地完成，給實(shí)際工作帶來(lái)許多不便。筆者發(fā)現(xiàn)，用EXCEL可以很直觀、方便地完成坐標(biāo)換算工作，不需要編制任何軟件，只需要在EXCEL的相應(yīng)單元格中輸入相應(yīng)的公式即可。下面以54系為例，介紹具體的計(jì)算方法。

完成經(jīng)緯度BL到平面直角坐標(biāo)XY的換算，在EXCEL中大約需要占用21列，當(dāng)然讀者可以通過(guò)簡(jiǎn)化計(jì)算公式或考慮直觀性，適當(dāng)增加或減少所占列數(shù)。在 EXCEL中，輸入公式的起始單元格不同，則反映出來(lái)的公式不同，以公式從第2行第1列（A2格）為起始單元格為例，各單元格的公式如下：

單元格
單元格內(nèi)容
說(shuō)明

A2
輸入中央子午線，以度.分秒形式輸入，如115度30分則輸入115.30
起算數(shù)據(jù)L0

B2
=INT(A2)+(INT(A2*100)-INT(A2)*100)/60+(A2*10000-INT(A2*100)*100)/3600
把L0化成度

C2
以度小數(shù)形式輸入緯度值，如38°14′20″則輸入38.1420
起算數(shù)據(jù)B

D2
以度小數(shù)形式輸入經(jīng)度值
起算數(shù)據(jù)L

E2
=INT(C2)+(INT(C2*100)-INT(C2)*100)/60+(C2*10000-INT(C2*100)*100)/3600
把B化成度

F2
=INT(D2)+(INT(D2*100)-INT(D2)*100)/60+(D2*10000-INT(D2*100)*100)/3600
把L化成度

G2
=F2-B2
L-L0

H2
=G2/57.2957795130823
化作弧度

I2
=TAN(RADIANS(E2))
Tan(B)

J2
=COS(RADIANS(E2))
COS(B)

K2
=0.006738525415*J2*J2


L2
=I2*I2


M2
=1+K2


N2
=6399698.9018/SQRT(M2)


O2
=H2*H2*J2*J2


P2
=I2*J2


Q2
=P2*P2


R2
=(32005.78006+Q2*(133.92133+Q2*0.7031))


S2
=6367558.49686*E2/57.29577951308-P2*J2*R2+((((L2-58)*L2+61)*