发布时间：2025-06-24 20:35:04  作者：北方职教升学中心  阅读量：559

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

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

import arcpyimport math# 常数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）公式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'> 椭球面积计算完成...')

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