如何用Python（优先Cartopy）计算椭球体经纬度网格单元面积？

优先方案：使用Cartopy的Geodesic类 #

Cartopy自带的cartopy.geodesic.Geodesic可直接基于自定义椭球体计算地理多边形面积，完美适配规则经纬度网格场景，精度高且实现简洁。

代码示例

import numpy as np
import cartopy.crs as ccrs
from cartopy.geodesic import Geodesic

# 自定义椭球体参数（替换为你的a、b值，单位：米）
a = 6378137.0
b = 6356752.314245

# 定义经纬度网格参数
lon_min, lon_max, d_lon = -180, 180, 1.0  # 经度范围、步长
lat_min, lat_max, d_lat = -90, 90, 1.0    # 纬度范围、步长

# 生成网格点数组
lons = np.arange(lon_min, lon_max + d_lon, d_lon)
lats = np.arange(lat_min, lat_max + d_lat, d_lat)

# 初始化自定义椭球体的Geodesic对象
geod = Geodesic(semi_major_axis=a, semi_minor_axis=b)

# 创建存储网格面积的数组（维度为网格行数×列数）
grid_areas = np.zeros((len(lats)-1, len(lons)-1))

# 遍历每个网格单元计算面积
for i in range(len(lats)-1):
    for j in range(len(lons)-1):
        # 获取当前网格的四个顶点经纬度（按顺时针顺序排列）
        vertices = [
            (lons[j], lats[i]),
            (lons[j+1], lats[i]),
            (lons[j+1], lats[i+1]),
            (lons[j], lats[i+1])
        ]
        # 计算面积（返回值单位：平方米），取绝对值确保为正
        area_sqm = abs(geod.geometry_area_perimeter(vertices)[0])
        # 转换为平方公里并存储
        grid_areas[i, j] = area_sqm / 1e6

关键注意事项

• 多边形顶点必须按连续的顺时针或逆时针顺序传入，否则面积计算会出错
• geometry_area_perimeter方法返回的第一个值是面积，第二个是周长，仅需提取面积部分
• 面积默认单位为平方米，可根据需求转换为平方公里（除以1e6）或其他单位

替代方案：使用pyproj的Geod类 #

若不想依赖Cartopy，可直接用pyproj计算（Cartopy底层也依赖pyproj），原理一致：

import numpy as np
from pyproj import Geod

# 初始化自定义椭球体的Geod对象
geod = Geod(a=a, b=b)

grid_areas = np.zeros((len(lats)-1, len(lons)-1))

for i in range(len(lats)-1):
    for j in range(len(lons)-1):
        # 提取网格顶点的经纬度数组
        poly_lons = [lons[j], lons[j+1], lons[j+1], lons[j]]
        poly_lats = [lats[i], lats[i], lats[i+1], lats[i+1]]
        # 计算面积
        area_sqm, _ = geod.polygon_area_perimeter(poly_lons, poly_lats)
        grid_areas[i, j] = abs(area_sqm) / 1e6

内容的提问来源于stack exchange，提问作者DingoTim