import numpy as np
def ellipsoidal_trapezoid_area(a, b, phi1_deg, phi2_deg, delta_L_deg, is_map_sheet=False):
"""
计算椭球面上梯形面积的通用公式和图幅公式
参数:
a: 椭球长半轴（米）
b: 椭球短半轴（米）
phi1_deg: 起始纬度（度）
phi2_deg: 终止纬度（度）
delta_L_deg: 经差（度）
is_map_sheet: 是否使用图幅公式（默认False，使用通用公式）
返回:
梯形面积（平方米）
"""
# 转换为弧度
phi1 = np.radians(phi1_deg)
phi2 = np.radians(phi2_deg)
# 计算第一偏心率平方 e²
e_sq = 1 - (b**2 / a**2)
e = np.sqrt(e_sq)
# 定义被积函数的两部分
def term1(phi):
return (np.sin(phi) * np.cos(phi)) / (1 - e_sq * np.sin(phi)**2)
def term2(phi):
return (1 / (2 * e)) * np.log((1 + e * np.sin(phi)) / (1 - e * np.sin(phi)))
# 计算积分结果
integral = (term1(phi2) + term2(phi2)) - (term1(phi1) + term2(phi1))
# 选择公式
if is_map_sheet:
# 图幅公式：经差需以“分”为单位，系数为 4πb²/(360*60)
delta_L_min = delta_L_deg * 60 # 度转分
area = (4 * np.pi * b**2 / (360 * 60)) * delta_L_min * integral
else:
# 通用公式：经差以弧度为单位
delta_L_rad = np.radians(delta_L_deg)
area = (b**2) * delta_L_rad * integral
return area
# 示例：WGS84 椭球参数
a = 6378137.0 # 长半轴（米）
b = 6356752.314245 # 短半轴（米）
# 案例1：计算纬度 20°-25°、经差 1° 的梯形面积
phi1_deg = 20.0
phi2_deg = 25.0
delta_L_deg = 1.0
# 通用公式计算
area_general = ellipsoidal_trapezoid_area(a, b, phi1_deg, phi2_deg, delta_L_deg)
# 图幅公式计算（经差需为分数，此处1°=60'）
area_map_sheet = ellipsoidal_trapezoid_area(a, b, phi1_deg, phi2_deg, delta_L_deg, is_map_sheet=True)
print("通用公式结果:")
print(f"面积 = {area_general:.2f} 平方米 | {area_general / 1e6:.2f} 平方公里")
print("\n图幅公式结果 (经差1°=60'):")
print(f"面积 = {area_map_sheet:.2f} 平方米 | {area_map_sheet / 1e6:.2f} 平方公里")
print(f"\n差值 = {abs(area_general - area_map_sheet):.2f} 平方米")
# 案例2：标准1:1万图幅（经差3.75'，纬差2.5'）
delta_L_deg_map = 3.75 / 60 # 3.75分转度
delta_phi_deg_map = 2.5 / 60 # 2.5分转度
phi_start = 30.0 # 起始纬度30°
phi_end = phi_start + delta_phi_deg_map
area_1_10000 = ellipsoidal_trapezoid_area(a, b, phi_start, phi_end, delta_L_deg_map, is_map_sheet=True)
print("\n1:1万标准图幅面积 (30°N, 经差3.75', 纬差2.5'):")
print(f"面积 = {area_1_10000:.2f} 平方米 | {area_1_10000 / 1e6:.4f} 平方公里")
