palaeomag a95 plots

[amer7632]

implemented plotting a95 circles from pmagpy (https://github.com/PmagPy/PmagPy)

-three functions added to gplately/plot.py -shoot (where the a95 circle is calculated) -equi (where the a95 is plotted) -plot pole (master function where equi is called)

-temporary data folder added ('./pmag_test_data') to store Merdith et al. (2021) reconstruction files and a .gpml of palaeomagnetic data from Torsvik et al. (2012) -temporary notebook added to show the plotting function ('./pmag_demo_plot')

[michaelchin]

Is this PR ready to merge? @brmather do you know anything about this PR?

[brmather]

@michaelchin not yet - I still need to go through the PR and change some things.

[brmather]

Hey @amer7632,

This might be a silly question, but can you achieve the same using matplotlib.patches.Circle? I've made a small example:

from matplotlib import patches as mpatches

proj = ccrs.Mollweide(central_longitude = 0)

# Set up a GeoAxis plot
fig = plt.figure(figsize=(16,12), dpi=100)
ax = fig.add_subplot(111, projection=proj)
ax.gridlines(color='k',linestyle='--', xlocs=np.arange(-180,180,15), ylocs=np.arange(-90,90,15))
plt.title('Subduction zones and mid-ocean ridges reconstructed to %i Ma' % (time))

# Plot shapefile features, subduction zones and MOR boundaries at 50 Ma
gplot.time = 50 # Ma
gplot.plot_continent_ocean_boundaries(ax, color='b', alpha=0.05)
gplot.plot_continents(ax, facecolor='palegoldenrod', alpha=0.2)
gplot.plot_coastlines(ax, color='DarkKhaki')
gplot.plot_ridges_and_transforms(ax, color='red')
gplot.plot_trenches(ax, color='k')
gplot.plot_subduction_teeth(ax, color='k')
ax.set_global()

# define pole lon/lat/a95
pole_lon = 45
pole_lat = 45
pole_a95 = 15

# shoot and equi
gplately.plot.plot_pole(ax, pole_lon, pole_lat, pole_a95)

# adding a patch
ax.add_patch(mpatches.Circle(xy=[pole_lon, pole_lat], radius=pole_a95,
                             color='red', alpha=0.3, transform=ccrs.PlateCarree(), zorder=30))

I define a point at 45,45 degrees (lon,lat) and an uncertainty bracket of 15 degrees (the gridlines are spaced every 15 degrees). The black circle Is plotted using the plot_pole function and the red circle using the circle patch. I expected the black circle to have a radius of 15 degrees but it is more like 15 height, 20 width. Both are consistent if there is very little map distortion e.g. 0,0 (lon,lat).

Do you know why this is the case? If matplotlib.patches.Circle does the job then we should use that because it is much more efficient and uses fewer lines of code.

[brmather]

Here is another map with more comparisons at different lat,lon coordinates.

[brmather]

And another using a Mercator projection.

[amer7632]

can you do a plot with orthographic projection to check?, possibly with some data that would overlap the pole as well!

I think that the original method calculates a series of points the centre following the azimith at X distance away (X being defined by the a95 uncertainty). Following a few of the links it seems that originally it was designed for basemap and then just adopted into cartopy. It is... quite chunky and a bit hard to follow!

All we need for the a95 is to draw a circle of X radius around a point so I think what you've done is fine (and way simpler)

a few things:

1/ can it be plotted unfilled (with different colours passed)?

2/ in some cases a95s are given as an ellipse with a width/height, is it possible also to pass in as a secondary option for mpatches.ellipse?

otherwise I think it's good, the final check will be to compare visually to GPlates display but I can do that afterwards

[brmather]

OK - I think it now makes sense why the great circle distances were different. The map below has an additional set of circles (in purple) which calculates the distance using an Orthographic projection. It's based on this Stack Overflow answer. This plots an error envelope similar to the shoot / equi method. The same will work for a mpatches.ellipse setting a width and height that is double the radius.

I'll go ahead and implement that now in this PR. I'm thinking something like gplot.plot_pole(ax, lon, lat, a95, **kwargs) where you can specify facecolors, edgecolors, linestyles, etc. using keyword-argument pairs like other matplotlib patches.

[amer7632]

ah, this looks perfect now, and yes, agreed, the suggested format should be perfect. might make a note in the documentation of the function that it plots pole lat/long, and that the site lat/long (i.e. the 'dot' in the continent) has to be plotted separately (which should be pretty simple with just a basic plt.scatter() call

[michaelchin]

@amer7632 @brmather

Just out of curiosity, what's the plan for this PR? Do you want this change in the next major release?

[brmather]

This has now been reimplemented in #249. Close this PR.