arcgis 橢球面積計(jì)算公式

1、圖幅理論面積與圖斑橢球面積計(jì)算公式及要求一、 圖幅理論面積計(jì)算公式 （1）式中：a橢球長(zhǎng)半軸(單位：米)，橢球扁率，b橢球短半軸(單位：米)。²(a²b²)/a²。A1（3/6）²（30/80）4（35/112）6（630/2304）8。B （1/6）²（15/80）4（21/112）6（420/2304）8。C （3/80）4 （7/112）6（180/2304）8。D （1/112）6 （45/2304）8。E （5/2304）8。L圖幅東西圖廓的經(jīng)差（單位：分）。(B2B1)圖幅南北圖廓的緯差（單位：弧度），Bm（B1B2）/2

2、。二、橢球面上任意梯形面積計(jì)算公式 （2）其中：A,B,C,D,E 為常數(shù)，按下式計(jì)算:²(a²b²)/a²A1（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橢球長(zhǎng)半軸(單位：米)，b橢球短半軸(單位：米)；L圖塊經(jīng)差(單位：弧度)； (B2B1)圖塊緯差(單位：弧度) Bm（B1B2）/2。三、高斯投影反解變換（）

3、模型 （若坐標(biāo)不帶帶號(hào)，則不需減去帶號(hào)×1000000；） +中央子午線經(jīng)度值（弧度） （3）式中： 公式說明：若坐標(biāo)為沒有帶號(hào)前綴格式，則不需減去帶號(hào)×1000000；若坐標(biāo)為有帶號(hào)前綴格式，則需減去帶號(hào)×1000000。四、計(jì)算用到的常數(shù)、橢球參數(shù)在計(jì)算圖幅理論面積與任意圖斑橢球面積時(shí)，有關(guān)常數(shù)及保留的位數(shù)按給定數(shù)值計(jì)算。常數(shù)： 206264.806247180橢球常數(shù)： = 6378140 = 1/ 298.257 = 6356755.29 = 6.69438499958795E-03 = 6.73950181947292E-03 = 6399596.651

4、98801相關(guān)常數(shù): k0 = 1.57048687472752E-07k1 = 5.05250559291393E-03k2 = 2.98473350966158E-05k3k4 = 2.22241909461273E-09五、計(jì)算中的取位及要求高斯投影反解變換后的B,L保留到秒后6位小數(shù)，四舍五入。 采用計(jì)算機(jī)計(jì)算時(shí)，所有變量數(shù)據(jù)類型均要定義為雙精度。 面積計(jì)算結(jié)果以平方米為單位，保留一位小數(shù)，四舍五入。 各種比例尺標(biāo)準(zhǔn)分幅圖經(jīng)差、緯差見表1。 在用大地坐標(biāo)生成標(biāo)準(zhǔn)分幅圖框時(shí)，要求在每條邊框線的整秒處插入加密點(diǎn)。表1 各種比例尺標(biāo)準(zhǔn)分幅圖經(jīng)差、緯差表比例尺1:100萬1:50萬1:25萬1

5、:10萬1:5萬1:2.5萬1:1萬1:5千經(jīng)差6º3º1º303015730345152.5緯差4º2º1º20105230115六、任意圖斑橢球面積計(jì)算方法任意封閉圖斑橢球面積計(jì)算的原理：將任意封閉圖斑高斯平面坐標(biāo)利用高斯投影反解變換模型，將高斯平面坐標(biāo)換算為相應(yīng)橢球的大地坐標(biāo)，再利用橢球面上任意梯形圖塊面積計(jì)算模型計(jì)算其橢球面積，從而得到任意封閉圖斑的橢球面積。1、計(jì)算方法：任意封閉區(qū)域總是可以分割成有限個(gè)任意小的梯形圖塊，因此，任意封閉區(qū)域的面積，式中Si為分割的任意小的梯形圖塊面積（i=1,2,n）用公式(2)計(jì)算。求封閉區(qū)

6、域（多邊形如圖1）ABCD的面積 ，其具體方法為：（1）對(duì)封閉區(qū)域（多邊形）的界址點(diǎn)連續(xù)編號(hào)（順時(shí)針或逆時(shí)針）ABCD，提取各界址點(diǎn)的高斯平面坐標(biāo)A(X1,Y1)，B(X2,Y2)，C(X3,Y3)，D(X4,Y4)；（2）利用高斯投影反解變換模型公式（3），將高斯平面坐標(biāo)換算為相應(yīng)橢球的大地坐標(biāo)A(B1,L1)，B(B2,L2)，C(B3,L3)，D(B4,L4)；（3）任意給定一經(jīng)線L0（如L060°），這樣多邊形ABCD的各邊AB、BC、CD、DA與L0就圍成了4個(gè)梯形圖塊（ABB1A1、BCC1B1、CDD1C1、DAA1D1）； （4）由于在橢球面上同一經(jīng)差隨著緯度升高，梯

