PyHydroGeophysX.data_processing.mt package

Submodules

PyHydroGeophysX.data_processing.mt.analysis module

What a site’s transfer functions say about the ground: dimensionality and depth.

• phase_tensor() - Caldwell et al. (2004): Phi = X^-1 Y for Z = X + iY is free of galvanic distortion; its principal phases, skew beta, angle alpha, strike alpha - beta and ellipticity. |beta| > 3 deg is the usual sign of a 3D response (Booker 2014).

• swift_skew(), bahr_skew(), swift_strike() - the classical rotational invariants (Swift 1967; Bahr 1988).

• induction_arrows() - the tipper as arrows, in the Parkinson convention (real arrows point towards conductors) or Wiese’s.

• niblett_bostick() - an approximate resistivity-depth profile from apparent resistivity and phase (Bostick 1977; Jones 1983).

• estimate_static_shift(), static_shift_from_layers(), apply_static_shift() - static shift: the frequency-independent factor small near-surface bodies multiply apparent resistivity by, found against a reference (often a TEM sounding at the site, Sternberg et al. 1988) and divided out.

Caldwell, T. G., Bibby, H. M. & Brown, C. (2004). The magnetotelluric phase tensor. Geophysical Journal International, 158(2), 457-469. https://doi.org/10.1111/j.1365-246X.2004.02281.x

Booker, J. R. (2014). The magnetotelluric phase tensor: a critical review. Surveys in Geophysics, 35(1), 7-40. https://doi.org/10.1007/s10712-013-9234-2

Swift, C. M. (1967). A magnetotelluric investigation of an electrical conductivity anomaly in the southwestern United States. PhD thesis, MIT.

Bahr, K. (1988). Interpretation of the magnetotelluric impedance tensor: regional induction and local telluric distortion. Journal of Geophysics, 62(2), 119-127.

Bostick, F. X. (1977). A simple almost exact method of MT analysis. Workshop on Electrical Methods in Geothermal Exploration, U.S. Geological Survey.

Jones, A. G. (1983). On the equivalence of the “Niblett” and “Bostick” transformations in the magnetotelluric method. Journal of Geophysics, 53(1), 72-73.

Sternberg, B. K., Washburne, J. C. & Pellerin, L. (1988). Correction for the static shift in magnetotellurics using transient electromagnetic soundings. Geophysics, 53(11), 1459-1468. https://doi.org/10.1190/1.1442426

PyHydroGeophysX.data_processing.mt.analysis.apply_static_shift(tf: TransferFunction, factor_x: float = 1.0, factor_y: float | None = None) → TransferFunction

Divide the apparent resistivities of the Ex row by factor_x and Ey’s by factor_y.

factor_y defaults to factor_x. The phases are unchanged; errors and covariances scale with the rows.

PyHydroGeophysX.data_processing.mt.analysis.bahr_skew(tf: TransferFunction) → ndarray

Bahr’s phase-sensitive skew; above about 0.3 the response is 3D (Bahr 1988).

PyHydroGeophysX.data_processing.mt.analysis.estimate_static_shift(tf: TransferFunction, reference_rho_a: Any, *, period_range: Tuple[float, float] | None = None) → Dict[str, float]

Static-shift factors of the Ex and Ey rows against a reference apparent resistivity.

reference_rho_a is an array at tf.frequency or a callable of frequency. The factor is the geometric mean of rho_a / reference over period_range (s), where the reference holds - at the short periods a TEM sounding sees. Returns {"x", "y", "x_spread", "y_spread", "n"}.

PyHydroGeophysX.data_processing.mt.analysis.induction_arrows(tf: TransferFunction, *, convention: str = 'parkinson') → Dict[str, ndarray]

Real and imaginary induction arrows: lengths and azimuths (deg clockwise from north).

"parkinson" reverses the real arrow so it points at current concentrations; "wiese" keeps the tipper’s sign.

PyHydroGeophysX.data_processing.mt.analysis.niblett_bostick(tf: TransferFunction, *, mode: str = 'det') → Dict[str, ndarray]

depth (m) and resistivity (ohm m) by the Niblett-Bostick transform.

mode is "det" (the determinant impedance), "xy" or "yx". rho_NB = rho_a (pi / (2 phi) - 1), depth = sqrt(rho_a / (omega mu0)).

PyHydroGeophysX.data_processing.mt.analysis.phase_tensor(tf: TransferFunction, *, n_realizations: int = 0, seed: int = 0) → Dict[str, ndarray]

The phase tensor and its invariants (degrees), per frequency.

Returns phi (n, 2, 2) (in the TF’s frame) and phi_max, phi_min, alpha, beta (skew), strike and ellipticity; alpha and strike are azimuths clockwise from north, modulo 180 (strike has the usual 90-degree ambiguity). With n_realizations, their standard errors (*_err) come from perturbing Z by its errors.

PyHydroGeophysX.data_processing.mt.analysis.static_shift_from_layers(tf: TransferFunction, resistivity: Any, thickness: Any, *, period_range: Tuple[float, float] | None = None) → Dict[str, float]

Static-shift factors against the MT response of a layered model, such as a TEM inversion’s.

The model’s apparent resistivity is computed at the site’s frequencies (forward1d.apparent_resistivity_1d()); period_range should keep to the periods whose skin depth the model resolves.

PyHydroGeophysX.data_processing.mt.analysis.swift_skew(tf: TransferFunction) → ndarray

|Zxx + Zyy| / |Zxy - Zyx|; near 0 for 1D and 2D (Swift 1967).

PyHydroGeophysX.data_processing.mt.analysis.swift_strike(tf: TransferFunction) → ndarray

The strike (deg clockwise from north, +-90 ambiguity) that minimizes |Zxx|^2 + |Zyy|^2.

PyHydroGeophysX.data_processing.mt.edi module

SEG MT/EMAP Electrical Data Interchange (EDI) files: read and write.

The format is Wight’s 1987 SEG standard. A file is a sequence of blocks, each a >KEYWORD, an option list NAME=value and, for data blocks, a data set //count followed by the values; >=KEYWORD opens a section and >!...! is a comment. Files written by real instruments and programs bend it, so this reader takes keywords and option names in any case, options with spaces after = (ID=    11.001), // 80 with a space, and values written with a Fortran D exponent.

What a file can hold, in the order the reader prefers it:

• the impedance >ZXXR/>ZXXI …, with >ZXX.VAR (or .R.VAR and .I.VAR) and the tipper >TXR.EXP … with >TXVAR.EXP;

• apparent resistivities and phases only, >RHOXY/>PHSXY with .ERR, turned into the impedance they imply. Programs commonly write PHSYX folded into the first quadrant; it is moved back to the third, and the site says so;

• power spectra (>=SPECTRASECT), from which the impedance and tipper are estimated, Z = <E R*> <H R*>^-1 with the section’s reference channels.

Impedances are in (mV/km)/nT and are converted to ohms. Values at or beyond the file’s EMPTY (1.0E32) are missing. The standard’s time dependence is e^{+i omega t}; a file whose phases sit in the quadrants of the opposite convention throughout is conjugated, and the site records it.

PyHydroGeophysX.data_processing.mt.edi.is_edi_file(path: Any) → bool

Whether path is an EDI file: a .edi name, or a text that opens with >HEAD.

PyHydroGeophysX.data_processing.mt.edi.read_edi(path: Any, *, sign_convention: str = 'auto') → TransferFunction

Read an EDI file into a TransferFunction.

sign_convention is "+" (the standard’s e^{+i omega t}), "-" (conjugate everything), or "auto": "+" unless the Zxy and Zyx phases lie in the quadrants of e^{-i omega t} at four frequencies in five.

PyHydroGeophysX.data_processing.mt.edi.write_edi(tf: TransferFunction, path: Any) → Path

Write tf as a standard EDI file, impedances in (mV/km)/nT, e^{+i omega t}.

The impedance blocks carry Z**.VAR as the squared standard error, the tipper the experimental T*.EXP blocks every program reads, and both rotation blocks the site’s rotation.

PyHydroGeophysX.data_processing.mt.emtf_xml module

EMTF XML transfer functions (Kelbert 2020): read and write.

The format of the EarthScope/IRIS SPUD archive. Each <Period> holds the impedance <Z>, the tipper <T>, their variances <Z.VAR>/<T.VAR> and, from EMTF, the full error model: the inverse signal power <Z.INVSIGCOV> (inputs Hx, Hy) and the residual covariance <Z.RESIDCOV> (outputs). The file states its sign convention (<SignConvention>exp(+ i\omega t)</SignConvention>) and its units ([mV/km]/[nT], [V/m]/[T] or ohms); both are converted to this package’s. The data are in the frame <Orientation angle_to_geographic_north=...>.

Kelbert, A. (2020). EMTF XML: New data interchange format and conversion tools for electromagnetic transfer functions. Geophysics, 85(1), F1-F17. https://doi.org/10.1190/geo2018-0679.1

Kelbert, A., Erofeeva, S., Trabant, C., Karstens, R., & Van Fossen, M. (2018). Taking magnetotelluric data out of the drawer. Eos, 99. https://doi.org/10.1029/2018EO112859

PyHydroGeophysX.data_processing.mt.emtf_xml.is_emtf_xml_file(path: Any) → bool

PyHydroGeophysX.data_processing.mt.emtf_xml.read_emtf_xml(path: Any) → TransferFunction

Read an EMTF XML file into a TransferFunction.

PyHydroGeophysX.data_processing.mt.emtf_xml.write_emtf_xml(tf: TransferFunction, path: Any) → Path

Write tf as an EMTF XML file, impedances in (mV/km)/nT, e^{+i omega t}.

The full error model is written when the site carries it; otherwise only the variances, err^2. One rotation angle is written for the whole site, so a site whose frequencies are in different frames is first rotated to the frame of its first frequency.

PyHydroGeophysX.data_processing.mt.forward1d module

The impedance of a layered earth: the 1D MT forward problem.

For layers of resistivity rho_k and thickness h_k over a half-space, the impedance at the surface follows from the bottom up (Wait’s recursion):

k_j = sqrt(i omega mu0 / rho_j), Z_N = i omega mu0 / k_N, Z_j = z_j (Z_{j+1} + z_j tanh(k_j h_j)) / (z_j + Z_{j+1} tanh(k_j h_j)),

with z_j = i omega mu0 / k_j the layer’s intrinsic impedance, in ohms and e^{+i omega t}, so a half-space has rho_a = rho and a phase of 45 degrees. The derivative of log rho_a and the phase with respect to log rho_j is the same recursion differentiated, for the inversion.

Wait, J. R. (1954). On the relation between telluric currents and the Earth’s magnetic field. Geophysics, 19(2), 281-289. https://doi.org/10.1190/1.1437994

PyHydroGeophysX.data_processing.mt.forward1d.apparent_resistivity_1d(resistivity: Any, thickness: Any, frequency: Any) → Tuple[ndarray, ndarray]

(rho_a, phase_deg) of a layered earth.

PyHydroGeophysX.data_processing.mt.forward1d.impedance_1d(resistivity: Any, thickness: Any, frequency: Any) → ndarray

Surface impedance (ohm, e^{+i omega t}) of a layered earth at each frequency.

PyHydroGeophysX.data_processing.mt.forward1d.sensitivity_1d(resistivity: Any, thickness: Any, frequency: Any) → Tuple[ndarray, ndarray, ndarray]

Z and the derivatives of log10 rho_a and phase (rad) to log rho_j.

Returns (Z, d_log10_rho_a, d_phase), the latter two (n_freq, n_layers).

PyHydroGeophysX.data_processing.mt.inversion1d module

Occam’s 1D inversion of MT soundings, alone or jointly with a TEM sounding.

The model is log10 resistivity in many thin layers - thicknesses growing geometrically from a twentieth of the shortest period’s skin depth to beyond the longest’s - and Occam’s inversion (Constable et al. 1987) looks for the smoothest such model (least first-difference roughness) that fits the data to a target RMS misfit. Each iteration linearizes the forward problem (forward1d.sensitivity_1d() gives the exact derivatives), solves the regularized normal equations for a sweep of trade-off parameters mu, and takes the mu with the least misfit while the target is out of reach, then the largest mu that still reaches it.

The data are log10 rho_a and phase (degrees) of the determinant, of Zxy, of Zyx, or of both modes against one model. A static shift per mode can be estimated with the model: a free offset of that mode’s log10 rho_a, held to zero by a Gaussian prior of static_shift_prior (log10 units); the phases fix the model, the TEM (or the prior) the offset.

Joined with a TEM sounding (tem=), the transient’s response is fitted by the same model through PyHydroGeophysX’s TDEM forward operator - the TEM resolves the near surface that sets the static shift, the MT the depths below (Meju 1996). That needs SimPEG, as TEM inversion does; MT alone needs only NumPy.

Constable, S. C., Parker, R. L. & Constable, C. G. (1987). Occam’s inversion: a practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics, 52(3), 289-300. https://doi.org/10.1190/1.1442303

Meju, M. A. (1996). Joint inversion of TEM and distorted MT soundings: some effective practical considerations. Geophysics, 61(1), 56-65. https://doi.org/10.1190/1.1443956

class PyHydroGeophysX.data_processing.mt.inversion1d.Anisotropic1DResult(strike: float, resistivity: ~numpy.ndarray, thickness: ~numpy.ndarray, rms: float, mu: float, static_shift: ~typing.Dict[str, float], predicted: ~typing.Dict[str, ~typing.Dict[str, ~numpy.ndarray]], sounding: ~PyHydroGeophysX.data_processing.mt.inversion1d.Sounding1D, iterations: int, history: ~typing.List[~typing.Dict[str, float]] = <factory>, converged: bool = False)

Bases: object

Occam’s azimuthally anisotropic layered model: two resistivities per layer on one strike.

resistivity[0] is each layer’s resistivity along the axis at azimuth strike (degrees clockwise from north), which Zxy fits in that frame; resistivity[1] along the axis a quarter turn on, which Zyx fits.

property anisotropy: ndarray

Each layer’s resistivity across the strike over along it.

axis(k: int) → Occam1DResult

One axis’s model as an Occam1DResult, for the plots and conversions that take one.

converged: bool = False

property depth: ndarray

Each layer’s top, m (the half-space last).

history: List[Dict[str, float]]

iterations: int

mu: float

predicted: Dict[str, Dict[str, ndarray]]

resistivity: ndarray

rms: float

sounding: Sounding1D

static_shift: Dict[str, float]

strike: float

thickness: ndarray

class PyHydroGeophysX.data_processing.mt.inversion1d.Occam1DResult(resistivity: 'np.ndarray', thickness: 'np.ndarray', rms: 'float', roughness: 'float', mu: 'float', static_shift: 'Dict[str, float]', predicted: 'Dict[str, np.ndarray]', sounding: 'Sounding1D', iterations: 'int', history: 'List[Dict[str, float]]'=<factory>, tem: 'Optional[Dict[str, Any]]'=None, converged: 'bool' = False)

Bases: object

converged: bool = False

property depth: ndarray

Each layer’s top, m (the half-space last).

depth_profile() → Tuple[ndarray, ndarray]

Step-plot arrays (depth, resistivity), the half-space drawn one layer deep.

history: List[Dict[str, float]]

iterations: int

mu: float

predicted: Dict[str, ndarray]

resistivity: ndarray

rms: float

roughness: float

sounding: Sounding1D

static_shift: Dict[str, float]

tem: Dict[str, Any] | None = None

thickness: ndarray

class PyHydroGeophysX.data_processing.mt.inversion1d.Sounding1D(frequency: ndarray, modes: List[str], log_rho: ndarray, phase: ndarray, log_rho_err: ndarray, phase_err: ndarray, station: str = '')

Bases: object

What a 1D inversion fits: per mode, log10 rho_a and phase with errors.

frequency: ndarray

log_rho: ndarray

log_rho_err: ndarray

modes: List[str]

property n_data: int

phase: ndarray

phase_err: ndarray

station: str = ''

PyHydroGeophysX.data_processing.mt.inversion1d.layer_thicknesses(frequency: Any, n_layers: int = 50, *, resistivity_guess: Any = 100.0, depth_range: Tuple[float, float] | None = None) → ndarray

Geometric thicknesses from a twentieth of the least skin depth to twice the greatest.

resistivity_guess is one resistivity or a pair: the apparent resistivity at the highest and at the lowest frequency, which set the shallowest and the deepest skin depth. One value for both made the top layer of a resistive site with a conductive cover tens of metres thick.

