13.5. Site Editing#

The pycsamt.site.edit module contains practical editing helpers for EDI-like objects and site collections. These helpers are designed for survey preparation: rotate impedance tensors, reduce frequency ranges, normalize station names, assign coordinates, fill missing arrays, and recompute derived apparent resistivity and phase quantities after edits.

The editing functions are deliberately tolerant. When an optional section is missing or has an incompatible shape, most functions skip that part and keep processing the rest of the object. This is helpful for real field data, where some stations may lack tipper, errors, or derived arrays.

For a complete survey-level cleanup that reads EDI files, applies edits, recomputes derived quantities, preserves line folders, and writes new pyCSAMT-authored EDI files, use EDI Recompute Workflow. The functions on this page are the lower-level building blocks used by that workflow.

The examples below use a small synthetic EDI-like station so the outputs can be reproduced exactly. In production, replace DemoSite(...) with pycsamt.seg.edi.EDIFile, pycsamt.site.base.Site, or a pycsamt.site.base.Sites collection loaded from survey files.

>>> import copy
>>> import numpy as np

>>> class Head:
...     def __init__(self, dataid, lat=np.nan, lon=np.nan, elev=np.nan):
...         self.dataid = dataid
...         self.station = dataid
...         self.lat = lat
...         self.lon = lon
...         self.elev = elev
...
>>> class ZBlock:
...     def __init__(self, freq, z):
...         self.freq = np.asarray(freq, float)
...         self.z = np.asarray(z, complex)
...         self.rho = None
...         self.phase = None
...
...     def compute_resistivity_phase(self):
...         mu0 = 4 * np.pi * 1e-7
...         self.rho = (
...             np.abs(self.z) ** 2
...             / (mu0 * 2 * np.pi * self.freq[:, None, None])
...         )
...         self.phase = np.degrees(np.angle(self.z))
...
>>> class TipBlock:
...     def __init__(self, tipper):
...         self.tipper = np.asarray(tipper, complex)
...
>>> class DemoSite:
...     def __init__(self, name, scale=1.0):
...         self.name = name
...         self.Head = Head(name)
...         freq = np.array([1.0, 3.0, 10.0, 30.0, 100.0, 300.0, 1000.0])
...         base = (1 + 0.2j) * np.sqrt(freq / 100.0) * scale
...         z = np.zeros((freq.size, 2, 2), dtype=complex)
...         z[:, 0, 0] = 0.05 * base
...         z[:, 1, 1] = -0.03 * base
...         z[:, 0, 1] = base * (1 + 0.05j)
...         z[:, 1, 0] = -0.8 * base * (1 - 0.04j)
...         z[2, 0, 1] = np.nan + 0j
...         self.Z = ZBlock(freq, z)
...         self.Tip = TipBlock(
...             np.column_stack([
...                 0.10 * scale * np.ones(freq.size),
...                 0.04j * np.linspace(1.0, 1.5, freq.size),
...             ])
...         )
...
...     def get_section(self, name):
...         if str(name).lower() == "head":
...             return self.Head
...         return getattr(self, name, None)
...
...     def __copy__(self):
...         new = type(self).__new__(type(self))
...         new.__dict__ = copy.deepcopy(self.__dict__)
...         return new
...

13.5.1. Editing Map#

Function

Scope

Main purpose

rotate()

One site

Rotate impedance tensors and tipper values by an azimuthal angle.

select_freq()

One site

Subset frequency-indexed arrays by range, indices, or boolean mask.

rename()

One site

Set an explicit station name or apply a naming policy.

set_coords()

One site

Update latitude, longitude, and elevation in the EDI header.

fill_missing()

One site

Fill or allocate missing Z and tipper arrays.

recompute_res_phase()

One site

Recompute apparent resistivity and phase from the impedance tensor.

rotate_all()

Collection

Rotate every site in a collection.

select_freq_all()

Collection

Apply frequency subsetting to every site.

rename_all()

Collection

Rename stations across a collection.

set_coords_all()

Collection

Assign coordinates from a callable, mapping, or table holder.

set_coords_from_table()

Collection

Assign coordinates from CSV, DataFrame, structured array, or row list.

set_coords_from_en()

One site

Project easting/northing coordinates and store lon/lat on a site.

13.5.2. Copy Versus In-Place Editing#

