hello，最近yogurt給大家的更新很頻繁哦~~今天要分享的內容是緊接着前面兩篇的內容做的擴展~~
我們不僅要求取某地區在地球橢球體這個三維空間中的面積，還要與該地區投影到二維空間后平面多邊形的面積進行對比。怎么求取二維平面多邊形的面積，大家可以看看我之前寫過的《求解多邊形面積2S= Σ【Xi （Yi+1-Yi-1）】,（i屬於1~n），公式解析及編程實現》http://www.cnblogs.com/to-sunshine/p/7642222.html。至於怎么進行蘭伯特投影墨卡托投影，大家可以參考我之前寫過的《.gen地圖文件的投影編程實現（以墨卡托投影和蘭伯特投影為例）》http://www.cnblogs.com/to-sunshine/p/6048438.html。
============================yogurt小課堂開課了===========================
下面yogurt要跟大家分享這次需要用到的知識點：
地球橢球體表面上的梯形面積：
取長度很小很小的一個長和一個寬，組成的梯形是非常小的，近似可以看作矩形。那么 S = 長 X 寬 ，如下圖：
（圖來自我老師的PPT）
我們知道 AB = CD ，是南北方向上的邊長，假設與任意一條經線平行。由上節的知識，我們不難知道： AB = CD = M X dφ （dφ是AD所在緯線與BC所在緯線的緯度間隔，很小）；
同樣，由於緯度間隔非常小，可以近似看作 BC ≈ AD ，是東西方向上的邊長，假設與任意一條緯線平行。由上節的知識，我們也不難知道： BC ≈ AD = r X dλ （dλ是AB所在經線與CD所在經線的經度間隔，很小），而 r = N X cosφ（φ看作AD所在緯線或者BC所在緯線的緯度，由於間隔很小，所以只要所有梯形統一用上邊的或者統一都用下邊的緯度即可）。因此，BC = AD = N X cosφ X dλ 。
綜上，便能夠得到地球橢球體上梯形的面積 dS = （M X dφ） X （N X cosφ X dλ）= M N cosφ dλ dφ ，那么 S = 
。
對於每一個小梯形，用前后兩個點的平均緯度作為 φ2，以該區域的最低緯度作為 φ1，dφ = φ2 - φ1 ；前后兩個點的經度對應 λ1、λ2，再利用積分的原理就可以計算得到橢球體上的某區域的面積啦！
=================================下課了================================
假設我們拿到的數據是某區域墨卡托投影后的二維平面上一系列的區域邊界點數據，那么我接下來的步驟將分為六步：1、計算二維平面上的墨卡托投影后的平面面積 S1 --> 2、墨卡托投影反算，得到每個點在地球橢球體上對應的經緯度坐標 --> 3、計算地球橢球體上該區域的面積 S2 --> 4、把該區域進行蘭伯特投影得到二維平面上又一系列的區域邊界點數據 --> 5、計算二維平面上的蘭伯特投影后的平面面積 S3 --> 6、對比三種面積的區別。
1、計算二維平面上的墨卡托投影后的平面面積 S1
程序如下：
2、墨卡托投影反算，得到每個點在地球橢球體上對應的經緯度坐標
先利用墨卡托投影反解公式，計算B、L，其中對於求解 L 需要用到 K ，K 的值由公式求出；對於求解 B，則需要設定一個初始值，然后進行迭代求解，直到前后兩次計算出的 B 之差小於0.00000000001，則認為后一次計算的 B 的值為最終解。
程序如下：
在ArcGIS中查看墨卡托反算前后的數據，對比顯示如下：
反算前 反算后
3、計算地球橢球體上該區域的面積 S2
因為ds足夠小，所以把梯形近似看做一個矩形來計算，矩形的長為東西方向的弧長，寬為南北方向的弧長。根據弧長計算公式：弧長=半徑*弧度，涉及到子午圈曲率半徑M和主法截面曲率半徑N的計算公式。通過查閱資料，可知：
南北方向上的弧長d=M*d；東西方向上的弧長d=N*d
對於每一個小梯形，用前后兩個點的平均緯度作為 φ2，以該區域的最低緯度作為 φ1，前后兩個點的經度對應 λ1、λ2，再利用積分的原理來計算得到橢球體上的該區域面積。
程序如下：
先聲明和賦值程序中將會用到的基本數據長半軸a、第一偏心率e、基准緯度（江蘇省最低緯度）B；並通過指向矢量文件的指針獲得前后兩點的緯度B1、B2和經度L1、L2。用TB來代替2-1，用AB來代替（1+2）/2：
計算小梯形的面積積分公式所需要用到的參數K、A、B、C、D，帶入公式進行計算面積：
4、把該區域進行蘭伯特投影得到二維平面上又一系列的區域邊界點數據
這里參考《.gen地圖文件的投影編程實現（以墨卡托投影和蘭伯特投影為例）》http://www.cnblogs.com/to-sunshine/p/6048438.html。
5、計算二維平面上的蘭伯特投影后的平面面積 S3
方法同第一步，只是帶入的數據不同。
6、對比三種面積的區別
整個程序主函數如下：
運行后結果如下：
可見，進行Albers等積投影之后，矢量面的面積變化誤差相比之整體面積來說較小，所以視為等積投影是成功的。