Instrument

From radiance to what a sensor records: the Instrument protocol, BandInstrument, JAX band integration with a differentiable spectral shift, and GaussianInstrument, Gaussian bands shifted analytically in cm⁻¹ or nm (spectral); spectral response functions and NumPy band integration (srf), point-spread and ground-sampling operators (psf); the NoiseModel protocol and its diagonal, SNR and NEΔT noise (noise).

instrument

From radiance to what a sensor records.

srf integrates hyperspectral radiance over spectral response functions (SpectralResponseFunction); psf holds the JAX point-spread and ground-sampling operators with their hand-written adjoints. spectral holds the Instrument protocol and BandInstrument, the JAX band integration with a differentiable spectral shift (radtrax#36), and GaussianInstrument, Gaussian bands shifted analytically in cm⁻¹ or nm (radtrax#96). noise holds the NoiseModel protocol, the measurement-noise covariance (radtrax#108), with NEdTNoise for thermal bands (radtrax#249).

Ported from plumax 0.1.2 at 28f257e (plumax.radtran.srf, plumax.radtran.instrument).

spectral

The Instrument protocol and a band instrument with a spectral shift.

An instrument turns radiance on the model's wavenumber grid into what the sensor records per band (radtrax · 3). v0's instrument integrates over each band's spectral response and lets the wavenumber scale shift:

.. math:: y_b = \sum_g f_{bg}\, L(ν_g + δ),

with f the row-normalised response matrix (B, G), L the radiance (G,) and δ the spectral shift [cm⁻¹]: a positive δ moves every band's response to higher wavenumber (shorter wavelength). A uniform shift in wavenumber is not one in wavelength (δλ = −λ² δν / 10⁷): across 1600–2500 nm the same δν moves the bands by amounts 2.4 times apart, so a calibration given per nm converts band by band.

L(ν + δ) is linear interpolation on the grid, so y is piecewise linear in δ with kinks at multiples of the grid step; on a node the derivative is the slope of the segment above it. On a 0.01 cm⁻¹ grid that derivative is within 0.2 % of an analytically shifted SRF's and a Gauss–Newton step from δ = 0 converges; ∂²y/∂δ² is zero, so a method that needs the curvature in δ needs the analytic shift. Beyond the grid's ends the radiance is held at its edge value, so keep each band's response well inside the grid. A changing FWHM and the PSF and GSD of :mod:radtrax.instrument.psf join it in v1 of the roadmap.

The response matrix comes from the NumPy :class:~radtrax.instrument.SpectralResponseFunction (exact Gaussian FWHM, trapezoidal weights) through :meth:BandInstrument.from_srf, so both integrate a band the same way; this module only moves it onto the gradient path. Ported from the radtrax · 3 design notebook's Instrument.

:class:GaussianInstrument (radtrax#96) builds the same Gaussian bands in JAX from their centres and FWHM in nm, and shifts the response rather than the spectrum:

.. math:: y_b = \frac{\sum_g w_g\, f_b(ν_g)\, L(ν_g)}{\sum_g w_g\, f_b(ν_g)},\qquad f_b(ν) = \exp\Big(-\frac{(10^7/(ν - δ) - δ_λ - c_b)^2}{2σ_b^2}\Big),

with w the trapezoidal weights on the wavelengths 10⁷/ν, c_b the centre, σ_b the FWHM over 2√(2 ln 2), and a shift δ in cm⁻¹ or δ_λ in nm. At zero shift it is BandInstrument.from_srf of the same bands; y is smooth in both shifts and in the FWHM, so ∂²y/∂δ² is the band's curvature, not zero, and a calibration per nm needs no conversion.

Shapes: G wavenumbers, B bands.

Instrument

Bases: Module

Protocol: radiance on a wavenumber grid in, one value per band out.

An abstract Equinox module, so filter_jit compiles each concrete instrument and a forward model can hold any of them.

The contract: an instrument is linear in radiance, y = A L with A its response matrix (B, G) (:meth:response_matrix), so I(a L₁ + b L₂) = a I(L₁) + b I(L₂). Its parameters (a shift, a FWHM) may enter A non-linearly. Band brightness temperature and NEΔT (:mod:radtrax.radiometry.brightness) rely on it.

Attributes:

Name	Type	Description
nu_cm	AbstractVar[Float[Array, ' G']]	The wavenumber grid the instrument reads [cm⁻¹], (G,).

response_matrix

response_matrix() -> Float[Array, 'B G']

The response matrix A = ∂y/∂L, (B, G).

By default jax.jacrev of the call at zero radiance (B reverse passes, exact for a linear instrument); an instrument that holds A in closed form overrides it.

Source code in src/radtrax/instrument/spectral.py

def response_matrix(self) -> Float[Array, "B G"]:
    """The response matrix ``A = ∂y/∂L``, ``(B, G)``.

    By default ``jax.jacrev`` of the call at zero radiance (``B``
    reverse passes, exact for a linear instrument); an instrument that
    holds ``A`` in closed form overrides it.
    """
    return jax.jacrev(self.__call__)(jnp.zeros_like(self.nu_cm))

BandInstrument

BandInstrument(nu_cm: Float[ArrayLike, ' G'], response: Float[ArrayLike, 'B G'], shift_cm: Float[ArrayLike, ''] = 0.0, band_names: tuple[str, ...] | None = None)

Bases: Instrument

Band integration with a spectral shift, y_b = Σ_g f_bg L(ν_g + δ).

Attributes:

Name	Type	Description
nu_cm	Float[Array, ' G']	Wavenumber grid [cm⁻¹], strictly increasing, (G,).
response	Float[Array, 'B G']	Spectral response per band on nu_cm, each row summing to 1, (B, G).
shift_cm	Float[Array, '']	Spectral shift δ [cm⁻¹], a differentiable leaf.
band_names	tuple[str, ...]	One label per band.

Examples:

>>> import jax.numpy as jnp
>>> instrument = BandInstrument(
...     jnp.array([4300.0, 4300.1, 4300.2]),
...     jnp.array([[0.5, 0.5, 0.0], [0.0, 0.5, 0.5]]),
... )
>>> [float(y) for y in instrument(jnp.array([1.0, 2.0, 4.0]))]
[1.5, 3.0]

Source code in src/radtrax/instrument/spectral.py

def __init__(
    self,
    nu_cm: Float[ArrayLike, " G"],
    response: Float[ArrayLike, "B G"],
    shift_cm: Float[ArrayLike, ""] = 0.0,
    band_names: tuple[str, ...] | None = None,
) -> None:
    self.nu_cm = jnp.asarray(nu_cm, dtype=float)
    self.response = jnp.asarray(response, dtype=float)
    self.shift_cm = jnp.asarray(shift_cm, dtype=float)
    n_bands = self.response.shape[0] if self.response.ndim == 2 else 0
    self.band_names = (
        tuple(f"B{b}" for b in range(n_bands))
        if band_names is None
        else tuple(band_names)
    )

n_bands property

n_bands: int

Number of bands B.

from_srf classmethod

from_srf(srf: SpectralResponseFunction, nu_cm: Float[ArrayLike, ' G'], *, shift_cm: Float[ArrayLike, ''] = 0.0, per_um: bool = False) -> BandInstrument

From a :class:~radtrax.instrument.SpectralResponseFunction built on the wavelengths of nu_cm.

srf.wavelengths_hr_nm must be 1e7 / nu_cm in ascending order (nu_cm reversed), so the response is the SRF's own matrix, its columns reversed onto the wavenumber grid.

With per_um=True the instrument takes radiance per cm⁻¹ and returns the band average per µm, which is what a thermal Level-1 product reports (Landsat TIRS, ECOSTRESS). The SRF's rows w are weights in wavelength (R Δλ, summing to 1), so the band average of \(L_\lambda = L_{\tilde\nu} \tilde\nu^2 / 10^4\) is

\[ \bar L_{\lambda,b} = \frac{\int R_b L_\lambda \, d\lambda}{\int R_b \, d\lambda} = \sum_g w_{bg} \frac{\tilde\nu_g^2}{10^4} L_{\tilde\nu}(\tilde\nu_g), \]

and the factor is folded into the response, \(f_{bg} = w_{bg} \tilde\nu_g^2 / 10^4\), whose rows then sum to the band's mean of \(\tilde\nu^2 / 10^4\) (about 84 at 10.9 µm), not to 1. Averaging \(L_{\tilde\nu}\) with w and converting afterwards is not the same: across a 0.9 µm band at 10.9 µm, \(\tilde\nu^2\) changes by about ±8 %. The factor is taken at the unshifted grid, so a shift \(\delta\) adds a relative error of about \(2\delta/\tilde\nu\) (2e-4 for 0.1 cm⁻¹ at 900 cm⁻¹).

Parameters:

Name	Type	Description	Default
srf	SpectralResponseFunction	The SRF, on wavelengths_hr_nm = 1e7 / nu_cm[::-1].	required
nu_cm	Float[ArrayLike, ' G']	Wavenumber grid [cm⁻¹], strictly increasing, (G,).	required
shift_cm	Float[ArrayLike, '']	Spectral shift δ [cm⁻¹].	0.0
per_um	bool	Return band radiance per µm [W m⁻² sr⁻¹ µm⁻¹] from radiance per cm⁻¹ [W m⁻² sr⁻¹ (cm⁻¹)⁻¹]; otherwise the band average is in the input's units.	False

Returns:

Type	Description
BandInstrument	The instrument, (G,) → (B,).

Raises:

Type	Description
ValueError	If the SRF is on a different grid.

Examples:

A flat spectrum per µm, \(L_{\tilde\nu} = c \cdot 10^4 / \tilde\nu^2\), reads c in every band:

>>> import numpy as np
>>> from radtrax.instrument import SpectralResponseFunction
>>> nu = np.linspace(800.0, 1000.0, 2001)
>>> srf = SpectralResponseFunction(
...     wavelengths_hr_nm=1e7 / nu[::-1],
...     band_centers_nm=np.array([10900.0, 11500.0]),
...     band_widths_nm=np.array([900.0, 900.0]),
...     band_names=("B10", "B11"),
... )
>>> band = BandInstrument.from_srf(srf, nu, per_um=True)
>>> [round(float(y), 12) for y in band(7.5 * 1e4 / nu**2)]
[7.5, 7.5]

Source code in src/radtrax/instrument/spectral.py

@classmethod
def from_srf(
    cls,
    srf: SpectralResponseFunction,
    nu_cm: Float[ArrayLike, " G"],
    *,
    shift_cm: Float[ArrayLike, ""] = 0.0,
    per_um: bool = False,
) -> BandInstrument:
    """From a :class:`~radtrax.instrument.SpectralResponseFunction` built
    on the wavelengths of ``nu_cm``.

    ``srf.wavelengths_hr_nm`` must be ``1e7 / nu_cm`` in ascending
    order (``nu_cm`` reversed), so the response is the SRF's own matrix,
    its columns reversed onto the wavenumber grid.

    With ``per_um=True`` the instrument takes radiance per cm⁻¹ and
    returns the band average per µm, which is what a thermal Level-1
    product reports (Landsat TIRS, ECOSTRESS). The SRF's rows ``w`` are
    weights in wavelength (``R Δλ``, summing to 1), so the band average
    of $L_\\lambda = L_{\\tilde\\nu} \\tilde\\nu^2 / 10^4$ is

    $$
    \\bar L_{\\lambda,b} = \\frac{\\int R_b L_\\lambda \\, d\\lambda}{\\int R_b \\, d\\lambda}
        = \\sum_g w_{bg} \\frac{\\tilde\\nu_g^2}{10^4} L_{\\tilde\\nu}(\\tilde\\nu_g),
    $$

    and the factor is folded into the response,
    $f_{bg} = w_{bg} \\tilde\\nu_g^2 / 10^4$, whose rows then sum to the
    band's mean of $\\tilde\\nu^2 / 10^4$ (about 84 at 10.9 µm), not to 1.
    Averaging $L_{\\tilde\\nu}$ with ``w`` and converting afterwards is not
    the same: across a 0.9 µm band at 10.9 µm, $\\tilde\\nu^2$ changes by
    about ±8 %. The factor is taken at the unshifted grid, so a shift
    $\\delta$ adds a relative error of about $2\\delta/\\tilde\\nu$
    (2e-4 for 0.1 cm⁻¹ at 900 cm⁻¹).

    Args:
        srf: The SRF, on ``wavelengths_hr_nm = 1e7 / nu_cm[::-1]``.
        nu_cm: Wavenumber grid [cm⁻¹], strictly increasing, ``(G,)``.
        shift_cm: Spectral shift ``δ`` [cm⁻¹].
        per_um: Return band radiance per µm [W m⁻² sr⁻¹ µm⁻¹] from
            radiance per cm⁻¹ [W m⁻² sr⁻¹ (cm⁻¹)⁻¹]; otherwise the band
            average is in the input's units.

    Returns:
        The instrument, ``(G,) → (B,)``.

    Raises:
        ValueError: If the SRF is on a different grid.

    Examples:
        A flat spectrum per µm, $L_{\\tilde\\nu} = c \\cdot 10^4 / \\tilde\\nu^2$,
        reads ``c`` in every band:

        >>> import numpy as np
        >>> from radtrax.instrument import SpectralResponseFunction
        >>> nu = np.linspace(800.0, 1000.0, 2001)
        >>> srf = SpectralResponseFunction(
        ...     wavelengths_hr_nm=1e7 / nu[::-1],
        ...     band_centers_nm=np.array([10900.0, 11500.0]),
        ...     band_widths_nm=np.array([900.0, 900.0]),
        ...     band_names=("B10", "B11"),
        ... )
        >>> band = BandInstrument.from_srf(srf, nu, per_um=True)
        >>> [round(float(y), 12) for y in band(7.5 * 1e4 / nu**2)]
        [7.5, 7.5]
    """
    nu = np.asarray(nu_cm, dtype=float)
    expected_nm = 1e7 / nu[::-1]
    if srf.wavelengths_hr_nm.shape != expected_nm.shape or not np.allclose(
        srf.wavelengths_hr_nm, expected_nm, rtol=1e-12, atol=0.0
    ):
        raise GridError(
            "BandInstrument.from_srf: build the SRF on wavelengths_hr_nm = "
            "1e7 / nu_cm[::-1] (the wavenumber grid in ascending wavelength)."
        )
    response = srf.matrix[:, ::-1]  # (B, G), rows summing to 1
    if per_um:  # f_bg = w_bg ν_g² / 10⁴, so y_b = Σ_g w_bg L_λ(ν_g)
        response = np.asarray(to_per_um(response, nu))
    return cls(nu, response, shift_cm, srf.band_names)

GaussianInstrument

GaussianInstrument(nu_cm: Float[ArrayLike, ' G'], center_nm: Float[ArrayLike, ' B'], fwhm_nm: Float[ArrayLike, ' B'], *, shift_cm: Float[ArrayLike, ''] = 0.0, shift_nm: Float[ArrayLike, ''] = 0.0, band_names: tuple[str, ...] | None = None)

Bases: Instrument

Gaussian bands built in JAX, shifted analytically in cm⁻¹ or in nm.

y_b is the trapezoidal average of the radiance over band b's response f_b(ν) = exp(−(10⁷/(ν − δ) − δ_λ − c_b)² / 2σ_b²); see the module docstring. Keep each band's response inside the grid: a band shifted off it raises equinox.EquinoxRuntimeError (a RuntimeError).

Attributes:

Name	Type	Description
nu_cm	Float[Array, ' G']	Wavenumber grid [cm⁻¹], strictly increasing, (G,).
center_nm	Float[Array, ' B']	Band centres [nm], (B,).
fwhm_nm	Float[Array, ' B']	Band FWHM [nm], (B,), a differentiable leaf.
shift_cm	Float[Array, '']	Spectral shift δ [cm⁻¹]; positive moves the bands to higher wavenumber, as in :class:BandInstrument.
shift_nm	Float[Array, '']	Spectral shift δ_λ [nm]; positive moves the bands to longer wavelength.
band_names	tuple[str, ...]	One label per band.

Examples:

>>> import jax.numpy as jnp
>>> nu = jnp.linspace(4250.0, 4350.0, 1001)
>>> instrument = GaussianInstrument(
...     nu, jnp.array([2325.0]), jnp.array([8.5])
... )
>>> round(float(instrument(jnp.ones(1001))[0]), 12)  # a flat spectrum
1.0

Source code in src/radtrax/instrument/spectral.py

def __init__(
    self,
    nu_cm: Float[ArrayLike, " G"],
    center_nm: Float[ArrayLike, " B"],
    fwhm_nm: Float[ArrayLike, " B"],
    *,
    shift_cm: Float[ArrayLike, ""] = 0.0,
    shift_nm: Float[ArrayLike, ""] = 0.0,
    band_names: tuple[str, ...] | None = None,
) -> None:
    self.nu_cm = jnp.asarray(nu_cm, dtype=float)
    self.center_nm = jnp.asarray(center_nm, dtype=float)
    self.fwhm_nm = jnp.asarray(fwhm_nm, dtype=float)
    self.shift_cm = jnp.asarray(shift_cm, dtype=float)
    self.shift_nm = jnp.asarray(shift_nm, dtype=float)
    n_bands = self.center_nm.shape[0] if self.center_nm.ndim == 1 else 0
    self.band_names = (
        tuple(f"B{b}" for b in range(n_bands))
        if band_names is None
        else tuple(band_names)
    )

n_bands property

n_bands: int

Number of bands B.

response

response() -> Float[Array, 'B G']

The shifted, trapezoid-weighted response, each row summing to 1, (B, G).

Source code in src/radtrax/instrument/spectral.py

def response(self) -> Float[Array, "B G"]:
    """The shifted, trapezoid-weighted response, each row summing to 1, ``(B, G)``."""
    wavelength_nm = 1e7 / (self.nu_cm - self.shift_cm) - self.shift_nm  # (G,)
    sigma = self.fwhm_nm[:, None] / FWHM_PER_SIGMA  # (B, 1)
    f = jnp.exp(-((wavelength_nm - self.center_nm[:, None]) ** 2) / (2 * sigma**2))
    # Trapezoidal weights on the (unshifted) wavelengths, as the NumPy SRF.
    half = jnp.abs(jnp.diff(1e7 / self.nu_cm)) / 2  # (G − 1,)
    weights = jnp.pad(half, (0, 1)) + jnp.pad(half, (1, 0))  # (G,)
    weighted = f * weights
    total = jnp.sum(weighted, axis=1, keepdims=True)  # (B, 1)
    # A band shifted off the grid has no response left: fail loudly (also
    # under jit), as the NumPy SRF does, rather than return 0/0.
    total = eqx.error_if(
        total,
        jnp.any(total <= 0),
        "GaussianInstrument: a band has no response on the grid.",
    )
    return weighted / total

response_matrix

response_matrix() -> Float[Array, 'B G']

The response matrix A = ∂y/∂L, (B, G): :meth:response.

Source code in src/radtrax/instrument/spectral.py

def response_matrix(self) -> Float[Array, "B G"]:
    """The response matrix ``A = ∂y/∂L``, ``(B, G)``: :meth:`response`."""
    return self.response()

from_srf classmethod

from_srf(srf: SpectralResponseFunction, nu_cm: Float[ArrayLike, ' G'], *, shift_cm: Float[ArrayLike, ''] = 0.0, shift_nm: Float[ArrayLike, ''] = 0.0) -> GaussianInstrument

The bands of a Gaussian :class:~radtrax.instrument.SpectralResponseFunction.

Raises:

Type	Description
ParameterError	If the SRF's bands are not Gaussian.

Source code in src/radtrax/instrument/spectral.py

@classmethod
def from_srf(
    cls,
    srf: SpectralResponseFunction,
    nu_cm: Float[ArrayLike, " G"],
    *,
    shift_cm: Float[ArrayLike, ""] = 0.0,
    shift_nm: Float[ArrayLike, ""] = 0.0,
) -> GaussianInstrument:
    """The bands of a Gaussian :class:`~radtrax.instrument.SpectralResponseFunction`.

    Raises:
        ParameterError: If the SRF's bands are not Gaussian.
    """
    if srf.srf_type != "gaussian":
        raise ParameterError(
            f"GaussianInstrument.from_srf: the SRF must be Gaussian, got {srf.srf_type!r}."
        )
    return cls(
        nu_cm,
        srf.band_centers_nm,
        srf.band_widths_nm,
        shift_cm=shift_cm,
        shift_nm=shift_nm,
        band_names=srf.band_names,
    )

srf

Spectral Response Function — band integration of hyperspectral radiance.

The HAPI LUTs in :mod:radtrax.spectroscopy.hitran tabulate σ(ν, T, P) on a fine line-by-line wavenumber grid. Real instruments observe band-integrated radiance

L_b = ∫ f_b(λ) · L(λ) dλ          / ∫ f_b(λ) dλ

where f_b(λ) is the instrument's Spectral Response Function for band b. This module builds the (n_bands, n_λ) SRF matrix once — each row is a normalised response profile — and applies it as a plain matrix multiplication against a high-resolution radiance spectrum (1-D) or radiance cube (3-D).

The class mirrors the API in the Orbio Project Eucalyptus tutorial (SpectralFilter + get_sensor_spectral_response_function) but exposes a linear-operator view: apply (forward), jacobian (tangent action — identical to forward since SRF is linear), and adjoint (Sᵀ for reverse-mode / VJP use in 3DVar). These three methods match the operator interface used by the observation-operator classes in the original jej_vc_snippets/methane_retrieval/lut_obs_op.py.

Three parametric SRF shapes are supported:

• 'gaussian' : exp[-(λ-λ_c)² / (2 σ²)] with FWHM = 2√(2 ln 2) σ.
• 'rectangular' : 1 inside [λ_c − W/2, λ_c + W/2], 0 outside.
• 'triangular' : linear ramp peaking at λ_c, zero at the edges.

Plus 'custom' which accepts an arbitrary (n_bands, n_λ) matrix (each row gets normalised to sum to 1).

Ported from plumax 0.1.2 at 28f257e (src/plumax/radtran/srf.py); only import paths changed.

SpectralResponseFunction dataclass

SpectralResponseFunction(wavelengths_hr_nm: ndarray, band_centers_nm: ndarray, band_widths_nm: ndarray, band_names: tuple[str, ...] = ('B0',), srf_type: SRFType = 'gaussian', custom_srfs: ndarray | None = None)

Band-integration operator for a multispectral / hyperspectral instrument.

Attributes

wavelengths_hr_nm : np.ndarray High-resolution wavelength grid [nm], shape (n_lambda,). Must be monotone increasing. band_centers_nm : np.ndarray Band centre wavelengths [nm], shape (n_bands,). band_widths_nm : np.ndarray Band FWHM (for 'gaussian'), full-width (for 'rectangular'/'triangular') [nm], shape (n_bands,). band_names : tuple of str Band labels, length n_bands. srf_type : str 'gaussian' | 'rectangular' | 'triangular' | 'custom'. custom_srfs : np.ndarray or None Used only when srf_type == 'custom'. Shape (n_bands, n_lambda).

Notes

Rows of :attr:matrix are L1-normalised so that a flat input spectrum L(λ) = c maps to L_b = c for every band; this matches the "dot-product-in-frequency" convention in the Eucalyptus tutorial.

n_bands property

n_bands: int

Number of bands.

n_lambda property

n_lambda: int

Number of high-resolution wavelength samples.

apply

apply(radiance_hr: Float[ArrayLike, '*batch Ghr']) -> Float[ndarray, '*batch B']

Forward operator: map HR spectrum/cube to band-integrated values.

Parameters

radiance_hr : np.ndarray Input radiance. Shape (n_lambda,) for a single spectrum, or (..., n_lambda) for a cube (spatial + spectral).

Returns

radiance_bands : np.ndarray Band-integrated radiance with the final axis replaced by n_bands.

Source code in src/radtrax/instrument/srf.py

def apply(
    self, radiance_hr: Float[ArrayLike, "*batch Ghr"]
) -> Float[np.ndarray, "*batch B"]:
    """Forward operator: map HR spectrum/cube to band-integrated values.

    Parameters
    ----------
    radiance_hr : np.ndarray
        Input radiance. Shape ``(n_lambda,)`` for a single spectrum, or
        ``(..., n_lambda)`` for a cube (spatial + spectral).

    Returns
    -------
    radiance_bands : np.ndarray
        Band-integrated radiance with the final axis replaced by
        ``n_bands``.
    """
    x = np.asarray(radiance_hr, dtype=float)
    if x.shape[-1] != self.n_lambda:
        raise ShapeError(
            f"SpectralResponseFunction.apply: last axis size {x.shape[-1]} "
            f"must equal n_lambda={self.n_lambda}"
        )
    # ``S @ xᵀ``: contract last axis of `x` (wavelength) with second axis
    # of `matrix`, so `...l, bl -> ...b`.
    return np.einsum("bl, ...l -> ...b", self.matrix, x)

jacobian

jacobian(tangent_hr: Float[ArrayLike, '*batch Ghr']) -> Float[ndarray, '*batch B']

Jacobian-vector product S · v. Linear operator → same as apply.

Source code in src/radtrax/instrument/srf.py

def jacobian(
    self, tangent_hr: Float[ArrayLike, "*batch Ghr"]
) -> Float[np.ndarray, "*batch B"]:
    """Jacobian-vector product ``S · v``. Linear operator → same as ``apply``."""
    return self.apply(tangent_hr)

adjoint

adjoint(cotangent_bands: Float[ArrayLike, '*batch B']) -> Float[ndarray, '*batch Ghr']

Adjoint operator Sᵀ: map band-space cotangents back to HR.

Used for reverse-mode / VJP gradient calculations.

Parameters

cotangent_bands : np.ndarray Band-space cotangent, shape (..., n_bands).

Returns

cotangent_hr : np.ndarray HR-space cotangent, shape (..., n_lambda).

Source code in src/radtrax/instrument/srf.py

def adjoint(
    self, cotangent_bands: Float[ArrayLike, "*batch B"]
) -> Float[np.ndarray, "*batch Ghr"]:
    """Adjoint operator ``Sᵀ``: map band-space cotangents back to HR.

    Used for reverse-mode / VJP gradient calculations.

    Parameters
    ----------
    cotangent_bands : np.ndarray
        Band-space cotangent, shape ``(..., n_bands)``.

    Returns
    -------
    cotangent_hr : np.ndarray
        HR-space cotangent, shape ``(..., n_lambda)``.
    """
    u = np.asarray(cotangent_bands, dtype=float)
    if u.shape[-1] != self.n_bands:
        raise ShapeError(
            f"SpectralResponseFunction.adjoint: last axis size {u.shape[-1]} "
            f"must equal n_bands={self.n_bands}"
        )
    return np.einsum("bl, ...b -> ...l", self.matrix, u)

check_srf_lut_coverage

check_srf_lut_coverage(srf: SpectralResponseFunction, lut_wavelengths_nm: Float[ArrayLike, ' G'], *, context: str, tol: float = 0.01) -> Float[ndarray, ' B']

Warn when SRF bands respond outside the LUT's wavelength coverage.

Resampling a cross-section / target spectrum onto the SRF grid zero-fills wavelengths beyond the LUT range (np.interp(..., left=0, right=0)), silently attenuating any band whose response reaches there — a band whose support lies wholly outside gets an exactly-zero (and misleading) value. This computes each band's fraction of SRF mass outside the LUT coverage and warns for bands exceeding tol.

Parameters

srf Instrument SRF. lut_wavelengths_nm Wavelengths [nm] covered by the LUT (any order); coverage is taken as [min, max]. context Caller name, prepended to the warning message. tol Out-of-coverage SRF-mass fraction above which a band is flagged.

Returns

frac_out : np.ndarray Per-band fraction of SRF response outside the LUT coverage.

Source code in src/radtrax/instrument/srf.py

def check_srf_lut_coverage(
    srf: SpectralResponseFunction,
    lut_wavelengths_nm: Float[ArrayLike, " G"],
    *,
    context: str,
    tol: float = 1e-2,
) -> Float[np.ndarray, " B"]:
    """Warn when SRF bands respond outside the LUT's wavelength coverage.

    Resampling a cross-section / target spectrum onto the SRF grid zero-fills
    wavelengths beyond the LUT range (``np.interp(..., left=0, right=0)``),
    silently attenuating any band whose response reaches there — a band whose
    support lies wholly outside gets an exactly-zero (and misleading) value.
    This computes each band's fraction of SRF mass *outside* the LUT coverage
    and warns for bands exceeding ``tol``.

    Parameters
    ----------
    srf
        Instrument SRF.
    lut_wavelengths_nm
        Wavelengths [nm] covered by the LUT (any order); coverage is taken as
        ``[min, max]``.
    context
        Caller name, prepended to the warning message.
    tol
        Out-of-coverage SRF-mass fraction above which a band is flagged.

    Returns
    -------
    frac_out : np.ndarray
        Per-band fraction of SRF response outside the LUT coverage.
    """
    lam = np.asarray(lut_wavelengths_nm, dtype=float)
    lam_min, lam_max = float(lam.min()), float(lam.max())
    hr = srf.wavelengths_hr_nm
    in_range = (hr >= lam_min) & (hr <= lam_max)
    total = srf.matrix.sum(axis=1)
    covered = srf.matrix[:, in_range].sum(axis=1)
    frac_out = np.where(total > 0.0, 1.0 - covered / total, 0.0)
    for b in np.nonzero(frac_out > tol)[0]:
        warnings.warn(
            f"{context}: band {srf.band_names[b]!r} has {frac_out[b] * 100:.1f}% "
            f"of its SRF response outside the LUT wavelength coverage "
            f"[{lam_min:.1f}, {lam_max:.1f}] nm; its band-integrated value is "
            "silently zero-filled there.",
            stacklevel=3,
        )
    return frac_out

psf

JAX point-spread + ground-sampling instrument operators.

These two operators model the deterministic part of an imaging-spectrometer forward chain that sits between the radiative-transfer output and the band-integration step:

HR radiance cube  ──PSF──►  blurred HR cube  ──GSD──►  LR cube  ──SRF──►  bands

The :class:~radtrax.instrument.srf.SpectralResponseFunction already covers the SRF stage; this module fills in the PSF (spatial blur) and GSD (detector down-sampling). Both operators here are:

• JAX-native — every method is jit-able and vmap-able; gradients flow through apply/jacobian/adjoint via reverse-mode AD without manual bookkeeping. This is the prerequisite for the variational retrieval module :mod:radtrax.retrieval.var3d.
• Linear — jacobian == apply (state-independent), and adjoint implements the Euclidean transpose so the inner-product identity ⟨A x, y⟩ = ⟨x, Aᵀ y⟩ holds to machine precision.
• Cube-shaped (band-last) — inputs are (ny, nx, n_lambda); channels are processed independently. This matches the convention in :mod:radtrax.instrument.srf and plays nicely with jax.vmap along the spectral axis.

The numpy reference for the PSF/GSD logic lives in jej_vc_snippets/methane_retrieval/lut_obs_op.py; this is the JAX rewrite that lets the same operators participate in autodiff-based 3D-Var solves.

Both operators are Equinox modules (radtrax#57): the PSF kernel is a pytree leaf, so jax.grad reaches it, jax.vmap maps over a batch of PSFs, and eqx.filter_jit compiles once per kernel shape rather than once per instance.

Ported from plumax 0.1.2 at 28f257e (src/plumax/radtran/instrument.py), where both were frozen dataclasses.

PointSpreadFunction

PointSpreadFunction(kernel: Real[ArrayLike, 'k k'])

Bases: Module

2D PSF blur applied independently per spectral channel.

The PSF is a fixed convolution kernel (kernel_size × kernel_size, L1-normalised) so the operator is linear and shift-invariant. The Jacobian is the operator itself; the adjoint is convolution with the 180°-flipped kernel (real-valued case).

The kernel is a pytree leaf: differentiate through it with eqx.filter_grad, map over a stack of PSFs with eqx.filter_vmap. Unlike the plumax dataclass, an instance is not hashable (its kernel is an array), so it cannot be a static jit argument or a dict key; pass it to eqx.filter_jit instead.

A kernel that sums to zero raises DegenerateInputError; a traced kernel is checked with eqx.error_if instead, which raises a RuntimeError at run time.

Attributes

kernel : Array L1-normalised 2-D kernel, shape (kernel_size, kernel_size).

Source code in src/radtrax/instrument/psf.py

def __init__(self, kernel: Real[ArrayLike, "k k"]) -> None:
    k = jnp.asarray(kernel, dtype=jnp.float64)
    if k.ndim != 2:
        raise ShapeError(
            f"PointSpreadFunction: kernel must be 2-D, got shape {k.shape}."
        )
    if k.shape[0] != k.shape[1]:
        raise ShapeError(
            f"PointSpreadFunction: kernel must be square, got {k.shape}."
        )
    if k.shape[0] % 2 != 1 or k.shape[0] < 3:
        raise ShapeError(
            f"PointSpreadFunction: kernel side length must be odd ≥ 3 (got {k.shape[0]})."
        )
    s = k.sum()
    try:
        sums_to_zero = float(s) == 0.0
    except jax.errors.ConcretizationTypeError:
        # A traced kernel has no value to check here; fail at run time.
        k = eqx.error_if(k, s == 0.0, "PointSpreadFunction: kernel sums to zero.")
    else:
        if sums_to_zero:
            raise DegenerateInputError(
                "PointSpreadFunction: kernel sums to zero; cannot normalise."
            )
    self.kernel = k / s

normalised_kernel property

normalised_kernel: Float[Array, 'k k']

The kernel divided by its sum, as apply and adjoint use it.

__init__ stores a normalised kernel, but an update that skips __init__ (an optimiser step, eqx.tree_at) need not keep the sum at 1; normalising at use keeps the blur radiance-conserving and makes the kernel's gradient tangent to the Σ K = 1 constraint.

gaussian classmethod

gaussian(fwhm_pixels: float, kernel_size: int = 9) -> PointSpreadFunction

Build an isotropic Gaussian PSF.

FWHM = 2√(2 ln 2) σ; the kernel is sampled on the discrete pixel grid and L1-normalised in __init__.

Source code in src/radtrax/instrument/psf.py

@classmethod
def gaussian(cls, fwhm_pixels: float, kernel_size: int = 9) -> PointSpreadFunction:
    """Build an isotropic Gaussian PSF.

    ``FWHM = 2√(2 ln 2) σ``; the kernel is sampled on the discrete pixel grid
    and L1-normalised in ``__init__``.
    """
    if fwhm_pixels <= 0.0:
        raise ParameterError(
            "PointSpreadFunction.gaussian: `fwhm_pixels` must be > 0."
        )
    c = kernel_size // 2
    y, x = np.ogrid[-c : c + 1, -c : c + 1]
    r2 = x * x + y * y
    sigma = fwhm_pixels / FWHM_PER_SIGMA
    k = np.exp(-r2 / (2.0 * sigma * sigma))
    return cls(kernel=k)

airy classmethod

airy(fwhm_pixels: float, kernel_size: int = 9) -> PointSpreadFunction

Build an Airy-disk-like PSF as a diffraction-ish surrogate.

fwhm_pixels is the kernel's full width at half maximum, matching :meth:gaussian. The Airy argument is scaled so the half-max radius (x = _AIRY_HALF_MAX_X) falls at r = fwhm_pixels / 2; the earlier 1.22·fwhm (first-zero) scaling produced a kernel ≈ 0.63× the requested FWHM.

Source code in src/radtrax/instrument/psf.py

@classmethod
def airy(cls, fwhm_pixels: float, kernel_size: int = 9) -> PointSpreadFunction:
    """Build an Airy-disk-like PSF as a diffraction-ish surrogate.

    ``fwhm_pixels`` is the kernel's full width at half maximum, matching
    :meth:`gaussian`. The Airy argument is scaled so the half-max radius
    (``x = _AIRY_HALF_MAX_X``) falls at ``r = fwhm_pixels / 2``; the earlier
    ``1.22·fwhm`` (first-zero) scaling produced a kernel ≈ 0.63× the
    requested FWHM.
    """
    from scipy.special import j1

    if fwhm_pixels <= 0.0:
        raise ParameterError("PointSpreadFunction.airy: `fwhm_pixels` must be > 0.")
    c = kernel_size // 2
    y, x = np.ogrid[-c : c + 1, -c : c + 1]
    r = np.sqrt(x * x + y * y)
    xarg = np.where(r == 0.0, 1e-12, 2.0 * _AIRY_HALF_MAX_X * r / fwhm_pixels)
    k = (2.0 * j1(xarg) / xarg) ** 2
    return cls(kernel=k)

apply

apply(radiance: HRCube) -> HRCube

Forward blur: y = K * x per channel.

Source code in src/radtrax/instrument/psf.py

def apply(self, radiance: HRCube) -> HRCube:
    """Forward blur: ``y = K * x`` per channel."""
    x = jnp.asarray(radiance)
    if x.ndim != 3:
        raise ShapeError(
            f"PointSpreadFunction.apply: input must be 3-D (ny, nx, n_lambda), "
            f"got shape {x.shape}."
        )
    return _convolve_channels(x, self.normalised_kernel)

jacobian

jacobian(tangent: HRCube) -> HRCube

JVP — same as :meth:apply because the operator is linear.

Source code in src/radtrax/instrument/psf.py

def jacobian(self, tangent: HRCube) -> HRCube:
    """JVP — same as :meth:`apply` because the operator is linear."""
    return self.apply(tangent)

adjoint

adjoint(cotangent: HRCube) -> HRCube

Adjoint: convolution with the flipped kernel.

For a real-valued L1-normalised kernel K the discrete adjoint of the per-channel cross-correlation y = K ⋆ x is Kᵀ ⋆ u = K[::-1] ⋆ u. Symmetric kernels (e.g. Gaussian) are self-adjoint.

Source code in src/radtrax/instrument/psf.py

def adjoint(self, cotangent: HRCube) -> HRCube:
    """Adjoint: convolution with the flipped kernel.

    For a real-valued L1-normalised kernel ``K`` the discrete adjoint of
    the per-channel cross-correlation `y = K ⋆ x` is `Kᵀ ⋆ u = K[::-1] ⋆ u`.
    Symmetric kernels (e.g. Gaussian) are self-adjoint.
    """
    u = jnp.asarray(cotangent)
    if u.ndim != 3:
        raise ShapeError(
            f"PointSpreadFunction.adjoint: input must be 3-D, got shape {u.shape}."
        )
    return _convolve_channels(u, self.normalised_kernel[::-1, ::-1])

GroundSamplingDistance

Bases: Module

Block-mean down-sampling operator for a fixed integer factor.

Maps an HR cube (ny, nx, n_lambda) to LR (ny // f, nx // f, n_lambda) by averaging non-overlapping f × f spatial blocks. If (ny, nx) is not divisible by f, the input is truncated to (ny // f * f, nx // f * f) (top-left crop).

The adjoint distributes each LR value uniformly across the corresponding HR block with weight 1 / f² so that ⟨A x, y⟩ = ⟨x, Aᵀ y⟩ holds exactly under the standard Euclidean inner product.

Attributes

downsample_factor : int Integer block size f. Static: it fixes the output shape, so a new factor is a new compiled function.

Notes

The classmethod :meth:from_optics derives f from the standard pinhole-camera GSD formula GSD [m/px] = sensor_width [mm] · altitude [m] / (focal_length [mm] · image_width [px]) and rounds to the nearest integer. For non-integer ratios you should pre-resample the HR cube — block-mean is only well-defined for integer factors.

from_optics classmethod

from_optics(sensor_width_mm: float, focal_length_mm: float, image_width_px: int, altitude_m: float, pixel_size_hr_m: float) -> GroundSamplingDistance

Build a GSD operator from camera + altitude + HR pixel size.

Dimensional analysis:

GSD [m/px] = (sensor_pixel_size [m/px]) · (altitude [m]) / (focal_length [m])
           = (sensor_width_mm / image_width_px) · altitude_m / focal_length_mm

The mm units on numerator and denominator cancel, so no unit-conversion factor is needed — the numpy reference snippet we ported from had a spurious × 1000 that inflated every realistic camera's GSD by a factor of 1000.

Source code in src/radtrax/instrument/psf.py

@classmethod
def from_optics(
    cls,
    sensor_width_mm: float,
    focal_length_mm: float,
    image_width_px: int,
    altitude_m: float,
    pixel_size_hr_m: float,
) -> GroundSamplingDistance:
    """Build a GSD operator from camera + altitude + HR pixel size.

    Dimensional analysis:

        GSD [m/px] = (sensor_pixel_size [m/px]) · (altitude [m]) / (focal_length [m])
                   = (sensor_width_mm / image_width_px) · altitude_m / focal_length_mm

    The ``mm`` units on numerator and denominator cancel, so no unit-conversion
    factor is needed — the numpy reference snippet we ported from had a spurious
    ``× 1000`` that inflated every realistic camera's GSD by a factor of 1000.
    """
    gsd_m_per_px = (sensor_width_mm * altitude_m) / (
        focal_length_mm * image_width_px
    )
    f = round(gsd_m_per_px / pixel_size_hr_m)
    return cls(downsample_factor=max(1, f))

lr_shape

lr_shape(hr_shape: tuple[int, int, int]) -> tuple[int, int, int]

LR shape implied by an HR shape and the current factor.

Source code in src/radtrax/instrument/psf.py

def lr_shape(self, hr_shape: tuple[int, int, int]) -> tuple[int, int, int]:
    """LR shape implied by an HR shape and the current factor."""
    ny_hr, nx_hr, n_lambda = hr_shape
    f = self.downsample_factor
    return (ny_hr // f, nx_hr // f, n_lambda)

apply

apply(radiance: HRCube) -> LRCube

Forward down-sample by block mean.

Source code in src/radtrax/instrument/psf.py

def apply(self, radiance: HRCube) -> LRCube:
    """Forward down-sample by block mean."""
    x = jnp.asarray(radiance)
    if x.ndim != 3:
        raise ShapeError(
            f"GroundSamplingDistance.apply: input must be 3-D, got shape {x.shape}."
        )
    ny_hr, nx_hr, n_lambda = x.shape
    f = self.downsample_factor
    ny_lr, nx_lr = ny_hr // f, nx_hr // f
    if ny_lr == 0 or nx_lr == 0:
        raise ShapeError(
            f"GroundSamplingDistance.apply: input shape {x.shape[:2]} too small "
            f"for factor {f}."
        )
    x = x[: ny_lr * f, : nx_lr * f, :]
    return x.reshape(ny_lr, f, nx_lr, f, n_lambda).mean(axis=(1, 3))

jacobian

jacobian(tangent: HRCube) -> LRCube

JVP — same as :meth:apply.

Source code in src/radtrax/instrument/psf.py

def jacobian(self, tangent: HRCube) -> LRCube:
    """JVP — same as :meth:`apply`."""
    return self.apply(tangent)

adjoint

adjoint(cotangent: LRCube, hr_shape: tuple[int, int, int]) -> HRCube

Adjoint: spread each LR value uniformly across its HR block.

Parameters

cotangent : Array LR cotangent, shape (ny_lr, nx_lr, n_lambda). hr_shape : tuple Target HR shape (ny_hr, nx_hr, n_lambda). Truncated entries from :meth:apply (when ny_hr or nx_hr are not multiples of f) are filled with zeros — the adjoint of the truncation is the inclusion (... → ny_lr·f rows).

Source code in src/radtrax/instrument/psf.py

def adjoint(self, cotangent: LRCube, hr_shape: tuple[int, int, int]) -> HRCube:
    """Adjoint: spread each LR value uniformly across its HR block.

    Parameters
    ----------
    cotangent : Array
        LR cotangent, shape ``(ny_lr, nx_lr, n_lambda)``.
    hr_shape : tuple
        Target HR shape ``(ny_hr, nx_hr, n_lambda)``. Truncated entries
        from :meth:`apply` (when ``ny_hr`` or ``nx_hr`` are not multiples
        of ``f``) are filled with zeros — the adjoint of the truncation
        is the inclusion ``(... → ny_lr·f rows)``.
    """
    u = jnp.asarray(cotangent)
    if u.ndim != 3:
        raise ShapeError(
            f"GroundSamplingDistance.adjoint: input must be 3-D, got shape {u.shape}."
        )
    ny_hr, nx_hr, n_lambda = hr_shape
    f = self.downsample_factor
    ny_lr, nx_lr = ny_hr // f, nx_hr // f
    if u.shape != (ny_lr, nx_lr, n_lambda):
        raise ShapeError(
            f"GroundSamplingDistance.adjoint: cotangent shape {u.shape} does not "
            f"match expected {(ny_lr, nx_lr, n_lambda)} for hr_shape={hr_shape}."
        )
    # Spread + scale: each LR value contributes 1/f² to every HR pixel
    # in its block (the gradient of a mean reduction).
    block = jnp.repeat(jnp.repeat(u, f, axis=0), f, axis=1) / (f * f)
    if ny_lr * f == ny_hr and nx_lr * f == nx_hr:
        return block
    # Pad with zeros along the truncated rows/cols.
    pad_y = ny_hr - ny_lr * f
    pad_x = nx_hr - nx_lr * f
    return jnp.pad(block, ((0, pad_y), (0, pad_x), (0, 0)))

noise

The NoiseModel protocol: the measurement-noise covariance S_ε.

Optimal estimation needs S_ε as a structured operator, and synthetic scenes need draws from it. A noise model returns S_ε(y) for the band values y it is evaluated at (signal-dependent noise needs y), as a lineax operator tagged positive semidefinite (so gx.solve, gx.logdet and lx.CG all take it), draws ε ~ 𝒩(0, S_ε(y)) from a jax.random key the caller passes, and scores an observation with log 𝒩(y_obs | y, S_ε(y)).

Three diagonal models:

• DiagonalNoise: a fixed per-band standard deviation, S_ε = diag(σ²).
• SNRNoise: a signal-to-noise ratio with a noise floor, σ_b = √((y_b / SNR_b)² + σ_floor,b²).
• NEdTNoise: a thermal band's noise-equivalent temperature difference turned into radiance, σ_b = NEΔT_b · ∂B̄_b/∂T at a reference temperature or at the scene's band brightness temperature (radtrax#249).

Shapes: B bands.

NoiseModel

Bases: Module

Protocol: the measurement-noise covariance S_ε(y) and draws from it.

An abstract Equinox module, so filter_jit compiles each concrete model and a retrieval can take any of them. A subclass implements covariance; sample draws through gx.sample_mvn, which dispatches on the operator's structure (√σ² · z for a diagonal).

covariance abstractmethod

covariance(y: Float[ArrayLike, ' B']) -> AbstractLinearOperator

S_ε at the band values y, a positive-definite (B, B) operator.

Parameters:

Name	Type	Description	Default
y	Float[ArrayLike, ' B']	Band values the noise is evaluated at, (B,).	required

Source code in src/radtrax/instrument/noise.py

@abc.abstractmethod
def covariance(self, y: Float[ArrayLike, " B"]) -> lx.AbstractLinearOperator:
    """``S_ε`` at the band values ``y``, a positive-definite ``(B, B)`` operator.

    Args:
        y: Band values the noise is evaluated at, ``(B,)``.
    """

sample

sample(key: PRNGKeyArray, y: Float[ArrayLike, ' B']) -> Float[Array, ' B']

One draw ε ~ 𝒩(0, S_ε(y)); add it to y for a noisy observation.

Parameters:

Name	Type	Description	Default
key	PRNGKeyArray	A jax.random key; the same key gives the same draw.	required
y	Float[ArrayLike, ' B']	Band values the noise is evaluated at, (B,).	required

Returns:

Type	Description
Float[Array, ' B']	The noise ε, (B,).

Source code in src/radtrax/instrument/noise.py

def sample(
    self, key: PRNGKeyArray, y: Float[ArrayLike, " B"]
) -> Float[Array, " B"]:
    """One draw ``ε ~ 𝒩(0, S_ε(y))``; add it to ``y`` for a noisy observation.

    Args:
        key: A ``jax.random`` key; the same key gives the same draw.
        y: Band values the noise is evaluated at, ``(B,)``.

    Returns:
        The noise ``ε``, ``(B,)``.
    """
    y = jnp.asarray(y, dtype=float)
    return gx.sample_mvn(jnp.zeros_like(y), self.covariance(y), key=key)[0]

log_prob

log_prob(y_obs: Float[ArrayLike, ' B'], y_model: Float[ArrayLike, ' B']) -> Float[Array, '']

log 𝒩(y_obs | y_model, S_ε(y_model)), through gx.gaussian_log_prob.

The noise is evaluated at the modelled values, so for signal-dependent noise the log |S_ε| term depends on them too and is kept: dropping it would favour brighter solutions.

Parameters:

Name	Type	Description	Default
y_obs	Float[ArrayLike, ' B']	Observed band values, (B,).	required
y_model	Float[ArrayLike, ' B']	Modelled band values, (B,).	required

Returns:

Type	Description
Float[Array, '']	The log-density, a scalar.

Source code in src/radtrax/instrument/noise.py

def log_prob(
    self, y_obs: Float[ArrayLike, " B"], y_model: Float[ArrayLike, " B"]
) -> Float[Array, ""]:
    """``log 𝒩(y_obs | y_model, S_ε(y_model))``, through ``gx.gaussian_log_prob``.

    The noise is evaluated at the modelled values, so for signal-dependent
    noise the ``log |S_ε|`` term depends on them too and is kept: dropping
    it would favour brighter solutions.

    Args:
        y_obs: Observed band values, ``(B,)``.
        y_model: Modelled band values, ``(B,)``.

    Returns:
        The log-density, a scalar.
    """
    y_model = jnp.asarray(y_model, dtype=float)
    return gx.gaussian_log_prob(
        y_model, self.covariance(y_model), jnp.asarray(y_obs, dtype=float)
    )

DiagonalNoise

DiagonalNoise(sigma: Real[ArrayLike, ' B'])

Bases: NoiseModel

Independent noise with a fixed standard deviation per band, S_ε = diag(σ²).

Attributes:

Name	Type	Description
sigma	Float[Array, ' B']	Standard deviation per band, in the units of y, (B,); each strictly positive.

Examples:

>>> import jax.numpy as jnp
>>> noise = DiagonalNoise(jnp.array([0.1, 0.2]))
>>> import gaussx as gx
>>> [round(v, 6) for v in gx.diag(noise.covariance(jnp.ones(2))).tolist()]
[0.01, 0.04]

Source code in src/radtrax/instrument/noise.py

def __init__(self, sigma: Real[ArrayLike, " B"]) -> None:
    self.sigma = jnp.asarray(sigma, dtype=float)
    _check_bands("DiagonalNoise", sigma=self.sigma)
    _check_positive("DiagonalNoise", sigma=self.sigma)

covariance

covariance(y: Float[ArrayLike, ' B']) -> AbstractLinearOperator

diag(σ²), whatever y; y sets only the band count.

Source code in src/radtrax/instrument/noise.py

def covariance(self, y: Float[ArrayLike, " B"]) -> lx.AbstractLinearOperator:
    """``diag(σ²)``, whatever ``y``; ``y`` sets only the band count."""
    _check_y("DiagonalNoise", y, self.sigma)
    return _diagonal_covariance(self.sigma**2)

SNRNoise

SNRNoise(snr: Real[ArrayLike, ' B'], floor: Real[ArrayLike, ' B'])

Bases: NoiseModel

Signal-dependent noise from a signal-to-noise ratio and a floor.

σ_b = √((y_b / SNR_b)² + σ_floor,b²), S_ε = diag(σ²): the floor keeps S_ε positive definite where the signal vanishes.

Attributes:

Name	Type	Description
snr	Float[Array, ' B']	Signal-to-noise ratio per band (dimensionless), (B,); each strictly positive.
floor	Float[Array, ' B']	Noise floor per band, in the units of y, (B,); each strictly positive.

Examples:

>>> import jax.numpy as jnp
>>> noise = SNRNoise(snr=jnp.array([100.0]), floor=jnp.array([0.03]))
>>> import gaussx as gx
>>> round(float(gx.diag(noise.covariance(jnp.array([4.0])))[0]), 6)
0.0025

Source code in src/radtrax/instrument/noise.py

def __init__(
    self, snr: Real[ArrayLike, " B"], floor: Real[ArrayLike, " B"]
) -> None:
    self.snr = jnp.asarray(snr, dtype=float)
    self.floor = jnp.asarray(floor, dtype=float)
    _check_bands("SNRNoise", snr=self.snr, floor=self.floor)
    _check_positive("SNRNoise", snr=self.snr, floor=self.floor)

covariance

covariance(y: Float[ArrayLike, ' B']) -> AbstractLinearOperator

diag((y / SNR)² + σ_floor²), differentiable in y.

Source code in src/radtrax/instrument/noise.py

def covariance(self, y: Float[ArrayLike, " B"]) -> lx.AbstractLinearOperator:
    """``diag((y / SNR)² + σ_floor²)``, differentiable in ``y``."""
    y = _check_y("SNRNoise", y, self.snr)
    return _diagonal_covariance((y / self.snr) ** 2 + self.floor**2)

NEdTNoise

NEdTNoise(netd_k: Real[ArrayLike, ' B'], instrument: Instrument, *, t_ref_k: Real[ArrayLike, ''] = 300.0, mode: Literal['reference', 'scene'] = 'reference')

Bases: NoiseModel

Thermal-band noise from a noise-equivalent temperature difference.

\[ S_ε = \mathrm{diag}(σ^2), \qquad σ_b = \mathrm{NEΔT}_b \, \frac{\partial \bar B_b}{\partial T}\Big|_{T}, \qquad \bar B_b(T) = \sum_g A_{bg} B(\tilde\nu_g, T), \]

with \(A\) the instrument's response matrix and \(B\) the Planck radiance (radtrax.radiometry.band_planck_radiance). Two conventions for \(T\):

• "reference" (the default): \(T = T_\mathrm{ref}\), the temperature the specification quotes NEΔT at (300 K for TIRS-2 and ECOSTRESS). \(σ\) is constant in y, as for a detector-noise-limited instrument; \(\partial \bar B / \partial T\) at a scalar \(T_\mathrm{ref}\) is one JVP.
• "scene": \(T = T_b(y)\), each band's brightness temperature of y (radtrax.radiometry.band_brightness_temperature), for "NEΔT at the scene temperature" taken literally. \(σ\) depends on y, so log_prob keeps \(\log |S_ε|\) (the NoiseModel contract), and a band with \(y_b \le 0\) has no brightness temperature: its \(σ\) is NaN.

\(σ\) is recomputed in every call (not cached at construction), so gradients reach the instrument's parameters; it is in the instrument's radiance units (W m⁻² sr⁻¹ (cm⁻¹)⁻¹ for rows that sum to 1).

Attributes:

Name	Type	Description
netd_k	Float[Array, ' B']	NEΔT per band [K], (B,); each strictly positive.
instrument	Instrument	The bands \(\bar B\) is integrated over.
t_ref_k	Float[Array, '']	Reference temperature [K], a scalar > 0 (used in the "reference" mode).
mode	Literal['reference', 'scene']	"reference" or "scene".

Examples:

>>> import jax.numpy as jnp
>>> from radtrax.instrument import GaussianInstrument
>>> nu = jnp.linspace(916.0, 919.5, 3501)  # a 1 nm band at 10.9 µm
>>> band = GaussianInstrument(nu, jnp.array([10900.0]), jnp.array([1.0]))
>>> noise = NEdTNoise(jnp.array([0.048]), band)
>>> f"{float(noise.radiance_sigma(jnp.zeros(1))[0]):.4e}"
'8.1485e-05'

Source code in src/radtrax/instrument/noise.py

def __init__(
    self,
    netd_k: Real[ArrayLike, " B"],
    instrument: Instrument,
    *,
    t_ref_k: Real[ArrayLike, ""] = 300.0,
    mode: Literal["reference", "scene"] = "reference",
) -> None:
    self.netd_k = jnp.asarray(netd_k, dtype=float)
    self.instrument = instrument
    self.t_ref_k = jnp.asarray(t_ref_k, dtype=float)
    self.mode = mode
    if mode not in ("reference", "scene"):
        raise ParameterError(
            f"NEdTNoise: mode must be 'reference' or 'scene', got {mode!r}."
        )
    _check_bands("NEdTNoise", netd_k=self.netd_k)
    if self.t_ref_k.ndim != 0:
        raise ShapeError(
            f"NEdTNoise: t_ref_k must be a scalar, got {self.t_ref_k.shape}."
        )
    _check_positive("NEdTNoise", netd_k=self.netd_k)
    _check_positive("NEdTNoise", t_ref_k=self.t_ref_k[None])
    n_bands = jax.eval_shape(instrument, instrument.nu_cm).shape
    if n_bands != self.netd_k.shape:
        raise ShapeError(
            f"NEdTNoise: netd_k is {self.netd_k.shape} for an instrument "
            f"with {n_bands} bands."
        )

radiance_sigma

radiance_sigma(y: Float[ArrayLike, ' B']) -> Float[Array, ' B']

\(σ_b = \mathrm{NEΔT}_b \, \partial \bar B_b / \partial T\), in the instrument's units.

Parameters:

Name	Type	Description	Default
y	Float[ArrayLike, ' B']	Band values the noise is evaluated at, (B,); read only in the "scene" mode.	required

Returns:

Type	Description
Float[Array, ' B']	\(σ\), (B,).

Source code in src/radtrax/instrument/noise.py

def radiance_sigma(self, y: Float[ArrayLike, " B"]) -> Float[Array, " B"]:
    """$σ_b = \\mathrm{NEΔT}_b \\, \\partial \\bar B_b / \\partial T$, in the instrument's units.

    Args:
        y: Band values the noise is evaluated at, ``(B,)``; read only in
            the ``"scene"`` mode.

    Returns:
        $σ$, ``(B,)``.
    """
    y = _check_y("NEdTNoise", y, self.netd_k)
    return self.netd_k * self._slope(self._temperature(y))

netd_at

netd_at(y: Float[ArrayLike, ' B']) -> Float[Array, ' B']

The NEΔT this noise amounts to at the scene: \(σ_b / \partial \bar B_b / \partial T|_{T_b(y)}\) [K].

The same radiance noise is a different NEΔT at another scene temperature: in the "reference" mode it equals netd_k only at \(y = \bar B(T_\mathrm{ref})\); in the "scene" mode it is netd_k wherever \(y > 0\).

Parameters:

Name	Type	Description	Default
y	Float[ArrayLike, ' B']	Band values, (B,), each > 0.	required

Returns:

Type	Description
Float[Array, ' B']	NEΔT per band [K], (B,).

Source code in src/radtrax/instrument/noise.py

def netd_at(self, y: Float[ArrayLike, " B"]) -> Float[Array, " B"]:
    """The NEΔT this noise amounts to at the scene: $σ_b / \\partial \\bar B_b / \\partial T|_{T_b(y)}$ [K].

    The same radiance noise is a different NEΔT at another scene
    temperature: in the ``"reference"`` mode it equals ``netd_k`` only at
    $y = \\bar B(T_\\mathrm{ref})$; in the ``"scene"`` mode it is
    ``netd_k`` wherever $y > 0$.

    Args:
        y: Band values, ``(B,)``, each > 0.

    Returns:
        NEΔT per band [K], ``(B,)``.
    """
    y = _check_y("NEdTNoise", y, self.netd_k)
    t_b = band_brightness_temperature(self.instrument, y)
    return self.radiance_sigma(y) / self._slope(t_b)

covariance

covariance(y: Float[ArrayLike, ' B']) -> AbstractLinearOperator

diag(σ²) with \(σ\) from radiance_sigma, differentiable in the instrument.

Source code in src/radtrax/instrument/noise.py

def covariance(self, y: Float[ArrayLike, " B"]) -> lx.AbstractLinearOperator:
    """``diag(σ²)`` with $σ$ from `radiance_sigma`, differentiable in the instrument."""
    return _diagonal_covariance(self.radiance_sigma(y) ** 2)