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

1.前言

作为一名从事规划行业的人来讲，计算图斑的椭球面积在工作中是必不可少的，熟悉ArcGIS的人一般会使用下面这些步骤来计算（我的软件版本为ArcGIS 10.2.2 for Desktop）：

1. 打开图层属性表；
2. 右键要计算椭球面积的字段，进入字段计算器；
3. 解析程序选择Python，并在下方输入【!shape.geodesicArea!】，如果需要保留两位小数则输入【round(!shape.geodesicArea!,2)】，点击确定即可计算出椭球面积。

但在实际工作中，难免会遇到面积很小的图斑，而这些面积很小的图斑用上述步骤计算出来可能会出现面积为0的情况。

在最初的时候，这种情况令我百思不得其解，不过我还是找到了一个解决办法，就是去ArcGIS Pro中计算几何（选择测地线面积），这种方法可以保证不管图斑有多小都会有数值计算出来，就在我以为问题解决了的时候，新的问题又出现了。如下：

上面的图片是用ArcGIS 10.2.2计算的，下面的图片是用ArcGIS Pro 3.0.0计算的，可以看到面积相差1.2平方米左右，经过核查，我发现这次是ArcGIS 10.2.2算的准，而ArcGIS Pro却算的不准了，而差的这1.2平方米在数据检查时是会报错的。

这种情况勾起了我的好奇心，想要知道这是为什么，因此我去网上查找有没有相关经验，结果是没有。我又去找有没有计算椭球面积的小工具，但结果是虽然有很多人分享工具和代码出来，但工具大多是要付费的，代码也不是我熟悉的Python语言。

所以这次只能自己写一个计算椭球面积的代码，前后加起来大概花了一天时间完成（在此感谢空间规划工具箱作者@规划酱大大对我的帮助与解惑），经过测试，椭球面积达到了准确无误，源代码免费分享出来，希望能帮到和我一样有困惑的人。

2.源代码分享

椭球面积计算依据：《TD/T 1055-2019 第三次全国国土调查技术规程》中的附录D-图斑椭球面积计算公式及要求。

编写环境：Anaconda3 (64-bit)、Python2.7、ArcGIS10.2.2。

最终做成一个ArcGIS的工具箱，代码如下（由于是我个人使用，函数和变量的命名方式不喜勿喷，我认为有必要加注释的地方都加了注释）：

import arcpy
import math
# 常数
pi = 3.14159265358979
ipi = pi / 180.0
p = 206264.8062471
# 中央经线弧度
L0 = 114.0 * ipi
# CGCS2000 椭球常数如下
a = 6378137.0 # 椭球长半轴
f = 1.0 / 298.257222101 #椭球扁率
b = 6356752.31414036 # 椭球短半轴
e1 = 0.0066943800229 # 椭球第一偏心率 e^2
ee = 0.00673949677548 # 椭球第二偏心率 e`^2
c = 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 * e1
KB = 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 * e1
KC = 3.0/80.0 * e1 * e1 + 7.0/112.0 * e1 * e1 * e1 + 180.0/2304.0 * e1 * e1 * e1 * e1
KD = 1.0/112.0 * e1 * e1 * e1 + 45.0/2304.0 * e1 * e1 * e1 * e1
KE = 5.0/2304.0 * e1 * e1 * e1 * e1
 
# D.3 高斯投影反解变换（x,y->B,L）公式
def XY2BL(x,y):
    y1 = y - 500000 - 38 * 1000000 #坐标带号
    E = K0 * x
    Bf  = 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 / b
    V = math.sqrt(1.0 + nn)
    N = C / V
    
    VVT = V * V * t
    yn = y1 / N
    ynn = (y1 / N) * (y1 / N)
    tt = t * t
    CBf = 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 * ynn
    L = 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 + L0
    return B,L
 
# D.1 椭球面上任一梯形图块面积计算公式
def tx_area(B1,L1,B2,L2):
    BM = (B1 + B2) / 2.0
    BC = B2 - B1
    LC = (L1 + L2) / 2.0
    S = 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))
    return S
 
