Computed Diagnostics#

The pycsamt.site.compute module provides lightweight station diagnostics for EDI-like sites and site collections. These functions do not modify the input data. They read impedance, frequency, and tipper arrays and return small numerical summaries that are useful during quality control, survey preparation, and interpretation.

Use computed diagnostics before heavier processing when you need to answer questions such as:

  • does each station have a usable impedance tensor?

  • what strike angle is suggested by the off-diagonal tensor structure?

  • what apparent resistivity is observed near a frequency of interest?

  • does phase vary smoothly across the band used for inversion?

  • is the tipper response small, large, missing, or spatially variable?

Input Contract#

The diagnostics accept either one EDI-like object or an iterable of EDI-like objects. In practice this means:

The object must expose the data needed by the selected diagnostic:

Diagnostic

Required arrays

Missing-data behaviour

strike_estimate()

Frequency vector and impedance tensor shaped (n_freq, 2, 2).

Returns NaN when tensor or frequency data are missing.

res_at_freq()

Frequency vector and impedance tensor shaped (n_freq, 2, 2).

Returns NaN resistivity values when required data are missing.

phase_slope()

Frequency vector and impedance tensor shaped (n_freq, 2, 2).

Returns NaN slopes when the requested band is empty or unusable.

tipper_magnitude()

Frequency vector and tipper array shaped (n_freq, 2).

Returns NaN summary values, or an empty per-frequency result.

For deterministic diagnostics on incomplete data, prepare the site first with pycsamt.site.edit.fill_missing() or filter the survey with pycsamt.site.selection.keep_finite_z().

Return Types#

The module follows one consistent rule:

Input

Return value

Single site

A scalar or dictionary.

Iterable of sites

A pandas.DataFrame.

Iterable of sites with api=True

A pyCSAMT APIFrame wrapper around the DataFrame, when API view support is enabled.

This makes the functions convenient both in notebooks and in automated pipeline steps.

 1from pycsamt.seg.edi import EDIFile
 2from pycsamt.site.compute import res_at_freq, strike_estimate
 3
 4site = EDIFile("data/edi/S01.edi")
 5
 6theta = strike_estimate(site)
 7rho = res_at_freq(site, 100.0)
 8
 9print(theta)
10print(rho["res_xy"], rho["res_yx"], rho["f_used"])

For a collection:

1from pycsamt.site import Sites
2from pycsamt.site.compute import phase_slope
3
4sites = Sites.from_path("data/edi")
5table = phase_slope(sites, band=(1.0, 1000.0))
6
7print(table.head())

Diagnostic Map#

Function

Output

Main use

strike_estimate()

theta_deg in degrees.

Quick geoelectric strike estimate for tensor rotation checks and 2-D assumption screening.

res_at_freq()

res_xy, res_yx, and f_used.

Compare stations at one frequency, select target frequency slices, or build compact QC tables.

phase_slope()

slope_xy and slope_yx in degrees per frequency decade.

Detect abrupt phase behaviour, band-edge problems, or unstable curves.

tipper_magnitude()

Mean, median, max, or per-frequency tipper magnitude.

Evaluate vertical magnetic transfer response and possible 3-D structure indicators.

Strike Estimate#

strike_estimate() estimates a 2-D geoelectric strike angle from the impedance tensor. It supports three method names:

"swift"

Grid search from 0 to 179 degrees. The selected angle minimizes diagonal tensor power after rotation.

"groom"

Currently an alias for "swift" in this lightweight diagnostic layer.

"phase_diff"

A coarse fallback returning either 0 or 90 degrees based on the relative median magnitude of \(Z_{xy}\) and \(Z_{yx}\).

For the Swift-style criterion, the tensor is rotated by a trial angle \(\theta\), then the diagonal power is summarized:

\[J(\theta) = \operatorname{median} \left( |Z'_{xx}|^2 + |Z'_{yy}|^2 \right).\]

The reported strike is the angle that minimizes \(J(\theta)\) over the 1-degree search grid.

 1from pycsamt.site import Sites
 2from pycsamt.site.compute import strike_estimate
 3from pycsamt.site.edit import fill_missing
 4
 5sites = Sites.from_path("data/edi")
 6
 7# Optional: make missing tensor entries explicit before diagnostics.
 8prepared = [
 9    fill_missing(site, how="zero", components=("Z",), inplace=False)
10    for site in sites
11]
12
13strike_table = strike_estimate(prepared, method="swift")
14print(strike_table[["station", "theta_deg"]])

Interpretation notes:

  • strike estimates are screening values, not final structural interpretation;

  • unstable or sparse tensors can produce misleading angles;

  • the result is in the interval \([0, 180)\);

  • compare strike values across neighbouring stations before rotating an entire survey line;

  • use pycsamt.site.edit.rotate() or pycsamt.site.edit.rotate_all() after deciding on a rotation angle.

Apparent Resistivity At One Frequency#

res_at_freq() evaluates apparent resistivity for the off-diagonal components \(Z_{xy}\) and \(Z_{yx}\) at a requested frequency.

The apparent resistivity is computed as:

\[\rho_a = \frac{|Z|^2}{\mu_0\,2\pi f},\]

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

Two frequency-selection modes are available:

how="nearest"

Select the nearest available native frequency and report that value in f_used.

how="interp"

Compute resistivity at all native frequencies, then interpolate to the requested frequency using linear interpolation in frequency.

 1from pycsamt.site import Sites
 2from pycsamt.site.compute import res_at_freq
 3
 4sites = Sites.from_path("data/edi")
 5
 6near = res_at_freq(sites, 150.0, how="nearest")
 7interp = res_at_freq(sites, 150.0, how="interp")
 8
 9print(near[["station", "res_xy", "res_yx", "f_used"]])
