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
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.
