Analysis¶
Module Distortion¶
mtpy/analysis/distortion.py
Contains functions for the determination of (galvanic) distortion of impedance tensors. The methods used follow Bibby et al 2005. As it has been pointed out in that paper, there are various possibilities for constraining the solution, esp. in the 2D case.
Here we just implement the ‘most basic’ variety for the calculation of the distortion tensor. Other methods can be implemented, but since the optimal assumptions and constraints depend on the application, the actual place for further functions is in an independent, personalised module.
Algorithm Details: Finding the distortion of a Z array. Using the phase tensor so, Z arrays are transformed into PTs first), following Bibby et al. 2005.
First, try to find periods that indicate 1D. From them determine D incl. the g-factor by calculatiing a weighted mean. The g is assumed in order to cater for the missing unknown in the system, it is here set to det(X)^0.5. After that is found, the function no_distortion from the Z module can be called to obtain the unperturbated regional impedance tensor.
Second, if there are no 1D sections: Find the strike angle, then rotate the Z to the principal axis. In order to do that, use the rotate(-strike) method of the Z module. Then take the real part of the rotated Z. As in the 1D case, we need an assumption to get rid of the (2) unknowns: set det(D) = P and det(D) = T, where P,T can be chosen. Common choice is to set one of P,T to an arbitrary value (e.g. 1). Then check, for which values of the other parameter S^2 = T^2+4*P*X_12*X_21/det(X) > 0 holds.
@UofA, 2013 (LK)
Edited by JP, 2016
-
mtpy.analysis.distortion.
find_1d_distortion
(z_object, include_non1d=False)[source]¶ find 1D distortion tensor from z object
ONly use the 1D part of the Z to determine D. Treat all frequencies as 1D, if “include_non1d = True”.
-
mtpy.analysis.distortion.
find_2d_distortion
(z_object, include_non2d=False)[source]¶ find 2D distortion tensor from z object
ONly use the 2D part of the Z to determine D. Treat all frequencies as 2D, if “include_non2d = True”.
-
mtpy.analysis.distortion.
find_distortion
(z_object, g='det', num_freq=None, lo_dims=None)[source]¶ find optimal distortion tensor from z object
automatically determine the dimensionality over all frequencies, then find the appropriate distortion tensor D
Parameters: **z_object** : mtpy.core.z object
- g : [ ‘det’ | ‘01’ | ‘10 ]
type of distortion correction default is ‘det’
- num_freq : int
number of frequencies to look for distortion from the index 0 default is None, meaning all frequencies are used
- lo_dims : list
list of dimensions for each frequency default is None, meaning calculated from data
Returns: **distortion** : np.ndarray(2, 2)
distortion array all real values
- distortion_err : np.ndarray(2, 2)
distortion error array
Examples
Estimate Distortion: >>> import mtpy.analysis.distortion as distortion >>> dis, dis_err = distortion.find_distortion(z_obj, num_freq=12)
-
mtpy.analysis.distortion.
remove_distortion
(z_array=None, z_object=None, num_freq=None, g='det')[source]¶ remove distortion from an impedance tensor using the method outlined by Bibby et al., [2005].
Parameters: **z_array** : np.ndarray((nf, 2, 2))
numpy array of impedance tensor default is None
- z_object : mtpy.core.z object
default is None
- num_freq : int
number of frequecies to look for distortion default is None, meaning look over all frequencies
- g : [ ‘det’ | ‘01’ | ‘10 ]
type of distortion to look for default is ‘det’
Returns: **distortion** : np.ndarray (2, 2)
distortion array
- new_z_obj : mtpy.core.z
z object with distortion removed and error calculated
Examples
Remove Distortion: >>> import mtpy.analysis.distortion as distortion >>> d, new_z = distortion.remove_distortion(z_object=z_obj)
Module Geometry¶
mtpy/mtpy/analysis/geometry.py
Contains classes and functions for handling geometry analysis of impedance tensors:
dimensionality, strike directions, alphas/skews/…
- 1d - 2d : excentricity of ellipses
- 2d - 3d : skew < threshold (to be given as argument)
- strike: frequency - depending angle (incl. 90degree ambiguity)
@UofA, 2013(LK)
Edited by JP, 2016
-
mtpy.analysis.geometry.
dimensionality
(z_array=None, z_object=None, pt_array=None, pt_object=None, skew_threshold=5, eccentricity_threshold=0.1)[source]¶ Esitmate dimensionality of an impedance tensor, frequency by frequency.
Dimensionality is estimated from the phase tensor given the threshold criteria on the skew angle and eccentricity following Bibby et al., 2005 and Booker, 2014.
Returns: **dimensions** : np.ndarray(nf, dtype=int)
an array of dimesions for each frequency the values are [ 1 | 2 | 3 ]
Examples
Estimate Dimesions: >>> import mtpy.analysis.geometry as geometry >>> dim = geometry.dimensionality(z_object=z_obj, >>> skew_threshold=3)
-
mtpy.analysis.geometry.
eccentricity
(z_array=None, z_object=None, pt_array=None, pt_object=None)[source]¶ Estimate eccentricy of a given impedance or phase tensor object
Returns: **eccentricity** : np.ndarray(nf)
eccentricity_err : np.ndarray(nf)
Examples
Estimate Dimesions: >>> import mtpy.analysis.geometry as geometry >>> ec, ec_err= geometry.eccentricity(z_object=z_obj)
-
mtpy.analysis.geometry.
strike_angle
(z_array=None, z_object=None, pt_array=None, pt_object=None, skew_threshold=5, eccentricity_threshold=0.1)[source]¶ Estimate strike angle from 2D parts of the impedance tensor given the skew and eccentricity thresholds
Returns: **strike** : np.ndarray(nf)
an array of strike angles in degrees for each frequency assuming 0 is north, and e is 90. There is a 90 degree ambiguity in the angle.
Examples
Estimate Dimesions: >>> import mtpy.analysis.geometry as geometry >>> strike = geometry.strike_angle(z_object=z_obj, >>> skew_threshold=3)
Module Phase Tensor¶
Following Caldwell et al, 2004
Residual Phase Tensor following Heise et al., [2008]
@UofA, 2013 (LK)
Revised by Peacock, 2016
-
class
mtpy.analysis.pt.
PhaseTensor
(pt_array=None, pt_err_array=None, z_array=None, z_err_array=None, z_object=None, freq=None, pt_rot=0.0)[source]¶ PhaseTensor class - generates a Phase Tensor (PT) object.
Methods include reading and writing from and to edi-objects, rotations combinations of Z instances, as well as calculation of invariants, inverse, amplitude/phase,…
PT is a complex array of the form (n_freq, 2, 2), with indices in the following order:
PTxx: (0,0) - PTxy: (0,1) - PTyx: (1,0) - PTyy: (1,1)- All internal methods are based on (Caldwell et al.,2004) and
- (Bibby et al.,2005), in which they use the canonical cartesian 2D
reference (x1, x2). However, all components, coordinates, and angles for in- and outputs are given in the geographical reference frame:
x-axis = North ; y-axis = East (; z-axis = Down)- Therefore, all results from using those methods are consistent
- (angles are referenced from North rather than x1).
Attributes Description freq array of frequencies associated with elements of impedance tensor. pt phase tensor array pt_err phase tensor error z impedance tensor z_err impedance error rotation_angle rotation angle in degrees Attributes
alpha
Return the principal axis angle (strike) of PT in degrees (incl. alpha_err
azimuth
Returns the azimuth angle related to geoelectric strike in degrees including uncertainties azimuth_err
beta
Return the 3D-dimensionality angle Beta of PT in degrees (incl. beta_err
det
Return the determinant of PT (incl. det_err
ellipticity
Returns the ellipticity of the phase tensor, related to dimesionality ellipticity_err
freq
freq array invariants
Return a dictionary of PT-invariants. only1d
only2d
phimax
Return the angle Phi_max of PT (incl. phimax_err
phimin
Return the angle Phi_min of PT (incl. phimin_err
pt
Phase tensor array pt_err
Phase tensor error array, must be same shape as pt skew
Return the skew of PT (incl. skew_err
trace
Return the trace of PT (incl. trace_err
Methods
rotate
(alpha)Rotate PT array. set_z_object
(z_object)Read in Z object and convert information into PhaseTensor object attributes. -
alpha
¶ - Return the principal axis angle (strike) of PT in degrees
- (incl. uncertainties).
Output: - Alpha - Numpy array - Error of Alpha - Numpy array
-
azimuth
¶ Returns the azimuth angle related to geoelectric strike in degrees including uncertainties
-
beta
¶ Return the 3D-dimensionality angle Beta of PT in degrees (incl. uncertainties).
Output: - Beta - Numpy array - Error of Beta - Numpy array
-
det
¶ Return the determinant of PT (incl. uncertainties).
Output: - Det(PT) - Numpy array - Error of Det(PT) - Numpy array
-
ellipticity
¶ Returns the ellipticity of the phase tensor, related to dimesionality
-
freq
¶ freq array
-
invariants
¶ Return a dictionary of PT-invariants.
Contains: trace, skew, det, phimax, phimin, beta
-
phimax
¶ Return the angle Phi_max of PT (incl. uncertainties).
Phi_max is calculated according to Bibby et al. 2005: Phi_max = Pi2 + Pi1
Output: - Phi_max - Numpy array - Error of Phi_max - Numpy array
-
phimin
¶ Return the angle Phi_min of PT (incl. uncertainties).
- Phi_min is calculated according to Bibby et al. 2005:
- Phi_min = Pi2 - Pi1
Output: - Phi_min - Numpy array - Error of Phi_min - Numpy array
-
pt
¶ Phase tensor array
-
pt_err
¶ Phase tensor error array, must be same shape as pt
-
rotate
(alpha)[source]¶ 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.
-
set_z_object
(z_object)[source]¶ Read in Z object and convert information into PhaseTensor object attributes.
-
skew
¶ Return the skew of PT (incl. uncertainties).
Output: - Skew(PT) - Numpy array - Error of Skew(PT) - Numpy array
-
trace
¶ Return the trace of PT (incl. uncertainties).
Output: - Trace(PT) - Numpy array - Error of Trace(PT) - Numpy array
-
class
mtpy.analysis.pt.
ResidualPhaseTensor
(pt_object1=None, pt_object2=None)[source]¶ PhaseTensor class - generates a Phase Tensor (PT) object DeltaPhi
DeltaPhi = 1 - Phi1^-1*Phi2
Methods
compute_residual_pt
(pt_o1, pt_o2)Read in two instance of the MTpy PhaseTensor class. read_pts
(pt1, pt2[, pt1err, pt2err])Read two PT arrays and calculate the ResPT array (incl. set_rpt
(rpt_array)Set the attribute ‘rpt’ (ResidualPhaseTensor array). set_rpt_err
(rpt_err_array)Set the attribute ‘rpt_err’ (ResidualPhaseTensor-error array). -
compute_residual_pt
(pt_o1, pt_o2)[source]¶ Read in two instance of the MTpy PhaseTensor class.
Update attributes: rpt, rpt_err, _pt1, _pt2, _pt1err, _pt2err
-
read_pts
(pt1, pt2, pt1err=None, pt2err=None)[source]¶ Read two PT arrays and calculate the ResPT array (incl. uncertainties).
Input: - 2x PT array
Optional: - 2x pt_error array
-
-
mtpy.analysis.pt.
edi_file2pt
(filename)[source]¶ Calculate Phase Tensor from Edi-file (incl. uncertainties)
Input: - Edi-file : full path to the Edi-file
Return: - PT object
Module Static Shift¶
module for estimating static shift
Created on Mon Aug 19 10:06:21 2013
@author: jpeacock
-
mtpy.analysis.staticshift.
estimate_static_spatial_median
(edi_fn, radius=1000.0, num_freq=20, freq_skip=4, shift_tol=0.15)[source]¶ Remove static shift from a station using a spatial median filter. This will look at all the edi files in the same directory as edi_fn and find those station within the given radius (meters). Then it will find the medain static shift for the x and y modes and remove it, given that it is larger than the shift tolerance away from 1. A new edi file will be written in a new folder called SS.
Returns: **shift_corrections** : (float, float)
static shift corrections for x and y modes
-
mtpy.analysis.staticshift.
remove_static_shift_spatial_filter
(edi_fn, radius=1000, num_freq=20, freq_skip=4, shift_tol=0.15, plot=False)[source]¶ Remove static shift from a station using a spatial median filter. This will look at all the edi files in the same directory as edi_fn and find those station within the given radius (meters). Then it will find the medain static shift for the x and y modes and remove it, given that it is larger than the shift tolerance away from 1. A new edi file will be written in a new folder called SS.
Returns: **new_edi_fn_ss** : string
new path to the edi file with static shift removed
- shift_corrections : (float, float)
static shift corrections for x and y modes
- plot_obj : mtplot.plot_multiple_mt_responses object
If plot is True a plot_obj is returned If plot is False None is returned
Module Z Invariants¶
Created on Wed May 08 09:40:42 2013
Interpreted from matlab code written by Stephan Thiel 2005
@author: jpeacock
-
class
mtpy.analysis.zinvariants.
Zinvariants
(z_object=None, z_array=None, z_err_array=None, freq=None, rot_z=0)[source]¶ calculates invariants from Weaver et al. [2000, 2003]. At the moment it does not calculate the error for each invariant, only the strike.
Attributes
**inv1** (real off diaganol part normalizing factor) inv2 : imaginary off diaganol normalizing factor inv3 : real anisotropy factor (range from [0,1]) inv4 : imaginary anisotropy factor (range from [0,1]) inv5 : suggests electric field twist inv6 : suggests in phase small scale distortion inv7 : suggests 3D structure strike : strike angle (deg) assuming positive clockwise 0=N strike_err : strike angle error (deg) q : dependent variable suggesting dimensionality Methods
compute_invariants
()Computes the invariants according to Weaver et al., [2000, 2003] rotate
(rot_z)Rotates the impedance tensor by the angle rot_z clockwise positive assuming 0 is North set_freq
(freq)set the freq array, needs to be the same length at z set_z
(z_array)set the z array. set_z_err
(z_err_array)set the z_err array. -
compute_invariants
()[source]¶ Computes the invariants according to Weaver et al., [2000, 2003]
Mostly used to plot Mohr’s circles
In a 1D case: rho = mu (inv1**2+inv2**2)/w & phi = arctan(inv2/inv1)
- Sets the invariants as attributes:
inv1 : real off diaganol part normalizing factor
inv2 : imaginary off diaganol normalizing factor
inv3 : real anisotropy factor (range from [0,1])
inv4 : imaginary anisotropy factor (range from [0,1])
inv5 : suggests electric field twist
inv6 : suggests in phase small scale distortion
inv7 : suggests 3D structure
strike : strike angle (deg) assuming positive clockwise 0=N
strike_err : strike angle error (deg)
q : dependent variable suggesting dimensionality
-
rotate
(rot_z)[source]¶ Rotates the impedance tensor by the angle rot_z clockwise positive assuming 0 is North
-