如何使用Cartopy结合AACGMv2绘制磁坐标系下的海岸线等地理要素

你最初的思路是完全可落地的，也是目前实现该需求的最优方案，整体实现成本低、精度可控，不需要依赖额外的第三方工具。

核心实现逻辑 #

先提取cartopy内置海岸线的地理坐标点，批量送入AACGMv2转换为磁经纬度，最后按原始的线段组合直接绘制转换后的坐标即可。

具体操作步骤 #

• 提取cartopy海岸线的所有坐标段：遍历cartopy.feature.COASTLINE的geometries接口，把所有线要素的(lon, lat)坐标对按原始分段提取出来，提前过滤±89°以外超出AACGMv2转换范围的无效点。
• 批量完成坐标转换：调用aacgmv2的convert()方法，批量传入提取到的纬度、经度数组，指定转换对应的年份（适配地磁漂移参数），返回对应的磁纬、磁经数组，按原来的线段分组保存。
• 绘制磁坐标系海岸线：如果不需要叠加其他地理坐标系要素，直接用普通matplotlib axes逐段绘制转换后的坐标即可；如果需要极射赤面投影等特殊展示效果，直接对转换后的磁坐标做常规的投影变换即可。

示例代码片段 #

import aacgmv2
import numpy as np
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
from datetime import datetime

# 提取所有海岸线坐标段
coast_segments = []
for geom in cfeature.COASTLINE.geometries():
    if geom.geom_type == 'LineString':
        coast_segments.append(np.array(list(geom.coords)))
    elif geom.geom_type == 'MultiLineString':
        for line in geom.geoms:
            coast_segments.append(np.array(list(line.coords)))

# 批量转换为AACGMv2磁坐标，指定转换年份为2020年
mag_segments = []
for seg in coast_segments:
    lon, lat = seg[:, 0], seg[:, 1]
    # 过滤超出转换范围的点
    valid_mask = (lat >= -89) & (lat <= 89)
    if np.sum(valid_mask) < 2:
        continue
    mlat, mlon, _ = aacgmv2.convert(lat[valid_mask], lon[valid_mask], 0, datetime(2020, 1, 1))
    mag_segments.append(np.column_stack([mlon, mlat]))

# 绘制磁坐标系海岸线
fig, ax = plt.subplots(figsize=(8, 8))
for seg in mag_segments:
    ax.plot(seg[:, 0], seg[:, 1], color='k', linewidth=0.5)
ax.set_xlabel('磁经度 (°)')
ax.set_ylabel('磁纬度 (°)')
plt.show()

可选高阶方案 #

如果需要保留cartopy的全部投影能力、叠加其他地理要素，可以自定义cartopy的CRS类，把AACGMv2的转换逻辑嵌入到坐标转换回调中，不过该方案实现复杂度较高，没有特殊需求不推荐使用。

内容的提问来源于stack exchange，提问作者Nathaniel F.