mtpy.processing package

Subpackages

Submodules

mtpy.processing.base module

Created on Tue Jul 30 17:13:07 2024

@author: jpeacock

class mtpy.processing.base.BaseProcessing(**kwargs)[source]

Bases: KernelDataset

Base processing class containing paths to various files.

config

Processing configuration object.

Type:

object or None

run_summary

Run summary object containing metadata about processing runs.

Type:

RunSummary or None

get_run_summary() → RunSummary[source]

Get the RunSummary object from MTH5 files.

Returns:

Run summary object created from MTH5 file list.

Return type:

RunSummary

has_kernel_dataset() → bool[source]

Check if kernel dataset is set.

Returns:

True if df attribute is not None, False otherwise.

Return type:

bool

has_run_summary() → bool[source]

Check if run summary is set.

Returns:

True if run_summary is not None, False otherwise.

Return type:

bool

property mth5_list: list[str]

Get list of MTH5 file paths.

Returns:

List of MTH5 file paths to get run summary from.

Return type:

list[str]

Raises:

ValueError – If no MTH5 file paths are set.

property run_summary: RunSummary | None

Get the run summary object.

Returns:

Run summary containing metadata about processing runs.

Return type:

RunSummary or None

mtpy.processing.birrp module

BIRRP

  • deals with inputs and outputs from BIRRP

Created on Tue Sep 20 14:33:20 2016

@author: jrpeacock

exception mtpy.processing.birrp.BIRRPParameterError[source]

Bases: Exception

class mtpy.processing.birrp.BIRRPParameters(ilev: int = 0, **kwargs)[source]

Bases: object

Class to hold and produce appropriate BIRRP processing parameters.

Parameters:

  • ilev (int, optional) – Processing level (0 for basic, 1 for advanced), by default 0.

  • **kwargs (dict) – Additional BIRRP parameters to set as attributes.

from_dict(birrp_dict: dict[str, int | float | str | list]) → None[source]

Set BIRRP parameters from dictionary.

Parameters:

birrp_dict (dict) – Dictionary containing BIRRP parameter key-value pairs.

read_config_file(birrp_config_fn: str) → None[source]

Read configuration file and set BIRRP parameters.

Parameters:

birrp_config_fn (str) – Full path to BIRRP configuration file.

to_dict() → dict[str, int | float | str | list[float]][source]

Get appropriate BIRRP parameters as dictionary.

Returns:

Dictionary of BIRRP parameters based on processing level.

Return type:

dict

write_config_file(save_fn: str) → None[source]

Write BIRRP parameters to configuration file.

Parameters:

save_fn (str) – Base filename for configuration file.

class mtpy.processing.birrp.J2Edi(**kwargs)[source]

Bases: object

Convert BIRRP .j output file to EDI format.

Parameters:

  • **birrp_dir** (string) – full path to directory where birrp outputs are

  • **station** (string) – station name

  • **survey_config_fn** (string) – full path to survey configuration file with information on location and site setup must have a key that is the same as station.

  • **birrp_config_fn** (string) – full path to configuration file that was used to process with (all the birrp parameters used). If None is input, the file is searched for, if it is not found, the processing parameters are used from the .j file.

  • **j_fn** – full path to j file. If none is input the .j file is searched for in birrp_dir.

Example

>>> import mtpy.proceessing.birrp as birrp
>>> j2edi_obj = birrp.J_To_Edi()
>>> j2edi_obj.birrp_dir = r"/home/data/mt01/BF/256"
>>> j2edi_obj.station = 'mt01'
>>> j2edi_obj.survey_config_fn = r"/home/data/2016_survey.cfg"
>>> j2edi_obj.write_edi_file()
get_birrp_config_fn() → None[source]

Locate BIRRP configuration file in processing directory.

Notes

Sets self.birrp_config_fn to the found configuration file path.

get_j_file(birrp_dir: str | None = None) → None[source]

Locate BIRRP .j output file.

Parameters:

birrp_dir (str or None, optional) – Directory containing BIRRP outputs, by default None.

Raises:

MTpyError_inputarguments – If birrp_dir does not exist.

read_birrp_config_fn(birrp_config_fn: str | None = None) → None[source]

Read BIRRP configuration file.

Parameters:

birrp_config_fn (str or None, optional) – Full path to BIRRP configuration file, by default None.

read_survey_config_fn(survey_config_fn: str | None = None) → None[source]

Read survey configuration file.

Parameters:

survey_config_fn (str or None, optional) – Full path to survey configuration file, by default None.

Raises:

MTpyError_inputarguments – If configuration file does not exist.

write_edi_file(station: str | None = None, birrp_dir: str | None = None, survey_config_fn: str | None = None, birrp_config_fn: str | None = None, copy_path: str | None = None) → str[source]

Convert BIRRP .j file to EDI format.

Parameters:

  • station (str or None, optional) – Station name, by default None.

  • birrp_dir (str or None, optional) – Full path to BIRRP output directory, by default None.

  • survey_config_fn (str or None, optional) – Full path to survey configuration file, by default None.

  • birrp_config_fn (str or None, optional) – Full path to BIRRP configuration file, by default None.

  • copy_path (str or None, optional) – Directory to copy EDI file to, by default None.

