arcgis 椭球面积计算公式

图幅理论面积与图斑椭球面积计算公式及要求一、 图幅理论面积计算公式 （1）式中：a—椭球长半轴(单位：米)，α—椭球扁率，b—椭球短半轴(单位：米)е²﹦(a²﹣b²)/a²A﹦1﹢（3/6）е²﹢（30/80）е4﹢（35/112）е6﹢（630/2304）е8B﹦ （1/6）е²﹢（15/80）е4﹢（21/112）е6﹢（420/2304）е8C﹦ （3/80）е4﹢ （7/112）е6﹢（180/2304）е8D﹦ （1/112）е6﹢ （45/2304）е8E﹦ （5/2304）е8ΔL—图幅东西图廓的经差（单位：分）B2﹣B1)—图幅南北图廓的纬差（单位：弧度），Bm﹦（B1﹢B2）/2二、椭球面上任意梯形面积计算公式 （2）其中：A,B,C,D,E 为常数，按下式计算:е²﹦(a²﹣b²)/a²A﹦1﹢（3/6）е²﹢（30/80）е4﹢（35/112）е6﹢（630/2304）е8B﹦ （1/6）е²﹢（15/80）е4﹢（21/112）е6﹢（420/2304）е8C﹦ （3/80）е4﹢ （7/112）е6﹢（180/2304）е8D﹦ （1/112）е6﹢（45/2304）е8E﹦ （5/2304）е8式中：a—椭球长半轴(单位：米)，b—椭球短半轴(单位：米)；ΔL—图块经差(单位：弧度)； (B2﹣B1)—图块纬差(单位：弧度) Bm﹦（B1﹢B2）/2。

三、高斯投影反解变换（）模型 （若坐标不带带号，则不需减去带号×1000000；） +中央子午线经度值（弧度） （3）式中： 公式说明：若坐标为没有带号前缀格式，则不需减去带号×1000000；若坐标为有带号前缀格式，则需减去带号×1000000四、计算用到的常数、椭球参数在计算图幅理论面积与任意图斑椭球面积时，有关常数及保留的位数按给定数值计算常数：π﹦3.14159265358979 206264.806247180椭球常数： = 6378140 = 1/ 298.257 = 6356755.29 = 6.69438499958795E-03 = 6.73950181947292E-03 = 6399596.65198801相关常数: k0 = 1.57048687472752E-07k1 = 5.05250559291393E-03k2 = 2.98473350966158E-05k3 = 2.41627215981336E-07k4 = 2.22241909461273E-09五、计算中的取位及要求①高斯投影反解变换后的B,L保留到秒后6位小数，四舍五入。

② 采用计算机计算时，所有变量数据类型均要定义为双精度③ 面积计算结果以平方米为单位，保留一位小数，四舍五入④ 各种比例尺标准分幅图经差、纬差见表1⑤ 在用大地坐标生成标准分幅图框时，要求在每条边框线的整秒处插入加密点表1 各种比例尺标准分幅图经差、纬差表比例尺1:100万1:50万1:25万1:10万1:5万1:2.5万1:1万1:5千经差6º3º1º30′30′15′7′30″3′45″1′52.5″纬差4º2º1º20′10′5′2′30″1′15″六、任意图斑椭球面积计算方法任意封闭图斑椭球面积计算的原理：将任意封闭图斑高斯平面坐标利用高斯投影反解变换模型，将高斯平面坐标换算为相应椭球的大地坐标，再利用椭球面上任意梯形图块面积计算模型计算其椭球面积，从而得到任意封闭图斑的椭球面积1、计算方法：任意封闭区域总是可以分割成有限个任意小的梯形图块，因此，任意封闭区域的面积，式中Si为分割的任意小的梯形图块面积（i=1,2,…n）用公式(2)计算求封闭区域（多边形如图1）ABCD的面积 ，其具体方法为：（1）对封闭区域（多边形）的界址点连续编号（顺时针或逆时针）ABCD，提取各界址点的高斯平面坐标A(X1,Y1)，B(X2,Y2)，C(X3,Y3)，D(X4,Y4)；（2）利用高斯投影反解变换模型公式（3），将高斯平面坐标换算为相应椭球的大地坐标A(B1,L1)，B(B2,L2)，C(B3,L3)，D(B4,L4)；（3）任意给定一经线L0（如L0＝60°），这样多边形ABCD的各边AB、BC、CD、DA与L0就围成了4个梯形图块（ABB1A1、BCC1B1、CDD1C1、DAA1D1）； （4）由于在椭球面上同一经差随着纬度升高，梯形图块的面积逐渐减小，而同一纬差上经差梯形图块的面积相等，所以，将梯形图块ABB1A1按纬差分割成许多个小梯形图块AEiFiA1，用公式（2）计算出各小梯形图块AEiFiA1的面积Si，然后累加Si就得到梯形图块ABB1A1的面积，同理，依次计算出梯形图块BCC1B1、CDD1C1、DAA1D1的面积（注：用公式（2）计算面积时，B1、B2分别取沿界址点编号方向的前一个、后一个界址点的大地纬度，ΔL为沿界址点编号方向的前一个、后一个界址点的大地经度的平均值与L0的差）；（5）多边形ABCD的面积就等于4个梯形图块（ABB1A1、BCC1B1、CDD1C1、DAA1D1）面积的代数和。

