In this example we will plot phase tensors from ModEM files. This example is a bit more complex than previous examples, as, unlike the previous examples, the basemap plotting functionality is not contained within MTPy. This has the benefit that it makes it easier to customise the plot. But it may mean it takes a bit longer to become familiar with the functionality.
The first step is to import the required modules needed. We only have one import from MTPy - PlotPTMpas. Then there is some standard matplotlib functionality and importantly the basemap module which creates coastlines and the nice borders.
In [1]:
from mtpy.modeling.modem import PlotPTMaps
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib as mpl
from mpl_toolkits.basemap import Basemap
from shapely.geometry import Polygon
from descartes import PolygonPatch
import numpy as np
The next step is to create a function that will draw an inset map showing the survey boundaries on Australia.
In [2]:
# function to draw a bounding box
def drawBBox( minLon, minLat, maxLon, maxLat, bm, **kwargs):
bblons = np.array([minLon, maxLon, maxLon, minLon, minLon])
bblats = np.array([minLat, minLat, maxLat, maxLat, minLat])
x, y = bm( bblons, bblats )
xy = zip(x,y)
poly = Polygon(xy)
bm.ax.add_patch(PolygonPatch(poly, **kwargs))
We now need to define our file paths for the response and data files
In [3]:
# define paths
data_fn = r'C:/mtpywin/mtpy/examples/model_files/ModEM_2/ModEM_Data.dat'
resp_fn = r'C:/mtpywin/mtpy/examples/model_files/ModEM_2/Modular_MPI_NLCG_004.dat'
# define extents
minLat = -22.5
maxLat = -18.5
minLon = 135.5
maxLon = 140.5
# define period index to plot
periodIdx = 10
# position of inset axes (bottom,left,width,height)
inset_ax_position = [0.6,0.2,0.3,0.2]
We can now create the plot!
In [8]:
# read in ModEM data to phase tensor object
plotPTM = PlotPTMaps(data_fn=data_fn, resp_fn=resp_fn)
# make a figure
fig, ax = plt.subplots(figsize=(10,10))
# make a basemap
m = Basemap(resolution='c', # c, l, i, h, f or None
ax=ax,
projection='merc',
lat_0=-20.5, lon_0=132, # central lat/lon for projection
llcrnrlon=minLon, llcrnrlat=minLat, urcrnrlon=maxLon, urcrnrlat=maxLat)
# draw lat-lon grids
m.drawparallels(np.linspace(minLat, maxLat, 5), labels=[1,1,0,0], linewidth=0.1)
m.drawmeridians(np.linspace(minLon, maxLon, 5), labels=[0,0,1,1], linewidth=0.1)
## draw shaded topographic relief map
## (need to install pil first before this will work - conda install pil)
# m.shadedrelief()
# plot inset map ==================================================================
insetAx = fig.add_axes(inset_ax_position)
mInset = Basemap(resolution='c', # c, l, i, h, f or None
ax=insetAx,
projection='merc',
lat_0=-20, lon_0=132,
llcrnrlon=110, llcrnrlat=-40, urcrnrlon=155, urcrnrlat=-10)
mInset.fillcontinents(color='lightgray')
mInset.drawstates(color="grey")
drawBBox(minLon, minLat, maxLon, maxLat, mInset, fill='True', facecolor='k')
# plot phase tensors =============================================================
# fetch attribute to color phase tensor ellipses with
cmapAttrib = plotPTM.get_period_attributes(periodIdx, 'phimin', ptarray='data')
sm = cm.ScalarMappable(norm=mpl.colors.Normalize(vmin=np.min(cmapAttrib),
vmax=np.max(cmapAttrib)),
cmap='gist_heat')
# extract color values from colormap
cvals = sm.cmap(sm.norm(cmapAttrib))
plotPTM.plot_on_axes(ax, m, periodIdx=periodIdx, ptarray='data',
cvals=cvals, ellipse_size_factor=2e4,
edgecolor='k')
# show colormap
cbax = fig.add_axes([0.25,0,0.525,.025])
cbar = mpl.colorbar.ColorbarBase(cbax, cmap=sm.cmap,
norm = sm.norm,
orientation='horizontal')
cbar.set_label('Phase tensors coloured by phimin')
plt.savefig('/tmp/a.pdf', dpi=300)