Returns:

Full path to created EDI file.

Return type:

str

Raises:

MTpyError_inputarguments – If station name or birrp_dir is not provided or invalid. .. note:: The survey_config_fn is a file that has the structure: [station] b_instrument_amplification = 1 b_instrument_type = coil b_logger_gain = 1 b_logger_type = zen b_xaxis_azimuth = 0 b_yaxis_azimuth = 90 box = 26 date = 2015/06/09 e_instrument_amplification = 1 e_instrument_type = Ag-Agcl electrodes e_logger_gain = 1 e_logger_type = zen e_xaxis_azimuth = 0 e_xaxis_length = 100 e_yaxis_azimuth = 90 e_yaxis_length = 100 elevation = 2113.2 hx = 2274 hy = 2284 hz = 2254 lat = 37.7074236995 location = Earth lon = -118.999542099 network = USGS notes = Generic config file rr_box = 25 rr_date = 2015/06/09 rr_hx = 2334 rr_hy = 2324 rr_lat = 37.6909139779 rr_lon = -119.028707542 rr_station = 302 sampling_interval = all save_path = homemtdatasurvey_01mt_01 station = 300 station_type = mt This file can be written using mtpy.utils.configfile:: >>> import mtpy.utils.configfile as mtcfg >>> station_dict = {} >>> station_dict[‘lat’] = 21.346 >>> station_dict[‘lon’] = 122.45654 >>> station_dict[‘elev’] = 123.43 >>> cfg_fn = r”homemtdatasurvey_01” >>> mtcfg.write_dict_to_configfile({station: station_dict}, cfg_fn)

class mtpy.processing.birrp.ScriptFile(script_fn: str | None = None, fn_arr: ndarray | None = None, **kwargs)[source]

Bases: BIRRPParameters

Class to read and write BIRRP script files.

Parameters:

  • script_fn (str or None, optional) – Full path to script file, by default None.

  • fn_arr (np.ndarray or None, optional) – Numpy structured array containing file information, by default None.

  • **kwargs (dict) – Additional BIRRP parameters.

\*\*fn_arr**

numpy.ndarray([[block 1], [block 2]])

Type:

numpy.ndarray

Note

[block n] is a numpy structured array with data type

Name

Description

Type

fn

file path/name

string

nread

number of points to read

int

nskip

number of points to skip

int

comp

component

[ ex | ey | hx hy | hz ]

calibration_fn

calibration file path/name

string

rr

a remote reference channel

[ True | False ]

rr_num

remote reference pair number

int

start

start time iso format

Timestamp

stop

stop time iso format

Timestamp

station

station name

string

sampling_rate

sampling rate

int

BIRRP Parameters:

parameter

description

ilev

processing mode 0 for basic and 1 for advanced RR-2 stage

nout

Number of Output time series (2 or 3-> for BZ)

ninp

Number of input time series for E-field (1,2,3)

nref

Number of reference channels (2 for MT)

nrr

bounded remote reference (0) or 2 stage bounded influence (1)

tbw

Time bandwidth for Sepian sequence

deltat

Sampling rate (+) for (s), (-) for (Hz)

nfft

Length of FFT (should be even)

nsctinc

section increment divisor (2 to divide by half)

nsctmax

Number of windows used in FFT

nf1

1st frequency to extract from FFT window (>=3)

nfinc

frequency extraction increment

nfsect

number of frequencies to extract

mfft

AR filter factor, window divisor (2 for half)

uin

Quantile factor determination

ainlin

Residual rejection factor low end (usually 0)

ainuin

Residual rejection factor high end (.95-.99)

c2threshb

Coherence threshold for magnetics (0 if undesired)

c2threshe

Coherence threshold for electrics (0 if undesired)

nz
Threshold for Bz (0=separate from E, 1=E threshold,

2=E and B)

Input if 3 B components else None

c2thresh1

Squared coherence for Bz, input if NZ=0, Nout=3

perlo

longest period to apply coherence threshold over

perhi

shortes period to apply coherence threshold over

ofil
Output file root(usually three letters, can add full

path)

nlev

Output files (0=Z; 1=Z,qq; 2=Z,qq,w; 3=Z,qq,w,d)

nprej

number of frequencies to reject

prej

frequencies to reject (+) for period, (-) for frequency

npcs

Number of independent data to be processed (1 for one segement)

nar

Prewhitening Filter (3< >15) or 0 if not desired’,

imode

