如何使用 cartopy 在地磁坐标中绘制数据？

Question

我是一名 Python 初学者，我使用 cartopy 将我的数据指向地理坐标。但是，我想在地磁坐标中绘制我的数据。我看到cartopy中有一个旋转杆选项。我不确定的是：如果我使用磁北极的位置（纬度和经度）来旋转两极，那会给我正确的结果吗？因为，北方和南极彼此不对称（由于太阳风的影响），但地理两极是对称的。据我了解，cartopy 中的旋转极点选项会根据给定的新极点位置对称地移动地理极点。

此外，地球上地磁北极的纬度和经度每年都在变化。所以，我想，磁北极的位置应该针对我们感兴趣的每个日期单独更新。这些值可以从多个地方获得，所以这不是问题。

来到我的实际问题：我应该在下面的代码中进行哪些更改才能在地磁坐标中正确绘制我的数据？在这种情况下，cartopy 中的旋转杆选项会起作用吗？也感谢任何有关处理此类问题的其他方法的建议。

这是我用来绘制地理坐标的：

# set the map coverage:
extent = [-90, -60, 30, 60]
central_lon = np.mean(extent[:2])
central_lat = np.mean(extent[2:])

fig2 = plt.figure(figsize=(7, 7))
ax2 = plt.axes(projection=ccrs.Orthographic(central_longitude=central_lon, central_latitude=central_lat))
ax2.set_extent(extent, ccrs.PlateCarree())
ax2.add_feature(cartopy.feature.OCEAN, color='white', alpha=1, zorder=0)
ax2.add_feature(cartopy.feature.LAND, edgecolor='white', 
                color='silver', alpha=0.3, zorder=10)
ax2.add_feature(cartopy.feature.LAKES, color='white', alpha=1, zorder=0)

# plot the pointing direction.
for i in range (0,ind1):
# downsampling to every 10th arrow ([::10])
    L=i*10
# Orientation (X body vector)
    if 360+lon[L] > 360+OLon: # if spacecraft is in the east 
        ax2.quiver(np.array([lon[L]]), np.array([lat[L]]), 
                np.array([-XV[L][2]]), np.array([XV[L][1]]), 
                color='black', scale=7, width=.006, axes=ax2,
                transform=ccrs.PlateCarree(), angles = "xy", zorder=20)
ax2.text(-0.09, 0.5, 'Geographic Latitude', fontsize=14, va='bottom', ha='center',
        rotation='vertical', rotation_mode='anchor',
        transform=ax2.transAxes)
ax2.text(0.5, -0.12, 'Geographic Longitude', fontsize=14, va='bottom', ha='center',
        rotation='horizontal', rotation_mode='anchor',
        transform=ax2.transAxes)

Answer 1

虽然旋转磁极的想法很有趣，但准确的地磁坐标不仅仅是一个不同的磁极位置，因为地球的磁场很奇怪（并且总是在变化）。 IGRF 或 AACGMv2（更精确）可用于从地磁坐标转换为地磁坐标，然后使用旋转磁极方法无需担心磁极不对称 - IGRF/AACGMv2 转换可以解决此问题。

数据（不是Cartopy-shapely对象）可以直接用[lat_mag, long_mag, _] = aacgmv2.convert_latlon_arr(lat, long, alt4mag, time4mag, method_code='G2A'); #converts from geographic to geomagnetic (AACGMv2)进行转换

根据https://github.com/SciTools/cartopy/issues/1945，以下代码将 Cartopy 海岸线从地理坐标转换为地磁坐标并绘制它们（我不太确定如何处理海洋/陆地/湖泊坐标，但应该是类似的方法）：

time4mag 是本地化到 UTC 时区的日期时间对象

alt4mag 是一个整数或浮点数，以千米为单位的高度

注意：对于 geo->mag 坐标以外的其他内容，请将 aacgmv2.convert_latlon_arr 替换为任何其他为特定应用调整纬度/经度坐标的函数

    import cartopy.crs as ccrs
    import cartopy.feature as cfeature
    cc = cfeature.NaturalEarthFeature('physical', 'coastline', '110m', color='xkcd:black', zorder=75); #get a feature to mess with
    geom_mag = []; #prep a list
    for geom in cc.geometries():
        geom_mag.append(convert_to_mag(geom, time4mag, alt4mag));
        # geom_mag.append(geom);
        # coords = list(geom.coords);
    #END FOR geom
    cc_mag = cfeature.ShapelyFeature(geom_mag, ccrs.PlateCarree() color='xkcd:black',zorder=75);
    for geom in cc_mag.geometries():
        ax.plot(*geom.coords.xy, color='xkcd:black', linewidth=1.0, zorder=75, transform=ccrs.PlateCarree());
    #END FOR geom

使用函数 convert_to_mag：

def convert_to_mag(geom, time4mag, alt4mag): #based on fabulously thorough code at https://gis.stackexchange.com/a/291293
    import aacgmv2 #install with: pip install aacgmv2 [need buildtools what a pain]
    # from math import isnan as isnan
    if geom.is_empty:
        return geom
    #END IF
    if geom.has_z:
        def convert_to_mag_doer(coords, time4mag, _):
            for long, lat, alt4mag in coords:
                [lat_mag, long_mag, alt_mag] = aacgmv2.convert_latlon(lat, long, alt4mag, time4mag, method_code='G2A'); #converts from geographic to geomagnetic (AACGMv2)
                yield (long_mag, lat_mag, alt_mag)
            #END FOR long, lat, alt_mag
        #END DEF
    else:
        def convert_to_mag_doer(coords, time4mag, alt4mag):
            for long, lat in coords:
                [lat_mag, long_mag, _] = aacgmv2.convert_latlon_arr(lat, long, alt4mag, time4mag, method_code='G2A'); #converts from geographic to geomagnetic (AACGMv2)
                # tryCntr = 0;
                # while( (isnan(long_mag.item()) | isnan(lat_mag.item())) & (tryCntr < 5) ):
                #     #recalc at higher altitude
                #     [long_mag, lat_mag, _] = aacgmv2.convert_latlon_arr(lat, long, alt4mag+200, time4mag, method_code='G2A'); #converts from geographic to geomagnetic (AACGMv2)
                #     tryCntr += 1 #increment
                # #END IF
                yield (long_mag.item(), lat_mag.item())
            #END FOR long, lat
        #END DEF
    #END IF
    # Process coordinates from each supported geometry type
    if geom.type in ('Point', 'LineString', 'LinearRing'):
        return type(geom)(list(convert_to_mag_doer(geom.coords, time4mag, alt4mag)))
    elif geom.type == 'Polygon':
        ring = geom.exterior
        shell = type(ring)(list(convert_to_mag_doer(ring.coords, time4mag, alt4mag)));
        holes = list(geom.interiors);
        for pos, ring in enumerate(holes):
            holes[pos] = type(ring)(list(convert_to_mag_doer(ring.coords, time4mag, alt4mag)));
        #END FOR pos, ring
        return type(geom)(shell, holes)
    elif geom.type.startswith('Multi') or geom.type == 'GeometryCollection':
        # Recursive call
        return type(geom)([convert_to_mag(part, time4mag, alt4mag) for part in geom.geoms])
    else:
        raise ValueError('Type %r not recognized' % geom.type)
    #END IF
#END DEF