Package 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)

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)

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)

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')

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)

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)

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)

Estimate strike angle from 2D parts of the phase 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)

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

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.

property alpha

Return the principal axis angle (strike) of PT in degrees

(incl. uncertainties).

Output: - Alpha - Numpy array - Error of Alpha - Numpy array

property azimuth

Returns the azimuth angle related to geoelectric strike in degrees including uncertainties

property beta

Return the 3D-dimensionality angle Beta of PT in degrees (incl. uncertainties).

Output: - Beta - Numpy array - Error of Beta - Numpy array

property det

Return the determinant of PT (incl. uncertainties).

Output: - Det(PT) - Numpy array - Error of Det(PT) - Numpy array

property ellipticity

Returns the ellipticity of the phase tensor, related to dimesionality

property freq

freq array

property invariants

Return a dictionary of PT-invariants.

Contains: trace, skew, det, phimax, phimin, beta

property 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

property 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

property pt

Phase tensor array

property pt_err

Phase tensor error array, must be same shape as pt

rotate(alpha)

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)

Read in Z object and convert information into PhaseTensor object attributes.

property skew

Return the skew of PT (incl. uncertainties).

Output: - Skew(PT) - Numpy array - Error of Skew(PT) - Numpy array

property 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, residualtype='heise')

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)

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)

Read two PT arrays and calculate the ResPT array (incl. uncertainties).

Input: - 2x PT array

Optional: - 2x pt_error array

set_rpt(rpt_array)

Set the attribute ‘rpt’ (ResidualPhaseTensor array).

Input: ResPT array

Test for shape, but no test for consistency!

set_rpt_err(rpt_err_array)

Set the attribute ‘rpt_err’ (ResidualPhaseTensor-error array).

Input: ResPT-error array

Test for shape, but no test for consistency!

mtpy.analysis.pt.edi_file2pt(filename)

Calculate Phase Tensor from Edi-file (incl. uncertainties)

Input: - Edi-file : full path to the Edi-file

Return: - PT object

mtpy.analysis.pt.z2pt(z_array, z_err_array=None)

Calculate Phase Tensor from Z array (incl. uncertainties)

Input: - Z : 2x2 complex valued Numpy array

Optional: - Z-error : 2x2 real valued Numpy array

Return: - PT : 2x2 real valued Numpy array - PT-error : 2x2 real valued Numpy array

mtpy.analysis.pt.z_object2pt(z_object)

Calculate Phase Tensor from Z object (incl. uncertainties)

Input: - Z-object : instance of the MTpy Z class

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)

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)

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)

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()

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)

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.

If the shape changes or the freq are changed need to input those as well.

set_z_err(z_err_array)

set the z_err array.

If the shape changes or the freq are changed need to input those as well.