How to use cartopy to plot data in geomagnetic coordinates?

Question

I am a python beginner and I use cartopy to point my data in geographic coordinates. However, I want to plot my data in geomagnetic coordinates. I saw that there is a rotated pole option in cartopy. What I am not sure of is: if I rotate the poles using the magnetic North pole's location (latitude and longitude), would that give me correct results? Because, the north and south magnetic poles are not symmetric to each other (because of the solar wind's influence), but geographic poles are symmetric. As far as I understand, the rotated pole option in cartopy shifts the geographic poles according to the given new pole location symmetrically.

Furthermore, the magnetic north pole's latitude and longitude on Earth changes from year to year. So, I guess, the magnetic north pole's location should be updated separately for each date we are interested in. These values can be obtained from several places, so this is not a problem.

Coming to my actual question: what should I change in my code below to plot my data correctly in geomagnetic coordinates? Would the rotated pole option in cartopy work in this case? Any suggestions about other ways of handling this kind of problem is also appreciated.

This is what I use for plotting in geographical coordinates:

# 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

While the rotated pole idea is interesting, accurate geomagnetic coordinates aren't just a different pole location since the Earth's magnetic field is funky (and always changing). IGRF or AACGMv2 (more precise) can be used to convert to geomagnetic coordinates from geomagnetic coordinates and then there's no worry about unsymmetrical poles with the rotated pole method - the IGRF/AACGMv2 conversion takes care of that.

Data (not Cartopy-shapely objects) can be converted directly with [lat_mag, long_mag, _] = aacgmv2.convert_latlon_arr(lat, long, alt4mag, time4mag, method_code='G2A'); #converts from geographic to geomagnetic (AACGMv2)

Per https://github.com/SciTools/cartopy/issues/1945 the following code converts Cartopy coastlines from geographic coordinates to geomagnetic coordinates and plots them (I'm not quite sure how to deal with ocean/land/lake coordinates but it should be a similar methodology):

time4mag is a datetime object localized to the UTC timezone

alt4mag is an integer or float that is an altitude in km

Note: for something other than geo->mag coords, replace aacgmv2.convert_latlon_arr with any other function that adjusts lat/long coords for the specific application

    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

with the function 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