pycsamt.zonge.proc_utils#

pycsamt.zonge.proc_utils#

This module provides the core filtering algorithms used for static shift correction, mirroring the methods described in the Zonge ASTATIC program manual.

Functions

ama(-> ~pandas.DataFrame  -> ~numpy.ndarray)

Apply an Adaptive Moving Average (AMA) filter.

flma(-> ~pandas.DataFrame  -> ~numpy.ndarray)

Apply a Fixed-Length Moving Average (FLMA) filter.

get_reference_frequency(df[, mode, ...])

Determine a suitable reference frequency for static shift.

get_skew(df)

Calculate Swift's skew for the impedance tensor.

get_strike(df)

Calculate the geoelectric strike angle.

interpolate_to_log_space(df, *[, freq_min, ...])

Interpolate AVG data onto a regular log-spaced grid.

prepare_strike_frame(*[, z_frame, df, ...])

Build a strike-ready frame with complex z.

smooth_rho_from_phase(df, *[, smoothing_factor])

Smooth apparent resistivity using the Hilbert transform.

tma(-> ~pandas.DataFrame  -> ~numpy.ndarray)

Apply a Trimmed Moving Average (TMA) filter.
pycsamt.zonge.proc_utils.tma(rho_profile: str, *, data: DataFrame, window_size: int = 5, trim_proportion: float = 0.2) → DataFrame[source]#
pycsamt.zonge.proc_utils.tma(rho_profile: ndarray, *, data: None = None, window_size: int = 5, trim_proportion: float = 0.2) → ndarray
pycsamt.zonge.proc_utils.tma(rho_profile: Series, *, data: None = None, window_size: int = 5, trim_proportion: float = 0.2) → Series

Apply a Trimmed Moving Average (TMA) filter.

This filter is designed to remove single-station static offsets while preserving broader geological trends. It can operate on a pandas Series, a NumPy array, or a column within a pandas DataFrame.

Parameters:

  • rho_profile (pd.Series, np.ndarray, or str) –

    The apparent resistivity data to be filtered. This can be: - A pandas Series of resistivity values. - A 1D NumPy array of resistivity values. - The string name of the column to use if data is

    provided.

  • data (pd.DataFrame, optional) – A DataFrame containing the resistivity data. Required if rho_profile is provided as a string.

  • window_size (int, default 5) – The size of the moving window for the filter. Must be an odd integer.

  • trim_proportion (float, default 0.2) – The proportion of observations to trim from each end of the window before computing the mean. For a window of 5, 0.2 trims 1 value from each end (the min and max).

Returns:

The smoothed resistivity profile, returned in the same format as the input (rho_profile). If a DataFrame was passed to data, the function returns a new DataFrame with an added column for the smoothed data.

Return type:

pd.Series, np.ndarray, or pd.DataFrame

Raises:

ValueError – If rho_profile is a string but data is not provided, or if window_size is not an odd integer.

pycsamt.zonge.proc_utils.flma(z_profile: str, stations: str, dipole_length: float, *, data: DataFrame, filter_width_dipoles: float = 5.0) → DataFrame[source]#
pycsamt.zonge.proc_utils.flma(z_profile: ndarray, stations: ndarray, dipole_length: float, *, data: None = None, filter_width_dipoles: float = 5.0) → ndarray
pycsamt.zonge.proc_utils.flma(z_profile: Series, stations: Series, dipole_length: float, *, data: None = None, filter_width_dipoles: float = 5.0) → Series

Apply a Fixed-Length Moving Average (FLMA) filter.

This filter smooths complex impedance data using a spatial Hanning window whose width is a fixed multiple of the receiver dipole length.

Parameters:

  • z_profile (pd.Series, np.ndarray, or str) – The complex impedance data (Z) to be filtered. Can be: - A pandas Series of complex impedance values. - A 1D NumPy array of complex impedance values. - The string name of the column if data is provided.

  • stations (pd.Series, np.ndarray, or str) – The station locations (coordinates). Must correspond to the z_profile data. Can be a Series, array, or column name if data is provided.

  • dipole_length (float) – The length of the E-field receiver dipole in the same units as the station locations.

  • data (pd.DataFrame, optional) – A DataFrame containing the impedance and station data. Required if z_profile or stations are strings.

  • filter_width_dipoles (float, default 5.0) – The total width of the Hanning window, expressed in multiples of the dipole_length.

Returns:

The smoothed complex impedance profile, returned in the same format as the input. If a DataFrame was passed to data, a new DataFrame with an added column for the smoothed data is returned.

Return type:

pd.Series, np.ndarray, or pd.DataFrame

pycsamt.zonge.proc_utils.ama(z_profile: str, stations: str, dipole_length: float, frequency: float, *, data: DataFrame, skin_depth_factor: float = 2.0, iterations: int = 3) → DataFrame[source]#
pycsamt.zonge.proc_utils.ama(z_profile: ndarray, stations: ndarray, dipole_length: float, frequency: float, *, data: None = None, skin_depth_factor: float = 2.0, iterations: int = 3) → ndarray
pycsamt.zonge.proc_utils.ama(z_profile: Series, stations: Series, dipole_length: float, frequency: float, *, data: None = None, skin_depth_factor: float = 2.0, iterations: int = 3) → Series