#检查线段是否大于70米
def check_len(start_x, start_y, end_x, end_y):
    # 计算原始线段的长度
    original_length = math.sqrt((end_x - start_x) ** 2 + (end_y - start_y) ** 2)
    # 如果线段长度小于等于70m，则无需分割，直接返回起点和终点
    if original_length <= 70.0:
        return [[start_x, start_y]]
    # 计算需要分割的段数（向上取整，因为最后一段即使小于70m也要单独存在）
    num_segments = int(math.ceil(original_length / 70.0))
    # 计算每段的长度
    segment_length = original_length / num_segments
    # 计算分割点的坐标
    points = [[start_x, start_y]]
    for i in range(1, num_segments):
        # 使用线性插值计算分割点的坐标
        t = float(i) / float(num_segments)
        x = start_x + t * (end_x - start_x)
        y = start_y + t * (end_y - start_y)
        points.append([x, y])
    return points
 
#坐标列表及面积计算
def tqmj(fc):
    with arcpy.da.UpdateCursor(fc, ["SHAPE@",'TQMJ']) as cursor:
        for row in cursor:
            #获取要素初始坐标点列表
            hz = []
            crip_lst = []
            partnum = 1
            for part in row[0]:
                jzd = 1
                for n,point in enumerate(part):
                    start_lst = []
                    if point:
                        start_lst.append(jzd)
                        start_lst.append(partnum)
                        start_lst.append(round(part[n].Y,4))
                        start_lst.append(round(part[n].X,4))
                        jzd += 1
                    else:
                        partnum += 1
                    if start_lst:
                        hz.append(start_lst)
            for n,i in enumerate(hz):
                if hz[n][1] != hz[n-1][1] and n != 0:
                    crip_lst.append(n)
                elif hz[n][0] == 1 and n != 0:
                    crip_lst.append(n)
            for i in hz:
                i.pop(0)
                i.pop(0)
            # 使用循环切分列表
            after_lst = []
            start_index = 0
            for end_index in crip_lst:  
                part = hz[start_index:end_index]
                after_lst.append(part)
                start_index = end_index
            # 最后一部分的长度可能不固定，所以需要单独处理
            if start_index < len(hz):
                after_lst.append(hz[start_index:])
            #检查线段长度是否大于70米并输出最终坐标列表
            end_lst = []
            for y in after_lst:
                end2_lst = []
                for n,i in enumerate(y):
                    if n != len(y) - 2 and n != len(y) - 1:
                        X1,Y1 = y[n][0],y[n][1]
                        X2,Y2 = y[n+1][0],y[n+1][1]
                        for m in check_len(X1,Y1,X2,Y2):
                            end2_lst.append(m)
                    elif n == len(y) - 2:
                        X1,Y1 = y[n][0],y[n][1]
                        X2,Y2 = y[n+1][0],y[n+1][1]
                        for m in check_len(X1,Y1,X2,Y2):
                            end2_lst.append(m)
                        end2_lst.append([X2,Y2])
                end_lst.append(end2_lst)
            #计算面积
            TQMJ = 0
            for y in end_lst:
                for n,i in enumerate(y):
                    if n != len(y) - 2 and n != len(y) - 1:
                        B1,L1 = XY2BL(y[n][0],y[n][1])
                        B2,L2 = XY2BL(y[n+1][0],y[n+1][1])
                        TQ = tx_area(B1,L1,B2,L2)
                        TQMJ += TQ
                    elif n == len(y) - 2:
                        B1,L1 = XY2BL(y[n][0],y[n][1])
                        B2,L2 = XY2BL(y[n+1][0],y[n+1][1])
                        TQ = tx_area(B1,L1,B2,L2)
                        TQMJ += TQ
            if sfxs is True:
                row[1] = round(abs(TQMJ),2)
                cursor.updateRow(row)
            else:
                row[1] = abs(TQMJ)
                cursor.updateRow(row)
if __name__ == '__main__':
    in_fc = arcpy.GetParameterAsText(0) #要素图层
    sfxs = arcpy.GetParameter(1) #是否保留两位小数
    fields = [i.name for i in arcpy.ListFields(in_fc)]
    if 'TQMJ' not in fields:
        arcpy.AddField_management (in_fc, 'TQMJ', field_type = "DOUBLE")
    tqmj(in_fc)
    arcpy.AddMessage (u'> 椭球面积计算完成...')

有疑问可后台私信或者评论。