基于ArcGIS和Python的图斑椭球面积计算方法

慧谷森林

1.前言作为一名从事规划行业的人来讲，计算图斑的椭球面积在工作中是必不可少的，熟悉ArcGIS的人一般会使用下面这些步骤来计算（我的软件版本为ArcGIS10.2.2forDesktop）：打开图层属性表；右键要计算椭球面积的字段，进入字段计算器；解析程序选择Python，并在下方输入【!shape.geodesicArea!】，如果需要保留两位小数则输入【round(!shape.geodesicArea!,2)】，点击确定即可计算出椭球面积。但在实际工作中，难免会遇到面积很小的图斑，而这些面积很小的图斑用上述步骤计算出来可能会出现面积为0的情况。在最初的时候，这种情况令我百思不得其解，不过我还是找到了一个解决办法，就是去ArcGISPro中计算几何（选择测地线面积），这种方法可以保证不管图斑有多小都会有数值计算出来，就在我以为问题解决了的时候，新的问题又出现了。如下：上面的图片是用ArcGIS10.2.2计算的，下面的图片是用ArcGISPro3.0.0计算的，可以看到面积相差1.2平方米左右，经过核查，我发现这次是ArcGIS10.2.2算的准，而ArcGISPro却算的不准了，而差的这1.2平方米在数据检查时是会报错的。这种情况勾起了我的好奇心，想要知道这是为什么，因此我去网上查找有没有相关经验，结果是没有。我又去找有没有计算椭球面积的小工具，但结果是虽然有很多人分享工具和代码出来，但工具大多是要付费的，代码也不是我熟悉的Python语言。所以这次只能自己写一个计算椭球面积的代码，前后加起来大概花了一天时间完成（在此感谢空间规划工具箱作者@规划酱大大对我的帮助与解惑），经过测试，椭球面积达到了准确无误，源代码免费分享出来，希望能帮到和我一样有困惑的人。2.源代码分享椭球面积计算依据：《TD/T1055-2019第三次全国国土调查技术规程》中的附录D-图斑椭球面积计算公式及要求。编写环境：Anaconda3(64-bit)、Python2.7、ArcGIS10.2.2。最终做成一个ArcGIS的工具箱，代码如下（由于是我个人使用，函数和变量的命名方式不喜勿喷，我认为有必要加注释的地方都加了注释）：importarcpyimportmath#常数pi=3.14159265358979ipi=pi/180.0p=206264.8062471#中央经线弧度L0=114.0*ipi#CGCS2000椭球常数如下a=6378137.0#椭球长半轴f=1.0/298.257222101#椭球扁率b=6356752.31414036#椭球短半轴e1=0.0066943800229#椭球第一偏心率e^2ee=0.00673949677548#椭球第二偏心率e`^2c=6399593.62586#极点子午圈曲率半径#D.3其他相关常数如下K0=1.57048761144159*math.pow(10.0,-7.0)K1=5.05250178820567*math.pow(10.0,-3.0)K2=2.98472900956587*math.pow(10.0,-5.0)K3=2.41626669230084*math.pow(10.0,-7.0)K4=2.22241238938534*math.pow(10.0,-9.0)#D.1其他相关常数如下KA=1.0+3.0/6.0*e1+30.0/80.0*e1*e1+35.0/112.0*e1*e1*e1+630.0/2304.0*e1*e1*e1*e1KB=1.0/6.0*e1+15.0/80.0*e1*e1+21.0/112.0*e1*e1*e1+420.0/2304.0*e1*e1*e1*e1KC=3.0/80.0*e1*e1+7.0/112.0*e1*e1*e1+180.0/2304.0*e1*e1*e1*e1KD=1.0/112.0*e1*e1*e1+45.0/2304.0*e1*e1*e1*e1KE=5.0/2304.0*e1*e1*e1*e1#D.3高斯投影反解变换（x,y->B,L）公式defXY2BL(x,y):y1=y-500000-38*1000000#坐标带号E=K0*xBf=E+math.cos(E)*(K1*math.sin(E)-K2*math.sin(E)*math.sin(E)*math.sin(E)+K3*math.sin(E)*math.sin(E)*math.sin(E)*math.sin(E)*math.sin(E)-K4*math.sin(E)*math.sin(E)*math.sin(E)*math.sin(E)*math.sin(E)*math.sin(E)*math.sin(E))t=math.tan(Bf)nn=ee*math.cos(Bf)*math.cos(Bf)C=a*a/bV=math.sqrt(1.0+nn)N=C/VVVT=V*V*tyn=y1/Nynn=(y1/N)*(y1/N)tt=t*tCBf=1.0/math.cos(Bf)B=Bf-1.0/2.0*VVT*ynn+1.0/24.0*(5.0+3.0*tt+nn-9.0*nn*tt)*VVT*ynn*ynn-1.0/720.0*(61.0+90.0*tt+45.0*tt*tt)*VVT*ynn*ynn*ynnL=CBf*yn-1.0/6.0*(1.0+2.0*tt+nn)*CBf*yn*yn*yn+1.0/120.0*(5.0+28.0*tt+24.0*tt*tt+6.0*nn+8.0*nn*tt)*CBf*yn*yn*yn*yn*yn+L0returnB,L#D.1椭球面上任一梯形图块面积计算公式deftx_area(B1,L1,B2,L2):BM=(B1+B2)/2.0BC=B2-B1LC=(L1+L2)/2.0S=2.0*b*b*LC*(KA*math.sin(0.5*BC)*math.cos(BM)-KB*math.sin(1.5*BC)*math.cos(3.0*BM)+KC*math.sin(2.5*BC)*math.cos(5.0*BM)-KD*math.sin(3.5*BC)*math.cos(7.0*BM)+KE*math.sin(4.5*BC)*math.cos(9.0*BM))returnS#检查线段是否大于70米defcheck_len(start_x,start_y,end_x,end_y):#计算原始线段的长度original_length=math.sqrt((end_x-start_x)**2+(end_y-start_y)**2)#如果线段长度小于等于70m，则无需分割，直接返回起点和终点iforiginal_length椭球面积计算完成...')有疑问可后台私信或者评论。