O2/30+(4*K2+5)*M2-L2)*O2/12+1)*N2*I2*O2/2
計算結果X

T2
=((((L2-18)*L2-(58*L2-14)*K2+5)*O2/20+M2-L2)*O2/6+1)*N2*(H2*J2)
計算結果Y




表中公式的來源及EXCEL軟件的操作方法，請參閱有關資料，這里不再贅述。按上面表格中的公式輸入到相應單元格后，就可方便地由經緯度求得平面直角坐標。當輸入完所有的經緯度后，用鼠標下拉即可得到所有的計算結果。表中的許多單元格公式為中間過程，可以用EXCEL的列隱藏功能把這些沒有必要顯示的列隱藏起來，表面上形成標準的計算報表，使整個計算表簡單明了。從理論上講，可計算的數據量是無限的，當第一次輸入公式后，相當于自己完成了一軟件的編制，可另存起來供今后重復使用，一勞永逸。



二、GPS坐標轉換方法與面積計算

GPS所采用的坐標系是美國國防部1984世界坐標系，簡稱WGS-84，它是一個協議地球參考系，坐標系原點在地球質心。GPS的測量結果與我國的54 系或80系坐標相差幾十米至一百多米，隨區域不同，差別也不同，經粗落統計，我國西部相差70米左右，東北部140米左右，南部75米左右，中部45米左右。由此可見，必須將WGS-84坐標進行坐標系轉換才能供標圖使用。坐標系之間的轉換一般采用七參數法或三參數法，其中七參數為X平移、Y平移、Z平移、X旋轉、Y旋轉、Z旋轉以及尺度比參數，若忽略旋轉參數和尺度比參數則為三參數方法，三參數法為七參數法的特例。這里的Z、Y、Z是空間大地直角坐標系坐標，為轉換過程的中間值。在實際工作中我們常用的是平面直角坐標，是否可以跳過空間直角坐標系，省略復雜的運算，進行簡單轉換呢？為此，筆者進行了長期的實踐，證明是可行的。其在原理是：不把GPS所測定的WGS-84坐標當作WGS-84坐標，而是當作具有一定系統性誤差的54系坐標值，然后通過國家已知點糾正，消除該系統誤差。我們暫把該方法稱作坐標改正法，下面以WGS-84坐標轉換成54系坐標為例，介紹數據處理方法：

首先，在測區附近選擇一國家已知點，在該已知點上用GPS測定WGPS-84坐標系經緯度B和L，把此坐標視為有誤差的54系坐標，利用54系EXCEL將經緯度BL轉換成平面直角坐標X’Y’，然后與已知坐標比較則可計算出偏移量：

△X=X－X’

△Y=Y－Y’

式中的X、Y為國家控制點的已知坐標，X’、Y’為測定坐標，△X和△Y為偏移量。

求得偏移量后，就可以用此偏移量糾正測區內的其他測量點了。把其他GPS測量點的經緯度測量值，轉換成平面坐標X’Y’，在此XY坐標值上直接加上偏移值就得到了轉換后的54系坐標：

X=X’+△X

Y=Y’+△Y

在上述EXCEL計算表的最后兩列，附加上求得的改正數并分別與計算出來的XY相加后，即得到轉換結果。若測量路線是一閉合區域的話，可把計算結果按路線順序排列起來，再輸入相應的計算公式，即可計算出該區域的面積。有關用坐標計算面積的原理與公式，這里不再敘述，讀者可參閱有關資料。需要說明的是，面積的計算精度基本上不受坐標轉換精度的影響，若只需要求算面積的話，可不進行坐標系轉換這一步，只需要把BL化成XY就行了。

就1：1萬比例尺成圖而言，在一般的縣行政區范圍內（如40Km×40Km），用此簡單的坐標改正法進行轉換與較復雜的七參數法沒有多大差別。能否滿足 1：1萬比例尺變更調查的要求，主要取決于GPS接收機本身的精度，與轉換方法的選擇關系不大。當面積較大時，使用該方法可能會使誤差增大，這時可考慮分區域轉換。



































#include <windows.h>
#include  <math.h>
#include <stdio.h>
#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()
{
}




























 轉載： 高斯投影變換_cpp
#include 
#include 
#include 
#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()
{
}



以下是11位的：＄ＧＰＲＭＣ