一：首先補充下理論知識（網上找的）

對於不同的橢球體之間的轉換要用到七引數的方法，這個不還沒有研究。本文不是講的這個。

因為本文章主講演算法，理論知識就不多說了。

分帶方法
1．我國採用6度分帶和3度分帶： 　　
         1∶2.5萬及1∶5萬的地形圖採用6度分帶投影，即經差為6度，從零度子午線開始，自西向東每個經差6度為一投影帶，全球共分60個帶，用1，2， 3，4，5，……表示．即東經0～6度為第一帶，其中央經線的經度為東經3度，東經6～12度為第二帶，其中央經線的經度為9度。 　　
         1∶1萬的地形圖採用3度分帶，從東經1.5度的經線開始，每隔3度為一帶，用1，2，3，……表示，全球共劃分120個投影帶，即東經1.5～ 4.5度為第1帶，其中央經線的經度為東經3度，東經4.5～7.5度為第2帶，其中央經線的經度為東經6度．比如福建省位於東經113度－東經120度之間， 跨第38、39、40共計3個帶，其中東經115.5度以西為第38帶，其中央經線為東經114度；東經115.5～118.5度為39帶，其中央經線為 東經117度；東經118.5度以東到山海關為40帶，其中央經線為東經120度。 地形圖上公里網橫座標前2位就是帶號，例如：1∶5萬地形圖上的橫座標為20345486，其中20即為帶號，345486為橫座標值。
2．當地中央經線經度的計算
           六度帶中央經線經度的計算：當地中央經線經度＝６°×當地帶號－３°，例如：地形圖上的橫座標為２０３４５，其所處的六度帶的中央經線經度為：６°×２０－３°＝１１７°（適用於１∶２．５萬和１∶５萬地形圖）；６°×49－180－３°＝111°（適用於１∶50萬以下地形圖）。
三度帶中央經線經度的計算：中央經線經度＝３°×當地帶號（適用於１∶１萬地形圖）。

我國x座標都是正的，y座標的最大值（在赤道上）約為330km。為了避免出現負的橫座標，可在橫座標上加上500 OOOm。此外還應在座標前面再冠以帶號。這種座標稱為國家統一座標。

二. 轉換演算法，用C++來實現的

 int zonewide=6;     //3度帶比較準，可以用6度帶
 int base_type=WGS84;   //投影基準面型別：北京54基準面為54，西安80基準面為80，WGS84基準面
 double L,B;//大地經度，緯度

 //----------初始化已知引數-----------
 L = longitude;//經度
 B = latitude;//緯度

 double l,L0;//投影帶中央線起算的經差，投影帶中央子午線
 double X,Y;//直角座標
 double f,e2,a,b,e12;////參考橢球體扁率,第一偏心率，長半軸,短半軸，第二偏心率
 double t;//儲存tan B
 double m;//儲存lcosB
 double q2;//卯酉圈的分量 的平方
 double M,N;//子午，卯酉圈曲半徑
 double S;//子午線弧長
 double A0,B0,C0,D0,E0;
 double n;//投影帶號
 double a_1_e2;//a*(1-e*e)的值

 if(WGS84 == base_type)
 {
  a = 6378137;
  b = 6356752.3142;
  f = 1/298.257223563;
  e2 = 0.0066943799013;//--第一偏心率
  e12 = 0.00673949674227;//--第二偏心率
 }
 else if(BJ54 == base_type)
 {
  a = 6378245;
  b = 6356863.0187730473;
  f = 1/298.3;
  e2 = 0.006693421622966;//--第一偏心率
  e12 = 0.006738525414683;//--第二偏心率
 }
 else if(XIAN80 == base_type)
 {
  a = 6378140;
  b = 6356755.2881575287;
  f = 1/298.257;
  e2 = 0.006694384999588;//--第一偏心率
  e12 = 0.006739501819473;//--第二偏心率
 }

 //求中央經線
 if (zonewide==6)
 {
  n=(int)(L/6)+1;
  L0=n*zonewide-3;
 }
 else
 {
  n=(int)((L-1.5)/3)+1;
  L0=n*3;
 }

 //---基本計算公式 a 長半軸，b 短半軸--------
 //f = (a-b)/a;//--橢球體扁率
 //e2 = (a*a-b*b)/(a*a);//--第一偏心率
 //e12 = (a*a-b*b)/(b*b);//--第二偏心率
 a_1_e2 = a*(1-e2);

 //---轉化為弧度,sin,cos 等運算都是以弧度為基礎的
 L0=L0*PI/180;
 L=L*PI/180;
 B=B*PI/180;
 l = L - L0;
 t = tan(B);
 q2 = e12*(cos(B)*cos(B));
 N = a/sqrt(1-e2*sin(B)*sin(B));////卯酉圈的曲率半徑
 m = l*cos(B);

 //--以下演算法是對原始公式的多項式進行了合併處理，儲存計算更準
 double m0 = a_1_e2;
 double m2 = 3.0 * e2 * m0 / 2.0;
 double m4 = 5.0 * e2 * m2 / 4.0;
 double m6 = 7.0 * e2 * m4 / 6.0;
 double m8 = 9.0 * e2 * m6 / 8.0;
 double a0 = m0 + m2 / 2.0 + 3.0 * m4 / 8.0 + 5.0 * m6 / 16.0 + 35.0 * m8 / 128;
 double a2 = m2 / 2.0 + m4 / 2.0 + 15.0 * m6 / 32.0 + 7.0 * m8 / 16.0;
 double a4 = m4 / 8.0 + 3.0 * m6 / 16.0 + 7.0 * m8 / 32.0;
 double a6 = m6 / 32.0 + m8 / 16.0;
 double a8 = m8 / 128.0;
 A0 = a0;
 B0 = a2;
 C0 = a4;
 D0 = a6;
 E0 = a8;
 S = (A0 * B - B0* sin(2.0*B)/2.0 + C0 * sin(4.0*B)/4.0 - D0 *sin(6.0*B)/6.0 + E0*sin(8.0*B)/8.0);

 //----高斯投影公式-------
 X = S + N * t *((0.5 + ((5.0 - t * t + 9 * q2 +4 * q2 * q2)/24 + (61.0 - 58.0 * t * t + pow(t,4)) * m * m /720.0) * m * m) * m * m);//pow(t,4)為t的4次方
 Y = N * ((1 + ((1 - t * t +q2)/6.0 +(5.0 - 18.0 * t * t + pow(t,4) + 14 * q2 - 58 * q2 * t * t ) * m * m / 120.0) * m * m ) * m);
 //----根據中國國情座標縱軸西移500公里當作起始軸
 Y = 1000000 * n + 500000 + Y;

三. 面積計算公式

注意本公式的首末結點都為（0,0），求之前要把所有的結點都減去第一個結點。

PS：這樣求出來的面積是有負值的，因為這是向量的，和你點的順序有關係。具體來說順時針則為負，逆時針則為正。這樣又可以多了一種判斷地物向量化方向的方法。

四. 結果檢測

這是ArcMap中的面積計算結果，我把各個座標輸入進去後得到的結果是654865.33362435829。

誤差大約在0.81%。

原因分析：

1. 在ArcMap中我選的投影為UTM投影（0.9996），並不是嚴格的高斯克裡格投影。

2.座標輸入的時候，並不是用專門的軟體讀取出來的，而是手動點選然後copy的。

3.由於簡單的結算採用的並不是七引數的轉換，所以...

 不過怎麼樣，證明我的轉換演算法和麵積計算公式都是可行的，哈哈。。。

改天把工程和資料打包上傳。