PyHydroGeophysX.data_processing.mt.inversion1d.occam1d(data: Any, *, mode: str = 'det', n_layers: int = 50, thickness: Sequence[float] | None = None, target_rms: float = 1.0, max_iterations: int = 30, starting_resistivity: float | None = None, static_shift: bool = False, static_shift_prior: float = 0.3, error_floor: float = 0.05, period_range: Tuple[float, float] | None = None, tem: Dict[str, Any] | None = None, tem_weight: float = 1.0, log: Callable[[str], None] | None = None, resistivity_bounds: Tuple[float, float] = (0.01, 100000.0)) → Occam1DResult

Occam’s smoothest model fitting a site’s 1D data (and a TEM sounding).

data is a TransferFunction (its mode is taken, see sounding_from_tf()) or a Sounding1D. thickness fixes the layers; otherwise n_layers are spread over the skin-depth range. static_shift adds one offset of log10 rho_a per mode. tem={"data": ..., "geometry": ..., "inversion": ...} is a TDEM sounding as PyHydroGeophysX.inversion.em1d.tdem_invert() takes it, fitted with weight tem_weight; the RMS target then applies to all data together.

PyHydroGeophysX.data_processing.mt.inversion1d.occam1d_anisotropic(data: TransferFunction, *, strike: float | None = None, n_layers: int = 50, thickness: Sequence[float] | None = None, target_rms: float = 1.0, max_iterations: int = 30, starting_resistivity: float | None = None, static_shift: bool = False, static_shift_prior: float = 0.3, anisotropy_weight: float = 1.0, error_floor: float = 0.05, period_range: Tuple[float, float] | None = None, log: Callable[[str], None] | None = None, resistivity_bounds: Tuple[float, float] = (0.01, 100000.0)) → Anisotropic1DResult

Occam’s smoothest azimuthally anisotropic layered model with one strike.

Each layer conducts differently along two horizontal axes, the same two in every layer: along strike (degrees clockwise from north; by default the phase tensor’s, averaged over the periods used) and across it. In that frame the layered earth’s impedance has no diagonal, and Zxy and Zyx are each a 1D isotropic response - Zxy of the resistivities along the strike, Zyx of those across it (Pek & Santos 2002, the case of horizontal principal axes sharing one azimuth). The two models share the layers and are found together: Occam’s roughness of each, plus anisotropy_weight times the difference of their log10 resistivities layer by layer, so a layer is made anisotropic only where the data ask for it.

Why the modes differ cannot be told from one site: a 2D structure with the same strike, or galvanic distortion, splits them too. A static shift moves rho_a and not the phase, so with static_shift the difference in level between the modes is put down to the shifts and the anisotropy comes from the phases’ split alone.

Pek, J. & Santos, F. A. M. (2002). Magnetotelluric impedances and parametric sensitivities for 1-D anisotropic layered media. Computers & Geosciences, 28(8), 939-950. https://doi.org/10.1016/S0098-3004(02)00014-6

PyHydroGeophysX.data_processing.mt.inversion1d.sounding_from_tf(tf: TransferFunction, *, mode: str = 'det', error_floor: float = 0.05, period_range: Tuple[float, float] | None = None) → Sounding1D

The 1D data of a site: apparent resistivity and phase of a mode, with errors.

error_floor is the least relative error of the impedance (5%: about 0.043 in log10 rho_a and 1.4 degrees of phase). Phases are folded into the first quadrant; values outside it are dropped.

PyHydroGeophysX.data_processing.mt.inversion1d.water_content_profile(result: Occam1DResult, *, rhos: Any, n: Any, porosity: Any, sigma_sur: Any = 0.0) → Dict[str, ndarray]

Volumetric water content of each layer from its resistivity (Waxman-Smits).

The petrophysical link the rest of PyHydroGeophysX uses (PyHydroGeophysX.petrophysics.resistivity_models.resistivity_to_water_content()): rhos the saturated resistivity, n the saturation exponent, porosity and sigma_sur the surface conductivity - each one value or one per layer.

PyHydroGeophysX.data_processing.mt.inversion2d module

2D MT forward modelling and inversion of a profile, on SimPEG’s NSEM simulations.

A profile of sites across a 2D structure is modelled on a tensor mesh in the vertical plane through it: the profile horizontal, the strike out of the plane, air above the surface. SimPEG (Cockett et al. 2015) solves Maxwell’s equations there for the two polarizations - TE, the electric field along strike, which is Zxy in the strike frame, and TM, the electric field across strike, Zyx - with 1D fields as the side boundary conditions.

invert_profile() rotates each site’s impedance to the strike, fits the real and imaginary parts of the chosen modes (errors: the impedance’s own, with a relative floor) for the logarithm of conductivity below the surface, with a smoothness regularization whose weight is cooled until the misfit reaches its target (Gauss-Newton, conjugate gradients). The mesh is laid out from the frequencies: cells a tenth of the least skin depth at the top, growing to three times the greatest below and beside the profile.

SimPEG works in east, north, up; here the strike is SimPEG’s y. Its e^{+i omega t} impedances agree with this package’s, and on a layered earth both modes reproduce forward1d.impedance_1d() to within the mesh’s discretization: about 1-2% over most of the band, rising to a few per cent at the highest frequency (9% for TE at 1 kHz over 30 ohm m on the default mesh; a finer top_cell halves that but not the TM error, which depends on how the surface fields are interpolated). On a synthetic block the inversion recovers 12 ohm m for 10 inside it and 96 for 100 beside it.

Cockett, R., Kang, S., Heagy, L. J., Pidlisecky, A. & Oldenburg, D. W. (2015). SimPEG: an open source framework for simulation and gradient based parameter estimation in geophysical applications. Computers & Geosciences, 85, 142-154. https://doi.org/10.1016/j.cageo.2015.09.015

deGroot-Hedlin, C. & Constable, S. (1990). Occam’s inversion to generate smooth, two-dimensional models from magnetotelluric data. Geophysics, 55(12), 1613-1624. https://doi.org/10.1190/1.1442813

PyHydroGeophysX.data_processing.mt.inversion2d.PROFILE_FREQUENCY_TOLERANCE = 0.05

A site frequency within 5 % of a profile frequency (|ln(f_site / f)| below this) is taken as that frequency, as it always was.

PyHydroGeophysX.data_processing.mt.inversion2d.PROFILE_INTERPOLATION_SPAN = 0.8109302162163288

Widest gap between two of a site’s own frequencies across which its impedance is interpolated to a profile frequency: a factor of 2.25, so a band sampled three times per decade still interpolates. Midway across a gap that wide the interpolation is off by at most 2.3 % on three test earths (100/10/300, 1000/1/1000 and 10/1000 ohm m), under half the 5 % error floor; at 5 per decade by 0.8 %, at 10 by 0.2 %.

class PyHydroGeophysX.data_processing.mt.inversion2d.ProfileInversionResult(profile: 'ProfileMesh', resistivity: 'np.ndarray', frequencies: 'np.ndarray', modes: 'List[str]', observed: 'Dict[str, np.ndarray]', predicted: 'Dict[str, np.ndarray]', uncertainty: 'Dict[str, np.ndarray]', rms: 'float', history: 'List[Dict[str, float]]'=<factory>, data_per_site: 'List[int]' = <factory>, warnings: 'List[str]' = <factory>)

Bases: object

data_per_site: List[int]

Impedance values each site gave the fit (one per frequency and mode).

frequencies: ndarray

history: List[Dict[str, float]]

modes: List[str]

observed: Dict[str, ndarray]

predicted: Dict[str, ndarray]

profile: ProfileMesh

resistivity: ndarray

rms: float

section() → Tuple[ndarray, ndarray, ndarray]

(x_nodes, z_nodes, resistivity grid) of the ground, for pcolormesh.

uncertainty: Dict[str, ndarray]

warnings: List[str]

A line for each site that gave fewer than the profile asks for.

class PyHydroGeophysX.data_processing.mt.inversion2d.ProfileMesh(mesh: Any, active: ndarray, station_x: ndarray)

Bases: object

A 2D tensor mesh for a profile: mesh (discretize) and its ground cells.

active: ndarray

property ground_centers: ndarray

mesh: Any

station_x: ndarray

PyHydroGeophysX.data_processing.mt.inversion2d.build_profile_mesh(station_x: Sequence[float], frequencies: Sequence[float], *, resistivity: float = 100.0, core_cell: float | None = None, top_cell: float | None = None, vertical_growth: float = 1.12, padding_growth: float = 1.4) → ProfileMesh

A mesh fitted to the profile and the band.

core_cell (default: half the least station spacing, at most a fifth of the least skin depth) spans the stations; top_cell (default: a tenth of the least skin depth) starts the ground, which grows by vertical_growth to three greatest skin depths; padding grows by padding_growth beyond the profile as far, and the air above as high.

PyHydroGeophysX.data_processing.mt.inversion2d.forward_profile(profile: ProfileMesh, resistivity: Any, frequencies: Sequence[float], *, station_x: Sequence[float] | None = None, modes: Sequence[str] = ('te', 'tm')) → Dict[str, ndarray]

TE and TM impedances (ohm, strike frame) at the stations: {mode: (n_freq, n_sta)}.

resistivity is one value per ground cell (profile.active), or a callable of the ground cells’ centres (x, z).

PyHydroGeophysX.data_processing.mt.inversion2d.invert_profile(tfs: Sequence[TransferFunction], station_x: Sequence[float], *, strike: float = 0.0, modes: Sequence[str] = ('te', 'tm'), frequencies: Sequence[float] | None = None, error_floor: float = 0.05, starting_resistivity: float | None = None, profile: ProfileMesh | None = None, max_iterations: int = 20, alpha_s: float = 0.0001, beta0_ratio: float = 10.0, cooling_factor: float = 2.0, target_chifact: float = 1.0, seed: int = 0, log: Callable[[str], None] | None = None) → ProfileInversionResult

Smooth 2D resistivity section that fits a profile’s TE and/or TM impedances.

tfs are the sites’ transfer functions and station_x their distances along the profile (m); strike (degrees clockwise from north) is the strike the impedances are rotated to. frequencies defaults to those of the first site. A site recorded at other frequencies is interpolated to them within its own band (_site_impedance()); a site that gives fewer values than the profile asks for is named in the result’s warnings, with data_per_site counting what each gave.

PyHydroGeophysX.data_processing.mt.inversion2d.station_distances(tfs: Sequence[TransferFunction]) → Tuple[ndarray, float]

Each site’s distance along its profile (m) and the profile’s azimuth (deg).

The sites’ latitudes and longitudes are projected on a local plane (equirectangular, about their mean); the profile is their principal direction, and distances count from the site at its start.

PyHydroGeophysX.data_processing.mt.io module

One entry point for MT transfer-function files and one for instrument time series.

PyHydroGeophysX.data_processing.mt.io.TIMESERIES_FORMATS = ('phoenix', 'phoenix_legacy', 'metronix', 'zonge', 'lemi424')

The instrument formats read_timeseries() recognizes.

PyHydroGeophysX.data_processing.mt.io.TRANSFER_FUNCTION_PATTERNS = ('*.edi', '*.EDI', '*.xml', '*.zmm', '*.zrr', '*.zss', '*.j')

File patterns read_transfer_functions() collects from a folder.

PyHydroGeophysX.data_processing.mt.io.is_transfer_function_file(path: Any) → bool

PyHydroGeophysX.data_processing.mt.io.read_timeseries(path: Any, *, azimuths: Mapping[str, float] | None = None, **options: Any) → List[TimeSeriesRun]

Read an instrument’s MT time series - a file or a recording folder - into runs.

The format is recognized from the files: Phoenix MTU-5C/5P/8A recordings (read_phoenix()), legacy Phoenix MTU-5A TS/TBL sites (read_phoenix_legacy()), Metronix ATS (read_metronix()), Zonge Z3D (read_zonge()) and LEMI-424 text files (read_lemi424()). options go to that reader; format forces one.

azimuths ({"ex": 180, "hy": 270, ...}, degrees clockwise from north) replaces the direction the files give those channels, in every run: for a layout that differs from what the receiver was told, such as a dipole laid out to the south or a coil pointing west. Processing expresses every channel in the frame of Hx, so a reversed channel needs nothing more.

PyHydroGeophysX.data_processing.mt.io.read_transfer_function(path: Any, *, sign_convention: str = 'auto', z_units: str = 'auto') → TransferFunction

Read one site’s transfer functions from an EDI, EMTF XML, Z- or J-file.

sign_convention applies to EDI and J-files, which do not state it (see read_edi()), and z_units to J-files (see read_jfile()); EMTF XML and Z-files state both.

PyHydroGeophysX.data_processing.mt.io.read_transfer_functions(sources: str | Path | Iterable[Any], **options: Any) → List[TransferFunction]

Read every transfer-function file of a folder, or of a list of paths, in name order.

PyHydroGeophysX.data_processing.mt.io.timeseries_format(path: Any) → str | None

Which instrument wrote path (a file or a folder), or None.

PyHydroGeophysX.data_processing.mt.jfile module

Jones’s J-format (.j) transfer functions, as BIRRP writes them: read.

Comment lines start with # and keyword lines with > (>LATITUDE =, >AZIMUTH =); the first other line is the station name. Then come blocks, each a component name - ZXX..``ZYY`` impedances, RXX..``RYY`` apparent resistivities and phases, TZX/TZY tipper - with an optional units word, a line with the number of rows, and the rows: period, real part, imaginary part, error, … for Z and T; period, resistivity, phase, … for R. Missing values are written -999; a negative period is a frequency.

The units word is unreliable: BIRRP labels field-unit impedances S.I.. When the file also has R blocks, the impedance units are those that reproduce the file’s own apparent resistivities; otherwise field units, (mV/km)/nT, are assumed unless z_units says otherwise.

Jones, A. G. (1994). J-format: a format for the exchange of magnetotelluric transfer functions. MTNet, https://www.mtnet.info/docs/jformat.txt

PyHydroGeophysX.data_processing.mt.jfile.is_jfile(path: Any) → bool

PyHydroGeophysX.data_processing.mt.jfile.read_jfile(path: Any, *, z_units: str = 'auto', sign_convention: str = 'auto') → TransferFunction

Read a J-format file; z_units is "auto", "field" or "ohm".

BIRRP writes e^{-i omega t}; sign_convention is resolved as for EDI files (resolve_sign_convention()), so such a file is conjugated.

PyHydroGeophysX.data_processing.mt.lemi module

LEMI-424 long-period MT stations: the logger’s 1 S/s text files.

Each line is one second: year, month, day, hour, minute, second (UTC); Bx, By, Bz in nT from the fluxgate; the electronics’ and the sensor’s temperatures; E1-E4 in mV (the potential differences of up to four dipoles); battery voltage; elevation; latitude and longitude as ddmm.mmmm with a hemisphere letter; satellites; fix flag; time difference. Files usually hold a day (or an hour) and are joined in time.

A fluxgate measures B itself, flat across the MT band, so the magnetic channels are already in nT. The files do not hold the dipole lengths: pass dipole_lengths (m) to bring E1/E2 to mV/km; E1 is taken as Ex and E2 as Ey, E3 and E4 keep their names unless components renames them.

LEMI LLC (2009). LEMI-424 long-period magnetotelluric station: user manual. Lviv Centre of Institute of Space Research, Lviv, Ukraine.

PyHydroGeophysX.data_processing.mt.lemi.is_lemi424_file(path: Any) → bool

PyHydroGeophysX.data_processing.mt.lemi.read_lemi424(path: Any, *, dipole_lengths: Dict[str, float] | None = None, components: Iterable[str] | None = None, station: str = '', calibration: Dict[str, Any] | None = None, start: Any = None, end: Any = None) → List[TimeSeriesRun]

Read LEMI-424 text files - one, a folder of them, or a list - into runs.

dipole_lengths maps "ex"/"ey" (or "e1"…) to metres. components lists the channels to keep, by LEMI name or component (default: hx, hy, hz, ex, ey). calibration maps a component to a text table (frequency, amplitude, phase in degrees) of its response.

PyHydroGeophysX.data_processing.mt.metronix module

Metronix ADU-06/07/08 recordings: ATS files, the measurement XML and coil calibrations.