10print(interp[["station", "res_xy", "res_yx", "f_used"]])

Use nearest when you want to preserve the exact sampled frequency axis. Use interp when you need one common comparison frequency across stations whose frequency grids are slightly different.

Phase Slope#

phase_slope() summarizes how phase changes across a frequency band. For each site, it computes phase for the off-diagonal components:

\[\phi(f) = \arg(Z(f)) \times 180 / \pi,\]

then fits a straight line against log-frequency:

\[\phi(f) \approx a \log_{10}(f) + b.\]

The reported slope \(a\) is measured in degrees per decade. The function returns one slope for \(Z_{xy}\) and one for \(Z_{yx}\).

 1from pycsamt.site import Sites
 2from pycsamt.site.compute import phase_slope
 3
 4sites = Sites.from_path("data/edi")
 5
 6slopes = phase_slope(sites, band=(1.0, 1000.0))
 7steep = slopes[
 8    slopes["slope_xy"].abs().gt(30.0)
 9    | slopes["slope_yx"].abs().gt(30.0)
10]
11
12print(steep)

Phase-slope diagnostics are useful for finding:

  • phase curves that change abruptly within the inversion band;

  • stations with unstable tensor estimates;

  • band selections that cross noisy frequency ranges;

  • components that may require visual inspection before inversion.

The function does not unwrap phase. If phase wraps are important for your dataset, inspect the curves directly before treating the slope as a physical trend.

Tipper Magnitude#

tipper_magnitude() computes the magnitude of the tipper vector:

\[\|\mathbf{T}\| = \sqrt{|T_x|^2 + |T_y|^2}.\]

By default it returns summary statistics. Set per_freq=True to return one row per frequency.

 1from pycsamt.site import Sites
 2from pycsamt.site.compute import tipper_magnitude
 3
 4sites = Sites.from_path("data/edi")
 5
 6summary = tipper_magnitude(sites)
 7long_table = tipper_magnitude(sites, per_freq=True)
 8
 9print(summary[["station", "mean", "median", "max"]])
10print(long_table.head())

Summary mode returns:

Column

Meaning

mean

Average tipper magnitude over available frequencies.

median

Median tipper magnitude, less sensitive to isolated spikes.

max

Maximum tipper magnitude, useful for finding extreme response bands.

Per-frequency mode returns station, freq, and mag columns for collection inputs. For one site, it returns a dictionary containing freq and mag arrays.

APIFrame Output#

For collection inputs, pass api=True when the result should carry pyCSAMT API metadata in addition to tabular values.

 1from pycsamt.site import Sites
 2from pycsamt.site.compute import (
 3    phase_slope,
 4    res_at_freq,
 5    strike_estimate,
 6    tipper_magnitude,
 7)
 8
 9sites = Sites.from_path("data/edi")
10
11strike = strike_estimate(sites, api=True)
12rho = res_at_freq(sites, 100.0, api=True)
13slopes = phase_slope(sites, (1.0, 1000.0), api=True)
14tipper = tipper_magnitude(sites, api=True)
15
16print(strike.kind)
17print(rho.kind)
18print(slopes.kind)
19print(tipper.kind)

This is useful when diagnostics are emitted by CLI commands, agents, or pipelines that preserve result provenance.

Quality-Control Workflow#

The following example combines selection, editing, and computed diagnostics.

 1from pycsamt.site import Sites
 2from pycsamt.site.compute import (
 3    phase_slope,
 4    res_at_freq,
 5    strike_estimate,
 6    tipper_magnitude,
 7)
 8from pycsamt.site.edit import fill_missing
 9from pycsamt.site.selection import by_freq, keep_finite_z
10
11sites = Sites.from_path("data/edi")
12sites = keep_finite_z(sites)
13sites = by_freq(sites, fmin=1.0, fmax=1000.0)
14
15prepared = [
16    fill_missing(site, how="zero", components=("Z",), inplace=False)
17    for site in sites
18]
19
20strike = strike_estimate(prepared)
21rho100 = res_at_freq(prepared, 100.0, how="nearest")
22slopes = phase_slope(prepared, band=(1.0, 1000.0))
23tipper = tipper_magnitude(prepared)
24
25qc = (
26    strike
27    .merge(rho100, on="station", how="outer")
28    .merge(slopes, on="station", how="outer")
29    .merge(tipper, on="station", how="outer")
30)
31
32print(qc)

Common Mistakes#

Using diagnostics as final interpretation

These functions are quick screening tools. Confirm important decisions with maps, pseudo-sections, tensor plots, and inversion sensitivity tests.

Mixing frequency axes without checking f_used

In nearest mode, different stations may use different native frequencies. Always inspect f_used before comparing values.

Ignoring missing arrays

Missing Z or tipper data are reported as NaN rather than raising hard failures. Filter or fill intentionally before trusting the table.

Over-reading phase slope

A large phase slope may indicate structure, noise, phase wrapping, or a bad frequency band. Use it as a prompt for inspection, not as a standalone classifier.

Next Pages#

Continue with:

  • Site Selection for selecting stations before diagnostics;

  • Site Editing for preparing tensors, frequencies, names, and coordinates;

  • Export And Reporting for writing diagnostic summaries into survey deliverables;

  • ../theory/impedance_tensor for the physical meaning of impedance, apparent resistivity, phase, and tensor components.