Most editing helpers accept inplace. The default is False:

>>> from pycsamt.site.edit import rename
>>> from pycsamt.site.utils import station_name

>>> edi = DemoSite("S01")

>>> edited = rename(edi, name="LINE01_S01", inplace=False)

>>> print(edited is edi)
False
>>> print(station_name(edi), station_name(edited))
S01 LINE01_S01

Use inplace=True when you intentionally want to mutate the provided object:

>>> _ = rename(edi, name="LINE01_S01", inplace=True)
>>> print(station_name(edi))
LINE01_S01

For collection helpers, inplace=False returns a new pycsamt.site.base.Sites wrapper. The original EDI objects are left untouched as far as the helper can manage. Use inplace=True only when the calling workflow owns the input objects.

13.5.3. Tensor Rotation#

rotate() rotates the impedance tensor and, when present, the tipper. For the impedance tensor, pyCSAMT applies a horizontal coordinate transform. At each frequency \(f\),

(1)#\[\begin{split}\mathbf{Z}'(f) = \mathbf{R}(\theta)\,\mathbf{Z}(f)\,\mathbf{R}^{-1}(\theta), \qquad \mathbf{R}(\theta) = \begin{bmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{bmatrix}.\end{split}\]

This keeps the electric and magnetic axes consistent while expressing the same tensor in a rotated field coordinate frame. The inverse appears on the right because the tensor maps horizontal magnetic-field components to horizontal electric-field components.

>>> import numpy as np
>>> from pycsamt.site.edit import rotate

>>> edi = DemoSite("S01")
>>> before = edi.Z.z[0].copy()
>>> rotated = rotate(edi, angle_deg=30.0, inplace=False)

>>> print("before:")
before:
>>> print(np.round(before, 3))
[[ 0.005+0.001j  0.099+0.025j]
 [-0.081-0.013j -0.003-0.001j]]
>>> print("after:")
after:
>>> print(np.round(rotated.Z.z[0], 3))
[[ 0.011+0.006j  0.091+0.021j]
 [-0.089-0.017j -0.009-0.005j]]
>>> print("original unchanged:", np.allclose(edi.Z.z[0], before))
original unchanged: True

Because equation (1) is a similarity transformation, it changes individual tensor components but preserves the trace and determinant at each frequency. Those invariants provide a compact numerical check that the operation changed coordinates rather than the underlying response:

>>> after = rotated.Z.z[0]
>>> print(np.allclose(np.trace(after), np.trace(before)))
True
>>> print(np.allclose(np.linalg.det(after), np.linalg.det(before)))
True

The helper checks common tensor attribute names such as z, impedance, and _z. Tipper arrays may live under T, TIP, Tip, or sometimes as a tipper-like attribute under Z.

Error arrays are rotated with a magnitude-only approximation using absolute values of the rotation matrices. If \(\boldsymbol{\sigma}_Z\) contains component-wise standard deviations, the approximation is

(2)#\[\boldsymbol{\sigma}_{Z'} \approx |\mathbf R|\,\boldsymbol{\sigma}_Z\,|\mathbf R^{-1}|.\]

Equation (2) ignores cross-component covariances. Treat it as practical uncertainty propagation, not as a full rotation of a covariance matrix.

Rotate a whole collection with rotate_all():

>>> from pycsamt.site.edit import rotate_all

>>> collection = [DemoSite("S01"), DemoSite("S02", scale=1.2)]
>>> rotated_sites = rotate_all(collection, angle_deg=30.0)
>>> print(len(rotated_sites.as_list()))
2
>>> print(np.round(rotated_sites.as_list()[0].Z.z[0], 3))
[[ 0.011+0.006j  0.091+0.021j]
 [-0.089-0.017j -0.009-0.005j]]

13.5.4. Frequency Subsetting#

select_freq() subsets a station along its frequency axis and slices aligned arrays together. It handles Z, Z errors, apparent resistivity, phase, tipper, tipper errors, and several common aliases.

For a frequency vector \(\mathbf f=(f_1,\ldots,f_N)\), range selection constructs one mask

(3)#\[m_k=(f_k\geq f_{\min})\land(f_k\leq f_{\max}),\]

and applies the same \(m_k\) to every frequency-indexed array. Reusing one mask is essential: slicing impedance and tipper independently could associate values with the wrong frequencies even when their final row counts match.

