Static-Shift Correction#
Static shift is a station-dependent, frequency-independent offset in
apparent resistivity. In CSAMT, AMT, and MT data it is usually caused
by shallow galvanic distortion: the curve keeps nearly the same shape,
but it is lifted or lowered by a multiplier. Because apparent
resistivity is proportional to |Z|**2, pyCSAMT estimates the shift in
log10(rho) and applies the square-root correction to the impedance
tensor Z.
The static-shift tools in pycsamt.emtools cover four related jobs:
Job |
Main tools |
Use when |
|---|---|---|
Estimate correction factors |
|
You want a per-station table of static-shift amplitudes before modifying data. |
Apply factors |
|
You have accepted a factor table and want corrected |
Check the correction visually |
|
You need before/after diagnostics rather than only a table. |
Separate static shift from near-surface effects |
|
The distortion may be frequency-dependent and therefore not removable by a conventional static-shift multiplier. |
The examples below use two-level imports from pycsamt.emtools because
the static-shift functions are exported as public user-facing symbols.
Load The Survey#
Start from the canonical loader. It accepts a directory, a glob-like
collection, or an existing Sites object and skips files without valid
impedance data.
1from pathlib import Path
2
3from pycsamt.emtools import ensure_sites
4
5edi_dir = Path("data/AMT/WILLY_DATA/L18PLT")
6sites = ensure_sites(edi_dir, recursive=True, verbose=0)
For static-shift work, station order matters. The AMA method estimates a
local spatial trend from neighbouring stations, so choose sort_by to
match the survey line direction:
|
Meaning |
|---|---|
|
Order stations by longitude. This is useful for an east-west line. |
|
Order stations by latitude. This is useful for a north-south line. |
|
Order stations lexically by station name. Use this only when the station names encode the real line order. |
Estimate AMA Factors#
estimate_ss_ama is the main estimator. AMA means adaptive
moving-average: for each station, pyCSAMT builds a weighted local
reference from neighbouring stations, compares the target station with
that reference in log10(rho_det), and reduces the per-frequency
residuals to one static-shift value.
1from pycsamt.emtools import estimate_ss_ama
2
3factors = estimate_ss_ama(
4 sites,
5 sort_by="lat",
6 half_window=3,
7 weights="tri",
8 pband=(1e-4, 10.0),
9 max_skew=45.0,
10 robust_freq="median",
11 robust_overall="median",
12)
13
14print(factors[["station", "delta_log10_rho", "fac_rho", "fac_z", "n_used"]])
station delta_log10_rho fac_rho fac_z n_used
0 18-001A 1.315686 0.048341 0.219865 22
1 18-002U 0.199660 0.631451 0.794639 32
2 18-003A -0.429163 2.686355 1.639010 38
3 18-004A 0.829915 0.147940 0.384630 34
4 18-005U 1.204691 0.062418 0.249836 31
5 18-006A -0.134703 1.363652 1.167755 35
6 18-007U 1.090404 0.081207 0.284969 31
7 18-008U 0.842759 0.143629 0.378984 30
8 18-009A 0.311615 0.487960 0.698542 43
9 18-010U 0.531531 0.294082 0.542293 40
10 18-011A 0.429333 0.372107 0.610005 39
11 18-012A 0.716283 0.192184 0.438388 32
12 18-013U 0.971147 0.106869 0.326909 37
13 18-014A 0.077767 0.836051 0.914358 35
14 18-015U -0.008743 1.020336 1.010117 41
15 18-016A 0.553432 0.279620 0.528791 39
16 18-017U 0.243962 0.570214 0.755125 40
17 18-018A 0.570181 0.269041 0.518692 8
18 18-019U -0.231489 1.704076 1.305403 11
19 18-020A 0.611525 0.244610 0.494581 25
20 18-021B 0.714440 0.193001 0.439319 22
21 18-021U 1.176125 0.066661 0.258189 23
22 18-022U -0.415305 2.601986 1.613067 18
23 18-022V 0.211366 0.614659 0.784002 21
24 18-023A -0.261413 1.825629 1.351158 16
25 18-023V 0.342848 0.454101 0.673870 11
26 18-024U 0.003641 0.991652 0.995817 9
27 18-025A -1.909562 81.201150 9.011168 17
The returned table has one row per station with these columns:
Column |
Interpretation |
|---|---|
|
Station identifier used to match the factor back to the survey. |
|
Estimated vertical offset of the station relative to the local spatial trend. Positive means the station’s apparent resistivity is above the trend. |
|
Resistivity correction factor, computed as
|
|
Impedance correction factor, computed as
|
|
Number of finite frequency samples that contributed to the station estimate after period-band and skew filtering. |
The phase-tensor skew filter is important. The default threshold is conservative, so a structurally complex survey may return a sparse table unless you choose a threshold appropriate for the line. A quick audit is usually better than accepting the first number silently:
1strict_factors = estimate_ss_ama(
2 sites,
3 sort_by="lat",
4 half_window=3,
5 max_skew=6.0,
6)
7
8survey_factors = estimate_ss_ama(
9 sites,
10 sort_by="lat",
11 half_window=3,
12 max_skew=45.0,
13)
14
15print("strict stations:", len(strict_factors))
16print("survey stations:", len(survey_factors))
17print("median samples:", survey_factors["n_used"].median())
strict stations: 26
survey stations: 28
median samples: 31.5
If n_used is low for many stations, the result may be technically
valid but weakly constrained. In that case, inspect the skew page,
restrict pband to the useful period range, or compare several
estimators before applying the correction.
Read Factor Signs Correctly#
A positive delta_log10_rho means the raw apparent resistivity is high
relative to the local trend, so the correction factor is less than one.
A negative delta_log10_rho means the raw apparent resistivity is low,
so the correction factor is greater than one.
1view = factors.assign(
2 direction=factors["delta_log10_rho"].map(
3 lambda value: "lower Z" if value > 0 else "raise Z"
4 )
5)
6
7print(
8 view[
9 ["station", "delta_log10_rho", "fac_z", "n_used", "direction"]
10 ].sort_values("delta_log10_rho")
11)
station delta_log10_rho fac_z n_used direction
27 18-025A -1.909562 9.011168 17 raise Z
2 18-003A -0.429163 1.639010 38 raise Z
22 18-022U -0.415305 1.613067 18 raise Z
24 18-023A -0.261413 1.351158 16 raise Z
18 18-019U -0.231489 1.305403 11 raise Z
5 18-006A -0.134703 1.167755 35 raise Z
14 18-015U -0.008743 1.010117 41 raise Z
26 18-024U 0.003641 0.995817 9 lower Z
13 18-014A 0.077767 0.914358 35 lower Z
1 18-002U 0.199660 0.794639 32 lower Z
23 18-022V 0.211366 0.784002 21 lower Z
16 18-017U 0.243962 0.755125 40 lower Z
8 18-009A 0.311615 0.698542 43 lower Z
25 18-023V 0.342848 0.673870 11 lower Z
10 18-011A 0.429333 0.610005 39 lower Z
9 18-010U 0.531531 0.542293 40 lower Z
15 18-016A 0.553432 0.528791 39 lower Z
17 18-018A 0.570181 0.518692 8 lower Z
19 18-020A 0.611525 0.494581 25 lower Z
20 18-021B 0.714440 0.439319 22 lower Z
11 18-012A 0.716283 0.438388 32 lower Z
3 18-004A 0.829915 0.384630 34 lower Z
7 18-008U 0.842759 0.378984 30 lower Z
12 18-013U 0.971147 0.326909 37 lower Z
6 18-007U 1.090404 0.284969 31 lower Z
21 18-021U 1.176125 0.258189 23 lower Z
4 18-005U 1.204691 0.249836 31 lower Z
0 18-001A 1.315686 0.219865 22 lower Z
fac_z must be finite and positive. apply_ss_factors guards the
impedance tensor against non-finite, zero, and negative factors by
leaving those stations unchanged, but your processing report should still
flag such rows because they indicate a failed or unreliable estimate.
1import numpy as np
2
3good = np.isfinite(factors["fac_z"]) & (factors["fac_z"] > 0.0)
4rejected = factors.loc[~good, ["station", "fac_z", "n_used"]]
5
6if not rejected.empty:
7 print("Rows that will not be applied:")
8 print(rejected)
Apply Factors#
Use apply_ss_factors when you want an explicit estimate-then-apply
workflow. This is the preferred pattern for production processing
because the table can be saved, reviewed, plotted, and reproduced.
1from pycsamt.emtools import apply_ss_factors
2
3corrected = apply_ss_factors(
4 sites,
5 factors,
6 key="fac_z",
7 inplace=False,
8)
Keep inplace=False while developing a processing flow. It returns a
corrected copy and leaves the original Sites object available for
before/after plots. Use inplace=True only when you deliberately want
to mutate the object already in memory.
For exploratory notebooks, correct_ss_ama combines estimation and
application in one call:
1from pycsamt.emtools import correct_ss_ama
2
3corrected = correct_ss_ama(
4 sites,
5 sort_by="lat",
6 half_window=3,
7 max_skew=45.0,
8 inplace=False,
9)
This convenience call is useful, but it does not show the factor table by
itself. When documenting a survey, keep the explicit factors table
as part of the processing record.
Compare Estimators#
Static-shift estimation is not a single universal recipe. The four estimators use different references:
Estimator |
Reference model |
Practical reading |
|---|---|---|
|
Weighted local median/mean from neighbouring stations. |
Good default for line data when station order is meaningful. |
|
Local polynomial regression across neighbouring stations. |
Useful when the background level changes smoothly along the line. |
|
Neighbour trend weighted by both distance and value similarity. |
Can protect sharp local contrasts, but may disagree at outliers. |
|
Global frequency-wise median reference. |
Useful as a broad check, less local than AMA or LOESS. |
Run several estimators before correcting a valuable line. Agreement is a stronger signal than a single factor table.
1from functools import reduce
2
3from pycsamt.emtools import (
4 estimate_ss_bilateral,
5 estimate_ss_loess,
6 estimate_ss_refmedian,
7)
8
9tables = {
10 "ama": estimate_ss_ama(
11 sites,
12 sort_by="lat",
13 half_window=3,
14 max_skew=45.0,
15 ),
16 "loess": estimate_ss_loess(
17 sites,
18 half_window=3,
19 max_skew=45.0,
20 ),
21 "bilateral": estimate_ss_bilateral(
22 sites,
23 half_window=4,
24 max_skew=45.0,
25 ),
26 "refmedian": estimate_ss_refmedian(
27 sites,
28 max_skew=45.0,
29 ),
30}
31
32compact = []
33for name, table in tables.items():
34 compact.append(
35 table[["station", "delta_log10_rho"]].rename(
36 columns={"delta_log10_rho": name}
37 )
38 )
39
40comparison = reduce(
41 lambda left, right: left.merge(right, on="station", how="inner"),
42 compact,
43)
44
45print(comparison)
46print(comparison.drop(columns="station").corr().round(2))
station ama loess bilateral refmedian
0 18-001A 1.332021 1.159768 0.404790 1.147046
1 18-002U 0.204137 0.027926 -0.024785 0.344798
2 18-003A -0.423116 -0.228044 -0.056300 -0.013361
3 18-004A 0.845136 0.544652 0.027856 0.593522
4 18-005U 1.212859 1.208384 0.350171 1.008823
5 18-006A -0.023122 -0.111229 -0.094693 0.382904
6 18-007U 1.108909 1.196616 0.022501 0.870229
7 18-008U 0.826022 0.576419 0.101546 0.799377
8 18-009A 0.308615 0.164553 0.064798 0.466872
9 18-010U 0.387382 0.575710 0.001690 0.436586
10 18-011A 0.316431 0.272670 -0.071976 0.329512
11 18-012A 0.709257 0.750474 0.524762 0.782521
12 18-013U 0.987698 0.823783 0.408062 0.983713
13 18-014A 0.060567 0.085557 -0.001016 0.328062
14 18-015U -0.027194 0.148288 0.011463 0.437894
15 18-016A 0.548499 0.682097 0.222652 0.387643
16 18-017U 0.207717 -0.072589 0.175048 -0.020784
17 18-018A 0.559294 0.622097 0.204397 0.442325
18 18-019U -0.787610 -0.750317 -0.310667 -0.562525
19 18-020A 0.724057 0.577375 0.219258 0.793404
20 18-021B 0.454139 0.194675 0.176764 0.632271
21 18-021U 1.304396 1.190724 0.551496 1.188714
22 18-022U -0.404453 -0.554035 -0.135679 -0.246598
23 18-022V 0.276068 0.284507 0.184105 0.341684
24 18-023A -0.261413 -0.354546 -0.189863 -0.345754
25 18-023V 0.258398 -0.256520 0.058670 -0.092407
26 18-024U 0.015435 0.353124 -0.039578 0.014523
27 18-025A -1.910982 -2.802634 1.019139 -1.579179
ama loess bilateral refmedian
ama 1.00 0.96 0.08 0.95
loess 0.96 1.00 -0.08 0.95
bilateral 0.08 -0.08 1.00 0.05
refmedian 0.95 0.95 0.05 1.00
Look for sign disagreements, not only large absolute differences. If AMA says a station should be lowered while bilateral filtering says it should be raised, that station deserves a manual plot before correction.
QC Plots In One Call#
The ss_qc_* wrappers estimate, apply, and plot in one call. They are
good first-look tools when you do not need to reuse the corrected object.
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import (
4 ss_qc_profile,
5 ss_qc_psection,
6 ss_qc_station_curves,
7)
8
9ss_qc_psection(
10 sites,
11 method="ama",
12 sort_by="lat",
13 half_window=3,
14 max_skew=45.0,
15 pband=(1e-4, 10.0),
16)
17
18ss_qc_profile(
19 sites,
20 method="ama",
21 sort_by="lat",
22 half_window=3,
23 max_skew=45.0,
24 robust="median",
25)
26
27ss_qc_station_curves(
28 sites,
29 method="ama",
30 station="18-016A",
31 sort_by="lat",
32 half_window=3,
33 max_skew=45.0,
34)
35
36for i, num in enumerate(plt.get_fignums(), start=1):
37 fig = plt.figure(num)
38 fig.savefig(f"ss_qc_plots_{i:02d}.png", dpi=200, bbox_inches="tight")
39 plt.close(fig)
ss_qc_psection shows the pointwise
log10(rho_after) - log10(rho_before) correction as a pseudosection.
For a true static-shift correction, each station tends to show a nearly
vertical band because the multiplier is constant with frequency.
ss_qc_profile compresses that difference to one value per station.
Use it to find stations with unusually large corrections.
ss_qc_station_curves overlays the before and after curves for one
station. It is the quickest way to confirm that the correction shifted
the level without changing the curve shape.
Before/After Plots From Existing Sites#
When you already have a corrected Sites object, call the lower-level
plotters directly. These functions do not re-estimate the correction;
they compare the two data sets you pass in.
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import (
4 plot_ss_delta_profile,
5 plot_ss_delta_psection,
6 plot_ss_station_curves,
7)
8
9plot_ss_delta_psection(
10 sites,
11 corrected,
12 axis_y="logperiod",
13 pband=(1e-4, 10.0),
14)
15
16plot_ss_delta_profile(
17 sites,
18 corrected,
19 robust="median",
20 pband=(1e-4, 10.0),
21)
22
23plot_ss_station_curves(
24 sites,
25 corrected,
26 station="18-016A",
27 pband=(1e-4, 10.0),
28)
29
30for i, num in enumerate(plt.get_fignums(), start=1):
31 fig = plt.figure(num)
32 fig.savefig(f"ss_before_after_{i:02d}.png", dpi=200, bbox_inches="tight")
33 plt.close(fig)
These calls are useful in reports because they make the provenance clear:
sites is the uncorrected input and corrected is the accepted
corrected output.
Publication Comparison Figures#
ss_comparison_psection is the high-level publication helper. It
estimates the correction, applies it, builds aligned before/after
log10(rho_det) arrays internally, and renders a shared-scale
pseudosection comparison.
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import ss_comparison_psection
4
5fig = ss_comparison_psection(
6 sites,
7 method="ama",
8 sort_by="lat",
9 half_window=3,
10 max_skew=45.0,
11 show_delta=True,
12 suptitle="Static-shift correction: AMA",
13)
14fig.savefig("ss_comparison_psection_ama.png", dpi=200, bbox_inches="tight")
15plt.close(fig)
Use show_delta=True when the figure must show both the corrected
resistivity field and the actual correction amplitude. The before and
after panels share colour limits, so station-level offsets remain
visible instead of being hidden by automatic rescaling.
The lower-level plot_ss_comparison_psection, plot_ss_1d_curves,
and plot_ss_summary functions accept precomputed arrays:
logRho_before with shape (n_stations, n_frequencies),
logRho_after with the same shape, and freqs in hertz. Use those
array-level functions when you have already built a custom resistivity
matrix outside the normal Sites workflow.
Radar View#
plot_ss_radar shows the xy and yx apparent-resistivity
components for one station on a polar grid. The angle around the circle
represents period or frequency, and the radius represents resistivity or
log10(resistivity).
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import plot_ss_radar
4
5plot_ss_radar(
6 sites,
7 station="18-016A",
8 rotate="none",
9 radial="log10rho",
10 theta_axis="logperiod",
11 fill_between=True,
12)
13
14plot_ss_radar(
15 sites,
16 station="18-016A",
17 rotate="pt",
18 rotate_stat="median",
19 radial="log10rho",
20)
21
22for i, num in enumerate(plt.get_fignums(), start=1):
23 fig = plt.figure(num)
24 fig.savefig(f"ss_radar_18-016A_{i:02d}.png", dpi=200, bbox_inches="tight")
25 plt.close(fig)
Use rotate="none" for a raw component view, rotate="deg" with
rotate_deg=... for a fixed rotation, and rotate="pt" for
phase-tensor based rotation. The radar plot is diagnostic: it helps you
see directional imbalance and frequency structure, but it is not itself a
static-shift estimator.
Near-Surface Versus Static Shift#
Conventional static-shift correction assumes the distortion is a constant multiplier. Some shallow effects are frequency-dependent: the high frequency part of the curve may become much noisier or steeper than the low frequency part. A single static multiplier cannot fix that behavior.
detect_near_surface classifies each station using three diagnostics:
Diagnostic |
Meaning |
|---|---|
|
Ratio of high-frequency residual spread to low-frequency residual
spread. Values above |
|
Difference between high-frequency and low-frequency log-log slopes. Larger values indicate a change in curve shape. |
|
AMA-like static-shift residual. Values above |
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import detect_near_surface, plot_ns_detection
4
5ns = detect_near_surface(
6 sites,
7 f_split=1.0,
8 ns_threshold=2.0,
9 ss_threshold=0.1,
10 sort_by="lat",
11 half_window=3,
12 max_skew=45.0,
13)
14
15print(ns[[
16 "station",
17 "ns_index",
18 "gradient_delta",
19 "ss_delta_log10",
20 "distortion_type",
21]])
22print(ns["distortion_type"].value_counts())
23
24plot_ns_detection(
25 sites,
26 f_split=1.0,
27 ns_threshold=2.0,
28 ss_threshold=0.1,
29 sort_by="lat",
30 half_window=3,
31 max_skew=45.0,
32 show_ss=True,
33)
34plt.gcf().savefig("ss_near_surface_detection.png", dpi=200, bbox_inches="tight")
35plt.close(plt.gcf())
station ns_index gradient_delta ss_delta_log10 distortion_type
0 18-001A NaN NaN 1.332021 static
1 18-002U NaN NaN 0.204137 static
2 18-003A NaN NaN -0.423116 static
3 18-004A NaN NaN 0.845136 static
4 18-005U NaN NaN 1.212859 static
5 18-006A NaN NaN -0.023122 clean
6 18-007U NaN NaN 1.108909 static
7 18-008U NaN NaN 0.826022 static
8 18-009A NaN NaN 0.308615 static
9 18-010U NaN NaN 0.387382 static
10 18-011A NaN NaN 0.316431 static
11 18-012A NaN NaN 0.709257 static
12 18-013U NaN NaN 0.987698 static
13 18-014A NaN NaN 0.060567 clean
14 18-015U NaN NaN -0.027194 clean
15 18-016A NaN NaN 0.548499 static
16 18-017U NaN NaN 0.207717 static
17 18-018A NaN NaN 0.559294 static
18 18-019U NaN NaN -0.787610 static
19 18-020A NaN NaN 0.724057 static
20 18-021B NaN NaN 0.454139 static
21 18-021U NaN NaN 1.304396 static
22 18-022U NaN NaN -0.404453 static
23 18-022V NaN NaN 0.276068 static
24 18-023A NaN NaN -0.261413 static
25 18-023V NaN NaN 0.258398 static
26 18-024U NaN NaN 0.015435 clean
27 18-025A NaN NaN -1.910982 static
distortion_type
static 24
clean 4
Name: count, dtype: int64
The distortion_type column has four possible values:
Type |
Meaning |
|---|---|
|
No strong near-surface index and no strong static-shift residual. |
|
Static-shift residual is large, but the high-frequency residual spread is not. A conventional static-shift correction may help. |
|
High-frequency residual spread is large, but static residual is not. Treat this as frequency-dependent distortion, not as a simple multiplier. |
|
Both diagnostics are large. Static correction may remove the constant part, but the remaining frequency-dependent effect needs separate geological and inversion judgment. |
Recommended Processing Pattern#
For a survey report, keep the static-shift workflow explicit:
1import matplotlib.pyplot as plt
2
3from pathlib import Path
4
5from pycsamt.emtools import (
6 apply_ss_factors,
7 detect_near_surface,
8 ensure_sites,
9 estimate_ss_ama,
10 ss_comparison_psection,
11)
12
13edi_dir = Path("data/AMT/WILLY_DATA/L18PLT")
14sites = ensure_sites(edi_dir, recursive=True)
15
16factors = estimate_ss_ama(
17 sites,
18 sort_by="lat",
19 half_window=3,
20 pband=(1e-4, 10.0),
21 max_skew=45.0,
22)
23
24factors.to_csv("static_shift_factors.csv", index=False)
25
26ns = detect_near_surface(
27 sites,
28 f_split=1.0,
29 sort_by="lat",
30 half_window=3,
31 max_skew=45.0,
32)
33ns.to_csv("near_surface_diagnostics.csv", index=False)
34
35accepted = factors.merge(
36 ns[["station", "distortion_type"]],
37 on="station",
38 how="left",
39)
40
41print(
42 accepted.groupby("distortion_type", dropna=False)
43 ["station"]
44 .count()
45 .rename("n_stations")
46)
47
48corrected = apply_ss_factors(
49 sites,
50 factors,
51 key="fac_z",
52 inplace=False,
53)
54
55fig = ss_comparison_psection(
56 sites,
57 method="ama",
58 sort_by="lat",
59 half_window=3,
60 pband=(1e-4, 10.0),
61 max_skew=45.0,
62 show_delta=True,
63)
64fig.savefig("ss_recommended_comparison.png", dpi=200, bbox_inches="tight")
65plt.close(fig)
distortion_type
clean 4
static 24
Name: n_stations, dtype: int64
The important decisions are visible in that script: loader, station ordering, period band, skew threshold, accepted estimator, saved factor table, and a near-surface check before interpreting the correction.
Common Pitfalls#
sort_by does not describe the coordinate you like best; it describes
the along-line order used by the neighbourhood estimator. A north-south
line usually wants sort_by="lat". An east-west line usually wants
sort_by="lon".
max_skew is a filter, not a quality score. A strict value can leave
too few samples for a complex survey. Check n_used before trusting
the factors.
Static shift changes level, not shape. If the before/after station plot would need a frequency-dependent correction to align the curves, use the near-surface diagnostics and interpret with caution.
Do not apply an empty or nearly empty factor table just because the code returns without error. Empty estimates are a valid no-op outcome for single-station inputs or surveys where no station has usable neighbours.
Prefer fac_z for impedance correction. fac_rho is the
resistivity-domain factor and is useful for interpretation, but
apply_ss_factors scales Z.
Worked Example#
The gallery example uses the L18PLT survey and walks through raw curve spread, AMA estimation, exact factor application, estimator comparison, QC wrappers, radar plots, and near-surface classification.
Open the rendered gallery page here: Static-shift estimation, correction, and QC (pycsamt.emtools.ss).