使用 Cartopy 在磁坐标中绘制海岸线和其他特征

Question

我想使用磁坐标而不是地理坐标来绘制海岸线形状。我想用于坐标转换的例程是AACGMv2，它将接收地理纬度和经度的向量，并返回磁经纬度的向量。

我希望截取 cartopy.feature.COASTLINE 的坐标，运行转换，然后绘制结果。但是，我还没有弄清楚如何做到这一点。

有人知道最好的方法（或任何方法）吗？

Answer 1

每https://github.com/SciTools/cartopy/issues/1945:

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