??? hello，最近yogurt給大家的更新很頻繁哦~~今天要分享的內(nèi)容是緊接著前面兩篇的內(nèi)容做的擴展~~

??? 我們不僅要求取某地區(qū)在地球橢球體這個三維空間中的面積，還要與該地區(qū)投影到二維空間后平面多邊形的面積進行對比。怎么求取二維平面多邊形的面積，大家可以看看我之前寫過的《求解多邊形面積2S= Σ【Xi （Yi+1-Yi-1）】,（i屬于1~n），公式解析及編程實現(xiàn)》http://www.cnblogs.com/to-sunshine/p/7642222.html。至于怎么進行蘭伯特投影墨卡托投影，大家可以參考我之前寫過的《.gen地圖文件的投影編程實現(xiàn)（以墨卡托投影和蘭伯特投影為例）》http://www.cnblogs.com/to-sunshine/p/6048438.html。

============================yogurt小課堂開課了===========================

??? 下面yogurt要跟大家分享這次需要用到的知識點：

??? 地球橢球體表面上的梯形面積：

??? 取長度很小很小的一個長和一個寬，組成的梯形是非常小的，近似可以看作矩形。那么 S = 長 X 寬 ，如下圖：

（圖來自我老師的PPT）

??? ??? 我們知道 AB = CD ，是南北方向上的邊長，假設(shè)與任意一條經(jīng)線平行。由上節(jié)的知識，我們不難知道： AB = CD = M X dφ （dφ是AD所在緯線與BC所在緯線的緯度間隔，很小）；

??? 同樣，由于緯度間隔非常小，可以近似看作 BC ≈ AD ，是東西方向上的邊長，假設(shè)與任意一條緯線平行。由上節(jié)的知識，我們也不難知道： BC ≈ AD = r X dλ （dλ是AB所在經(jīng)線與CD所在經(jīng)線的經(jīng)度間隔，很小），而 r = N X cosφ（φ看作AD所在緯線或者BC所在緯線的緯度，由于間隔很小，所以只要所有梯形統(tǒng)一用上邊的或者統(tǒng)一都用下邊的緯度即可）。因此，BC = AD = N X cosφ X dλ 。

??? 綜上，便能夠得到地球橢球體上梯形的面積 dS = （M X dφ） X （N X cosφ X dλ）= M N cosφ dλ dφ ，那么 S =?。

??? 對于每一個小梯形，用前后兩個點的平均緯度作為 φ2，以該區(qū)域的最低緯度作為 φ1，dφ = φ2 - φ1 ；前后兩個點的經(jīng)度對應(yīng) λ1、λ2，再利用積分的原理就可以計算得到橢球體上的某區(qū)域的面積啦！

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

??? 假設(shè)我們拿到的數(shù)據(jù)是某區(qū)域墨卡托投影后的二維平面上一系列的區(qū)域邊界點數(shù)據(jù)，那么我接下來的步驟將分為六步：1、計算二維平面上的墨卡托投影后的平面面積 S1? -->? 2、墨卡托投影反算，得到每個點在地球橢球體上對應(yīng)的經(jīng)緯度坐標? -->? 3、計算地球橢球體上該區(qū)域的面積 S2 --> 4、把該區(qū)域進行蘭伯特投影得到二維平面上又一系列的區(qū)域邊界點數(shù)據(jù)? -->? 5、計算二維平面上的蘭伯特投影后的平面面積 S3 -->? 6、對比三種面積的區(qū)別。

?????

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

程序如下：

?

?

2、墨卡托投影反算，得到每個點在地球橢球體上對應(yīng)的經(jīng)緯度坐標

??? 先利用墨卡托投影反解公式，計算B、L，其中對于求解 L 需要用到 K ，K 的值由公式求出；對于求解 B，則需要設(shè)定一個初始值，然后進行迭代求解，直到前后兩次計算出的 B 之差小于0.00000000001，則認為后一次計算的 B 的值為最終解。

程序如下：

?

?

在ArcGIS中查看墨卡托反算前后的數(shù)據(jù)，對比顯示如下：

????????? ? ???? 反算前?????????????????????????? ? ? ? ?? 反算后

?

3、計算地球橢球體上該區(qū)域的面積 S2

??? 因為ds足夠小，所以把梯形近似看做一個矩形來計算，矩形的長為東西方向的弧長，寬為南北方向的弧長。根據(jù)弧長計算公式：弧長=半徑*弧度，涉及到子午圈曲率半徑M和主法截面曲率半徑N的計算公式。通過查閱資料，可知：

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

??? 對于每一個小梯形，用前后兩個點的平均緯度作為 φ2，以該區(qū)域的最低緯度作為 φ1，前后兩個點的經(jīng)度對應(yīng) λ1、λ2，再利用積分的原理來計算得到橢球體上的該區(qū)域面積。

程序如下：

先聲明和賦值程序中將會用到的基本數(shù)據(jù)長半軸a、第一偏心率e、基準緯度（江蘇省最低緯度）B；并通過指向矢量文件的指針獲得前后兩點的緯度B1、B2和經(jīng)度L1、L2。用TB來代替2-1，用AB來代替（1+2）/2：

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

?

4、把該區(qū)域進行蘭伯特投影得到二維平面上又一系列的區(qū)域邊界點數(shù)據(jù)

這里參考《.gen地圖文件的投影編程實現(xiàn)（以墨卡托投影和蘭伯特投影為例）》http://www.cnblogs.com/to-sunshine/p/6048438.html。

?

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

方法同第一步，只是帶入的數(shù)據(jù)不同。

?

6、對比三種面積的區(qū)別

整個程序主函數(shù)如下：

運行后結(jié)果如下：

可見，進行Albers等積投影之后，矢量面的面積變化誤差相比之整體面積來說較小，所以視為等積投影是成功的。

?

轉(zhuǎn)載于:https://www.cnblogs.com/to-sunshine/p/7655432.html

總結(jié)

以上是生活随笔為你收集整理的椭球体上某区域面积的求算，及该区域兰伯特投影与墨卡托投影到二维平面后面积对比...的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。