Keep a frequency range:

>>> from pycsamt.site.edit import select_freq

>>> edi = DemoSite("S01")
>>> band = select_freq(
...     edi,
...     fmin=10.0,
...     fmax=300.0,
...     inplace=False,
... )
>>> print(band.Z.freq.tolist())
[10.0, 30.0, 100.0, 300.0]
>>> print(band.Z.z.shape, band.Tip.tipper.shape)
(4, 2, 2) (4, 2)

The four retained frequencies therefore correspond to four rows in both arrays; this shape check is a quick audit of the alignment promised by equation (3).

Keep explicit row indices:

>>> edges = select_freq(edi, keep=[0, -1])
>>> print(edges.Z.freq.tolist())
[1.0, 1000.0]

Use a boolean mask:

>>> import numpy as np

>>> freq = np.asarray(edi.Z.freq)
>>> mask = freq >= 30.0
>>> high = select_freq(edi, keep=mask)
>>> print(high.Z.freq.tolist())
[30.0, 100.0, 300.0, 1000.0]

When keep is provided, fmin and fmax are ignored. Use select_freq_all() for collections:

>>> from pycsamt.site.edit import select_freq_all

>>> trimmed = select_freq_all(
...     collection,
...     fmin=10.0,
...     fmax=300.0,
... )
>>> print(trimmed.as_list()[0].Z.freq.tolist())
[10.0, 30.0, 100.0, 300.0]

The same synthetic station can also be plotted before and after rotation and frequency subsetting. The selected points are drawn in one panel so the frequency mask is easy to audit visually.

>>> import matplotlib.pyplot as plt
>>> from pycsamt.site.edit import fill_missing, rotate, select_freq

>>> raw = DemoSite("S01")
>>> rotated = rotate(raw, angle_deg=30.0, inplace=False)
>>> raw_filled = fill_missing(raw, how="zero", components=("Z",), inplace=False)
>>> band = select_freq(raw_filled, fmin=10.0, fmax=300.0, inplace=False)

>>> fig, ax = plt.subplots(1, 2, figsize=(8, 3.2), constrained_layout=True)

>>> _ = ax[0].plot(raw.Z.freq, np.abs(raw.Z.z[:, 0, 1]), marker="o", label="raw Zxy")
>>> _ = ax[0].plot(
...     rotated.Z.freq,
...     np.abs(rotated.Z.z[:, 0, 1]),
...     marker="s",
...     label="rotated Zxy",
... )
>>> _ = ax[0].set_xscale("log")
>>> _ = ax[0].set_xlabel("frequency (Hz)")
>>> _ = ax[0].set_ylabel("|Zxy|")
>>> _ = ax[0].set_title("Rotation effect")
>>> _ = ax[0].legend(frameon=False)

>>> _ = ax[1].plot(
...     raw_filled.Z.freq,
...     np.abs(raw_filled.Z.z[:, 0, 1]),
...     marker="o",
...     label="all rows",
... )
>>> _ = ax[1].scatter(
...     band.Z.freq,
...     np.abs(band.Z.z[:, 0, 1]),
...     s=80,
...     label="selected band",
... )
>>> _ = ax[1].set_xscale("log")
>>> _ = ax[1].set_xlabel("frequency (Hz)")
>>> _ = ax[1].set_title("Frequency subset")
>>> _ = ax[1].legend(frameon=False)

>>> for axis in ax:
...     axis.grid(True, alpha=0.25)
...
>>> fig.savefig("editing_rotation_frequency.png", dpi=160)
Two-panel plot showing rotated impedance magnitude and selected frequency rows for a synthetic station.

A compact check of the two edits that most often affect numerical interpretation: tensor rotation and frequency subsetting.#

The left panel shows a small change in \(|Z_{xy}|\) after rotation. This is redistribution among tensor components, not a change in the measured electromagnetic response, as the invariant checks above confirm. In the right panel, the selected band contains the 10–300 Hz rows only. The point at zero near 10 Hz comes from the deliberately zero-filled missing value; it is a processing artifact and must not be interpreted as a measured vanishing impedance.

13.5.5. Station Renaming#

rename() updates common station identifiers in the EDI header and mirrors the result to edi.name when possible. This makes later lookup through pycsamt.site.base.Site and pycsamt.site.base.Sites more consistent.

