椭球体上某区域面积的求算，及该区域兰伯特投影与墨卡托投影到二维平面后面积对比

    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等积投影以后，矢量面的面积变化偏差相比之总体面积来讲较小，因此视为等积投影是成功的。