15.2. Transfer Functions and Estimates#
TransferFunction stores one electromagnetic transfer function
as a matrix – shape (n_period, n_output_channels, n_input_channels)
– with its axes tied to named channels rather than hard-coded
Zxx/Zxy/Zyx/Zyy components.
StatisticalEstimate attaches variance or covariance to that
matrix in the form it was actually supplied, instead of coercing every
kind of uncertainty into the legacy z_err standard-error model.
Both are pure scientific objects: no EDI section, no XML element, is
ever stored inside them.
15.2.1. Matrix Shape and Validation#
TransferFunction accepts scalar, 1-D, 2-D, or 3-D data and
normalizes all of them to the canonical 3-D shape, using the declared
channel counts to decide what a lower-dimensional array must mean:
>>> from pycsamt.emtf import TransferFunction
>>> tf = TransferFunction(
... name="z", data=[[1.0, 2.0], [3.0, 4.0]],
... input_channels=("Hx", "Hy"), output_channels=("Ex", "Ey"),
... )
>>> print(tf.shape)
(1, 2, 2)
>>> print(tf.data.dtype)
complex128
>>> print(tf.data[0])
[[1.+0.j 2.+0.j]
[3.+0.j 4.+0.j]]
A bare 2-D array is treated as one period’s matrix and promoted to
(1, n_output, n_input) – real, plain floats are silently upcast
to complex because "z" normalizes to the impedance datatype,
whose declared data_kind is "complex" (more on that below). A
shape that does not match the declared channels raises immediately
with the exact mismatch, not a generic broadcasting error:
>>> try:
... TransferFunction(
... name="impedance", data=[[1, 2], [3, 4], [5, 6]],
... input_channels=("Hx",), output_channels=("Ex", "Ey"),
... )
... except ValueError as exc:
... print(type(exc).__name__, exc)
ValueError 2-D data are interpreted as one TF matrix; expected (2, 1), got (3, 2)
>>> try:
... TransferFunction(
... name="impedance", data=[[1, 2], [3, 4]],
... input_channels=("Hx", "Hx"), output_channels=("Ex", "Ey"),
... )
... except ValueError as exc:
... print(type(exc).__name__, exc)
ValueError duplicate channel name: 'Hx'
A scalar quantity – no meaningful input/output channels, like phase
tensor skew – uses empty channel tuples, and a 1-D array is then
promoted to (n_period, 1, 1) instead of a matrix:
>>> skew = TransferFunction(
... name="impedance_skew", data=[0.1, 0.2, 0.3],
... input_channels=(), output_channels=(),
... )
>>> print(skew.shape, skew.n_periods)
(3, 1, 1) 3
15.2.2. Name Normalization and Units#
The name passed to the constructor does not have to be the exact
tag: "z" above normalized to "impedance" because the datatype
registry (covered in full later on this page) resolves it. Once
resolved, the definition also fills in units automatically when
none was supplied, and enforces the declared numeric kind:
>>> print(tf.name, tf.units)
impedance [mV/km]/[nT]
>>> try:
... TransferFunction(
... name="impedance_skew", data=[0.1 + 0.5j, 0.2, 0.3],
... input_channels=(), output_channels=(),
... )
... except TypeError as exc:
... print(type(exc).__name__, exc)
TypeError impedance_skew is defined as real but complex data with non-zero imaginary values were supplied
>>> zero_imag = TransferFunction(
... name="impedance_skew", data=[0.1 + 0j, 0.2 + 0j, 0.3 + 0j],
... input_channels=(), output_channels=(),
... )
>>> print(zero_imag.data.dtype, zero_imag.data.ravel())
float64 [0.1 0.2 0.3]
A complex array with a genuinely non-zero imaginary part is rejected
for a "real"-kind type – that data almost certainly represents a
mistake, not a real quantity, and pycsamt refuses to discard
information silently. A complex array whose imaginary part happens to
be exactly zero is accepted and quietly narrowed to real, since no
information is actually lost.
15.2.3. The Missing-Value Policy#
The EMTF Document already showed that a real EDI-derived tipper can
contain genuine NaN entries, not zeros, wherever the source
station’s response was not estimated. That pattern is visible across
every component at once by comparing the real gv100 impedance and
tipper together:
import numpy as np
import matplotlib.pyplot as plt
z_valid = np.isfinite(imp.data.real) & np.isfinite(imp.data.imag)
t_valid = np.isfinite(tip.data.real) & np.isfinite(tip.data.imag)
rows, labels = [], []
for i, name in enumerate(["Zxx", "Zxy", "Zyx", "Zyy"]):
r, c = divmod(i, 2)
rows.append(z_valid[:, r, c])
labels.append(name)
for j, name in enumerate(["Tx", "Ty"]):
rows.append(t_valid[:, 0, j])
labels.append(name)
grid = np.array(rows, dtype=float)
fig, ax = plt.subplots(figsize=(9, 3))
im = ax.imshow(grid, aspect="auto", cmap="RdYlGn", vmin=0, vmax=1, interpolation="none")
ax.set_yticks(range(len(labels)))
ax.set_yticklabels(labels)
ax.set_xlabel("Period index (short -> long)")
ax.set_title("gv100: finite (green) vs NaN (red) per component and period")
fig.colorbar(im, ax=ax, ticks=[0, 1], label="finite", shrink=0.6)
fig.tight_layout()
All four real gv100 impedance components are finite at every
period. Tipper (Tx/Ty) is genuinely missing – not
zero – at 14 of 48 periods, split between the shortest and
longest ends of the band.#
Nothing in TransferFunction fills, interpolates, or masks
these gaps – it stores exactly the array it was given. Any downstream
code (plotting, inversion, QC) has to decide deliberately what to do
with a NaN; Processing and Quality’s
DataQuality is the tool that turns
this same pattern into a coverage percentage and a quality flag.
15.2.4. Statistical Estimates#
StatisticalEstimate attaches to a specific
TransferFunction through add_estimate(),
which validates that the estimate’s period axis actually matches the
parent matrix before accepting it:
>>> import numpy as np
>>> from pycsamt.emtf import StatisticalEstimate
>>> periods = np.array([1.0, 10.0, 100.0])
>>> z = np.ones((3, 2, 2), dtype=complex)
>>> tf3 = TransferFunction(
... name="impedance", data=z,
... input_channels=("Hx", "Hy"), output_channels=("Ex", "Ey"),
... periods=periods,
... )
>>> var = StatisticalEstimate(name="var", data=np.ones((3, 2, 2)), kind="variance")
>>> tf3.add_estimate(var)
TransferFunction(name='impedance', data=ndarray(shape=(3, 2, 2), dtype=complex128), input_channels=tuple(['Hx', 'Hy']), output_channels=tuple(['Ex', 'Ey']), units='[mV/km]/[nT]', periods=ndarray(shape=(3,), dtype=float64), estimates=dict(len=1, keys=[variance]), attrs=dict(len=0, keys=[]))
>>> print(tf3.get_estimate("VAR") is tf3.get_estimate("variance"))
True
Like add_transfer_function(),
add_estimate returns self and echoes the whole
TransferFunction repr when called bare, for the same
chaining reason discussed in The EMTF Document.
get_estimate() matches on the stored
.name ("VAR") or the semantic .kind ("variance")
interchangeably.
Attempting to attach a second estimate under the same key, or one whose period axis does not match, both fail with a specific message rather than corrupting the existing set:
>>> var2 = StatisticalEstimate(name="var", data=np.ones((3, 2, 2)) * 2, kind="variance")
>>> try:
... tf3.add_estimate(var2)
... except ValueError as exc:
... print(type(exc).__name__, exc)
ValueError estimate already exists: variance
>>> _ = tf3.add_estimate(var2, replace=True)
>>> print(tf3.get_estimate("variance").data[0, 0, 0])
2.0
>>> bad_est = StatisticalEstimate(name="bad", data=np.ones((5, 2, 2)), kind="variance")
>>> try:
... tf3.add_estimate(bad_est, key="bad")
... except ValueError as exc:
... print(type(exc).__name__, exc)
ValueError estimate period axis does not match TF data: 5 != 3
>>> print(tf3.get_estimate("does-not-exist"))
None
15.2.5. The Datatype Registry#
DataTypeDefinition describes one EMTF data type: its short
FCU-style code ("Z"), semantic tag ("impedance"), numeric
kind, channel families, units, and whether it is a primary
(measured/estimated directly) or derived (computed from another
type) quantity. list_emtf_datatypes() filters the registry by
that distinction:
>>> from pycsamt.emtf import get_emtf_datatype, list_emtf_datatypes
>>> z_def = get_emtf_datatype("Z")
>>> print(z_def.is_primary, z_def.is_derived)
True False
>>> rho_def = get_emtf_datatype("RHO")
>>> print(rho_def.is_primary, rho_def.is_derived, rho_def.derived_from)
False True impedance
>>> print(sorted(list_emtf_datatypes(intention="primary").keys()))
['impedance', 'interstation_impedance', 'interstation_transfer_functions', 'off_diagonal_impedance', 'tipper']
>>> print(len(list_emtf_datatypes(intention="derived")))
15
The primary set is deliberately small – just impedance, tipper, and interstation variants – while the fifteen derived entries cover phase tensor, apparent resistivity/phase, determinant and effective impedance, strike/skew/ellipticity, and induction-arrow quantities, each recording exactly which primary type it derives from.
aliases exist because the 2020 EMTF paper and EMTF FCU v4.1 do not
always agree on short names – interstation impedance and interstation
transfer functions are the clearest real example, registered under
"ZI"/"TI" but historically also known as "Q"/"P":
>>> zi = get_emtf_datatype("Q")
>>> print(zi.name, zi.tag, zi.aliases)
ZI interstation_impedance ('Q',)
>>> ti = get_emtf_datatype("interstation_transfer_functions")
>>> print(ti.name, ti.tag, ti.aliases)
TI interstation_transfer_functions ('P',)
Both lookups above – one starting from the historical alias, one
from the semantic tag – reach the same
DataTypeDefinition, which is exactly how "Z" and
"impedance" reached the same object in The EMTF Document.
15.2.6. Registering a Custom Datatype#
register_emtf_datatype() extends the registry without modifying
pycsamt itself – the mechanism the committed pycsamt.airborne
package (see the top-level user guide) would use to register an
admittance-style quantity that has no equivalent in ground MT:
>>> from pycsamt.emtf import DataTypeDefinition, register_emtf_datatype
>>> admittance_def = DataTypeDefinition(
... name="Y", tag="admittance", data_kind="complex",
... input_kind="E", output_kind="H",
... units="[nT]/[mV/km]", intention="derived",
... description="Airborne-style admittance (inverse impedance convention)",
... derived_from="impedance",
... )
>>> register_emtf_datatype(admittance_def)
DataTypeDefinition(name='Y', tag='admittance', data_kind='complex', input_kind='E', output_kind='H', units='[nT]/[mV/km]', intention='derived', description='Airborne-style admittance (inverse impedance ...', ...)
>>> periods_y = np.array([25.0, 385.0])
>>> y = np.ones((2, 2, 2), dtype=complex) * 0.01
>>> tf_y = TransferFunction(
... name="Y", data=y,
... input_channels=("Ex", "Ey"), output_channels=("Hx", "Hy"),
... periods=periods_y,
... )
>>> from pycsamt.emtf import EMTF
>>> doc = EMTF()
>>> _ = doc.add_transfer_function(tf_y)
>>> print(sorted(doc.transfer_functions.keys()))
['admittance']
>>> print(doc.get_transfer_function("Y") is tf_y)
True
>>> try:
... register_emtf_datatype(admittance_def)
... except ValueError as exc:
... print(type(exc).__name__, exc)
ValueError EMTF data type already registered: admittance
The registration itself is process-global – once "admittance" is
registered, any TransferFunction named "Y" or
"admittance" resolves to it for the rest of the process, which is
exactly what lets get_transfer_function("Y") find tf_y above.
Registering the same tag twice fails unless overwrite=True is
passed explicitly, so a project cannot accidentally redefine a type
another part of the same process already depends on.
15.2.7. Choosing the Right Object#
Need |
Object/function |
Notes |
|---|---|---|
One matrix-valued transfer function |
Shape normalized from scalar/1-D/2-D/3-D input automatically. |
|
Variance or covariance attached to a TF |
Stored as supplied; never coerced into the legacy |
|
Look up what a short code or tag means |
Matches tag, primary code, or any registered alias. |
|
List primary or derived types |
|
|
Add a project-specific data type |
Global for the process; requires |
15.2.8. Next Steps#
The EMTF Document covers the
EMTFcontainer these transfer functions attach to;EDI Interoperability covers how a real EDI file’s
Z/TIPblocks becomeTransferFunctionobjects in the first place;Rotation and Covariance covers how a
TransferFunctionand its attached estimates transform together under rotation.