?? source.txt
字號:
這段代碼是偶找來的，里面存在一些錯誤的代碼偶沒有改正，哈哈，希望大家自己去找出來，因為代碼不是偶自己寫的，如何調用也不要問偶，如果需要的話偶可以直接提供dll（偶不想寫注釋），保證絕對正確，計算出來的數值和國家正在使用的軟件系統中的轉換保持一致。偶不是搞地理信息系統的，其他問題不要找我。呵呵。
代碼如下：
//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; // 中央子午線經度
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 */
//高斯投影正、反算
//////6度帶寬 54年北京坐標系
//高斯投影由經緯度(Unit:DD)反算大地坐標(含帶號，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年北京坐標系參數
////a=6378140.0; f=1/298.257; //80年西安坐標系參數
ProjNo = (int)(longitude / ZoneWide) ;
longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
longitude0 = longitude0 * iPI ;
latitude0=0;
longitude1 = longitude * iPI ; //經度轉換為弧度
latitude1 = latitude * iPI ; //緯度轉換為弧度
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;
}
//高斯投影由大地坐標(Unit:Metres)反算經緯度(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年北京坐標系參數
////a=6378140.0; f=1/298.257; //80年西安坐標系參數
ZoneWide = 6; ////6度帶寬
ProjNo = (int)(X/1000000L) ; //查找帶號
longitude0 = (ProjNo-1) * ZoneWide + ZoneWide / 2;
longitude0 = longitude0 * iPI ; //中央經線
X0 = ProjNo*1000000L+500000L;
Y0 = 0;
xval = X-X0; yval = Y-Y0; //帶內大地坐標
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;
//計算經度(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);
//轉換為度 DD
*longitude = longitude1 / iPI;
*latitude = latitude1 / iPI;
}
對利用EXCEL電子表格進行高斯投影換算的方法進行了較詳細的介紹，對如何進行GPS坐標系轉換進行了分析，提出了一種簡單實用的坐標改正轉換方法，介紹了用EXCEL完成轉換的思路。
[關鍵字] 電子表格；GPS；坐標轉換
作為尖端技術GPS，能方便快捷性地測定出點位坐標，無論是操作上還是精度上，比全站儀等其他常規測量設備有明顯的優越性。隨著我國各地GPS差分臺站的不斷建立以及美國SA政策的取消，使得單機定位的精度大大提高，有的已經達到了亞米級精度，能夠滿足國土資源調查、土地利用更新、遙感監測、海域使用權清查等工作的應用。在一般情況下，我們使用的是1954年北京坐標系或1980年西安坐標系（以下分別簡稱54系和80系），而GPS測定的坐標是WGS- 84坐標系坐標，需要進行坐標系轉換。對于非測量專業的工作人員來說，雖然GPS定位操作非常容易，但坐標轉換則難以掌握，EXCEL是比較普及的電子表格軟件，能夠處理較復雜的數學運算，用它來進行GPS坐標轉換、面積計算會非常輕松自如。要進行坐標系轉換，離不開高斯投影換算，下面分別介紹用 EXCEL進行換算的方法和GPS坐標轉換方法。
一、用EXCEL進行高斯投影換算
從經緯度BL換算到高斯平面直角坐標XY（高斯投影正算），或從XY換算成BL（高斯投影反算），一般需要專用計算機軟件完成，在目前流行的換算軟件中，存在一個共同的不足之處，就是靈活性較差，大都需要一個點一個點地進行，不能成批量地完成，給實際工作帶來許多不便。筆者發現，用EXCEL可以很直觀、方便地完成坐標換算工作，不需要編制任何軟件，只需要在EXCEL的相應單元格中輸入相應的公式即可。下面以54系為例，介紹具體的計算方法。
完成經緯度BL到平面直角坐標XY的換算，在EXCEL中大約需要占用21列，當然讀者可以通過簡化計算公式或考慮直觀性，適當增加或減少所占列數。在 EXCEL中，輸入公式的起始單元格不同，則反映出來的公式不同，以公式從第2行第1列（A2格）為起始單元格為例，各單元格的公式如下：
單元格
單元格內容
說明
A2
輸入中央子午線，以度.分秒形式輸入，如115度30分則輸入115.30
起算數據L0
B2
=INT(A2)+(INT(A2*100)-INT(A2)*100)/60+(A2*10000-INT(A2*100)*100)/3600
把L0化成度
C2
以度小數形式輸入緯度值，如38°14′20″則輸入38.1420
起算數據B
D2
以度小數形式輸入經度值
起算數據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)*
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -