Phase Tensor And Impedance Tensor Tools#
The tensor page is one of the central pages in the pyCSAMT user guide
because it connects interpretation and preprocessing. The same
impedance tensor Z supports apparent resistivity, phase, skew,
geoelectric strike, static-shift checks, dimensionality checks, and 2-D
rotation decisions. The pycsamt.emtools tensor tools expose two
main workflows:
Workflow |
Main tools |
Purpose |
|---|---|---|
Phase-tensor interpretation |
|
Diagnose dimensionality, strike, skew, ellipticity, and spatial coherence without being dominated by static-shift amplitude. |
Impedance-tensor editing |
|
Apply deliberate preprocessing operations to |
The phase tensor follows the Caldwell et al. style decomposition. For
each frequency, pyCSAMT splits the impedance tensor into real and
imaginary parts, then computes the phase tensor
Phi = real(Z)^-1 imag(Z) and its invariants.
The examples in this guide use public two-level imports from
pycsamt.emtools. One name needs special attention:
pycsamt.emtools.rotate_to_strike belongs to the strike module, while
the tensor module’s own rotation-to-strike helper is exported as
rotate_z_to_strike.
Load Data#
Start with the canonical loader. It returns a Sites object and skips
stations 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(
7 edi_dir,
8 recursive=True,
9 on_dup="replace",
10 strict=False,
11 verbose=0,
12)
Keep sites as the unmodified reference object. Tensor editing
functions can return corrected copies when inplace=False, so it is
easy to compare original and edited data.
Build The Phase-Tensor Table#
build_phase_tensor_table is the foundation for the plotting tools.
It returns one row per station and frequency.
1from pycsamt.emtools import build_phase_tensor_table
2
3pt = build_phase_tensor_table(
4 sites,
5 recursive=False,
6)
7
8print(pt.head())
9print(pt[["station", "freq", "period", "theta", "skew", "ellipt"]])
station freq period ... beta skew ellipt
0 18-001A 10400.0 0.000096 ... -56.700714 -56.700714 0.194909
1 18-001A 8707.0 0.000115 ... -54.693184 -54.693184 0.210163
2 18-001A 7289.0 0.000137 ... -51.452210 -51.452210 0.213894
3 18-001A 6102.0 0.000164 ... -61.983725 -61.983725 0.212655
4 18-001A 5108.0 0.000196 ... -60.874439 -60.874439 0.425503
[5 rows x 10 columns]
station freq period theta skew ellipt
0 18-001A 10400.000 0.000096 120.687698 -56.700714 0.194909
1 18-001A 8707.000 0.000115 123.342495 -54.693184 0.210163
2 18-001A 7289.000 0.000137 126.743525 -51.452210 0.213894
3 18-001A 6102.000 0.000164 116.947420 -61.983725 0.212655
4 18-001A 5108.000 0.000196 126.074830 -60.874439 0.425503
... ... ... ... ... ... ...
1479 18-025A 2.052 0.487329 -34.147397 -54.498901 0.559680
1480 18-025A 1.718 0.582072 4.537986 -41.477337 0.196258
1481 18-025A 1.438 0.695410 40.621212 -32.840732 0.216635
1482 18-025A 1.204 0.830565 -125.802958 -33.743310 0.527659
1483 18-025A 1.008 0.992063 -107.576221 -21.335294 0.931748
[1484 rows x 6 columns]
The table contains these core columns:
Column |
Meaning |
|---|---|
|
Station identifier. |
|
Frequency in hertz. |
|
Period in seconds, computed as |
|
Phase-tensor principal values. They control ellipse major and minor axes. |
|
Principal-axis angle, interpreted as phase-tensor strike. It is axial, so directions separated by 180 degrees are equivalent. |
|
Phase-tensor coordinate angle. |
|
Skew angle. |
|
Ellipticity, computed from the principal values. Values near zero are closer to circular; larger values indicate stronger 2-D anisotropy of the phase tensor. |
The table is also the best place to audit a survey numerically before plotting:
1summary = pt.groupby("station").agg(
2 n=("freq", "count"),
3 median_abs_skew=("skew", lambda values: values.abs().median()),
4 median_ellipt=("ellipt", "median"),
5 theta_iqr=("theta", lambda values: values.quantile(0.75) - values.quantile(0.25)),
6)
7
8print(summary.sort_values("median_abs_skew", ascending=False))
n median_abs_skew median_ellipt theta_iqr
station
18-023A 53 67.022970 0.776955 332.406202
18-022V 53 66.787861 0.734789 33.283984
18-018A 53 66.547818 0.747821 19.123004
18-022U 53 65.349813 0.765743 330.236042
18-024U 53 63.853268 0.696343 69.080189
18-019U 53 61.945234 0.610832 44.257571
18-023V 53 59.306877 0.568990 39.420996
18-025A 53 55.354566 0.666770 174.077409
18-021U 53 55.017266 0.975587 5.728713
18-021B 53 52.388774 0.917700 65.673584
18-001A 53 50.326802 0.423331 15.689190
18-020A 53 45.332864 0.890491 7.151592
18-008U 53 40.862593 0.459678 30.040274
18-005U 53 36.404849 0.598509 16.449050
18-002U 53 36.059416 0.535202 16.329767
18-014A 53 35.458520 0.735137 273.677122
18-007U 53 34.174772 0.492508 29.303937
18-012A 53 32.307319 0.710183 37.311060
18-003A 53 31.245824 0.674498 23.567376
18-004A 53 31.005169 0.683143 18.628263
18-006A 53 30.359830 0.574245 24.700266
18-013U 53 29.705384 0.591064 41.602228
18-011A 53 27.534040 0.692296 261.601465
18-010U 53 26.006304 0.676668 257.925336
18-009A 53 25.288856 0.477635 19.936171
18-016A 53 23.525350 0.732322 302.853539
18-017U 53 22.912833 0.660533 269.958366
18-015U 53 22.459269 0.721964 316.822496
Large median_abs_skew means the station is not behaving like a clean
1-D or 2-D response in that period range. Large theta_iqr means
phase-tensor strike changes strongly with frequency, so one rotation
angle may be a poor summary.
Filter By Period#
Most tensor interpretation should be tied to a period band. A shallow band and a deeper band can show different strike, skew, or ellipticity.
1period_band = (0.001, 10.0)
2band_pt = pt[
3 (pt["period"] >= period_band[0])
4 & (pt["period"] <= period_band[1])
5]
6
7print("rows in band:", len(band_pt))
8print("stations in band:", band_pt["station"].nunique())
9print("median |skew|:", band_pt["skew"].abs().median())
rows in band: 1092
stations in band: 28
median |skew|: 39.133827836405146
When you report a tensor result, always report the period band. A map at
period=1.0 second and a summary over 0.001 to 10.0 seconds
answer different questions.
Read Dimensionality From Skew And Ellipticity#
A simple phase-tensor dimensionality rule uses skew and ellipticity:
Class |
Rule |
Interpretation |
|---|---|---|
1-D |
|
Low skew and nearly circular phase tensor. |
2-D |
|
Low skew but elongated phase tensor. |
3-D |
|
High skew; strike and 2-D rotation should be treated with caution. |
1import numpy as np
2
3skew_threshold = 3.0
4ellipt_threshold = 0.2
5
6work = band_pt.copy()
7abs_skew = work["skew"].abs()
8abs_ellipt = work["ellipt"].abs()
9
10work["dimensionality"] = np.select(
11 [
12 (abs_skew <= skew_threshold) & (abs_ellipt <= ellipt_threshold),
13 (abs_skew <= skew_threshold) & (abs_ellipt > ellipt_threshold),
14 ],
15 ["1D", "2D"],
16 default="3D",
17)
18
19print(work["dimensionality"].value_counts(normalize=True))
dimensionality
3D 0.975275
2D 0.023810
1D 0.000916
Name: proportion, dtype: float64
The default 3 degree skew threshold is strict. It is useful as a
textbook 1-D/2-D screen, but it can classify many real field samples as
3-D. That is not a failure of the function; it is a warning about the
data and the 2-D assumption.
Simple Phase-Tensor Views#
Use the simpler plots before the full ellipse plot. They make it easy to identify which invariant is causing the interpretation.
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import (
4 plot_dimensionality_grid,
5 plot_dimensionality_psection,
6 plot_ellipticity_psection,
7 plot_phase_tensor_skewmap,
8 plot_theta_vs_period,
9)
10
11plot_theta_vs_period(sites, recursive=False)
12plt.gcf().savefig("tensor_simple_views_01.png", dpi=200, bbox_inches="tight")
13plt.close()
14
15plot_phase_tensor_skewmap(sites, recursive=False, axis_y="logperiod")
16plt.gcf().savefig("tensor_simple_views_02.png", dpi=200, bbox_inches="tight")
17plt.close()
18
19plot_ellipticity_psection(sites, recursive=False)
20plt.gcf().savefig("tensor_simple_views_03.png", dpi=200, bbox_inches="tight")
21plt.close()
22
23plot_dimensionality_psection(
24 sites,
25 skew_th=3.0,
26 ellipt_th=0.2,
27 recursive=False,
28)
29plt.gcf().savefig("tensor_simple_views_04.png", dpi=200, bbox_inches="tight")
30plt.close()
31
32plot_dimensionality_grid(
33 sites,
34 skew_th=3.0,
35 ellipt_th=0.2,
36 recursive=False,
37)
38plt.gcf().savefig("tensor_simple_views_05.png", dpi=200, bbox_inches="tight")
39plt.close()
plot_theta_vs_period is a scatter view of strike angle by period.
It is quick, but it puts an axial angle on a linear y-axis. Treat jumps
near the wrap boundary carefully.
plot_phase_tensor_skewmap and plot_ellipticity_psection show
station-by-period heatmaps. They are good for finding period bands that
are consistently low-skew or strongly elongated.
plot_dimensionality_psection and plot_dimensionality_grid apply
the simple skew/ellipticity classification to every station-period cell.
Phase-Tensor Ellipse Pseudosection#
plot_phase_tensor_psection is the main phase-tensor figure. Each
cell is an ellipse:
Visual element |
Meaning |
|---|---|
Major axis |
|
Minor axis |
|
Orientation |
|
Fill color |
Controlled by |
Thick border |
Optional marker for cells where |
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import plot_phase_tensor_psection
4
5plot_phase_tensor_psection(
6 sites,
7 stations=None,
8 period_range=(0.001, 10.0),
9 axis_y="logperiod",
10 period_up=True,
11 c_by="skew",
12 skew_threshold=3.0,
13 mark_3d=True,
14 normalise_by="cell",
15 recursive=False,
16)
17plt.gcf().savefig(
18 "tensor_phase_tensor_psection.png",
19 dpi=200,
20 bbox_inches="tight",
21)
22plt.close()
Useful c_by values include "skew", "beta", "theta",
"ellipt", "s1", "s2", "|skew|", "phi_mean",
"phi_max", and "phi_min".
Use normalise_by="cell" for most survey plots because it scales the
ellipses to the local plotting grid. Use normalise_by="unity" when
you want the 45-degree, 1-D reference to have an explicit visual meaning.
Use normalise_by="abs" only when absolute ellipse sizes in data units
are intentional.
Strike As A Director Field#
theta is axial. A director field is often easier to interpret than
a scatter plot because the glyph has no arrow head and therefore
respects the 180-degree ambiguity.
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import plot_strike_director_field
4
5plot_strike_director_field(
6 sites,
7 color_by="skew",
8 length_by="ellipt",
9 skew_max=6.0,
10 streamlines=True,
11 period_subsample=40,
12 recursive=False,
13)
14plt.gcf().savefig(
15 "tensor_strike_director_field.png",
16 dpi=200,
17 bbox_inches="tight",
18)
19plt.close()
Interpret long, aligned, low-skew directors as a more coherent 2-D strike signal. Short directors mean the phase tensor is close to circular and strike is poorly defined. High-skew directors mean the response is more 3-D or distorted, even if the direction looks visually organized.
Rose And Stability Plots#
Rose plots are useful for summarizing axial direction. Stability plots are useful for seeing whether that summary hides period dependence.
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import (
4 plot_phase_tensor_rose,
5 plot_theta_rose_grid,
6 plot_theta_stability_stripe,
7)
8
9plot_phase_tensor_rose(
10 sites,
11 band=(0.001, 10.0),
12 bins=36,
13 recursive=False,
14)
15plt.gcf().savefig("tensor_phase_tensor_rose.png", dpi=200, bbox_inches="tight")
16plt.close()
17
18plot_theta_rose_grid(
19 sites,
20 n_bands=6,
21 bins=24,
22 recursive=False,
23)
24plt.gcf().savefig("tensor_theta_rose_grid.png", dpi=200, bbox_inches="tight")
25plt.close()
26
27plot_theta_stability_stripe(
28 sites,
29 win=5,
30 recursive=False,
31)
32plt.gcf().savefig(
33 "tensor_theta_stability_stripe.png",
34 dpi=200,
35 bbox_inches="tight",
36)
37plt.close()
plot_phase_tensor_rose folds all selected theta values into one
axial histogram. plot_theta_rose_grid splits the period range into
equal log-width bands. plot_theta_stability_stripe uses hue for
theta and saturation for local stability.
If the rose grid changes direction from one period band to another, avoid selecting a single strike for the entire data set.
Skew-Ellipticity Density#
plot_skew_ellipt_density shows the joint distribution of
abs(beta) and abs(ellipt). It is a compact way to see whether
the data cluster in a 1-D, 2-D, or 3-D region.
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import plot_skew_ellipt_density
4
5plot_skew_ellipt_density(
6 sites,
7 band=(0.001, 10.0),
8 gridsize=40,
9 recursive=False,
10)
11plt.gcf().savefig(
12 "tensor_skew_ellipt_density.png",
13 dpi=200,
14 bbox_inches="tight",
15)
16plt.close()
Use this plot with the dimensionality grid. The grid tells you where problem cells occur; the density plot tells you how the full population is distributed.
Summary Figure#
plot_phase_tensor_summary combines the main ellipse pseudosection,
the dimensionality grid, and the skew-ellipticity density into one
figure.
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import plot_phase_tensor_summary
4
5fig = plot_phase_tensor_summary(
6 sites,
7 stations=None,
8 period_range=(0.001, 10.0),
9 c_by="skew",
10 skew_threshold=3.0,
11 ellipt_threshold=0.2,
12 recursive=False,
13)
14fig.savefig("tensor_phase_tensor_summary.png", dpi=200, bbox_inches="tight")
15plt.close(fig)
This is a good report figure when the audience needs the whole phase-tensor story in one place: what the ellipses look like, how much of the band is 1-D/2-D/3-D, and where the skew/ellipticity population sits.
Geographic Phase-Tensor Map#
plot_phase_tensor_map draws one phase-tensor ellipse per station at
the period nearest a requested target period. It can also overlay tipper
arrows when vertical magnetic transfer functions are present.
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import plot_phase_tensor_map
4
5plot_phase_tensor_map(
6 sites,
7 period=1.0,
8 c_by="skew",
9 show_tipper=True,
10 tipper_convention="parkinson",
11 station_labels=True,
12 recursive=False,
13)
14plt.gcf().savefig("tensor_phase_tensor_map.png", dpi=200, bbox_inches="tight")
15plt.close()
Use period as a target; each station uses its nearest available
period. If the EDI headers do not provide usable coordinates, pass an
explicit coords dictionary:
1import matplotlib.pyplot as plt
2
3coords = {
4 "18-001A": (7.312, -5.218),
5 "18-002U": (7.318, -5.211),
6 "18-003A": (7.324, -5.204),
7}
8
9plot_phase_tensor_map(
10 sites,
11 period=1.0,
12 coords=coords,
13 show_tipper=False,
14 recursive=False,
15)
16plt.gcf().savefig(
17 "tensor_phase_tensor_map_custom_coords.png",
18 dpi=200,
19 bbox_inches="tight",
20)
21plt.close()
The coordinate tuple is (lat, lon). A map with no coordinates is not
a tensor failure; it is a metadata problem. Use pseudosections and
profiles until coordinates are supplied.
Per-Station Ellipse Strips#
plot_phase_tensor_strip draws one station as an ellipse sequence
through period. plot_phase_tensor_strip_grid tiles several such
strips, optionally grouped by profile.
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import (
4 plot_phase_tensor_strip,
5 plot_phase_tensor_strip_grid,
6)
7
8plot_phase_tensor_strip(
9 sites,
10 station="18-016A",
11 period_range=(0.001, 10.0),
12 c_by="skew",
13 recursive=False,
14)
15plt.gcf().savefig("tensor_phase_tensor_strip.png", dpi=200, bbox_inches="tight")
16plt.close()
17
18groups = {
19 "L18PLT": ["18-001A", "18-002U", "18-003A", "18-004A"],
20}
21
22plot_phase_tensor_strip_grid(
23 sites,
24 groups,
25 period_range=(0.001, 10.0),
26 c_by="skew",
27 recursive=False,
28)
29plt.gcf().savefig(
30 "tensor_phase_tensor_strip_grid.png",
31 dpi=200,
32 bbox_inches="tight",
33)
34plt.close()
Use strips when a single station deserves close inspection. Use the grid when comparing stations or profile groups without the compression of a full pseudosection.
Standalone Legend#
phase_tensor_legend draws a reference ellipse that can be composed
into custom figures.
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import phase_tensor_legend
4
5phase_tensor_legend(size=1.0)
6plt.gcf().savefig("tensor_phase_tensor_legend.png", dpi=200, bbox_inches="tight")
7plt.close()
It is not a diagnostic by itself; it is a small plotting component for figures where the phase-tensor ellipse convention needs to be explained.
Impedance-Tensor Editing#
Tensor editing functions change Z. They should be treated as
processing operations, not as harmless plots. Keep inplace=False
unless you deliberately want to mutate the object in memory.
Function |
What it does |
Typical use |
|---|---|---|
|
Rotates all station impedance tensors by one fixed angle. |
Apply a known regional strike or sensor-frame correction. |
|
Rotates each station by a value from a station-angle dictionary. |
Apply station-specific rotations from an external interpretation. |
|
Estimates a tensor strike for each station and rotates by it. |
Exploratory tensor-level strike rotation. |
|
Enforces off-diagonal antisymmetry. |
Prepare data for methods that assume |
|
Corrects electric and magnetic sensor orientations. |
Fix known sensor azimuth errors. |
|
Sets outlying |
Remove isolated spikes before plotting or inversion. |
|
Balances off-diagonal magnitudes. |
Exploratory symmetry conditioning. |
|
Replaces each 2-by-2 impedance tensor by its matrix inverse. |
Advanced workflows that explicitly need admittance-like tensors. |
Fixed-Angle Rotation#
Use rotate when every station should receive the same angle.
1from pycsamt.emtools import rotate
2
3rotated_30 = rotate(
4 sites,
5 30.0,
6 inplace=False,
7 recursive=False,
8)
Use this only when the angle is justified by strike analysis, field geometry, or a documented processing convention.
Station-Specific Rotation#
Use rotate_by_map when each station needs its own angle.
1from pycsamt.emtools import rotate_by_map
2
3angle_by_station = {
4 "18-001A": 25.0,
5 "18-002A": 27.5,
6 "18-003A": 24.0,
7}
8
9rotated_stationwise = rotate_by_map(
10 sites,
11 angle_by_station,
12 inplace=False,
13 recursive=False,
14)
Stations missing from the dictionary receive a 0 degree rotation.
For production work, check that every intended station is present in the
map before applying it.
Tensor Rotation To Strike#
The tensor module’s rotation-to-strike function is exported as
rotate_z_to_strike at the pycsamt.emtools level. This avoids a
name collision with the strike module’s rotate_to_strike.
1from pycsamt.emtools import rotate_z_to_strike
2
3rotated_to_tensor_strike = rotate_z_to_strike(
4 sites,
5 method="swift",
6 inplace=False,
7 recursive=False,
8)
Use this as a tensor editing step, not as the main strike-analysis
interface. For detailed strike estimation, use the dedicated
strike.rst workflow and its station-level tables.
Antisymmetrize And Balance#
antisymmetrize enforces off-diagonal antisymmetry. balance_offdiag
balances |Zxy| and |Zyx| while preserving phase.
1from pycsamt.emtools import antisymmetrize, balance_offdiag
2
3antisymmetric = antisymmetrize(
4 sites,
5 how="rms",
6 inplace=False,
7 recursive=False,
8)
9
10balanced = balance_offdiag(
11 sites,
12 mode="avgabs",
13 inplace=False,
14 recursive=False,
15)
These operations can be useful for controlled experiments or for preparing data for algorithms with strong 2-D assumptions. They can also hide real 3-D information if applied blindly. Always keep an uncorrected copy.
Sensor Orientation#
Use orient_from_sensors when field notes show that the electric and
magnetic sensors were not aligned with the assumed axes.
1from pycsamt.emtools import orient_from_sensors
2
3oriented = orient_from_sensors(
4 sites,
5 ex=5.0,
6 ey=95.0,
7 bx=5.0,
8 by=95.0,
9 degrees=True,
10 inplace=False,
11 recursive=False,
12)
Angles are degrees by default. Pass degrees=False only when your
inputs are radians. This operation is meaningful only when the sensor
orientation metadata or field notes are trustworthy.
Sigma Clip Outliers#
sigma_clip_z flags outlying entries in the complex impedance tensor
and sets them to NaN.
1from pycsamt.emtools import sigma_clip_z
2
3clipped = sigma_clip_z(
4 sites,
5 sigma=3.0,
6 inplace=False,
7 recursive=False,
8)
After clipping, rerun coverage and QC checks. A small number of clipped entries may remove isolated spikes; a large number means the survey or threshold needs review.
Invert Tensor#
invert applies the 2-by-2 matrix inverse frequency by frequency.
1from pycsamt.emtools import invert
2
3admittance_like = invert(
4 sites,
5 inplace=False,
6 recursive=False,
7)
This is an advanced operation. Do not use it as a generic noise correction. It changes the physical quantity represented by the tensor.
Audit Edits With Phase-Tensor Tables#
A safe editing workflow compares before and after diagnostics.
1import matplotlib.pyplot as plt
2
3from pycsamt.emtools import (
4 build_phase_tensor_table,
5 plot_phase_tensor_summary,
6 rotate,
7)
8
9before = build_phase_tensor_table(sites, recursive=False)
10
11rotated = rotate(
12 sites,
13 30.0,
14 inplace=False,
15 recursive=False,
16)
17
18after = build_phase_tensor_table(rotated, recursive=False)
19
20audit = before[["station", "period", "theta", "skew", "ellipt"]].merge(
21 after[["station", "period", "theta", "skew", "ellipt"]],
22 on=["station", "period"],
23 suffixes=("_before", "_after"),
24)
25
26audit["theta_change"] = (
27 (audit["theta_after"] - audit["theta_before"] + 90.0) % 180.0
28) - 90.0
29
30print(audit[["station", "period", "theta_change", "skew_before", "skew_after"]])
31
32fig = plot_phase_tensor_summary(
33 rotated,
34 period_range=(0.001, 10.0),
35 recursive=False,
36)
37fig.savefig("tensor_audit_summary.png", dpi=200, bbox_inches="tight")
38plt.close(fig)
station period theta_change skew_before skew_after
0 18-001A 0.000096 -30.0 -56.700714 -86.700714
1 18-001A 0.000115 -30.0 -54.693184 -84.693184
2 18-001A 0.000137 -30.0 -51.452210 -81.452210
3 18-001A 0.000164 -30.0 -61.983725 88.016275
4 18-001A 0.000196 -30.0 -60.874439 89.125561
... ... ... ... ... ...
1479 18-025A 0.487329 -30.0 -54.498901 -84.498901
1480 18-025A 0.582072 -30.0 -41.477337 -71.477337
1481 18-025A 0.695410 -30.0 -32.840732 -62.840732
1482 18-025A 0.830565 -30.0 -33.743310 -63.743310
1483 18-025A 0.992063 -30.0 -21.335294 -51.335294
[1484 rows x 5 columns]
The axial difference formula keeps theta comparisons honest across
the 180-degree wrap boundary.
Recommended Interpretation Workflow#
For a survey report, keep the phase-tensor interpretation explicit:
1from pathlib import Path
2
3import matplotlib.pyplot as plt
4
5from pycsamt.emtools import (
6 build_phase_tensor_table,
7 ensure_sites,
8 plot_phase_tensor_psection,
9 plot_phase_tensor_summary,
10 plot_skew_ellipt_density,
11 plot_theta_rose_grid,
12)
13
14sites = ensure_sites(
15 Path("data/AMT/WILLY_DATA/L18PLT"),
16 recursive=True,
17)
18
19period_range = (0.001, 10.0)
20skew_threshold = 3.0
21ellipt_threshold = 0.2
22
23pt = build_phase_tensor_table(sites, recursive=False)
24pt_band = pt[
25 (pt["period"] >= period_range[0])
26 & (pt["period"] <= period_range[1])
27]
28pt_band.to_csv("phase_tensor_table.csv", index=False)
29
30plot_phase_tensor_psection(
31 sites,
32 period_range=period_range,
33 c_by="skew",
34 skew_threshold=skew_threshold,
35 mark_3d=True,
36 recursive=False,
37)
38plt.gcf().savefig(
39 "tensor_recommended_psection.png",
40 dpi=200,
41 bbox_inches="tight",
42)
43plt.close()
44
45plot_theta_rose_grid(
46 sites,
47 n_bands=6,
48 recursive=False,
49)
50plt.gcf().savefig(
51 "tensor_recommended_rose_grid.png",
52 dpi=200,
53 bbox_inches="tight",
54)
55plt.close()
56
57plot_skew_ellipt_density(
58 sites,
59 band=period_range,
60 recursive=False,
61)
62plt.gcf().savefig(
63 "tensor_recommended_density.png",
64 dpi=200,
65 bbox_inches="tight",
66)
67plt.close()
68
69fig = plot_phase_tensor_summary(
70 sites,
71 period_range=period_range,
72 skew_threshold=skew_threshold,
73 ellipt_threshold=ellipt_threshold,
74 recursive=False,
75)
76fig.savefig("tensor_recommended_summary.png", dpi=200, bbox_inches="tight")
77plt.close(fig)
This workflow saves the table, plots the core ellipse section, checks strike by band, checks skew/ellipticity distribution, and creates a summary figure with the same thresholds.
Common Pitfalls#
Do not treat theta as an ordinary linear angle. It is axial. Use an
axial difference formula when comparing directions.
Do not rotate data only because a strike value exists. High skew,
unstable theta, or strong period dependence can make the rotation
misleading.
Do not hide 3-D behavior by editing Z too early. Plot the raw phase
tensor before applying antisymmetrization, balancing, clipping, or
rotation.
Do not interpret a missing map as missing tensor data. A map also needs
coordinates. If coordinates are missing, use pseudosections or pass an
explicit coords mapping.
Do not mix the two rotation-to-strike functions. Use
rotate_z_to_strike for the tensor editing helper exported from
pycsamt.emtools. Use the dedicated strike page for the strike module
workflow.
Worked Example#
The gallery example computes the phase-tensor table, builds the main pseudosections and rose diagrams, demonstrates the map and strip views, and audits impedance-tensor editing operations on real data.
Open the rendered gallery page here: Phase-tensor analysis and tensor editing (pycsamt.emtools.tensor).