An ADU writes one .ats file per channel - 085_V01_C02_R001_THx_BL_2048H.ats, system 085, channel 2, run 1, Hx, band 2048 Hz - into a meas_<start> folder with a measurement XML. An ATS file (header versions 80 and 81) is a 1024-byte little-endian header followed by int32 samples, which times the header’s least significant bit give millivolts at the channel input (the gains are already in it).

The header fields read here: header length (0), version (2), samples (4), sample rate (8, float32), start in Unix seconds (12), LSB in mV (16, double), UTC offset (24), system serial (32), channel number (36), chopper (37), channel type (38, Ex…), sensor type (40) and serial (46), electrode positions x1, y1, z1, x2, y2, z2 in m (48-68, x north, y east), dipole length and angle (72, 76), latitude and longitude in milliseconds of arc (96, 100), elevation in cm (104) and the system type (132).

Coils (MFS-06e, MFS-07e) are calibrated as V/(nT Hz) against frequency, with the chopper on and off; the table for the channel’s chopper setting is used, from the measurement XML’s calibration_sensors or a Metronix calibration text file (MFS06E0005.txt), given as calibration or found in a cal folder next to the measurements.

Metronix Geophysics (2009). ADU-07 system manual, appendix: ATS file format.

PyHydroGeophysX.data_processing.mt.metronix.is_ats_file(path: Any) → bool

PyHydroGeophysX.data_processing.mt.metronix.read_ats_header(path: Any) → Dict[str, Any]

PyHydroGeophysX.data_processing.mt.metronix.read_metronix(path: Any, *, components: Iterable[str] | None = None, calibration: Any = None, start: Any = None, end: Any = None) → List[TimeSeriesRun]

Read Metronix ATS files - one file, a meas_ folder or a site folder - into runs.

Files of one measurement folder with the same run number and sample rate form a run. calibration is a coil calibration text file, a folder of them, or a list; cal folders beside and above the measurements are searched too.

PyHydroGeophysX.data_processing.mt.metronix.read_metronix_calibration(path: Any) → Dict[str, Any]

A Metronix coil calibration text file: {"sensor", "serial", "on", "off"} stages.

PyHydroGeophysX.data_processing.mt.phoenix module

Phoenix Geophysics MTU-5C, MTU-5P, MTU-8A, RXU-8A and MTU-2C recordings.

A recording is a folder <serial>_<yyyy-mm-dd-hhmmss> holding recmeta.json and config.json and one numbered folder per channel; the channel map in recmeta.json says which connector (E1, E2, H1, H2, H3…) each folder holds. Each folder has the channel’s files, fragmented in time:

• .bin - the native 24 kS/s stream: 64-byte frames of twenty 24-bit big-endian samples and a little-endian footer whose 28-bit counter numbers the frame from the start of the recording;

• .td_150, .td_30 - decimated continuous float32 volts;

• .td_24k, .td_2400 and other rates - decimated segmented files, bursts each with its own 32-byte subheader (GPS time stamp, sample count).

All of them start with the 128-byte header of Phoenix’s specification, which this reader follows; the hardware-configuration bits that set a board’s gains, which the specification leaves undocumented, are decoded in _board_gain().

Samples come out in volts at the instrument input: native counts are scaled by the A/D converter’s +-5 V over 2^23 and divided by the board’s total gain, which is what the decimated files already hold. The responses come from Phoenix’s calibration files: <receiver>_rxcal.json (each channel’s anti-alias filter, one table per low-pass setting) and <coil>_*.scal.json (each coil, mV/nT).

Phoenix time stamps count GPS seconds from 1970; they are converted to UTC with the receiver’s own leap-second count. Files written before firmware 2.0 (native version < 4, decimated < 3) stamp one second early, which is undone. The stamps do not include the decimation filters’ delay (about 15 ms at 150 S/s against the 24 kS/s bursts): every channel of a stream shares it, so it cancels in transfer functions, but streams of different rates are offset by it.

Phoenix Geophysics (2021). Time series file specifications for MTU-8A, RXU-8A, MTU-5C, MTU-2C and MTU-5D. Document DAA09, version 210915.

PyHydroGeophysX.data_processing.mt.phoenix.is_phoenix_recording(path: Any) → bool

PyHydroGeophysX.data_processing.mt.phoenix.phoenix_coil_calibrations(path: Any, calibration: Any = None) → List[Dict[str, Any]]

The coil calibrations a Phoenix recording would use: {"serial", "model", "file"} each.

calibration is as read_phoenix() takes it. One entry per serial, its newest file. Files that do not read as calibrations are left out.

PyHydroGeophysX.data_processing.mt.phoenix.read_phoenix(path: Any, *, sample_rate: None | float | Sequence[float] = None, components: Iterable[str] | None = None, start: Any = None, end: Any = None, calibration: Any = None, coils: Dict[str, str] | None = None) → List[TimeSeriesRun]

Read a Phoenix MTU-5C-family recording into runs.

path is the recording folder, one of its channel folders, or a file (which reads that file’s stream, every channel of it). sample_rate picks the stream(s) - 24000 for the native .bin or a decimated rate; by default every rate present is read, highest first. start and end (UTC) limit what is read, which matters for long native recordings. calibration is a calibration file, a folder of them (searched through) or a list; by default *rxcal.json, *scal.json and binary *.scal files are looked for in the recording folder and the two above it, and in any folder there named for calibrations.

coils says which coil was on which magnetic channel, {"H1": "4684", ...} (or "hx"…), by serial or calibration file. A receiver whose coils were not detected records serial 0, so without it their responses are not known and the magnetic channels stay in volts.

Continuous streams come out as one run per gap-free stretch, segmented ones as one run per burst.

PyHydroGeophysX.data_processing.mt.phoenix.read_phoenix_calibration(path: Any) → Dict[str, Any]

A Phoenix rxcal.json, scal.json or binary .scal: each channel’s response tables.

Returns {"kind": "receiver" | "sensor", "instrument_serial", "sensor_serial", "model", "channels": {tag: [ResponseStage, ...]}}; a receiver file has one table per low-pass setting.

PyHydroGeophysX.data_processing.mt.phoenix.read_phoenix_header(path: Any) → Dict[str, Any]

The 128-byte header of a Phoenix time-series file, decoded.

PyHydroGeophysX.data_processing.mt.phoenix_legacy module

Phoenix Geophysics legacy MTU-5, MTU-5A and V8 recordings: TS2-TS5 and TBL.

A site’s recording is a table file <name>.TBL and one time-series file per band, <name>.TS2 … <name>.TS5 (MTU-5A: 2400, 150 and 15 S/s for TS3-TS5; TS2 for AMT). A TS file is a sequence of records, each a 32-byte tag - the UTC time of its first scan (second, minute, hour, day, month, year, weekday, century), serial number, scans in the record, channels per scan, tag length, status and saturation flags, sample length, sample rate - followed by the scans: per channel a 24-bit signed little-endian integer. Records that follow each other in time form one run; the high bands are bursts.

The TBL file holds 25-byte records: a 12-byte name and a 13-byte value (integer, double, text or time). Its gains turn counts into volts at the box input, counts * FSCV / 2^23 / gain with gain EGN for E and HGN * HATT for H, and its EXLN/EYLN dipole lengths and HNOM coil sensitivity (mV/nT) give a flat response; Phoenix’s per-frequency box and coil calibrations (CLB/CLC) are not decoded, but any table can be given as calibration (see read_response_table()). Without them the magnetic channels are right only in the coil’s flat band, and only up to the coil’s phase there: on Phoenix’s sample site the impedance phases come out 180 degrees from their usual quadrants and the apparent resistivity climbs at long periods, where an induction coil’s output falls with frequency.

Phoenix Geophysics (2005). Instrument and sensor calibration: concepts and utility programs. Phoenix Geophysics Ltd., Toronto.

PyHydroGeophysX.data_processing.mt.phoenix_legacy.is_phoenix_legacy_file(path: Any) → bool

PyHydroGeophysX.data_processing.mt.phoenix_legacy.read_phoenix_legacy(path: Any, *, bands: Iterable[str] | None = None, components: Iterable[str] | None = None, start: Any = None, end: Any = None, calibration: Dict[str, Any] | None = None) → List[TimeSeriesRun]

Read a legacy Phoenix site (.TBL with .TS2-.TS5) into runs.

path is one of the site’s files. bands picks "TS3", "TS4"… (default: every TS file of the site). calibration maps a component to a ResponseStage or a text table of frequency, amplitude (mV/nT for coils), phase in degrees, which replaces the flat coil sensitivity.

PyHydroGeophysX.data_processing.mt.phoenix_legacy.read_phoenix_table(path: Any) → Dict[str, Any]

The records of a legacy Phoenix .TBL file, decoded by name.

PyHydroGeophysX.data_processing.mt.processing module

From time series to transfer functions: the impedance and tipper of a site.

process_mt() takes a site’s runs - any instrument, any number, any sample rates - and, optionally, a remote reference’s, and estimates the impedance E = Z H and the tipper Hz = T H band by band:

1. runs are grouped by sample rate; with a remote reference, each run is cut to the stretch the remote site recorded at the same rate;

2. spectra gives each run’s calibrated Fourier coefficients at the bands’ decimation levels, the remote’s on the same windows;

3. the coefficients are expressed in the frame of the Hx azimuth (an electric or magnetic pair that is not orthogonal, or turned from the other, is resolved into it);

4. for each band, robust.robust_regression() estimates each output’s row from every window’s coefficients at the band’s harmonics; the error covariance is scaled by spectra.variance_inflation(), because a band’s neighbouring harmonics are correlated through the taper (with tapers above one, every Slepian taper’s coefficients are points);

5. bands from different sample rates are merged, the better estimate kept where two overlap.

Before step 1 every run can be rid of power-line hum (powerline) and of interference repeating on the receiver’s clock (periodic_noise).

The defaults are EMTF’s (Egbert 1997): 128-point windows overlapping by 32, decimation by 4, its band set, Huber then redescending weights.

Egbert, G. D. (1997). Robust multiple-station magnetotelluric data processing. Geophysical Journal International, 130(2), 475-496. https://doi.org/10.1111/j.1365-246X.1997.tb05663.x

Egbert, G. D. & Livelybrooks, D. W. (1996). Single station magnetotelluric impedance estimation: Coherence weighting and the regression M-estimate. Geophysics, 61(4), 964-970. https://doi.org/10.1190/1.1444045 (presorting)

Ernst, T., Sokolova, E. Yu., Varentsov, Iv. M. & Golubev, N. G. (2001). Comparison of two techniques for magnetotelluric data processing using synthetic data sets. Acta Geophysica Polonica, 49(2), 213-243 (multi-window stacking)

class PyHydroGeophysX.data_processing.mt.processing.ProcessingConfig(window_length: int = 128, overlap: int = 32, decimation_factor: int = 4, n_levels: int | None = None, min_windows: int = 5, bands: List[Band] | None = None, bands_per_decade: int | None = None, taper: str = 'hamming', prewhiten: bool = True, huber: float = 1.5, redescend: float = 2.8, max_iterations: int = 10, redescend_iterations: int = 2, tolerance: float = 0.005, leverage_cutoff: float | None = None, min_points: int = 10, allow_uncalibrated: bool = False, tipper: bool = True, periodic_noise: Tuple[float, ...] = (), periodic_cycles: int = 31, periodic_where_found: bool = True, skip_notched: bool = True, widen_above: float = 0.0, powerline: float | None = None, drop_saturated: bool = True, bounded_influence: bool = False, jackknife: bool = False, min_coherence: float = 0.0, presort_windows: int = 10, min_reference_coherence: float = 0.0, tapers: int = 1, multitaper_above: float = 0.0, multitaper_tapers: int = 3, multitaper_periods: Tuple[float, float] | None = None, multi_window: bool = False, noise_classes: int = 0)

Bases: object

Every choice of process_mt(); the defaults are EMTF’s.

allow_uncalibrated: bool = False

bands: List[Band] | None = None

bands_per_decade: int | None = None

Overlapping bands, this many to a decade (spectra.log_bands()), in place of EMTF’s band set: continuous curves, errors correlated between neighbours.

bounded_influence: bool = False

Weight points by their leverage too, every iteration (robust): noise in the magnetic channels then cannot pull the estimate.

decimation_factor: int = 4

drop_saturated: bool = True

Leave out the windows holding a sample at the recorder’s full scale (timeseries.saturated()), where the reader knows it.

huber: float = 1.5

jackknife: bool = False

Errors from the delete-a-group jackknife (20 groups) instead of the parametric model (robust).

leverage_cutoff: float | None = None

max_iterations: int = 10

min_coherence: float = 0.0

in each band, stretches of presort_windows consecutive windows whose multiple coherence (squared) of the output with Hx, Hy is below this are left out of that output’s regression.

Type:

Coherence presorting

min_points: int = 10

min_reference_coherence: float = 0.0

where a band’s Hx or Hy is coherent with the remote’s same channel below this (squared coherence), the reference holds almost no natural signal; the elements estimated from that input (Hx: Zxx, Zyx; Hy: Zxy, Zyy; and the tipper’s) are then a ratio of cross-spectra near zero, and the jackknife understates their error. Their relative errors are raised to at least 100 %. 0 is off.

Type:

Remote reference

min_windows: int = 5

multi_window: bool = False

each band also takes the coefficients at its frequencies of the windows one level up or down - windows decimation_factor times longer or shorter - in one regression. Those overlap the band’s own data in time, so the errors are the jackknife’s.

Type:

Multi-window stacking (Ernst et al. 2001)

multitaper_above: float = 0.0

An output (Ex for Zxy, Ey for Zyx, Hz) whose relative error in a band is above this with one taper is estimated again there with multitaper_tapers Slepian tapers, and that estimate kept. Where the field is weakest the tapers’ extra points cut the scatter; where the errors are a few per cent, the tapers’ smoothing over neighbouring harmonics (a bias of a few per cent where the curve is steep, which the errors do not show) would be the larger error, so those bands keep one taper. Where the two estimates differ, which is right is not known: the difference is added to the errors (in quadrature), so the curve is the steadier one and its errors hold the other. 0 is off.

multitaper_periods: Tuple[float, float] | None = None

The periods (s, lowest and highest) whose bands multitaper_above may estimate again; None is every band. Tapers help where few windows are the limit, as in a weak dead band; in short bursts at high rates they can make the curve rougher.

multitaper_tapers: int = 3

n_levels: int | None = None

noise_classes: int = 0

up to this many classes of noise that grows with the signal, each normalised before the M-estimate; 0 is off.

Type:

Noise classes (robust)

overlap: int = 32

periodic_cycles: int = 31

fewer follow interference whose shape drifts, each notch then a little wider.

Type:

Cycles each interference waveform is averaged over

periodic_noise: Tuple[float, ...] = ()

Periods (s) of interference repeating exactly on the receiver’s clock, subtracted from every run before its spectra (timeseries.remove_periodic()).

periodic_where_found: bool = True

Remove a period from a run only where timeseries.periodic_interference() finds it in that run (when the run holds 30 cycles or more, enough to tell). The removal also takes the field at every harmonic of the period, and where a level’s windows cannot resolve those notches it biases the band: so one site’s interference is not taken out of another site’s data. Bursts too short to fold are cleaned for every period, as before.

powerline: float | None = None

Grid frequency (50 or 60 Hz) whose hum and harmonics are notched out of every run before its spectra (timeseries.remove_powerline()); None is off.

presort_windows: int = 10

prewhiten: bool = True

redescend: float = 2.8

redescend_iterations: int = 2

regression() → RegressionConfig

skip_notched: bool = True

Leave out of the bands the harmonics inside those notches, where windows are long enough to resolve them (twice the period or more): the notch takes the field there too, and what is left is the interference’s drift.

spectra() → SpectraConfig

taper: str = 'hamming'

tapers: int = 1

Slepian tapers per window (multitaper, spectra); 1 is taper alone.

tipper: bool = True

tolerance: float = 0.005

widen_above: float = 0.0

A band whose Zxy and Zyx relative errors average above this is estimated again from a band about twice as wide (Q 1.5 instead of 3, same level and centre) and the wider one kept if its errors are smaller: more data where the field is weakest, at the price of resolution. 0 is off.

window_length: int = 128

PyHydroGeophysX.data_processing.mt.processing.merge_estimates(estimates: List[Dict[str, Any]]) → TransferFunction

One TransferFunction from the bands of every sample rate.

