# -*- coding: utf-8 -*-
"""
Created on Thu May 30 18:28:24 2013
@author: jpeacock-pr
"""
#==============================================================================
import os
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.ticker import MultipleLocator
from mtpy.analysis.zinvariants import Zinvariants
import mtpy.imaging.mtplottools as mtpl
#==============================================================================
[docs]class PlotStrike(object):
"""
PlotStrike will plot the strike estimated from the invariants, phase tensor
and the tipper in either a rose diagram of xy plot
plots the strike angle as determined by phase tensor azimuth (Caldwell et
al. [2004]) and invariants of the impedance tensor (Weaver et al. [2003]).
The data is split into decades where the histogram for each is plotted in
the form of a rose diagram with a range of 0 to 180 degrees.
Where 0 is North and 90 is East. The median angle of the period band is
set in polar diagram. The top row is the strike estimated from
the invariants of the impedance tensor. The bottom row is the azimuth
estimated from the phase tensor. If tipper is 'y' then the 3rd row is the
strike determined from the tipper, which is orthogonal to the induction
arrow direction.
Arguments
----------
**fn_list** : list of strings
full paths to .edi files to plot
**z_object** : class mtpy.core.z.Z
object of mtpy.core.z. If this is input be sure the
attribute z.freq is filled. *default* is None
**mt_object** : class mtpy.imaging.mtplot.MTplot
object of mtpy.imaging.mtplot.MTplot
*default* is None
**fignum** : int
figure number to be plotted. *Default* is 1
**fs** : float
font size for labels of plotting. *Default* is 10
**dpi** : int
dots-per-inch resolution of figure, 300 is needed for
publications. *Default* is 300
**rot_z** : float
angle of rotation clockwise positive. *Default* is 0
**period_tolerance** : float
Tolerance level to match periods from different edi files.
*Default* is 0.05
**text_dict** : dictionary
*'pad' : float
padding of the angle label at the bottom of each
polar diagram. *Default* is 1.65
*'size' : float
font size
**plot_range** : [ 'data' | (period_min,period_max) ]
period range to estimate the strike angle. Options are:
* *'data'* for estimating the strike for all periods
in the data.
* (pmin,pmax) for period min and period max, input as
(log10(pmin),log10(pmax))
**plot_type** : [ 1 | 2 ]
-*1* to plot individual decades in one plot
-*2* to plot all period ranges into one polar diagram
for each strike angle estimation
**plot_tipper** : [ 'y' | 'n' ]
-*'y'* to plot the tipper strike
-*'n'* to not plot tipper strike
**pt_error_floor** : float
Maximum error in degrees that is allowed to estimate strike.
*Default* is None allowing all estimates to be used.
**fold** : [ True | False ]
*True to plot only from 0 to 180
*False to plot from 0 to 360
:Example: ::
>>> import os
>>> import mtpy.imaging.mtplot as mtplot
>>> edipath = r"/home/EDIFiles"
>>> edilist = [os.path.join(edipath,edi) for edi in os.listdir(edipath)
>>> ... if edi.find('.edi')>0]
>>> #---plot rose plots in decades with tipper and an error floor on pt
>>> strike = mtplot.plot_strike(fn_list=edilist, plot_type=1,\
pt_error_floor=5)
>>> #---plot all decades into one rose plot for each estimation---
>>> strike.plot_type = 2
>>> strike.redraw_plot()
>>> #---save the plot---
>>> strike.save_plot(r"/home/Figures")
'Figure saved to /home/Figures/StrikeAnalysis_.pdf'
Attributes
-----------
-axhinv matplotlib.axes instance for invariant strike
-axhpt matplotlib.axes instance for phase tensor strike
-axhtip matplotlib.axes instance for tipper strike
-barinv matplotlib.axes.bar instance for invariant strike
-barpt matplotlib.axes.bar instance for pt strike
-bartr matplotlib.axes.bar instance for tipper strike
-bin_width width of histogram bins in degrees
-fig matplotlib.figure instance of plot
-fig_dpi dots-per-inch resolution of figure
-fig_num number of figure being plotted
-fig_size size of figure in inches
-fold boolean to fold angles to range from [0,180] or
[0,360]
-font_size font size of axes tick labels
-mt_list list of mtplot.MTplot instances containing all
the important information for each station
-period_tolerance tolerance to look for periods being plotted
-plot_range range of periods to plot
-plot_tipper string to tell program to plot induction arrows
-plot_type string to tell program how to plot strike angles
-plot_yn plot strike on instance creation
-pt_error_floor error floor to plot phase tensor strike, anything
above this error will not be plotted
-text_pad padding between text and rose diagram
-text_size font size of text labeling the mode of the histogram
-title_dict title dictionary
Methods
--------
-plot plots the pseudo section
-redraw_plot on call redraws the plot from scratch
-save_figure saves figure to a file of given format
-update_plot updates the plot while still active
-export_pt_params_to_file writes parameters of the phase tensor and tipper
to text files.
"""
def __init__(self, **kwargs):
fn_list = kwargs.pop('fn_list', None)
z_object_list = kwargs.pop('z_object_list', None)
tipper_object_list = kwargs.pop('tipper_object_list', None)
mt_object_list = kwargs.pop('mt_object_list', None)
#------Set attributes of the class-----------------
#--> get the inputs into a list of mt objects
self.mt_list = mtpl.get_mtlist(fn_list=fn_list,
z_object_list=z_object_list,
tipper_object_list=tipper_object_list,
mt_object_list=mt_object_list)
self._rot_z = kwargs.pop('rot_z', 0)
if isinstance(self._rot_z, float) or isinstance(self._rot_z, int):
self._rot_z = np.array([self._rot_z] * len(self.mt_list))
# if the rotation angle is an array for rotation of different
# freq than repeat that rotation array to the len(mt_list)
elif isinstance(self._rot_z, np.ndarray):
if self._rot_z.shape[0] != len(self.mt_list):
self._rot_z = np.repeat(self._rot_z, len(self.mt_list))
else:
pass
#--> set plot properties
self.fig_num = kwargs.pop('fig_num', 1)
self.fig_dpi = kwargs.pop('fig_dpi', 300)
self.fig_size = kwargs.pop('fig_size', [7, 5])
self.plot_num = kwargs.pop('plot_num', 1)
self.plot_type = kwargs.pop('plot_type', 2)
self.plot_title = kwargs.pop('plot_title', None)
self.plot_range = kwargs.pop('plot_range', 'data')
self.plot_tipper = kwargs.pop('plot_tipper', 'n')
self.period_tolerance = kwargs.pop('period_tolerance', .05)
self.pt_error_floor = kwargs.pop('pt_error_floor', None)
self.fold = kwargs.pop('fold', True)
self.show_ptphimin = kwargs.pop('show_ptphimin',False)
self.bin_width = kwargs.pop('bin_width', 5)
self.font_size = kwargs.pop('font_size', 7)
text_dict = kwargs.pop('text_dict', {})
try:
self.text_pad = text_dict['pad']
except KeyError:
self.text_pad = 0.6
try:
self.text_size = text_dict['size']
except KeyError:
self.text_size = self.font_size
# make a dictionary for plotting titles
self.title_dict = {}
self.title_dict[-5] = '10$^{-5}$--10$^{-4}$s'
self.title_dict[-4] = '10$^{-4}$--10$^{-3}$s'
self.title_dict[-3] = '10$^{-3}$--10$^{-2}$s'
self.title_dict[-2] = '10$^{-2}$--10$^{-1}$s'
self.title_dict[-1] = '10$^{-1}$--10$^{0}$s'
self.title_dict[0] = '10$^{0}$--10$^{1}$s'
self.title_dict[1] = '10$^{1}$--10$^{2}$s'
self.title_dict[2] = '10$^{2}$--10$^{3}$s'
self.title_dict[3] = '10$^{3}$--10$^{4}$s'
self.title_dict[4] = '10$^{4}$--10$^{5}$s'
self.title_dict[5] = '10$^{5}$--10$^{6}$s'
self.plot_yn = kwargs.pop('plot_yn', 'y')
if self.plot_yn == 'y':
self.plot()
#---need to rotate data on setting rotz
def _set_rot_z(self, rot_z):
"""
need to rotate data when setting z
"""
# if rotation angle is an int or float make an array the length of
# mt_list for plotting purposes
if isinstance(rot_z, float) or isinstance(rot_z, int):
self._rot_z = np.array([rot_z] * len(self.mt_list))
# if the rotation angle is an array for rotation of different
# freq than repeat that rotation array to the len(mt_list)
elif isinstance(rot_z, np.ndarray):
if rot_z.shape[0] != len(self.mt_list):
self._rot_z = np.repeat(rot_z, len(self.mt_list))
else:
pass
for ii, mt in enumerate(self.mt_list):
mt.rot_z = self._rot_z[ii]
def _get_rot_z(self):
return self._rot_z
rot_z = property(fget=_get_rot_z, fset=_set_rot_z,
doc="""rotation angle(s)""")
def plot(self, show=True):
plt.rcParams['font.size'] = self.font_size
plt.rcParams['figure.subplot.left'] = .07
plt.rcParams['figure.subplot.right'] = .98
plt.rcParams['figure.subplot.bottom'] = .09
plt.rcParams['figure.subplot.top'] = .90
plt.rcParams['figure.subplot.wspace'] = .2
plt.rcParams['figure.subplot.hspace'] = .4
bw = self.bin_width
if self.fold == True:
histrange = (-180, 180)
elif self.fold == False:
histrange = (0, 360)
# set empty lists that will hold dictionaries with keys as the period
invlist = []
ptlist = []
tiprlist = []
# initialize some parameters
nc = len(self.mt_list)
nt = 0
kk = 0
for dd, mt in enumerate(self.mt_list):
#-----------get strike angle from invariants-----------------------
zinv = Zinvariants(mt.Z)
# add 90 degrees because invariants assume 0 is north, but plotting
# assumes that 90 is north and measures clockwise, thus the negative
# because the strike angle from invariants is measured
# counter-clockwise
zs = 90 - zinv.strike
# fold so the angle goes from 0 to 180
if self.fold == True:
# for plotting put the NW angles into the SE quadrant
zs[np.where(zs > 90)] -= 180
zs[np.where(zs < -90)] += 180
# leave as the total unit circle 0 to 360
elif self.fold == False:
zs[np.where(zs < 0)] += 360
#zs = 360-zs
# make a dictionary of strikes with keys as period
mdictinv = dict([(ff, jj) for ff, jj in zip(mt.period, zs)])
invlist.append(mdictinv)
#------------get strike from phase tensor strike angle-------------
pt = mt.pt
az = 90 - pt.azimuth
az_err = pt.azimuth_err
az[pt.phimax == 0] = np.nan
# need to add 90 because pt assumes 0 is north and
# negative because measures clockwise.
# put an error max on the estimation of strike angle
if self.pt_error_floor:
az[np.where(az_err > self.pt_error_floor)] = 0.0
# fold so the angle goes from 0 to 180
if self.fold:
az[np.where(az > 90)] -= 180
az[np.where(az < -90)] += 180
# leave as the total unit circle 0 to 360
elif self.fold == False:
az[np.where(az < 0)] += 360
# make a dictionary of strikes with keys as period
mdictpt = dict([(ff, jj) for ff, jj in zip(mt.period, az)])
ptlist.append(mdictpt)
#-----------get tipper strike------------------------------------
tip = mt.Tipper
if tip.tipper is None:
tip.tipper = np.zeros((len(mt.period), 1, 2),
dtype='complex')
tip.compute_components()
# needs to be negative because measures clockwise
tipr = -tip.angle_real
tipr[np.where(tipr == 180.)] = 0.0
# fold so the angle goes from 0 to 180
if self.fold == True:
tipr[np.where(tipr > 90)] -= 180
tipr[np.where(tipr < -90)] += 180
# leave as the total unit circle 0 to 360
elif self.fold == False:
tipr[np.where(tipr < 0)] += 360
tipr[np.where(tipr == 360.0)] = 0.0
# make a dictionary of strikes with keys as period
tiprdict = dict([(ff, jj) for ff, jj in zip(mt.period, tipr)])
tiprlist.append(tiprdict)
# get unique periods from all data and min and max values
plist = np.unique(np.hstack([mm.keys() for mm in ptlist]))
minper = np.min(plist)
maxper = np.max(plist)
nt = len(plist)
pdict = dict([(ii, jj) for jj, ii in enumerate(plist)])
# make empty arrays to put data into for easy manipulation
medinv = np.zeros((nt, nc))
medpt = np.zeros((nt, nc))
medtipr = np.zeros((nt, nc))
self._plist = plist
# put data into arrays
for ii, mm in enumerate(invlist):
mperiod = list(mm.keys())
for jj, mp in enumerate(mperiod):
for kk in list(pdict.keys()):
if mp > kk * (1 - self.period_tolerance) and \
mp < kk * (1 + self.period_tolerance):
ll = pdict[kk]
medinv[ll, ii] = invlist[ii][mp]
medpt[ll, ii] = ptlist[ii][mp]
medtipr[ll, ii] = tiprlist[ii][mp]
else:
pass
# medpt = np.hstack([medpt,medpt+180])
# make the arrays local variables
self._medinv = medinv
self._medpt = medpt
self._medtp = medtipr
#-----Plot Histograms of the strike angles-----------------------------
if self.plot_range == 'data':
brange = np.arange(np.floor(np.log10(minper)),
np.ceil(np.log10(maxper)), 1)
else:
brange = np.arange(np.floor(self.plot_range[0]),
np.ceil(self.plot_range[1]), 1)
self._brange = brange
# font dictionary
fd = {'size': self.font_size, 'weight': 'normal'}
#------------------plot indivdual decades------------------------------
if self.plot_type == 1:
# plot specs
plt.rcParams['figure.subplot.hspace'] = .3
plt.rcParams['figure.subplot.wspace'] = .3
self.fig = plt.figure(self.fig_num, dpi=self.fig_dpi)
plt.clf()
nb = len(brange)
for jj, bb in enumerate(brange, 1):
# make subplots for invariants and phase tensor azimuths
if self.plot_tipper == 'n':
self.axhinv = self.fig.add_subplot(2, nb, jj, polar=True)
self.axhpt = self.fig.add_subplot(
2, nb, jj + nb, polar=True)
axlist = [self.axhinv, self.axhpt]
if self.plot_tipper == 'y':
self.axhinv = self.fig.add_subplot(3, nb, jj, polar=True)
self.axhpt = self.fig.add_subplot(
3, nb, jj + nb, polar=True)
self.axhtip = self.fig.add_subplot(3, nb, jj + 2 * nb,
polar=True)
axlist = [self.axhinv, self.axhpt, self.axhtip]
# make a list of indicies for each decades
binlist = []
for ii, ff in enumerate(plist):
if ff > 10**bb and ff < 10**(bb + 1):
binlist.append(ii)
# extract just the subset for each decade
hh = medinv[binlist, :]
gg = medpt[binlist, :]
ptplotdata = gg.flatten()
ptplotdata = ptplotdata[np.nonzero(ptplotdata)]
ptplotdata = ptplotdata[np.isfinite(ptplotdata)]
if self.plot_tipper == 'y':
tr = medtipr[binlist, :]
# compute the historgram for the tipper strike
trhist = np.histogram(tr[np.nonzero(tr)].flatten(),
bins=int(360/bw),
range=histrange)
# make a bar graph with each bar being width of bw degrees
bartr = self.axhtip.bar((trhist[1][:-1]) * np.pi / 180,
trhist[0],
width=bw * np.pi / 180)
# set color of the bars according to the number in that bin
# tipper goes from dark blue (low) to light blue (high)
for cc, bar in enumerate(bartr):
try:
fc = float(trhist[0][cc]) / trhist[0].max() * .9
bar.set_facecolor((0, 1 - fc / 2, fc))
except ZeroDivisionError:
pass
plotinvdata = hh.flatten()
plotinvdata = plotinvdata[np.nonzero(plotinvdata)]
plotinvdata = plotinvdata[np.isfinite(plotinvdata)]
invhist = np.histogram(plotinvdata,
bins= int(360/bw),
range=histrange)
pthist = np.histogram(ptplotdata,
bins=int(360/bw),
range=histrange)
# plot the histograms
invhistlocs = np.mean([invhist[1][:-1],invhist[1][1:]],axis=0)
pthistlocs = np.mean([pthist[1][:-1],pthist[1][1:]],axis=0)
self.barinv = self.axhinv.bar(invhistlocs * np.pi/180,
invhist[0],
width=bw * np.pi / 180)
self.barpt = self.axhpt.bar(pthistlocs * np.pi / 180,
pthist[0],
width=bw * np.pi / 180)
# set the color of the bars according to the number in that bin
# invariants go from purple (low) to red (high)
for cc, bar in enumerate(self.barinv):
try:
fc = float(invhist[0][cc]) / invhist[0].max() * .8
bar.set_facecolor((fc, 0, 1 - fc))
except ZeroDivisionError:
pass
# pt goes from green (low) to orange (high)
for cc, bar in enumerate(self.barpt):
try:
fc = float(pthist[0][cc]) / pthist[0].max() * .8
bar.set_facecolor((fc, 1 - fc, 0))
except ZeroDivisionError:
pass
# make axis look correct with N to the top at 90.
for aa, axh in enumerate(axlist):
# set multiple locator to be every 15 degrees
axh.xaxis.set_major_locator(
MultipleLocator(30 * np.pi / 180))
# set labels on the correct axis
axh.xaxis.set_ticks(np.arange(0.,360.,90.)*np.pi/180.)
axh.xaxis.set_ticklabels(['E','N','W','S'])
axh.tick_params(pad=-5)
# make a light grid
axh.grid(alpha=.25)
# properties for the invariants
if aa == 0:
# limits need to be rotate 90 counter clockwise because
# we already rotated by 90 degrees so the range is
# from -90 to 270 with -90 being east
axh.set_xlim(-90 * np.pi / 180, 270 * np.pi / 180)
# label the plot with the mode value of strike
# need to subtract 90 again because the histogram is
# for ploting 0 east, 90 north measuring
# counter-clockwise
invmode = 90 - invhistlocs[np.where(
invhist[0] == invhist[0].max())[0][0]]
if invmode < 0:
invmode += 360
axh.text(np.pi, axh.get_ylim()[1] * self.text_pad,
'{0:.1f}$^o$'.format(invmode),
horizontalalignment='center',
verticalalignment='baseline',
fontdict={'size': self.text_size},
bbox={'facecolor': (.9, 0, .1), 'alpha': .25})
# print out the statistics of the strike angles
invmedian = 90 - np.median(plotinvdata)
if invmedian < 0:
invmedian += 360
invmean = 90 - np.mean(plotinvdata)
if invmean < 0:
invmean += 360
print('-----Period Range {0:.3g} to {1:.3g} (s)-----'.format(10**bb,
10**(bb + 1)))
print(' *Z-Invariants: median={0:.1f} mode={1:.1f} mean={2:.1f}'.format(
invmedian,
invmode,
invmode))
#--> set title of subplot
axh.set_title(self.title_dict[bb], fontdict=fd,
bbox={'facecolor': 'white', 'alpha': .25})
#--> set the title offset
axh.titleOffsetTrans._t = (0, .1)
# set pt axes properties
elif aa == 1:
# limits go from -180 to 180 as that is how the angle
# is calculated
axh.set_xlim(-180 * np.pi / 180, 180 * np.pi / 180)
# label plot with the mode of the strike angle
ptmode = 90 - pthist[1][np.where(
pthist[0] == pthist[0].max())[0][0]]
if ptmode < 0:
ptmode += 360
ptmedian = 90 - np.median(ptplotdata)
if ptmedian < 0:
ptmedian += 360
ptmean = 90 - np.mean(ptplotdata)
if ptmean < 0:
ptmean += 360
axh.text(np.pi, axh.get_ylim()[1] * self.text_pad,
'{0:.1f}$^o$'.format(ptmode),
horizontalalignment='center',
verticalalignment='baseline',
fontdict={'size': self.text_size},
bbox={'facecolor': (.9, .9, 0), 'alpha': .25})
# print out the results for the strike angles
print(' *PT Strike: median={0:.1f} mode={1:.1f} mean={2:.1f}'.format(
ptmedian,
ptmode,
ptmean))
if self.plot_tipper != 'y':
print('\n')
if nb > 5:
axh.set_title(self.title_dict[bb], fontdict=fd,
bbox={'facecolor': 'white', 'alpha': .25})
# set tipper axes properties
elif aa == 2:
# limits go from -180 to 180
axh.set_xlim(-180 * np.pi / 180, 180 * np.pi / 180)
# label plot with mode
tpmode = 90 - trhist[1][np.where(
trhist[0] == trhist[0].max())[0][0]]
if tpmode < 0:
tpmode += 360.
tpmedian = 90 - np.median(tr[np.nonzero(tr)])
if tpmedian < 0:
tpmedian += 360
tpmean = 90 - np.mean(tr[np.nonzero(tr)])
if tpmean < 0:
tpmean += 360
axh.text(np.pi, axh.get_ylim()[1] * self.text_pad,
'{0:.1f}$^o$'.format(tpmode),
horizontalalignment='center',
verticalalignment='baseline',
fontdict={'size': self.text_size},
bbox={'facecolor': (0, .1, .9), 'alpha': .25})
# print out statistics for strike angle
print(' *Tipper Strike: median={0:.1f} mode={1:.1f} mean={2:.1f}'.format(
tpmedian,
tpmode,
tpmode))
print('\n')
if nb > 5:
axh.set_title(self.title_dict[bb], fontdict=fd,
bbox={'facecolor': 'white', 'alpha': .25})
# set plot labels
if jj == 1:
if aa == 0:
axh.set_ylabel('Strike (Z)', fontdict=fd,
labelpad=self.font_size,
bbox={'facecolor': (.9, 0, .1),
'alpha': 0.25})
elif aa == 1:
axh.set_ylabel('PT Azimuth', fontdict=fd,
labelpad=self.font_size,
bbox={'facecolor': (.9, .9, 0),
'alpha': .25})
elif aa == 2:
axh.set_ylabel('Tipper Strike', fd,
labelpad=self.font_size,
bbox={'facecolor': (0, .1, .9),
'alpha': 0.25})
plt.setp(axh.yaxis.get_ticklabels(), visible=False)
print('Note: North is assumed to be 0 and the strike angle is measured' +\
'clockwise positive.')
if show:
plt.show()
#------------------Plot strike angles for all period ranges------------
elif self.plot_type == 2:
# plot specs
plt.rcParams['figure.subplot.left'] = .07
plt.rcParams['figure.subplot.right'] = .98
plt.rcParams['figure.subplot.bottom'] = .100
plt.rcParams['figure.subplot.top'] = .88
plt.rcParams['figure.subplot.hspace'] = .3
plt.rcParams['figure.subplot.wspace'] = .2
self.fig = plt.figure(self.fig_num,
self.fig_size,
dpi=self.fig_dpi)
plt.clf()
# make subplots for invariants and phase tensor azimuths
if self.plot_tipper == 'n':
self.axhinv = self.fig.add_subplot(1, 2, 1, projection='polar') # polar=True)
self.axhpt = self.fig.add_subplot(1, 2, 2, projection='polar') # polar=True)
axlist = [self.axhinv, self.axhpt]
else:
self.axhinv = self.fig.add_subplot(1, 3, 1, polar=True)
self.axhpt = self.fig.add_subplot(1, 3, 2, polar=True)
self.axhtip = self.fig.add_subplot(1, 3, 3, polar=True)
axlist = [self.axhinv, self.axhpt, self.axhtip]
# make a list of indicies for each decades
binlist = [pdict[ff] for ff in plist
if ff > 10**brange.min() and ff < 10**brange.max()]
# extract just the subset for each decade
hh = medinv[binlist, :]
gg = medpt[binlist, :]
plotinvdata = hh.flatten()
plotinvdata = plotinvdata[np.nonzero(plotinvdata)]
plotinvdata = plotinvdata[np.isfinite(plotinvdata)]
# estimate the histogram for the decade for invariants and pt
invhist = np.histogram(hh[np.nonzero(hh)].flatten(),
bins=int(360/bw),
range=histrange)
ptplotdata = gg.flatten()
ptplotdata = ptplotdata[np.nonzero(ptplotdata)]
ptplotdata = ptplotdata[np.isfinite(ptplotdata)]
if not self.fold:
ptplotdata = np.hstack([ptplotdata,ptplotdata+180])
if self.show_ptphimin:
if self.fold:
# we are working in matplotlib coordinates (counterclockwise from east)
# therefore range of values is -90 to 90 for a range of 0 to 180 east of north
# if 0 to 90 need to subtract 90 to get angle of phimin
# if -90 to 0 need to add 90 to get angle of phimin.
zero_to_90 = np.where(np.all([ptplotdata > 0, ptplotdata <= 90],axis=0))
m90_to_zero = np.where(np.all([ptplotdata > -90, ptplotdata < 0],axis=0))
ptplotdata = np.hstack([ptplotdata,
ptplotdata[zero_to_90] - 90,
ptplotdata[m90_to_zero] + 90])
else:
ptplotdata = np.hstack([ptplotdata,(ptplotdata+90)%360])
pthist = np.histogram(ptplotdata,
bins=int(360/bw),
range=histrange
)
# centre points for the bins
invhistlocs = np.mean([invhist[1][:-1],invhist[1][1:]],axis=0)
pthistlocs = np.mean([pthist[1][:-1],pthist[1][1:]],axis=0)
# plot the histograms
self.barinv = self.axhinv.bar(invhistlocs * np.pi / 180,
invhist[0],
width=bw * np.pi / 180)
self.barpt = self.axhpt.bar(pthistlocs * np.pi / 180,
pthist[0],
width=bw * np.pi / 180)
# set color of invariants from purple (low) to red (high count)
for cc, bar in enumerate(self.barinv):
fc = float(invhist[0][cc]) / invhist[0].max() * .8
bar.set_facecolor((fc, 0, 1 - fc))
# set color of pt from green (low) to orange (high count)
for cc, bar in enumerate(self.barpt):
fc = float(pthist[0][cc]) / pthist[0].max() * .8
bar.set_facecolor((fc, 1 - fc, 0))
# plot tipper if desired
if self.plot_tipper == 'y':
tr = self._medtp[binlist, :]
trhist = np.histogram(tr[np.nonzero(tr)].flatten(),
bins= int(360/bw),
range=histrange
)
self.bartr = self.axhtip.bar((trhist[1][:-1]) * np.pi / 180,
trhist[0],
width=bw * np.pi / 180)
# set tipper color from dark blue (low) to light blue (high)
for cc, bar in enumerate(self.bartr):
try:
fc = float(trhist[0][cc]) / trhist[0].max() * .9
bar.set_facecolor((0, 1 - fc / 2, fc))
except ZeroDivisionError:
pass
# make axis look correct with N to the top at 90.
for aa, axh in enumerate(axlist):
# set major ticks to be every 30 degrees
axh.xaxis.set_major_locator(MultipleLocator(2 * np.pi / 12))
# set a light grid
axh.grid(alpha=0.25) # works in 2.0.2 not 2.1.0
# set tick labels to be invisible
plt.setp(axh.yaxis.get_ticklabels(), visible=False) # works in 2.0.2 not 2.1.0
# place the correct label at the cardinal directions
axh.xaxis.set_ticks(np.arange(0.,360.,90.)*np.pi/180.)
axh.xaxis.set_ticklabels(['E','N','W','S'])
axh.tick_params(pad=-5)
# set invariant axes propertiesz
if aa == 0:
axh.set_ylim(0, invhist[0].max())
invmode = 90 - invhist[1][np.where(
invhist[0] == invhist[0].max())[0][0]]
if invmode < 0:
invmode += 360
axh.text(170 * np.pi / 180, axh.get_ylim()[1] * .65,
'{0:.1f}$^o$'.format(invmode),
horizontalalignment='center',
verticalalignment='baseline',
fontdict={'size': self.font_size},
bbox={'facecolor': (.9, 0, .1), 'alpha': .25})
# print out the statistics of the strike angles
invmedian = 90 - np.median(plotinvdata)
if invmedian < 0:
invmedian += 360
invmean = 90 - np.mean(plotinvdata)
if invmean < 0:
invmean += 360
print('-----Period Range {0:.3g} to {1:.3g} (s)-----'.format(10**brange[0],
10**brange[-1]))
print(' *Z-Invariants: median={0:.1f} mode={1:.1f} mean={2:.1f}'.format(
invmedian,
invmode,
invmode))
axh.set_title('Strike (Z)', fontdict=fd,
bbox={'facecolor': (.9, 0, .1), 'alpha': 0.25})
# set pt axes properties
elif aa == 1:
axh.set_ylim(0, pthist[0].max())
ptmode = 90 - pthist[1][np.where(
pthist[0] == pthist[0].max())[0][0]]
if ptmode < 0:
ptmode += 360
if ptmode < 0:
ptmode += 360
ptmedian = 90 - np.median(ptplotdata)
if ptmedian < 0:
ptmedian += 360
ptmean = 90 - np.mean(ptplotdata)
if ptmean < 0:
ptmean += 90
axh.text(170 * np.pi / 180, axh.get_ylim()[1] * .65,
'{0:.1f}$^o$'.format(ptmode),
horizontalalignment='center',
verticalalignment='baseline',
fontdict={'size': self.text_size},
bbox={'facecolor': (.9, .9, 0), 'alpha': 0.25})
# print results of strike analysis for pt
print(' *PT Strike: median={0:.1f} mode={1:.1f} mean={2:.1f}'.format(
ptmedian,
ptmode,
ptmean))
if self.plot_tipper != 'y':
print('\n')
axh.set_title('PT Azimuth', fontdict=fd,
bbox={'facecolor': (.9, .9, 0), 'alpha': 0.25})
# set tipper axes properties
elif aa == 2:
axh.set_ylim(0, trhist[0].max())
tpmode = 90 - trhist[1][np.where(
trhist[0] == trhist[0].max())[0][0]]
if tpmode < 0:
tpmode += 360
tpmedian = 90 - np.median(tr[np.nonzero(tr)])
if tpmedian < 0:
tpmedian += 360
tpmean = 90 - np.mean(tr[np.nonzero(tr)])
if tpmean < 0:
tpmean += 360
axh.text(170 * np.pi / 180, axh.get_ylim()[1] * .65,
'{0:.1f}$^o$'.format(tpmode),
horizontalalignment='center',
verticalalignment='baseline',
fontdict={'size': self.text_size},
bbox={'facecolor': (0, .1, .9), 'alpha': 0.25})
print(' *Tipper Stike: median={0:.1f} mode={1:.1f} mean={2:.1f}\n'.format(
tpmedian,
tpmode,
tpmean))
axh.set_title('Tipper Strike', fontdict=fd,
bbox={'facecolor': (0, .1, .9), 'alpha': 0.25})
# move title up a little to make room for labels
axh.titleOffsetTrans._t = (0, .15)
# remind the user what the assumptions of the strike angle are
print('Note: North is assumed to be 0 and the strike angle is ' +\
'measured clockwise positive.')
if show:
plt.show()
return ptplotdata,plotinvdata
[docs] def save_plot(self, save_fn, file_format='pdf',
orientation='portrait', fig_dpi=None, close_plot='y'):
"""
save_plot will save the figure to save_fn.
Arguments
-----------
**save_fn** : string
full path to save figure to, can be input as
- directory path -> the directory path to save to
in which the file will be saved as
save_fn/station_name_ResPhase.file_format
- full path -> file will be save to the given
path. If you use this option then the format
will be assumed to be provided by the path
**file_format** : [ pdf | eps | jpg | png | svg ]
file type of saved figure pdf,svg,eps...
**orientation** : [ landscape | portrait ]
orientation in which the file will be saved
*default* is portrait
**fig_dpi** : int
The resolution in dots-per-inch the file will be
saved. If None then the dpi will be that at
which the figure was made. I don't think that
it can be larger than dpi of the figure.
**close_plot** : [ y | n ]
-'y' will close the plot after saving.
-'n' will leave plot open
Examples
---------
:Example: ::
>>> # to save plot as jpg
>>> import mtpy.imaging.mtplottools as mtplot
>>> p1 = mtplot.PlotPhaseTensorMaps(edilist,freqspot=10)
>>> p1.save_plot(r'/home/MT', file_format='jpg')
'Figure saved to /home/MT/PTMaps/PTmap_phimin_10Hz.jpg'
"""
# get rid of . in file format as it will be added later
if file_format is not None:
file_format = file_format.replace('.','')
if fig_dpi is None:
fig_dpi = self.fig_dpi
if os.path.isdir(save_fn) == False:
file_format = save_fn[-3:]
self.fig.savefig(save_fn, dpi=fig_dpi, format=file_format,
orientation=orientation)
# plt.clf()
# plt.close(self.fig)
else:
if not os.path.exists(save_fn):
os.mkdir(save_fn)
save_fn = os.path.join(save_fn, 'StrikeAnalysis.' + file_format)
self.fig.savefig(save_fn, dpi=fig_dpi, format=file_format,
orientation=orientation)
if close_plot == 'y':
plt.clf()
plt.close(self.fig)
else:
pass
self.fig_fn = save_fn
print('Saved figure to: ' + self.fig_fn)
[docs] def update_plot(self):
"""
update any parameters that where changed using the built-in draw from
canvas.
Use this if you change an of the .fig or axes properties
:Example: ::
>>> # to change the grid lines to only be on the major ticks
>>> import mtpy.imaging.mtplottools as mtplot
>>> p1 = mtplot.PlotResPhase(r'/home/MT/mt01.edi')
>>> [ax.grid(True, which='major') for ax in [p1.axr,p1.axp]]
>>> p1.update_plot()
"""
self.fig.canvas.draw()
[docs] def redraw_plot(self):
"""
use this function if you updated some attributes and want to re-plot.
:Example: ::
>>> # change the color and marker of the xy components
>>> import mtpy.imaging.mtplottools as mtplot
>>> p1 = mtplot.PlotResPhase(r'/home/MT/mt01.edi')
>>> p1.xy_color = (.5,.5,.9)
>>> p1.xy_marker = '*'
>>> p1.redraw_plot()
"""
plt.close(self.fig)
self.plot()
def __str__(self):
"""
rewrite the string builtin to give a useful message
"""
return "Plots phase tensor maps for one freq"
[docs] def writeTextFiles(self, save_path=None):
"""
Saves the strike information as a text file.
"""
# check to see if the strikes have been calculated
try:
self.bin_width
except AttributeError:
self.plot()
# get the path to save the file to
if save_path is None:
try:
svpath = os.path.dirname(self.mt_list[0].fn)
except TypeError:
raise IOError('Need to input save_path, could not find path')
else:
svpath = save_path
# set
if self.fold == True:
histrange = (-180, 180)
elif self.fold == False:
histrange = (0, 360)
# set the bin width
bw = self.bin_width
slistinv = [['station']]
slistpt = [['station']]
slisttip = [['station']]
# calculate the strikes for the different period bands
for jj, bb in enumerate(self._brange):
tstr = self.title_dict[bb].replace('$', '')
tstr = tstr.replace('{', '').replace('}', '').replace('^', 'e')
tstr = tstr.replace('s', '(s)')
slistinv[0].append(tstr)
slistpt[0].append(tstr)
slisttip[0].append(tstr)
# calculate the strike for the different period bands per station
for kk, mt in enumerate(self.mt_list, 1):
if jj == 0:
slistinv.append([mt.station])
slistpt.append([mt.station])
slisttip.append([mt.station])
zinv = mt.Z.invariants
pt = mt.pt
tp = mt.Tipper
bnlist = []
for nn, per in enumerate(mt.period):
if per > 10**bb and per < 10**(bb + 1):
bnlist.append(nn)
#---> strike from invariants
zs = 90 - zinv.strike[bnlist]
# fold so the angle goes from 0 to 180
if self.fold == True:
# for plotting put the NW angles into the SE quadrant
zs[np.where(zs > 90)] = zs[np.where(zs > 90)] - 180
zs[np.where(zs < -90)] = zs[np.where(zs < -90)] + 180
# leave as the total unit circle 0 to 360
elif self.fold == False:
pass
zshist = np.histogram(zs[np.nonzero(zs)].flatten(),
bins=int(360/bw),
range=histrange)
#==============================================================
# For putting the values into a useful text file
# need to subtract 90 from the values to put them into
# coordinates where north is 0 and east is 90, which is
# different from plotting where east in the plotting function
# is equal to 0 and north is 90, measuring counter-clockwise
#==============================================================
#==> compute mean
invmean = 90 - zs.mean()
if invmean < 0:
invmean += 360
invmed = 90 - np.median(zs)
#==> compute median
if invmed < 0:
invmed += 360
#==> compute mode
invmode = 90 - zshist[1][np.where(
zshist[0] == zshist[0].max())[0][0]]
if invmode < 0:
invmode += 360
#==> append to list
slistinv[kk].append((invmean,
invmed,
invmode))
#---> strike from phase tensor
az = pt.azimuth[0][bnlist]
# fold so the angle goes from 0 to 180
if self.fold == True:
az[np.where(az > 90)] = az[np.where(az > 90)] - 180
az[np.where(az < -90)] = az[np.where(az < -90)] + 180
# leave as the total unit circle 0 to 360
elif self.fold == False:
az[np.where(az < 0)] = az[np.where(az < 0)] + 360
# == > compute mean
ptmean1 = 90 - az.mean()
if ptmean1 < 0:
ptmean1 += 360
# == > compute median
ptmed1 = 90 - np.median(az)
if ptmed1 < 0:
ptmed1 += 360
# == > compute mode
azhist = np.histogram(az[np.nonzero(az)].flatten(),
bins=int(360/bw),
range=histrange)
ptmode1 = 90 - azhist[1][np.where(
azhist[0] == azhist[0].max())[0][0]]
if ptmode1 < 0:
ptmode1 += 360
slistpt[kk].append((ptmean1,
ptmed1,
ptmode1))
#---> strike from tipper
# needs to be negative because measures clockwise
if tp._Tipper.tipper is None:
tp._Tipper.tipper = np.zeros((len(mt.period), 1, 2),
dtype='complex')
tp.compute_components()
tipr = -tp.angle_real[bnlist]
# fold so the angle goes from 0 to 180
if self.fold == True:
tipr[np.where(tipr > 90)] = tipr[np.where(tipr > 90)] - 180
tipr[np.where(tipr < -90)
] = tipr[np.where(tipr < -90)] + 180
# leave as the total unit circle 0 to 360
elif self.fold == False:
tipr[np.where(tipr < 0)] = tipr[np.where(tipr < 0)] + 360
tphist = np.histogram(tipr[np.nonzero(tipr)].flatten(),
bins=int(360/bw),
range=histrange)
#==> compute mean
tpmean1 = 90 - tipr.mean()
if tpmean1 < 0:
tpmean1 += 360
#==> compute median
tpmed1 = 90 - np.median(tipr)
if tpmed1 < 0:
tpmed1 += 360
#==> compute mode
tpmode1 = 90 - tphist[1][np.where(
tphist[0] == tphist[0].max())[0][0]]
if tpmode1 < 0:
tpmode1 += 360
#--> append statistics to list
slisttip[kk].append((tpmean1,
tpmed1,
tpmode1))
# make a list of indicies for each decades
binlist = []
for ii, ff in enumerate(self._plist):
if ff > 10**bb and ff < 10**(bb + 1):
binlist.append(ii)
# extract just the subset for each decade
hh = self._medinv[binlist, :]
gg = self._medpt[binlist, :]
ptplotdata = gg[np.nonzero(gg)].flatten()
ptplotdata = ptplotdata[np.isfinite(ptplotdata)]
tr = self._medtp[binlist, :]
# estimate the histogram for the decade for invariants and pt
invhist = np.histogram(hh[np.nonzero(hh)].flatten(),
bins=int(360/bw),
range=histrange)
pthist = np.histogram(ptplotdata,
bins=int(360/bw),
range=histrange)
trhist = np.histogram(tr[np.nonzero(tr)].flatten(),
bins=int(360/bw),
range=histrange)
#--> include the row for mean, median and mode for each parameter
if jj == 0:
slistinv.append(['mean'])
slistinv.append(['median'])
slistinv.append(['mode'])
slistpt.append(['mean'])
slistpt.append(['median'])
slistpt.append(['mode'])
slisttip.append(['mean'])
slisttip.append(['median'])
slisttip.append(['mode'])
#--> compute mean, median and mode for invariants
# == > mean
imean = 90 - np.mean(hh[np.nonzero(hh)])
if imean < 0:
imean += 360
# == > median
imed = 90 - np.median(hh[np.nonzero(hh)])
if imed < 0:
imed += 360
# == > mode
imode = 90 - invhist[1][np.where(
invhist[0] == invhist[0].max())[0][0]]
if imode < 0:
imode += 360
#--> add them to the list of estimates
slistinv[kk + 1].append(imean)
slistinv[kk + 2].append(imed)
slistinv[kk + 3].append(imode)
#--> compute pt statistics
# == > mean
ptmean = 90 - np.mean(ptplotdata)
if ptmean < 0:
ptmean = np.mean(ptplotdata)
# == > median
ptmed = 90 - np.median(ptplotdata)
if ptmed < 0:
ptmed += 360
# == > mode
ptmode = 90 - pthist[1][np.where(
pthist[0] == pthist[0].max())[0][0]]
if ptmode < 0:
ptmode += 360
#--> add the statistics to the parameter list
slistpt[kk + 1].append(ptmean)
slistpt[kk + 2].append(ptmed)
slistpt[kk + 3].append(ptmode)
#--> compute tipper statistics
# == > mean
tpmean = 90 - np.mean(tipr[np.nonzero(tipr)])
if tpmean < 0:
tpmean += 360
# == > median
tpmed = 90 - np.median(tipr[np.nonzero(tipr)])
if tpmed < 0:
tpmed += 360
# == > mode
tpmode = 90 - trhist[1][np.where(
trhist[0] == trhist[0].max())[0][0]]
if tpmode < 0:
tpmode += 360
#--> add the statistics to parameter list
slisttip[kk + 1].append(tpmean)
slisttip[kk + 2].append(tpmed)
slisttip[kk + 3].append(tpmode)
invfid = file(os.path.join(svpath, 'Strike.invariants'), 'w')
ptfid = file(os.path.join(svpath, 'Strike.pt'), 'w')
tpfid = file(os.path.join(svpath, 'Strike.tipper'), 'w')
#---> write strike from the invariants
# == > mean
invfid.write('-' * 20 + 'MEAN' + '-' * 20 + '\n')
for ii, l1 in enumerate(slistinv):
for jj, l2 in enumerate(l1):
if ii == 0:
invfid.write('{0:^16}'.format(l2))
else:
if jj == 0:
invfid.write('{0:>16}'.format(l2 + ' ' * 6))
else:
try:
invfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2[0])))
except IndexError:
invfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2)))
invfid.write('\n')
# == > median
invfid.write('-' * 20 + 'MEDIAN' + '-' * 20 + '\n')
for ii, l1 in enumerate(slistinv):
for jj, l2 in enumerate(l1):
if ii == 0:
invfid.write('{0:^16}'.format(l2))
else:
if jj == 0:
invfid.write('{0:>16}'.format(l2 + ' ' * 6))
else:
try:
invfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2[1])))
except IndexError:
invfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2)))
invfid.write('\n')
# == > mode
invfid.write('-' * 20 + 'MODE' + '-' * 20 + '\n')
for ii, l1 in enumerate(slistinv):
for jj, l2 in enumerate(l1):
if ii == 0:
invfid.write('{0:^16}'.format(l2))
else:
if jj == 0:
invfid.write('{0:>16}'.format(l2 + ' ' * 6))
else:
try:
invfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2[2])))
except IndexError:
invfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2)))
invfid.write('\n')
invfid.close()
#---> write the phase tensor text files
ptfid.write('-' * 20 + 'MEAN' + '-' * 20 + '\n')
for ii, l1 in enumerate(slistpt):
for jj, l2 in enumerate(l1):
if ii == 0:
ptfid.write('{0:^16}'.format(l2))
else:
if jj == 0:
ptfid.write('{0:>16}'.format(l2 + ' ' * 6))
else:
try:
ptfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2[0])))
except IndexError:
ptfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2)))
ptfid.write('\n')
ptfid.write('-' * 20 + 'MEDIAN' + '-' * 20 + '\n')
for ii, l1 in enumerate(slistpt):
for jj, l2 in enumerate(l1):
if ii == 0:
ptfid.write('{0:^16}'.format(l2))
else:
if jj == 0:
ptfid.write('{0:>16}'.format(l2 + ' ' * 6))
else:
try:
ptfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2[1])))
except IndexError:
ptfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2)))
ptfid.write('\n')
ptfid.write('-' * 20 + 'MODE' + '-' * 20 + '\n')
for ii, l1 in enumerate(slistpt):
for jj, l2 in enumerate(l1):
if ii == 0:
ptfid.write('{0:^16}'.format(l2))
else:
if jj == 0:
ptfid.write('{0:>16}'.format(l2 + ' ' * 6))
else:
try:
ptfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2[2])))
except IndexError:
ptfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2)))
ptfid.write('\n')
ptfid.close()
#---> write the tipper text files
tpfid.write('-' * 20 + 'MEAN' + '-' * 20 + '\n')
for ii, l1 in enumerate(slisttip):
for jj, l2 in enumerate(l1):
if ii == 0:
tpfid.write('{0:^16}'.format(l2))
else:
if jj == 0:
tpfid.write('{0:>16}'.format(l2 + ' ' * 6))
else:
try:
tpfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2[0])))
except IndexError:
tpfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2)))
tpfid.write('\n')
tpfid.write('-' * 20 + 'MEDIAN' + '-' * 20 + '\n')
for ii, l1 in enumerate(slisttip):
for jj, l2 in enumerate(l1):
if ii == 0:
tpfid.write('{0:^16}'.format(l2))
else:
if jj == 0:
tpfid.write('{0:>16}'.format(l2 + ' ' * 6))
else:
try:
tpfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2[1])))
except IndexError:
tpfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2)))
tpfid.write('\n')
tpfid.write('-' * 20 + 'MODE' + '-' * 20 + '\n')
for ii, l1 in enumerate(slisttip):
for jj, l2 in enumerate(l1):
if ii == 0:
tpfid.write('{0:^16}'.format(l2))
else:
if jj == 0:
tpfid.write('{0:>16}'.format(l2 + ' ' * 6))
else:
try:
tpfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2[2])))
except IndexError:
tpfid.write('{0:^16}'.format(
'{0: .2f}'.format(l2)))
tpfid.write('\n')
tpfid.close()