Error in plotting coastlines in magnetic co-ordinates

39 Views Asked by At

SO I was trying to plot coastlines in both geographical and geomagnetic co-ordinates. I took the latitude and longitude values of the geographical co-ordinate system, and converted the co-ordinates to geomagnetic using AACGMv2 as mentioned here (How to use cartopy to plot data in geomagnetic coordinates? and here - Plotting Coastlines and Other Features in Magnetic Coordinates with Cartopy). However, they had wrongly used aacgmv2.convert_latlon_arr() to calculate magnetic co-ordinates instead of aacgmv2.get_aacgm_coord_arr() so I fixed it and tried from my end.

But I am getting some "weird" polygon shape inside my map when plotting in geomagnetic co-ordinates (as shown by the red arrow). Does anyone know why?

weird polygon shape inside the map as marked by a red arrow

Fso0L

This is the code which I had tried to get the plots:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import datetime as dt
import numpy as np
import aacgmv2

fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(1,2,1, projection=ccrs.SouthPolarStereo())  # for geographic lat, lon plot
ax1 = fig.add_subplot(1,2,2, projection=ccrs.SouthPolarStereo()) # for geomagnetic lat, lon plot


cc_geo = cfeature.NaturalEarthFeature('physical', 'coastline', '110m', color='xkcd:black', zorder=75); #get a feature to mess with
geom_mag = []; #prep a list

time4mag = dt.datetime(2013,11,7,2,30)
alt4mag = 300

for geom in cc_geo.geometries():
    geom_mag.append(convert_to_mag(geom, time4mag, alt4mag));


cc_mag = cfeature.ShapelyFeature(geom_mag, ccrs.SouthPolarStereo(), color='xkcd:black',zorder=75);


for geo in cc_geo.geometries():
    ax.plot(*geo.coords.xy, color='xkcd:black', linewidth=1.0, zorder=75, transform=ccrs.PlateCarree());
    ax.set_extent([-180, 180, -90, -60], ccrs.PlateCarree())
    ax.set_title('Coastlines in geographical lat and long')
    ax.title.set_size(20)

for geom in cc_mag.geometries():
    ax1.plot(*geom.coords.xy, color='xkcd:black', linewidth=1.0, zorder=75, transform=ccrs.PlateCarree());
    ax1.set_extent([-180, 180, -90, -60], ccrs.PlateCarree())
    ax1.set_title('Coastlines in geomagnetic lat and long')
    ax1.title.set_size(20)

where the function convert_to_mag() is given as:

def convert_to_mag(geom, time4mag, alt4mag): 
    import aacgmv2 
    
    if geom.is_empty:
        return geom

   
    def convert_to_mag_doer(coords, time4mag, alt4mag):
        for long, lat in coords:
            [mlat, mlon, mlt] = aacgmv2.get_aacgm_coord_arr(
                lat, long, alt4mag, time4mag
            ) # convert geographic lat, lon to geomagnetic lat, lon, magnetic local time
       
            yield (mlon.item(), mlat.item())

    
    # 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)

The coastline plots finally produced are not looking satisfactory though.

0

There are 0 best solutions below