Where two rates cover the same frequencies, the one whose impedance errors are smaller there, as a whole, takes the stretch - so the curve does not hop between rates band by band. Bands that still fall within 10% in frequency keep the one with the smaller errors.

PyHydroGeophysX.data_processing.mt.processing.process_mt(runs: TimeSeriesRun | Sequence[TimeSeriesRun], *, remote: TimeSeriesRun | Sequence[TimeSeriesRun] | None = None, config: ProcessingConfig | None = None, station: str | None = None, log: Callable[[str], None] | None = None, **options: Any) → TransferFunction

Estimate a site’s impedance (and tipper, when it has Hz) from its runs.

remote holds the reference site’s runs (its Hx and Hy are used). config or keyword options set ProcessingConfig fields; log receives a line per sample rate.

PyHydroGeophysX.data_processing.mt.robust module

Robust regression of one output’s Fourier coefficients on two inputs.

The transfer function b in y = X b + r - y one band’s Ex (or Ey, Hz) coefficients, X the Hx, Hy ones - is estimated as the instrumental variable b = (R^H W X)^{-1} R^H W y, with R = X for a single site and the remote site’s Hx, Hy for a remote reference (Gamble et al. 1979), and W diagonal weights that a regression M-estimator sets from the residuals (Egbert & Booker 1986; Chave et al. 1987):

1. least squares, then iterations with Huber weights min(1, c / |r_i|/s) (c = 1.5) until b changes by less than tolerance; the scale s is re-estimated every time from the winsorized squared residuals, divided by 1 - exp(-c^2), their expectation for complex Gaussian residuals;

2. a few iterations with the redescending weights exp(-exp(u0 (|r_i|/s - u0))) (u0 = 2.8), which reject outliers outright, with the scale correction computed for those weights;

3. optionally, bounded influence, after Chave & Thomson (2004): the M-estimate bounds the influence of a large residual but not of a large input - a burst of noise in the magnetic channels gives its points a large leverage and they pull the estimate however small their residuals. With bounded_influence, every iteration also weights each point by its leverage, the hat value h_i = w_i x_i (X^H W X)^{-1} x_i^H (x_i a row): points above the leverage_probability point of the hat values’ distribution for complex Gaussian inputs, Beta(p, n - p), get (cutoff / h_i)^2. The leverage weights accumulate from iteration to iteration, so a point once downweighted stays so and none keeps more leverage than the cutoff. leverage_cutoff is the older, single-pass form: hat values above that many times their mean p/n downweighted in proportion, after the M-estimate.

4. optionally, noise classes (the noise-weighted estimate of Chen 2020): when the noise grows with the signal - the residual power |r_i|^2 rising with the input power q_i = |x_i|^2, from more than one source - the residuals are heavy-tailed for every M-estimate. The points are split into k classes by k-medians on (q_i, |r_i|^2); in each class a bisquare line gives the noise power a + b q_i, and every row of the regression is divided by its square root, so the noise is the same everywhere; the M-estimate is then made again. k grows from 1 until the normalised residuals’ stable index alpha (McCulloch’s quantile estimate; 2 is Gaussian) reaches noise_alpha; failing that, the k with the largest alpha is kept, and none if no k beats the unweighted residuals.

The error model is EMTF’s: var(b_j) = sigma^2 S_jj with the residual variance sigma^2 and S = (R^H W X)^{-1} (R^H W^2 R) (X^H W R)^{-1}, the inverse signal power ((X^H X)^{-1} for least squares). With jackknife_groups, the variance is also estimated without that model, by the delete-a-group jackknife (Chave & Thomson 1989): the points’ windows, in time order, are cut into that many groups, the estimate is made again without each group at the final weights, and the spread of those estimates is the variance. A group holds whole windows, so the correlation between a window’s harmonics - which the parametric error needs inflated for - is carried by it. Both variances are of the complex coefficient, E|db|^2: the sum of its real and imaginary parts’ variances.

Chave, A. D. & Thomson, D. J. (1989). Some comments on magnetotelluric response function estimation. Journal of Geophysical Research, 94(B10), 14215-14225. https://doi.org/10.1029/JB094iB10p14215

Gamble, T. D., Goubau, W. M. & Clarke, J. (1979). Magnetotellurics with a remote magnetic reference. Geophysics, 44(1), 53-68. https://doi.org/10.1190/1.1440923

Chave, A. D., Thomson, D. J. & Ander, M. E. (1987). On the robust estimation of power spectra, coherences, and transfer functions. Journal of Geophysical Research, 92(B1), 633-648. https://doi.org/10.1029/JB092iB01p00633

Chave, A. D. & Thomson, D. J. (2004). Bounded influence magnetotelluric response function estimation. Geophysical Journal International, 157(3), 988-1006. https://doi.org/10.1111/j.1365-246X.2004.02203.x

Chen, H. (2020). Study on the processing method of magnetotelluric impedance estimation for non-Gaussian noise. Undergraduate thesis, School of Geosciences and Info-Physics, Central South University (advisor R. Guo).

McCulloch, J. H. (1986). Simple consistent estimators of stable distribution parameters. Communications in Statistics - Simulation and Computation, 15(4), 1109-1136. https://doi.org/10.1080/03610918608812563

class PyHydroGeophysX.data_processing.mt.robust.RegressionConfig(huber: 'float' = 1.5, redescend: 'float' = 2.8, max_iterations: 'int' = 10, redescend_iterations: 'int' = 2, tolerance: 'float' = 0.005, leverage_cutoff: 'Optional[float]' = None, bounded_influence: 'bool' = False, leverage_probability: 'float' = 0.999, jackknife_groups: 'int' = 0, noise_classes: 'int' = 0, noise_alpha: 'float' = 1.8)

Bases: object

bounded_influence: bool = False

huber: float = 1.5

jackknife_groups: int = 0

leverage_cutoff: float | None = None

leverage_probability: float = 0.999

max_iterations: int = 10

noise_alpha: float = 1.8

noise_classes: int = 0

none); see the module notes.

Type:

