橢球體上某區域面積的求算，及該區域蘭伯特投影與墨卡託投影到二維平面後面積對比

    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 ；先後兩個點的經度對應 λ一、λ2，再利用積分的原理就能夠計算獲得橢球體上的某區域的面積啦！

=================================下課了================================

    假設咱們拿到的數據是某區域墨卡託投影后的二維平面上一系列的區域邊界點數據，那麼我接下來的步驟將分爲六步：一、計算二維平面上的墨卡託投影后的平面面積 S1  -->  二、墨卡託投影反算，獲得每一個點在地球橢球體上對應的經緯度座標  -->  三、計算地球橢球體上該區域的面積 S2 --> 四、把該區域進行蘭伯特投影獲得二維平面上又一系列的區域邊界點數據  -->  五、計算二維平面上的蘭伯特投影后的平面面積 S3 -->  六、對比三種面積的區別。

     

一、計算二維平面上的墨卡託投影后的平面面積 S1

程序以下：

 

 

二、墨卡託投影反算，獲得每一個點在地球橢球體上對應的經緯度座標

    先利用墨卡託投影反解公式，計算B、L，其中對於求解 L 須要用到 K ，K 的值由公式求出；對於求解 B，則須要設定一個初始值，而後進行迭代求解，直到先後兩次計算出的 B 之差小於0.00000000001，則認爲後一次計算的 B 的值爲最終解。

程序以下：

 

 

在ArcGIS中查看墨卡託反算先後的數據，對比顯示以下：

                 反算前                                    反算後

 

三、計算地球橢球體上該區域的面積 S2

    由於ds足夠小，因此把梯形近似看作一個矩形來計算，矩形的長爲東西方向的弧長，寬爲南北方向的弧長。根據弧長計算公式：弧長=半徑*弧度，涉及到子午圈曲率半徑M和主法截面曲率半徑N的計算公式。經過查閱資料，可知：

南北方向上的弧長d=M*d；東西方向上的弧長d=N*d

    對於每個小梯形，用先後兩個點的平均緯度做爲 φ2，以該區域的最低緯度做爲 φ1，先後兩個點的經度對應 λ一、λ2，再利用積分的原理來計算獲得橢球體上的該區域面積。

程序以下：

先聲明和賦值程序中將會用到的基本數據長半軸a、第一偏愛率e、基準緯度（江蘇省最低緯度）B；並經過指向矢量文件的指針得到先後兩點的緯度B一、B2和經度L一、L2。用TB來代替2-1，用AB來代替（1+2）/2：

計算小梯形的面積積分公式所須要用到的參數K、A、B、C、D，帶入公式進行計算面積：

 

四、把該區域進行蘭伯特投影獲得二維平面上又一系列的區域邊界點數據

這裏參考《.gen地圖文件的投影編程實現（以墨卡託投影和蘭伯特投影爲例）》http://www.cnblogs.com/to-sunshine/p/6048438.html。

 

五、計算二維平面上的蘭伯特投影后的平面面積 S3

方法同第一步，只是帶入的數據不一樣。

 

六、對比三種面積的區別

整個程序主函數以下：

運行後結果以下：

可見，進行Albers等積投影以後，矢量面的面積變化偏差相比之總體面積來講較小，因此視爲等積投影是成功的。