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)
?? ?
??? # 計算第一偏心率平方 e2
??? 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πb2/(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} 平方公里")