Up to this many noise classes (0

redescend: float = 2.8

redescend_iterations: int = 2

tolerance: float = 0.005

class PyHydroGeophysX.data_processing.mt.robust.RegressionResult(coefficients: 'np.ndarray', weights: 'np.ndarray', residual_variance: 'float', inverse_signal_power: 'np.ndarray', residuals: 'np.ndarray', coherence: 'float', iterations: 'int', jackknife_variance: 'Optional[np.ndarray]' = None, alpha: 'float' = nan, noise_classes: 'int' = 0)

Bases: object

alpha: float = nan

The residuals’ stable index (2 is Gaussian), and the noise classes used.

coefficients: ndarray

coherence: float

inverse_signal_power: ndarray

iterations: int

jackknife_variance: ndarray | None = None

E|b_j - b|^2 from the delete-a-group jackknife, when it was asked for.

noise_classes: int = 0

residual_variance: float

residuals: ndarray

weights: ndarray

PyHydroGeophysX.data_processing.mt.robust.robust_regression(y: ndarray, X: ndarray, R: ndarray | None = None, config: RegressionConfig | None = None, blocks: ndarray | None = None) → RegressionResult

Estimate b in y = X b robustly; R is the remote reference (default X).

blocks labels each point with its window (its start time, say): the jackknife then deletes whole windows, in groups of neighbours in label order, instead of runs of points that may split a window’s harmonics.

PyHydroGeophysX.data_processing.mt.robust.stable_alpha(residuals: ndarray) → float

The stable index of complex residuals (real and imaginary parts pooled), McCulloch’s way.

PyHydroGeophysX.data_processing.mt.spectra module

Fourier coefficients of MT time series, level by level, in field units.

The cascade EMTF made standard (Egbert & Booker 1986; Egbert 1997): level 0 is the recorded series; each further level low-passes and decimates the one before by decimation_factor (a zero-phase FIR, so every channel keeps its timing), and drops the samples at either end that the filter computed partly from beyond the record; the first and last windows of each decimated level are left out too. At every level the series is cut into windows of window_length samples overlapping by overlap; each window is detrended, first differenced (prewhitening, undone after the transform), tapered and transformed, and every Fourier coefficient is divided by its channel’s response at that frequency - so the coefficients are mV/km and nT, whatever the instrument. Windows sit on a grid counted from each level’s first sample, so two runs cut to the same span - a site and its remote reference - share their windows.

A frequency band is a set of harmonics of one level’s windows, as in an EMTF band-setup file (level first last, level counted from 1).

With tapers above one, each window is instead tapered by that many Slepian sequences (Thomson’s multitaper method, time-bandwidth (tapers + 1) / 2): their coefficients at a harmonic are uncorrelated estimates of the same spectrum, so a window gives a band tapers points per harmonic, where one taper would give one - most worth it in the narrow bands at the long-period end of each level. The price is resolution: a coefficient averages the spectrum over ±(tapers + 1) / 2 harmonics.

Egbert, G. D. & Booker, J. R. (1986). Robust estimation of geomagnetic transfer functions. Geophysical Journal of the Royal Astronomical Society, 87(1), 173-194. https://doi.org/10.1111/j.1365-246X.1986.tb04552.x

Thomson, D. J. (1982). Spectrum estimation and harmonic analysis. Proceedings of the IEEE, 70(9), 1055-1096. https://doi.org/10.1109/PROC.1982.12433

class PyHydroGeophysX.data_processing.mt.spectra.Band(level: int, first: int, last: int, skip: Tuple[int, ...] = ())

Bases: object

Harmonics first..last of the windows of decimation level (0-based), less any in skip.

center(sample_rate: float, window_length: int) → float

The band’s frequency: the mean of its harmonics’ (EMTF’s convention).

first: int

frequencies(sample_rate: float, window_length: int) → ndarray

harmonics() → ndarray

last: int

level: int

skip: Tuple[int, ...] = ()

class PyHydroGeophysX.data_processing.mt.spectra.LevelSpectra(level: int, sample_rate: float, window_starts: ndarray, harmonics: ndarray, coefficients: Dict[str, ~numpy.ndarray]=<factory>)

Bases: object

One level’s Fourier coefficients: coefficients[c] is (windows, harmonics).

With several tapers, harmonics holds each harmonic once per taper (taper by taper), and so do the coefficients’ columns.

coefficients: Dict[str, ndarray]

harmonics: ndarray

level: int

sample_rate: float

window_starts: ndarray

class PyHydroGeophysX.data_processing.mt.spectra.SpectraConfig(window_length: 'int' = 128, overlap: 'int' = 32, decimation_factor: 'int' = 4, taper: 'str' = 'hamming', prewhiten: 'bool' = True, min_windows: 'int' = 5, tapers: 'int' = 1)

Bases: object

decimation_factor: int = 4

min_windows: int = 5

overlap: int = 32

prewhiten: bool = True

taper: str = 'hamming'

tapers: int = 1

Slepian tapers per window; above one they replace taper (multitaper).

window_length: int = 128

PyHydroGeophysX.data_processing.mt.spectra.band_data(level_spectra: Sequence[LevelSpectra], band: Band, components: Sequence[str]) → Tuple[Dict[str, ndarray], ndarray]

One band’s coefficients, flattened over windows and harmonics, with window times.

PyHydroGeophysX.data_processing.mt.spectra.band_harmonics(bands: Sequence[Band]) → Dict[int, ndarray]

PyHydroGeophysX.data_processing.mt.spectra.default_bands(n_levels: int, window_length: int = 128) → List[Band]

EMTF’s standard band set for n_levels levels of 128-point windows.

For other window lengths the harmonics are scaled by window_length / 128.

PyHydroGeophysX.data_processing.mt.spectra.log_bands(sample_rate: float, n_levels: int, window_length: int = 128, decimation_factor: int = 4, per_decade: int = 10, width: float = 0.3333333333333333) → List[Band]

Overlapping bands centred on frequencies evenly spaced in log, per_decade to a decade.

After Chen (2020): a band holds the harmonics within 1 +- width of its centre (width = 1/Q, Q = 3), so it spans a factor of about 2 and neighbouring bands share most of their harmonics - the curve is continuous, as in commercial processing that smooths over frequency, at the price of resolution and of errors correlated between neighbours. Each centre is taken at the deepest level that holds it at harmonic 20 (of 128) or below - the longest windows that resolve it - and at level 0 up to harmonic 30; never below harmonic 5. The centres are the same for every sample rate, so two rates meet on the same grid.

PyHydroGeophysX.data_processing.mt.spectra.n_levels_for(n_samples: int, config: SpectraConfig) → int

How many levels leave at least min_windows windows that are used.

Counted as run_spectra() makes them: each decimation keeps (n - 1) // D + 1 samples less the _EDGE trimmed at each end, and a decimated level leaves out its first and last window.

PyHydroGeophysX.data_processing.mt.spectra.read_band_setup(path: Any) → List[Band]

An EMTF band-setup file: the number of bands, then level first last rows.

PyHydroGeophysX.data_processing.mt.spectra.run_spectra(run: TimeSeriesRun, components: Sequence[str], harmonics: Dict[int, ndarray], config: SpectraConfig) → List[LevelSpectra]

Calibrated Fourier coefficients of components at the levels in harmonics.

harmonics maps a level to the harmonic indices its bands use. Windows holding a non-finite sample are dropped.

PyHydroGeophysX.data_processing.mt.spectra.shared_windows(local: LevelSpectra, remote: LevelSpectra) → Tuple[LevelSpectra, LevelSpectra]

The two levels cut to the windows they share, matched by start time.

A remote window is the local one’s when it starts within half a sample of it at that level: start times rounded to the nanosecond, or a remote clock a fraction of a sample off, still pair (a constant offset of the remote changes the reference’s phase only, which cancels).

PyHydroGeophysX.data_processing.mt.spectra.variance_inflation(config: SpectraConfig, n_harmonics: int) → float

How much a band of n_harmonics adjacent harmonics overstates its information.

The taper makes the coefficients of neighbouring harmonics - and of overlapping windows - correlated, for the noise and the fields alike, so a regression on them has (1/m) sum_kl |rho_(k-l)|^2 times the variance it would have on independent points, rho_k the taper’s correlation at a lag of k harmonics; overlapping windows add the same sum over the pairs of points one window and the next share.

With several tapers the sums run over every pair of points, rho then the correlation of taper j at one harmonic with taper l at another: zero between two tapers at the same harmonic, which is why they add information. (A Monte Carlo of white noise agrees within a few per cent, for one taper and for several.)

PyHydroGeophysX.data_processing.mt.timeseries module

MT time series as the instruments record them, with their responses.

Every instrument reader returns the same objects:

• TimeSeriesRun - one stretch of simultaneous, gap-free samples of several channels at one sample rate: a continuous recording, or one burst of a segmented (scheduled) one. start is UTC as numpy.datetime64.

• Channel - one component’s samples in the units the logger stored them in (volts, millivolts, counts, or nT for a fluxgate), with what is needed to reach the field: its ResponseStage chain, orientation, dipole length and sensor.

Response convention. Each stage maps its input to its output, physical side first: a coil is nT -> mV with a table of mV/nT, a dipole mV/km -> mV with gain L / 1000, a receiver’s anti-alias filter V -> V. The product of the stages is what the recorded spectrum is per unit of field, so the field’s spectrum is the recorded one divided by Channel.response(). Electric fields come out in mV/km and magnetic ones in nT, the units every MT code and file uses.

Precision. Data stay in the precision of the source: 24-bit samples and float32 files as float32 (exact for 24-bit counts), 32-bit counts and text as float64. A fluxgate’s total field (~50 000 nT at pT resolution) needs float64.

PyHydroGeophysX.data_processing.mt.timeseries.CLIP_LEVEL = 0.995

a decimated stream’s filter rounds the flat top of a clipped stretch off a little below it.

Type:

Share of a channel’s full scale read as clipped

PyHydroGeophysX.data_processing.mt.timeseries.COMPONENT_ALIASES = {'bx': 'hx', 'by': 'hy', 'bz': 'hz', 'e1': 'ex', 'e2': 'ey', 'h1': 'hx', 'h2': 'hy', 'h3': 'hz'}

E1 is Ex, H3 is Hz…

Type:

What a logger’s channel name means, lower-cased

class PyHydroGeophysX.data_processing.mt.timeseries.Channel(component: str, data: ndarray, units: str, response: List[ResponseStage] = <factory>, azimuth: float = nan, tilt: float = 0.0, dipole_length: float = nan, sensor: str = '', sensor_serial: str = '', metadata: Dict[str, ~typing.Any]=<factory>)

Bases: object

One component’s samples, in units, with the chain that reaches the field.

metadata["full_scale"], where the reader knows it, is the level in units at which the recorder clips (Phoenix: the A/D’s 5 V over the gain).

azimuth: float = nan

component: str

copy(data: ndarray | None = None) → Channel

data: ndarray

dipole_length: float = nan

property is_calibrated: bool

Whether removing the response gives mV/km (E) or nT (H).

property kind: str

metadata: Dict[str, Any]

property physical_units: str

The units of the field once the response is removed.

response: List[ResponseStage]

response_at(frequency: Any) → ndarray

Recorded units per physical unit at frequency: the stages’ product.

sensor: str = ''

sensor_serial: str = ''

tilt: float = 0.0

units: str

class PyHydroGeophysX.data_processing.mt.timeseries.ResponseStage(name: str, input_units: str, output_units: str, gain: float = 1.0, frequency: ndarray | None = None, values: ndarray | None = None, metadata: Dict[str, ~typing.Any]=<factory>)

Bases: object

One stage of an instrument’s response: a gain, a table, or both.

frequency (Hz) and complex values tabulate the stage; between frequencies the amplitude is interpolated in log-log and the unwrapped phase in log frequency. Outside the table the amplitude follows the slope of its last three points in log-log (a coil’s response falls as f below its corner) and the phase is held; frequency_range is what was measured. gain multiplies the table (or stands alone).

frequency: ndarray | None = None

property frequency_range: Tuple[float, float]

The tabulated band in Hz; a scalar stage holds at every frequency.

classmethod from_amplitude_phase(name: str, input_units: str, output_units: str, frequency: Any, amplitude: Any, phase: Any, *, phase_units: str = 'deg', gain: float = 1.0, metadata: Dict[str, Any] | None = None) → ResponseStage

gain: float = 1.0

input_units: str

property is_table: bool

metadata: Dict[str, Any]

name: str

output_units: str

values: ndarray | None = None

class PyHydroGeophysX.data_processing.mt.timeseries.TimeSeriesRun(channels: ~typing.Dict[str, ~PyHydroGeophysX.data_processing.mt.timeseries.Channel], sample_rate: float, start: ~typing.Any, station: str = '', latitude: float = nan, longitude: float = nan, elevation: float = nan, declination: float = 0.0, instrument: str = '', serial: str = '', metadata: ~typing.Dict[str, ~typing.Any] = <factory>)

Bases: object

Simultaneous, gap-free samples of several channels at one rate.

array(components: Sequence[str] | None = None, dtype: Any = <class 'float'>) → ndarray

The chosen channels stacked as (n_channels, n_samples).

calibrated(components: Sequence[str] | None = None, *, band: Tuple[float, float] | None = None) → TimeSeriesRun

The channels in field units (mV/km, nT), the response removed by FFT.

Each channel is demeaned and divided by its response at every Fourier frequency; band (Hz) zeroes the spectrum outside it, which a coil’s vanishing low-frequency response needs. This holds the whole series in memory as complex numbers: for long, fast recordings select a window first.

channels: Dict[str, Channel]

property components: List[str]

declination: float = 0.0

property duration: float

Seconds the samples span, n / sample_rate.

elevation: float = nan

property end: datetime64

The time just after the last sample.

instrument: str = ''

latitude: float = nan

longitude: float = nan

metadata: Dict[str, Any]

property n_samples: int

sample_rate: float

select(components: Iterable[str]) → TimeSeriesRun

serial: str = ''

slice(first: int, last: int | None = None) → TimeSeriesRun

Samples first:last (views, not copies).

start: Any

station: str = ''

summary() → str

times() → ndarray

Each sample’s UTC time as datetime64[ns].

window(start: Any = None, end: Any = None) → TimeSeriesRun

The samples between two UTC times (either may be None).

PyHydroGeophysX.data_processing.mt.timeseries.as_datetime64(value: Any) → datetime64 | None

value (string, datetime, datetime64) as UTC datetime64[ns]; None stays None.

PyHydroGeophysX.data_processing.mt.timeseries.assemble_runs(pieces: Dict[str, List[Tuple[datetime64, ndarray]]], sample_rate: float, make_channel, *, min_samples: int = 1, **run_fields: Any) → List[TimeSeriesRun]

Runs from each component’s contiguous segments, where all of them overlap.

pieces maps a component to its (start, data) segments; make_channel(component, data) builds the Channel. Overlaps are cut to whole samples on a common grid; offsets of a fraction of a sample between channels are reported in the run’s notes.

PyHydroGeophysX.data_processing.mt.timeseries.contiguous_blocks(blocks: Sequence[Tuple[datetime64, ndarray]], sample_rate: float, *, tolerance: float = 0.25) → List[Tuple[datetime64, ndarray]]

Join time-stamped blocks that follow each other without a gap.

A block continues the previous one when it starts within tolerance samples of where that one ended; otherwise a new segment starts.

PyHydroGeophysX.data_processing.mt.timeseries.dipole_stage(length_m: float) → ResponseStage

The dipole: mV/km -> mV, gain L / 1000.

PyHydroGeophysX.data_processing.mt.timeseries.gps_leap_seconds(when: Any) → int

GPS - UTC at a date, from the leap-second table (13 s before 1999).

PyHydroGeophysX.data_processing.mt.timeseries.gps_seconds_to_utc(seconds: Any, *, epoch: str = 'unix', leap_seconds: int | None = None) → datetime64

UTC of a GPS-scale time stamp.

epoch="unix" reads seconds as GPS time counted from 1970-01-01, as Phoenix writes it; "gps" counts from 1980-01-06. leap_seconds is GPS - UTC; by default it comes from the table for that date.

PyHydroGeophysX.data_processing.mt.timeseries.load_runs(source: Any) → List[TimeSeriesRun]

The runs a save_runs() container holds.

PyHydroGeophysX.data_processing.mt.timeseries.normalize_component(name: str) → str

PyHydroGeophysX.data_processing.mt.timeseries.periodic_interference(run: TimeSeriesRun, periods: Sequence[float], *, threshold: float = 10.0) → Dict[float, List[str]]

Which channels carry something repeating exactly every one of periods seconds.

Each channel is band-passed to 0.1-20 Hz, so neither the field’s long periods nor power-line hum swamps a pulse that lasts a fraction of a second, and folded on the period; a mean waveform standing more than threshold times above what chance leaves (about 4) is interference. The periods are tried shortest first, each one found removed before the next is tried: folded on 30 s, a pulse every 4 s still stands out (at half its height, since 30 s is not a multiple of 4), and would pass for a second interference. Returns {period: [channel, ...]} for the periods found.

PyHydroGeophysX.data_processing.mt.timeseries.powerline_frequency(run: TimeSeriesRun, candidates: Sequence[float] = (50.0, 60.0), *, threshold: float = 100.0) → float | None

The grid frequency whose hum run carries, of candidates: the one whose line stands highest above its floor, if more than threshold times; None when none does or none lies below the Nyquist.

PyHydroGeophysX.data_processing.mt.timeseries.read_response_table(path: Any, *, name: str | None = None, input_units: str = 'nT', output_units: str = 'mV', phase_units: str = 'deg', amplitude_per_hz: bool = False) → ResponseStage

A calibration table from a text file of frequency, amplitude, phase rows.

Lines that do not start with three numbers are skipped, so headers and comments may be anywhere. amplitude_per_hz multiplies the amplitude by frequency, for tables normalized that way (Metronix V/(nT Hz)).

PyHydroGeophysX.data_processing.mt.timeseries.remove_periodic(run: TimeSeriesRun, period: float, *, cycles: int = 31) → TimeSeriesRun

run less whatever repeats exactly every period seconds - a receiver’s own interference.

Each channel is first high-passed by subtracting its centred moving mean over one period, which passes the harmonics k / period whole and takes out everything slower; the waveform of each cycle is then the mean of that over the cycles cycles around it, and is subtracted from the channel. The whole is one time-invariant, zero-phase comb filter, the same on every channel, so an impedance is left as it was: notches about 1 / (cycles * period) wide at the harmonics, the frequencies between them kept. (Taking each cycle’s own mean off instead turned the field’s slow drift into a staircase with steps at every cycle - a near-real “impedance” carrying the site’s static shift.) Interference timed by the receiver’s GPS clock repeats exactly; the field does not. period must be a whole number of samples; a run shorter than three periods is returned as it is.

PyHydroGeophysX.data_processing.mt.timeseries.remove_periodic_bursts(runs: Sequence[TimeSeriesRun], period: float, *, cycles: int = 31) → List[TimeSeriesRun]

Bursts of a segmented recording less what repeats every period seconds.

A burst shorter than a few periods cannot be folded on itself, but bursts that start at the same phase of the period - and have the same rate and length - carry the interference at the same place. Each such burst has the mean of the other bursts of that group among the cycles around it subtracted (its own field is left out of its template), less its own mean; the field, which does not repeat, averages out of it. A group of fewer than three bursts is left as recorded. This also takes off whatever every burst of the group starts with alike.

PyHydroGeophysX.data_processing.mt.timeseries.remove_powerline(run: TimeSeriesRun, frequency: float = 60.0, *, threshold: float = 10.0) → TimeSeriesRun

run with the power-line hum notched out: frequency and its harmonics below the Nyquist.

Hum is many orders of magnitude stronger than the field, and in a record whose Nyquist lies above it - 60 Hz in 150 Hz data - it leaks through the taper’s sidelobes into bands well away from it. There it passes for an impedance with apparent resistivity growing as the period and a phase of 0 or 180 degrees, and as the grid is the same at a remote site, remote referencing keeps it. Each line is found where the channels’ power, smoothed over 0.01 Hz, stands threshold times above the floor 1-3 Hz either side, and every bin of it - as wide as the grid’s frequency wandered over the run - is set to zero, in one transform of the whole run: a zero-phase notch that leaves the frequencies between the harmonics as they were. A harmonic that does not stand out is left alone.

A run shorter than a minute - a burst of a segmented recording - has bins too wide for that: hum that is not on a bin leaks across many of them and stays. There each line is fitted instead, as a sinusoid at its own frequency (the peak of the four-times zero-padded spectrum) with an amplitude that may change linearly over the run, and subtracted.

PyHydroGeophysX.data_processing.mt.timeseries.runs_from_payload(payload: Sequence[Dict[str, Any]]) → List[TimeSeriesRun]

The runs runs_to_payload() took apart.

PyHydroGeophysX.data_processing.mt.timeseries.runs_to_payload(runs: Sequence[TimeSeriesRun]) → List[Dict[str, Any]]

Runs as plain values and arrays, for PyHydroGeophysX.data_processing.run_inputs.save_container().

PyHydroGeophysX.data_processing.mt.timeseries.saturated(run: TimeSeriesRun, components: Sequence[str] | None = None, *, level: float = 0.995, guard: int = 2) → ndarray | None

Which samples of run have one of components at the recorder’s full scale.

Only channels whose reader knows the full scale (metadata["full_scale"]) count; None when none does. Each clipped sample is widened by guard samples either side. A clipped stretch is the recorder’s limit, not the field: windows holding it give an impedance that is plausible and wrong.

PyHydroGeophysX.data_processing.mt.timeseries.save_runs(destination: Any, runs: Sequence[TimeSeriesRun])

Write runs as one compressed container (.npz); returns its path.

PyHydroGeophysX.data_processing.mt.timeseries.seconds_between(start: datetime64, end: datetime64) → float

PyHydroGeophysX.data_processing.mt.timeseries.shifted(start: datetime64, seconds: float) → datetime64

PyHydroGeophysX.data_processing.mt.timeseries.unit_stage(input_units: str, output_units: str) → ResponseStage

PyHydroGeophysX.data_processing.mt.timeseries.without_samples(run: TimeSeriesRun, mask: ndarray) → TimeSeriesRun

run with the samples mask marks set to NaN in every channel (their windows are then dropped).

PyHydroGeophysX.data_processing.mt.transfer_function module

The magnetotelluric transfer functions of one site: impedance and tipper.

One container serves every reader, writer and processing step, so a site read from an EDI file, an EMTF XML file or a Z-file, or estimated from time series, is the same object. Its conventions are fixed, and every reader converts to them:

• Impedance in ohms (SI), E [V/m] / H [A/m]. Files usually carry field units, (mV/km)/nT, which are 1 / (4 pi 1e-4) times larger.

• Time dependence e^{+i omega t}, the SEG EDI standard’s: a half-space gives a phase of +45 degrees for Zxy and -135 degrees for Zyx.

• Rotation is the azimuth of the x axis, in degrees clockwise from north, for every frequency. Rotating Z' = R Z R^T and T' = T R^T.

• Errors are one standard error of the real part and of the imaginary part alike - the square root of the EDI Z**.VAR, Gamble’s average of the two variances. When the full covariance is known (EMTF XML, Z-files and the robust estimator here: inverse signal power and residual covariance), it is kept, and errors are derived from it, so they rotate correctly.

The apparent resistivity is |Z|^2 / (omega mu0), and the phase atan2(Im Z, Re Z) in degrees; their errors follow from the impedance error to first order, drho / rho = 2 dZ / |Z| and dphi = dZ / |Z| radians.

PyHydroGeophysX.data_processing.mt.transfer_function.FIELD_TO_OHM = 0.0012566370614359172

(mV/km)/nT -> ohm. E 1e-6 V/m over H = B / mu0 with B 1e-9 T.

class PyHydroGeophysX.data_processing.mt.transfer_function.TransferFunction(frequency: ndarray, z: ndarray | None = None, z_err: ndarray | None = None, tipper: ndarray | None = None, tipper_err: ndarray | None = None, rotation: Any = 0.0, inverse_signal_power: ndarray | None = None, z_residual_covariance: ndarray | None = None, tipper_residual_covariance: ndarray | None = None, station: str = '', latitude: float = nan, longitude: float = nan, elevation: float = nan, x: float = nan, y: float = nan, declination: float = nan, metadata: Dict[str, ~typing.Any]=<factory>)

Bases: object

Impedance and tipper of one MT site, against frequency.

z is (n, 2, 2) complex ohms, tipper (n, 1, 2) complex and dimensionless; either may be None. z_err and tipper_err are the standard errors of each real and imaginary part. inverse_signal_power (n, 2, 2) and the residual covariances - z_residual_covariance (n, 2, 2) for (Ex, Ey) and tipper_residual_covariance (n, 1, 1) for Hz - are the full error model where a file or the estimator supplies it: var(Z_ij) = residual_ii S_jj. Values that are missing are NaN.

property angular_frequency: ndarray

apparent_resistivity() → ndarray

|Z|^2 / (omega mu0) in ohm-m, (n, 2, 2).

apparent_resistivity_err() → ndarray

First-order standard error of the apparent resistivity, ohm-m.

conjugated() → TransferFunction

The site under the opposite time convention (e^{-i omega t} <-> e^{+i omega t}).

copy() → TransferFunction

declination: float = nan

determinant() → ndarray

The rotation-invariant determinant impedance sqrt(Zxx Zyy - Zxy Zyx), (n,).

The branch is the one with the phase of Zxy, so a half-space gives +45.

elevation: float = nan

frequency: ndarray

property has_impedance: bool

property has_tipper: bool

inverse_signal_power: ndarray | None = None

latitude: float = nan

longitude: float = nan

metadata: Dict[str, Any]

property n_frequencies: int

property period: ndarray

phase() → ndarray

Impedance phase in degrees, (n, 2, 2); Zyx sits near -135 for a half-space.

phase_err() → ndarray

First-order standard error of the phase, degrees.

rotated(angle: Any) → TransferFunction

The site with its x axis at azimuth angle (degrees from north, clockwise).

angle is one value or one per frequency. A missing element (Zxx and Zyy of a site given as Zxy and Zyx only) leaves missing only what it enters: a rotation by 0 or 90 degrees keeps the rest.

The errors are recomputed from the rotated covariance when the site’s errors are the covariance’s (EMTF XML, Z-files); otherwise - errors of their own, such as the jackknife’s - each rotated variance is sum (R_ik R_jl)^2 var_kl from the errors in the frame the site was read in. That assumes the elements independent, an approximation that does not depend on the path the rotations took and leaves the errors as they were at the frame they came in.

rotation: Any = 0.0

select(mask: Any) → TransferFunction

The frequencies a boolean mask or an index array selects.

sorted_by_period() → TransferFunction

Frequencies in order of increasing period.

station: str = ''

summary() → Dict[str, Any]

A JSON-safe description of the site, for logs and workflow results.

tipper: ndarray | None = None

tipper_err: ndarray | None = None

tipper_residual_covariance: ndarray | None = None

x: float = nan

Local coordinates (m), when the site belongs to a projected survey.

y: float = nan

z: ndarray | None = None

z_err: ndarray | None = None

property z_field_units: ndarray | None

The impedance in (mV/km)/nT.

z_residual_covariance: ndarray | None = None

PyHydroGeophysX.data_processing.mt.transfer_function.Z_COMPONENTS = ('xx', 'xy', 'yx', 'yy')

Which impedance element holds which output/input pair.

PyHydroGeophysX.data_processing.mt.transfer_function.from_apparent_resistivity(frequency: Sequence[float], rho: Dict[str, Sequence[float]], phase_deg: Dict[str, Sequence[float]], *, rho_err: Dict[str, Sequence[float]] | None = None, phase_err_deg: Dict[str, Sequence[float]] | None = None, **site: Any) → TransferFunction

A site from apparent resistivities and phases, keyed "xy", "yx", …

|Z| = sqrt(rho omega mu0), Z = |Z| e^{i phi}. The impedance error is taken from the resistivity error, dZ = |Z| drho / (2 rho), or else from the phase error, dZ = |Z| dphi. Elements not given are NaN.

PyHydroGeophysX.data_processing.mt.transfer_function.resolve_sign_convention(tf: TransferFunction, convention: str) → TransferFunction

tf in e^{+i omega t}, given the convention it was read in.

"+" keeps it, "-" conjugates it, and "auto" conjugates it only when, at four frequencies in five (of three or more), the phases of Zxy and Zyx lie closer to (-45, +135) - a half-space under e^{-i omega t} - than to (+45, -135). Distances rather than quadrants, because the phases of real data stray past 90 degrees at long periods. The choice is recorded in metadata["sign_convention"] and, when it changed the data, in the notes.

PyHydroGeophysX.data_processing.mt.transfer_function.rotation_matrix(angle_deg: Any) → ndarray

R for a frame rotated clockwise by angle_deg: shape (..., 2, 2).

A vector’s components in the rotated frame are R v: the new x axis, at azimuth theta, is (cos theta, sin theta) in the old (north, east) frame.

PyHydroGeophysX.data_processing.mt.zfile module

Egbert’s EMTF Z-files (.zmm, .zrr, .zss): read.

The text output of Egbert’s EMTF processing: a header with the site coordinates, the declination and each channel’s azimuth and tilt, then, per period, the transfer functions of every output channel (Hz, Ex, Ey, …) on the inputs Hx, Hy, the inverse coherent signal power matrix S (2 x 2) and the residual covariance (outputs x outputs), both stored as lower triangles of complex numbers. The variance of an element is residual_ii S_jj; this reader keeps S and the covariance, so errors rotate correctly.

The transfer functions are in measurement coordinates, field units ((mV/km)/nT) and e^{+i omega t}. Channels laid out orthogonally keep that frame, with the rotation the Hx azimuth; otherwise the tensor is transformed to the orthogonal frame of the azimuths’ north, Z = U_e^-1 Z_m U_h.

Egbert, G. D. (1997). Robust multiple-station magnetotelluric data processing. Geophysical Journal International, 130(2), 475-496. https://doi.org/10.1111/j.1365-246X.1997.tb05663.x

PyHydroGeophysX.data_processing.mt.zfile.is_zfile(path: Any) → bool

PyHydroGeophysX.data_processing.mt.zfile.read_zfile(path: Any) → TransferFunction

Read an EMTF Z-file into a TransferFunction.

PyHydroGeophysX.data_processing.mt.zonge module

Zonge International ZEN receivers: Z3D files, one per channel.

A Z3D file starts with 512-byte text records: the GPS board’s header (A/D Rate, Ch.Factor, GpsWeek, Lat/Long in radians…), the schedule, and metadata records of |KEY=VALUE| pairs - CH.CMP, CH.OFFSET.XYZ1/XYZ2 (electrode offsets, m), CH.ANTSN (coil serial), RX.XAZIMUTH and, in magnetic files, the coil calibration CAL.ANT (frequency in Hz : amplitude in mV/nT : phase in mrad - read as radians per second instead, the frequencies give impedance phases outside their quadrants on Zonge’s sample site, as Hz they give a smooth sounding whose 4096 and 256 S/s runs agree where they overlap). Then come int32 samples, one second per block, each block led by a 64-byte GPS stamp that starts with the words 0x7FFFFFFF, 0x80000000 and gives the time in 1/1024 s of the GPS week. Samples times Ch.Factor are volts at the input (divided by the channel gain).

Files of one station, sample rate and schedule start form a run; the first two seconds are dropped, because the ZEN’s SD-card buffer can corrupt them, and a block that does not hold one second of samples ends a run.

Zonge International (2014). ZEN receiver and Z3D file format notes, in ZenACQ documentation. Zonge International, Tucson, AZ.

PyHydroGeophysX.data_processing.mt.zonge.is_z3d_file(path: Any) → bool

PyHydroGeophysX.data_processing.mt.zonge.read_z3d_header(path: Any) → Dict[str, Any]

The text records of a Z3D file: header, schedule and metadata, and the data offset.

PyHydroGeophysX.data_processing.mt.zonge.read_zonge(path: Any, *, components: Iterable[str] | None = None, calibration: Dict[str, Any] | None = None, start: Any = None, end: Any = None, skip_seconds: int = 2) → List[TimeSeriesRun]

Read Z3D files - one file or a folder - into runs.

calibration maps a component (or a coil serial) to a ResponseStage or a text table (frequency Hz, amplitude mV/nT, phase in degrees) that replaces the coil table stored in the file.

Module contents

Magnetotelluric (MT) data: transfer functions, instrument time series and processing.

Everything here is written for PyHydroGeophysX on NumPy and SciPy alone.

• transfer_function - TransferFunction, one site’s impedance and tipper with their errors, in SI units and e^{+i omega t};

• edi, emtf_xml, zfile, jfile - the exchange formats; read_transfer_function() picks the reader;

• timeseries - TimeSeriesRun and Channel, recorded samples with their instrument responses;

• phoenix, phoenix_legacy, metronix, zonge, lemi - the instruments’ own files (Phoenix MTU-5C/5P/8A and MTU-5A/V8, Metronix ADU, Zonge ZEN, LEMI-424); read_timeseries() picks the reader;

• processing (with spectra and robust) - process_mt(), robust single-site and remote-reference estimation of the impedance and tipper from time series;

• analysis - phase tensor, skews, strike, induction arrows, Niblett-Bostick depths and static shift;

• forward1d - the layered-earth impedance and its derivatives;

• inversion1d - Occam’s 1D inversion, with static shift and jointly with a TEM sounding, an azimuthally anisotropic variant, and the water content of the result;

• inversion2d - 2D forward modelling and inversion of a profile on SimPEG (an optional dependency, imported when used).

class PyHydroGeophysX.data_processing.mt.Anisotropic1DResult(strike: float, resistivity: ~numpy.ndarray, thickness: ~numpy.ndarray, rms: float, mu: float, static_shift: ~typing.Dict[str, float], predicted: ~typing.Dict[str, ~typing.Dict[str, ~numpy.ndarray]], sounding: ~PyHydroGeophysX.data_processing.mt.inversion1d.Sounding1D, iterations: int, history: ~typing.List[~typing.Dict[str, float]] = <factory>, converged: bool = False)

Bases: object

Occam’s azimuthally anisotropic layered model: two resistivities per layer on one strike.

resistivity[0] is each layer’s resistivity along the axis at azimuth strike (degrees clockwise from north), which Zxy fits in that frame; resistivity[1] along the axis a quarter turn on, which Zyx fits.

property anisotropy: ndarray

Each layer’s resistivity across the strike over along it.

axis(k: int) → Occam1DResult

One axis’s model as an Occam1DResult, for the plots and conversions that take one.

converged: bool = False

property depth: ndarray

Each layer’s top, m (the half-space last).

history: List[Dict[str, float]]

iterations: int

mu: float

predicted: Dict[str, Dict[str, ndarray]]

resistivity: ndarray

rms: float

sounding: Sounding1D

static_shift: Dict[str, float]

strike: float

thickness: ndarray

class PyHydroGeophysX.data_processing.mt.Band(level: int, first: int, last: int, skip: Tuple[int, ...] = ())

Bases: object

Harmonics first..last of the windows of decimation level (0-based), less any in skip.

center(sample_rate: float, window_length: int) → float

The band’s frequency: the mean of its harmonics’ (EMTF’s convention).

first: int

frequencies(sample_rate: float, window_length: int) → ndarray

harmonics() → ndarray

last: int

level: int

skip: Tuple[int, ...] = ()

class PyHydroGeophysX.data_processing.mt.Channel(component: str, data: ndarray, units: str, response: List[ResponseStage] = <factory>, azimuth: float = nan, tilt: float = 0.0, dipole_length: float = nan, sensor: str = '', sensor_serial: str = '', metadata: Dict[str, ~typing.Any]=<factory>)

Bases: object

One component’s samples, in units, with the chain that reaches the field.

metadata["full_scale"], where the reader knows it, is the level in units at which the recorder clips (Phoenix: the A/D’s 5 V over the gain).

azimuth: float = nan

component: str

copy(data: ndarray | None = None) → Channel

data: ndarray

dipole_length: float = nan

property is_calibrated: bool

Whether removing the response gives mV/km (E) or nT (H).

property kind: str

metadata: Dict[str, Any]

property physical_units: str

The units of the field once the response is removed.

response: List[ResponseStage]

response_at(frequency: Any) → ndarray

Recorded units per physical unit at frequency: the stages’ product.

sensor: str = ''

sensor_serial: str = ''

tilt: float = 0.0

units: str

class PyHydroGeophysX.data_processing.mt.Occam1DResult(resistivity: 'np.ndarray', thickness: 'np.ndarray', rms: 'float', roughness: 'float', mu: 'float', static_shift: 'Dict[str, float]', predicted: 'Dict[str, np.ndarray]', sounding: 'Sounding1D', iterations: 'int', history: 'List[Dict[str, float]]'=<factory>, tem: 'Optional[Dict[str, Any]]'=None, converged: 'bool' = False)

Bases: object

converged: bool = False

property depth: ndarray

Each layer’s top, m (the half-space last).

depth_profile() → Tuple[ndarray, ndarray]

Step-plot arrays (depth, resistivity), the half-space drawn one layer deep.

history: List[Dict[str, float]]

iterations: int

mu: float

predicted: Dict[str, ndarray]

resistivity: ndarray

rms: float

roughness: float

sounding: Sounding1D

static_shift: Dict[str, float]

tem: Dict[str, Any] | None = None

thickness: ndarray

class PyHydroGeophysX.data_processing.mt.ProcessingConfig(window_length: int = 128, overlap: int = 32, decimation_factor: int = 4, n_levels: int | None = None, min_windows: int = 5, bands: List[Band] | None = None, bands_per_decade: int | None = None, taper: str = 'hamming', prewhiten: bool = True, huber: float = 1.5, redescend: float = 2.8, max_iterations: int = 10, redescend_iterations: int = 2, tolerance: float = 0.005, leverage_cutoff: float | None = None, min_points: int = 10, allow_uncalibrated: bool = False, tipper: bool = True, periodic_noise: Tuple[float, ...] = (), periodic_cycles: int = 31, periodic_where_found: bool = True, skip_notched: bool = True, widen_above: float = 0.0, powerline: float | None = None, drop_saturated: bool = True, bounded_influence: bool = False, jackknife: bool = False, min_coherence: float = 0.0, presort_windows: int = 10, min_reference_coherence: float = 0.0, tapers: int = 1, multitaper_above: float = 0.0, multitaper_tapers: int = 3, multitaper_periods: Tuple[float, float] | None = None, multi_window: bool = False, noise_classes: int = 0)

Bases: object

Every choice of process_mt(); the defaults are EMTF’s.

allow_uncalibrated: bool = False

bands: List[Band] | None = None

bands_per_decade: int | None = None

Overlapping bands, this many to a decade (spectra.log_bands()), in place of EMTF’s band set: continuous curves, errors correlated between neighbours.

bounded_influence: bool = False

Weight points by their leverage too, every iteration (robust): noise in the magnetic channels then cannot pull the estimate.

decimation_factor: int = 4

drop_saturated: bool = True

Leave out the windows holding a sample at the recorder’s full scale (timeseries.saturated()), where the reader knows it.

huber: float = 1.5

jackknife: bool = False

Errors from the delete-a-group jackknife (20 groups) instead of the parametric model (robust).

leverage_cutoff: float | None = None

max_iterations: int = 10

min_coherence: float = 0.0

in each band, stretches of presort_windows consecutive windows whose multiple coherence (squared) of the output with Hx, Hy is below this are left out of that output’s regression.

Type:

Coherence presorting

min_points: int = 10

min_reference_coherence: float = 0.0

where a band’s Hx or Hy is coherent with the remote’s same channel below this (squared coherence), the reference holds almost no natural signal; the elements estimated from that input (Hx: Zxx, Zyx; Hy: Zxy, Zyy; and the tipper’s) are then a ratio of cross-spectra near zero, and the jackknife understates their error. Their relative errors are raised to at least 100 %. 0 is off.

Type:

Remote reference

min_windows: int = 5

multi_window: bool = False

each band also takes the coefficients at its frequencies of the windows one level up or down - windows decimation_factor times longer or shorter - in one regression. Those overlap the band’s own data in time, so the errors are the jackknife’s.

Type:

Multi-window stacking (Ernst et al. 2001)

multitaper_above: float = 0.0

An output (Ex for Zxy, Ey for Zyx, Hz) whose relative error in a band is above this with one taper is estimated again there with multitaper_tapers Slepian tapers, and that estimate kept. Where the field is weakest the tapers’ extra points cut the scatter; where the errors are a few per cent, the tapers’ smoothing over neighbouring harmonics (a bias of a few per cent where the curve is steep, which the errors do not show) would be the larger error, so those bands keep one taper. Where the two estimates differ, which is right is not known: the difference is added to the errors (in quadrature), so the curve is the steadier one and its errors hold the other. 0 is off.

multitaper_periods: Tuple[float, float] | None = None

The periods (s, lowest and highest) whose bands multitaper_above may estimate again; None is every band. Tapers help where few windows are the limit, as in a weak dead band; in short bursts at high rates they can make the curve rougher.

multitaper_tapers: int = 3

n_levels: int | None = None

noise_classes: int = 0

up to this many classes of noise that grows with the signal, each normalised before the M-estimate; 0 is off.

Type:

Noise classes (robust)

overlap: int = 32

periodic_cycles: int = 31

fewer follow interference whose shape drifts, each notch then a little wider.

Type:

Cycles each interference waveform is averaged over

periodic_noise: Tuple[float, ...] = ()

Periods (s) of interference repeating exactly on the receiver’s clock, subtracted from every run before its spectra (timeseries.remove_periodic()).

periodic_where_found: bool = True

Remove a period from a run only where timeseries.periodic_interference() finds it in that run (when the run holds 30 cycles or more, enough to tell). The removal also takes the field at every harmonic of the period, and where a level’s windows cannot resolve those notches it biases the band: so one site’s interference is not taken out of another site’s data. Bursts too short to fold are cleaned for every period, as before.

powerline: float | None = None

Grid frequency (50 or 60 Hz) whose hum and harmonics are notched out of every run before its spectra (timeseries.remove_powerline()); None is off.

presort_windows: int = 10

prewhiten: bool = True

redescend: float = 2.8

redescend_iterations: int = 2

regression() → RegressionConfig

skip_notched: bool = True

Leave out of the bands the harmonics inside those notches, where windows are long enough to resolve them (twice the period or more): the notch takes the field there too, and what is left is the interference’s drift.

spectra() → SpectraConfig

taper: str = 'hamming'

tapers: int = 1

Slepian tapers per window (multitaper, spectra); 1 is taper alone.

tipper: bool = True

tolerance: float = 0.005

widen_above: float = 0.0

A band whose Zxy and Zyx relative errors average above this is estimated again from a band about twice as wide (Q 1.5 instead of 3, same level and centre) and the wider one kept if its errors are smaller: more data where the field is weakest, at the price of resolution. 0 is off.

window_length: int = 128

class PyHydroGeophysX.data_processing.mt.RegressionConfig(huber: 'float' = 1.5, redescend: 'float' = 2.8, max_iterations: 'int' = 10, redescend_iterations: 'int' = 2, tolerance: 'float' = 0.005, leverage_cutoff: 'Optional[float]' = None, bounded_influence: 'bool' = False, leverage_probability: 'float' = 0.999, jackknife_groups: 'int' = 0, noise_classes: 'int' = 0, noise_alpha: 'float' = 1.8)

Bases: object

bounded_influence: bool = False

huber: float = 1.5

jackknife_groups: int = 0

leverage_cutoff: float | None = None

leverage_probability: float = 0.999

max_iterations: int = 10

noise_alpha: float = 1.8

noise_classes: int = 0

none); see the module notes.

