From a8bf78afc9db0dd6b491a08095b71ed7f48454a8 Mon Sep 17 00:00:00 2001 From: JP Date: Sun, 3 May 2020 21:48:35 -0700 Subject: [PATCH 01/20] set rotation matrix to be formulated as counterclockwise positive --- mtpy/utils/calculator.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/mtpy/utils/calculator.py b/mtpy/utils/calculator.py index 585bdaddb..1286656da 100644 --- a/mtpy/utils/calculator.py +++ b/mtpy/utils/calculator.py @@ -493,8 +493,8 @@ def rotatematrix_incl_errors(inmatrix, angle, inmatrix_err = None) : # counter clockwise, I cannot find a good reason for this except that # when you plot the strike and phase tensors the look correct with this # formulation. - rotmat = np.array([[ cphi, sphi], [-sphi, cphi]]) - # rotmat = np.array([[ cphi, -sphi], [sphi, cphi]]) + # rotmat = np.array([[ cphi, sphi], [-sphi, cphi]]) + rotmat = np.array([[ cphi, -sphi], [sphi, cphi]]) rotated_matrix = np.dot(np.dot(rotmat, inmatrix), np.linalg.inv(rotmat)) errmat = None @@ -546,8 +546,8 @@ def rotatevector_incl_errors(invector, angle, invector_err = None): # counter clockwise, I cannot find a good reason for this except that # when you plot the strike and phase tensors the look correct with this # formulation. - rotmat = np.array([[ cphi, sphi],[-sphi, cphi]]) - # rotmat = np.array([[ cphi, -sphi],[sphi, cphi]]) + # rotmat = np.array([[ cphi, sphi],[-sphi, cphi]]) + rotmat = np.array([[ cphi, -sphi],[sphi, cphi]]) if invector.shape == (1, 2): rotated_vector = np.dot( invector, np.linalg.inv(rotmat) ) From 054eec1b4b7acb5bc7ea804205736f4ce2f50dad Mon Sep 17 00:00:00 2001 From: JP Date: Sun, 3 May 2020 23:05:38 -0700 Subject: [PATCH 02/20] set rotation matrix to be formulated as counterclockwise positive --- mtpy/analysis/pt.py | 2 +- mtpy/core/z.py | 9 --------- mtpy/utils/calculator.py | 8 ++++---- 3 files changed, 5 insertions(+), 14 deletions(-) diff --git a/mtpy/analysis/pt.py b/mtpy/analysis/pt.py index 3d62b8f05..7a4f99106 100644 --- a/mtpy/analysis/pt.py +++ b/mtpy/analysis/pt.py @@ -571,7 +571,7 @@ def azimuth(self): @property def azimuth_err(self): if self.pt_err is not None: - az_err = np.sqrt(self.alpha+self.beta) + az_err = np.sqrt(self.alpha + self.beta) else: az_err = None diff --git a/mtpy/core/z.py b/mtpy/core/z.py index 60208c1da..43e53c3ca 100644 --- a/mtpy/core/z.py +++ b/mtpy/core/z.py @@ -111,21 +111,12 @@ def compute_resistivity_phase(self, z_array=None, z_err_array=None, for idx_f in range(self.freq.size): for ii in range(2): for jj in range(2): -# r_err, phi_err = MTcc.z_error2r_phi_error( -# np.real(self._z[idx_f, ii, jj]), -# self._z_err[idx_f, ii, jj], -# np.imag(self._z[idx_f, ii, jj]), -# self._z_err[idx_f, ii, jj]) - r_err, phi_err = MTcc.z_error2r_phi_error( self._z[idx_f, ii, jj].real, self._z[idx_f, ii, jj].imag, self._z_err[idx_f, ii, jj]) self._resistivity_err[idx_f, ii, jj] = \ self._resistivity[idx_f, ii, jj] * r_err -# self._resistivity_err[idx_f, ii, jj] = \ -# 0.4 * np.abs(self._z[idx_f, ii, jj]) / \ -# self.freq[idx_f] * r_err self._phase_err[idx_f, ii, jj] = phi_err def set_res_phase(self, res_array, phase_array, freq, res_err_array=None, diff --git a/mtpy/utils/calculator.py b/mtpy/utils/calculator.py index 1286656da..585bdaddb 100644 --- a/mtpy/utils/calculator.py +++ b/mtpy/utils/calculator.py @@ -493,8 +493,8 @@ def rotatematrix_incl_errors(inmatrix, angle, inmatrix_err = None) : # counter clockwise, I cannot find a good reason for this except that # when you plot the strike and phase tensors the look correct with this # formulation. - # rotmat = np.array([[ cphi, sphi], [-sphi, cphi]]) - rotmat = np.array([[ cphi, -sphi], [sphi, cphi]]) + rotmat = np.array([[ cphi, sphi], [-sphi, cphi]]) + # rotmat = np.array([[ cphi, -sphi], [sphi, cphi]]) rotated_matrix = np.dot(np.dot(rotmat, inmatrix), np.linalg.inv(rotmat)) errmat = None @@ -546,8 +546,8 @@ def rotatevector_incl_errors(invector, angle, invector_err = None): # counter clockwise, I cannot find a good reason for this except that # when you plot the strike and phase tensors the look correct with this # formulation. - # rotmat = np.array([[ cphi, sphi],[-sphi, cphi]]) - rotmat = np.array([[ cphi, -sphi],[sphi, cphi]]) + rotmat = np.array([[ cphi, sphi],[-sphi, cphi]]) + # rotmat = np.array([[ cphi, -sphi],[sphi, cphi]]) if invector.shape == (1, 2): rotated_vector = np.dot( invector, np.linalg.inv(rotmat) ) From e4d726514222a152c8e12849b140f3cdaef71f03 Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 4 May 2020 11:55:03 -0700 Subject: [PATCH 03/20] fixing plotting phase tensor ellipses in correct coordinate system --- mtpy/imaging/phase_tensor_pseudosection.py | 58 ++++------------------ mtpy/imaging/plotstrike.py | 28 ++++++----- mtpy/utils/calculator.py | 12 +---- 3 files changed, 25 insertions(+), 73 deletions(-) diff --git a/mtpy/imaging/phase_tensor_pseudosection.py b/mtpy/imaging/phase_tensor_pseudosection.py index 6ac0c25f7..f13cdec15 100644 --- a/mtpy/imaging/phase_tensor_pseudosection.py +++ b/mtpy/imaging/phase_tensor_pseudosection.py @@ -112,16 +112,9 @@ class PlotPhaseTensorPseudoSection(mtpl.PlotSettings): start and stop of station name indicies. ex: for MT01dr station_id=(0,4) will be MT01 - **rotz** : float or np.ndarray + **rotation_angle** : float or np.ndarray angle in degrees to rotate the data clockwise positive. - Can be an array of angle to individually rotate stations or - periods or both. - - If rotating each station by a constant - angle the array needs to have a shape of - (# of stations) - - If rotating by period needs to have shape - # of periods - - If rotating both individually shape=(ns, nf) + a single number for all input files. *Default* is 0 **title** : string @@ -131,7 +124,7 @@ class PlotPhaseTensorPseudoSection(mtpl.PlotSettings): dots per inch of the resolution. *default* is 300 - **fignum** : int + **fig_num** : int figure number. *Default* is 1 **plot_tipper** : [ 'yri' | 'yr' | 'yi' | 'n' ] @@ -394,19 +387,6 @@ def __init__(self, **kwargs): if 'text_offset_y' not in list(self.scale_arrow_dict.keys()): self.scale_arrow_dict['text_offset_y'] = 0. - 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 induction arrow properties ------------------------------- self.plot_tipper = kwargs.pop('plot_tipper', 'n') @@ -434,9 +414,6 @@ def rotation_angle(self, value): only a single value is allowed """ for ii, mt in enumerate(self.mt_list): - # JP: need to set the rotation angle negative for plotting - # I think its because the way polar plots work by measuring - # counter clockwise mt.rotation_angle = value self._rotation_angle = value @@ -532,7 +509,7 @@ def plot(self, show=True): azimuth = pt.azimuth[::-1] # if there are induction arrows, flip them as pt - if self.plot_tipper.find('y') == 0: + if 'y' in self.plot_tipper: tip = mt.Tipper if tip.mag_real is not None: tmr = tip.mag_real[::-1] @@ -596,14 +573,17 @@ def plot(self, show=True): # create an ellipse scaled by phimin and phimax and orient # the ellipse so that north is up and east is right - # need to add 90 to do so instead of subtracting + # Ellipse patch assumes the angle measure counterclockwise + # therefore need to multiply the azimuth by -1 because + # it is measured clockwise positive and subtract 90 to put + # it in a coordinate system of N=0, E=90 ellipd = patches.Ellipse((offset * self.xstretch, np.log10(ff) * self.ystretch), width=ewidth, height=eheight, edgecolor='k', lw=0.5, - angle=azimuth[jj] + 90) + angle=90 - azimuth[jj]) # get ellipse color if cmap.find('seg') > 0: @@ -683,10 +663,6 @@ def plot(self, show=True): self._plot_periodlist = plot_periodlist n = len(plot_periodlist) - # calculate minimum period and maximum period with a stretch factor - # pmin = np.log10(plot_periodlist.min())*self.ystretch - # pmax = np.log10(plot_periodlist.max())*self.ystretch - pmin = int(np.floor(np.log10(plot_periodlist.min()))) pmax = int(np.ceil(np.log10(plot_periodlist.max()))) @@ -698,28 +674,17 @@ def plot(self, show=True): self.offsetlist = offset_sort['offset'] self.stationlist = offset_sort['station'] - # if self.offsetlist[0] > 0: - # print 'rotating' - # print self.stationlist - # self.stationlist = self.stationlist[::-1] # set y-ticklabels if self.tscale == 'period': yticklabels = [mtpl.labeldict[ii] for ii in range(pmin, pmax + 1, 1)] - # yticklabels = ['{0:>4}'.format('{0: .1e}'.format(plot_period_list[ll])) - # for ll in np.arange(0, n, self.ystep)]+\ - # ['{0:>4}'.format('{0: .1e}'.format(plot_period_list[-1]))] self.ax.set_ylabel('Period (s)', fontsize=self.font_size + 2, fontweight='bold') elif self.tscale == 'frequency': - # yticklabels = ['{0:>4}'.format('{0: .1e}'.format(1./plot_period_list[ll])) - # for ll in np.arange(0, n, self.ystep)]+\ - # ['{0:>4}'.format('{0: .1e}'.format(1./plot_period_list[-1]))] - # yticklabels = [mtpl.labeldict[-ii] for ii in range(pmin, pmax + 1, 1)] self.ax.set_ylabel('Frequency (Hz)', @@ -732,15 +697,10 @@ def plot(self, show=True): # --> set tick locations and labels # set y-axis major ticks - # self.ax.yaxis.set_ticks([np.log10(plot_periodlist[ll])*self.ystretch - # for ll in np.arange(0, n, self.ystep)]): self.ax.yaxis.set_ticks(np.arange(pmin * self.ystretch, (pmax + 1) * self.ystretch, self.ystretch)) - # set y-axis minor ticks - # self.ax.yaxis.set_ticks([np.log10(plot_periodlist[ll])*self.ystretch - # for ll in np.arange(0, n, 1)],minor=True) # set y-axis tick labels self.ax.set_yticklabels(yticklabels) diff --git a/mtpy/imaging/plotstrike.py b/mtpy/imaging/plotstrike.py index cd8f5b242..8596176cc 100644 --- a/mtpy/imaging/plotstrike.py +++ b/mtpy/imaging/plotstrike.py @@ -240,15 +240,12 @@ def make_strike_array(self): if mt.period.size > nt: nt = mt.period.size - #-----------get strike angle from invariants----------------------- + #-----------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 + # subtract 90 because polar plot assumes 0 is on the x an 90 is + # on the y + zs = zinv.strike - 90 # fold so the angle goes from 0 to 180 if self.fold == True: @@ -264,9 +261,11 @@ def make_strike_array(self): mdictinv = dict([(ff, jj) for ff, jj in zip(mt.period, zs)]) inv_list.append(mdictinv) - #------------get strike from phase tensor strike angle------------- + #------------get strike from phase tensor strike angle------------ + # subtract 90 because polar plot assumes 0 is on the x an 90 is + # on the y pt = mt.pt - az = 90 - pt.azimuth + az = pt.azimuth - 90 az_err = pt.azimuth_err az[pt.phimax == 0] = np.nan @@ -296,8 +295,9 @@ def make_strike_array(self): dtype='complex') tip.compute_components() - # needs to be negative because measures clockwise - tipr = -tip.angle_real + # # subtract 90 because polar plot assumes 0 is on the x an 90 is + # on the y + tipr = tip.angle_real + 90 tipr[np.where(tipr == 180.)] = 0.0 @@ -316,8 +316,10 @@ def make_strike_array(self): tip_list.append(tiprdict) #--> get min and max period - self.max_per = np.amax([np.max(list(mm.keys())) for mm in inv_list], axis=0) - self.min_per = np.amin([np.min(list(mm.keys())) for mm in pt_list], axis=0) + self.max_per = np.amax([np.max(list(mm.keys())) for mm in inv_list], + axis=0) + self.min_per = np.amin([np.min(list(mm.keys())) for mm in pt_list], + axis=0) # make empty arrays to put data into for easy manipulation medinv = np.zeros((nt, nc)) diff --git a/mtpy/utils/calculator.py b/mtpy/utils/calculator.py index 585bdaddb..6a78037ff 100644 --- a/mtpy/utils/calculator.py +++ b/mtpy/utils/calculator.py @@ -489,12 +489,7 @@ def rotatematrix_incl_errors(inmatrix, angle, inmatrix_err = None) : cphi = np.cos(phi) sphi = np.sin(phi) - # JP: Changed the rotation matrix to be formulated to rotate - # counter clockwise, I cannot find a good reason for this except that - # when you plot the strike and phase tensors the look correct with this - # formulation. rotmat = np.array([[ cphi, sphi], [-sphi, cphi]]) - # rotmat = np.array([[ cphi, -sphi], [sphi, cphi]]) rotated_matrix = np.dot(np.dot(rotmat, inmatrix), np.linalg.inv(rotmat)) errmat = None @@ -541,13 +536,8 @@ def rotatevector_incl_errors(invector, angle, invector_err = None): cphi = np.cos(phi) sphi = np.sin(phi) - - # JP: Changed the rotation matrix to be formulated to rotate - # counter clockwise, I cannot find a good reason for this except that - # when you plot the strike and phase tensors the look correct with this - # formulation. + rotmat = np.array([[ cphi, sphi],[-sphi, cphi]]) - # rotmat = np.array([[ cphi, -sphi],[sphi, cphi]]) if invector.shape == (1, 2): rotated_vector = np.dot( invector, np.linalg.inv(rotmat) ) From abbf1e071e7fa95e50004661a8f6dbbb8872a01a Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 4 May 2020 14:04:44 -0700 Subject: [PATCH 04/20] set rotation angle to be opposite of input in pt.rotate to be consistent with other rotation angles --- mtpy/analysis/pt.py | 22 ++++---- mtpy/imaging/plotpt.py | 120 +++++++++++++++++++---------------------- 2 files changed, 67 insertions(+), 75 deletions(-) diff --git a/mtpy/analysis/pt.py b/mtpy/analysis/pt.py index 7a4f99106..5065e9dcf 100644 --- a/mtpy/analysis/pt.py +++ b/mtpy/analysis/pt.py @@ -741,14 +741,13 @@ def phimax(self): if self.pt is None: return None -# return self._pi2()[0] + self._pi1()[0] return np.degrees(np.arctan(self._pi2()[0] + self._pi1()[0])) @property def phimax_err(self): phimaxerr = None if self.pt_err is not None: - phimaxerr = np.sqrt(self._pi2()[1]**2+self._pi1()[1]**2) + phimaxerr = np.sqrt(self._pi2()[1]**2 + self._pi1()[1]**2) return np.degrees(np.arctan(phimaxerr)) else: @@ -756,14 +755,12 @@ def phimax_err(self): def rotate(self, alpha): """ - Rotate PT array. Change the rotation angles attribute respectively. + Rotate PT array. Change the rotation angles attribute respectively. - Rotation angle must be given in degrees. All angles are referenced to - geographic North, positive in clockwise direction. - (Mathematically negative!) - - In non-rotated state, X refs to North and Y to East direction. + Rotation angle must be given in degrees. All angles are referenced to + North, positive in clockwise direction. (Mathematically negative!) + In non-rotated state, X refs to North and Y to East direction. """ @@ -820,9 +817,14 @@ def rotate(self, alpha): angle = 0. if self.pt_err is not None: - pt_rot[idx_freq], pt_err_rot[idx_freq] = MTcc.rotatematrix_incl_errors(self.pt[idx_freq,:,:], angle, self.pt_err[idx_freq,:,:]) + pt_rot[idx_freq], pt_err_rot[idx_freq] = \ + MTcc.rotatematrix_incl_errors(self.pt[idx_freq,:,:], + -angle, + self.pt_err[idx_freq,:,:]) else: - pt_rot[idx_freq], pt_err_rot = MTcc.rotatematrix_incl_errors(self.pt[idx_freq,:,:], angle) + pt_rot[idx_freq], pt_err_rot = \ + MTcc.rotatematrix_incl_errors(self.pt[idx_freq,:,:], + -angle) # --> set the rotated tensors as the current attributes self._pt = pt_rot diff --git a/mtpy/imaging/plotpt.py b/mtpy/imaging/plotpt.py index 7e7bee9ca..1ee3a2537 100755 --- a/mtpy/imaging/plotpt.py +++ b/mtpy/imaging/plotpt.py @@ -167,19 +167,25 @@ class PlotPhaseTensor(mtpl.MTEllipse): """ def __init__(self, **kwargs): + super().__init__() fn = kwargs.pop('fn', None) z_object = kwargs.pop('z_object', None) mt_object = kwargs.pop('mt_object', None) pt_object = kwargs.pop('pt_object', None) + + self._rotation_angle = 0 #--> get mt object if fn is not None: self._mt = mtpl.MTplot(fn=fn) + self.pt = self._mt.pt elif z_object is not None: self._mt = mtpl.MTplot(z_object=z_object) + self.pt = self._mt.pt elif mt_object is not None: self._mt = mt_object + self.pt = self._mt.pt elif pt_object is not None: self.pt = pt_object self._mt = mtpl.MTplot() @@ -222,11 +228,9 @@ def __init__(self, **kwargs): # read ellipse dict ellipse_dict = kwargs.pop('ellipse_dict', None) if ellipse_dict is None: - self._ellipse_dict = {'size': .25} - else: - self._ellipse_dict = ellipse_dict + ellipse_dict = {'size': .25} - self._read_ellipse_dict() + self._read_ellipse_dict(ellipse_dict) self.ellipse_spacing = kwargs.pop('ellipse_spacing', 1) @@ -234,6 +238,21 @@ def __init__(self, **kwargs): if self.plot_yn == 'y': self.plot() + + #---need to rotate data on setting rotz + @property + def rotation_angle(self): + return self._rotation_angle + + @rotation_angle.setter + def rotation_angle(self, value): + """ + only a single value is allowed + """ + + self.pt.rotation_angle = value + + self._rotation_angle = value def plot(self): """ @@ -256,16 +275,6 @@ def plot(self): self.fig = plt.figure(self.fig_num, self.fig_size, dpi=self.fig_dpi) plt.clf() - # get phase tensor instance - try: - self.pt - self.pt.rotate(self.rot_z) - except AttributeError: - self.pt = self._mt.pt - self.pt.rotate(self.rot_z) - self.zinv = self._mt.Z.invariants - self.zinv.rotate(self.rot_z) - cmap = self.ellipse_cmap ckmin = self.ellipse_range[0] ckmax = self.ellipse_range[1] @@ -281,29 +290,30 @@ def plot(self): # get the properties to color the ellipses by if self.ellipse_colorby == 'phiminang' or \ self.ellipse_colorby == 'phimin': - colorarray = self.pt.phimin[0] + colorarray = self.pt.phimin elif self.ellipse_colorby == 'phidet': - colorarray = np.sqrt(abs(self.pt.det[0])) * (180 / np.pi) + colorarray = np.sqrt(abs(self.pt.det)) * (180 / np.pi) elif self.ellipse_colorby == 'skew' or\ self.ellipse_colorby == 'skew_seg': - colorarray = self.pt.beta[0] + colorarray = self.pt.beta elif self.ellipse_colorby == 'ellipticity': - colorarray = self.pt.ellipticity[0] + colorarray = self.pt.ellipticity + elif self.ellipse_colorby in ['strike', 'azimuth']: + colorarray = (self.pt.azimuth % 180) else: raise NameError(self.ellipse_colorby + ' is not supported') #-------------plotPhaseTensor----------------------------------- self.ax1 = self.fig.add_subplot(3, 1, 1, aspect='equal') - self._mt._period = 1. / self._mt.freq for ii, ff in enumerate(self._mt.period): # make sure the ellipses will be visable - if self.pt.phimax[0][ii] != 0: - eheight = self.pt.phimin[0][ii] / self.pt.phimax[0][ii] *\ + if self.pt.phimax[ii] != 0: + eheight = self.pt.phimin[ii] / self.pt.phimax[ii] *\ self.ellipse_size ewidth = self.ellipse_size @@ -313,22 +323,13 @@ def plot(self): eheight = 0.01 * self.ellipse_size ewidth = 0.01 * self.ellipse_size - # ewidth = self.pt.phimax[0][ii]/self.pt.phimax[0][ii]*\ - # self.ellipse_size - - # alternative scaling - # eheight = self.pt.phimin[0][ii]/max(np.abs(self.pt.phimax[0]))*\ - # self.ellipse_size - # ewidth = self.pt.phimax[0][ii]/max(np.abs(self.pt.phimax[0]))*\ - # self.ellipse_size - # create an ellipse scaled by phimin and phimax and oriented along # the azimuth which is calculated as clockwise but needs to # be plotted counter-clockwise hence the negative sign. ellipd = patches.Ellipse((np.log10(ff) * self.ellipse_spacing, 0), width=ewidth, height=eheight, - angle=90 - self.pt.azimuth[0][ii]) + angle=90 - self.pt.azimuth[ii]) self.ax1.add_patch(ellipd) @@ -352,12 +353,13 @@ def plot(self): xlimits = (np.floor(np.log10(self._mt.period[0])), np.ceil(np.log10(self._mt.period[-1]))) - self.ax1.set_xlim(xlimits) + self.ax1.set_xlim((xlimits[0] * self.ellipse_spacing, + xlimits[1] * self.ellipse_spacing)) tklabels = [] xticks = [] for tk in self.ax1.get_xticks(): try: - tklabels.append(mtpl.labeldict[tk]) + tklabels.append(mtpl.labeldict[tk/self.ellipse_spacing]) xticks.append(tk) except KeyError: pass @@ -415,12 +417,8 @@ def plot(self): #---------------plotStrikeAngle----------------------------------- self.ax2 = self.fig.add_subplot(3, 2, 3) - az = self.pt.azimuth[0] - az_err = self.pt.azimuth[1] - - # put the strike into a coordinate system that goes from -90 to 90 - az[np.where(az > 90)] -= 180 - az[np.where(az < -90)] += 180 + az = self.pt.azimuth % 180 + az_err = self.pt.azimuth_err stlist = [] stlabel = [] @@ -442,12 +440,8 @@ def plot(self): stlist.append(ps2[0]) stlabel.append('PT') try: - strike = self.zinv.strike - strikeerr = np.nan_to_num(self.zinv.strike_err) - # put the strike into a coordinate system that goes from -90 to 90 - strike[np.where(strike > 90)] = strike[np.where(strike > 90)] - 180 - strike[np.where(strike < -90) - ] = strike[np.where(strike < -90)] + 180 + strike = self.zinv.strike % 180 + strike_err = np.nan_to_num(self.zinv.strike_err) # plot invariant strike erxy = self.ax2.errorbar(self._mt.period, @@ -458,7 +452,7 @@ def plot(self): mec=self.strike_inv_color, mew=self.marker_lw, ls='none', - yerr=strikeerr, + yerr=strike_err, ecolor=self.strike_inv_color, capsize=self.marker_size, elinewidth=self.marker_lw) @@ -468,14 +462,10 @@ def plot(self): except AttributeError: print('Could not get z_invariants from pt, input z if desired.') - if self._mt.tipper is not None: + if self._mt.Tipper.tipper is not None: # strike from tipper tp = self._mt.Tipper - s3 = tp.angle_real + 90 - - # fold to go from -90 to 90 - s3[np.where(s3 > 90)] = s3[np.where(s3 > 90)] - 180 - s3[np.where(s3 < -90)] = s3[np.where(s3 < -90)] + 180 + s3 = (tp.angle_real % 180) # plot strike with error bars ps3 = self.ax2.errorbar(self._mt.period, @@ -510,7 +500,7 @@ def plot(self): plt.setp(ltext, fontsize=6) # the legend text fontsize if self.strike_limits is None: - self.strike_limits = (-89.99, 89.99) + self.strike_limits = (-4, 184) self.ax2.set_yscale('linear') self.ax2.set_xscale('log', nonposx='clip') @@ -524,10 +514,10 @@ def plot(self): self.ax2.set_title('Strike', fontdict=font_dictt) #---------plot Min & Max Phase----------------------------------------- - minphi = self.pt.phimin[0] - minphierr = self.pt.phimin[1] - maxphi = self.pt.phimax[0] - maxphierr = self.pt.phimax[1] + minphi = self.pt.phimin + minphierr = self.pt.phimin_err + maxphi = self.pt.phimax + maxphierr = self.pt.phimax_err self.ax3 = self.fig.add_subplot(3, 2, 4, sharex=self.ax2) @@ -558,10 +548,10 @@ def plot(self): elinewidth=self.marker_lw) if self.pt_limits is None: - self.pt_limits = [min([self.pt.phimax[0].min(), - self.pt.phimin[0].min()]) - 3, - max([self.pt.phimax[0].max(), - self.pt.phimin[0].max()]) + 3] + self.pt_limits = [min([self.pt.phimax.min(), + self.pt.phimin.min()]) - 3, + max([self.pt.phimax.max(), + self.pt.phimin.max()]) + 3] if self.pt_limits[0] < -10: self.pt_limits[0] = -9.9 if self.pt_limits[1] > 100: @@ -595,8 +585,8 @@ def plot(self): #-----------------------plotSkew--------------------------------------- - skew = self.pt.beta[0] - skewerr = self.pt.beta[1] + skew = self.pt.beta + skewerr = self.pt.beta_err self.ax4 = self.fig.add_subplot(3, 2, 5, sharex=self.ax2) erskew = self.ax4.errorbar(self._mt.period, @@ -639,8 +629,8 @@ def plot(self): self.ax4.set_title('Skew Angle', fontdict=font_dictt) #----------------------plotEllipticity-------------------------------- - ellipticity = self.pt.ellipticity[0] - ellipticityerr = self.pt.ellipticity[1] + ellipticity = self.pt.ellipticity + ellipticityerr = self.pt.ellipticity_err self.ax5 = self.fig.add_subplot(3, 2, 6, sharex=self.ax2) erskew = self.ax5.errorbar(self._mt.period, From d32d084e322bc7f78d136af3f321018c8ad5691a Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 4 May 2020 14:27:01 -0700 Subject: [PATCH 05/20] set rotation angle to be opposite of input in pt.rotate to be consistent with other rotation angles --- mtpy/core/mt.py | 2 +- mtpy/imaging/plotpt.py | 24 +++++++++++++++++------- 2 files changed, 18 insertions(+), 8 deletions(-) diff --git a/mtpy/core/mt.py b/mtpy/core/mt.py index bde853484..b97b316fa 100644 --- a/mtpy/core/mt.py +++ b/mtpy/core/mt.py @@ -302,7 +302,7 @@ def rotation_angle(self, theta_r): self._Tipper.rotate(theta_r) self.pt.rotate(theta_r) - print(("Rotated Z, Tipper, Phase Tensor and Zinvariants by" + print(("Rotated Z, Tipper, Phase Tensor and Zinvariants by " "{0:.3f} degrees".format(self._rotation_angle))) @Z.setter diff --git a/mtpy/imaging/plotpt.py b/mtpy/imaging/plotpt.py index 1ee3a2537..bc4c62cae 100755 --- a/mtpy/imaging/plotpt.py +++ b/mtpy/imaging/plotpt.py @@ -250,10 +250,20 @@ def rotation_angle(self, value): only a single value is allowed """ - self.pt.rotation_angle = value + self.pt.rotate(value) + self._mt.rotation_angle = value self._rotation_angle = value + def fold_strike(self, strike): + """ + + """ + strike = strike % 180 + strike[np.where(strike > 90)] = 90 - strike[np.where(strike > 90)] + + return strike + def plot(self): """ plots the phase tensor elements @@ -302,7 +312,7 @@ def plot(self): elif self.ellipse_colorby == 'ellipticity': colorarray = self.pt.ellipticity elif self.ellipse_colorby in ['strike', 'azimuth']: - colorarray = (self.pt.azimuth % 180) + colorarray = self.fold_strike(self.pt.azimuth) else: raise NameError(self.ellipse_colorby + ' is not supported') @@ -374,7 +384,7 @@ def plot(self): which='major', color=(.25, .25, .25), lw=.25) - + self.ax1.set_axisbelow(True) plt.setp(self.ax1.get_yticklabels(), visible=False) # add colorbar for PT self.cbax = self.fig.add_axes(self.cb_position) @@ -417,7 +427,7 @@ def plot(self): #---------------plotStrikeAngle----------------------------------- self.ax2 = self.fig.add_subplot(3, 2, 3) - az = self.pt.azimuth % 180 + az = self.fold_strike(self.pt.azimuth) az_err = self.pt.azimuth_err stlist = [] @@ -440,7 +450,7 @@ def plot(self): stlist.append(ps2[0]) stlabel.append('PT') try: - strike = self.zinv.strike % 180 + strike = self.fold_strike(self.zinv.strike) strike_err = np.nan_to_num(self.zinv.strike_err) # plot invariant strike @@ -465,7 +475,7 @@ def plot(self): if self._mt.Tipper.tipper is not None: # strike from tipper tp = self._mt.Tipper - s3 = (tp.angle_real % 180) + s3 = self.fold_strike(tp.angle_real) # plot strike with error bars ps3 = self.ax2.errorbar(self._mt.period, @@ -500,7 +510,7 @@ def plot(self): plt.setp(ltext, fontsize=6) # the legend text fontsize if self.strike_limits is None: - self.strike_limits = (-4, 184) + self.strike_limits = (-94, 94) self.ax2.set_yscale('linear') self.ax2.set_xscale('log', nonposx='clip') From 6f3bcffb632d39a4daa78960ba6d14e5e4202798 Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 11 May 2020 15:59:36 -0700 Subject: [PATCH 06/20] making rotation matrix postitive clockwise and fixing plotting schemes --- mtpy/analysis/pt.py | 4 +- mtpy/core/z.py | 8 +- mtpy/imaging/plotpt.py | 2 +- mtpy/imaging/plotstrike.py | 23 ++- mtpy/utils/calculator.py | 330 +++++++++++++++++++------------------ 5 files changed, 194 insertions(+), 173 deletions(-) diff --git a/mtpy/analysis/pt.py b/mtpy/analysis/pt.py index 5065e9dcf..ca57d2c02 100644 --- a/mtpy/analysis/pt.py +++ b/mtpy/analysis/pt.py @@ -818,12 +818,12 @@ def rotate(self, alpha): if self.pt_err is not None: pt_rot[idx_freq], pt_err_rot[idx_freq] = \ - MTcc.rotatematrix_incl_errors(self.pt[idx_freq,:,:], + MTcc.rotate_matrix_with_errors(self.pt[idx_freq,:,:], -angle, self.pt_err[idx_freq,:,:]) else: pt_rot[idx_freq], pt_err_rot = \ - MTcc.rotatematrix_incl_errors(self.pt[idx_freq,:,:], + MTcc.rotate_matrix_with_errors(self.pt[idx_freq,:,:], -angle) # --> set the rotated tensors as the current attributes diff --git a/mtpy/core/z.py b/mtpy/core/z.py index 43e53c3ca..cd22cc8cb 100644 --- a/mtpy/core/z.py +++ b/mtpy/core/z.py @@ -604,12 +604,12 @@ def rotate(self, alpha): if self.z_err is not None: z_rot[idx_freq], z_err_rot[idx_freq] = \ - MTcc.rotatematrix_incl_errors(self.z[idx_freq, :, :], + MTcc.rotate_matrix_with_errors(self.z[idx_freq, :, :], angle, self.z_err[idx_freq, :, :]) else: z_rot[idx_freq], z_err_rot = \ - MTcc.rotatematrix_incl_errors(self.z[idx_freq, :, :], + MTcc.rotate_matrix_with_errors(self.z[idx_freq, :, :], angle) self.z = z_rot @@ -1512,12 +1512,12 @@ def rotate(self, alpha): if self.tipper_err is not None: tipper_rot[idx_freq], tipper_err_rot[idx_freq] = \ - MTcc.rotatevector_incl_errors(self.tipper[idx_freq, :, :], + MTcc.rotate_vector_with_errors(self.tipper[idx_freq, :, :], angle, self.tipper_err[idx_freq, :, :]) else: tipper_rot[idx_freq], tipper_err_rot = \ - MTcc.rotatevector_incl_errors(self.tipper[idx_freq, :, :], + MTcc.rotate_vector_with_errors(self.tipper[idx_freq, :, :], angle) self.tipper = tipper_rot diff --git a/mtpy/imaging/plotpt.py b/mtpy/imaging/plotpt.py index bc4c62cae..04fe4ade5 100755 --- a/mtpy/imaging/plotpt.py +++ b/mtpy/imaging/plotpt.py @@ -260,7 +260,7 @@ def fold_strike(self, strike): """ strike = strike % 180 - strike[np.where(strike > 90)] = 90 - strike[np.where(strike > 90)] + strike[np.where(strike > 90)] -= 180 return strike diff --git a/mtpy/imaging/plotstrike.py b/mtpy/imaging/plotstrike.py index 8596176cc..f95324239 100644 --- a/mtpy/imaging/plotstrike.py +++ b/mtpy/imaging/plotstrike.py @@ -184,6 +184,7 @@ def __init__(self, **kwargs): self.font_size = 7 self.text_pad = 0.6 self.text_size = self.font_size + self.polar_limits = (np.deg2rad(-180), np.deg2rad(180)) for key, value in kwargs.items(): setattr(self, key, value) @@ -227,6 +228,12 @@ def rotation_angle(self, value): def make_strike_array(self): """ make strike array + + ..note:: Polar plots assume the azimuth is an angle measured + counterclockwise positive from x = 0. Therefore all angles + are calculated as 90 - angle to make them conform to the + polar plot convention. + """ inv_list = [] pt_list = [] @@ -245,7 +252,7 @@ def make_strike_array(self): # subtract 90 because polar plot assumes 0 is on the x an 90 is # on the y - zs = zinv.strike - 90 + zs = 90 - zinv.strike # fold so the angle goes from 0 to 180 if self.fold == True: @@ -265,7 +272,7 @@ def make_strike_array(self): # subtract 90 because polar plot assumes 0 is on the x an 90 is # on the y pt = mt.pt - az = pt.azimuth - 90 + az = 90 - pt.azimuth az_err = pt.azimuth_err az[pt.phimax == 0] = np.nan @@ -297,7 +304,7 @@ def make_strike_array(self): # # subtract 90 because polar plot assumes 0 is on the x an 90 is # on the y - tipr = tip.angle_real + 90 + tipr = 90 - tip.angle_real tipr[np.where(tipr == 180.)] = 0.0 @@ -352,8 +359,7 @@ def make_strike_array(self): self.med_inv = medinv self.med_pt = medpt self.med_tip = medtipr - - + def get_mean(self, st_array): """ get mean value @@ -589,13 +595,14 @@ def plot(self, show=True): # make a light grid axh.grid(alpha=.25, zorder=0) + axh.set_xlim(self.polar_limits) # 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 @@ -632,7 +639,7 @@ def plot(self, show=True): 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) + #axh.set_xlim(-180 * np.pi / 180, 180 * np.pi / 180) pt_median, pt_mode, pt_mean = self.get_stats(plot_pt, pt_hist, @@ -650,7 +657,7 @@ def plot(self, show=True): # set tipper axes properties elif aa == 2: # limits go from -180 to 180 - axh.set_xlim(-180 * np.pi / 180, 180 * np.pi / 180) + #axh.set_xlim(-180 * np.pi / 180, 180 * np.pi / 180) tr_median, tr_mode, tr_mean = self.get_stats(tr, tr_hist, diff --git a/mtpy/utils/calculator.py b/mtpy/utils/calculator.py index 6a78037ff..7528bc6c5 100644 --- a/mtpy/utils/calculator.py +++ b/mtpy/utils/calculator.py @@ -42,13 +42,14 @@ def roundsf(number, sf): round a number to a specified number of significant figures (sf) """ # can't have < 1 s.f. - sf = max(sf,1.) + sf = max(sf, 1.) rounding = int(np.ceil(-np.log10(number) + sf - 1.)) return np.round(number, rounding) -def get_period_list(period_min,period_max,periods_per_decade,include_outside_range=True): +def get_period_list(period_min, period_max, periods_per_decade, + include_outside_range=True): """ get a list of values (e.g. periods), evenly spaced in log space and including values on multiples of 10 @@ -118,7 +119,8 @@ def nearest_index(val,array): return np.where(diff==min(diff))[0][0] -def make_log_increasing_array(z1_layer, target_depth, n_layers, increment_factor=0.999): +def make_log_increasing_array(z1_layer, target_depth, n_layers, + increment_factor=0.999): """ create depth array with log increasing cells, down to target depth, inputs are z1_layer thickness, target depth, number of layers (n_layers) @@ -189,17 +191,17 @@ def invertmatrix_incl_errors(inmatrix, inmatrix_err=None): for l in range(2): #each entry has 4 summands - err += np.abs (- inv_matrix[i,k] * inv_matrix[l,j] * inmatrix_err[k,l]) - + err += np.abs(- inv_matrix[i,k] * inv_matrix[l,j] *\ + inmatrix_err[k,l]) inv_matrix_err[i,j] = err - return inv_matrix, inv_matrix_err def rhophi2z(rho, phi, freq): """ - Convert impedance-style information given in Rho/Phi format into complex valued Z. + Convert impedance-style information given in Rho/Phi format into + complex valued Z. Input: rho - 2x2 array (real) - in Ohm m @@ -228,7 +230,8 @@ def rhophi2z(rho, phi, freq): return z -def compute_determinant_error(z_array, z_err_array, method = 'theoretical', repeats=1000): +def compute_determinant_error(z_array, z_err_array, method='theoretical', + repeats=1000): """ compute the error of the determinant of z using a stochastic method seed random z arrays with a normal distribution around the input array @@ -257,33 +260,44 @@ def compute_determinant_error(z_array, z_err_array, method = 'theoretical', repe error = np.std(detlist,axis=0) else: - error = np.abs(z_err_array[:,0,0]*np.abs(z_array[:,1,1]) + z_err_array[:,1,1]*np.abs(z_array[:,0,0]) \ - - z_err_array[:,0,1]*np.abs(z_array[:,1,0]) - z_err_array[:,1,0]*np.abs(z_array[:,0,1])) + error = np.abs(z_err_array[:, 0, 0] * np.abs(z_array[:, 1, 1]) +\ + z_err_array[:, 1, 1] * np.abs(z_array[:, 0, 0]) \ + - z_err_array[:, 0, 1]*np.abs(z_array[:, 1, 0]) - \ + z_err_array[:, 1, 0]*np.abs(z_array[:, 0, 1])) return error - - - def propagate_error_polar2rect(r,r_error,phi, phi_error): """ - Find error estimations for the transformation from polar to cartesian coordinates. - - Uncertainties in polar representation define a section of an annulus. Find the 4 corners of this section and additionally the outer boundary point, which is defined by phi = phi0, rho = rho0 + sigma rho. - The cartesian "box" defining the uncertainties in x,y is the outer bound around the annulus section, defined by the four outermost points. So check the four corners as well as the outer boundary edge of the section to find the extrema in x znd y. These give you the sigma_x/y. + Find error estimations for the transformation from polar to cartesian + coordinates. + + Uncertainties in polar representation define a section of an annulus. + Find the 4 corners of this section and additionally the outer boundary + point, which is defined by phi = phi0, rho = rho0 + sigma rho. + The cartesian "box" defining the uncertainties in x,y is the outer + bound around the annulus section, defined by the four outermost + points. So check the four corners as well as the outer boundary edge + of the section to find the extrema in x znd y. These give you the + sigma_x/y. """ - corners = [ ( np.real(cmath.rect(r-r_error, phi-phi_error)), np.imag(cmath.rect(r-r_error, phi-phi_error))),\ - ( np.real(cmath.rect(r+r_error, phi-phi_error)), np.imag(cmath.rect(r+r_error, phi-phi_error))),\ - ( np.real(cmath.rect(r+r_error, phi+phi_error)), np.imag(cmath.rect(r+r_error, phi+phi_error))),\ - ( np.real(cmath.rect(r-r_error, phi+phi_error)), np.imag(cmath.rect(r-r_error, phi+phi_error))),\ - ( np.real(cmath.rect(r+r_error, phi)), np.imag(cmath.rect(r+r_error, phi))) ] + corners = [( np.real(cmath.rect(r-r_error, phi-phi_error)), + np.imag(cmath.rect(r-r_error, phi-phi_error))), + (np.real(cmath.rect(r+r_error, phi-phi_error)), + np.imag(cmath.rect(r+r_error, phi-phi_error))), + (np.real(cmath.rect(r+r_error, phi+phi_error)), + np.imag(cmath.rect(r+r_error, phi+phi_error))), + (np.real(cmath.rect(r-r_error, phi+phi_error)), + np.imag(cmath.rect(r-r_error, phi+phi_error))), + (np.real(cmath.rect(r+r_error, phi)), + np.imag(cmath.rect(r+r_error, phi))) ] lo_x = [i[0] for i in corners] lo_y = [i[1] for i in corners] - point = (np.real(cmath.rect(r, phi)), np.imag(cmath.rect(r, phi)) ) + point = (np.real(cmath.rect(r, phi)), np.imag(cmath.rect(r, phi))) lo_xdiffs = [ abs(point[0] - i) for i in lo_x] lo_ydiffs = [ abs(point[1] - i) for i in lo_y] @@ -292,25 +306,24 @@ def propagate_error_polar2rect(r,r_error,phi, phi_error): return xerr, yerr - - - def propagate_error_rect2polar(x,x_error,y, y_error): # x_error, y_error define a rectangular uncertainty box - # rho error is the difference between the closest and furthest point of the box (w.r.t. the origin) + # rho error is the difference between the closest and furthest point + # of the box (w.r.t. the origin) # approximation: just take corners and midpoint of edges - lo_points = [ (x + x_error, y), (x - x_error, y), (x, y - y_error ), (x, y + y_error ),\ - (x - x_error, y - y_error) ,(x + x_error, y - y_error) ,(x + x_error, y + y_error) ,(x - x_error, y + y_error) ] - + lo_points = [(x + x_error, y), (x - x_error, y), + (x, y - y_error ), (x, y + y_error ), + (x - x_error, y - y_error), (x + x_error, y - y_error), + (x + x_error, y + y_error), (x - x_error, y + y_error)] #check, if origin is within the box: origin_in_box = False if x_error >= np.abs(x) and y_error >= np.abs(y): origin_in_box = True - lo_polar_points = [ cmath.polar(np.complex(*i)) for i in lo_points ] + lo_polar_points = [cmath.polar(np.complex(*i)) for i in lo_points] lo_rho = [i[0] for i in lo_polar_points ] lo_phi = [math.degrees(i[1])%360 for i in lo_polar_points ] @@ -319,14 +332,13 @@ def propagate_error_rect2polar(x,x_error,y, y_error): phi_err = 0.5*(max(lo_phi) - min(lo_phi)) if (270 < max(lo_phi) < 360 ) and (0 < min(lo_phi) < 90 ): - tmp1 = [ i for i in lo_phi if (0 < i < 90) ] - tmp4 = [ i for i in lo_phi if (270 < i < 360) ] - phi_err = 0.5*((max(tmp1) - min(tmp4))%360) + tmp1 = [i for i in lo_phi if (0 < i < 90)] + tmp4 = [i for i in lo_phi if (270 < i < 360)] + phi_err = 0.5 * ((max(tmp1) - min(tmp4)) % 360) if phi_err > 180: - #print phi_err,' -> ',(-phi_err)%360 - phi_err = (-phi_err)%360 + phi_err = (-phi_err) % 360 if origin_in_box is True: #largest rho: @@ -336,8 +348,6 @@ def propagate_error_rect2polar(x,x_error,y, y_error): return rho_err, phi_err - - def z_error2r_phi_error(z_real, z_imag, error): """ Error estimation from rectangular to polar coordinates. @@ -366,9 +376,9 @@ def z_error2r_phi_error(z_real, z_imag, error): res_rel_err = 2.*z_rel_err - #if the relative error of the amplitude is >=100% that means that the relative - #error of the resistivity is 200% - that is then equivalent to an uncertainty - #in the phase angle of 90 degrees: + # if the relative error of the amplitude is >=100% that means that the + # relative error of the resistivity is 200% - that is then equivalent + # to an uncertainty in the phase angle of 90 degrees if np.iterable(z_real): phi_err = np.degrees(np.arctan(z_rel_err)) phi_err[res_rel_err > 1.] = 90. @@ -378,58 +388,50 @@ def z_error2r_phi_error(z_real, z_imag, error): phi_err = 90 else: phi_err = np.degrees(np.arctan(z_rel_err)) - - + return res_rel_err, phi_err - def old_z_error2r_phi_error(x,x_error,y, y_error): """ - Error estimation from rect to polar, but with small variation needed for - MT: the so called 'relative phase error' is NOT the relative phase error, - but the ABSOLUTE uncertainty in the angle that corresponds to the relative - error in the amplitude. - - So, here we calculate the transformation from rect to polar coordinates, - esp. the absolute/length of the value. Then we find the uncertainty in - this length and calculate the relative error of this. The relative error of - the resistivity will be double this value, because it's calculated by taking - the square of this length. - - The relative uncertainty in length defines a circle around (x,y) - (APPROXIMATION!). The uncertainty in phi is now the absolute of the - angle beween the vector to (x,y) and the origin-vector tangential to the - circle. - BUT....since the phase angle uncertainty is interpreted with regard to - the resistivity and not the Z-amplitude, we have to look at the square of - the length, i.e. the relative error in question has to be halfed to get - the correct relationship between resistivity and phase errors!!. + Error estimation from rect to polar, but with small variation needed for + MT: the so called 'relative phase error' is NOT the relative phase error, + but the ABSOLUTE uncertainty in the angle that corresponds to the relative + error in the amplitude. + + So, here we calculate the transformation from rect to polar coordinates, + esp. the absolute/length of the value. Then we find the uncertainty in + this length and calculate the relative error of this. The relative error of + the resistivity will be double this value, because it's calculated by taking + the square of this length. + + The relative uncertainty in length defines a circle around (x,y) + (APPROXIMATION!). The uncertainty in phi is now the absolute of the + angle beween the vector to (x,y) and the origin-vector tangential to the + circle. + BUT....since the phase angle uncertainty is interpreted with regard to + the resistivity and not the Z-amplitude, we have to look at the square of + the length, i.e. the relative error in question has to be halfed to get + the correct relationship between resistivity and phase errors!!. """ # x_error, y_error define a rectangular uncertainty box - # rho error is the difference between the closest and furthest point of the box (w.r.t. the origin) - # approximation: just take corners and midpoint of edges - lo_points = [ (x + x_error, y), (x - x_error, y), (x, y - y_error ), - (x, y + y_error ), (x - x_error, y - y_error) , - (x + x_error, y - y_error) ,(x + x_error, y + y_error) , - (x - x_error, y + y_error) ] - + # rho error is the difference between the closest and furthest point + # of the box (w.r.t. the origin) approximation: just take corners and + # midpoint of edges + lo_points = [(x + x_error, y), (x - x_error, y), + (x, y - y_error ), (x, y + y_error ), + (x - x_error, y - y_error), (x + x_error, y - y_error), + (x + x_error, y + y_error), (x - x_error, y + y_error)] - #check, if origin is within the box: - origin_in_box = False - if x_error >= np.abs(x) and y_error >= np.abs(y): - origin_in_box = True + lo_polar_points = [cmath.polar(np.complex(*i)) for i in lo_points] - lo_polar_points = [ cmath.polar(np.complex(*i)) for i in lo_points ] - - lo_rho = [i[0] for i in lo_polar_points ] - lo_phi = [math.degrees(i[1])%360 for i in lo_polar_points ] + lo_rho = [i[0] for i in lo_polar_points] #uncertainty in amplitude is defined by half the diameter of the box around x,y - rho_err = 0.5*(max(lo_rho) - min(lo_rho) ) + rho_err = 0.5 * (max(lo_rho) - min(lo_rho) ) rho = cmath.polar(np.complex(x,y))[0] try: @@ -443,14 +445,11 @@ def old_z_error2r_phi_error(x,x_error,y, y_error): if rel_error_rho > 1.: phi_err = 90 else: - phi_err = np.degrees(np.arcsin( rel_error_rho)) + phi_err = np.degrees(np.arcsin(rel_error_rho)) return rho_err, phi_err - - - #rotation: #1. rotation positive in clockwise direction #2. orientation of new X-axis X' given by rotation angle @@ -470,96 +469,111 @@ def old_z_error2r_phi_error(x,x_error,y, y_error): # That is NOT the same as the rotated error matrix Zerr (although the result is similar) -def rotatematrix_incl_errors(inmatrix, angle, inmatrix_err = None) : - - if inmatrix is None : - raise MTex.MTpyError_inputarguments('Matrix AND eror matrix must be defined') +def rotate_matrix_with_errors(in_matrix, angle, error=None): + """ + + Rotate a matrix including errors clockwise given an angle in degrees. + + + :param inmatrix: DESCRIPTION + :type inmatrix: TYPE + :param angle: DESCRIPTION + :type angle: TYPE + :param inmatrix_err: DESCRIPTION, defaults to None + :type inmatrix_err: TYPE, optional + :raises MTex: DESCRIPTION + :return: DESCRIPTION + :rtype: TYPE - if (inmatrix_err is not None) and (inmatrix.shape != inmatrix_err.shape): - raise MTex.MTpyError_inputarguments('Matrix and err-matrix shapes do not match: %s - %s'%(str(inmatrix.shape), str(inmatrix_err.shape))) + """ + + if in_matrix is None: + raise MTex.MTpyError_inputarguments('Matrix must be defined') + if (error is not None) and (in_matrix.shape != error.shape): + msg = 'matricies are not the same shape in_matrix={0}, err={1}'.format( + in_matrix.shape, error.shape) + raise MTex.MTpyError_inputarguments(msg) try: - degreeangle = angle % 360 - except: - raise MTex.MTpyError_inputarguments('"Angle" must be a valid number (in degrees)') - - phi = math.radians(degreeangle) + phi = np.deg2rad(float(angle) % 360) + except TypeError: + raise MTex.MTpyError_inputarguments('"Angle" must be a float') cphi = np.cos(phi) sphi = np.sin(phi) - rotmat = np.array([[ cphi, sphi], [-sphi, cphi]]) - rotated_matrix = np.dot(np.dot(rotmat, inmatrix), np.linalg.inv(rotmat)) + # clockwise rotation matrix is given by [[cos, -sin], [sin, cos]] + rot_mat = np.array([[ cphi, -sphi], [sphi, cphi]]) + rotated_matrix = np.dot(np.dot(rot_mat, in_matrix), np.linalg.inv(rot_mat)) - errmat = None - if (inmatrix_err is not None) : - err_orig = np.real(inmatrix_err) - errmat = np.zeros_like(inmatrix_err) + err_mat = None + if (error is not None) : + err_orig = np.real(error) + err_mat = np.zeros_like(error) # standard propagation of errors: - errmat[0,0] = np.sqrt( (cphi**2 * err_orig[0,0])**2 + \ - (cphi * sphi * err_orig[0,1])**2 + \ - (cphi * sphi * err_orig[1,0])**2 + \ - (sphi**2 * err_orig[1,1])**2 ) - errmat[0,1] = np.sqrt( (cphi**2 * err_orig[0,1])**2 + \ - (cphi * sphi * err_orig[1,1])**2 + \ - (cphi * sphi * err_orig[0,0])**2 + \ - (sphi**2 * err_orig[1,0])**2 ) - errmat[1,0] = np.sqrt( (cphi**2 * err_orig[1,0])**2 + \ - (cphi * sphi * err_orig[1,1])**2 +\ - (cphi * sphi * err_orig[0,0])**2 + \ - (sphi**2 * err_orig[0,1])**2 ) - errmat[1,1] = np.sqrt( (cphi**2 * err_orig[1,1])**2 + \ - (cphi * sphi * err_orig[0,1])**2 + \ - (cphi * sphi * err_orig[1,0])**2 + \ - (sphi**2 * err_orig[0,0])**2 ) - - return rotated_matrix, errmat - - -def rotatevector_incl_errors(invector, angle, invector_err = None): + err_mat[0,0] = np.sqrt((cphi**2 * err_orig[0, 0])**2 + \ + (cphi * sphi * err_orig[0, 1])**2 + \ + (cphi * sphi * err_orig[1, 0])**2 + \ + (sphi**2 * err_orig[1, 1])**2) + err_mat[0,1] = np.sqrt((cphi**2 * err_orig[0, 1])**2 + \ + (cphi * sphi * err_orig[1, 1])**2 + \ + (cphi * sphi * err_orig[0, 0])**2 + \ + (sphi**2 * err_orig[1, 0])**2) + err_mat[1,0] = np.sqrt((cphi**2 * err_orig[1, 0])**2 + \ + (cphi * sphi * err_orig[1, 1])**2 +\ + (cphi * sphi * err_orig[0, 0])**2 + \ + (sphi**2 * err_orig[0, 1])**2) + err_mat[1,1] = np.sqrt((cphi**2 * err_orig[1, 1])**2 + \ + (cphi * sphi * err_orig[0, 1])**2 + \ + (cphi * sphi * err_orig[1, 0])**2 + \ + (sphi**2 * err_orig[0, 0])**2) + + return rotated_matrix, err_mat + + +def rotate_vector_with_errors(in_vector, angle, error=None): #check for row or column vector - if invector is None : - raise MTex.MTpyError_inputarguments('Vector AND error-vector must be defined') - - if (invector_err is not None) and (invector.shape != invector_err.shape): - raise MTex.MTpyError_inputarguments('Vector and errror-vector shapes do not match: %s - %s'%(str(invector.shape), str(invector_err.shape))) + if in_vector is None: + raise MTex.MTpyError_inputarguments('Vector AND error-vector must'+ + ' be defined') + if (error is not None) and (in_vector.shape != error.shape): + msg = 'matricies are not the same shape in_matrix={0}, err={1}'.format( + in_vector.shape, error.shape) + raise MTex.MTpyError_inputarguments(msg) + try: - degreeangle = angle%360 - except: - raise MTex.MTpyError_inputarguments('"Angle" must be a valid number (in degrees)') + phi = np.deg2rad(float(angle) % 360) + except TypeError: + raise MTex.MTpyError_inputarguments('"Angle" must be a float') - phi = math.radians(degreeangle) - cphi = np.cos(phi) sphi = np.sin(phi) - rotmat = np.array([[ cphi, sphi],[-sphi, cphi]]) + rot_mat = np.array([[ cphi, -sphi],[sphi, cphi]]) - if invector.shape == (1, 2): - rotated_vector = np.dot( invector, np.linalg.inv(rotmat) ) + if in_vector.shape == (1, 2): + rotated_vector = np.dot(in_vector, np.linalg.inv(rot_mat) ) else: - rotated_vector = np.dot( rotmat, invector ) - - - errvec = None - if (invector_err is not None) : - errvec = np.zeros_like(invector_err) - - if invector_err.shape == (1, 2): - errvec = np.dot(invector_err, np.abs(np.linalg.inv(rotmat))) - else: - errvec = np.dot(np.abs(rotmat), invector_err ) + rotated_vector = np.dot(rot_mat, in_vector) + err_vector = None + if (error is not None): + err_vector = np.zeros_like(error) - return rotated_vector, errvec + if error.shape == (1, 2): + err_vector = np.dot(error, np.abs(np.linalg.inv(rot_mat))) + else: + err_vector = np.dot(np.abs(rot_mat), error) + return rotated_vector, err_vector -def multiplymatrices_incl_errors(inmatrix1, inmatrix2, inmatrix1_err = None,inmatrix2_err = None ): +def multiplymatrices_incl_errors(inmatrix1, inmatrix2, + inmatrix1_err=None, inmatrix2_err=None): if inmatrix1 is None or inmatrix2 is None: raise MTex.MTpyError_inputarguments('ERROR - two 2x2 arrays needed as input') @@ -591,20 +605,20 @@ def multiplymatrices_incl_errors(inmatrix1, inmatrix2, inmatrix1_err = None,inma def reorient_data2D(x_values, y_values, x_sensor_angle = 0 , y_sensor_angle = 90): """ - Re-orient time series data of a sensor pair, which has not been in default (x=0, y=90) orientation. + Re-orient time series data of a sensor pair, which has not been in default (x=0, y=90) orientation. - Input: - - x-values - Numpy array - - y-values - Numpy array - Note: same length for both! - If not, the shorter length is taken + Input: + - x-values - Numpy array + - y-values - Numpy array + Note: same length for both! - If not, the shorter length is taken - Optional: - - Angle of the x-sensor - measured in degrees, clockwise from North (0) - - Angle of the y-sensor - measured in degrees, clockwise from North (0) + Optional: + - Angle of the x-sensor - measured in degrees, clockwise from North (0) + - Angle of the y-sensor - measured in degrees, clockwise from North (0) - Output: - - corrected x-values (North) - - corrected y-values (East) + Output: + - corrected x-values (North) + - corrected y-values (East) """ x_values = np.array(x_values) From 5dd664f05021e64389768c2b508ca55cf7daf94b Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 11 May 2020 16:15:59 -0700 Subject: [PATCH 07/20] fixed plot_pt, plot_pt_pseudosection --- mtpy/imaging/phase_tensor_maps.py | 38 +++++++------------- mtpy/imaging/phase_tensor_pseudosection.py | 42 +++++++++++----------- mtpy/imaging/plotstrike.py | 8 ++--- 3 files changed, 38 insertions(+), 50 deletions(-) diff --git a/mtpy/imaging/phase_tensor_maps.py b/mtpy/imaging/phase_tensor_maps.py index b4c6af295..113a1f87d 100644 --- a/mtpy/imaging/phase_tensor_maps.py +++ b/mtpy/imaging/phase_tensor_maps.py @@ -518,34 +518,22 @@ def __init__(self, **kwargs): self.plot() # self.save_figure(self.save_fn, file_format='png') - # ---need to rotate data on setting rotz - def _set_rot_z(self, rot_z): + #---need to rotate data on setting rotz + @property + def rotation_angle(self): + return self._rotation_angle + + @rotation_angle.setter + def rotation_angle(self, value): """ - need to rotate data when setting z + only a single value is allowed """ - - # 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)""") + mt.rotation_angle = value + + self._rotation_angle = value + + self.make_strike_array() # ----------------------------------------------- # The main plot method for this module diff --git a/mtpy/imaging/phase_tensor_pseudosection.py b/mtpy/imaging/phase_tensor_pseudosection.py index f13cdec15..9337074cb 100644 --- a/mtpy/imaging/phase_tensor_pseudosection.py +++ b/mtpy/imaging/phase_tensor_pseudosection.py @@ -435,6 +435,7 @@ def plot(self, show=True): # create a plot instance self.fig = plt.figure(self.fig_num, self.fig_size, dpi=self.fig_dpi) + self.fig.clf() self.ax = self.fig.add_subplot(1, 1, 1, aspect='equal') # FZ: control tick rotation=30 not that good @@ -528,26 +529,27 @@ def plot(self, show=True): # get the properties to color the ellipses by if self.ellipse_colorby == 'phimin': - colorarray = pt.phimin[::-1] + color_array = pt.phimin[::-1] elif self.ellipse_colorby == 'phimax': - colorarray = pt.phimin[::-1] + color_array = pt.phimin[::-1] elif self.ellipse_colorby == 'phidet': - colorarray = np.sqrt(abs(pt.det[::-1])) * (180 / np.pi) + color_array = np.sqrt(abs(pt.det[::-1])) * (180 / np.pi) elif self.ellipse_colorby == 'skew' or \ self.ellipse_colorby == 'skew_seg': - colorarray = pt.beta[::-1] + color_array = pt.beta[::-1] elif self.ellipse_colorby == 'normalized_skew' or \ self.ellipse_colorby == 'normalized_skew_seg': - colorarray = 2 * pt.beta[::-1] + color_array = 2 * pt.beta[::-1] elif self.ellipse_colorby == 'ellipticity': - colorarray = pt.ellipticity[::-1] + color_array = pt.ellipticity[::-1] elif self.ellipse_colorby in ['strike', 'azimuth']: - colorarray = pt.azimuth[::-1] % 180 + color_array = pt.azimuth[::-1] % 180 + color_array[color_array > 90] -= 180 else: raise NameError(self.ellipse_colorby + ' is not supported') @@ -562,8 +564,8 @@ def plot(self, show=True): plot_periodlist = periodlist # get min and max of the color array for scaling later - minlist.append(min(colorarray)) - maxlist.append(max(colorarray)) + minlist.append(min(color_array)) + maxlist.append(max(color_array)) for jj, ff in enumerate(periodlist): @@ -587,14 +589,14 @@ def plot(self, show=True): # get ellipse color if cmap.find('seg') > 0: - ellipd.set_facecolor(mtcl.get_plot_color(colorarray[jj], + ellipd.set_facecolor(mtcl.get_plot_color(color_array[jj], self.ellipse_colorby, cmap, ckmin, ckmax, bounds=bounds)) else: - ellipd.set_facecolor(mtcl.get_plot_color(colorarray[jj], + ellipd.set_facecolor(mtcl.get_plot_color(color_array[jj], self.ellipse_colorby, cmap, ckmin, @@ -1159,45 +1161,45 @@ def writeTextFiles(self, save_path=None, ptol=0.10): tpiazlines.append(tpiazline + '\n') # write files - skfid = file(os.path.join(svpath, 'PseudoSection.skew'), 'w') + skfid = open(os.path.join(svpath, 'PseudoSection.skew'), 'w') skfid.writelines(sklines) skfid.close() - phiminfid = file(os.path.join(svpath, 'PseudoSection.phimin'), 'w') + phiminfid = open(os.path.join(svpath, 'PseudoSection.phimin'), 'w') phiminfid.writelines(phiminlines) phiminfid.close() - phimaxfid = file(os.path.join(svpath, 'PseudoSection.phimax'), + phimaxfid = open(os.path.join(svpath, 'PseudoSection.phimax'), 'w') phimaxfid.writelines(phimaxlines) phimaxfid.close() - ellipfid = file(os.path.join(svpath, 'PseudoSection.ellipticity'), + ellipfid = open(os.path.join(svpath, 'PseudoSection.ellipticity'), 'w') ellipfid.writelines(elliplines) ellipfid.close() - azfid = file(os.path.join(svpath, 'PseudoSection.azimuth'), + azfid = open(os.path.join(svpath, 'PseudoSection.azimuth'), 'w') azfid.writelines(azimlines) azfid.close() - tprfid = file(os.path.join(svpath, 'PseudoSection.tipper_mag_real'), + tprfid = open(os.path.join(svpath, 'PseudoSection.tipper_mag_real'), 'w') tprfid.writelines(tprlines) tprfid.close() - tprazfid = file(os.path.join(svpath, 'PseudoSection.tipper_ang_real'), + tprazfid = open(os.path.join(svpath, 'PseudoSection.tipper_ang_real'), 'w') tprazfid.writelines(tprazlines) tprazfid.close() - tpifid = file(os.path.join(svpath, 'PseudoSection.tipper_mag_imag'), + tpifid = open(os.path.join(svpath, 'PseudoSection.tipper_mag_imag'), 'w') tpifid.writelines(tpilines) tpifid.close() - tpiazfid = file(os.path.join(svpath, 'PseudoSection.tipper_ang_imag'), + tpiazfid = open(os.path.join(svpath, 'PseudoSection.tipper_ang_imag'), 'w') tpiazfid.writelines(tpiazlines) tpiazfid.close() diff --git a/mtpy/imaging/plotstrike.py b/mtpy/imaging/plotstrike.py index f95324239..4a13038a8 100644 --- a/mtpy/imaging/plotstrike.py +++ b/mtpy/imaging/plotstrike.py @@ -16,8 +16,6 @@ import mtpy.imaging.mtplottools as mtpl #============================================================================== - - class PlotStrike(object): """ PlotStrike will plot the strike estimated from the invariants, phase tensor @@ -1288,9 +1286,9 @@ def writeTextFiles(self, save_path=None): 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') + invfid = open(os.path.join(svpath, 'Strike.invariants'), 'w') + ptfid = open(os.path.join(svpath, 'Strike.pt'), 'w') + tpfid = open(os.path.join(svpath, 'Strike.tipper'), 'w') #---> write strike from the invariants # == > mean From 4ee148c76c8f72aed13c484bc4e7bee471a35639 Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 11 May 2020 16:28:18 -0700 Subject: [PATCH 08/20] fixed plot_pt, plot_pt_maps --- mtpy/imaging/phase_tensor_maps.py | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/mtpy/imaging/phase_tensor_maps.py b/mtpy/imaging/phase_tensor_maps.py index 113a1f87d..c48ea67fa 100644 --- a/mtpy/imaging/phase_tensor_maps.py +++ b/mtpy/imaging/phase_tensor_maps.py @@ -532,8 +532,15 @@ def rotation_angle(self, value): mt.rotation_angle = value self._rotation_angle = value - - self.make_strike_array() + + def fold_strike(self, strike): + """ + + """ + strike = strike % 180 + strike[np.where(strike > 90)] -= 180 + + return strike # ----------------------------------------------- # The main plot method for this module @@ -565,14 +572,6 @@ def plot(self, fig=None, save_path=None, show=True, plt.rcParams['figure.subplot.top'] = .93 plt.rcParams['figure.subplot.wspace'] = .55 plt.rcParams['figure.subplot.hspace'] = .70 - # FZ: tweaks to make plot positioned better - # plt.rcParams['font.size']=self.font_size - # plt.rcParams['figure.subplot.left']=.1 - # plt.rcParams['figure.subplot.right']=.90 - # plt.rcParams['figure.subplot.bottom']=.2 - # plt.rcParams['figure.subplot.top']=.90 - # plt.rcParams['figure.subplot.wspace']=.70 - # plt.rcParams['figure.subplot.hspace']=.70 lpfig = None lpax = None @@ -804,6 +803,9 @@ def plot(self, fig=None, save_path=None, show=True, elif self.ellipse_colorby == 'ellipticity': colorarray = pt.ellipticity[jj] + + elif self.ellipse_colorby in ['strike', 'azimuth']: + colorarray = self.fold_strike(self.pt.azimuth) else: raise NameError(self.ellipse_colorby + ' is not supported') @@ -1256,7 +1258,7 @@ def redraw_plot(self): use this function if you updated some attributes and want to re-plot. """ - plt.close(self.fig) + self.fig.clf() self.plot() def export_params_to_file(self, save_path=None): From e303ec3c6836e7d59e76e6044ca1871bb39441e2 Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 11 May 2020 18:01:17 -0700 Subject: [PATCH 09/20] fixed plot_pt, plot_pt_maps --- mtpy/analysis/pt.py | 199 ++++++------------------------ mtpy/analysis/zinvariants.py | 5 +- mtpy/imaging/phase_tensor_maps.py | 7 +- mtpy/imaging/plot_mt_response.py | 41 +++--- mtpy/imaging/plotresponse.py | 12 +- mtpy/imaging/plotstrike.py | 3 +- 6 files changed, 75 insertions(+), 192 deletions(-) diff --git a/mtpy/analysis/pt.py b/mtpy/analysis/pt.py index ca57d2c02..7e8a58cfd 100644 --- a/mtpy/analysis/pt.py +++ b/mtpy/analysis/pt.py @@ -108,7 +108,12 @@ def __init__(self, pt_array=None, pt_err_array=None, z_array=None, # define get/set functions and properties # ========================================================================== # ---phase tensor array---------------------------------------- - def _set_pt(self, pt_array): + @property + def pt(self): + return self._pt + + @pt.setter + def pt(self, pt_array): """ Set the attribute 'pt'. @@ -176,13 +181,13 @@ def _set_pt(self, pt_array): else: pass - def _get_pt(self): + # ---phase tensor Error----------------------------------------------------- + @property + def pt_err(self): return self._pt - pt = property(_get_pt, _set_pt, doc="Phase tensor array") - - # ---phase tensor Error----------------------------------------------------- - def _set_pt_err(self, pt_err_array): + @pt_err.setter + def pt_err(self, pt_err_array): """ Set the attribute 'pt_err'. @@ -224,14 +229,13 @@ def _set_pt_err(self, pt_err_array): else: pass - def _get_pt_err(self): - return self._pt_err - - pt_err = property(_get_pt_err, _set_pt_err, - doc='Phase tensor error array, must be same shape as pt') - # ---freq------------------------------------------------------------ - def _set_freq(self, lo_freq): + @property + def freq(self): + return self._freq + + @freq.setter + def freq(self, lo_freq): """ Set array of freq. @@ -253,11 +257,6 @@ def _set_freq(self, lo_freq): except: self._freq = None - def _get_freq(self): - return self._freq - - freq = property(_get_freq, _set_freq, doc="freq array") - # ---z_object--------------------------------------------------------------- def set_z_object(self, z_object): @@ -278,12 +277,7 @@ def set_z_object(self, z_object): self._pt[idx_f], self._pt_err[idx_f] = z2pt(self._z[idx_f], self._z_err[idx_f]) except MTex.MTpyError_PT: - try: - print('Singular Matrix at {0:.5g} Hz'.format( - self._freq[idx_f])) - except AttributeError: - print('Computed singular matrix') - print(' --> pt[{0}]=np.zeros((2,2))'.format(idx_f)) + pass # --> if there is not error to the impedance tensor else: @@ -291,30 +285,19 @@ def set_z_object(self, z_object): try: self._pt[idx_f] = z2pt(self._z[idx_f])[0] except MTex.MTpyError_PT: - try: - print('Singular Matrix at {0:.5g}'.format( - self._freq[idx_f])) - except AttributeError: - print('Computed singular matrix') - print(' --> pt[{0}]=np.zeros((2,2))'.format(idx_f)) + pass self.rotation_angle = z_object.rotation_angle - # def _get_z_object(self): - # z_object = MTz.Z(z_array=self._z, z_err_array=self._z_err) - # z_object.freq = self._freq - # z_object.rotation_angle = self.rotation_angle - - # return z_object - - # _z_object = property(_get_z_object, _set_z_object, - # doc="class mtpy.core.z.Z") - - # ---z array--------------------------------------------------------------- - def _set_z(self, z_array): + @property + def z(self): + return self._z + + @z.setter + def z(self, z_array): """ - Set Z array as PhaseTensor object attribute. + Set Z array as PhaseTensor object attribute. """ self._z = z_array @@ -327,12 +310,7 @@ def _set_z(self, z_array): self._pt[idx_f], self._pt_err[idx_f] = z2pt(self._z[idx_f], self._z_err[idx_f]) except MTex.MTpyError_PT: - try: - print('Singular Matrix at {0:.5g} Hz'.format( - self._freq[idx_f])) - except AttributeError: - print('Computed singular matrix') - print(' --> pt[{0}]=np.zeros((2,2))'.format(idx_f)) + pass # --> if there is not error to the impedance tensor elif self._z is not None: @@ -340,20 +318,14 @@ def _set_z(self, z_array): try: self._pt[idx_f] = z2pt(self._z[idx_f])[0] except MTex.MTpyError_PT: - try: - print('Singular Matrix at {0:.5g}'.format( - self._freq[idx_f])) - except AttributeError: - print('Computed singular matrix') - print(' --> pt[{0}]=np.zeros((2,2))'.format(idx_f)) - - # def _get_z(self): - # return self._z + pass - # z = property(_get_z, _set_z, - # doc="impedance tensor numpy.array((nf, 2, 2))") # ---Z Error array--------------------------------------------------------------- + def z_err(self): + return self._z_err + + @z_err.setter def _set_z_err(self, z_err_array): """ Set Z-error array as PhaseTensor object attribute. @@ -373,12 +345,7 @@ def _set_z_err(self, z_err_array): self.pt[idx_f], self.pt_err[idx_f] = z2pt(self._z[idx_f], self._z_err[idx_f]) except MTex.MTpyError_PT: - try: - print('Singular Matrix at {0:.5g} Hz'.format( - self._freq[idx_f])) - except AttributeError: - print('Computed singular matrix') - print(' --> pt[{0}]=np.zeros((2,2))'.format(idx_f)) + pass # --> if there is not error to the impedance tensor else: @@ -386,20 +353,7 @@ def _set_z_err(self, z_err_array): try: self._pt[idx_f] = z2pt(self._z[idx_f])[0] except MTex.MTpyError_PT: - try: - print('Singular Matrix at {0:.5g}'.format( - self._freq[idx_f])) - except AttributeError: - print('Computed singular matrix') - print(' --> pt[{0}]=np.zeros((2,2))'.format(idx_f)) - - # def _get_z_err(self): - # return self._z_err - - # z_err = property(_get_z_err, _set_z_err, - # doc="impedance tensor numpy.array((nf, 2, 2))") - - + pass # ========================================================================== # define get methods for read only properties @@ -1309,87 +1263,4 @@ def z2pt(z_array, z_err_array=None): np.abs( -pt_array[idx_f,1,1] * realz[0,0] * z_err_array[1,1]) + \ np.abs( realz[0,0] * z_err_array[1,1]) + np.abs( -realz[1,0] * z_err_array[0,1]) ) - return pt_array, pt_err_array - - -def z_object2pt(z_object): - """ - Calculate Phase Tensor from Z object (incl. uncertainties) - - Input: - - Z-object : instance of the MTpy Z class - - - Return: - - PT object - """ - # - PT : nx2x2 real valued Numpy array - # - PT-error : nx2x2 real valued Numpy array - - # """ - - try: - p = PhaseTensor(z_object=z_object) - except: - raise MTex.MTpyError_Z('Input argument is not a valid instance of the Z class') - - # pt_array = p.pt - # pterr_array = p.pterr - - # return pt_array, pterr_array - return p - - -def _edi_object2pt(edi_object): - """ - Calculate Phase Tensor from Edi object (incl. uncertainties) - - Input: - - Edi-object : instance of the MTpy Edi class - - Return: - - PT : nx2x2 real valued Numpy array - - PT-error : nx2x2 real valued Numpy array - - """ - - if not isinstance(edi_object, MTedi.Edi): - raise MTex.MTpyError_EDI('Input argument is not an instance of the Edi class') - p = PhaseTensor(edi_object=edi_object) - - pt_array = p.pt - - pterr_array = p.pterr - - return pt_array, pterr_array - - -def edi_file2pt(filename): - """ - Calculate Phase Tensor from Edi-file (incl. uncertainties) - - Input: - - Edi-file : full path to the Edi-file - - Return: - - PT object - - """ - # Return: - # - PT : nx2x2 real valued Numpy array - # - PT-error : nx2x2 real valued Numpy array - - # """ - - e = MTedi.Edi() - e.readfile(filename) - - p = PhaseTensor(z_object=e.Z) - - # pt_array = p.pt - - # pterr_array = p.pterr - - # return pt_array, pterr_array - - return p + return pt_array, pt_err_array \ No newline at end of file diff --git a/mtpy/analysis/zinvariants.py b/mtpy/analysis/zinvariants.py index eace6360f..ef0eb4859 100644 --- a/mtpy/analysis/zinvariants.py +++ b/mtpy/analysis/zinvariants.py @@ -151,7 +151,6 @@ def compute_invariants(self): **q** : dependent variable suggesting dimensionality """ - print("computing invariants") # get the length of z to initialize some empty arrays nz = self.z.shape[0] @@ -184,8 +183,8 @@ def compute_invariants(self): ex = x1 * e1 - x2 * e2 - x3 * e3 + x4 * e4 if ex == 0.0: - print('Could not compute invariants for {0:5e} Hz'.format( - self.freq[ii])) + # print('Could not compute invariants for {0:5e} Hz'.format( + # self.freq[ii])) self.inv1[ii] = np.nan self.inv2[ii] = np.nan self.inv3[ii] = np.nan diff --git a/mtpy/imaging/phase_tensor_maps.py b/mtpy/imaging/phase_tensor_maps.py index c48ea67fa..34f719ed4 100644 --- a/mtpy/imaging/phase_tensor_maps.py +++ b/mtpy/imaging/phase_tensor_maps.py @@ -805,7 +805,12 @@ def plot(self, fig=None, save_path=None, show=True, colorarray = pt.ellipticity[jj] elif self.ellipse_colorby in ['strike', 'azimuth']: - colorarray = self.fold_strike(self.pt.azimuth) + colorarray = pt.azimuth[jj] % 180 + if colorarray > 90: + colorarray -= 180 + self.ellipse_range = (-90, 90) + ckmin = float(self.ellipse_range[0]) + ckmax = float(self.ellipse_range[1]) else: raise NameError(self.ellipse_colorby + ' is not supported') diff --git a/mtpy/imaging/plot_mt_response.py b/mtpy/imaging/plot_mt_response.py index 53dc99c3a..f28e0df7c 100644 --- a/mtpy/imaging/plot_mt_response.py +++ b/mtpy/imaging/plot_mt_response.py @@ -476,7 +476,9 @@ def plot(self, show=True, overlay_mt_obj=None): 'ellipticity': r'Ellipticity', 'skew_seg': r'Skew (deg)', 'normalized_skew_seg': r'Normalized Skew (deg)', - 'geometric_mean': r'$\sqrt{\Phi_{min} \cdot \Phi_{max}}$'} + 'geometric_mean': r'$\sqrt{\Phi_{min} \cdot \Phi_{max}}$', + 'strike': r"Strike (deg)", + 'azimuth': r"Strike (deg)"} if self.plot_tipper.find('y') == 0: if self.Tipper is None or np.all(self.Tipper.tipper == 0 + 0j): @@ -542,6 +544,7 @@ def plot(self, show=True, overlay_mt_obj=None): # make figure instance self.fig = plt.figure(self.fig_num, self.fig_size, dpi=self.fig_dpi) + self.fig.clf() # --> make figure for xy,yx components if self.plot_num == 1 or self.plot_num == 3: @@ -762,15 +765,15 @@ def plot(self, show=True, overlay_mt_obj=None): # -----plot tipper---------------------------------------------------- if self.plot_tipper.find('y') == 0: - txr = self.Tipper.mag_real * np.sin(self.Tipper.angle_real * np.pi / 180 + \ - np.pi * self.arrow_direction) - tyr = self.Tipper.mag_real * np.cos(self.Tipper.angle_real * np.pi / 180 + \ - np.pi * self.arrow_direction) + txr = self.Tipper.mag_real * \ + np.cos(np.deg2rad(self.Tipper.angle_real)) + tyr = self.Tipper.mag_real * \ + np.sin(np.deg2rad(self.Tipper.angle_real)) - txi = self.Tipper.mag_imag * np.sin(self.Tipper.angle_imag * np.pi / 180 + \ - np.pi * self.arrow_direction) - tyi = self.Tipper.mag_imag * np.cos(self.Tipper.angle_imag * np.pi / 180 + \ - np.pi * self.arrow_direction) + txi = self.Tipper.mag_imag * \ + np.cos(np.deg2rad(self.Tipper.angle_imag)) + tyi = self.Tipper.mag_imag * \ + np.sin(np.deg2rad(self.Tipper.angle_imag)) nt = len(txr) @@ -895,23 +898,27 @@ def plot(self, show=True, overlay_mt_obj=None): # get the properties to color the ellipses by if self.ellipse_colorby == 'phiminang' or \ self.ellipse_colorby == 'phimin': - colorarray = self.pt.phimin + color_array = self.pt.phimin elif self.ellipse_colorby == 'phimaxang' or \ self.ellipse_colorby == 'phimax': - colorarray = self.pt.phimax + color_array = self.pt.phimax elif self.ellipse_colorby == 'phidet': - colorarray = np.sqrt(abs(self.pt.det)) * (180 / np.pi) + color_array = np.sqrt(abs(self.pt.det)) * (180 / np.pi) elif self.ellipse_colorby == 'skew' or \ self.ellipse_colorby == 'skew_seg': - colorarray = self.pt.beta + color_array = self.pt.beta elif self.ellipse_colorby == 'ellipticity': - colorarray = self.pt.ellipticity + color_array = self.pt.ellipticity + + elif self.ellipse_colorby in ['strike', 'azimuth']: + color_array = self.pt.azimuth % 180 + color_array[np.where(color_array > 90)] -= 180 else: raise NameError(self.ellipse_colorby + ' is not supported') @@ -937,14 +944,14 @@ def plot(self, show=True, overlay_mt_obj=None): # get ellipse color if cmap.find('seg') > 0: - ellipd.set_facecolor(mtcl.get_plot_color(colorarray[ii], + ellipd.set_facecolor(mtcl.get_plot_color(color_array[ii], self.ellipse_colorby, cmap, ckmin, ckmax, bounds=bounds)) else: - ellipd.set_facecolor(mtcl.get_plot_color(colorarray[ii], + ellipd.set_facecolor(mtcl.get_plot_color(color_array[ii], self.ellipse_colorby, cmap, ckmin, @@ -1364,7 +1371,7 @@ def redraw_plot(self): >>> p1.redraw_plot() """ - plt.close(self.fig) + self.fig.clf() self.plot() def __str__(self): diff --git a/mtpy/imaging/plotresponse.py b/mtpy/imaging/plotresponse.py index 318101496..415c859dc 100644 --- a/mtpy/imaging/plotresponse.py +++ b/mtpy/imaging/plotresponse.py @@ -615,7 +615,7 @@ def plot(self): #phase_yx: Note add 180 to place it in same quadrant as phase_xy self.ebyxp = mtpl.plot_errorbar(self.axp, self.mt.period, - self.mt.Z.phase_yx+180*self.phase_quadrant, + self.mt.Z.phase_yx + 180*self.phase_quadrant, marker=self.yx_marker, ms=self.marker_size, color=self.yx_color, @@ -626,14 +626,14 @@ def plot(self): #check the phase to see if any point are outside of [0:90] if self.phase_limits == None: - if min(self.mt.Z.phase_xy)<0 or min(self.mt.Z.phase_yx)<0: + if min(self.mt.Z.phase_xy) < 0 or min(self.mt.Z.phase_yx) < 0: pymin = min([min(self.mt.Z.phase_xy), min(self.mt.Z.phase_yx)]) if pymin > 0: pymin = 0 else: pymin = 0 - if max(self.mt.Z.phase_xy)>90 or max(self.mt.Z.phase_yx)>90: + if max(self.mt.Z.phase_xy) > 90 or max(self.mt.Z.phase_yx) > 90: pymax = min([max(self.mt.Z.phase_xy), max(self.mt.Z.phase_yx)]) if pymax < 91: pymax = 89.9 @@ -817,8 +817,8 @@ def plot(self): s2 = self.mt.pt.azimuth s2_err = self.mt.pt.azimuth_err #fold angles to go from -90 to 90 - s2[np.where(s2>90)] -= 180 - s2[np.where(s2<-90)] += 180 + s2[np.where(s2 > 90)] -= 180 + s2[np.where(s2 < -90)] += 180 #plot strike with error bars ps2 = mtpl.plot_errorbar(self.axst, @@ -839,7 +839,7 @@ def plot(self): if self._plot_strike.find('t') > 0: #strike from tipper - s3 = self.mt.Tipper.angle_real+90 + s3 = self.mt.Tipper.angle_real + 90 #fold to go from -90 to 90 s3[np.where(s3 > 90)] -= 180 diff --git a/mtpy/imaging/plotstrike.py b/mtpy/imaging/plotstrike.py index 4a13038a8..d16c9b62f 100644 --- a/mtpy/imaging/plotstrike.py +++ b/mtpy/imaging/plotstrike.py @@ -304,7 +304,8 @@ def make_strike_array(self): # on the y tipr = 90 - tip.angle_real - tipr[np.where(tipr == 180.)] = 0.0 + tipr[np.where(abs(tipr) == 180.)] = np.nan + tipr[np.where(abs(tipr) == 0)] = np.nan # fold so the angle goes from 0 to 180 if self.fold == True: From 73429f205e23643af31778d5a994669321f2f071 Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 11 May 2020 18:08:57 -0700 Subject: [PATCH 10/20] added documentation to utils.calculator --- mtpy/utils/calculator.py | 53 +++++++++++++++++++++++++++++++--------- 1 file changed, 42 insertions(+), 11 deletions(-) diff --git a/mtpy/utils/calculator.py b/mtpy/utils/calculator.py index 7528bc6c5..b59bf7c27 100644 --- a/mtpy/utils/calculator.py +++ b/mtpy/utils/calculator.py @@ -474,16 +474,24 @@ def rotate_matrix_with_errors(in_matrix, angle, error=None): Rotate a matrix including errors clockwise given an angle in degrees. + :param in_matrix: A n x 2 x 2 matrix to rotate + :type inmatrix: np.ndarray - :param inmatrix: DESCRIPTION - :type inmatrix: TYPE - :param angle: DESCRIPTION - :type angle: TYPE - :param inmatrix_err: DESCRIPTION, defaults to None - :type inmatrix_err: TYPE, optional - :raises MTex: DESCRIPTION - :return: DESCRIPTION - :rtype: TYPE + :param angle: Angle to rotate by assuming clockwise positive from + 0 = north + :type angle: float + + :param error: A n x 2 x 2 matrix of associated errors, + defaults to None + :type error: np.ndarray, optional + + :raises MTex: If input array is incorrect + + :return: rotated matrix + :rtype: np.ndarray + + :return: rotated matrix errors + :rtype: np.ndarray """ @@ -534,14 +542,37 @@ def rotate_matrix_with_errors(in_matrix, angle, error=None): def rotate_vector_with_errors(in_vector, angle, error=None): - #check for row or column vector + """ + + Rotate a vector including errors clockwise given an angle in degrees. + + :param in_matrix: A n x 1 x 2 vector to rotate + :type invector: np.ndarray + + :param angle: Angle to rotate by assuming clockwise positive from + 0 = north + :type angle: float + + :param error: A n x 1 x 2 vector of associated errors, + defaults to None + :type error: np.ndarray, optional + + :raises MTex: If input array is incorrect + + :return: rotated vector + :rtype: np.ndarray + + :return: rotated vector errors + :rtype: np.ndarray + + """ if in_vector is None: raise MTex.MTpyError_inputarguments('Vector AND error-vector must'+ ' be defined') if (error is not None) and (in_vector.shape != error.shape): - msg = 'matricies are not the same shape in_matrix={0}, err={1}'.format( + msg = 'matricies are not the same shape in_vector={0}, err={1}'.format( in_vector.shape, error.shape) raise MTex.MTpyError_inputarguments(msg) From dbed2d0740ff161affd1434cf36f67b0ef0ef8e3 Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 11 May 2020 18:18:10 -0700 Subject: [PATCH 11/20] fixed plotnresponses --- mtpy/imaging/plotnresponses.py | 185 +++------------------------------ 1 file changed, 16 insertions(+), 169 deletions(-) diff --git a/mtpy/imaging/plotnresponses.py b/mtpy/imaging/plotnresponses.py index 98319bdc0..d456da4db 100644 --- a/mtpy/imaging/plotnresponses.py +++ b/mtpy/imaging/plotnresponses.py @@ -309,23 +309,6 @@ def __init__(self, **kwargs): self.plot_style = kwargs.pop('plot_style', '1') self.plot_title = kwargs.pop('plot_title', None) - # if rotation angle is an int or float make an array the length of - # mt_list for plotting purposes - 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 - - self._set_rot_z(self._rot_z) - # set plot limits self.xlimits = kwargs.pop('xlimits', None) self.res_limits = kwargs.pop('res_limits', None) @@ -385,34 +368,20 @@ def __init__(self, **kwargs): if self.plot_yn == 'y': self.plot() - # ---rotate data on setting rot_z - def _set_rot_z(self, rot_z): + #---need to rotate data on setting rotz + @property + def rotation_angle(self): + return self._rotation_angle + + @rotation_angle.setter + def rotation_angle(self, value): """ - need to rotate data when setting z + only a single value is allowed """ - - # 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)""") + mt.rotation_angle = value + + self._rotation_angle = value # --> on setting plot_ make sure to update the order and list def _set_plot_tipper(self, plot_tipper): @@ -501,86 +470,6 @@ def plot(self, show=True): # set height ratios of the subplots hr = [2, 1.5] + [1] * (len(list(pdict.keys())) - 2) - # if self.plot_style == '1': - # self.plotlist = [] - # - # #--> plot from edi's if given, don't need to rotate because - # # data has already been rotated by the funcion _set_rot_z - ## if self.fig_size is None: - ## self.fig_size = [6, 6] - # for ii, mt in enumerate(self.mt_list, 1): - # p1 = plotresponse(mt_object=mt, - # fig_num=ii, - # fig_size=self.fig_size, - # plot_num=self.plot_num, - # fig_dpi=self.fig_dpi, - # plot_yn='n', - # plot_tipper=self._plot_tipper, - # plot_strike=self._plot_strike, - # plot_skew=self._plot_skew, - # plot_pt=self._plot_pt) - # - # #make sure all the properties are set to match the users - # #line style between points - # p1.xy_ls = self.xy_ls - # p1.yx_ls = self.yx_ls - # p1.det_ls = self.det_ls - # - # #outline color - # p1.xy_color = self.xy_color - # p1.yx_color = self.yx_color - # p1.det_color = self.det_color - # - # #face color - # p1.xy_mfc = self.xy_mfc - # p1.yx_mfc = self.yx_mfc - # p1.det_mfc = self.det_mfc - # - # #maker - # p1.xy_marker = self.xy_marker - # p1.yx_marker = self.yx_marker - # p1.det_marker = self.det_marker - # - # #size - # p1.marker_size = 2 - # - # #set plot limits - # p1.xlimits = self.xlimits - # p1.res_limits = self.res_limits - # p1.phase_limits = self.phase_limits - # - # #set font parameters - # p1.font_size = self.font_size - # - # #set arrow properties - # p1.arrow_lw = self.arrow_lw - # p1.arrow_head_width = self.arrow_head_width - # p1.arrow_head_length = self.arrow_head_length - # p1.arrow_color_real = self.arrow_color_real - # p1.arrow_color_imag = self.arrow_color_imag - # p1.arrow_direction = self.arrow_direction - # p1.tipper_limits = self.tipper_limits - # - # #skew properties - # p1.skew_color = self.skew_color - # p1.skew_marker = self.skew_marker - # - # #strike properties - # p1.strike_inv_marker = self.strike_inv_marker - # p1.strike_inv_color = self.strike_inv_color - # - # p1.strike_pt_marker = self.strike_pt_marker - # p1.strike_pt_color = self.strike_pt_color - # - # p1.strike_tip_marker = self.strike_tip_marker - # p1.strike_tip_color = self.strike_tip_color - # - # #--> plot the apparent resistivity and phase - # self.plotlist.append(p1) - # - # p1.plot() - # - # -----Plot All in one figure with each plot as a subfigure------------ if self.plot_style == 'all': @@ -802,28 +691,6 @@ def plot(self, show=True): pymin = min(0, min([min(rp.phase_xy), min(rp.phase_yx)])) pymax = max(89.9, max([max(rp.phase_xy), max(rp.phase_yx)])) self.phase_limits = (pymin, pymax) - # self.phase_limits = (pymin, pymax) - # else: - # self.phase_limits = (min(self.phase_limits[0], pymin), - # max(self.phase_limits[1], pymax)) - # if self.phase_limits is None: - # if min(rp.phasexy) < 0 or min(rp.phase_yx) < 0: - # pymin = min([min(rp.phase_xy), - # min(rp.phase_yx)]) - # if pymin > 0: - # pymin = 0 - # else: - # pymin = 0 - # - # if max(rp.phasexy) > 90 or max(rp.phase_yx) > 90: - # pymax = min([max(rp.phase_xy), # YG: should use max instead ?? - # max(rp.phase_yx)]) - # if pymax < 91: - # pymax = 89.9 # YG: why?? - # else: - # pymax = 89.9 - # - # self.phase_limits = (pymin, pymax) # --> set axes properties if ii == 0: @@ -1182,6 +1049,11 @@ def plot(self, show=True): elif self.ellipse_colorby == 'ellipticity': colorarray = pt.ellipticity + elif self.ellipse_colorby in ['strike', 'azimuth']: + colorarray = self.fold_strike(pt.azimuth) + self.ellipse_range = (-90, 90) + ckmin = self.ellipse_range[0] + ckmax = self.ellipse_range[1] else: raise NameError(self.ellipse_colorby + ' is not supported') @@ -1888,31 +1760,6 @@ def plot(self, show=True): # ------plot strike angles---------------------------------------------- if self._plot_strike.find('y') == 0: - # if self._plot_strike.find('i') > 0: - # #strike from invariants - # zinv = mt.Z.invariants - # s1 = zinv.strike - # - # #fold angles so go from -90 to 90 - # s1[np.where(s1>90)] -= -180 - # s1[np.where(s1<-90)] += 180 - # - # #plot strike with error bars - # ps1 = self.axst.errorbar(mt.period, - # s1, - # marker=mxy[ii % len(mxy)], - # ms=self.marker_size, - # mfc=cst[ii], - # mec=cst[ii], - # mew=self.marker_lw, - # ls='none', - # yerr=zinv.strike_err, - # ecolor=cst[ii], - # capsize=self.marker_size, - # elinewidth=self.marker_lw) - # - # stlist.append(ps1[0]) - if self._plot_strike.find('p') > 0: # strike from phase tensor s2 = mt.pt.azimuth From 0287d7015274842611d8ebad2c5b71bc763ee7c6 Mon Sep 17 00:00:00 2001 From: JP Date: Thu, 21 May 2020 10:00:21 -0700 Subject: [PATCH 12/20] fixed z_err.setter issue, set z_err as a property --- mtpy/analysis/pt.py | 1 + 1 file changed, 1 insertion(+) diff --git a/mtpy/analysis/pt.py b/mtpy/analysis/pt.py index 7e8a58cfd..776fa7257 100644 --- a/mtpy/analysis/pt.py +++ b/mtpy/analysis/pt.py @@ -322,6 +322,7 @@ def z(self, z_array): # ---Z Error array--------------------------------------------------------------- + @property def z_err(self): return self._z_err From 17238b7bd9675adfa91bd171548090a377d90bc5 Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 25 May 2020 22:57:39 -0700 Subject: [PATCH 13/20] refactored plot_pt_maps --- examples/scripts/plot_phase_tensor_map.py | 14 +- mtpy/imaging/phase_tensor_maps.py | 862 +++++++++++----------- 2 files changed, 418 insertions(+), 458 deletions(-) diff --git a/examples/scripts/plot_phase_tensor_map.py b/examples/scripts/plot_phase_tensor_map.py index 094917d14..92b88ec76 100644 --- a/examples/scripts/plot_phase_tensor_map.py +++ b/examples/scripts/plot_phase_tensor_map.py @@ -15,11 +15,10 @@ import mtpy.imaging.phase_tensor_maps as pptmaps # directory containing edis -#edipath = r'C:\mtpywin\mtpy\examples\data\edi2' -edipath = '/home/bren/data_mtpy/phase_tensor_map/MT086_Edited_EDIs' +edipath = r'C:\mtpywin\mtpy\examples\data\edi2' # whether or not to save the figure to file -save = True +save = False # full path to file to save to savepath = '/tmp' @@ -28,8 +27,9 @@ plot_freq = 1e-2 # value to color ellipses by, options are phimin,phimax,skew -colorby = 'skew' +colorby = 'strike' ellipse_range = [-9, 9] +rotation_angle = 0 image_fn = 'phase_tensor_map%1is_' % (int(1. / plot_freq)) + colorby + '.png' @@ -57,11 +57,11 @@ arrow_head_width=0.002, arrow_head_length=0.002, #ellipse_cmap='mt_seg_bl2wh2rd' - station_dict={'id': (5, 7)}, - background_image='/home/bren/data_mtpy/phase_tensor_map/gravity_Sept2017_Complete_Bouguer_focus_area1_10000_residual_he_ness.tif' + station_id=(5, 7), + #background_image='/home/bren/data_mtpy/phase_tensor_map/gravity_Sept2017_Complete_Bouguer_focus_area1_10000_residual_he_ness.tif' # bimg_band=1, # Optional, set to None by default #bimg_cmap='viridis' # Optional, set to 'viridis' by default - ) + rotation_angle=rotation_angle) if save: diff --git a/mtpy/imaging/phase_tensor_maps.py b/mtpy/imaging/phase_tensor_maps.py index 34f719ed4..57cfab6b8 100644 --- a/mtpy/imaging/phase_tensor_maps.py +++ b/mtpy/imaging/phase_tensor_maps.py @@ -367,6 +367,10 @@ def __init__(self, **kwargs): :param kwargs: keyword-value pairs """ super(PlotPhaseTensorMaps, self).__init__(**kwargs) + + self.fig = None + self.ax_ptm = None + self.ax_cb = None fn_list = kwargs.pop('fn_list', None) z_object_list = kwargs.pop('z_object_list', None) @@ -382,141 +386,103 @@ def __init__(self, **kwargs): mt_object_list=mt_object_list) # set the freq to plot - self.plot_freq = kwargs.pop('plot_freq', 1.0) - self.ftol = kwargs.pop('ftol', 0.1) - self.interpolate = kwargs.pop('interpolate', True) + self.plot_freq = 1.0 + self.ftol = 0.1 + self.interpolate = True # read in map scale - self.mapscale = kwargs.pop('mapscale', 'deg') + self.mapscale = 'deg' # map background image - self.background_image = kwargs.pop('background_image', None) - self.bimg_band = kwargs.pop('bimg_band', None) - self.bimg_cmap = kwargs.pop('bimg_cmap', 'viridis') + self.background_image = None + self.bimg_band = None + self.bimg_cmap = 'viridis' # --> set the ellipse properties ------------------- # set default size to 2 if self.mapscale == 'deg': - self.ellipse_size = kwargs.pop('ellipse_size', .05) - self.arrow_size = kwargs.pop('arrow_size', .05) - self.arrow_head_length = kwargs.pop('arrow_head_length', .005) - self.arrow_head_width = kwargs.pop('arrow_head_width', .005) - self.arrow_lw = kwargs.pop('arrow_lw', .0005) - self.xpad = kwargs.pop('xpad', .05) - self.ypad = kwargs.pop('xpad', .05) + self.ellipse_size = .05 + self.arrow_size = .05 + self.arrow_head_length = .005 + self.arrow_head_width = .005 + self.arrow_lw = .0005 + self.xpad = .05 + self.ypad = .05 elif self.mapscale == 'm': - self.ellipse_size = kwargs.pop('ellipse_size', 500) - self.arrow_size = kwargs.pop('arrow_size', 500) - self.arrow_head_length = kwargs.pop('arrow_head_length', 50) - self.arrow_head_width = kwargs.pop('arrow_head_width', 50) - self.arrow_lw = kwargs.pop('arrow_lw', 5) - self.xpad = kwargs.pop('xpad', 500) - self.ypad = kwargs.pop('xpad', 500) + self.ellipse_size = 500 + self.arrow_size = 500 + self.arrow_head_length = 50 + self.arrow_head_width = 50 + self.arrow_lw = 5 + self.xpad = 500 + self.ypad = 500 elif self.mapscale == 'km': - self.ellipse_size = kwargs.pop('ellipse_size', .5) - self.arrow_size = kwargs.pop('arrow_size', .5) - self.arrow_head_length = kwargs.pop('arrow_head_length', .05) - self.arrow_head_width = kwargs.pop('arrow_head_width', .05) - self.arrow_lw = kwargs.pop('arrow_lw', .005) - self.xpad = kwargs.pop('xpad', .5) - self.ypad = kwargs.pop('xpad', .5) + self.ellipse_size = .5 + self.arrow_size = .5 + self.arrow_head_length = .05 + self.arrow_head_width = .05 + self.arrow_lw = .005 + self.xpad = .5 + self.ypad = .5 - self.minorticks_on = kwargs.pop('minorticks_on', True) + self.minorticks_on = True # --> set colorbar properties--------------------------------- # set orientation to horizontal - cb_dict = kwargs.pop('cb_dict', {}) - try: - self.cb_orientation = cb_dict['orientation'] - except KeyError: - self.cb_orientation = 'vertical' - - # set the position to middle outside the plot - try: - self.cb_position = cb_dict['position'] - except KeyError: - self.cb_position = None + self.cb_orientation = 'vertical' + self.cb_position = None # --> set plot properties ------------------------------ # set some of the properties as attributes much to Lars' discontent - self.fig_num = kwargs.pop('fig_num', 1) - self.plot_num = kwargs.pop('plot_num', 1) - self.plot_style = kwargs.pop('plot_style', '1') - self.plot_title = kwargs.pop('plot_title', None) - self.fig_dpi = kwargs.pop('fig_dpi', 300) + self.fig_num = 1 + self.plot_title = 'Phase Tensor Map for ' + self.fig_dpi = 300 - self.tscale = kwargs.pop('tscale', 'period') - self.fig_size = kwargs.pop('fig_size', [8, 8]) - self.tickstrfmt = kwargs.pop('tickstrfmt', '%.2f') + self.tscale = 'period' + self.fig_size = None + self.tickstrfmt = '%.2f' - self.font_size = kwargs.pop('font_size', 7) + self.font_size = 7 - self.ref_ax_loc = kwargs.pop('ref_ax_loc', (.85, .1, .1, .1)) + self.ref_ax_loc = (.85, .1, .1, .1) # if rotation angle is an int or float make an array the length of # mt_list for plotting purposes - 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 induction arrow properties ------------------------------- - self.plot_tipper = kwargs.pop('plot_tipper', 'n') + self.plot_tipper = 'n' # --> set arrow legend properties ------------------------------- - arrow_legend_dict = kwargs.pop('arrow_legend_dict', {}) - self.arrow_legend_position = arrow_legend_dict.pop('position', 'lower right') - - # set x-border pad - self.arrow_legend_xborderpad = arrow_legend_dict.pop('xborderpad', .2) - - # set y-border pad - self.arrow_legend_yborderpad = arrow_legend_dict.pop('yborderpad', .2) - - # set font pad - self.arrow_legend_fontpad = arrow_legend_dict.pop('fontpad', .05) - - # set font properties - self.arrow_legend_fontdict = arrow_legend_dict.pop('fontdict', {'size': self.font_size, - 'weight': 'bold'}) + self.arrow_legend_position = 'lower right' + self.arrow_legend_xborderpad = .2 + self.arrow_legend_yborderpad = .2 + self.arrow_legend_fontpad = .05 + self.arrow_legend_fontdict = {'size': self.font_size, + 'weight': 'bold'} # --> set a central reference point - self.plot_reference_point = kwargs.pop('reference_point', (0, 0)) + self.plot_reference_point = (0, 0) # --> set station name properties - station_dict = kwargs.pop('station_dict', None) - if station_dict is not None: - self.station_id = station_dict.pop('id', (0, 2)) + self.plot_station = False + self.station_id = (0, 2) + - # set spacing of station name and ellipse - self.station_pad = station_dict.pop('pad', .0005) + # set spacing of station name and ellipse + self.station_pad = .0005 - # set font properties of the station label - self.station_font_dict = station_dict.pop('font_dict', {'size': self.font_size, - 'weight': 'bold'}) + # set font properties of the station label + self.station_font_dict = {'size': self.font_size, 'weight': 'bold'} - self.plot_yn = kwargs.pop('plot_yn', 'y') - self.save_fn = kwargs.pop('save_fn', "/c/tmp") + self.plot_yn = 'y' + self.save_fn = "/c/tmp" - # By this stage all keyword arguments meant to be set as class properties will have - # been processed. Popping all class properties that still exist in kwargs - self.kwargs = kwargs - for key in vars(self): - self.kwargs.pop(key, None) + for key, value in kwargs.items(): + setattr(self, key, value) # --> plot if desired ------------------------ if self.plot_yn == 'y': self.plot() - # self.save_figure(self.save_fn, file_format='png') #---need to rotate data on setting rotz @property @@ -528,8 +494,10 @@ def rotation_angle(self, value): """ only a single value is allowed """ - for ii, mt in enumerate(self.mt_list): - mt.rotation_angle = value + if value == 0: + return + for ii, mt_obj in enumerate(self.mt_list): + mt_obj.rotation_angle = value self._rotation_angle = value @@ -541,61 +509,11 @@ def fold_strike(self, strike): strike[np.where(strike > 90)] -= 180 return strike - - # ----------------------------------------------- - # The main plot method for this module - # ----------------------------------------------- - def plot(self, fig=None, save_path=None, show=True, - raster_dict={'lons': [], 'lats': [], - 'vals': [], 'levels': 50, 'cmap': 'rainbow', - 'cbar_title': 'Arbitrary units', - 'cbar_position': None}): + + def add_raster(self, raster_dict): """ - Plots the phase tensor map. - :param fig: optional figure object - :param save_path: path to folder for saving plots - :param show: show plots if True - :param raster_dict: Plotting of raster data is currently only supported when mapscale='deg'. - This parameter is a dictionary of parameters for plotting raster data, - on top of which phase tensor data are plotted. 'lons', 'lats' and 'vals' - are one dimensional lists (or numpy arrays) for longitudes, latitudes - and corresponding values, respectively. 'levels', 'cmap' and 'cbar_title' - are the number of levels to be used in the colormap, the colormap and its - title, respectively. """ - - # set position properties for the plot - plt.rcParams['font.size'] = self.font_size - plt.rcParams['figure.subplot.left'] = .1 - plt.rcParams['figure.subplot.right'] = .98 - plt.rcParams['figure.subplot.bottom'] = .1 - plt.rcParams['figure.subplot.top'] = .93 - plt.rcParams['figure.subplot.wspace'] = .55 - plt.rcParams['figure.subplot.hspace'] = .70 - - lpfig = None - lpax = None - lpax2 = None - # make figure instance - if(fig is None): - self.fig = plt.figure(self.fig_num, figsize=self.fig_size, dpi=self.fig_dpi) - # self.fig = plt.figure(self.fig_num, dpi=self.fig_dpi) - - # clear the figure if there is already one up - plt.clf() - lpfig = self.fig - else: - lpfig = fig - # end if - - lpax = lpfig.add_subplot(1, 1, 1, aspect='equal') - - # plt.locator_params(axis='x', nbins=3) # control number of ticks in axis (nbins ticks) - plt.xticks(rotation='vertical') # FZ: control tick rotation=30 not that good - - # get the reference point - refpoint = self.plot_reference_point - + # plot raster data if provided and if mapscale is 'deg' if(len(raster_dict['lons']) and self.mapscale == 'deg'): lons = np.array(raster_dict['lons']) @@ -607,24 +525,28 @@ def plot(self, fig=None, save_path=None, show=True, else: vals = np.array(raster_dict['vals']) - assert len(lons) == len(lats) == len(vals), 'Lons, Lats and Vals must all have the same length' + if not len(lons) == len(lats) == len(vals): + msg= 'Raster Lons, Lats and Vals must all have the same length' + _logger.error(msg) - lons -= refpoint[0] - lats -= refpoint[1] + lons -= self.plot_reference_point[0] + lats -= self.plot_reference_point[1] levels = raster_dict.pop('levels', 50) cmap = raster_dict.pop('cmap', 'rainbow') cbar_title = raster_dict.pop('cbar_title', 'Arbitrary Units') triangulation = tri.Triangulation(lons, lats) - cbinfo = lpax.tricontourf(triangulation, vals, - levels=np.linspace(vals.min(), vals.max(), levels), + cbinfo = self.ax_ptm.tricontourf(triangulation, vals, + levels=np.linspace(vals.min(), + vals.max(), + levels), cmap=cmap) if raster_dict['cbar_position'] is not None: cbax = self.fig.add_axes(raster_dict['cbar_position']) else: - cbax, kw = mcb.make_axes(lpax, + cbax, kw = mcb.make_axes(self.ax_ptm, orientation=self.cb_orientation, shrink=.35) - cbar = lpfig.colorbar(cbinfo, cbax) + cbar = self.fig.colorbar(cbinfo, cbax) if(self.cb_orientation == 'horizontal'): cbar.ax.set_xlabel(cbar_title) @@ -639,25 +561,257 @@ def plot(self, fig=None, save_path=None, show=True, cbar.ax.tick_params(axis='y', direction='in') # end if + + def add_tipper(self, tipper_obj, plot_x, plot_y): + """ + + :param tipper: DESCRIPTION + :type tipper: TYPE + :return: DESCRIPTION + :rtype: TYPE + + """ + + # make some local parameters for easier typing + ascale = self.arrow_size + adir = self.arrow_direction * np.pi + + # plot real tipper + ti = tipper_obj + if ti is None: + return + + if self.plot_tipper == 'yri' or self.plot_tipper == 'yr': + if ti.mag_real[self.jj] <= self.arrow_threshold: + txr = ti.mag_real[self.jj] * ascale * \ + np.sin((ti.angle_real[self.jj]) * np.pi / 180 + adir) + tyr = ti.mag_real[self.jj] * ascale * \ + np.cos((ti.angle_real[self.jj]) * np.pi / 180 + adir) + + self.ax_ptm.arrow(plot_x, + plot_y, + txr, + tyr, + width=self.arrow_lw, + facecolor=self.arrow_color_real, + edgecolor=self.arrow_color_real, + length_includes_head=False, + head_width=self.arrow_head_width, + head_length=self.arrow_head_length) + else: + pass + + # plot imaginary tipper + if self.plot_tipper == 'yri' or self.plot_tipper == 'yi': + if ti.mag_imag[self.jj] <= self.arrow_threshold: + txi = ti.mag_imag[self.jj] * ascale * \ + np.sin((ti.angle_imag[self.jj]) * np.pi / 180 + adir) + tyi = ti.mag_imag[self.jj] * ascale * \ + np.cos((ti.angle_imag[self.jj]) * np.pi / 180 + adir) + + self.ax_ptm.arrow(plot_x, + plot_y, + txi, + tyi, + width=self.arrow_lw, + facecolor=self.arrow_color_imag, + edgecolor=self.arrow_color_imag, + length_includes_head=False, + head_width=self.arrow_head_width, + head_length=self.arrow_head_length) + + def get_xy(self, mt_obj): + """ + """ + # if map scale is lat lon set parameters + if self.mapscale == 'deg': + plot_x = mt_obj.lon - self.plot_reference_point[0] + plot_y = mt_obj.lat - self.plot_reference_point[1] + + # if map scale is in meters easting and northing + elif self.mapscale in ['km', 'm']: + if self.mapscale == 'km': + scale = 1./1000 + else: + scale = 1. + east, north, zone = gis_tools.project_point_ll2utm(mt_obj.lat, + mt_obj.lon) + + # set the first point read in as a refernce other points + if not hasattr(self._utm_zone_01): + self._utm_zone_01 = zone + plot_x = east - self.plot_reference_point[0] + plot_y = north - self.plot_reference_point[1] + + else: + # check to make sure the zone is the same this needs + # to be more rigorously done + if self._utm_zone_01 != zone: + print('Zone change at station ' + mt_obj.station) + if zone1[0:2] == zone[0:2]: + pass + elif int(zone1[0:2]) < int(zone[0:2]): + east += 500000 + else: + east -= -500000 + plot_x = east - self.plot_reference_point[0] + plot_y = north - self.plot_reference_point[1] + else: + plot_x = east - self.plot_reference_point[0] + plot_y = north - self.plot_reference_point[1] + + plot_x *= scale + plot_y *= scale + + else: + raise NameError('mapscale not recognized') + + return plot_x, plot_y + + def get_pt_ellipse(self, pt, plot_x, plot_y): + """ + """ + + # --> set local variables + phimin = np.nan_to_num(pt.phimin[self.jj]) + phimax = np.nan_to_num(pt.phimax[self.jj]) + eangle = np.nan_to_num(pt.azimuth[self.jj]) + + if self.ellipse_cmap == 'mt_seg_bl2wh2rd': + bounds = np.arange(self.ellipse_range[0], + self.ellipse_range[1] + self.ellipse_range[2], + self.ellipse_range[2]) + nseg = float((self.ellipse_range[1] - self.ellipse_range[0]) / \ + (2 * self.ellipse_range[2])) + + # get the properties to color the ellipses by + if self.ellipse_colorby == 'phiminang' or \ + self.ellipse_colorby == 'phimin': + colorarray = pt.phimin[self.jj] + + elif self.ellipse_colorby == 'phimax': + colorarray = pt.phimax[self.jj] + + elif self.ellipse_colorby == 'phidet': + colorarray = np.sqrt(abs(pt.det[self.jj])) * (180 / np.pi) + + elif self.ellipse_colorby == 'skew' or \ + self.ellipse_colorby == 'skew_seg': + colorarray = pt.beta[self.jj] + + elif self.ellipse_colorby == 'normalized_skew' or \ + self.ellipse_colorby == 'normalized_skew_seg': + colorarray = 2 * pt.beta[self.jj] + + elif self.ellipse_colorby == 'ellipticity': + colorarray = pt.ellipticity[self.jj] + + elif self.ellipse_colorby in ['strike', 'azimuth']: + colorarray = pt.azimuth[self.jj] % 180 + if colorarray > 90: + colorarray -= 180 + self.ellipse_range = (-90, 90) + + else: + raise NameError(self.ellipse_colorby + ' is not supported') + + # --> get ellipse properties + # if the ellipse size is not physically correct make it a dot + if phimax == 0 or phimax > 100 or phimin == 0 or phimin > 100: + eheight = .0000001 * self.ellipse_size + ewidth = .0000001 * self.ellipse_size + else: + scaling = self.ellipse_size / phimax + eheight = phimin * scaling + ewidth = phimax * scaling + + # make an ellipse + ellipd = patches.Ellipse((plot_x, plot_y), + width=ewidth, + height=eheight, + angle=90 - eangle, + lw=self.lw) + + # get ellipse color + if self.ellipse_cmap.find('seg') > 0: + ellipd.set_facecolor(mtcl.get_plot_color(colorarray, + self.ellipse_colorby, + self.ellipse_cmap, + self.ellipse_range[0], + self.ellipse_range[1], + bounds=bounds)) + else: + ellipd.set_facecolor(mtcl.get_plot_color(colorarray, + self.ellipse_colorby, + self.ellipse_cmap, + self.ellipse_range[0], + self.ellipse_range[1])) + + return ellipd + + + # ----------------------------------------------- + # The main plot method for this module + # ----------------------------------------------- + def plot(self, fig=None, save_path=None, show=True, + raster_dict={'lons': [], 'lats': [], + 'vals': [], 'levels': 50, 'cmap': 'rainbow', + 'cbar_title': 'Arbitrary units', + 'cbar_position': None}): + """ + Plots the phase tensor map. + :param fig: optional figure object + :param save_path: path to folder for saving plots + :param show: show plots if True + :param raster_dict: Plotting of raster data is currently only + supported when mapscale='deg'. + This parameter is a dictionary of parameters for + plotting raster data, + on top of which phase tensor data are plotted. + 'lons', 'lats' and 'vals' are one dimensional + lists (or numpy arrays) for longitudes, latitudes + and corresponding values, respectively. 'levels', + 'cmap' and 'cbar_title' are the number of levels + to be used in the colormap, the colormap and its + title, respectively. + """ + + # set position properties for the plot + plt.rcParams['font.size'] = self.font_size + + # make figure instance + if fig is not None: + self.fig = fig + else: + self.fig = plt.figure(self.fig_num, figsize=self.fig_size, + dpi=self.fig_dpi) + plt.clf() + + self.ax_ptm = self.fig.add_subplot(1, 1, 1, aspect='equal') + + # plt.locator_params(axis='x', nbins=3) + # control number of ticks in axis (nbins ticks) + # FZ: control tick rotation=30 not that good + plt.xticks(rotation='vertical') + + self.add_raster(raster_dict) - # set some local parameters - es = float(self.ellipse_size) - cmap = self.ellipse_cmap - ckmin = float(self.ellipse_range[0]) - ckmax = float(self.ellipse_range[1]) try: - ckstep = float(self.ellipse_range[2]) + step = float(self.ellipse_range[2]) except IndexError: - if cmap == 'mt_seg_bl2wh2rd': + if self.ellipse_cmap == 'mt_seg_bl2wh2rd': raise ValueError('Need to input range as (min, max, step)') else: - ckstep = 3 - nseg = float((ckmax - ckmin) / (2 * ckstep)) + step = 3 + nseg = float((self.ellipse_range[1] - self.ellipse_range[0]) / \ + (2 * step)) ck = self.ellipse_colorby # --> set the bounds on the segmented colormap - if cmap == 'mt_seg_bl2wh2rd': - bounds = np.arange(ckmin, ckmax + ckstep, ckstep) + if self.ellipse_cmap == 'mt_seg_bl2wh2rd': + bounds = np.arange(self.ellipse_range[0], + self.ellipse_range[1] + step, + step) # set tick parameters depending on the mapscale if self.mapscale == 'deg': @@ -668,310 +822,119 @@ def plot(self, fig=None, save_path=None, show=True, # make some empty arrays elliplist = [] - latlist = np.zeros(len(self.mt_list)) - lonlist = np.zeros(len(self.mt_list)) self.plot_xarr = np.zeros(len(self.mt_list)) self.plot_yarr = np.zeros(len(self.mt_list)) - for ii, mt in enumerate(self.mt_list): + for ii, mt_obj in enumerate(self.mt_list): - newZ = None - newTipper = None + new_z_obj = None + new_tipper_obj = None fidx = 0 - if(self.interpolate): - newZ, newTipper = mt.interpolate([self.plot_freq], bounds_error=False) + if self.interpolate: + new_z_obj, new_tipper_obj = mt_obj.interpolate([self.plot_freq], + bounds_error=False) else: - fidx = np.argmin(np.fabs(mt.Z.freq - self.plot_freq)) + fidx = np.argmin(np.fabs(mt_obj.Z.freq - self.plot_freq)) - if((not self.interpolate and np.fabs(mt.Z.freq[fidx] - self.plot_freq) < self.ftol) or - (self.interpolate)): + if((not self.interpolate and \ + np.fabs(mt_obj.Z.freq[fidx] - self.plot_freq) < self.ftol) or \ + (self.interpolate)): self.jj = fidx jj = fidx # get phase tensor if(not self.interpolate): - pt = mt.pt + pt = mt_obj.pt else: - newZ.compute_resistivity_phase() - pt = MTpt.PhaseTensor(z_object=newZ) - - # if map scale is lat lon set parameters - if self.mapscale == 'deg': - latlist[ii] = mt.lat - lonlist[ii] = mt.lon - plotx = mt.lon - refpoint[0] - ploty = mt.lat - refpoint[1] - - # if map scale is in meters easting and northing - elif self.mapscale == 'm': - east, north, zone = gis_tools.project_point_ll2utm(mt.lat, - mt.lon) - - # set the first point read in as a refernce other points - if ii == 0: - zone1 = zone - plotx = east - refpoint[0] - ploty = north - refpoint[1] - - # read in all the other point - else: - # check to make sure the zone is the same this needs - # to be more rigorously done - if zone1 != zone: - print('Zone change at station ' + mt.station) - if zone1[0:2] == zone[0:2]: - pass - elif int(zone1[0:2]) < int(zone[0:2]): - east += 500000 - else: - east -= -500000 - latlist[ii] = north - refpoint[1] - lonlist[ii] = east - refpoint[0] - plotx = east - refpoint[0] - ploty = north - refpoint[1] - else: - latlist[ii] = north - refpoint[1] - lonlist[ii] = east - refpoint[0] - plotx = east - refpoint[0] - ploty = north - refpoint[1] - - # if mapscale is in km easting and northing - elif self.mapscale == 'km': - east, north, zone = gis_tools.project_point_ll2utm(mt.lat, - mt.lon) - if ii == 0: - zone1 = zone - plotx = (east - refpoint[0]) / 1000. - ploty = (north - refpoint[1]) / 1000. - - else: - if zone1 != zone: - print('Zone change at station ' + mt.station) - if zone1[0:2] == zone[0:2]: - pass - elif int(zone1[0:2]) < int(zone[0:2]): - east += 500000 - else: - east -= 500000 - latlist[ii] = (north - refpoint[1]) / 1000. - lonlist[ii] = (east - refpoint[0]) / 1000. - plotx = (east - refpoint[0]) / 1000. - ploty = (north - refpoint[1]) / 1000. - else: - latlist[ii] = (north - refpoint[1]) / 1000. - lonlist[ii] = (east - refpoint[0]) / 1000. - plotx = (east - refpoint[0]) / 1000. - ploty = (north - refpoint[1]) / 1000. - else: - raise NameError('mapscale not recognized') + new_z_obj.compute_resistivity_phase() + pt = MTpt.PhaseTensor(z_object=new_z_obj) + plot_x, plot_y = self.get_xy(mt_obj) # put the location of each ellipse into an array in x and y - self.plot_xarr[ii] = plotx - self.plot_yarr[ii] = ploty - - # --> set local variables - phimin = np.nan_to_num(pt.phimin[jj]) - phimax = np.nan_to_num(pt.phimax[jj]) - eangle = np.nan_to_num(pt.azimuth[jj]) - - # output to csv file: - # print('OUTCSV', mt.station, plotx, ploty, phimin, phimax, eangle) - - if cmap == 'mt_seg_bl2wh2rd': - bounds = np.arange(ckmin, ckmax + ckstep, ckstep) - nseg = float((ckmax - ckmin) / (2 * ckstep)) - - # get the properties to color the ellipses by - if self.ellipse_colorby == 'phiminang' or \ - self.ellipse_colorby == 'phimin': - colorarray = pt.phimin[jj] - - elif self.ellipse_colorby == 'phimax': - colorarray = pt.phimax[jj] - - elif self.ellipse_colorby == 'phidet': - colorarray = np.sqrt(abs(pt.det[jj])) * (180 / np.pi) - - elif self.ellipse_colorby == 'skew' or \ - self.ellipse_colorby == 'skew_seg': - colorarray = pt.beta[jj] - - elif self.ellipse_colorby == 'normalized_skew' or \ - self.ellipse_colorby == 'normalized_skew_seg': - colorarray = 2 * pt.beta[jj] - - elif self.ellipse_colorby == 'ellipticity': - colorarray = pt.ellipticity[jj] - - elif self.ellipse_colorby in ['strike', 'azimuth']: - colorarray = pt.azimuth[jj] % 180 - if colorarray > 90: - colorarray -= 180 - self.ellipse_range = (-90, 90) - ckmin = float(self.ellipse_range[0]) - ckmax = float(self.ellipse_range[1]) - - else: - raise NameError(self.ellipse_colorby + ' is not supported') - - # --> get ellipse properties - # if the ellipse size is not physically correct make it a dot - if phimax == 0 or phimax > 100 or phimin == 0 or phimin > 100: - eheight = .0000001 * es - ewidth = .0000001 * es - else: - scaling = es / phimax - eheight = phimin * scaling - ewidth = phimax * scaling - - # make an ellipse - ellipd = patches.Ellipse((plotx, ploty), - width=ewidth, - height=eheight, - angle=90 - eangle, - lw=self.lw, - **self.kwargs) - - # get ellipse color - if cmap.find('seg') > 0: - ellipd.set_facecolor(mtcl.get_plot_color(colorarray, - self.ellipse_colorby, - cmap, - ckmin, - ckmax, - bounds=bounds)) - else: - ellipd.set_facecolor(mtcl.get_plot_color(colorarray, - self.ellipse_colorby, - cmap, - ckmin, - ckmax)) + self.plot_xarr[ii] = plot_x + self.plot_yarr[ii] = plot_y + + ellipd = self.get_pt_ellipse(pt, plot_x, plot_y) # ==> add ellipse to the plot elliplist.append(ellipd) - lpax.add_artist(ellipd) + self.ax_ptm.add_artist(ellipd) # -----------Plot Induction Arrows--------------------------- if self.plot_tipper.find('y') == 0: - - # get tipper - - if mt.Tipper.tipper is None: - mt.Tipper.tipper = np.zeros((len(mt.period), 1, 2), - dtype='complex') - - # make some local parameters for easier typing - ascale = self.arrow_size - adir = self.arrow_direction * np.pi - - # plot real tipper - ti = None - if(not self.interpolate): - ti = mt.Tipper + # get phase tensor + if not self.interpolate: + new_tipper_obj = mt_obj.Tipper else: - ti = newTipper - # end if - - if self.plot_tipper == 'yri' or self.plot_tipper == 'yr': - if ti.mag_real[jj] <= self.arrow_threshold: - txr = ti.mag_real[jj] * ascale * \ - np.sin((ti.angle_real[jj]) * np.pi / 180 + adir) - tyr = ti.mag_real[jj] * ascale * \ - np.cos((ti.angle_real[jj]) * np.pi / 180 + adir) - - lpax.arrow(plotx, - ploty, - txr, - tyr, - width=self.arrow_lw, - facecolor=self.arrow_color_real, - edgecolor=self.arrow_color_real, - length_includes_head=False, - head_width=self.arrow_head_width, - head_length=self.arrow_head_length) - else: - pass - - # plot imaginary tipper - if self.plot_tipper == 'yri' or self.plot_tipper == 'yi': - if ti.mag_imag[jj] <= self.arrow_threshold: - txi = ti.mag_imag[jj] * ascale * \ - np.sin((ti.angle_imag[jj]) * np.pi / 180 + adir) - tyi = ti.mag_imag[jj] * ascale * \ - np.cos((ti.angle_imag[jj]) * np.pi / 180 + adir) - - lpax.arrow(plotx, - ploty, - txi, - tyi, - width=self.arrow_lw, - facecolor=self.arrow_color_imag, - edgecolor=self.arrow_color_imag, - length_includes_head=False, - head_width=self.arrow_head_width, - head_length=self.arrow_head_length) + new_tipper_obj.compute_mag_direction() + + if mt_obj.Tipper.tipper is None: + continue + self.add_tipper(new_tipper_obj, plot_x, plot_y) # ------------Plot station name------------------------------ - try: - lpax.text(plotx, - ploty + self.station_pad, - mt.station[self.station_id[0]:self.station_id[1]], - horizontalalignment='center', - verticalalignment='baseline', - fontdict=self.station_font_dict) - except AttributeError: - pass + if self.plot_station: + try: + name = mt_obj.station[self.station_id[0]:self.station_id[1]] + self.ax_ptm.text(plot_x, + plot_y + self.station_pad, + name, + horizontalalignment='center', + verticalalignment='baseline', + fontdict=self.station_font_dict) + except AttributeError: + pass # ==> print a message if couldn't find the freq else: - _logger.warn('Did not find {0:.5g} Hz for station {1}'.format(self.plot_freq, mt.station)) + _logger.warn('Did not find {0:.5g} Hz for station {1}'.format(self.plot_freq, mt_obj.station)) # --> set axes properties depending on map scale------------------------ if self.mapscale == 'deg': - lpax.set_xlabel('Longitude', + self.ax_ptm.set_xlabel('Longitude', fontsize=self.font_size, # +2, fontweight='bold') - lpax.set_ylabel('Latitude', + self.ax_ptm.set_ylabel('Latitude', fontsize=self.font_size, # +2, fontweight='bold') elif self.mapscale == 'm': - lpax.set_xlabel('Easting (m)', + self.ax_ptm.set_xlabel('Easting (m)', fontsize=self.font_size, # +2, fontweight='bold') - lpax.set_ylabel('Northing (m)', + self.ax_ptm.set_ylabel('Northing (m)', fontsize=self.font_size, # +2, fontweight='bold') elif self.mapscale == 'km': - lpax.set_xlabel('Easting (km)', + self.ax_ptm.set_xlabel('Easting (km)', fontsize=self.font_size, # +2, fontweight='bold') - lpax.set_ylabel('Northing (km)', + self.ax_ptm.set_ylabel('Northing (km)', fontsize=self.font_size, # +2, fontweight='bold') # --> set plot limits # need to exclude zero values from the calculation of min/max!!!! - lpax.set_xlim(self.plot_xarr[self.plot_xarr != 0.].min() - self.xpad, + self.ax_ptm.set_xlim(self.plot_xarr[self.plot_xarr != 0.].min() - self.xpad, self.plot_xarr[self.plot_xarr != 0.].max() + self.xpad) - lpax.set_ylim(self.plot_yarr[self.plot_yarr != 0.].min() - self.xpad, + self.ax_ptm.set_ylim(self.plot_yarr[self.plot_yarr != 0.].min() - self.xpad, self.plot_yarr[self.plot_xarr != 0.].max() + self.xpad) # BM: Now that we have the bounds of the axis, we can plot a # background image on the map. if self.background_image: - plot_geotiff_on_axes(self.background_image, lpax, + plot_geotiff_on_axes(self.background_image, self.ax_ptm, epsg_code=4326, band_number=self.bimg_band, cmap=self.bimg_cmap) # --> set tick label format - lpax.xaxis.set_major_formatter(FormatStrFormatter(self.tickstrfmt)) - lpax.yaxis.set_major_formatter(FormatStrFormatter(self.tickstrfmt)) -# lpax.set_xticklabels(np.round(self.plot_xarr, decimals=2), + self.ax_ptm.xaxis.set_major_formatter(FormatStrFormatter(self.tickstrfmt)) + self.ax_ptm.yaxis.set_major_formatter(FormatStrFormatter(self.tickstrfmt)) +# self.ax_ptm.set_xticklabels(np.round(self.plot_xarr, decimals=2), # rotation=45) - plt.setp(lpax.get_xticklabels(), rotation=45) + plt.setp(self.ax_ptm.get_xticklabels(), rotation=45) # --> set title in period or freq if self.tscale == 'period': @@ -979,17 +942,13 @@ def plot(self, fig=None, save_path=None, show=True, else: titlefreq = '{0:.5g} (Hz)'.format(self.plot_freq) - if not self.plot_title: - lpax.set_title('Phase Tensor Map for ' + titlefreq, - fontsize=self.font_size + 2, fontweight='bold') - else: - lpax.set_title(self.plot_title + titlefreq, - fontsize=self.font_size + 2, fontweight='bold') + self.ax_ptm.set_title('{0} {1}'.format(self.plot_title, titlefreq), + fontsize=self.font_size + 2, fontweight='bold') # --> plot induction arrow scale bar ----------------------------------- if self.plot_tipper.find('y') == 0: - parrx = lpax.get_xlim() - parry = lpax.get_ylim() + parrx = self.ax_ptm.get_xlim() + parry = self.ax_ptm.get_ylim() try: axpad = self.arrow_legend_xborderpad except AttributeError: @@ -1003,7 +962,7 @@ def plot(self, fig=None, save_path=None, show=True, try: txtpad = self.arrow_legend_fontpad except AttributeError: - txtpad = .25 * es + txtpad = .25 * self.ellipse_size # make arrow legend postion and arrows coordinates if self.arrow_legend_position == 'lower right': @@ -1045,19 +1004,19 @@ def plot(self, fig=None, save_path=None, show=True, pty = 0 # txy = pay + txtpad - lpax.arrow(pax, - pay, - ptx, - pty, - width=self.arrow_lw, - facecolor=self.arrow_color_real, - edgecolor=self.arrow_color_real, - length_includes_head=False, - head_width=self.arrow_head_width, - head_length=self.arrow_head_length) + self.ax_ptm.arrow(pax, + pay, + ptx, + pty, + width=self.arrow_lw, + facecolor=self.arrow_color_real, + edgecolor=self.arrow_color_real, + length_includes_head=False, + head_width=self.arrow_head_width, + head_length=self.arrow_head_length) # FZ: what is this '|T|=1'? and the horizontal line? - # lpax.text(txa, + # self.ax_ptm.text(txa, # txy, # '|T|=1', # horizontalalignment='center', @@ -1066,13 +1025,14 @@ def plot(self, fig=None, save_path=None, show=True, # END: if self.plot_tipper.find('yes') == 0 --------------------------- # make a grid with color lines - lpax.grid(True, alpha=.3, which='both', color=(0.5, 0.5, 0.5)) + self.ax_ptm.grid(True, alpha=.3, which='both', color=(0.5, 0.5, 0.5)) + self.ax_ptm.set_axisbelow(True) if self.minorticks_on: plt.minorticks_on() # turn on minor ticks automatically # ==> make a colorbar with appropriate colors if self.cb_position is None: - lpax2, kw = mcb.make_axes(lpax, + self.ax_cb, kw = mcb.make_axes(self.ax_ptm, orientation=self.cb_orientation, shrink=.35) # FZ: try to fix colorbar h-position @@ -1084,9 +1044,9 @@ def plot(self, fig=None, save_path=None, show=True, # self.ax2 = divider.append_axes("right", size="5%", pad=0.05) else: - lpax2 = lpfig.add_axes(self.cb_position) + self.ax_cb = self.fig.add_axes(self.cb_position) - if cmap == 'mt_seg_bl2wh2rd': + if self.ellipse_cmap == 'mt_seg_bl2wh2rd': # make a color list self.clist = [(cc, cc, 1) for cc in np.arange(0, 1 + 1. / (nseg), 1. / (nseg))] + \ @@ -1097,26 +1057,26 @@ def plot(self, fig=None, save_path=None, show=True, mt_seg_bl2wh2rd = colors.ListedColormap(self.clist) # make bounds so that the middle is white - bounds = np.arange(ckmin - ckstep, ckmax + 2 * ckstep, ckstep) + bounds = np.arange(self.ellipse_range[0] - self.ellipse_range[2], self.ellipse_range[1] + 2 * self.ellipse_range[2], self.ellipse_range[2]) # normalize the colors norms = colors.BoundaryNorm(bounds, mt_seg_bl2wh2rd.N) # make the colorbar - self.cb = mcb.ColorbarBase(lpax2, + self.cb = mcb.ColorbarBase(self.ax_cb, cmap=mt_seg_bl2wh2rd, norm=norms, orientation=self.cb_orientation, ticks=bounds[1:-1]) else: - if cmap in list(mtcl.cmapdict.keys()): - cmap_input = mtcl.cmapdict[cmap] + if self.ellipse_cmap in list(mtcl.cmapdict.keys()): + cmap_input = mtcl.cmapdict[self.ellipse_cmap] else: - cmap_input = mtcl.cm.get_cmap(cmap) - self.cb = mcb.ColorbarBase(lpax2, + cmap_input = mtcl.cm.get_cmap(self.ellipse_cmap) + self.cb = mcb.ColorbarBase(self.ax_cb, cmap=cmap_input, # mtcl.cmapdict[cmap], - norm=colors.Normalize(vmin=ckmin, - vmax=ckmax), + norm=colors.Normalize(vmin=self.ellipse_range[0], + vmax=self.ellipse_range[1]), orientation=self.cb_orientation) # label the color bar accordingly @@ -1139,19 +1099,19 @@ def plot(self, fig=None, save_path=None, show=True, show_phi = False # JingMingDuan does not want to show the black circle - it's not useful if show_phi is True: ref_ellip = patches.Ellipse((0, .0), - width=es, - height=es, + width=self.ellipse_size, + height=self.ellipse_size, angle=0) ref_ellip.set_facecolor((0, 0, 0)) - ref_ax_loc = list(lpax2.get_position().bounds) + ref_ax_loc = list(self.ax_cb.get_position().bounds) ref_ax_loc[0] *= .95 ref_ax_loc[1] -= .17 ref_ax_loc[2] = .1 ref_ax_loc[3] = .1 - self.ref_ax = lpfig.add_axes(ref_ax_loc, aspect='equal') + self.ref_ax = self.fig.add_axes(ref_ax_loc, aspect='equal') self.ref_ax.add_artist(ref_ellip) - self.ref_ax.set_xlim(-es / 2. * 1.05, es / 2. * 1.05) - self.ref_ax.set_ylim(-es / 2. * 1.05, es / 2. * 1.05) + self.ref_ax.set_xlim(-self.ellipse_size / 2. * 1.05, self.ellipse_size / 2. * 1.05) + self.ref_ax.set_ylim(-self.ellipse_size / 2. * 1.05, self.ellipse_size / 2. * 1.05) plt.setp(self.ref_ax.xaxis.get_ticklabels(), visible=False) plt.setp(self.ref_ax.yaxis.get_ticklabels(), visible=False) self.ref_ax.set_title(r'$\Phi$ = 1') @@ -1334,7 +1294,7 @@ def export_params_to_file(self, save_path=None): tiazmap = np.zeros((xlist.shape[0], ylist.shape[0])) stationmap = np.zeros((xlist.shape[0], ylist.shape[0]), dtype='|S8') - station_location = {} # a dict to store all MT stations (lan lat) + station_location = {} # a dict to store all mt_obj stations (lan lat) # put the information into the zeroed arrays for ii in range(nx): From b7388cb12b5ad8377ad131eb58a1e6feda4d4e2b Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 25 May 2020 23:34:13 -0700 Subject: [PATCH 14/20] added to plot_strike_rose_plot --- examples/scripts/plot_strike_roseplot.py | 68 ++++++++++++++++++------ mtpy/imaging/plotstrike.py | 16 +++--- 2 files changed, 63 insertions(+), 21 deletions(-) diff --git a/examples/scripts/plot_strike_roseplot.py b/examples/scripts/plot_strike_roseplot.py index 4316fd230..b6460d32f 100644 --- a/examples/scripts/plot_strike_roseplot.py +++ b/examples/scripts/plot_strike_roseplot.py @@ -9,31 +9,69 @@ """ import os.path as op import os -os.chdir(r'C:\mtpywin\mtpy') # change to path where mtpy is installed - +# os.chdir(r'C:\mtpywin\mtpy') # change to path where mtpy is installed +os.chdir(r'c:\Users\jpeacock\Documents\GitHub\mtpy\data\edifiles2') # change to path where mtpy is installed from mtpy.imaging.plotstrike import PlotStrike # directory containing edis -edipath = r'C:\mtpywin\mtpy\examples\data\edi_files_2' +edipath = r'c:\Users\jpeacock\Documents\GitHub\mtpy\data\edifiles2' -# full path to file to save to -savepath = r'C:\mtpywin\mtpy\examples\plots\edi_plots' +# # full path to file to save to +# savepath = r'C:\mtpywin\mtpy\examples\plots\edi_plots' # gets edi file names as a list elst = [op.join(edipath,f) for f in os.listdir(edipath) if (f.endswith('.edi'))]# and f.startswith('GL') -strikeplot = PlotStrike(fn_list=elst, - fold=False, - show_ptphimin=False, - plot_type=2 # 1 means divide into separate plots for different decades - # 2 means combine all data into one rose plot - ) -#strikeplot.save_plot(savepath, -# file_format='png', -# fig_dpi=400 -# ) \ No newline at end of file +### this will plot the estimated strike duplicated across quadrants and +# plot the orthogonal component. Plot type 2 will plot all estimates of +# strike into one ploe +strike_plot = PlotStrike(fn_list=elst, + fold=False, + show_ptphimin=False, + plot_type=2) + +# If you want to remove the orthogonal component +strike_plot.plot_orthogonal = False +strike_plot.fig_num = 2 +strike_plot.plot() + +# if you want to plot Tipper strike estimates as well +strike_plot.plot_tipper = True +strike_plot.fig_num = 3 +strike_plot.plot() + +# if you want to rotate the data +# note this will rotate the data N30W +strike_plot.rotation_angle = -30 +strike_plot.fig_num = 4 +strike_plot.plot() + +# to rotate back +# strike_plot.rotation_angle = 30 + +# # if you want to plot per decade plots +# strike_plot.plot_type = 1 +# strike_plot.fig_num = 5 +# strike_plot.text_pad = 1.4 +# strike_plot.plot() + +# # if you want to only look at a few period ranges +# # not the range is given in log10 of the period +# strike_plot.plot_range = [-2, 0] +# strike_plot.fig_num = 6 +# strike_plot.plot() + +# # if you want a vertical orientation instead of horizontal +# strike_plot.plot_orientation = 'v' +# strike_plot.plot_range = 'data' +# strike_plot.fig_num = 7 +# strike_plot.plot() + + + + diff --git a/mtpy/imaging/plotstrike.py b/mtpy/imaging/plotstrike.py index d16c9b62f..1f14ee271 100644 --- a/mtpy/imaging/plotstrike.py +++ b/mtpy/imaging/plotstrike.py @@ -483,9 +483,11 @@ def plot(self, show=True): self.ax_pt = self.fig.add_subplot(n_subplots, nb, jj + nb, polar=True) elif 'v' in self.plot_orientation: - self.ax_inv = self.fig.add_subplot(nb, n_subplots, 2*jj-1, + self.ax_inv = self.fig.add_subplot(nb, n_subplots, + jj * n_subplots - 2, polar=True) - self.ax_pt = self.fig.add_subplot(nb, n_subplots, 2*jj, + self.ax_pt = self.fig.add_subplot(nb, n_subplots, + jj * n_subplots - 1, polar=True) ax_list = [self.ax_inv, self.ax_pt] # vertical orientation @@ -495,8 +497,8 @@ def plot(self, show=True): jj + 2 * nb, polar=True) elif 'v' in self.plot_orientation: - self.ax_tip = self.fig.add_subplot(n_subplots, nb, - 2*jj + 2, + self.ax_tip = self.fig.add_subplot(nb, n_subplots, + jj * n_subplots, polar=True) ax_list.append(self.ax_tip) @@ -623,13 +625,15 @@ def plot(self, show=True): #--> set title of subplot if 'h' in self.plot_orientation: axh.set_title(self.title_dict[bb], fontdict=fd, - bbox={'facecolor': 'white', 'alpha': .25}) + bbox={'facecolor': 'white', + 'alpha': .25}) #--> set the title offset axh.titleOffsetTrans._t = (0, .1) elif 'v' in self.plot_orientation: axh.set_ylabel(self.title_dict[bb], fontdict=fd, - bbox={'facecolor': 'white', 'alpha': .25}, + bbox={'facecolor': 'white', + 'alpha': .25}, rotation=0, labelpad=50) axh.yaxis.set_label_position("right") From ec89752543f569ec05e2e1568fbe00130a81676a Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 25 May 2020 23:36:21 -0700 Subject: [PATCH 15/20] set plot_orthogonal = False as default in plot_strike --- examples/scripts/plot_strike_roseplot.py | 11 +++++------ mtpy/imaging/plotstrike.py | 2 +- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/examples/scripts/plot_strike_roseplot.py b/examples/scripts/plot_strike_roseplot.py index b6460d32f..a8dbc3a62 100644 --- a/examples/scripts/plot_strike_roseplot.py +++ b/examples/scripts/plot_strike_roseplot.py @@ -27,20 +27,19 @@ elst = [op.join(edipath,f) for f in os.listdir(edipath) if (f.endswith('.edi'))]# and f.startswith('GL') -### this will plot the estimated strike duplicated across quadrants and -# plot the orthogonal component. Plot type 2 will plot all estimates of -# strike into one ploe +### this will plot the estimated strike duplicated across quadrants +# Plot type 2 will plot all estimates of strike into one ploe strike_plot = PlotStrike(fn_list=elst, fold=False, - show_ptphimin=False, plot_type=2) -# If you want to remove the orthogonal component -strike_plot.plot_orthogonal = False +# If you want to plot the orthogonal estimation +strike_plot.plot_orthogonal = True strike_plot.fig_num = 2 strike_plot.plot() # if you want to plot Tipper strike estimates as well +strike_plot.plot_orthogonal = False strike_plot.plot_tipper = True strike_plot.fig_num = 3 strike_plot.plot() diff --git a/mtpy/imaging/plotstrike.py b/mtpy/imaging/plotstrike.py index 1f14ee271..28018b968 100644 --- a/mtpy/imaging/plotstrike.py +++ b/mtpy/imaging/plotstrike.py @@ -177,7 +177,7 @@ def __init__(self, **kwargs): self.color_tip = (.2, .65, .2) self.ring_spacing = 10 self.ring_limits = None - self.plot_orthogonal = True + self.plot_orthogonal = False self.font_size = 7 self.text_pad = 0.6 From 46b3adbbb32bc628b05a45d0ec234ed78f3143a8 Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 25 May 2020 23:42:08 -0700 Subject: [PATCH 16/20] added to plot_strike_rose_plot --- examples/scripts/plot_strike_roseplot.py | 47 ++++++++++++------------ 1 file changed, 23 insertions(+), 24 deletions(-) diff --git a/examples/scripts/plot_strike_roseplot.py b/examples/scripts/plot_strike_roseplot.py index a8dbc3a62..574b6a6b5 100644 --- a/examples/scripts/plot_strike_roseplot.py +++ b/examples/scripts/plot_strike_roseplot.py @@ -9,18 +9,17 @@ """ import os.path as op import os -# os.chdir(r'C:\mtpywin\mtpy') # change to path where mtpy is installed -os.chdir(r'c:\Users\jpeacock\Documents\GitHub\mtpy\data\edifiles2') # change to path where mtpy is installed +os.chdir(r'C:\mtpywin\mtpy') # change to path where mtpy is installed from mtpy.imaging.plotstrike import PlotStrike # directory containing edis -edipath = r'c:\Users\jpeacock\Documents\GitHub\mtpy\data\edifiles2' +edipath = r'C:\mtpywin\mtpy\data\edifiles2' -# # full path to file to save to -# savepath = r'C:\mtpywin\mtpy\examples\plots\edi_plots' +# full path to file to save to +savepath = r'C:\mtpywin\mtpy\examples\plots\edi_plots' # gets edi file names as a list @@ -51,25 +50,25 @@ strike_plot.plot() # to rotate back -# strike_plot.rotation_angle = 30 - -# # if you want to plot per decade plots -# strike_plot.plot_type = 1 -# strike_plot.fig_num = 5 -# strike_plot.text_pad = 1.4 -# strike_plot.plot() - -# # if you want to only look at a few period ranges -# # not the range is given in log10 of the period -# strike_plot.plot_range = [-2, 0] -# strike_plot.fig_num = 6 -# strike_plot.plot() - -# # if you want a vertical orientation instead of horizontal -# strike_plot.plot_orientation = 'v' -# strike_plot.plot_range = 'data' -# strike_plot.fig_num = 7 -# strike_plot.plot() +strike_plot.rotation_angle = 30 + +# if you want to plot per decade plots +strike_plot.plot_type = 1 +strike_plot.fig_num = 5 +strike_plot.text_pad = 1.4 +strike_plot.plot() + +# if you want to only look at a few period ranges +# not the range is given in log10 of the period +strike_plot.plot_range = [-2, 0] +strike_plot.fig_num = 6 +strike_plot.plot() + +# if you want a vertical orientation instead of horizontal +strike_plot.plot_orientation = 'v' +strike_plot.plot_range = 'data' +strike_plot.fig_num = 7 +strike_plot.plot() From 1ba3aab9d8e88247be11a7c8d0b156c5065c0fcf Mon Sep 17 00:00:00 2001 From: JP Date: Thu, 28 May 2020 08:09:42 -0700 Subject: [PATCH 17/20] fixed error in plot_tensor_pseudosection line 404, needed to be for key, value in kwargs.items() --- examples/scripts/plot_phase_tensor_section.py | 2 +- mtpy/imaging/phase_tensor_pseudosection.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/scripts/plot_phase_tensor_section.py b/examples/scripts/plot_phase_tensor_section.py index 979632cb7..9b7f1bb44 100644 --- a/examples/scripts/plot_phase_tensor_section.py +++ b/examples/scripts/plot_phase_tensor_section.py @@ -28,6 +28,7 @@ plot_tipper = 'yri', # plot tipper ('y') + 'ri' means real+imag font_size=5, lw=0.5, + rotation_angle=20, ellipse_dict = {'ellipse_colorby':'skew_seg',# option to colour by phimin, phimax, skew, skew_seg 'ellipse_range':[-12, 12, 3]} # set color limits - default 0,90 for phimin or max, # [-12,12] for skew. If plotting skew_seg need to provide @@ -37,7 +38,6 @@ # update ellipse size (tweak for your dataset) plotObj.ellipse_size = 2.5 - plotObj.plot() #plotObj.save_figure(save_fn = op.join(savepath,'PhaseTensorSection.png'), diff --git a/mtpy/imaging/phase_tensor_pseudosection.py b/mtpy/imaging/phase_tensor_pseudosection.py index 9337074cb..cd4729a05 100644 --- a/mtpy/imaging/phase_tensor_pseudosection.py +++ b/mtpy/imaging/phase_tensor_pseudosection.py @@ -400,7 +400,8 @@ def __init__(self, **kwargs): self.subplot_wspace = .05 self.subplot_hspace = .05 - for key, value in kwargs: + print(kwargs) + for key, value in kwargs.items(): setattr(self, key, value) #---need to rotate data on setting rotz From e76b74d0fdd5c55bb36952b1952d999e802158f2 Mon Sep 17 00:00:00 2001 From: JP Date: Thu, 28 May 2020 08:11:49 -0700 Subject: [PATCH 18/20] fixed error in plot_tensor_pseudosection line 404, needed to be for key, value in kwargs.items() --- examples/scripts/plot_phase_tensor_section.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/scripts/plot_phase_tensor_section.py b/examples/scripts/plot_phase_tensor_section.py index 9b7f1bb44..3efc9025a 100644 --- a/examples/scripts/plot_phase_tensor_section.py +++ b/examples/scripts/plot_phase_tensor_section.py @@ -41,4 +41,4 @@ plotObj.plot() #plotObj.save_figure(save_fn = op.join(savepath,'PhaseTensorSection.png'), -# fig_dpi=400) # change to your preferred file resolution \ No newline at end of file +# fig_dpi=400) # change to your preferred file resolution From 22be963998bd180466d671b8e20170145b3f6d2c Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 1 Jun 2020 13:27:12 -0700 Subject: [PATCH 19/20] changed test_z.test_rotate to positive clockwise rotation matrix --- tests/core/test_z.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/core/test_z.py b/tests/core/test_z.py index 1119a9b24..ba640fea8 100644 --- a/tests/core/test_z.py +++ b/tests/core/test_z.py @@ -18,14 +18,15 @@ def setUp(self): def test_rotate(self): alpharad = np.deg2rad(self.rotation_angle) - rotation_matrix = np.matrix([[np.cos(alpharad),np.sin(alpharad)], - [-np.sin(alpharad),np.cos(alpharad)]]) + rotation_matrix = np.matrix([[np.cos(alpharad), -np.sin(alpharad)], + [np.sin(alpharad), np.cos(alpharad)]]) z_array = self.MT.Z.z.copy() ztest = z_array[0] - ztest_rot = np.ma.dot(np.ma.dot(rotation_matrix,ztest),rotation_matrix.T) + ztest_rot = np.ma.dot(np.ma.dot(rotation_matrix,ztest), + rotation_matrix.T) self.MT.Z.rotate(self.rotation_angle) From 0a555358513fa94eb417287d286949a783a42668 Mon Sep 17 00:00:00 2001 From: JP Date: Mon, 1 Jun 2020 13:38:01 -0700 Subject: [PATCH 20/20] fixed test_geometry.test_eccentricity, had the wrong values for control --- tests/analysis/test_geometry.py | 38 +++++++++++++++++---------------- 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/tests/analysis/test_geometry.py b/tests/analysis/test_geometry.py index 62ae3279b..d8dd69600 100644 --- a/tests/analysis/test_geometry.py +++ b/tests/analysis/test_geometry.py @@ -33,24 +33,26 @@ def test_get_dimensionality_from_edi_file(self): def test_get_eccentricity_from_edi_file(self): mt_obj = MT(os.path.normpath(os.path.join(TEST_MTPY_ROOT, "examples/data/edi_files/pb42c.edi"))) - eccentricity_pb42c = (np.array([0.01675639, 0.01038589, 0.00527011, 0.00638819, 0.01483804, - 0.00385233, 0.00513294, 0.00403781, 0.02862114, 0.02689821, - 0.01425044, 0.05686524, 0.05742524, 0.02696736, 0.0275285, - 0.03647819, 0.04721932, 0.06336521, 0.12789841, 0.16409303, - 0.20630821, 0.34261225, 0.3967886, 0.51629705, 0.56645987, - 0.52558696, 0.46954261, 0.48028767, 0.47490701, 0.44927612, - 0.45185046, 0.44143159, 0.43570377, 0.41537978, 0.40546014, - 0.38785478, 0.37174031, 0.34534557, 0.35510941, 0.32282644, - 0.28501461, 0.22463964, 0.20683855]), - np.array([0.01648335, 0.01850986, 0.02356282, 0.02351531, 0.0245181 , - 0.0295284 , 0.03050838, 0.03262934, 0.04868243, 0.04905988, - 0.04742146, 0.04365671, 0.06595206, 0.09502736, 0.11485881, - 0.1330897 , 0.16488362, 0.18365985, 0.2233313 , 0.32718984, - 0.35370171, 0.38087385, 0.55555244, 0.59755086, 0.60360704, - 0.82315127, 0.75538697, 0.73747297, 0.46115785, 0.85300128, - 0.95439195, 0.4485694 , 0.42871072, 0.43533434, 0.43440266, - 0.4950619 , 0.56743593, 0.64693597, 0.74768197, 0.90040195, - 1.01527094, 1.15623132, 1.34441602])) + eccentricity_pb42c = (np.array( + [0.01675639, 0.01038589, 0.00527011, 0.00638819, 0.01483804, + 0.00385233, 0.00513294, 0.00403781, 0.02862114, 0.02689821, + 0.01425044, 0.05686524, 0.05742524, 0.02696736, 0.0275285 , + 0.03647819, 0.04721932, 0.06336521, 0.12789841, 0.16409303, + 0.20630821, 0.34261225, 0.3967886 , 0.51629705, 0.56645987, + 0.52558696, 0.46954261, 0.48028767, 0.47490701, 0.44927612, + 0.45185046, 0.44143159, 0.43570377, 0.41537978, 0.40546014, + 0.38785478, 0.37174031, 0.34534557, 0.35510941, 0.32282644, + 0.28501461, 0.22463964, 0.20683855]), + np.array( + [1.09157801, 1.88513021, 2.52461839, 2.82831998, 2.48686079, + 2.78766939, 2.37250674, 1.16767867, 2.82659526, 1.41909434, + 1.57463881, 2.80007766, 2.73439937, 2.71321983, 2.82682667, + 2.25051091, 1.435991 , 0.43231073, 0.81874286, 1.62593021, + 2.4864899 , 2.91852175, 3.11049976, 3.48724331, 3.6630486 , + 3.54297703, 3.40317309, 3.4175019 , 3.3910377 , 3.28570849, + 3.33766706, 3.31489955, 3.32710981, 3.29740796, 3.28179559, + 3.24496124, 3.2052965 , 3.13248919, 3.11445279, 2.98672687, + 2.85242595, 2.79044042, 2.73406595])) eccentricity = mtg.eccentricity(z_object=mt_obj.Z) for i in range(2):