C(B3,L3)D(B4,L4)B (B2,L2)A(B1,L1)LL0C1D1A1FiBEi(Bi,Li)图1 椭球面上任意多边形计算面积则任意多边形ABCD的面积P为：P=ABCD= BCC1B1+ CDD1C1+ DAA1D1- ABB1A12、计算要求① 利用图形坐标点将高斯坐标系下的几何图形反算投影到大地坐标系，进行投影变换② 任意指定一条经线L0，从选定多边形几何形状的起始点开始，沿顺时针方向依次计算相邻两点构成的线段，以及两点到指定经线的平行线构成的梯形面积③ 计算过程中应顺同一方向依坐标点逐个计算相邻两点连线与任意经线构成的梯形面积，坐标点不得有遗漏若多边形包含内多边形（洞），则该多边形面积为外多边形面积减去所有内多边形面积之和④ 计算所有梯形面积的代数和即为该多边形的面积七、算法伪代码描述为了确保编程使用的参数、算法一致，保证不同软件计算的椭球面积一致，我们用算法伪代码描述的方法对编程进行统一，在利用计算机编制椭球面积计算软件时，计算参数与计算顺序应严格按照以下代码执行1、概述计算规则：两个绝对值很大的数或两个绝对值很小的数相乘时，不能用幂计算符，并且尽量不使用连乘，应使用下述的方法：A*B*A*B；数据类型：当使用.net环境时，用Decimal数据类型代替Double类型；2、参数说明双精度类型： 圆周率值：PI = 3.14159265358979 中央经线:CenterL RHO = 206264.8062471 A:ParamA B:ParamB C:ParamC D:ParamD E:ParamE Const ZERO As Double = 0.000000000001 80椭球常数 椭球长半轴：aRadius = 6378140 椭球短半轴：bRadius = 6356755.29 椭球扁率：ParaAF = 1/ 298.257 椭球第一偏心率：ParaE1 = 6.69438499958795E-03 椭球第二偏心率：ParaE2 = 6.73950181947292E-03 极点子午圈曲率半径：ParaC = 6399596.65198801 k0:Parak0 = 1.57048687472752E-07 k1:Parak1 = 5.05250559291393E-03 k2:Parak2 = 2.98473350966158E-05 k3:Parak3 = 2.41627215981336E-07 k4:Parak4 = 2.22241909461273E-09 3、算法描述 初始化参数 Double e; Double a; e = ParaE1; ParamA = 1 + (3 / 6) * e + (30 / 80) * Power(e, 2) + (35 / 112) * Power(e, 3) + (630 / 2304) * Power(e, 4); ParamB = (1 / 6) * e + (15 / 80) * Power(e, 2) + (21 / 112) * Power(e, 3) + (420 / 2304) * Power(e, 4); ParamC = (3 / 80) * Power(e, 2) + (7 / 112) * Power(e, 3) + (180 / 2304) * Power(e, 4); ParamD = (1 / 112) * Power(e, 3) + (45 / 2304) * Power(e, 4); ParamE = (5 / 2304) * Power(e, 4); 参数初始化结束 中央经线转换为弧度 CenterL = TransDegreeToArc(CenterL)选定本初子午线为参考经线 StandardLat = 0 For 起始点 To 倒数第二点 由高斯坐标反解计算经纬度值 ComputeXYGeo (PntColl.Point(i).y, PntColl.Point(i).x, B, L, CenterL) ComputeXYGeo (PntColl.Point(i + 1).y, PntColl.Point(i + 1).x, B1, L1, CenterL) 将经纬度转换为弧度值 B = B / RHO L = L / RHO B1 = B1 / RHO L1 = L1 / RHO 计算梯形面积 Double AreaVal;//梯形面积值 Double lDiference ;//经差 Double bDiference; //纬差 Double bSum;//纬度和 Double ItemValue(5);//计算变量 bDiference = (B1 - B0) / 2； bSum = (B1 + B0) / 2； lDiference = (L1 + L) / 2; //按照以下计算顺序：短半径*经差*短半径 RadDiffVal = 2 * bRadius * lDiference * bRadius cosVal = Cos(bSum) sinVal = Sin(bDiference) ItemValue(0) = RadDiffVal * ParamA * cosVal * sinVal ItemValue(1) = RadDiffVal * ParamB * Sin(3 * bDiference) * Cos(3 * bSum) ItemValue(2) = RadDiffVal * ParamC * Sin(5 * bDiference) * Cos(5 * bSum) ItemValue(3) = RadDiffVal * ParamD * Sin(7 * bDiference) * Cos(7 * bSum) ItemValue(4) = RadDiffVal * ParamE * Sin(9 * bDiference) * Cos(9 * bSum) AreaVal = ItemValue(0) - ItemValue(1) + ItemValue(2) - ItemValue(3) + ItemValue(4) areaSum = areaSum + AreaVal; Next End Sub 4、高斯坐标反解算法Public Sub ComputeXYGeo(x As Double, y As Double, B As Double, L As Double, center As Double) Dim y1 As Double Dim bf As Double y1 = y - 500000 Dim e As Double e = Parak0 * x Dim se As Double se = Sin(e)bf = e + Cos(e) * (Parak1 * se - Parak2 * Power(se, 3) + Parak3 * Power(se, 5) - Parak4 * Power(se, 7)) Dim v As Double Dim t As Double Dim N As Double Dim nl As Double Dim vt As Double Dim yn As Double Dim t2 As Double Dim g As Double g = 1 t = Tan(bf) nl = ParaE2 * Power(Cos(bf), 2) v = Sqr(1 + nl) N = ParaC / v yn = y1 / N vt = Power(v, 2) * t t2 = Power(t, 2) B = bf - vt * Power(yn, 2) / 2 + (5 + 3 * t2 + nl - 9 * nl * t2) * vt * Power(yn, 4) / 24 - (61 + 90 * t2 + 45 * Power(t2, 2)) * vt * Power(yn, 6) / 720 B = TransArcToDegree(B) Dim cbf As Double cbf = 1 / Cos(bf) L = cbf * yn - (1 + 2 * t2 + nl) * cbf * Power(yn, 3) / 6 + (5 + 28 * t2 + 24 * Power(t2, 2) + 6 * nl + 8 * nl * t2) * cbf * Power(yn, 5) / 120 + center L = TransArcToDegree(L) End Sub 弧度转换为秒Public Function TransArcToDegree(arc As Double) As Double Dim sec as Double sec = arc * RHO //秒保留到小数点后6位，四舍五入 sec = Format(sec, "####.000000") TransArcToDegree = sec End Function9。