Apply an Adaptive Moving Average (AMA) filter.

This filter smooths complex impedance data using a spatial Hanning window whose width is adapted iteratively based on the local skin depth.

Parameters:

  • z_profile (pd.Series, np.ndarray, or str) – The complex impedance data (Z) to be filtered.

  • stations (pd.Series, np.ndarray, or str) – The station locations (coordinates).

  • dipole_length (float) – The length of the E-field receiver dipole.

  • frequency (float) – The frequency of the data in Hz.

  • data (pd.DataFrame, optional) – A DataFrame containing the impedance and station data. Required if z_profile or stations are strings.

  • skin_depth_factor (float, default 2.0) – The filter width, expressed as a multiple of the local skin depth.

  • iterations (int, default 3) – The number of iterations to perform to allow the filter width to adapt to the smoothed data.

Returns:

The smoothed complex impedance profile, returned in the same format as the input.

Return type:

pd.Series, np.ndarray, or pd.DataFrame

pycsamt.zonge.proc_utils.interpolate_to_log_space(df, *, freq_min=None, freq_max=None, num_points=50, interp_kind='cubic')[source]#

Interpolate AVG data onto a regular log-spaced grid.

This function resamples the apparent resistivity and phase data for each station and component onto a new, logarithmically spaced frequency axis. This is a common preprocessing step for inversion and modeling.

Parameters:

  • df (pandas.DataFrame) – The input DataFrame containing the AVG data. Must include ‘station’, ‘freq’, ‘rho’, and ‘phase’ columns.

  • freq_min (float, optional) – The minimum frequency for the new grid. If None, it is inferred from the minimum frequency in the data.

  • freq_max (float, optional) – The maximum frequency for the new grid. If None, it is inferred from the maximum frequency in the data.

  • num_points (int, default 50) – The number of logarithmically spaced points to create between freq_min and freq_max.

  • interp_kind (str, default 'cubic') – The kind of interpolation to perform, passed to scipy.interpolate.interp1d. Common options are ‘slinear’, ‘quadratic’, and ‘cubic’.

Returns:

A new DataFrame with the data interpolated onto the specified logarithmic frequency grid.

Return type:

pandas.DataFrame

pycsamt.zonge.proc_utils.smooth_rho_from_phase(df, *, smoothing_factor=0.1)[source]#

Smooth apparent resistivity using the Hilbert transform.

This function reconstructs a smoother apparent resistivity curve from the impedance phase data, leveraging the causal nature of magnetotelluric data.

Parameters:

  • df (pandas.DataFrame) – The input DataFrame. Must include ‘station’, ‘freq’, ‘rho’, and ‘phase’ columns.

  • smoothing_factor (float, default 0.1) – The smoothing factor for the UnivariateSpline applied to the phase data before the Hilbert transform. A smaller value results in less smoothing.

Returns:

A new DataFrame with a ‘rho_smoothed’ column containing the reconstructed, smoother apparent resistivity values.

Return type:

pandas.DataFrame

pycsamt.zonge.proc_utils.get_reference_frequency(df, mode='auto', *, qc_column='pc_rho', qc_threshold=20.0)[source]#

Determine a suitable reference frequency for static shift.

This utility automatically selects a reference frequency based on the principle of using the “highest frequency with clean data,” as recommended in the ASTATIC manual.

Parameters:

  • df (pandas.DataFrame) – The input DataFrame containing the AVG data. Must include ‘station’, ‘freq’, and a quality control column.

  • mode ({'auto'} or float, default 'auto') –

    The method for determining the reference frequency. - ‘auto’: Automatically select the frequency based on QC

    metrics.

    • float: A specific frequency value to use, which will be returned directly.

  • qc_column (str, default 'pc_rho') – The canonical name of the quality control column to use for filtering “clean” data points.

  • qc_threshold (float, default 20.0) – The threshold for the qc_column. Data points where the value in this column is below the threshold are considered clean.

Returns:

The determined reference frequency in Hz.

Return type:

float

Notes

The ‘auto’ mode follows a robust, data-driven approach: 1. It filters the dataset to include only “clean” data points

(where qc_column < qc_threshold).

  1. For each station, it finds the highest frequency that has clean data.

  2. It then returns the median of these highest frequencies. The median is used to provide a stable estimate that is robust against outlier stations.

If no data points meet the clean criteria, the function will issue a warning and fall back to using the absolute maximum frequency present in the entire dataset.

pycsamt.zonge.proc_utils.get_skew(df)[source]#

Calculate Swift’s skew for the impedance tensor.

Skew is a dimensionless parameter that quantifies the three-dimensionality of the subsurface conductivity structure. It is a crucial diagnostic tool for determining whether a 1D or 2D interpretation is appropriate for the data.

