2.8.1.2. pycsamt.site.compute#

Functions

phase_slope(obj, band, *[, api])

Estimate phase slopes within a frequency band.

res_at_freq(obj, freq, *[, how, api])

Evaluate apparent resistivity at a target frequency.

strike_estimate(obj, *[, method, api])

Estimate a strike angle from impedance tensors.

tipper_magnitude(obj, *[, per_freq, api])

Summarize or tabulate tipper magnitudes.
pycsamt.site.compute.strike_estimate(obj, *, method='swift', api=None)[source]

Estimate a strike angle from impedance tensors.

This computes a 2D geoelectric strike angle in degrees. The routine supports either a single site or many sites. For a single site a scalar angle is returned. For multiple sites a pandas.DataFrame is returned with one row per site.

Parameters:

  • obj (Any) – A single EDI-like object (e.g. EDIFile) or an iterable of EDI-like objects. Each item must expose a .Z section of shape (n_freq, 2, 2) or be convertible to such.

  • method (str, optional) –

    Strike method. Allowed values are:

    • "swift" (default): grid search over 0..179 degrees that minimizes the diagonal power after rotation.

    • "groom": alias of "swift" in this lightweight mode.

    • "phase_diff": heuristic that returns 0 or 90 degrees based on the relative magnitude of off-diagonals.

  • api (bool | None)

Returns:

If obj is a single site, returns a float angle in degrees within [0, 180). If obj is iterable, returns a DataFrame with columns:

station, method, theta_deg.

Return type:

float or pandas.DataFrame

Notes

The Swift-style criterion rotates the impedance tensor \(Z\) by a test angle \(\\theta\) and minimizes

\[\begin{split}J(\\theta) = \lvert Z'_{xx} \rvert^2 + \lvert Z'_{yy} \rvert^2 ,\end{split}\]