Set an explicit name:

>>> from pycsamt.site.edit import rename

>>> edi = DemoSite("S01")
>>> renamed = rename(edi, name="LINE01_S01")
>>> print(station_name(renamed))
LINE01_S01

Apply a policy to the current station name:

>>> renamed = rename(
...     edi,
...     policy=lambda old: f"LINE01_{old}",
... )
>>> print(station_name(renamed), renamed.Head.dataid)
LINE01_S01 LINE01_S01

If both name and policy are provided, the explicit name wins. Renaming does not change the on-disk filename. Use export templates later if you want output filenames to follow the new station names.

For collections, use rename_all():

>>> from pycsamt.site.edit import rename_all

>>> renamed_sites = rename_all(
...     collection,
...     name_fn=lambda edi: f"L01_{station_name(edi)}",
... )
>>> print([station_name(edi) for edi in renamed_sites.as_list()])
['L01_S01', 'L01_S02']

Prefer name_fn when the original EDI files have duplicate station names. It receives the EDI object and can use file paths, metadata, or external state to generate unique identifiers.

13.5.6. Coordinate Editing#

set_coords() updates latitude, longitude, and elevation on one site. Only the values provided are changed.

>>> from pycsamt.site.edit import set_coords

>>> moved = set_coords(
...     edi,
...     lat=35.125,
...     lon=12.750,
...     elev=1234.0,
... )
>>> print(moved.Head.lat, moved.Head.lon, moved.Head.elev)
35.125 12.75 1234.0

To apply coordinates to many sites, use set_coords_all().

From a mapping keyed by station name:

>>> from pycsamt.site.edit import set_coords_all

>>> collection = [DemoSite("S01"), DemoSite("S02", scale=1.2)]
>>> coords = {
...     "S01": (35.125, 12.750, 1234.0),
...     "S02": (35.200, 12.900, 1180.0),
... }

>>> updated = set_coords_all(collection, coords)
>>> for site in updated.as_list():
...     print(site.name, site.Head.lat, site.Head.lon, site.Head.elev)
...
S01 35.125 12.75 1234.0
S02 35.2 12.9 1180.0

From a callable:

>>> collection = [DemoSite("S01"), DemoSite("S02", scale=1.2)]

>>> def lookup(edi):
...     name = getattr(edi, "name", "")
...     if name == "S01":
...         return 35.125, 12.750, 1234.0
...     return None
...
>>> updated = set_coords_all(collection, lookup)
>>> for site in updated.as_list():
...     print(site.name, site.Head.lat, site.Head.lon, site.Head.elev)
...
S01 35.125 12.75 1234.0
S02 nan nan nan

From an object exposing a .frame attribute:

>>> import pandas as pd

>>> class CoordinateTable:
...     def __init__(self, frame):
...         self.frame = frame
...
>>> frame = pd.DataFrame(
...     {
...         "station": ["S01", "S02"],
...         "lat": [35.125, 35.200],
...         "lon": [12.750, 12.900],
...         "elev": [1234.0, 1180.0],
...     }
... )
>>> updated = set_coords_all(collection, CoordinateTable(frame))
>>> for site in updated.as_list():
...     print(site.name, site.Head.lat, site.Head.lon, site.Head.elev)
...
S01 35.125 12.75 1234.0
S02 35.2 12.9 1180.0

13.5.7. Coordinate Tables#

set_coords_from_table() is the high-level table loader. It accepts:

  • path to a CSV file;

  • path to a whitespace-delimited text file;

  • pandas.DataFrame;

  • NumPy structured array;

  • list of dictionaries;

  • list of row-like tuples that pandas can convert to a frame.

With standard geographic columns:

>>> import pandas as pd
>>> from pycsamt.site.edit import set_coords_from_table

>>> collection = [DemoSite("S01"), DemoSite("S02", scale=1.2)]
>>> table = pd.DataFrame(
...     {
...         "station": ["S01", "S02"],
...         "lat": [35.125, 35.200],
...         "lon": [12.750, 12.900],
...         "elev": [1234.0, 1180.0],
...     }
... )