Parameters:

df (pandas.DataFrame) – A tidy DataFrame containing the impedance tensor components. Must include ‘station’, ‘freq’, ‘comp’, and a complex ‘z’ column.

Returns:

A DataFrame with columns for ‘station’, ‘freq’, and the calculated ‘skew’.

Return type:

pandas.DataFrame

Notes

The function calculates the conventional skew definition proposed by Swift [1]_. The formula is:

\[\text{skew} = \frac{|Z_{xx} + Z_{yy}|}{|Z_{xy} - Z_{yx}|}\]

Generally, skew values are interpreted as follows: - skew < 0.1: Data can be considered 1D. - 0.1 <= skew <= 0.3: Data is likely 2D. - skew > 0.3: Data is likely 3D or affected by significant

local distortion.

Examples

>>> from pycsamt.zonge import AMTAVG
>>> from pycsamt.zonge.proc_utils import get_skew
>>> avg = AMTAVG.from_file('data/avg/K2.avg')
>>> skew_df = get_skew(avg.z.frame)
>>> print(skew_df.head())

References

See also

get_strike

Calculate the geoelectric strike angle.

pycsamt.zonge.proc_utils.get_strike(df)[source]#

Calculate the geoelectric strike angle.

This function determines the principal axis direction of the impedance tensor for each station and frequency. The strike is the angle to which the coordinate system must be rotated to minimize the diagonal components of the impedance tensor (\(Z_{xx}\), \(Z_{yy}\)), which is a common assumption for 2D geological structures.

Parameters:

df (pandas.DataFrame) – A tidy DataFrame containing the impedance tensor components. Must include ‘station’, ‘freq’, ‘comp’, and a complex ‘z’ column from which the tensor elements can be pivoted.

Returns:

A DataFrame with columns for ‘station’, ‘freq’, and the calculated ‘strike_angle’ in degrees clockwise from North.

Return type:

pandas.DataFrame

Notes

The calculation is performed using a common tensor decomposition method based on the real parts of the impedance tensor components [1]_. The formula used is:

\[\theta_s = \frac{1}{2} \arctan\left( \frac{\text{Re}(Z_{xy} + Z_{yx})} {\text{Re}(Z_{xx} - Z_{yy})} \right)\]

This angle represents the orientation of the primary 2D geological structure.

Examples

>>> from pycsamt.zonge import AMTAVG
>>> from pycsamt.zonge.proc_utils import get_strike
>>> avg = AMTAVG.from_file('data/avg/K2.avg')
>>> # The Z component's frame already has the required columns
>>> strike_df = get_strike(avg.z.frame)
>>> print(strike_df.head())

References

See also

get_skew

Calculate the dimensionality indicator (skew).

pycsamt.zonge.avg.AMTAVG.rotate

Rotate the tensor by a given angle.

pycsamt.zonge.proc_utils.prepare_strike_frame(*, z_frame=None, df=None, prefer='z', station_col='station', freq_col='freq', comp_col='comp', z_col='z', rho_col='rho', phase_col='phase', phase_unit='auto', drop_na=True, na_policy='any', components=None, ensure_sorted=True, copy=True, mu0=1.2566370614359173e-06)[source]#

Build a strike-ready frame with complex z.

Prefers an input frame that already contains a complex impedance column. If not available, derives z from rho, phase and freq:

\[|Z| = \sqrt{\mu_0 \, \omega \, \rho_a},\quad Z = |Z|\,[\cos\phi + i\sin\phi]\]
Parameters:

  • z_frame (DataFrame, optional) – Frame that already holds a complex column z.

  • df (DataFrame, optional) – Long frame with columns for apparent resistivity, phase and frequency (rho, phase, freq).

  • prefer ({"z","df"}, default "z") – Source priority when both frames are given.

  • station_col (str) – Column names in provided frames.

  • freq_col (str) – Column names in provided frames.

  • comp_col (str) – Column names in provided frames.

  • z_col (str) – Column names in provided frames.

  • rho_col (str) – Column names in provided frames.

  • phase_col (str) – Column names in provided frames.

  • phase_unit ({"auto","deg","mrad","rad"}, default "auto") – Unit for phase. "auto" guesses: - |phase|max <= π*1.5 → rad - else if <= 180 → deg - else → mrad

  • drop_na (bool, default True) – Drop NA rows before computing.

  • na_policy ({"any","all"}, default "any") – Row-drop policy when drop_na=True.

  • components (iterable of str, optional) – If given, keep only these components.

  • ensure_sorted (bool, default True) – Sort output by (station, freq, comp).

  • copy (bool, default True) – Work on a copy of the source frame(s).

  • mu0 (float, default 4π·1e−7) – Magnetic permeability of free space.

Returns:

Columns: station, freq, comp, z.

Return type:

DataFrame

Raises:

ProcessingError – If required columns are missing or inputs are not provided.