where \(Z'\) is the rotated tensor. The returned angle is the argmin over a 1 degree grid in 0..179.

The "phase_diff" fallback returns 0 if median \(\\lvert Z_{xy} \rvert \ge \lvert Z_{yx} \rvert\), else 90. It is intended for degraded or sparse data.

This function does not alter data. If you need deterministic behavior for incomplete arrays, consider preparing tensors with pycsamt.site.edit.fill_missing().

Examples

Single site, Swift estimate:

>>> from pycsamt.seg.edi import EDIFile
>>> from pycsamt.site import compute as cmp, edit as ed
>>> edf = EDIFile("S01.edi")
>>> edf = ed.fill_missing(
...     edf, how="zero", components=("Z",), inplace=False
... )
>>> ang = cmp.strike_estimate(edf, method="swift")
>>> 0.0 <= ang < 180.0
True

Many sites, returning a DataFrame:

>>> e1 = EDIFile("S01.edi")
>>> e2 = EDIFile("S02.edi")
>>> df = cmp.strike_estimate([e1, e2], method="phase_diff")
...
>>> list(df.columns)
['station', 'method', 'theta_deg']

See also

pycsamt.site.edit.rotate

Rotate site tensors by a user angle.

pycsamt.site.compute.phase_slope

Phase slope diagnostic over a frequency band.

References

[strike-estimate-1]

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

[strike-estimate-2]

Groom, R. W., and R. C. Bailey, 1989. Decomposition of magnetotelluric impedance tensors in the presence of local three dimensional galvanic distortion. JGR.

pycsamt.site.compute.res_at_freq(obj, freq, *, how='nearest', api=None)[source]

Evaluate apparent resistivity at a target frequency.

Computes apparent resistivity for the \(Z_{xy}\) and \(Z_{yx}\) components at a requested frequency. Works with a single site or a collection. For a single site, a dict is returned. For multiple sites, a pandas.DataFrame is returned.

Parameters:

  • obj (Any) – A single EDI-like object (e.g. EDIFile) or an iterable of such objects. Each item must expose a .Z section of shape (n_freq, 2, 2) and a frequency vector.

  • freq (float) – Query frequency in Hz.

  • how (str, optional) –

    Selection mode:

    • "nearest" (default): choose the nearest available frequency in the site data and report that value.

    • "interp": linearly interpolate resistivity versus frequency using numpy.interp. Interpolation occurs on linear frequency, not log frequency.

  • api (bool | None)

Returns:

If obj is a single site, returns a dictionary with keys "res_xy", "res_yx", "f_used". If obj is iterable, returns a DataFrame with columns station, res_xy, res_yx, f_used.

Return type:

dict or pandas.DataFrame

Notes

Apparent resistivity \(\\rho_a\) is computed as

\[\begin{split}\rho_a = \frac{\lvert Z \rvert^2} {\mu_0\,2\pi\\,f} ,\end{split}\]

where \(Z\) is the complex impedance for the selected component, \(\\mu_0\) is the magnetic permeability of free space, and \(f\) is frequency in Hz.

When how="interp", the function first computes \(\rho_a\) at all native frequencies, then interpolates the result to the query frequency using linear interpolation in frequency. If the frequency vector or impedance is missing, NaN values are returned.

Examples

Single site, nearest selection:

>>> from pycsamt.seg.edi import EDIFile
>>> from pycsamt.site import compute as cmp, edit as ed
>>> edf = EDIFile("S01.edi")
>>> edf = ed.fill_missing(
...     edf, how="zero", components=("Z",), inplace=False
... )
>>> out = cmp.res_at_freq(edf, 150.0, how="nearest")
...
>>> set(out.keys()) == {"res_xy", "res_yx", "f_used"}
...
True

Single site, interpolated:

>>> out = cmp.res_at_freq(edf, 150.0, how="interp")
...
>>> out["f_used"]
150.0

Many sites, DataFrame:

>>> e1 = EDIFile("S01.edi")
>>> e2 = EDIFile("S02.edi")
>>> df = cmp.res_at_freq([e1, e2], 1.0, how="interp")
...
>>> list(df.columns)
['station', 'res_xy', 'res_yx', 'f_used']

See also

pycsamt.site.compute.strike_estimate

Estimate 2D strike angle from Z.

pycsamt.site.edit.select_freq

Subset site data by frequency criteria.

References

[res-at-freq-1]

Vozoff, K., 1991. The magnetotelluric method. In Electromagnetic methods in applied geophysics.

[res-at-freq-2]

Simpson, F., and K. Bahr, 2005. Practical Magnetotellurics. Cambridge University Press.

pycsamt.site.compute.phase_slope(obj, band, *, api=None)[source]

Estimate phase slopes within a frequency band.

For each site, this computes the least-squares slope of phase (degrees) versus \(\log_{10}(f)\) over the requested band. Two slopes are reported, one for \(Z_{xy}\) and one for \(Z_{yx}\).

If a single site is provided, a dictionary is returned. If an iterable of sites is provided, a pandas.DataFrame is returned with one row per station.

Parameters:

  • obj (Any) – A single EDI-like object (e.g. EDIFile) or an iterable of such objects.

  • band (tuple of float) – Inclusive frequency band as (fmin, fmax) in Hz. The order does not matter; the function uses the numeric min and max.

  • api (bool | None)

Returns:

Single site -> {"slope_xy": float, "slope_yx": float}. Multi-site -> DataFrame with columns ["station", "slope_xy", "slope_yx"].

Return type:

dict or pandas.DataFrame

Notes

The phase series for each off-diagonal component is computed as

\[\begin{split}\phi(f) = \operatorname{angle}(Z(f)) \times 180/\\pi ,\end{split}\]

then a straight line is fit

\[\begin{split}\phi(f) \approx a\\,\\log_{10}(f) + b\end{split}\]

using numpy.polyfit(x, y, 1) where \(x=\\log_{10}(f)\). The reported slope is \(a\) in units of degrees per decade.

Rows or sites with missing data in the band are reported as NaN. The function does not unwrap phase.

Examples

Single site:

>>> from pycsamt.seg.edi import EDIFile
>>> from pycsamt.site import compute as cmp, edit as ed
>>> edf = EDIFile("S01.edi")
>>> edf = ed.fill_missing(
...     edf, how="zero", components=("Z",), inplace=False
... )
>>> out = cmp.phase_slope(edf, band=(1.0, 1000.0))
...
>>> set(out.keys()) == {"slope_xy", "slope_yx"}
...
True

Many sites:

>>> e1 = EDIFile("S01.edi")
>>> e2 = EDIFile("S02.edi")
>>> df = cmp.phase_slope([e1, e2], band=(0.1, 10.0))
...
>>> list(df.columns)
['station', 'slope_xy', 'slope_yx']

See also

pycsamt.site.compute.strike_estimate

Strike angle by Swift-style criterion.

pycsamt.site.compute.res_at_freq

Apparent resistivity at a target frequency.

References

[phase-slope-1]

Simpson, F., and K. Bahr, 2005. Practical Magnetotellurics. Cambridge University Press.

[phase-slope-2]

Vozoff, K., 1991. The magnetotelluric method. In Electromagnetic methods in applied geophysics.

pycsamt.site.compute.tipper_magnitude(obj, *, per_freq=False, api=None)[source]

Summarize or tabulate tipper magnitudes.

Computes the magnitude of the tipper vector per frequency as

\[\lVert \mathbf{T} \rVert = \sqrt{\lvert T_x \rvert^2 + \lvert T_y \rvert^2} ,\]

where \(T_x, T_y\) are the complex tipper components. The result can be returned as per-frequency values or summarized statistics.

For a single site, returns a dict. For an iterable of sites, returns a pandas.DataFrame.

Parameters:

  • obj (Any) – A single EDI-like object (e.g. EDIFile) or an iterable of such objects. The tipper may be attached as ed.Tip, ed.T, or ed.TIP and must expose a 2-component array shaped (n_freq, 2) or (n_freq, 1, 2).

  • per_freq (bool, optional) – If False (default), return summary statistics (mean, median, max). If True, return per-frequency values.

  • api (bool | None)

Returns:

Single site:

  • per_freq=False -> {"mean", "median", "max"}

  • per_freq=True -> {"freq", "mag"}

Multi-site:

  • per_freq=False -> DataFrame with columns ["station", "mean", "median", "max"]

  • per_freq=True -> DataFrame with columns ["station", "freq", "mag"]

Return type:

dict or pandas.DataFrame

Notes

If the site has no tipper section, summary statistics are NaN and per-frequency mode yields an empty result for that site. To initialize missing arrays, consider pycsamt.site.edit.fill_missing() with components=("Tip",).

Frequencies are reported from the site frequency vector. The function assumes the tipper array and frequency vector are aligned along their first dimension.

Examples

Single site, summary stats:

>>> from pycsamt.seg.edi import EDIFile
>>> from pycsamt.site import compute as cmp, edit as ed
>>> edf = EDIFile("S01.edi")
>>> edf = ed.fill_missing(
...     edf, how="zero", components=("Tip",), inplace=False
... )
>>> s = cmp.tipper_magnitude(edf, per_freq=False)
...
>>> set(s.keys()) == {"mean", "median", "max"}
...
True

Single site, per-frequency:

>>> out = cmp.tipper_magnitude(edf, per_freq=True)
...
>>> list(out.keys())
['freq', 'mag']

Many sites, summary:

>>> e1 = EDIFile("S01.edi")
>>> e2 = EDIFile("S02.edi")
>>> df = cmp.tipper_magnitude([e1, e2], per_freq=False)
...
>>> list(df.columns)
['station', 'mean', 'median', 'max']

See also

pycsamt.site.edit.fill_missing

Initialize or sanitize Z/Tip arrays in a site.

pycsamt.site.compute.res_at_freq

Apparent resistivity at a target frequency.

References

[tipper-magnitude-1]

Simpson, F., and K. Bahr, 2005. Practical Magnetotellurics. Cambridge University Press.