Type:

Up to this many noise classes (0

redescend: float = 2.8

redescend_iterations: int = 2

tolerance: float = 0.005

class PyHydroGeophysX.data_processing.mt.ResponseStage(name: str, input_units: str, output_units: str, gain: float = 1.0, frequency: ndarray | None = None, values: ndarray | None = None, metadata: Dict[str, ~typing.Any]=<factory>)

Bases: object

One stage of an instrument’s response: a gain, a table, or both.

frequency (Hz) and complex values tabulate the stage; between frequencies the amplitude is interpolated in log-log and the unwrapped phase in log frequency. Outside the table the amplitude follows the slope of its last three points in log-log (a coil’s response falls as f below its corner) and the phase is held; frequency_range is what was measured. gain multiplies the table (or stands alone).

frequency: ndarray | None = None

property frequency_range: Tuple[float, float]

The tabulated band in Hz; a scalar stage holds at every frequency.

classmethod from_amplitude_phase(name: str, input_units: str, output_units: str, frequency: Any, amplitude: Any, phase: Any, *, phase_units: str = 'deg', gain: float = 1.0, metadata: Dict[str, Any] | None = None) → ResponseStage

gain: float = 1.0

input_units: str

property is_table: bool

metadata: Dict[str, Any]

name: str

output_units: str

values: ndarray | None = None

class PyHydroGeophysX.data_processing.mt.Sounding1D(frequency: ndarray, modes: List[str], log_rho: ndarray, phase: ndarray, log_rho_err: ndarray, phase_err: ndarray, station: str = '')

Bases: object

What a 1D inversion fits: per mode, log10 rho_a and phase with errors.

frequency: ndarray

log_rho: ndarray

log_rho_err: ndarray

modes: List[str]

property n_data: int

phase: ndarray

phase_err: ndarray

station: str = ''

class PyHydroGeophysX.data_processing.mt.TimeSeriesRun(channels: ~typing.Dict[str, ~PyHydroGeophysX.data_processing.mt.timeseries.Channel], sample_rate: float, start: ~typing.Any, station: str = '', latitude: float = nan, longitude: float = nan, elevation: float = nan, declination: float = 0.0, instrument: str = '', serial: str = '', metadata: ~typing.Dict[str, ~typing.Any] = <factory>)

Bases: object

Simultaneous, gap-free samples of several channels at one rate.

array(components: Sequence[str] | None = None, dtype: Any = <class 'float'>) → ndarray

The chosen channels stacked as (n_channels, n_samples).

calibrated(components: Sequence[str] | None = None, *, band: Tuple[float, float] | None = None) → TimeSeriesRun

The channels in field units (mV/km, nT), the response removed by FFT.

Each channel is demeaned and divided by its response at every Fourier frequency; band (Hz) zeroes the spectrum outside it, which a coil’s vanishing low-frequency response needs. This holds the whole series in memory as complex numbers: for long, fast recordings select a window first.

channels: Dict[str, Channel]

property components: List[str]

declination: float = 0.0

property duration: float

Seconds the samples span, n / sample_rate.

elevation: float = nan

property end: datetime64

The time just after the last sample.

instrument: str = ''

latitude: float = nan

longitude: float = nan

metadata: Dict[str, Any]

property n_samples: int

sample_rate: float

select(components: Iterable[str]) → TimeSeriesRun

serial: str = ''

slice(first: int, last: int | None = None) → TimeSeriesRun

Samples first:last (views, not copies).

start: Any

station: str = ''

summary() → str

times() → ndarray

Each sample’s UTC time as datetime64[ns].

window(start: Any = None, end: Any = None) → TimeSeriesRun

The samples between two UTC times (either may be None).

class PyHydroGeophysX.data_processing.mt.TransferFunction(frequency: ndarray, z: ndarray | None = None, z_err: ndarray | None = None, tipper: ndarray | None = None, tipper_err: ndarray | None = None, rotation: Any = 0.0, inverse_signal_power: ndarray | None = None, z_residual_covariance: ndarray | None = None, tipper_residual_covariance: ndarray | None = None, station: str = '', latitude: float = nan, longitude: float = nan, elevation: float = nan, x: float = nan, y: float = nan, declination: float = nan, metadata: Dict[str, ~typing.Any]=<factory>)

Bases: object

Impedance and tipper of one MT site, against frequency.

z is (n, 2, 2) complex ohms, tipper (n, 1, 2) complex and dimensionless; either may be None. z_err and tipper_err are the standard errors of each real and imaginary part. inverse_signal_power (n, 2, 2) and the residual covariances - z_residual_covariance (n, 2, 2) for (Ex, Ey) and tipper_residual_covariance (n, 1, 1) for Hz - are the full error model where a file or the estimator supplies it: var(Z_ij) = residual_ii S_jj. Values that are missing are NaN.

property angular_frequency: ndarray

apparent_resistivity() → ndarray

|Z|^2 / (omega mu0) in ohm-m, (n, 2, 2).

apparent_resistivity_err() → ndarray

First-order standard error of the apparent resistivity, ohm-m.

conjugated() → TransferFunction

The site under the opposite time convention (e^{-i omega t} <-> e^{+i omega t}).

copy() → TransferFunction

declination: float = nan

determinant() → ndarray

The rotation-invariant determinant impedance sqrt(Zxx Zyy - Zxy Zyx), (n,).

The branch is the one with the phase of Zxy, so a half-space gives +45.

elevation: float = nan

frequency: ndarray

property has_impedance: bool

property has_tipper: bool

inverse_signal_power: ndarray | None = None

latitude: float = nan

longitude: float = nan

metadata: Dict[str, Any]

property n_frequencies: int

property period: ndarray

phase() → ndarray

Impedance phase in degrees, (n, 2, 2); Zyx sits near -135 for a half-space.

phase_err() → ndarray

First-order standard error of the phase, degrees.

rotated(angle: Any) → TransferFunction

The site with its x axis at azimuth angle (degrees from north, clockwise).

angle is one value or one per frequency. A missing element (Zxx and Zyy of a site given as Zxy and Zyx only) leaves missing only what it enters: a rotation by 0 or 90 degrees keeps the rest.

The errors are recomputed from the rotated covariance when the site’s errors are the covariance’s (EMTF XML, Z-files); otherwise - errors of their own, such as the jackknife’s - each rotated variance is sum (R_ik R_jl)^2 var_kl from the errors in the frame the site was read in. That assumes the elements independent, an approximation that does not depend on the path the rotations took and leaves the errors as they were at the frame they came in.

rotation: Any = 0.0

select(mask: Any) → TransferFunction

The frequencies a boolean mask or an index array selects.

sorted_by_period() → TransferFunction

Frequencies in order of increasing period.

station: str = ''

summary() → Dict[str, Any]

A JSON-safe description of the site, for logs and workflow results.

tipper: ndarray | None = None

tipper_err: ndarray | None = None

tipper_residual_covariance: ndarray | None = None

x: float = nan

Local coordinates (m), when the site belongs to a projected survey.

y: float = nan

z: ndarray | None = None

z_err: ndarray | None = None

property z_field_units: ndarray | None

The impedance in (mV/km)/nT.

z_residual_covariance: ndarray | None = None

PyHydroGeophysX.data_processing.mt.apparent_resistivity_1d(resistivity: Any, thickness: Any, frequency: Any) → Tuple[ndarray, ndarray]

(rho_a, phase_deg) of a layered earth.

PyHydroGeophysX.data_processing.mt.apply_static_shift(tf: TransferFunction, factor_x: float = 1.0, factor_y: float | None = None) → TransferFunction

Divide the apparent resistivities of the Ex row by factor_x and Ey’s by factor_y.

factor_y defaults to factor_x. The phases are unchanged; errors and covariances scale with the rows.

PyHydroGeophysX.data_processing.mt.bahr_skew(tf: TransferFunction) → ndarray

Bahr’s phase-sensitive skew; above about 0.3 the response is 3D (Bahr 1988).

PyHydroGeophysX.data_processing.mt.default_bands(n_levels: int, window_length: int = 128) → List[Band]

EMTF’s standard band set for n_levels levels of 128-point windows.

For other window lengths the harmonics are scaled by window_length / 128.

PyHydroGeophysX.data_processing.mt.estimate_static_shift(tf: TransferFunction, reference_rho_a: Any, *, period_range: Tuple[float, float] | None = None) → Dict[str, float]

Static-shift factors of the Ex and Ey rows against a reference apparent resistivity.

reference_rho_a is an array at tf.frequency or a callable of frequency. The factor is the geometric mean of rho_a / reference over period_range (s), where the reference holds - at the short periods a TEM sounding sees. Returns {"x", "y", "x_spread", "y_spread", "n"}.

PyHydroGeophysX.data_processing.mt.from_apparent_resistivity(frequency: Sequence[float], rho: Dict[str, Sequence[float]], phase_deg: Dict[str, Sequence[float]], *, rho_err: Dict[str, Sequence[float]] | None = None, phase_err_deg: Dict[str, Sequence[float]] | None = None, **site: Any) → TransferFunction

A site from apparent resistivities and phases, keyed "xy", "yx", …

|Z| = sqrt(rho omega mu0), Z = |Z| e^{i phi}. The impedance error is taken from the resistivity error, dZ = |Z| drho / (2 rho), or else from the phase error, dZ = |Z| dphi. Elements not given are NaN.

PyHydroGeophysX.data_processing.mt.gps_leap_seconds(when: Any) → int

GPS - UTC at a date, from the leap-second table (13 s before 1999).

PyHydroGeophysX.data_processing.mt.gps_seconds_to_utc(seconds: Any, *, epoch: str = 'unix', leap_seconds: int | None = None) → datetime64

UTC of a GPS-scale time stamp.

epoch="unix" reads seconds as GPS time counted from 1970-01-01, as Phoenix writes it; "gps" counts from 1980-01-06. leap_seconds is GPS - UTC; by default it comes from the table for that date.

PyHydroGeophysX.data_processing.mt.impedance_1d(resistivity: Any, thickness: Any, frequency: Any) → ndarray

Surface impedance (ohm, e^{+i omega t}) of a layered earth at each frequency.

PyHydroGeophysX.data_processing.mt.induction_arrows(tf: TransferFunction, *, convention: str = 'parkinson') → Dict[str, ndarray]

Real and imaginary induction arrows: lengths and azimuths (deg clockwise from north).

"parkinson" reverses the real arrow so it points at current concentrations; "wiese" keeps the tipper’s sign.

PyHydroGeophysX.data_processing.mt.is_edi_file(path: Any) → bool

Whether path is an EDI file: a .edi name, or a text that opens with >HEAD.

PyHydroGeophysX.data_processing.mt.is_emtf_xml_file(path: Any) → bool

PyHydroGeophysX.data_processing.mt.is_jfile(path: Any) → bool

PyHydroGeophysX.data_processing.mt.is_transfer_function_file(path: Any) → bool

PyHydroGeophysX.data_processing.mt.is_zfile(path: Any) → bool

PyHydroGeophysX.data_processing.mt.load_runs(source: Any) → List[TimeSeriesRun]

The runs a save_runs() container holds.

PyHydroGeophysX.data_processing.mt.niblett_bostick(tf: TransferFunction, *, mode: str = 'det') → Dict[str, ndarray]

depth (m) and resistivity (ohm m) by the Niblett-Bostick transform.

mode is "det" (the determinant impedance), "xy" or "yx". rho_NB = rho_a (pi / (2 phi) - 1), depth = sqrt(rho_a / (omega mu0)).

PyHydroGeophysX.data_processing.mt.occam1d(data: Any, *, mode: str = 'det', n_layers: int = 50, thickness: Sequence[float] | None = None, target_rms: float = 1.0, max_iterations: int = 30, starting_resistivity: float | None = None, static_shift: bool = False, static_shift_prior: float = 0.3, error_floor: float = 0.05, period_range: Tuple[float, float] | None = None, tem: Dict[str, Any] | None = None, tem_weight: float = 1.0, log: Callable[[str], None] | None = None, resistivity_bounds: Tuple[float, float] = (0.01, 100000.0)) → Occam1DResult

Occam’s smoothest model fitting a site’s 1D data (and a TEM sounding).

data is a TransferFunction (its mode is taken, see sounding_from_tf()) or a Sounding1D. thickness fixes the layers; otherwise n_layers are spread over the skin-depth range. static_shift adds one offset of log10 rho_a per mode. tem={"data": ..., "geometry": ..., "inversion": ...} is a TDEM sounding as PyHydroGeophysX.inversion.em1d.tdem_invert() takes it, fitted with weight tem_weight; the RMS target then applies to all data together.

PyHydroGeophysX.data_processing.mt.occam1d_anisotropic(data: TransferFunction, *, strike: float | None = None, n_layers: int = 50, thickness: Sequence[float] | None = None, target_rms: float = 1.0, max_iterations: int = 30, starting_resistivity: float | None = None, static_shift: bool = False, static_shift_prior: float = 0.3, anisotropy_weight: float = 1.0, error_floor: float = 0.05, period_range: Tuple[float, float] | None = None, log: Callable[[str], None] | None = None, resistivity_bounds: Tuple[float, float] = (0.01, 100000.0)) → Anisotropic1DResult

Occam’s smoothest azimuthally anisotropic layered model with one strike.

Each layer conducts differently along two horizontal axes, the same two in every layer: along strike (degrees clockwise from north; by default the phase tensor’s, averaged over the periods used) and across it. In that frame the layered earth’s impedance has no diagonal, and Zxy and Zyx are each a 1D isotropic response - Zxy of the resistivities along the strike, Zyx of those across it (Pek & Santos 2002, the case of horizontal principal axes sharing one azimuth). The two models share the layers and are found together: Occam’s roughness of each, plus anisotropy_weight times the difference of their log10 resistivities layer by layer, so a layer is made anisotropic only where the data ask for it.

Why the modes differ cannot be told from one site: a 2D structure with the same strike, or galvanic distortion, splits them too. A static shift moves rho_a and not the phase, so with static_shift the difference in level between the modes is put down to the shifts and the anisotropy comes from the phases’ split alone.

Pek, J. & Santos, F. A. M. (2002). Magnetotelluric impedances and parametric sensitivities for 1-D anisotropic layered media. Computers & Geosciences, 28(8), 939-950. https://doi.org/10.1016/S0098-3004(02)00014-6

PyHydroGeophysX.data_processing.mt.periodic_interference(run: TimeSeriesRun, periods: Sequence[float], *, threshold: float = 10.0) → Dict[float, List[str]]

Which channels carry something repeating exactly every one of periods seconds.

Each channel is band-passed to 0.1-20 Hz, so neither the field’s long periods nor power-line hum swamps a pulse that lasts a fraction of a second, and folded on the period; a mean waveform standing more than threshold times above what chance leaves (about 4) is interference. The periods are tried shortest first, each one found removed before the next is tried: folded on 30 s, a pulse every 4 s still stands out (at half its height, since 30 s is not a multiple of 4), and would pass for a second interference. Returns {period: [channel, ...]} for the periods found.

PyHydroGeophysX.data_processing.mt.phase_tensor(tf: TransferFunction, *, n_realizations: int = 0, seed: int = 0) → Dict[str, ndarray]

The phase tensor and its invariants (degrees), per frequency.

Returns phi (n, 2, 2) (in the TF’s frame) and phi_max, phi_min, alpha, beta (skew), strike and ellipticity; alpha and strike are azimuths clockwise from north, modulo 180 (strike has the usual 90-degree ambiguity). With n_realizations, their standard errors (*_err) come from perturbing Z by its errors.

PyHydroGeophysX.data_processing.mt.phoenix_coil_calibrations(path: Any, calibration: Any = None) → List[Dict[str, Any]]

The coil calibrations a Phoenix recording would use: {"serial", "model", "file"} each.

calibration is as read_phoenix() takes it. One entry per serial, its newest file. Files that do not read as calibrations are left out.

PyHydroGeophysX.data_processing.mt.powerline_frequency(run: TimeSeriesRun, candidates: Sequence[float] = (50.0, 60.0), *, threshold: float = 100.0) → float | None

The grid frequency whose hum run carries, of candidates: the one whose line stands highest above its floor, if more than threshold times; None when none does or none lies below the Nyquist.

PyHydroGeophysX.data_processing.mt.process_mt(runs: TimeSeriesRun | Sequence[TimeSeriesRun], *, remote: TimeSeriesRun | Sequence[TimeSeriesRun] | None = None, config: ProcessingConfig | None = None, station: str | None = None, log: Callable[[str], None] | None = None, **options: Any) → TransferFunction

Estimate a site’s impedance (and tipper, when it has Hz) from its runs.

remote holds the reference site’s runs (its Hx and Hy are used). config or keyword options set ProcessingConfig fields; log receives a line per sample rate.

PyHydroGeophysX.data_processing.mt.read_ats_header(path: Any) → Dict[str, Any]

PyHydroGeophysX.data_processing.mt.read_band_setup(path: Any) → List[Band]

An EMTF band-setup file: the number of bands, then level first last rows.

PyHydroGeophysX.data_processing.mt.read_edi(path: Any, *, sign_convention: str = 'auto') → TransferFunction

Read an EDI file into a TransferFunction.

sign_convention is "+" (the standard’s e^{+i omega t}), "-" (conjugate everything), or "auto": "+" unless the Zxy and Zyx phases lie in the quadrants of e^{-i omega t} at four frequencies in five.

PyHydroGeophysX.data_processing.mt.read_emtf_xml(path: Any) → TransferFunction

Read an EMTF XML file into a TransferFunction.

PyHydroGeophysX.data_processing.mt.read_jfile(path: Any, *, z_units: str = 'auto', sign_convention: str = 'auto') → TransferFunction

Read a J-format file; z_units is "auto", "field" or "ohm".

BIRRP writes e^{-i omega t}; sign_convention is resolved as for EDI files (resolve_sign_convention()), so such a file is conjugated.

PyHydroGeophysX.data_processing.mt.read_lemi424(path: Any, *, dipole_lengths: Dict[str, float] | None = None, components: Iterable[str] | None = None, station: str = '', calibration: Dict[str, Any] | None = None, start: Any = None, end: Any = None) → List[TimeSeriesRun]

Read LEMI-424 text files - one, a folder of them, or a list - into runs.

dipole_lengths maps "ex"/"ey" (or "e1"…) to metres. components lists the channels to keep, by LEMI name or component (default: hx, hy, hz, ex, ey). calibration maps a component to a text table (frequency, amplitude, phase in degrees) of its response.

PyHydroGeophysX.data_processing.mt.read_metronix(path: Any, *, components: Iterable[str] | None = None, calibration: Any = None, start: Any = None, end: Any = None) → List[TimeSeriesRun]

Read Metronix ATS files - one file, a meas_ folder or a site folder - into runs.

Files of one measurement folder with the same run number and sample rate form a run. calibration is a coil calibration text file, a folder of them, or a list; cal folders beside and above the measurements are searched too.

PyHydroGeophysX.data_processing.mt.read_metronix_calibration(path: Any) → Dict[str, Any]

A Metronix coil calibration text file: {"sensor", "serial", "on", "off"} stages.

PyHydroGeophysX.data_processing.mt.read_phoenix(path: Any, *, sample_rate: None | float | Sequence[float] = None, components: Iterable[str] | None = None, start: Any = None, end: Any = None, calibration: Any = None, coils: Dict[str, str] | None = None) → List[TimeSeriesRun]

Read a Phoenix MTU-5C-family recording into runs.

path is the recording folder, one of its channel folders, or a file (which reads that file’s stream, every channel of it). sample_rate picks the stream(s) - 24000 for the native .bin or a decimated rate; by default every rate present is read, highest first. start and end (UTC) limit what is read, which matters for long native recordings. calibration is a calibration file, a folder of them (searched through) or a list; by default *rxcal.json, *scal.json and binary *.scal files are looked for in the recording folder and the two above it, and in any folder there named for calibrations.

coils says which coil was on which magnetic channel, {"H1": "4684", ...} (or "hx"…), by serial or calibration file. A receiver whose coils were not detected records serial 0, so without it their responses are not known and the magnetic channels stay in volts.

Continuous streams come out as one run per gap-free stretch, segmented ones as one run per burst.

PyHydroGeophysX.data_processing.mt.read_phoenix_calibration(path: Any) → Dict[str, Any]

A Phoenix rxcal.json, scal.json or binary .scal: each channel’s response tables.

Returns {"kind": "receiver" | "sensor", "instrument_serial", "sensor_serial", "model", "channels": {tag: [ResponseStage, ...]}}; a receiver file has one table per low-pass setting.

PyHydroGeophysX.data_processing.mt.read_phoenix_header(path: Any) → Dict[str, Any]

The 128-byte header of a Phoenix time-series file, decoded.

PyHydroGeophysX.data_processing.mt.read_phoenix_legacy(path: Any, *, bands: Iterable[str] | None = None, components: Iterable[str] | None = None, start: Any = None, end: Any = None, calibration: Dict[str, Any] | None = None) → List[TimeSeriesRun]

Read a legacy Phoenix site (.TBL with .TS2-.TS5) into runs.

path is one of the site’s files. bands picks "TS3", "TS4"… (default: every TS file of the site). calibration maps a component to a ResponseStage or a text table of frequency, amplitude (mV/nT for coils), phase in degrees, which replaces the flat coil sensitivity.

PyHydroGeophysX.data_processing.mt.read_phoenix_table(path: Any) → Dict[str, Any]

The records of a legacy Phoenix .TBL file, decoded by name.

PyHydroGeophysX.data_processing.mt.read_response_table(path: Any, *, name: str | None = None, input_units: str = 'nT', output_units: str = 'mV', phase_units: str = 'deg', amplitude_per_hz: bool = False) → ResponseStage

A calibration table from a text file of frequency, amplitude, phase rows.

Lines that do not start with three numbers are skipped, so headers and comments may be anywhere. amplitude_per_hz multiplies the amplitude by frequency, for tables normalized that way (Metronix V/(nT Hz)).

PyHydroGeophysX.data_processing.mt.read_timeseries(path: Any, *, azimuths: Mapping[str, float] | None = None, **options: Any) → List[TimeSeriesRun]

Read an instrument’s MT time series - a file or a recording folder - into runs.

The format is recognized from the files: Phoenix MTU-5C/5P/8A recordings (read_phoenix()), legacy Phoenix MTU-5A TS/TBL sites (read_phoenix_legacy()), Metronix ATS (read_metronix()), Zonge Z3D (read_zonge()) and LEMI-424 text files (read_lemi424()). options go to that reader; format forces one.

azimuths ({"ex": 180, "hy": 270, ...}, degrees clockwise from north) replaces the direction the files give those channels, in every run: for a layout that differs from what the receiver was told, such as a dipole laid out to the south or a coil pointing west. Processing expresses every channel in the frame of Hx, so a reversed channel needs nothing more.

PyHydroGeophysX.data_processing.mt.read_transfer_function(path: Any, *, sign_convention: str = 'auto', z_units: str = 'auto') → TransferFunction

Read one site’s transfer functions from an EDI, EMTF XML, Z- or J-file.

sign_convention applies to EDI and J-files, which do not state it (see read_edi()), and z_units to J-files (see read_jfile()); EMTF XML and Z-files state both.

PyHydroGeophysX.data_processing.mt.read_transfer_functions(sources: str | Path | Iterable[Any], **options: Any) → List[TransferFunction]

Read every transfer-function file of a folder, or of a list of paths, in name order.

PyHydroGeophysX.data_processing.mt.read_z3d_header(path: Any) → Dict[str, Any]

The text records of a Z3D file: header, schedule and metadata, and the data offset.

PyHydroGeophysX.data_processing.mt.read_zfile(path: Any) → TransferFunction

Read an EMTF Z-file into a TransferFunction.

PyHydroGeophysX.data_processing.mt.read_zonge(path: Any, *, components: Iterable[str] | None = None, calibration: Dict[str, Any] | None = None, start: Any = None, end: Any = None, skip_seconds: int = 2) → List[TimeSeriesRun]

Read Z3D files - one file or a folder - into runs.

calibration maps a component (or a coil serial) to a ResponseStage or a text table (frequency Hz, amplitude mV/nT, phase in degrees) that replaces the coil table stored in the file.

PyHydroGeophysX.data_processing.mt.remove_periodic(run: TimeSeriesRun, period: float, *, cycles: int = 31) → TimeSeriesRun

run less whatever repeats exactly every period seconds - a receiver’s own interference.

Each channel is first high-passed by subtracting its centred moving mean over one period, which passes the harmonics k / period whole and takes out everything slower; the waveform of each cycle is then the mean of that over the cycles cycles around it, and is subtracted from the channel. The whole is one time-invariant, zero-phase comb filter, the same on every channel, so an impedance is left as it was: notches about 1 / (cycles * period) wide at the harmonics, the frequencies between them kept. (Taking each cycle’s own mean off instead turned the field’s slow drift into a staircase with steps at every cycle - a near-real “impedance” carrying the site’s static shift.) Interference timed by the receiver’s GPS clock repeats exactly; the field does not. period must be a whole number of samples; a run shorter than three periods is returned as it is.

PyHydroGeophysX.data_processing.mt.remove_periodic_bursts(runs: Sequence[TimeSeriesRun], period: float, *, cycles: int = 31) → List[TimeSeriesRun]

Bursts of a segmented recording less what repeats every period seconds.

A burst shorter than a few periods cannot be folded on itself, but bursts that start at the same phase of the period - and have the same rate and length - carry the interference at the same place. Each such burst has the mean of the other bursts of that group among the cycles around it subtracted (its own field is left out of its template), less its own mean; the field, which does not repeat, averages out of it. A group of fewer than three bursts is left as recorded. This also takes off whatever every burst of the group starts with alike.

PyHydroGeophysX.data_processing.mt.remove_powerline(run: TimeSeriesRun, frequency: float = 60.0, *, threshold: float = 10.0) → TimeSeriesRun

run with the power-line hum notched out: frequency and its harmonics below the Nyquist.

Hum is many orders of magnitude stronger than the field, and in a record whose Nyquist lies above it - 60 Hz in 150 Hz data - it leaks through the taper’s sidelobes into bands well away from it. There it passes for an impedance with apparent resistivity growing as the period and a phase of 0 or 180 degrees, and as the grid is the same at a remote site, remote referencing keeps it. Each line is found where the channels’ power, smoothed over 0.01 Hz, stands threshold times above the floor 1-3 Hz either side, and every bin of it - as wide as the grid’s frequency wandered over the run - is set to zero, in one transform of the whole run: a zero-phase notch that leaves the frequencies between the harmonics as they were. A harmonic that does not stand out is left alone.

A run shorter than a minute - a burst of a segmented recording - has bins too wide for that: hum that is not on a bin leaks across many of them and stays. There each line is fitted instead, as a sinusoid at its own frequency (the peak of the four-times zero-padded spectrum) with an amplitude that may change linearly over the run, and subtracted.

PyHydroGeophysX.data_processing.mt.robust_regression(y: ndarray, X: ndarray, R: ndarray | None = None, config: RegressionConfig | None = None, blocks: ndarray | None = None) → RegressionResult

Estimate b in y = X b robustly; R is the remote reference (default X).

blocks labels each point with its window (its start time, say): the jackknife then deletes whole windows, in groups of neighbours in label order, instead of runs of points that may split a window’s harmonics.

PyHydroGeophysX.data_processing.mt.rotation_matrix(angle_deg: Any) → ndarray

R for a frame rotated clockwise by angle_deg: shape (..., 2, 2).

A vector’s components in the rotated frame are R v: the new x axis, at azimuth theta, is (cos theta, sin theta) in the old (north, east) frame.

PyHydroGeophysX.data_processing.mt.runs_from_payload(payload: Sequence[Dict[str, Any]]) → List[TimeSeriesRun]

The runs runs_to_payload() took apart.

PyHydroGeophysX.data_processing.mt.runs_to_payload(runs: Sequence[TimeSeriesRun]) → List[Dict[str, Any]]

Runs as plain values and arrays, for PyHydroGeophysX.data_processing.run_inputs.save_container().

PyHydroGeophysX.data_processing.mt.saturated(run: TimeSeriesRun, components: Sequence[str] | None = None, *, level: float = 0.995, guard: int = 2) → ndarray | None

Which samples of run have one of components at the recorder’s full scale.

Only channels whose reader knows the full scale (metadata["full_scale"]) count; None when none does. Each clipped sample is widened by guard samples either side. A clipped stretch is the recorder’s limit, not the field: windows holding it give an impedance that is plausible and wrong.

PyHydroGeophysX.data_processing.mt.save_runs(destination: Any, runs: Sequence[TimeSeriesRun])

Write runs as one compressed container (.npz); returns its path.

PyHydroGeophysX.data_processing.mt.sensitivity_1d(resistivity: Any, thickness: Any, frequency: Any) → Tuple[ndarray, ndarray, ndarray]

Z and the derivatives of log10 rho_a and phase (rad) to log rho_j.

Returns (Z, d_log10_rho_a, d_phase), the latter two (n_freq, n_layers).

PyHydroGeophysX.data_processing.mt.sounding_from_tf(tf: TransferFunction, *, mode: str = 'det', error_floor: float = 0.05, period_range: Tuple[float, float] | None = None) → Sounding1D

The 1D data of a site: apparent resistivity and phase of a mode, with errors.

error_floor is the least relative error of the impedance (5%: about 0.043 in log10 rho_a and 1.4 degrees of phase). Phases are folded into the first quadrant; values outside it are dropped.

PyHydroGeophysX.data_processing.mt.static_shift_from_layers(tf: TransferFunction, resistivity: Any, thickness: Any, *, period_range: Tuple[float, float] | None = None) → Dict[str, float]

Static-shift factors against the MT response of a layered model, such as a TEM inversion’s.

The model’s apparent resistivity is computed at the site’s frequencies (forward1d.apparent_resistivity_1d()); period_range should keep to the periods whose skin depth the model resolves.

PyHydroGeophysX.data_processing.mt.swift_skew(tf: TransferFunction) → ndarray

|Zxx + Zyy| / |Zxy - Zyx|; near 0 for 1D and 2D (Swift 1967).

PyHydroGeophysX.data_processing.mt.swift_strike(tf: TransferFunction) → ndarray

The strike (deg clockwise from north, +-90 ambiguity) that minimizes |Zxx|^2 + |Zyy|^2.

PyHydroGeophysX.data_processing.mt.timeseries_format(path: Any) → str | None

Which instrument wrote path (a file or a folder), or None.

PyHydroGeophysX.data_processing.mt.water_content_profile(result: Occam1DResult, *, rhos: Any, n: Any, porosity: Any, sigma_sur: Any = 0.0) → Dict[str, ndarray]

Volumetric water content of each layer from its resistivity (Waxman-Smits).

The petrophysical link the rest of PyHydroGeophysX uses (PyHydroGeophysX.petrophysics.resistivity_models.resistivity_to_water_content()): rhos the saturated resistivity, n the saturation exponent, porosity and sigma_sur the surface conductivity - each one value or one per layer.

PyHydroGeophysX.data_processing.mt.write_edi(tf: TransferFunction, path: Any) → Path

Write tf as a standard EDI file, impedances in (mV/km)/nT, e^{+i omega t}.

The impedance blocks carry Z**.VAR as the squared standard error, the tipper the experimental T*.EXP blocks every program reads, and both rotation blocks the site’s rotation.

PyHydroGeophysX.data_processing.mt.write_emtf_xml(tf: TransferFunction, path: Any) → Path

Write tf as an EMTF XML file, impedances in (mV/km)/nT, e^{+i omega t}.

The full error model is written when the site carries it; otherwise only the variances, err^2. One rotation angle is written for the whole site, so a site whose frequencies are in different frames is first rotated to the frame of its first frequency.