7、形圖塊的面積逐漸減小，而同一緯差上經(jīng)差梯形圖塊的面積相等，所以，將梯形圖塊ABB1A1按緯差分割成許多個(gè)小梯形圖塊AEiFiA1，用公式（2）計(jì)算出各小梯形圖塊AEiFiA1的面積Si，然后累加Si就得到梯形圖塊ABB1A1的面積，同理，依次計(jì)算出梯形圖塊BCC1B1、CDD1C1、DAA1D1的面積（注：用公式（2）計(jì)算面積時(shí)，B1、B2分別取沿界址點(diǎn)編號(hào)方向的前一個(gè)、后一個(gè)界址點(diǎn)的大地緯度，L為沿界址點(diǎn)編號(hào)方向的前一個(gè)、后一個(gè)界址點(diǎn)的大地經(jīng)度的平均值與L0的差）；（5）多邊形ABCD的面積就等于4個(gè)梯形圖塊（ABB1A1、BCC1B1、CDD1C1、DAA1D1）面積的代數(shù)和。C(B3,L

8、3)D(B4,L4)B (B2,L2)A(B1,L1)LL0C1D1A1FiBEi(Bi,Li)圖1 橢球面上任意多邊形計(jì)算面積則任意多邊形ABCD的面積P為：P=ABCD= BCC1B1+ CDD1C1+ DAA1D1- ABB1A12、計(jì)算要求 利用圖形坐標(biāo)點(diǎn)將高斯坐標(biāo)系下的幾何圖形反算投影到大地坐標(biāo)系，進(jìn)行投影變換。 任意指定一條經(jīng)線L0，從選定多邊形幾何形狀的起始點(diǎn)開始，沿順時(shí)針方向依次計(jì)算相鄰兩點(diǎn)構(gòu)成的線段，以及兩點(diǎn)到指定經(jīng)線的平行線構(gòu)成的梯形面積。 計(jì)算過程中應(yīng)順同一方向依坐標(biāo)點(diǎn)逐個(gè)計(jì)算相鄰兩點(diǎn)連線與任意經(jīng)線構(gòu)成的梯形面積，坐標(biāo)點(diǎn)不得有遺漏。若多邊形包含內(nèi)多邊形（洞），則該多邊形面

9、積為外多邊形面積減去所有內(nèi)多邊形面積之和。 計(jì)算所有梯形面積的代數(shù)和即為該多邊形的面積。七、算法偽代碼描述為了確保編程使用的參數(shù)、算法一致，保證不同軟件計(jì)算的橢球面積一致，我們用算法偽代碼描述的方法對(duì)編程進(jìn)行統(tǒng)一，在利用計(jì)算機(jī)編制橢球面積計(jì)算軟件時(shí)，計(jì)算參數(shù)與計(jì)算順序應(yīng)嚴(yán)格按照以下代碼執(zhí)行。1、概述計(jì)算規(guī)則：兩個(gè)絕對(duì)值很大的數(shù)或兩個(gè)絕對(duì)值很小的數(shù)相乘時(shí)，不能用冪計(jì)算符，并且盡量不使用連乘，應(yīng)使用下述的方法：A*B*A*B；數(shù)據(jù)類型：當(dāng)使用.net環(huán)境時(shí)，用Decimal數(shù)據(jù)類型代替Double類型；2、參數(shù)說明雙精度類型：中央經(jīng)線:CenterL RHO = 206264.8062471A:P

10、aramAB:ParamBC:ParamCD:ParamDE:ParamEConst ZERO As Double = 0.000000000001 80橢球常數(shù)橢球長(zhǎng)半軸：aRadius = 6378140 橢球短半軸：bRadius = 6356755.29橢球扁率：ParaAF = 1/ 298.257橢球第一偏心率：ParaE1 = 6.69438499958795E-03橢球第二偏心率：ParaE2 = 6.73950181947292E-03極點(diǎn)子午圈曲率半徑：ParaC = 6399596.65198801 k0:Parak0 = 1.57048687472752E-07k1:P

11、arak1 = 5.05250559291393E-03k2:Parak2 = 2.98473350966158E-05k4:Parak4 = 2.22241909461273E-09 3、算法描述初始化參數(shù) 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 / 1

12、12) * 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); 參數(shù)初始化結(jié)束中央經(jīng)線轉(zhuǎn)換為弧度 CenterL = TransDegreeToArc(CenterL)選定本初子午線為參考經(jīng)線 Sta

13、ndardLat = 0 For 起始點(diǎn) To 倒數(shù)第二點(diǎn) 由高斯坐標(biāo)反解計(jì)算經(jīng)緯度值 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) 將經(jīng)緯度轉(zhuǎn)換為弧度值 B = B / RHO L = L / RHO B1 = B1 / RHO L1 = L1 / RHO 計(jì)算梯形面積 Double AreaVal;/梯形面積值Double lDiference

14、 ;/經(jīng)差 Double bDiference; /緯差 Double bSum;/緯度和 Double ItemValue(5);/計(jì)算變量 bDiference = (B1 - B0) / 2； bSum = (B1 + B0) / 2； lDiference = (L1 + L) / 2; /按照以下計(jì)算順序：短半徑*經(jīng)差*短半徑 RadDiffVal = 2 * bRadius * lDiference * bRadius cosVal = Cos(bSum) sinVal = Sin(bDiference)ItemValue(0) = RadDiffVal * ParamA * co

15、sVal * sinValItemValue(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 * bS

16、um) AreaVal = ItemValue(0) - ItemValue(1) + ItemValue(2) - ItemValue(3) + ItemValue(4) areaSum = areaSum + AreaVal; Next End Sub 4、高斯坐標(biāo)反解算法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

17、 = 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 =

18、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