Output file mode (0=ascii; 1=binary; 2=headerless ascii; 3=ascii in TS mode’,

jmode
input file mode (0=user defined; 1=sconvert2tart time

YYYY-MM-DD HH:MM:SS)’,

nread

Number of points to be read for each data set (if segments>1 -> npts1,npts2…)’,

nfil
Filter parameters (0=none; >0=input parameters;

<0=filename)

nskip
Skip number of points in time series (0) if no skip,

(if segements >1 -> input1,input2…)’,

nskipr

Number of points to skip over (0) if none, (if segements >1 -> input1,input2…)’,

thetae

Rotation angles for electrics (relative to geomagnetic North)(N,E,rot)’,

thetab

Rotation angles for magnetics (relative to geomagnetic North)(N,E,rot)’,

thetar

Rotation angles for calculation (relative to geomagnetic North)(N,E,rot)’

Note

Currently only supports jmode = 0 and imode = 0

See also

BIRRP Manual and publications by Chave and Thomson for more details on the parameters found at:

http://www.whoi.edu/science/AOPE/people/achave/Site/Next1.html

property comp_list: list[str]

Get list of component names.

Returns:

List of component names (ex, ey, hx, hy, hz, rrhx_*, rrhy_*).

Return type:

list[str]

Raises:

ValueError – If number of components is invalid.

property deltat: float

Get sampling interval (negative for Hz).

Returns:

Sampling rate (negative value for Hz).

Return type:

float

make_fn_lines_block_00(fn_arr: ndarray) → list[str][source]

Make file lines for first processing block in script file.

Parameters:

fn_arr (np.ndarray) – Structured array with file information.

Returns:

Lines containing nread, filter_fn, filename, and nskip.

Return type:

list[str]

make_fn_lines_block_n(fn_arr: ndarray) → list[str][source]

Make file lines for subsequent processing blocks in script file.

Parameters:

fn_arr (np.ndarray) – Structured array with file information.

Returns:

Lines containing filename and nskip.

Return type:

list[str]

property nout: int

Get number of output time series.

Returns:

Number of output channels (excluding remote reference).

Return type:

int

property npcs: int

Get number of processing blocks.

Returns:

Number of independent data blocks to process.

Return type:

int

property nref: int

Get number of reference channels.

Returns:

Number of remote reference channels.

Return type:

int

write_script_file(script_fn: str | None = None, ofil: str | None = None) → None[source]

Write BIRRP script file.

Parameters:

  • script_fn (str or None, optional) – Full path to script file, by default None.

  • ofil (str or None, optional) – Output file root name, by default None.

exception mtpy.processing.birrp.ScriptFileError[source]

Bases: Exception

mtpy.processing.birrp.run(birrp_exe: str, script_file: str) → bytes[source]

Run BIRRP script file from command line via subprocess.

Parameters:

  • birrp_exe (str) – Full path to compiled BIRRP executable.

  • script_file (str) – Full path to BIRRP script file.

Returns:

Output from BIRRP process.

Return type:

bytes

Raises:

MTpyError_inputarguments – If birrp_exe does not exist.

See also

BIRRP

mtpy.processing.filter module

mtpy/processing/filter.py

Functions for the frequency filtering of raw time series.

@UofA, 2013 (LK)

Revised 2017 JP

mtpy.processing.filter.adaptive_notch_filter(bx, df=100, notches=[50, 100], notchradius=0.5, freqrad=0.9, rp=0.1, dbstop_limit=5.0)[source]

Adaptive_notch_filter(bx, df, notches=[50,100], notchradius=.3, freqrad=.9) will apply a notch filter to the array bx by finding the nearest peak around the supplied notch locations. The filter is a zero-phase Chebyshev type 1 bandstop filter with minimal ripples.

Arguments::
bxnp.ndarray(len_time_series)

time series to filter

dffloat

sampling frequency in Hz

notches : list of frequencies (Hz) to filter

notchradiusfloat

radius of the notch in frequency domain (Hz)

freqradfloat

radius to searching for peak about notch from notches

rpfloat

ripple of Chebyshev type 1 filter, lower numbers means less ripples

dbstop_limitfloat (in decibels)

limits the difference between the peak at the notch and surrounding spectra. Any difference above dbstop_limit will be filtered, anything less will not

Outputs:

    **bx** : np.ndarray(len_time_series)
             filtered array

    **filtlst** : list
                  location of notches and power difference between peak of
                  notch and average power.

..Example: ::

    >>> import RemovePeriodicNoise_Kate as rmp
    >>> # make a variable for the file to load in
    >>> fn = r"/home/MT/mt01_20130101_000000.BX"
    >>> # load in file, if the time series is not an ascii file
    >>> # might need to add keywords to np.loadtxt or use another
    >>> # method to read in the file
    >>> bx = np.loadtxt(fn)
    >>> # create a list of frequencies to filter out
    >>> freq_notches = [50, 150, 200]
    >>> # filter data
    >>> bx_filt, filt_lst = rmp.adaptiveNotchFilter(bx, df=100.
    >>> ...                                         notches=freq_notches)
    >>> #save the filtered data into a file
    >>> np.savetxt(r"/home/MT/Filtered/mt01_20130101_000000.BX", bx_filt)
Notes::

Most of the time the default parameters work well, the only thing you need to change is the notches and perhaps the radius. I would test it out with a few time series to find the optimum parameters. Then make a loop over all you time series data. Something like

>>> import os
>>> dirpath = r"/home/MT"
>>> #make a director to save filtered time series
>>> save_path = r"/home/MT/Filtered"
>>> if not os.path.exists(save_path):
>>>     os.mkdir(save_path)
>>> for fn in os.listdir(dirpath):
>>>     bx = np.loadtxt(os.path.join(dirpath, fn)
>>>     bx_filt, filt_lst = rmp.adaptiveNotchFilter(bx, df=100.
>>>     ...                                         notches=freq_notches)
>>>     np.savetxt(os.path.join(save_path, fn), bx_filt)
mtpy.processing.filter.butter_bandpass(lowcut, highcut, samplingrate, order=4)[source]

Butter bandpass.

mtpy.processing.filter.butter_bandpass_filter(data, lowcut, highcut, samplingrate, order=4)[source]

Butter bandpass filter.

mtpy.processing.filter.low_pass(f, low_pass_freq, cutoff_freq, sampling_rate)[source]

Low pass.

mtpy.processing.filter.remove_periodic_noise(filename, dt, noiseperiods, save='n')[source]

RemovePeriodicNoise will take a window of length noise period and compute the median of signal for as many windows that can fit within the data. This median window is convolved with a series of delta functions at each window location to create a noise time series. This is then subtracted from the data to get a ‘noise free’ time series.

Arguments::
filenamestring (full path to file) or array

name of file to have periodic noise removed from can be an array

dtfloat

time sample rate (s)

noiseperiodslist

a list of estimated periods with a range of values to look around [[noiseperiod1,df1]…] where df1 is a fraction value find the peak about noiseperiod1 must be less than 1. (0 is a good start, but if you’re periodic noise drifts, might need to adjust df1 to .2 or something)

save[ ‘y’ | ‘n’ ]
  • ‘y’ to save file to:

    os.path.join(os.path.dirname(filename), ‘Filtered’, fn)

  • ‘n’ to return the filtered time series

Outputs:

    **bxnf** : np.ndarray
               filtered time series

    **pn** : np.ndarray
            periodic noise time series

    **fitlst** : list
                 list of peaks found in time series

..Example: ::

    >>> import RemovePeriodicNoise_Kate as rmp
    >>> # make a variable for the file to load in
    >>> fn = r"/home/MT/mt01_20130101_000000.BX"
    >>> # filter data assuming a 12 second period in noise and save data
    >>> rmp.remove_periodic_noise(fn, 100., [[12,0]], save='y')
Notes::

Test out the periodic noise period at first to see if it drifts. Then loop over files

>>> import os
>>> dirpath = r"/home/MT"
>>> for fn in os.listdir(dirpath):
>>>     rmp.remove_periodic_noise(fn, 100., [[12,0]], save='y')
mtpy.processing.filter.tukey(window_length, alpha=0.2)[source]

The Tukey window, also known as the tapered cosine window, can be regarded as a cosine lobe of width alpha * N / 2 that is convolved with a rectangle window of width (1 - alpha / 2). At alpha = 0 it becomes rectangular, and at alpha = 1 it becomes a Hann window.

output

Reference:

mtpy.processing.filter.zero_pad(input_array, power=2, pad_fill=0)[source]

Pad the input array with pad_fill to the next power of power.

For faster fft computation pad the array to the next power of 2 with zeros

Arguments::
input_arraynp.ndarray (only 1-d arrays are supported at the

moment)

power[ 2 | 10 ]

power look for

pad_fillfloat or int

pad the array with this

Output::

pad_array : np.ndarray padded with pad_fill

mtpy.processing.tf module

mtpy/processing/tf.py

Functions for the time-frequency analysis of time series data.

Output can be visualised with the help of mtpy/imaging/spectrogram.py

JP, 2013

mtpy.processing.tf.dctrend(f)[source]

dctrend(f) will remove a dc trend from the function f.

Arguments:

fnp.ndarray()

array to remove dc trend from

Returns:

fdcnp.ndarray()

array f with dc component removed

mtpy.processing.tf.decimate(f, m, window_function='hanning')[source]

resamples the data at the interval m so that the returned array is len(f)/m samples long

Arguments:

fnp.ndarray

array to be decimated

mint

decimation factor

window_functionwindowing function to apply to the data

to make sure the is no Gibbs ringing or aliasing see scipy.signal.window for all the options

Returns:

fdecnp.ndarray()

array f decimated by factor m

mtpy.processing.tf.dwindow(window)[source]

Calculates the derivative of the given window. Used for reassignment methods

Arguments:

windownp.ndarray

some sort of windowed array

Returns:

dwinnp.ndarray

derivative of window

mtpy.processing.tf.gausswin(window_len, alpha=2.5)[source]

gausswin will compute a gaussian window of length winlen with a variance of alpha

Arguments:

window_len: int

length of desired window

alphafloat

1/standard deviation of window, ie full width half max of window

Returns:

gauss_windownp.array

gaussian window

mtpy.processing.tf.modifiedb(fx, tstep=32, nfbins=1024, df=1.0, nh=255, beta=0.2)[source]

Calculates the modified b distribution as defined by cosh(n)^-2 beta for an array fx. Supposed to remove cross terms in the WVD.

Arguments:

fxlist or np.ndarray

the function to have a spectrogram computed for for cross-correlation input as [fx1, fx2]

nhint (should be odd)

window length for each time step default is None and window is calculated automatically

tstepint

number of sample between short windows default is 2**5 = 32

dffloat

sampling frequency

nfbinsint (should be power of 2 and equal or larger than nh)

number of frequency bins

betafloat

smoothing coefficient ussully between [0, 1]

Returns:

tfarraynp.ndarray(nfbins/2, len(fx)/tstep)

SPWVD spectrogram in units of amplitude

tlstnp.array()

array of time instances for each window calculated

flstnp.ndarray(nfbins/2)

frequency array containing only positive frequencies where the Fourier coeffients were calculated

mtpy.processing.tf.normalize_L2(f)[source]

normalize_L2(f) returns the array f normalized by the L2 norm -> f/(sqrt(sum(abs(x_i)^2))).

Arguments
fnp.ndarray()

array to be normalized

Returns:

fnormnp.ndarray()

array f normalized in L2 sense

mtpy.processing.tf.padzeros(f, npad=None, pad_pattern=None)[source]

padzeros(f) returns a function that is padded with zeros to the next power of 2 for faster processing for fft or to length npad if given.

Arguments:

fnp.ndarray(m, n)

array to pad

npadint

length to pad to if None finds next power of 2

pad_patternint or float

pattern to pad with if None set zeros

Returns:

fpadnp.ndarray(m, npad)

array f padded to length npad with pad_pattern

Example:

:: >>> x_array = np.sin(np.arange(0, 2, .01)*np.pi/3) >>> print len(x_array) >>> x_array_pad = padzeros(x_array) >>> print len(x_array_pad)

mtpy.processing.tf.reassigned_smethod(fx, nh=127, tstep=16, nfbins=512, df=1.0, alpha=4, thresh=0.01, L=5)[source]

Calulates the reassigned S-method as described by Djurovic[1999] by using the spectrogram to estimate the reassignment.

Arguments:

for cross-correlation input as [fx1, fx2]

Lint (should be odd)

length of window for S-method calculation, higher numbers tend toward WVD

nhint (should be power of 2)

window length for each time step default is 2**8 = 256

alphafloat

inverse of full-width half max of gaussian window, smaller numbers mean broader windows.

threshfloat

threshold for reassignment, lower numbers more points reassigned, higer numbers less points reassigned

tstepint

number of sample between short windows default is 2**7 = 128

dffloat

sampling frequency

nfbinsint (should be power of 2 and equal or larger than nh)

number of frequency bins

Returns:

tfarraynp.ndarray(nfbins/2, len(fx)/tstep)

S-method spectrogram in units of amplitude

tlstnp.array()

array of time instances for each window calculated

flstnp.ndarray(nfbins/2)

frequency array containing only positive frequencies where the Fourier coeffients were calculated

smnp.ndarray(nfbins/2, len(fx)/tstep)

S-method spectrogram in units of amplitude

mtpy.processing.tf.reassigned_stft(fx, nh=63, tstep=32, nfbins=1024, df=1.0, alpha=4, threshold=None)[source]

Computes the reassigned spectrogram by estimating the center of gravity of the signal and condensing dispersed energy back to that location. Works well for data with minimal noise and strong spectral structure.

Arguments:

fxnp.ndarray

time series to be analyzed

nhint(should be odd)

length of gaussian window that is applied to the short time intervals default is 127

tstepint

time step for each window calculation default is 64

nfbinsint (should be a power of 2 and larger or equal to nh

number of frequency bins to calculate, note result will be length nfbins/2 default is 1024

dffloat or int

sampling frequency (Hz)

alphafloat

reciprocal of full width half max of gaussian window default is 4

thresholdfloat

threshold value for reassignment If None the threshold is automatically calculated default is None

returns:
np.ndarray(nfbins/2, len(fx)/tstep)

reassigned spectrogram in units of amplitude

tlstnp.array()

array of time instances for each window calculated

flstnp.ndarray(nfbins/2)

frequency array containing only positive frequencies where the Fourier coeffients were calculated

stftnp.ndarray(nfbins/2, len(fx)/tstep)

standard spectrogram calculated from stft in units of amplitude

rtype:

rtfarray

mtpy.processing.tf.robust_smethod(fx, L=5, nh=128, tstep=32, nfbins=1024, df=1.0, robusttype='median', sigmaL=None, alpha=0.325)[source]

Computes the robust Smethod via the robust spectrogram.

Arguments:

fxlist or np.ndarray

the function to have a spectrogram computed for for cross-correlation input as [fx1, fx2]

Lint (should be odd)

length of window for S-method calculation, higher numbers tend toward WVD

nhint (should be power of 2)

window length for each time step default is 2**8 = 256

ngint (should be odd)

length of smoothing window along frequency plane

tstepint

number of sample between short windows default is 2**7 = 128

dffloat

sampling frequency

nfbinsint (should be power of 2 and equal or larger than nh)

number of frequency bins

robusttype[ ‘median’ | ‘L’ ]

type of robust STFT to compute. default is ‘median’

simgaLfloat

full-width half max of gaussian window applied in frequency

Returns:

tfarraynp.ndarray(nfbins/2, len(fx)/tstep)

S-method spectrogram in units of amplitude

tlstnp.array()

array of time instances for each window calculated

flstnp.ndarray(nfbins/2)

frequency array containing only positive frequencies where the Fourier coeffients were calculated

pxxnp.ndarray(nfbins/2, len(fx)/tstep)

STFT spectrogram in units of amplitude

mtpy.processing.tf.robust_stft_L(fx, alpha=0.325, nh=256, tstep=32, df=1.0, nfbins=1024)[source]

Calculates the robust spectrogram by estimating the vector median and summing terms estimated by alpha coefficients.

Arguments:

fxlist or np.ndarray

the function to have a spectrogram computed for for cross-correlation input as [fx1, fx2]

alphafloat

robust parameter [0,.5] -> 0 gives spectrogram, 0.5 gives median stft default is 0.325

nhint (should be power of 2)

window length for each time step default is 2**8 = 256

tstepint

number of sample between short windows default is 2**7 = 128

dffloat

sampling frequency

nfbinsint (should be power of 2 and equal or larger than nh)

number of frequency bins

Returns:

tfarraynp.ndarray(nfbins/2, len(fx)/tstep)

spectrogram in units of amplitude

tlstnp.array()

array of time instances for each window calculated

flstnp.ndarray(nfbins/2)

frequency array containing only positive frequencies where the Fourier coeffients were calculated

mtpy.processing.tf.robust_stft_median(fx, nh=256, tstep=32, df=1.0, nfbins=1024)[source]

Calculates the robust spectrogram using the vector median simplification.

Arguments:

fxlist or np.ndarray

the function to have a spectrogram computed for for cross-correlation input as [fx1, fx2]

nhint (should be power of 2)

window length for each time step default is 2**8 = 256

tstepint

number of sample between short windows default is 2**7 = 128

dffloat

sampling frequency (Hz)

nfbinsint (should be power of 2 and equal or larger than nh)

number of frequency bins

Returns:

tfarraynp.ndarray(nfbins/2, len(fx)/tstep)

spectrogram in units of amplitude

tlstnp.array()

array of time instances for each window calculated

flstnp.ndarray(nfbins/2)

frequency array containing only positive frequencies where the Fourier coeffients were calculated

mtpy.processing.tf.robust_wvd(fx, nh=127, ng=15, tstep=16, nfbins=256, df=1.0, sigmat=None, sigmaf=None)[source]

Calculate the robust Wigner-Ville distribution for an array fx. Smoothed with Gaussians windows to get best localization.

Arguments:

fxlist or np.ndarray

the function to have a spectrogram computed for for cross-correlation input as [fx1, fx2]

nhint (should be power of 2)

window length for each time step default is 2**8 = 256

tstepint

number of sample between short windows default is 2**7 = 128

ngint (should be odd)

length of smoothing window along frequency plane default is 2**4-1 = 15

dffloat

sampling frequency

nfbinsint (should be power of 2 and equal or larger than nh)

number of frequency bins

sigmatfloat

std of window h, ie full width half max of gaussian default is None and sigmat is calculate automatically

sigmaffloat

std of window g, ie full width half max of gaussian default is None and sigmaf is calculate automatically

Returns:

tfarraynp.ndarray(nfbins/2, len(fx)/tstep)

spectrogram in units of amplitude

tlstnp.array()

array of time instances for each window calculated

flstnp.ndarray(nfbins/2)

frequency array containing only positive frequencies where the Fourier coeffients were calculated

mtpy.processing.tf.sinc_filter(f, fcutoff=10.0, w=10.0, dt=0.001)[source]

Applies a sinc filter of width w to the function f by multipling in the frequency domain.

Arguments:

fnp.ndarray()

array to filter

fcuttofffloat

cutoff frequency of sinc filter

wfloat

length of filter

dtfloat

sampling rate in time (s)

Returns:

f_filtnp.ndarray()

f with sinc filter applied

mtpy.processing.tf.smethod(fx, L=11, nh=256, tstep=128, ng=1, df=1.0, nfbins=1024, sigmaL=None)[source]

Calculates the smethod by estimating the STFT first and computing the WV of window length L in the frequency domain.

For larger L more of WV estimation, if L=0 get back STFT

Arguments:

fxlist or np.ndarray

the function to have a spectrogram computed for for cross-correlation input as [fx1, fx2]

Lint (should be odd)

length of window for S-method calculation, higher numbers tend toward WVD

nhint (should be power of 2)

window length for each time step default is 2**8 = 256

ngint (should be odd)

length of smoothing window along frequency plane

tstepint

number of sample between short windows default is 2**7 = 128

dffloat

sampling frequency

nfbinsint (should be power of 2 and equal or larger than nh)

number of frequency bins

Returns:

tfarraynp.ndarray(nfbins/2, len(fx)/tstep)

S-method spectrogram in units of amplitude

tlstnp.array()

array of time instances for each window calculated

flstnp.ndarray(nfbins/2)

frequency array containing only positive frequencies where the Fourier coeffients were calculated

pxxnp.ndarray(nfbins/2, len(fx)/tstep)

STFT spectrogram in units of amplitude

mtpy.processing.tf.specwv(fx, tstep=32, nfbins=1024, nhs=256, nhwv=511, ngwv=7, df=1.0)[source]

Calculates the Wigner-Ville distribution mulitplied by the STFT windowed by the common gaussian window h for an array f. Handy for removing cross terms in the wvd.

Arguments:

fxlist or np.ndarray

the function to have a spectrogram computed for

tstepint

number of sample between short windows default is 2**7 = 128

nhsint (should be power of 2)

window length for each time step to caclulate STFT default is 2**8 = 256 and window is calculated automatically

nhwvint (should be odd)

length of smoothing window for each time step to calculate WVD. default is 2**9-1 = 511

ngwvint (should be odd)

length of frequency smoothing window for each time step to calculate WVD. default is 2**3-1 = 7

dffloat

sampling frequency

nfbinsint (should be power of 2 and equal or larger than nh)

number of frequency bins

Returns:

tfarraynp.ndarray(nfbins/2, len(fx)/tstep)

SPWVD spectrogram in units of amplitude

tlstnp.array()

array of time instances for each window calculated

flstnp.ndarray(nfbins/2)

frequency array containing only positive frequencies where the Fourier coeffients were calculated

mtpy.processing.tf.spwvd(fx, tstep=32, nfbins=1024, df=1.0, nh=None, ng=None, sigmat=None, sigmaf=None)[source]

Calculates the smoothed pseudo Wigner-Ville distribution for an array fx. Smoothed with Gaussians windows to get best localization.

Can be input as [fx1, fx2] to compute cross spectra.

Arguments:

fxlist or np.ndarray

the function to have a spectrogram computed for for cross-correlation input as [fx1, fx2]

nhint (should be odd)

window length for each time step default is None and window is calculated automatically

tstepint

number of sample between short windows default is 2**7 = 128

ngint (should be odd)

length of smoothing window along frequency plane

dffloat

sampling frequency

nfbinsint (should be power of 2 and equal or larger than nh)

number of frequency bins

sigmatfloat

std of window h, ie full width half max of gaussian default is None and sigmat is calculate automatically

sigmaffloat

std of window g, ie full width half max of gaussian default is None and sigmaf is calculate automatically

Returns:

tfarraynp.ndarray(nfbins/2, len(fx)/tstep)

SPWVD spectrogram in units of amplitude

tlstnp.array()

array of time instances for each window calculated

flstnp.ndarray(nfbins/2)

frequency array containing only positive frequencies where the Fourier coeffients were calculated

mtpy.processing.tf.stfbss(X, nsources=5, ng=31, nh=511, tstep=63, df=1.0, nfbins=1024, tftol=1e-08, L=7, normalize=True, tftype='spwvd', alpha=0.38)[source]
btfssX,nsources=5,ng=2**5-1,nh=2**9-1,tstep=2**6-1,df=1.0,nfbins=2**10,

tftol=1.E-8,normalize=True)

estimates sources using a blind source algorithm based on spatial time-frequency distributions. At the moment this algorithm uses the SPWVD to estimate TF distributions.

Arguments
X = m x n array of time series, where m is number of time series and n

is length of each time series

nsources = number of estimated sources ng = frequency window length nh = time window length tstep = time step increment df = sampling frequency (Hz) nfbins = number of frequencies tftol = tolerance for a time-frequency point to be estimated as a cross

term or as an auto term, the higher the number the more auto terms.

normalization = True or False, True to normalize, False if already

normalized

Returns:

Se = estimated individual signals up to a permutation and scale Ae = estimated mixing matrix as X=A*S

mtpy.processing.tf.stft(fx, nh=256, tstep=128, ng=1, df=1.0, nfbins=1024)[source]

calculate the spectrogam of the given function by calculating the fft of a window of length nh at each time instance with an interval of tstep. The frequency resolution is nfbins.

Can compute the cross STFT by inputting fx as [fx1, fx2]

Arguments:

fxlist or np.ndarray

the function to have a spectrogram computed for for cross-correlation input as [fx1, fx2]

nhint (should be power of 2)

window length for each time step default is 2**8 = 256

tstepint

number of sample between short windows default is 2**7 = 128

ngint (should be odd)

length of smoothing window along frequency plane

dffloat

sampling frequency

nfbinsint (should be power of 2 and equal or larger than nh)

number of frequency bins

Returns:

tfarraynp.ndarray(nfbins/2, len(fx)/tstep)

spectrogram in units of amplitude

tlstnp.array()

array of time instances for each window calculated

flstnp.ndarray(nfbins/2)

frequency array containing only positive frequencies where the Fourier coeffients were calculated

mtpy.processing.tf.wvd(fx, nh=255, tstep=32, nfbins=1024, df=1.0)[source]

calculates the Wigner-Ville distribution of f.

Can compute the cross spectra by inputting fx as [fx1,fx2]

Arguments:

fxlist or np.ndarray

the function to have a spectrogram computed for for cross-correlation input as [fx1, fx2]

nhint (should be odd)

window length for each time step default is 2**8-1 = 255

tstepint

number of sample between short windows default is 2**7 = 128

dffloat

sampling frequency

nfbinsint (should be power of 2 and equal or larger than nh)

number of frequency bins

Returns:

tfarraynp.ndarray(nfbins/2, len(fx)/tstep)

spectrogram in units of amplitude

tlstnp.array()

array of time instances for each window calculated

flstnp.ndarray(nfbins/2)

frequency array containing only positive frequencies where the Fourier coeffients were calculated

mtpy.processing.tf.wvd_analytic_signal(fx)[source]

Computes the analytic signal for WVVD as defined by J. M. O’ Toole, M. Mesbah, and B. Boashash, (2008), “A New Discrete Analytic Signal for Reducing Aliasing in the

Discrete Wigner-Ville Distribution”, IEEE Trans. on Signal Processing,

Argument:

fxnp.ndarray()

signal to compute anlytic signal for with length N

Returns:

fxanp.ndarray()

analytic signal of fx with length 2*N

Module contents

class mtpy.processing.AuroraProcessing(**kwargs)[source]

Bases: BaseProcessing

Convenience class to process with Aurora

from mtpy.processing.aurora.process_aurora import AuroraProcessing

ap = AuroraProcessing()

# set local station and path to MTH5
ap.local_station_id = "mt01"
ap.local_mth5_path = "/path/to/local_mth5.h5"

# set remote station and path to MTH5
ap.remote_station_id = "rr01"
ap.remote_mth5_path = "/path/to/remote_mth5.h5"

# process single sample rate
tf_obj = ap.process_single_sample_rate(sample_rate=1)

# process multiple sample rates, merge them all together and
# save transfer functions to the local MTH5
tf_processed_dict = ap.process(
    sample_rates=[4096, 1],
    merge=True,
    save_to_mth5=True
    ).
add_simple_coherence_weights(channel_list=[('ex', ['ex', 'hy']), ('ey', ['ey', 'hx']), ('hz', ['hx', 'hx'])], **kwargs) → list[ChannelWeightSpec][source]

Add coherence weights using the channel weight spec.

Parameters:

**kwargs (dict) – Additional keyword arguments (currently unused).

Returns:

List of channel weight specifications with coherence features.

Return type:

list[ChannelWeightSpec]

create_config(kernel_dataset: KernelDataset | None = None, decimation_kwargs: dict = {}, add_coherence_weights: bool = False, **kwargs) → Processing[source]

Create Aurora processing configuration.

Parameters:

  • kernel_dataset (KernelDataset or None, optional) – Kernel dataset defining processing runs, by default None.

  • decimation_kwargs (dict, optional) – Decimation parameters including window settings, by default {}.

  • add_coherence_weights (bool, optional) – Whether to add coherence-based weights, by default True.

  • **kwargs (dict) – Additional configuration parameters.

Returns:

Aurora configuration object.

Return type:

Processing

Raises:

ValueError – If kernel_dataset is None and no kernel dataset exists.

create_kernel_dataset(run_summary: RunSummary | None = None, local_station_id: str | None = None, remote_station_id: str | None = None, sample_rate: float | None = None) → KernelDataset[source]

Build and return a KernelDataset.

Parameters:

  • run_summary (RunSummary or None, optional) – Run summary to use, by default None (creates from MTH5).

  • local_station_id (str or None, optional) – Local station identifier, by default None.

  • remote_station_id (str or None, optional) – Remote reference station identifier, by default None.

  • sample_rate (float or None, optional) – Sample rate to filter runs, by default None.

Returns:

Kernel dataset defining processing configuration.

Return type:

KernelDataset

merge_transfer_functions(tf_dict: dict[float, dict[str, bool | MT]]) → MT[source]

Merge multiple transfer functions according to merge_dict.

Parameters:

tf_dict (dict) – Dictionary of transfer functions to merge.

Returns:

Merged transfer function combining all sample rates.

Return type:

MT

process(sample_rates: float | list[float] | None = None, processing_dict: dict[float, dict[str, Processing | KernelDataset]] | None = None, merge: bool = True, save_to_mth5: bool = True, plot: bool = False) → dict[float | str, dict[str, bool | MT]][source]

Process magnetotelluric data at multiple sample rates.

Parameters:

  • sample_rates (float, list of float, or None, optional) – Sample rate(s) to process, by default None.

  • processing_dict (dict or None, optional) – Dictionary mapping sample rates to config and kernel_dataset. Format: {sample_rate: {‘config’: Processing, ‘kernel_dataset’: KernelDataset}} By default None.

  • merge (bool, optional) – Whether to merge all sample rates into a single transfer function according to merge_dict, by default True.

  • save_to_mth5 (bool, optional) – Whether to save transfer functions to local MTH5 file, by default True.

Returns:

Dictionary with sample rates and ‘combined’ as keys, each containing {‘processed’: bool, ‘tf’: MT or None}.

Return type:

dict[float or str, dict[str, bool or MT]]

Raises:

  • ValueError – If neither sample_rates nor processing_dict is provided.

  • TypeError – If sample_rates or processing_dict is not the correct format.

Notes

If merge is True and multiple sample rates are processed, a ‘combined’ key is added with the merged transfer function.

Examples

>>> ap = AuroraProcessing()
>>> ap.local_station_id = "mt01"
>>> ap.local_mth5_path = "data.h5"
>>> results = ap.process(sample_rates=[1, 4], merge=True)
process_single_sample_rate(sample_rate: float, config: Processing | None = None, kernel_dataset: KernelDataset | None = None, plot: bool = False) → MT | None[source]

Process a single sample rate to generate transfer functions.

Parameters:

  • sample_rate (float) – Sample rate of time series data to process.

  • config (Processing or None, optional) – Aurora configuration object, by default None (creates from kernel_dataset).

  • kernel_dataset (KernelDataset or None, optional) – Kernel dataset defining processing runs, by default None (creates from run summary).

Returns:

Transfer function object, or None if processing fails.

Return type:

MT or None