>>> updated = set_coords_from_table(collection, table)
>>> for site in updated.as_list():
...     print(site.name, site.Head.lat, site.Head.lon, site.Head.elev)
...
S01 35.125 12.75 1234.0
S02 35.2 12.9 1180.0

For a loaded pycsamt.site.base.Sites collection, this updates each matching station header and returns a new Sites wrapper unless inplace=True is used.

The resolver understands common aliases:

Canonical field

Accepted names

station

station, name, site, id.

lat

lat, latitude.

lon

lon, long, longitude.

elev

elev, elevation, z.

easting

easting, x.

northing

northing, y.

Use an explicit column map when a field table uses non-standard names:

>>> import pandas as pd

>>> collection = [DemoSite("S01"), DemoSite("S02", scale=1.2)]
>>> alias_table = pd.DataFrame(
...     {
...         "name": ["S01", "S02"],
...         "latitude": [35.125, 35.200],
...         "long": [12.750, 12.900],
...         "elevation": [1234.0, 1180.0],
...     }
... )

>>> updated = set_coords_from_table(
...     collection,
...     alias_table,
...     columns={
...         "station": "name",
...         "lat": "latitude",
...         "lon": "long",
...         "elev": "elevation",
...     },
... )
>>> for site in updated.as_list():
...     print(site.name, site.Head.lat, site.Head.lon, site.Head.elev)
...
S01 35.125 12.75 1234.0
S02 35.2 12.9 1180.0

When both geographic coordinates and easting/northing are present, latitude and longitude are preferred. The coordinate table is a reproducibility artifact: keep it with the processing notes so reviewers know exactly which station positions were written into the edited files.

13.5.8. Easting/Northing Conversion#

If a table provides projected coordinates, pass crs_from:

>>> table = pd.DataFrame(
...     {
...         "station": ["S10", "S11"],
...         "easting": [400000.0, 401250.0],
...         "northing": [5750000.0, 5750400.0],
...         "elev": [250.0, 252.0],
...     }
... )

>>> updated = set_coords_from_table(
...     collection,
...     table,
...     crs_from="EPSG:32631",
... )

Projection uses pyproj. If pyproj is not installed and projected coordinates must be converted, the helper raises ImportError. The coordinate reference system controls the numerical meaning of easting and northing values; the same pair of numbers in a different CRS can project to a different geographic position.

For one site, use set_coords_from_en():

>>> from pycsamt.site.edit import set_coords_from_en

>>> projected = set_coords_from_en(
...     edi,
...     easting=400000.0,
...     northing=5750000.0,
...     crs_from="EPSG:32631",
...     elev=250.0,
... )
>>> print(round(projected.Head.lat, 3), round(projected.Head.lon, 3))
51.892 1.547

13.5.9. Missing Data Filling#

fill_missing() fills or allocates missing Z and tipper arrays. It is often used before diagnostics or before recomputing derived quantities.

>>> from pycsamt.site.edit import fill_missing

>>> filled = fill_missing(
...     edi,
...     how="zero",
...     components=("Z",),
... )
>>> print(int(np.isnan(filled.Z.z).sum()))
0
>>> print(filled.Z.z[2, 0, 1])
0j

Available fill policies are:

how="zero"

Replace non-finite values with numeric zeros.

how="nan"

Replace non-finite values with NaN.

The components argument is case-insensitive and accepts "Z" and "Tip". Z arrays are expected to have shape (n_freq, 2, 2). Tipper arrays are expected to have shape (n_freq, 2). The number of rows is inferred from the frequency vector.

Use this function carefully. Zero-filled tensors are convenient for tests and some robust workflows, but zeros are not measured values. Mathematically, the policy replaces each non-finite array entry \(x_i\) by

(4)#\[\begin{split}x_i' = \begin{cases} x_i, & x_i \text{ is finite},\\ 0, & x_i \text{ is not finite and } \texttt{how="zero"},\\ \mathrm{NaN}, & x_i \text{ is not finite and } \texttt{how="nan"}. \end{cases}\end{split}\]

Equation (4) makes clear that how="nan" does not infer a value—it preserves the gap explicitly. This can be verified on the same synthetic station:

>>> gap_preserved = fill_missing(edi, how="nan", components=("Z",))
>>> print(int(np.isnan(gap_preserved.Z.z).sum()))
1
>>> print(np.isnan(gap_preserved.Z.z[2, 0, 1]))
True

Keep the chosen policy in the processing log when missing data have been filled. For inversion or statistical estimation, preserving NaN and letting an explicit mask reject the row is generally safer than silently turning the gap into a zero-amplitude observation.

13.5.10. Recomputing Resistivity And Phase#

recompute_res_phase() calls the available Z-section compute_resistivity_phase() method. It is best used after changing frequency rows, tensor values, or masks.

>>> from pycsamt.site.edit import recompute_res_phase

>>> edited = select_freq(edi, fmin=1.0, fmax=1000.0)
>>> edited = fill_missing(edited, how="zero", components=("Z",))
>>> edited = recompute_res_phase(edited)
>>> print(edited.Z.rho.shape)
(7, 2, 2)
>>> print(round(float(edited.Z.phase[4, 0, 1]), 3))
14.172

The function is best-effort. If the Z section is missing, incompatible, or does not expose a recomputation method, the object is returned unchanged. For each component, a typical recomputation stores

(5)#\[\rho_a(f) = \frac{|Z(f)|^2}{\mu_0\,2\pi f}, \qquad \phi(f) = \arg(Z(f))\,\frac{180}{\pi},\]

so apparent resistivity and phase remain synchronized with the current frequency rows and tensor values. The \(Z\) in equation (5) is in SI impedance units; an EDI reader must apply the corresponding field-unit conversion before using this form.

13.5.11. Practical Preparation Workflow#

The following pattern is common before diagnostics or inversion preparation. For a real survey, sites would normally come from EDICollection.from_sources("data/raw_edi") or Sites.from_path(...); the synthetic list keeps the example output reproducible.

>>> import pandas as pd
>>> from pycsamt.site.edit import (
...     fill_missing,
...     recompute_res_phase,
...     rename_all,
...     rotate_all,
...     select_freq_all,
...     set_coords_from_table,
... )
>>> from pycsamt.site.utils import station_name

>>> sites = [
...     DemoSite("S01"),
...     DemoSite("S02", scale=1.2),
...     DemoSite("S03", scale=0.8),
... ]

>>> sites = rename_all(
...     sites,
...     name_fn=lambda edi: f"L01_{station_name(edi)}",
... )
>>> coords = pd.DataFrame(
...     {
...         "station": ["L01_S01", "L01_S02", "L01_S03"],
...         "lat": [35.1, 35.2, 35.3],
...         "lon": [12.7, 12.8, 12.9],
...         "elev": [100.0, 110.0, 120.0],
...     }
... )
>>> sites = set_coords_from_table(sites, coords)
>>> sites = rotate_all(sites, angle_deg=30.0)
>>> sites = select_freq_all(sites, fmin=10.0, fmax=300.0)

>>> rows = []
>>> for edi in sites.as_list():
...     edi = fill_missing(edi, how="nan", components=("Z",))
...     edi = recompute_res_phase(edi)
...     rows.append(
...         {
...             "name": station_name(edi),
...             "freq_rows": len(edi.Z.freq),
...             "nan_z": int(np.isnan(edi.Z.z).sum()),
...         }
...     )
...
>>> print(pd.DataFrame(rows).to_string(index=False))
   name  freq_rows  nan_z
L01_S01          4      4
L01_S02          4      4
L01_S03          4      4

Each stage should be recorded in a processing log, especially rotation angles, frequency bands, coordinate sources, and missing-data fill policies.

13.5.12. Common Mistakes#

Using zero fill as if it were measured data

fill_missing(..., how="zero") is useful for deterministic tests and some defensive algorithms. It should not silently replace failed field measurements in scientific interpretation.

Renaming without checking duplicates

A policy such as lambda name: "X_" + name preserves duplicates if the original station names were duplicated. Use name_fn with file stems or a station table when uniqueness matters.

Mixing frequency masks from different stations

keep=[...] and boolean masks are applied row-wise. A mask built from one station may not represent the same frequencies on another station if their frequency axes differ.

Forgetting to recompute derived arrays

After changing Z or frequency rows, apparent resistivity and phase may be stale. Run recompute_res_phase() when downstream steps use derived arrays.

Ignoring CRS metadata

Easting/northing values are meaningless without their source CRS. Always provide the correct crs_from when using projected coordinates.

13.5.13. Next Pages#

Continue with: