Geo

Sensor data and GeoTensor operators, behind the [geo] extra (pip install "radtrax[geo]"). Today: Sentinel-2 SRFs from the agency table georeader bundles, integrated by radtrax's SpectralResponseFunction; a Sentinel-2 L1C image read by georeader as a BandInstrument, a per-pixel ViewingGeometry and TOA radiance; Landsat 8/9 OLI B6 and B7 and TIRS B10 and B11 from the USGS spectral responses radtrax bundles, with radiance from Level-1 digital numbers and the MTL rescaling and a per-pixel ViewingGeometry from the Collection 2 angle bands (no extra needed); an EMIT L1B image read by georeader's EMITImage as a BandInstrument with Gaussian bands of its FWHM and a per-pixel ViewingGeometry from its obs file; and the matched filter as pipekit operators on GeoTensors, around the one core in radtrax.retrieval.matched_filter, with targets from the single-slab observation model or, through UnitAbsorptionTarget, from v0's ForwardModel, and MBMP, the Sentinel-2 multi-band multi-pass signal (geotoolz's former plume.SBMP).

Landsat's thermal bands report radiance per µm, while the forward model works per cm⁻¹: landsat_instrument(satellite, nu_cm, bands=("B10", "B11"), per_um=True) takes radiance per cm⁻¹ and returns the band average per µm, with the conversion inside the band integral. landsat_noise turns the median on-orbit NEΔT at 300 K (LANDSAT_NETD_K: 0.048 K and 0.054 K for TIRS, 0.048 K and 0.066 K for TIRS-2) into radiance noise through that instrument; LANDSAT_CALIBRATION_K holds TIRS-2's vicarious calibration residual (0.16 K and 0.61 K). Both tables come from J. Barsi's presentation at the ESA Living Planet Symposium 2022 (session H2-02), as quoted in radtrax#225. The TIRS responses are the USGS band-average RSR workbooks (L8_TIRS_RSR.xlsx, L9_TIRS2_RSR.xlsx), recorded with their SHA-256 in the bundled table's header.

ECOSTRESS's five thermal bands, B1 to B5 at about 8.3, 8.8, 9.2, 10.5 and 12.1 µm, follow the same pattern from plain inputs, no extra needed: ecostress_instrument(nu_cm, bands=...) integrates the science team's spectral responses (revision 3, as the EUMETSAT NWP SAF redistributes them, bundled with their SHA-256 in the table's header) and returns band radiance per µm by default, the L1B unit. Bands 1 and 3 hold fill values from 15 May 2019 until five-band mode resumed on 28 April 2023 (ECOSTRESS_3BAND_PERIOD), and the two modes have been interleaved since, so ecostress_bands reads a granule's band set from its BandSpecification metadata (the date alone decides only before 28 April 2023); the three-band set is ECOSTRESS_3BAND, bands 2, 4 and 5. ecostress_noise uses ECOSTRESS_NETD_K, 0.1 K per band at 300 K: that is the requirement (NEΔT < 0.1 K, ECOSTRESS L2 LSTE ATBD, 2018), an upper bound until per-band on-orbit values are sourced. ecostress_radiance stacks the L1B radiance_<n> datasets with their fill values as NaN, and ecostress_geometry builds a per-pixel ViewingGeometry from the L1B geolocation angles.

geo

Sensor data and GeoTensor operators, behind the [geo] extra.

• :func:sentinel2_srf reads the Sentinel-2 SRF table georeader bundles into radtrax's SRF (radtrax#28); :func:sentinel2_band_responses gives its raw B11 and B12 responses on a wavelength grid (radtrax#115).
• :func:sentinel2_instrument, :func:sentinel2_geometry and :func:sentinel2_radiance (:mod:radtrax.geo.sentinel2): a Sentinel-2 L1C image read by georeader as a BandInstrument, a per-pixel ViewingGeometry and TOA radiance (radtrax#120).
• :func:landsat_srf, :func:landsat_instrument, :func:landsat_radiance and :func:landsat_geometry (:mod:radtrax.geo.landsat): Landsat 8/9 OLI B6 and B7 from the USGS spectral responses radtrax bundles, radiance from Level-1 digital numbers and the MTL rescaling, and a per-pixel ViewingGeometry from the Collection 2 angle bands (radtrax#121); the TIRS bands B10 and B11 with band radiance per µm, and their NEΔT as a noise model, :func:landsat_noise (radtrax#225); no extra needed.
• :func:ecostress_bands, :func:ecostress_srf, :func:ecostress_instrument, :func:ecostress_noise, :func:ecostress_radiance and :func:ecostress_geometry (:mod:radtrax.geo.ecostress): ECOSTRESS's five thermal bands and its three-band mode, from the science team's spectral responses radtrax bundles, band radiance per µm, NEΔT as a noise model, L1B radiance and a per-pixel ViewingGeometry from the geolocation angles (radtrax#226); no extra needed.
• :func:emit_instrument and :func:emit_geometry (:mod:radtrax.geo.emit): an EMIT L1B image read by georeader's EMITImage as a BandInstrument with Gaussian bands of its FWHM, and a per-pixel ViewingGeometry from its obs file (radtrax#114).
• :mod:radtrax.geo.matched_filter: pipekit operators on GeoTensors around the one matched-filter core (radtrax#29), ported from geotoolz; :class:UnitAbsorptionTarget gives them v0's ForwardModel target (radtrax#98).
• :class:MBMP (:mod:radtrax.geo.mbmp): the Sentinel-2 multi-band multi-pass signal, geotoolz's former plume.SBMP (radtrax#30), with optional per-pass scale factors; :class:SBMP and :class:MBSP complete radtrax · 2's multi-pass family (radtrax#111).

Importing this subpackage does not need the extra; the operators load on first use, and a missing extra raises an ImportError naming it.

srf

Sentinel-2 spectral response functions from georeader's bundled table.

georeader keeps the agency SRF tables as sensor data (COPE-GSEG-EOPG-TN-15-0007 v4.0, bundled, so no network); radtrax owns the band integration (:class:radtrax.instrument.SpectralResponseFunction). This module reads the one into the other (radtrax#28). It needs the [geo] extra.

sentinel2_srf

sentinel2_srf(satellite: str = 'S2A', *, wavelengths_hr_nm: Sequence[float] | Float[ArrayLike, ' Ghr'] | None = None, bands: Sequence[str] | None = None) -> SpectralResponseFunction

Build a Sentinel-2 SRF from the agency table georeader bundles.

georeader drops the table's rows where every band is at most 1e-6, so the table is first restored onto its full 1 nm grid with zeros there; interpolating across a dropped gap would invent a response.

The 1 nm responses are then projected onto wavelengths_hr_nm, not sampled at it. With φ_j the linear-interpolation (hat) basis on the grid, band b's weight on grid point j is W[b, j] = Σ_k R_b(λ_k) φ_j(λ_k) over the table's wavelengths. Band integration Σ_j W[b, j] L(λ_j) / Σ_j W[b, j] is then exact for any radiance that is linear between grid points, and every band keeps its full response on a coarse grid. Sampling instead loses most of a narrow band at 5–10 nm, and a band can vanish entirely from about 25 nm. On the table's own 1 nm grid the two coincide.

Parameters:

Name	Type	Description	Default
satellite	str	"S2A", "S2B" or "S2C".	'S2A'
wavelengths_hr_nm	Sequence[float] | Float[ArrayLike, ' Ghr'] | None	The high-resolution grid [nm] to integrate on, strictly increasing. Defaults to the table's 1 nm grid.	None
bands	Sequence[str] | None	Band names to keep, in order (e.g. ("B11", "B12")). Defaults to all 13.	None

Returns:

Type	Description
SpectralResponseFunction	A custom SRF. Its band centres are the response-weighted mean
SpectralResponseFunction	wavelengths and its widths the equivalent widths ∫R dλ / max R,
SpectralResponseFunction	both computed on the 1 nm table. Its custom_srfs hold the
SpectralResponseFunction	projected weights divided by the grid's trapezoid bin widths, which
SpectralResponseFunction	SpectralResponseFunction multiplies back before normalising.

Warns:

Type	Description
UserWarning	If the grid misses more than 0.1 % of a band's response; the band is then normalised over the part it covers, which biases its value.

Raises:

Type	Description
ImportError	If the [geo] extra is not installed.
KeyError	If a requested band is not in the table.

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

def sentinel2_srf(
    satellite: str = "S2A",
    *,
    wavelengths_hr_nm: Sequence[float] | Float[ArrayLike, " Ghr"] | None = None,
    bands: Sequence[str] | None = None,
) -> SpectralResponseFunction:
    """Build a Sentinel-2 SRF from the agency table georeader bundles.

    georeader drops the table's rows where every band is at most 1e-6, so
    the table is first restored onto its full 1 nm grid with zeros there;
    interpolating across a dropped gap would invent a response.

    The 1 nm responses are then *projected* onto ``wavelengths_hr_nm``, not
    sampled at it. With ``φ_j`` the linear-interpolation (hat) basis on the
    grid, band ``b``'s weight on grid point ``j`` is
    ``W[b, j] = Σ_k R_b(λ_k) φ_j(λ_k)`` over the table's wavelengths. Band
    integration ``Σ_j W[b, j] L(λ_j) / Σ_j W[b, j]`` is then exact for any
    radiance that is linear between grid points, and every band keeps its full
    response on a coarse grid. Sampling instead loses most of a narrow band at
    5–10 nm, and a band can vanish entirely from about 25 nm. On the table's
    own 1 nm grid the two coincide.

    Args:
        satellite: ``"S2A"``, ``"S2B"`` or ``"S2C"``.
        wavelengths_hr_nm: The high-resolution grid [nm] to integrate on,
            strictly increasing. Defaults to the table's 1 nm grid.
        bands: Band names to keep, in order (e.g. ``("B11", "B12")``).
            Defaults to all 13.

    Returns:
        A ``custom`` SRF. Its band centres are the response-weighted mean
        wavelengths and its widths the equivalent widths ``∫R dλ / max R``,
        both computed on the 1 nm table. Its ``custom_srfs`` hold the
        projected weights divided by the grid's trapezoid bin widths, which
        ``SpectralResponseFunction`` multiplies back before normalising.

    Warns:
        UserWarning: If the grid misses more than ``0.1 %`` of a band's
            response; the band is then normalised over the part it covers,
            which biases its value.

    Raises:
        ImportError: If the ``[geo]`` extra is not installed.
        KeyError: If a requested band is not in the table.
    """
    table = _read_srf_table(satellite)
    if bands is not None:
        missing = [b for b in bands if b not in table.columns]
        if missing:
            raise MissingKeyError(f"sentinel2_srf: no band(s) {missing} in {satellite}")
        table = table[list(bands)]
    wl_table = table.index.to_numpy(dtype=float)
    grid = np.arange(wl_table[0], wl_table[-1] + 1.0)  # the full 1 nm grid
    full = table.reindex(grid, fill_value=0.0).to_numpy(dtype=float).T  # (B, N)
    names = tuple(str(c) for c in table.columns)
    return _projected_srf("sentinel2_srf", grid, full, names, wavelengths_hr_nm)

sentinel2_band_responses

sentinel2_band_responses(satellite: str = 'S2A', wavelengths_nm: Float[ArrayLike, ' G'] | None = None) -> tuple[Float[ndarray, ' G'], Float[ndarray, ' G']]

Sentinel-2 band 11 and 12 responses, interpolated onto a wavelength grid.

The raw agency responses (not normalised), zero outside the table, as :func:radtrax.spectroscopy.synthetic.methane_cross_section takes them (radtrax · 1 and · 2).

Parameters:

Name	Type	Description	Default
satellite	str	"S2A", "S2B" or "S2C".	'S2A'
wavelengths_nm	Float[ArrayLike, ' G'] | None	The grid [nm], (G,); by default :data:radtrax.spectroscopy.synthetic.WL_NM.	None

Returns:

Type	Description
tuple[Float[ndarray, ' G'], Float[ndarray, ' G']]	(f11, f12), each (G,).

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

def sentinel2_band_responses(
    satellite: str = "S2A", wavelengths_nm: Float[ArrayLike, " G"] | None = None
) -> tuple[Float[np.ndarray, " G"], Float[np.ndarray, " G"]]:
    """Sentinel-2 band 11 and 12 responses, interpolated onto a wavelength grid.

    The raw agency responses (not normalised), zero outside the table, as
    :func:`radtrax.spectroscopy.synthetic.methane_cross_section` takes them
    (*radtrax · 1* and *· 2*).

    Args:
        satellite: ``"S2A"``, ``"S2B"`` or ``"S2C"``.
        wavelengths_nm: The grid [nm], ``(G,)``; by default
            :data:`radtrax.spectroscopy.synthetic.WL_NM`.

    Returns:
        ``(f11, f12)``, each ``(G,)``.
    """
    if wavelengths_nm is None:
        from radtrax.spectroscopy.synthetic import WL_NM

        wavelengths_nm = WL_NM
    wl = np.asarray(wavelengths_nm, dtype=float)
    table = _read_srf_table(satellite)
    f11 = np.interp(wl, table.index.values, table["B11"].values, left=0, right=0)
    f12 = np.interp(wl, table.index.values, table["B12"].values, left=0, right=0)
    return f11, f12

sentinel2

Sentinel-2 L1C on radtrax's terms: instrument, geometry and radiance.

georeader's S2ImageL1C reads the product: TOA reflectance as digital numbers, the sun and per-band viewing angles on the 5 km tie-point grid of the tile metadata, and the band solar irradiance. This module turns them into the three things a forward model takes (radtrax#120):

• :func:sentinel2_instrument: a :class:~radtrax.instrument.BandInstrument for the SWIR pair B11 and B12 on the model's wavenumber grid, from the agency SRFs (:func:~radtrax.geo.sentinel2_srf);
• :func:sentinel2_geometry: a :class:~radtrax.atmosphere.ViewingGeometry per pixel, its angles interpolated from the tie points onto the image grid;
• :func:sentinel2_radiance: TOA radiance from TOA reflectance,

\[ L_b = \frac{\rho_b\, E_b \cos\theta_s}{\pi\, d^2}, \]

with E_b the band solar irradiance, d the Earth–Sun distance [AU] (the product's U = 1/d²) and θ_s the solar zenith angle of the pixel.

The tie-point grid. The L1C metadata gives each angle at the nodes of a 5 km grid whose first node is the tile's upper-left corner. georeader's sza/saa/vza/vaa GeoTensors place that first value at the centre of a 5 km cell instead, 2.5 km off in both directions; its own anglesULXY is the origin that puts the nodes on the corner. The angles are re-georeferenced on anglesULXY here, so bilinear interpolation onto the image grid uses the nodes where the metadata puts them, and a tile's pixels all lie inside the node grid (23 nodes span 110 km, a tile 109.8 km).

The functions need the [geo] extra; this module imports without it.

sentinel2_instrument

sentinel2_instrument(satellite: str, nu_cm: Float[ArrayLike, ' G'], *, bands: Sequence[str] = ('B11', 'B12'), shift_cm: float = 0.0) -> BandInstrument

Sentinel-2 bands as a :class:~radtrax.instrument.BandInstrument.

The agency SRFs (:func:~radtrax.geo.sentinel2_srf) projected onto the wavelengths of the model grid, λ_g = 10⁷ / ν_g, and turned onto the wavenumber grid by :meth:BandInstrument.from_srf <radtrax.instrument.BandInstrument.from_srf>:

\[ y_b = \sum_g f_{bg}\, L(\nu_g + \delta), \]

with f the SRF's row-normalised matrix, columns reversed. At zero shift it is :meth:SpectralResponseFunction.apply <radtrax.instrument.SpectralResponseFunction.apply> of the same SRF on the spectrum in ascending wavelength.

Parameters:

Name	Type	Description	Default
satellite	str	"S2A", "S2B" or "S2C".	required
nu_cm	Float[ArrayLike, ' G']	The model's wavenumber grid [cm⁻¹], strictly increasing, (G,); it should cover each band's response.	required
bands	Sequence[str]	Band names, in order.	('B11', 'B12')
shift_cm	float	Spectral shift δ [cm⁻¹], a differentiable leaf.	0.0

Returns:

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

Warns:

Type	Description
UserWarning	If the grid misses more than 0.1 % of a band's response (from :func:~radtrax.geo.sentinel2_srf).

Raises:

Type	Description
ImportError	If the [geo] extra is not installed.
KeyError	If a band is not in the SRF table.
ValueError	If nu_cm is not 1-D and strictly increasing.

Examples:

>>> import numpy as np
>>> nu = np.arange(4150.0, 6700.0)  # 1493–2410 nm, B11 and B12
>>> instrument = sentinel2_instrument("S2A", nu)
>>> instrument.band_names
('B11', 'B12')
>>> [round(float(y), 12) for y in instrument(np.ones(nu.size))]
[1.0, 1.0]

Source code in src/radtrax/geo/sentinel2.py

def sentinel2_instrument(
    satellite: str,
    nu_cm: Float[ArrayLike, " G"],
    *,
    bands: Sequence[str] = ("B11", "B12"),
    shift_cm: float = 0.0,
) -> BandInstrument:
    """Sentinel-2 bands as a :class:`~radtrax.instrument.BandInstrument`.

    The agency SRFs (:func:`~radtrax.geo.sentinel2_srf`) projected onto the
    wavelengths of the model grid, ``λ_g = 10⁷ / ν_g``, and turned onto the
    wavenumber grid by :meth:`BandInstrument.from_srf
    <radtrax.instrument.BandInstrument.from_srf>`:

    $$
    y_b = \\sum_g f_{bg}\\, L(\\nu_g + \\delta),
    $$

    with ``f`` the SRF's row-normalised matrix, columns reversed. At zero
    shift it is :meth:`SpectralResponseFunction.apply
    <radtrax.instrument.SpectralResponseFunction.apply>` of the same SRF on
    the spectrum in ascending wavelength.

    Args:
        satellite: ``"S2A"``, ``"S2B"`` or ``"S2C"``.
        nu_cm: The model's wavenumber grid [cm⁻¹], strictly increasing,
            ``(G,)``; it should cover each band's response.
        bands: Band names, in order.
        shift_cm: Spectral shift ``δ`` [cm⁻¹], a differentiable leaf.

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

    Warns:
        UserWarning: If the grid misses more than ``0.1 %`` of a band's
            response (from :func:`~radtrax.geo.sentinel2_srf`).

    Raises:
        ImportError: If the ``[geo]`` extra is not installed.
        KeyError: If a band is not in the SRF table.
        ValueError: If ``nu_cm`` is not 1-D and strictly increasing.

    Examples:
        >>> import numpy as np
        >>> nu = np.arange(4150.0, 6700.0)  # 1493–2410 nm, B11 and B12
        >>> instrument = sentinel2_instrument("S2A", nu)
        >>> instrument.band_names
        ('B11', 'B12')
        >>> [round(float(y), 12) for y in instrument(np.ones(nu.size))]
        [1.0, 1.0]
    """
    nu = np.asarray(nu_cm, dtype=float)
    if nu.ndim != 1 or nu.size < 2 or np.any(np.diff(nu) <= 0.0):
        raise GridError(
            "sentinel2_instrument: `nu_cm` must be 1-D, strictly increasing, "
            "with at least two points."
        )
    srf = sentinel2_srf(satellite, wavelengths_hr_nm=1e7 / nu[::-1], bands=bands)
    return BandInstrument.from_srf(srf, nu, shift_cm=shift_cm)

sentinel2_geometry

sentinel2_geometry(image: Any, band: str = 'B12') -> ViewingGeometry

Per-pixel viewing geometry of an L1C image, on its grid (H, W).

The solar zenith θ_s, the viewing zenith θ_v of band and the relative azimuth

\[ \phi = \arccos\big(\cos(\varphi_v - \varphi_s)\big) \in [0, \pi], \]

sensor minus sun (φ = 0 is the backscatter hotspot, as :mod:radtrax.surface.ross_li reads it), are taken from the tie-point grid (see the module docstring) and interpolated bilinearly onto the image's pixel centres. φ is formed on the tie points, before the interpolation, so azimuths that wrap through 0°/360° do not average across the cut; φ itself is continuous.

Parameters:

Name	Type	Description	Default
image	Any	A georeader S2ImageL1C; its transform, crs and shape set the grid (out_res and window_focus).	required
band	str	The band whose viewing angles to use; Sentinel-2's detectors see each band from a slightly different direction.	'B12'

Returns:

Type	Description
ViewingGeometry	A ViewingGeometry whose leaves are (H, W): the geometry of
ViewingGeometry	every pixel, built under eqx.filter_vmap, so
ViewingGeometry	eqx.filter_vmap over it (twice, or after a reshape to (N,))
ViewingGeometry	hands one pixel's geometry to a forward model.

Raises:

Type	Description
ImportError	If the [geo] extra is not installed.
KeyError	If the metadata has no viewing angles for band.
ValueError	If the sun and view azimuth grids differ.
RuntimeError	If an interpolated cosine is outside (0, 1], NaN included (equinox.EquinoxRuntimeError, from :class:~radtrax.atmosphere.ViewingGeometry).

Source code in src/radtrax/geo/sentinel2.py

def sentinel2_geometry(image: Any, band: str = "B12") -> ViewingGeometry:
    """Per-pixel viewing geometry of an L1C image, on its grid ``(H, W)``.

    The solar zenith ``θ_s``, the viewing zenith ``θ_v`` of ``band`` and the
    relative azimuth

    $$
    \\phi = \\arccos\\big(\\cos(\\varphi_v - \\varphi_s)\\big) \\in [0, \\pi],
    $$

    sensor minus sun (``φ = 0`` is the backscatter hotspot, as
    :mod:`radtrax.surface.ross_li` reads it), are taken from the tie-point
    grid (see the module docstring) and interpolated bilinearly onto the
    image's pixel centres. ``φ`` is formed on the tie points, before the
    interpolation, so azimuths that wrap through 0°/360° do not average
    across the cut; ``φ`` itself is continuous.

    Args:
        image: A georeader ``S2ImageL1C``; its ``transform``, ``crs`` and
            ``shape`` set the grid (``out_res`` and ``window_focus``).
        band: The band whose viewing angles to use; Sentinel-2's detectors
            see each band from a slightly different direction.

    Returns:
        A ``ViewingGeometry`` whose leaves are ``(H, W)``: the geometry of
        every pixel, built under ``eqx.filter_vmap``, so
        ``eqx.filter_vmap`` over it (twice, or after a reshape to ``(N,)``)
        hands one pixel's geometry to a forward model.

    Raises:
        ImportError: If the ``[geo]`` extra is not installed.
        KeyError: If the metadata has no viewing angles for ``band``.
        ValueError: If the sun and view azimuth grids differ.
        RuntimeError: If an interpolated cosine is outside ``(0, 1]``,
            NaN included (``equinox.EquinoxRuntimeError``, from
            :class:`~radtrax.atmosphere.ViewingGeometry`).
    """
    image.read_metadata_tl()
    if band not in image.vza or band not in image.vaa:
        raise MissingKeyError(
            f"sentinel2_geometry: the metadata has no viewing angles for {band!r}; "
            f"it has {sorted(image.vza)}."
        )
    sza = _on_grid(image, image.sza, image)  # (H, W) [deg]
    vza = _on_grid(image, image.vza[band], image)
    raa = _on_grid(image, _relative_azimuth(image, band), image)
    from_angles = eqx.filter_vmap(eqx.filter_vmap(ViewingGeometry.from_angles))
    return from_angles(sza, vza, raa)

sentinel2_radiance

sentinel2_radiance(image: Any, bands: Sequence[str] = ('B11', 'B12')) -> Any

TOA radiance from an L1C image's reflectance, per band and pixel.

\[ \rho_b = \frac{\mathrm{DN}_b}{Q},\qquad L_b = \frac{\rho_b\, E_b \cos\theta_s}{\pi\, d^2} = \frac{\rho_b\, E_b\, U \cos\theta_s}{\pi}, \]

with DN_b the digital numbers georeader reads (the radiometric offset of processing baseline 04.00 already added), Q the quantification value (10 000), E_b the band solar irradiance [W m⁻² nm⁻¹] and U = 1/d² from the product metadata, and θ_s each pixel's solar zenith angle, interpolated as in :func:sentinel2_geometry. ESA computes the L1C reflectance with the per-pixel angle, so this inverts it exactly; georeader's DN_to_radiance uses the tile's mean angle instead.

Parameters:

Name	Type	Description	Default
image	Any	A georeader S2ImageL1C.	required
bands	Sequence[str]	Band names, in order.	('B11', 'B12')

Returns:

Type	Description
Any	A GeoTensor (B, H, W) of radiance [W m⁻² sr⁻¹ nm⁻¹], float64,
Any	on the image grid. Fill pixels (DN = 0) stay 0, the
Any	fill_value_default georeader reads them with.

Raises:

Type	Description
ImportError	If the [geo] extra is not installed.

Source code in src/radtrax/geo/sentinel2.py

def sentinel2_radiance(image: Any, bands: Sequence[str] = ("B11", "B12")) -> Any:
    """TOA radiance from an L1C image's reflectance, per band and pixel.

    $$
    \\rho_b = \\frac{\\mathrm{DN}_b}{Q},\\qquad
    L_b = \\frac{\\rho_b\\, E_b \\cos\\theta_s}{\\pi\\, d^2}
        = \\frac{\\rho_b\\, E_b\\, U \\cos\\theta_s}{\\pi},
    $$

    with ``DN_b`` the digital numbers georeader reads (the radiometric
    offset of processing baseline 04.00 already added), ``Q`` the
    quantification value (10 000), ``E_b`` the band solar irradiance
    [W m⁻² nm⁻¹] and ``U = 1/d²`` from the product metadata, and ``θ_s``
    each pixel's solar zenith angle, interpolated as in
    :func:`sentinel2_geometry`. ESA computes the L1C reflectance with the
    per-pixel angle, so this inverts it exactly; georeader's
    ``DN_to_radiance`` uses the tile's mean angle instead.

    Args:
        image: A georeader ``S2ImageL1C``.
        bands: Band names, in order.

    Returns:
        A ``GeoTensor`` ``(B, H, W)`` of radiance [W m⁻² sr⁻¹ nm⁻¹], float64,
        on the image grid. Fill pixels (``DN = 0``) stay 0, the
        ``fill_value_default`` georeader reads them with.

    Raises:
        ImportError: If the ``[geo]`` extra is not installed.
    """
    GeoTensor = _georeader().GeoTensor
    sub = image.read_from_band_names(list(bands))
    dn = sub.load()  # (B, H, W)
    rho = np.asarray(dn.values, dtype=float) / sub.quantification_value()
    irradiance = sub.solar_irradiance()  # W m⁻² nm⁻¹ per band
    e_b = np.array([irradiance[b] for b in sub.bands])[:, None, None]  # (B, 1, 1)
    image.read_metadata_tl()
    sza = _on_grid(image, image.sza, dn)  # (H, W) [deg]
    radiance = rho * e_b * sub.scale_factor_U() * np.cos(np.deg2rad(sza)) / np.pi
    return GeoTensor(
        radiance,
        transform=dn.transform,
        crs=dn.crs,
        fill_value_default=dn.fill_value_default,
    )

emit

EMIT L1B on radtrax's terms: instrument and per-pixel geometry.

georeader's EMITImage reads an EMIT L1B radiance granule: the per-band centre wavelengths and fwhm [nm] of its current band selection, and, through the companion observation (obs) file, the sun and sensor angles of every pixel, georeferenced through the granule's geographic lookup table (GLT) by EMITImage.observation(name). This module turns them into what a forward model takes (radtrax#114):

• :func:emit_instrument: a :class:~radtrax.instrument.BandInstrument on the model's wavenumber grid, with Gaussian bands of the image's FWHM (exact FWHM, :data:~radtrax.instrument.FWHM_PER_SIGMA, radtrax#28);
• :func:emit_geometry: a :class:~radtrax.atmosphere.ViewingGeometry per pixel of the image grid, from the obs file's four angle bands.

Neither function imports georeader: they read the image's attributes and call its methods, so they need no extra beyond whatever produced the image (EMITImage itself comes with the [geo] extra).

emit_instrument

emit_instrument(image: Any, nu_cm: Float[ArrayLike, ' G'], *, bands: Sequence[int] | None = None, shift_cm: float = 0.0) -> BandInstrument

EMIT bands as a :class:~radtrax.instrument.BandInstrument on a wavenumber grid.

Band b is a Gaussian in wavelength with the image's centre \(c_b\) and FWHM \(w_b\),

\[ f_b(\lambda) = \exp\!\left(-\frac{(\lambda - c_b)^2}{2\sigma_b^2}\right), \qquad \sigma_b = \frac{w_b}{2\sqrt{2\ln 2}}, \]

sampled at \(\lambda_g = 10^7 / \nu_g\) by :class:~radtrax.instrument.SpectralResponseFunction (trapezoidal weights, rows normalised to 1) and turned onto the wavenumber grid by :meth:BandInstrument.from_srf <radtrax.instrument.BandInstrument.from_srf>.

A band of which the grid holds less than 99.9 % of the Gaussian's area, \(\Phi((\lambda_{\max} - c_b)/\sigma_b) - \Phi((\lambda_{\min} - c_b)/\sigma_b)\), is dropped, with one warning naming every dropped band: a band is kept when its centre lies at least \(3.09\,\sigma_b \approx 1.31\,w_b\) inside the grid's wavelength range (a one-sided tail of 0.1 %).

Parameters:

Name	Type	Description	Default
image	Any	A georeader EMITImage (its wavelengths and fwhm, [nm], for the current band selection), or anything with those two attributes.	required
nu_cm	Float[ArrayLike, ' G']	The model's wavenumber grid [cm⁻¹], strictly increasing, (G,).	required
bands	Sequence[int] | None	Indices into image.wavelengths to keep, in order; all bands when None.	None
shift_cm	float	Spectral shift δ [cm⁻¹], a differentiable leaf.	0.0

Returns:

Type	Description
BandInstrument	The instrument, (G,) → (B,). Band i of image.wavelengths
BandInstrument	is named "B<i>" (0-based, as georeader's band selection counts),
BandInstrument	so the names say which bands were kept.

Warns:

Type	Description
UserWarning	If bands are dropped because the grid does not cover them.

Raises:

Type	Description
GridError	If nu_cm is not 1-D and strictly increasing, or the grid covers none of the requested bands.
ShapeError	If wavelengths and fwhm differ in shape.
ParameterError	If a band index is out of range, or a FWHM is not positive.

Examples:

>>> import numpy as np
>>> from types import SimpleNamespace
>>> image = SimpleNamespace(
...     wavelengths=np.array([1000.0, 2300.0, 2310.0]),
...     fwhm=np.array([8.0, 8.5, 8.5]),
... )
>>> nu = np.arange(4250.0, 4400.0, 0.1)  # 2273–2353 nm
>>> instrument = emit_instrument(image, nu, bands=[1, 2])
>>> instrument.band_names
('B1', 'B2')
>>> [round(float(y), 12) for y in instrument(np.ones(nu.size))]
[1.0, 1.0]

Source code in src/radtrax/geo/emit.py

def emit_instrument(
    image: Any,
    nu_cm: Float[ArrayLike, " G"],
    *,
    bands: Sequence[int] | None = None,
    shift_cm: float = 0.0,
) -> BandInstrument:
    """EMIT bands as a :class:`~radtrax.instrument.BandInstrument` on a wavenumber grid.

    Band ``b`` is a Gaussian in wavelength with the image's centre $c_b$
    and FWHM $w_b$,

    $$
    f_b(\\lambda) = \\exp\\!\\left(-\\frac{(\\lambda - c_b)^2}{2\\sigma_b^2}\\right),
    \\qquad \\sigma_b = \\frac{w_b}{2\\sqrt{2\\ln 2}},
    $$

    sampled at $\\lambda_g = 10^7 / \\nu_g$ by
    :class:`~radtrax.instrument.SpectralResponseFunction` (trapezoidal
    weights, rows normalised to 1) and turned onto the wavenumber grid by
    :meth:`BandInstrument.from_srf <radtrax.instrument.BandInstrument.from_srf>`.

    A band of which the grid holds less than 99.9 % of the Gaussian's area,
    $\\Phi((\\lambda_{\\max} - c_b)/\\sigma_b) - \\Phi((\\lambda_{\\min} - c_b)/\\sigma_b)$,
    is dropped, with one warning naming every dropped band: a band is kept
    when its centre lies at least $3.09\\,\\sigma_b \\approx 1.31\\,w_b$ inside
    the grid's wavelength range (a one-sided tail of 0.1 %).

    Args:
        image: A georeader ``EMITImage`` (its ``wavelengths`` and ``fwhm``,
            [nm], for the current band selection), or anything with those
            two attributes.
        nu_cm: The model's wavenumber grid [cm⁻¹], strictly increasing,
            ``(G,)``.
        bands: Indices into ``image.wavelengths`` to keep, in order; all
            bands when ``None``.
        shift_cm: Spectral shift ``δ`` [cm⁻¹], a differentiable leaf.

    Returns:
        The instrument, ``(G,) → (B,)``. Band ``i`` of ``image.wavelengths``
        is named ``"B<i>"`` (0-based, as georeader's band selection counts),
        so the names say which bands were kept.

    Warns:
        UserWarning: If bands are dropped because the grid does not cover
            them.

    Raises:
        GridError: If ``nu_cm`` is not 1-D and strictly increasing, or the
            grid covers none of the requested bands.
        ShapeError: If ``wavelengths`` and ``fwhm`` differ in shape.
        ParameterError: If a band index is out of range, or a FWHM is not
            positive.

    Examples:
        >>> import numpy as np
        >>> from types import SimpleNamespace
        >>> image = SimpleNamespace(
        ...     wavelengths=np.array([1000.0, 2300.0, 2310.0]),
        ...     fwhm=np.array([8.0, 8.5, 8.5]),
        ... )
        >>> nu = np.arange(4250.0, 4400.0, 0.1)  # 2273–2353 nm
        >>> instrument = emit_instrument(image, nu, bands=[1, 2])
        >>> instrument.band_names
        ('B1', 'B2')
        >>> [round(float(y), 12) for y in instrument(np.ones(nu.size))]
        [1.0, 1.0]
    """
    nu = np.asarray(nu_cm, dtype=float)
    if nu.ndim != 1 or nu.size < 2 or np.any(np.diff(nu) <= 0.0):
        raise GridError(
            "emit_instrument: `nu_cm` must be 1-D, strictly increasing, with at "
            "least two points."
        )
    centres = np.asarray(image.wavelengths, dtype=float).reshape(-1)
    fwhm = np.asarray(image.fwhm, dtype=float).reshape(-1)
    if centres.shape != fwhm.shape:
        raise ShapeError(
            f"emit_instrument: {centres.size} wavelengths but {fwhm.size} FWHMs."
        )
    index = np.arange(centres.size) if bands is None else np.asarray(bands, int)
    if index.ndim != 1 or np.any((index < 0) | (index >= centres.size)):
        raise ParameterError(
            f"emit_instrument: band indices must be in [0, {centres.size}), "
            f"got {list(np.atleast_1d(index))}."
        )
    if np.any(~(fwhm[index] > 0.0)):
        raise ParameterError("emit_instrument: every FWHM must be positive.")

    wl = 1e7 / nu[::-1]  # ascending [nm]
    sigma = fwhm[index] / FWHM_PER_SIGMA
    coverage = np.array(
        [
            _normal_cdf((wl[-1] - c) / s) - _normal_cdf((wl[0] - c) / s)
            for c, s in zip(centres[index], sigma, strict=True)
        ]
    )
    keep = coverage >= _MIN_COVERAGE
    if not np.any(keep):
        raise GridError(
            f"emit_instrument: the grid [{wl[0]:g}, {wl[-1]:g}] nm covers none of "
            "the requested bands."
        )
    if not np.all(keep):
        dropped = ", ".join(f"B{i} ({centres[i]:.2f} nm)" for i in index[~keep])
        warnings.warn(
            f"emit_instrument: the grid [{wl[0]:g}, {wl[-1]:g}] nm holds less than "
            f"{100 * _MIN_COVERAGE:g} % of the response of {dropped}; dropped.",
            UserWarning,
            stacklevel=2,
        )
    kept = index[keep]
    srf = SpectralResponseFunction(
        wavelengths_hr_nm=wl,
        band_centers_nm=centres[kept],
        band_widths_nm=fwhm[kept],
        band_names=tuple(f"B{i}" for i in kept),
        srf_type="gaussian",
    )
    return BandInstrument.from_srf(srf, nu, shift_cm=shift_cm)

emit_geometry

emit_geometry(image: Any) -> ViewingGeometry

Per-pixel viewing geometry of an EMIT image, on its grid (H, W).

The solar zenith \(\theta_s\), the sensor zenith \(\theta_v\) and the relative azimuth

\[ \phi = \arccos\big(\cos(\varphi_v - \varphi_s)\big) \in [0, \pi], \]

sensor minus sun, both azimuths clockwise from north toward the sun and toward the sensor (φ = 0 is the backscatter hotspot, as :func:~radtrax.geo.sentinel2_geometry forms it), from the obs file's bands as image.observation(name) georeferences them (nearest sensor pixel through the GLT, no interpolation).

Pixels outside the swath (the GLT's no-data) and obs fill values get the mean of the valid pixels' angles, so that the whole grid is a valid ViewingGeometry; mask them with the radiance's fill value or image.valid_glt, as for any other per-pixel product.

Parameters:

Name	Type	Description	Default
image	Any	A georeader EMITImage whose obs file is available (image.nc_ds_obs downloads it next to the radiance file when it is missing, which needs Earthdata credentials).	required

Returns:

Type	Description
ViewingGeometry	A ViewingGeometry whose leaves are (H, W), built under
ViewingGeometry	eqx.filter_vmap as :func:~radtrax.geo.sentinel2_geometry's.

Raises:

Type	Description
MissingKeyError	If the obs file lacks one of the four angle bands.
ShapeError	If the four angle grids differ in shape.
DegenerateInputError	If no pixel has valid angles.
RuntimeError	If a valid zenith angle is 90° or more (equinox.EquinoxRuntimeError, from :class:~radtrax.atmosphere.ViewingGeometry).

Source code in src/radtrax/geo/emit.py

def emit_geometry(image: Any) -> ViewingGeometry:
    """Per-pixel viewing geometry of an EMIT image, on its grid ``(H, W)``.

    The solar zenith $\\theta_s$, the sensor zenith $\\theta_v$ and the
    relative azimuth

    $$
    \\phi = \\arccos\\big(\\cos(\\varphi_v - \\varphi_s)\\big) \\in [0, \\pi],
    $$

    sensor minus sun, both azimuths clockwise from north toward the sun and
    toward the sensor (``φ = 0`` is the backscatter hotspot, as
    :func:`~radtrax.geo.sentinel2_geometry` forms it), from the obs file's
    bands as ``image.observation(name)`` georeferences them (nearest
    sensor pixel through the GLT, no interpolation).

    Pixels outside the swath (the GLT's no-data) and obs fill values get the
    mean of the valid pixels' angles, so that the whole grid is a valid
    ``ViewingGeometry``; mask them with the radiance's fill value or
    ``image.valid_glt``, as for any other per-pixel product.

    Args:
        image: A georeader ``EMITImage`` whose obs file is available
            (``image.nc_ds_obs`` downloads it next to the radiance file
            when it is missing, which needs Earthdata credentials).

    Returns:
        A ``ViewingGeometry`` whose leaves are ``(H, W)``, built under
        ``eqx.filter_vmap`` as :func:`~radtrax.geo.sentinel2_geometry`'s.

    Raises:
        MissingKeyError: If the obs file lacks one of the four angle bands.
        ShapeError: If the four angle grids differ in shape.
        DegenerateInputError: If no pixel has valid angles.
        RuntimeError: If a valid zenith angle is 90° or more
            (``equinox.EquinoxRuntimeError``, from
            :class:`~radtrax.atmosphere.ViewingGeometry`).
    """
    names = [str(n) for n in image.observation_bands]
    angles = [
        _observation(image, names, key)
        for key in (_SUN_ZENITH, _SUN_AZIMUTH, _SENSOR_ZENITH, _SENSOR_AZIMUTH)
    ]
    if len({a.shape for a in angles}) != 1:
        raise ShapeError(
            f"emit_geometry: the angle grids differ in shape {[a.shape for a in angles]}."
        )
    valid = np.logical_and.reduce([np.isfinite(a) for a in angles])
    if not np.any(valid):
        raise DegenerateInputError("emit_geometry: no pixel has valid angles.")
    s_zen, s_az, v_zen, v_az = (np.where(valid, a, np.mean(a[valid])) for a in angles)
    # The azimuths' cosine difference is formed per pixel, so the fill mean of
    # an azimuth that wraps through 0°/360° only affects masked pixels.
    raa = np.rad2deg(np.arccos(np.cos(np.deg2rad(v_az - s_az))))
    from_angles = eqx.filter_vmap(eqx.filter_vmap(ViewingGeometry.from_angles))
    return from_angles(s_zen, v_zen, raa)

landsat

Landsat 8/9 OLI and TIRS: SRF, instrument, noise, radiance and geometry.

OLI (Landsat 8) and OLI-2 (Landsat 9) carry a SWIR band pair, B6 at 1.6 µm and B7 at 2.2 µm, like Sentinel-2's B11 and B12, so the multi-pass family and the band-table inversion apply unchanged (radtrax#121). TIRS (Landsat 8) and TIRS-2 (Landsat 9) add the thermal pair, B10 at about 10.9 µm and B11 at about 12 µm (radtrax#225). georeader has neither a Collection 2 Level-1 reader (spaceml-org/georeader#99) nor the Landsat spectral responses (spaceml-org/georeader#100), so this module works on plain inputs:

• :func:landsat_srf and :func:landsat_instrument: the USGS band-average relative spectral responses, bundled with radtrax (geo/data/landsat_oli_rsr.csv, B6, B7, B10 and B11 of each satellite; its header records the source workbooks, their SHA-256 and the retrieval date), as a :class:~radtrax.instrument.SpectralResponseFunction and a :class:~radtrax.instrument.BandInstrument;
• :func:landsat_noise: the thermal bands' NEΔT (LANDSAT_NETD_K) as an :class:~radtrax.instrument.NEdTNoise; the TIRS-2 calibration residual is LANDSAT_CALIBRATION_K;
• :func:landsat_radiance: TOA radiance from the Level-1 digital numbers and the MTL's rescaling, L = M Q + A;
• :func:landsat_geometry: a per-pixel :class:~radtrax.atmosphere.ViewingGeometry from the four Collection 2 angle bands (SZA, SAA, VZA, VAA).

Units. Level-1 radiance is per µm (W m⁻² sr⁻¹ µm⁻¹, the MTL's unit); the forward model's is per cm⁻¹. landsat_instrument(..., per_um=True) takes radiance per cm⁻¹ and returns band averages per µm, the conversion inside the band integral (:meth:BandInstrument.from_srf <radtrax.instrument.BandInstrument.from_srf>), so its output, the noise :func:landsat_noise builds on it and :func:landsat_radiance share one unit.

The table is long-format (satellite, band, wavelength, response). None of these functions needs the [geo] extra; :func:landsat_radiance returns a GeoTensor when it is given one.

landsat_srf

landsat_srf(satellite: str = 'L8', *, bands: Sequence[str] = ('B6', 'B7'), wavelengths_hr_nm: Sequence[float] | Float[ArrayLike, ' Ghr'] | None = None) -> SpectralResponseFunction

Build a Landsat 8/9 SRF from the bundled USGS band-average RSRs.

Each band's response is placed on one 1 nm grid spanning the requested bands, zero where it was not measured, and projected onto wavelengths_hr_nm with the linear-interpolation (hat) basis, as :func:~radtrax.geo.sentinel2_srf does: band integration is then exact for any radiance linear between grid points. OLI's responses are published every 1 nm and used as they are; TIRS's every 50 nm, and are taken as piecewise linear between the published values (interpolated onto the 1 nm grid). The values are used as published, including the few small negative values in Landsat 8 OLI's tails and the floor of about 1e-3 Landsat 8 TIRS reports across 9–14 µm (so a grid narrower than 9–14 µm misses more than 0.1 % of an L8 thermal band, and warns). The responses are relative to radiance (energy), as the Landsat 8 TIRS workbook states, so the band average is ∫R L dλ / ∫R dλ.

Parameters:

Name	Type	Description	Default
satellite	str	"L8" (OLI, TIRS) or "L9" (OLI-2, TIRS-2).	'L8'
bands	Sequence[str]	Band names, in order; "B6" (SWIR1), "B7" (SWIR2), "B10" (~10.9 µm) and "B11" (~12 µm) are bundled.	('B6', 'B7')
wavelengths_hr_nm	Sequence[float] | Float[ArrayLike, ' Ghr'] | None	The high-resolution grid [nm] to integrate on, strictly increasing. Defaults to the table's 1 nm grid.	None

Returns:

Type	Description
SpectralResponseFunction	A custom SRF whose band centres are the response-weighted mean
SpectralResponseFunction	wavelengths and widths the equivalent widths ∫R dλ / max R, both
SpectralResponseFunction	computed on the 1 nm table.

Warns:

Type	Description
UserWarning	If the grid misses more than 0.1 % of a band's response; the band is then normalised over the part it covers.

Raises:

Type	Description
ValueError	If satellite is unknown or the grid is not 1-D and strictly increasing.
KeyError	If a band is not bundled.

Examples:

>>> srf = landsat_srf("L9")
>>> srf.band_names
('B6', 'B7')
>>> [round(float(c)) for c in srf.band_centers_nm]
[1608, 2201]
>>> tirs = landsat_srf("L9", bands=("B10", "B11"))
>>> [round(float(c)) for c in tirs.band_centers_nm]
[10837, 12025]

Source code in src/radtrax/geo/landsat.py

def landsat_srf(
    satellite: str = "L8",
    *,
    bands: Sequence[str] = ("B6", "B7"),
    wavelengths_hr_nm: Sequence[float] | Float[ArrayLike, " Ghr"] | None = None,
) -> SpectralResponseFunction:
    """Build a Landsat 8/9 SRF from the bundled USGS band-average RSRs.

    Each band's response is placed on one 1 nm grid spanning the requested
    bands, zero where it was not measured, and projected onto
    ``wavelengths_hr_nm`` with the linear-interpolation (hat) basis, as
    :func:`~radtrax.geo.sentinel2_srf` does: band integration is then exact
    for any radiance linear between grid points. OLI's responses are
    published every 1 nm and used as they are; TIRS's every 50 nm, and are
    taken as piecewise linear between the published values (interpolated
    onto the 1 nm grid). The values are used as published, including the
    few small negative values in Landsat 8 OLI's tails and the floor of
    about 1e-3 Landsat 8 TIRS reports across 9–14 µm (so a grid narrower
    than 9–14 µm misses more than 0.1 % of an L8 thermal band, and warns).
    The responses are relative to radiance (energy), as the Landsat 8 TIRS
    workbook states, so the band average is ``∫R L dλ / ∫R dλ``.

    Args:
        satellite: ``"L8"`` (OLI, TIRS) or ``"L9"`` (OLI-2, TIRS-2).
        bands: Band names, in order; ``"B6"`` (SWIR1), ``"B7"`` (SWIR2),
            ``"B10"`` (~10.9 µm) and ``"B11"`` (~12 µm) are bundled.
        wavelengths_hr_nm: The high-resolution grid [nm] to integrate on,
            strictly increasing. Defaults to the table's 1 nm grid.

    Returns:
        A ``custom`` SRF whose band centres are the response-weighted mean
        wavelengths and widths the equivalent widths ``∫R dλ / max R``, both
        computed on the 1 nm table.

    Warns:
        UserWarning: If the grid misses more than ``0.1 %`` of a band's
            response; the band is then normalised over the part it covers.

    Raises:
        ValueError: If ``satellite`` is unknown or the grid is not 1-D and
            strictly increasing.
        KeyError: If a band is not bundled.

    Examples:
        >>> srf = landsat_srf("L9")
        >>> srf.band_names
        ('B6', 'B7')
        >>> [round(float(c)) for c in srf.band_centers_nm]
        [1608, 2201]
        >>> tirs = landsat_srf("L9", bands=("B10", "B11"))
        >>> [round(float(c)) for c in tirs.band_centers_nm]
        [10837, 12025]
    """
    table = _rsr_table()
    satellites = sorted({sat for sat, _ in table})
    if satellite not in satellites:
        raise ParameterError(
            f"landsat_srf: unknown satellite {satellite!r}; use one of {satellites}."
        )
    names = tuple(bands)
    missing = [b for b in names if (satellite, b) not in table]
    if missing or not names:
        bundled = sorted(b for sat, b in table if sat == satellite)
        raise MissingKeyError(
            f"landsat_srf: no band(s) {missing} for {satellite}; bundled: {bundled}."
        )
    samples = [table[satellite, b] for b in names]  # every 1 nm (OLI), 50 nm (TIRS)
    return _tabulated_srf("landsat_srf", samples, names, wavelengths_hr_nm)

landsat_instrument

landsat_instrument(satellite: str, nu_cm: Float[ArrayLike, ' G'], *, bands: Sequence[str] = ('B6', 'B7'), shift_cm: float = 0.0, per_um: bool = False) -> BandInstrument

Landsat 8/9 bands as a :class:~radtrax.instrument.BandInstrument.

The RSRs (:func:landsat_srf) projected onto the wavelengths of the model grid, λ_g = 10⁷ / ν_g, and turned onto the wavenumber grid by :meth:BandInstrument.from_srf <radtrax.instrument.BandInstrument.from_srf>:

\[ y_b = \sum_g f_{bg}\, L(\nu_g + \delta), \]

with f the SRF's row-normalised matrix, columns reversed. At zero shift it is :meth:SpectralResponseFunction.apply <radtrax.instrument.SpectralResponseFunction.apply> of the same SRF on the spectrum in ascending wavelength.

For the thermal bands pass per_um=True: the instrument then takes radiance per cm⁻¹ [W m⁻² sr⁻¹ (cm⁻¹)⁻¹] and returns the band average per µm [W m⁻² sr⁻¹ µm⁻¹], the unit of :func:landsat_radiance, with f_bg = w_bg ν_g² / 10⁴ (rows no longer summing to 1).

Parameters:

Name	Type	Description	Default
satellite	str	"L8" or "L9".	required
nu_cm	Float[ArrayLike, ' G']	The model's wavenumber grid [cm⁻¹], strictly increasing, (G,); it should cover each band's response.	required
bands	Sequence[str]	Band names, in order.	('B6', 'B7')
shift_cm	float	Spectral shift δ [cm⁻¹], a differentiable leaf.	0.0
per_um	bool	Band radiance per µm from radiance per cm⁻¹.	False

Returns:

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

Warns:

Type	Description
UserWarning	If the grid misses more than 0.1 % of a band's response (from :func:landsat_srf).

Raises:

Type	Description
KeyError	If a band is not bundled.
ValueError	If nu_cm is not 1-D and strictly increasing, or the satellite is unknown.

Examples:

>>> import numpy as np
>>> nu = np.arange(4200.0, 6650.0)  # 1504–2381 nm, B6 and B7
>>> instrument = landsat_instrument("L8", nu)
>>> instrument.band_names
('B6', 'B7')
>>> [round(float(y), 6) for y in instrument(np.ones(nu.size))]
[1.0, 1.0]
>>> nu_tir = np.arange(700.0, 1120.0, 0.5)  # 8.9–14.3 µm, B10 and B11
>>> tirs = landsat_instrument(
...     "L9", nu_tir, bands=("B10", "B11"), per_um=True
... )
>>> [round(float(y), 9) for y in tirs(2.0 * 1e4 / nu_tir**2)]  # 2 per µm
[2.0, 2.0]

Source code in src/radtrax/geo/landsat.py

def landsat_instrument(
    satellite: str,
    nu_cm: Float[ArrayLike, " G"],
    *,
    bands: Sequence[str] = ("B6", "B7"),
    shift_cm: float = 0.0,
    per_um: bool = False,
) -> BandInstrument:
    """Landsat 8/9 bands as a :class:`~radtrax.instrument.BandInstrument`.

    The RSRs (:func:`landsat_srf`) projected onto the wavelengths of the
    model grid, ``λ_g = 10⁷ / ν_g``, and turned onto the wavenumber grid by
    :meth:`BandInstrument.from_srf
    <radtrax.instrument.BandInstrument.from_srf>`:

    $$
    y_b = \\sum_g f_{bg}\\, L(\\nu_g + \\delta),
    $$

    with ``f`` the SRF's row-normalised matrix, columns reversed. At zero
    shift it is :meth:`SpectralResponseFunction.apply
    <radtrax.instrument.SpectralResponseFunction.apply>` of the same SRF on
    the spectrum in ascending wavelength.

    For the thermal bands pass ``per_um=True``: the instrument then takes
    radiance per cm⁻¹ [W m⁻² sr⁻¹ (cm⁻¹)⁻¹] and returns the band average per
    µm [W m⁻² sr⁻¹ µm⁻¹], the unit of :func:`landsat_radiance`, with
    ``f_bg = w_bg ν_g² / 10⁴`` (rows no longer summing to 1).

    Args:
        satellite: ``"L8"`` or ``"L9"``.
        nu_cm: The model's wavenumber grid [cm⁻¹], strictly increasing,
            ``(G,)``; it should cover each band's response.
        bands: Band names, in order.
        shift_cm: Spectral shift ``δ`` [cm⁻¹], a differentiable leaf.
        per_um: Band radiance per µm from radiance per cm⁻¹.

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

    Warns:
        UserWarning: If the grid misses more than ``0.1 %`` of a band's
            response (from :func:`landsat_srf`).

    Raises:
        KeyError: If a band is not bundled.
        ValueError: If ``nu_cm`` is not 1-D and strictly increasing, or the
            satellite is unknown.

    Examples:
        >>> import numpy as np
        >>> nu = np.arange(4200.0, 6650.0)  # 1504–2381 nm, B6 and B7
        >>> instrument = landsat_instrument("L8", nu)
        >>> instrument.band_names
        ('B6', 'B7')
        >>> [round(float(y), 6) for y in instrument(np.ones(nu.size))]
        [1.0, 1.0]
        >>> nu_tir = np.arange(700.0, 1120.0, 0.5)  # 8.9–14.3 µm, B10 and B11
        >>> tirs = landsat_instrument(
        ...     "L9", nu_tir, bands=("B10", "B11"), per_um=True
        ... )
        >>> [round(float(y), 9) for y in tirs(2.0 * 1e4 / nu_tir**2)]  # 2 per µm
        [2.0, 2.0]
    """
    nu = np.asarray(nu_cm, dtype=float)
    if nu.ndim != 1 or nu.size < 2 or np.any(np.diff(nu) <= 0.0):
        raise GridError(
            "landsat_instrument: `nu_cm` must be 1-D, strictly increasing, "
            "with at least two points."
        )
    srf = landsat_srf(satellite, bands=bands, wavelengths_hr_nm=1e7 / nu[::-1])
    return BandInstrument.from_srf(srf, nu, shift_cm=shift_cm, per_um=per_um)

landsat_noise

landsat_noise(satellite: str, instrument: BandInstrument, *, t_ref_k: Real[ArrayLike, ''] = 300.0) -> NEdTNoise

The thermal bands' NEΔT as a noise model on instrument's bands.

\[ σ_b = \mathrm{NEΔT}_b \, \frac{\partial \bar B_b}{\partial T}\Big|_{T_\mathrm{ref}}, \qquad S_ε = \mathrm{diag}(σ^2), \]

with \(\bar B_b\) the Planck radiance integrated by instrument and NEΔT the median on-orbit value at 300 K (LANDSAT_NETD_K). \(σ\) is in the instrument's units: build it with per_um=True to compare \(σ\) with :func:landsat_radiance. The same \(σ\) is a different NEΔT at another scene temperature (:meth:NEdTNoise.netd_at <radtrax.instrument.NEdTNoise.netd_at>), larger in a colder scene. The calibration residual (LANDSAT_CALIBRATION_K) is systematic and is not part of this noise.

Parameters:

Name	Type	Description	Default
satellite	str	"L8" or "L9".	required
instrument	BandInstrument	The thermal bands, from :func:landsat_instrument; its band_names pick the NEΔT values.	required
t_ref_k	Real[ArrayLike, '']	The temperature NEΔT is quoted at [K], a scalar.	300.0

Returns:

Name	Type	Description
An	NEdTNoise	class:~radtrax.instrument.NEdTNoise in the "reference" mode.

Raises:

Type	Description
ValueError	If satellite is unknown.
KeyError	If a band has no tabulated NEΔT (OLI's B6 and B7).

Examples:

>>> import numpy as np
>>> nu = np.arange(700.0, 1120.0, 0.5)
>>> tirs = landsat_instrument("L9", nu, bands=("B10", "B11"), per_um=True)
>>> sigma = landsat_noise("L9", tirs).radiance_sigma(np.zeros(2))
>>> [f"{float(s):.2e}" for s in sigma]  # W m⁻² sr⁻¹ µm⁻¹
['6.92e-03', '8.00e-03']

Source code in src/radtrax/geo/landsat.py

def landsat_noise(
    satellite: str, instrument: BandInstrument, *, t_ref_k: Real[ArrayLike, ""] = 300.0
) -> NEdTNoise:
    """The thermal bands' NEΔT as a noise model on ``instrument``'s bands.

    $$
    σ_b = \\mathrm{NEΔT}_b \\, \\frac{\\partial \\bar B_b}{\\partial T}\\Big|_{T_\\mathrm{ref}},
    \\qquad S_ε = \\mathrm{diag}(σ^2),
    $$

    with $\\bar B_b$ the Planck radiance integrated by ``instrument`` and
    NEΔT the median on-orbit value at 300 K (``LANDSAT_NETD_K``). $σ$ is in
    the instrument's units: build it with ``per_um=True`` to compare $σ$
    with :func:`landsat_radiance`. The same $σ$ is a different NEΔT at
    another scene temperature (:meth:`NEdTNoise.netd_at
    <radtrax.instrument.NEdTNoise.netd_at>`), larger in a colder scene. The
    calibration residual (``LANDSAT_CALIBRATION_K``) is systematic and is
    not part of this noise.

    Args:
        satellite: ``"L8"`` or ``"L9"``.
        instrument: The thermal bands, from :func:`landsat_instrument`; its
            ``band_names`` pick the NEΔT values.
        t_ref_k: The temperature NEΔT is quoted at [K], a scalar.

    Returns:
        An :class:`~radtrax.instrument.NEdTNoise` in the ``"reference"`` mode.

    Raises:
        ValueError: If ``satellite`` is unknown.
        KeyError: If a band has no tabulated NEΔT (OLI's B6 and B7).

    Examples:
        >>> import numpy as np
        >>> nu = np.arange(700.0, 1120.0, 0.5)
        >>> tirs = landsat_instrument("L9", nu, bands=("B10", "B11"), per_um=True)
        >>> sigma = landsat_noise("L9", tirs).radiance_sigma(np.zeros(2))
        >>> [f"{float(s):.2e}" for s in sigma]  # W m⁻² sr⁻¹ µm⁻¹
        ['6.92e-03', '8.00e-03']
    """
    if satellite not in LANDSAT_NETD_K:
        raise ParameterError(
            f"landsat_noise: unknown satellite {satellite!r}; use one of "
            f"{sorted(LANDSAT_NETD_K)}."
        )
    table = LANDSAT_NETD_K[satellite]
    missing = [b for b in instrument.band_names if b not in table]
    if missing:
        raise MissingKeyError(
            f"landsat_noise: no NEΔT for band(s) {missing} of {satellite}; "
            f"tabulated: {sorted(table)}."
        )
    netd = np.array([table[b] for b in instrument.band_names])
    return NEdTNoise(netd, instrument, t_ref_k=t_ref_k)

landsat_radiance

landsat_radiance(dn: Real[ArrayLike, ' B *shape'], *, bands: Sequence[str] = ('B6', 'B7'), mtl: Mapping[str, Any] | None = None, mult: Sequence[float] | Float[ArrayLike, ' B'] | None = None, add: Sequence[float] | Float[ArrayLike, ' B'] | None = None) -> Float[ndarray, ' B *shape']

TOA radiance from Level-1 digital numbers, per band and pixel.

\[ L_b = M_b\, Q_b + A_b, \]

with Q_b the calibrated digital numbers and M_b, A_b the MTL's RADIANCE_MULT_BAND_n and RADIANCE_ADD_BAND_n (group LEVEL1_RADIOMETRIC_RESCALING). Q = 0 is fill (no data) in Collection 2 and stays 0. The thermal bands B10 and B11 rescale the same way; their radiance is a band average per µm, so model it with landsat_instrument(..., per_um=True). The product resamples TIRS's 100 m samples onto the 30 m grid; no point-spread or ground-sampling operator is applied here.

Parameters:

Name	Type	Description	Default
dn	Real[ArrayLike, ' B *shape']	The digital numbers, (B, ...), one band per entry of bands; a georeader GeoTensor gives a GeoTensor.	required
bands	Sequence[str]	Band names, in order, "B<n>".	('B6', 'B7')
mtl	Mapping[str, Any] | None	A parsed MTL (the _MTL.json loaded, or any mapping), nested or flat; the rescaling keys are looked up at any depth and may be strings.	None
mult	Sequence[float] | Float[ArrayLike, ' B'] | None	M_b [W m⁻² sr⁻¹ µm⁻¹ per DN], one per band, instead of mtl.	None
add	Sequence[float] | Float[ArrayLike, ' B'] | None	A_b [W m⁻² sr⁻¹ µm⁻¹], one per band, with mult.	None

Returns:

Type	Description
Float[ndarray, ' B *shape']	Radiance [W m⁻² sr⁻¹ µm⁻¹, the unit of the MTL rescaling], float64,
Float[ndarray, ' B *shape']	(B, ...); a GeoTensor on dn's grid if dn is one.

Raises:

Type	Description
ValueError	If neither or both of mtl and mult/add are given, or the shapes do not match bands.
KeyError	If the MTL has no rescaling for a band.

Examples:

>>> import numpy as np
>>> dn = np.array([[[0, 10_000]], [[20_000, 30_000]]], dtype=np.uint16)
>>> landsat_radiance(dn, mult=[1e-3, 2e-4], add=[-5.0, -1.0]).tolist()
[[[0.0, 5.0]], [[3.0, 5.0]]]

Source code in src/radtrax/geo/landsat.py

def landsat_radiance(
    dn: Real[ArrayLike, " B *shape"],
    *,
    bands: Sequence[str] = ("B6", "B7"),
    mtl: Mapping[str, Any] | None = None,
    mult: Sequence[float] | Float[ArrayLike, " B"] | None = None,
    add: Sequence[float] | Float[ArrayLike, " B"] | None = None,
) -> Float[np.ndarray, " B *shape"]:
    """TOA radiance from Level-1 digital numbers, per band and pixel.

    $$
    L_b = M_b\\, Q_b + A_b,
    $$

    with ``Q_b`` the calibrated digital numbers and ``M_b``, ``A_b`` the
    MTL's ``RADIANCE_MULT_BAND_n`` and ``RADIANCE_ADD_BAND_n`` (group
    ``LEVEL1_RADIOMETRIC_RESCALING``). ``Q = 0`` is fill (no data) in
    Collection 2 and stays 0. The thermal bands B10 and B11 rescale the same
    way; their radiance is a band average per µm, so model it with
    ``landsat_instrument(..., per_um=True)``. The product resamples TIRS's
    100 m samples onto the 30 m grid; no point-spread or ground-sampling
    operator is applied here.

    Args:
        dn: The digital numbers, ``(B, ...)``, one band per entry of
            ``bands``; a georeader ``GeoTensor`` gives a ``GeoTensor``.
        bands: Band names, in order, ``"B<n>"``.
        mtl: A parsed MTL (the ``_MTL.json`` loaded, or any mapping), nested
            or flat; the rescaling keys are looked up at any depth and may
            be strings.
        mult: ``M_b`` [W m⁻² sr⁻¹ µm⁻¹ per DN], one per band, instead of
            ``mtl``.
        add: ``A_b`` [W m⁻² sr⁻¹ µm⁻¹], one per band, with ``mult``.

    Returns:
        Radiance [W m⁻² sr⁻¹ µm⁻¹, the unit of the MTL rescaling], float64,
        ``(B, ...)``; a ``GeoTensor`` on ``dn``'s grid if ``dn`` is one.

    Raises:
        ValueError: If neither or both of ``mtl`` and ``mult``/``add`` are
            given, or the shapes do not match ``bands``.
        KeyError: If the MTL has no rescaling for a band.

    Examples:
        >>> import numpy as np
        >>> dn = np.array([[[0, 10_000]], [[20_000, 30_000]]], dtype=np.uint16)
        >>> landsat_radiance(dn, mult=[1e-3, 2e-4], add=[-5.0, -1.0]).tolist()
        [[[0.0, 5.0]], [[3.0, 5.0]]]
    """
    names = tuple(bands)
    if (mtl is None) == (mult is None and add is None):
        raise ParameterError(
            "landsat_radiance: pass either `mtl` or both `mult` and `add`."
        )
    if mtl is not None:
        m = np.array([_mtl_value(mtl, f"RADIANCE_MULT_BAND_{_n(b)}") for b in names])
        a = np.array([_mtl_value(mtl, f"RADIANCE_ADD_BAND_{_n(b)}") for b in names])
    elif mult is None or add is None:
        raise ParameterError("landsat_radiance: pass both `mult` and `add`.")
    else:
        m, a = np.asarray(mult, dtype=float), np.asarray(add, dtype=float)
    q = np.asarray(dn)
    if m.shape != (len(names),) or a.shape != m.shape or q.shape[:1] != m.shape:
        raise ShapeError(
            f"landsat_radiance: {len(names)} band(s) {names}, but `dn` is "
            f"{q.shape} and the rescaling {m.shape}, {a.shape}."
        )
    shape = (-1,) + (1,) * (q.ndim - 1)
    q = q.astype(float)
    radiance = np.where(q == 0.0, 0.0, m.reshape(shape) * q + a.reshape(shape))
    if hasattr(dn, "transform") and hasattr(dn, "crs"):  # a GeoTensor
        from georeader.geotensor import GeoTensor

        return GeoTensor(
            radiance, transform=dn.transform, crs=dn.crs, fill_value_default=0.0
        )
    return radiance

landsat_geometry

landsat_geometry(sza: Real[ArrayLike, ' H W'], saa: Real[ArrayLike, ' H W'], vza: Real[ArrayLike, ' H W'], vaa: Real[ArrayLike, ' H W'], *, scale_deg: float = _ANGLE_SCALE_DEG) -> ViewingGeometry

Per-pixel viewing geometry from the Collection 2 angle bands.

The four Level-1 angle bands (_SZA, _SAA, _VZA, _VAA, computed for band 4) are signed 16-bit integers in hundredths of a degree (LSDS-1822 v7.0, section 3.4), hence scale_deg = 0.01; pass scale_deg=1.0 for angles already in degrees. The angle bands are computed for band 4 (OLI) only, and serve the thermal bands B10 and B11 here as well; a difference between OLI's and TIRS's lines of sight, if any, is not modelled. The relative azimuth

\[ \phi = \arccos\big(\cos(\varphi_v - \varphi_s)\big) \in [0, \pi], \]

sensor minus sun, as :func:~radtrax.geo.sentinel2_geometry forms it, does not depend on whether the azimuths are stored in [−180°, 180°] or [0°, 360°).

Parameters:

Name	Type	Description	Default
sza	Real[ArrayLike, ' H W']	Solar zenith angle, (H, W), in units of scale_deg.	required
saa	Real[ArrayLike, ' H W']	Solar azimuth angle, (H, W).	required
vza	Real[ArrayLike, ' H W']	Viewing (sensor) zenith angle, (H, W).	required
vaa	Real[ArrayLike, ' H W']	Viewing (sensor) azimuth angle, (H, W).	required
scale_deg	float	Degrees per stored unit.	_ANGLE_SCALE_DEG

Returns:

Type	Description
ViewingGeometry	A ViewingGeometry whose leaves are (H, W), built under
ViewingGeometry	eqx.filter_vmap as :func:~radtrax.geo.sentinel2_geometry's.

Raises:

Type	Description
ValueError	If the four bands differ in shape.
RuntimeError	If a cosine is outside (0, 1] (a zenith angle at or beyond 90°, or NaN; equinox.EquinoxRuntimeError, from :class:~radtrax.atmosphere.ViewingGeometry).

Examples:

>>> import numpy as np
>>> sza, saa = np.full((1, 2), 3000), np.array([[15000, -17000]])
>>> vza, vaa = np.zeros((1, 2)), np.array([[-6000, 10000]])
>>> g = landsat_geometry(sza, saa, vza, vaa)
>>> np.rad2deg(np.asarray(g.relative_azimuth_rad)).round(3).tolist()
[[150.0, 90.0]]

Source code in src/radtrax/geo/landsat.py

def landsat_geometry(
    sza: Real[ArrayLike, " H W"],
    saa: Real[ArrayLike, " H W"],
    vza: Real[ArrayLike, " H W"],
    vaa: Real[ArrayLike, " H W"],
    *,
    scale_deg: float = _ANGLE_SCALE_DEG,
) -> ViewingGeometry:
    """Per-pixel viewing geometry from the Collection 2 angle bands.

    The four Level-1 angle bands (``_SZA``, ``_SAA``, ``_VZA``, ``_VAA``,
    computed for band 4) are signed 16-bit integers in hundredths of a
    degree (LSDS-1822 v7.0, section 3.4), hence ``scale_deg = 0.01``; pass
    ``scale_deg=1.0`` for angles already in degrees. The angle bands are
    computed for band 4 (OLI) only, and serve the thermal bands B10 and B11
    here as well; a difference between OLI's and TIRS's lines of sight, if
    any, is not modelled. The relative azimuth

    $$
    \\phi = \\arccos\\big(\\cos(\\varphi_v - \\varphi_s)\\big) \\in [0, \\pi],
    $$

    sensor minus sun, as :func:`~radtrax.geo.sentinel2_geometry` forms it,
    does not depend on whether the azimuths are stored in ``[−180°, 180°]``
    or ``[0°, 360°)``.

    Args:
        sza: Solar zenith angle, ``(H, W)``, in units of ``scale_deg``.
        saa: Solar azimuth angle, ``(H, W)``.
        vza: Viewing (sensor) zenith angle, ``(H, W)``.
        vaa: Viewing (sensor) azimuth angle, ``(H, W)``.
        scale_deg: Degrees per stored unit.

    Returns:
        A ``ViewingGeometry`` whose leaves are ``(H, W)``, built under
        ``eqx.filter_vmap`` as :func:`~radtrax.geo.sentinel2_geometry`'s.

    Raises:
        ValueError: If the four bands differ in shape.
        RuntimeError: If a cosine is outside ``(0, 1]`` (a zenith angle at
            or beyond 90°, or NaN; ``equinox.EquinoxRuntimeError``, from
            :class:`~radtrax.atmosphere.ViewingGeometry`).

    Examples:
        >>> import numpy as np
        >>> sza, saa = np.full((1, 2), 3000), np.array([[15000, -17000]])
        >>> vza, vaa = np.zeros((1, 2)), np.array([[-6000, 10000]])
        >>> g = landsat_geometry(sza, saa, vza, vaa)
        >>> np.rad2deg(np.asarray(g.relative_azimuth_rad)).round(3).tolist()
        [[150.0, 90.0]]
    """
    angles = [np.asarray(x, dtype=float) * scale_deg for x in (sza, saa, vza, vaa)]
    return _geometry_from_degrees("landsat_geometry", *angles)

ecostress

ECOSTRESS: band set, SRF, instrument, noise, radiance and geometry.

ECOSTRESS (PHyTIR, on the ISS) measures five thermal bands, B1 to B5 at about 8.3, 8.8, 9.2, 10.5 and 12.1 µm (radtrax#226). From 15 May 2019 to 28 April 2023 bands 1 and 3 hold fill values only (the three-band mode, ECOSTRESS_3BAND, after the mass-storage anomaly); since 28 April 2023 five-band acquisitions are interleaved with three-band ones, so the band set is a property of the granule (its BandSpecification metadata), not of the date alone. georeader has no ECOSTRESS reader, so this module works on plain inputs, the arrays and metadata an HDF5 reader returns:

• :func:ecostress_bands: the granule's bands from BandSpecification, or from the acquisition date inside the three-band period;
• :func:ecostress_srf and :func:ecostress_instrument: the science team's spectral responses (revision 3), bundled with radtrax (geo/data/ecostress_srf.csv; its header records the source, the SHA-256 of each file and the retrieval date), as a :class:~radtrax.instrument.SpectralResponseFunction and a :class:~radtrax.instrument.BandInstrument;
• :func:ecostress_noise: NEΔT (ECOSTRESS_NETD_K, the < 0.1 K requirement as an upper bound) as an :class:~radtrax.instrument.NEdTNoise;
• :func:ecostress_radiance: the L1B radiance_<n> datasets stacked per band, fill values as NaN;
• :func:ecostress_geometry: a per-pixel :class:~radtrax.atmosphere.ViewingGeometry from the L1B geolocation product's solar and view angles.

Units. L1B radiance is per µm (W m⁻² sr⁻¹ µm⁻¹, ECOSTRESS L1 User Guide, table 3.1.1); the forward model's is per cm⁻¹. ecostress_instrument therefore defaults to per_um=True: it takes radiance per cm⁻¹ and returns band averages per µm, the conversion inside the band integral (:meth:BandInstrument.from_srf <radtrax.instrument.BandInstrument.from_srf>), so its output, the noise built on it and :func:ecostress_radiance share one unit. None of these functions needs the [geo] extra.

ecostress_bands

ecostress_bands(band_specification: Sequence[float] | Float[ArrayLike, ' C'] | None = None, *, acquired: date | str | None = None) -> tuple[str, ...]

The TIR bands with data in a granule.

The L1 products' BandSpecification metadata (group L1B_RADMetadata or L1GEOMetadata) lists six wavelengths [µm], the SWIR band then B1 to B5, with 0 for a band that holds fill values (ECOSTRESS L1 User Guide, D-57151). A band is valid when its entry is not 0. Without the metadata, the acquisition date decides inside the three-band period (ECOSTRESS_3BAND_PERIOD) and before it (five bands); after it the two modes are interleaved, so the date is not enough.

Parameters:

Name	Type	Description	Default
band_specification	Sequence[float] | Float[ArrayLike, ' C'] | None	The BandSpecification values, six numbers (SWIR, B1, ..., B5) [µm], 0 for fill; takes precedence over acquired.	None
acquired	date | str | None	The acquisition date, a date, datetime or ISO string ("2021-07-01").	None

Returns:

Type	Description
str	The valid bands, in order, from ECOSTRESS_BANDS; pass them as
...	bands= to :func:ecostress_instrument.

Raises:

Type	Description
ParameterError	If neither argument is given, or only a date on or after 28 April 2023.
ShapeError	If band_specification does not have six entries.
DegenerateInputError	If the metadata marks every TIR band as fill.

Examples:

>>> ecostress_bands([0.0, 0.0, 8.7, 0.0, 10.5, 12.0])
('B2', 'B4', 'B5')
>>> ecostress_bands(acquired="2018-08-01")
('B1', 'B2', 'B3', 'B4', 'B5')
>>> ecostress_bands(acquired="2021-07-01")
('B2', 'B4', 'B5')

Source code in src/radtrax/geo/ecostress.py

def ecostress_bands(
    band_specification: Sequence[float] | Float[ArrayLike, " C"] | None = None,
    *,
    acquired: dt.date | str | None = None,
) -> tuple[str, ...]:
    """The TIR bands with data in a granule.

    The L1 products' ``BandSpecification`` metadata (group
    ``L1B_RADMetadata`` or ``L1GEOMetadata``) lists six wavelengths [µm],
    the SWIR band then B1 to B5, with 0 for a band that holds fill values
    (ECOSTRESS L1 User Guide, D-57151). A band is valid when its entry is
    not 0. Without the metadata, the acquisition date decides inside the
    three-band period (``ECOSTRESS_3BAND_PERIOD``) and before it (five
    bands); after it the two modes are interleaved, so the date is not
    enough.

    Args:
        band_specification: The ``BandSpecification`` values, six numbers
            (SWIR, B1, ..., B5) [µm], 0 for fill; takes precedence over
            ``acquired``.
        acquired: The acquisition date, a ``date``, ``datetime`` or ISO
            string (``"2021-07-01"``).

    Returns:
        The valid bands, in order, from ``ECOSTRESS_BANDS``; pass them as
        ``bands=`` to :func:`ecostress_instrument`.

    Raises:
        ParameterError: If neither argument is given, or only a date on or
            after 28 April 2023.
        ShapeError: If ``band_specification`` does not have six entries.
        DegenerateInputError: If the metadata marks every TIR band as fill.

    Examples:
        >>> ecostress_bands([0.0, 0.0, 8.7, 0.0, 10.5, 12.0])
        ('B2', 'B4', 'B5')
        >>> ecostress_bands(acquired="2018-08-01")
        ('B1', 'B2', 'B3', 'B4', 'B5')
        >>> ecostress_bands(acquired="2021-07-01")
        ('B2', 'B4', 'B5')
    """
    if band_specification is not None:
        spec = np.asarray(band_specification, dtype=float).reshape(-1)
        if spec.size != 1 + len(ECOSTRESS_BANDS):
            raise ShapeError(
                "ecostress_bands: `band_specification` has six entries (SWIR, "
                f"B1, ..., B5), got {spec.size}."
            )
        bands = tuple(b for b, w in zip(ECOSTRESS_BANDS, spec[1:], strict=True) if w)
        if not bands:
            raise DegenerateInputError(
                "ecostress_bands: every TIR band is fill in `band_specification`."
            )
        return bands
    if acquired is None:
        raise ParameterError(
            "ecostress_bands: pass the granule's `band_specification` or its "
            "`acquired` date."
        )
    day = (
        dt.date.fromisoformat(acquired[:10]) if isinstance(acquired, str) else acquired
    )
    if isinstance(day, dt.datetime):
        day = day.date()
    start, end = ECOSTRESS_3BAND_PERIOD
    if day < start:
        return ECOSTRESS_BANDS
    if day < end:
        return ECOSTRESS_3BAND
    raise ParameterError(
        f"ecostress_bands: from {end.isoformat()} three- and five-band "
        "acquisitions are interleaved; pass the granule's `band_specification`."
    )

ecostress_srf

ecostress_srf(*, bands: Sequence[str] | None = None, wavelengths_hr_nm: Sequence[float] | Float[ArrayLike, ' Ghr'] | None = None) -> SpectralResponseFunction

Build the ECOSTRESS SRF from the bundled science-team responses.

The responses (revision 3, geo/data/ecostress_srf.csv) are published on a wavenumber grid about 9.85 nm apart in wavelength; each band is taken as piecewise linear in wavelength between the samples, placed on one 1 nm grid and projected onto wavelengths_hr_nm with the linear-interpolation (hat) basis, as :func:~radtrax.geo.landsat_srf does. The band average is ∫S L dλ / ∫S dλ, as in the L2 LSTE ATBD (2018, eq. 3).

Parameters:

Name	Type	Description	Default
bands	Sequence[str] | None	Band names, in order, from ECOSTRESS_BANDS; all five by default.	None
wavelengths_hr_nm	Sequence[float] | Float[ArrayLike, ' Ghr'] | None	The high-resolution grid [nm] to integrate on, strictly increasing. Defaults to the 1 nm grid.	None

Returns:

Type	Description
SpectralResponseFunction	A custom SRF whose band centres are the response-weighted mean
SpectralResponseFunction	wavelengths and widths the equivalent widths ∫S dλ / max S.

Warns:

Type	Description
UserWarning	If the grid misses more than 0.1 % of a band's response; the band is then normalised over the part it covers.

Raises:

Type	Description
KeyError	If a band is not one of ECOSTRESS_BANDS.
ValueError	If the grid is not 1-D and strictly increasing.

Examples:

>>> srf = ecostress_srf()
>>> [round(float(c) / 1e3, 2) for c in srf.band_centers_nm]  # µm
[8.29, 8.79, 9.2, 10.49, 12.1]

Source code in src/radtrax/geo/ecostress.py

def ecostress_srf(
    *,
    bands: Sequence[str] | None = None,
    wavelengths_hr_nm: Sequence[float] | Float[ArrayLike, " Ghr"] | None = None,
) -> SpectralResponseFunction:
    """Build the ECOSTRESS SRF from the bundled science-team responses.

    The responses (revision 3, ``geo/data/ecostress_srf.csv``) are
    published on a wavenumber grid about 9.85 nm apart in wavelength; each
    band is taken as piecewise linear in wavelength between the samples,
    placed on one 1 nm grid and projected onto ``wavelengths_hr_nm`` with
    the linear-interpolation (hat) basis, as
    :func:`~radtrax.geo.landsat_srf` does. The band average is
    ``∫S L dλ / ∫S dλ``, as in the L2 LSTE ATBD (2018, eq. 3).

    Args:
        bands: Band names, in order, from ``ECOSTRESS_BANDS``; all five by
            default.
        wavelengths_hr_nm: The high-resolution grid [nm] to integrate on,
            strictly increasing. Defaults to the 1 nm grid.

    Returns:
        A ``custom`` SRF whose band centres are the response-weighted mean
        wavelengths and widths the equivalent widths ``∫S dλ / max S``.

    Warns:
        UserWarning: If the grid misses more than ``0.1 %`` of a band's
            response; the band is then normalised over the part it covers.

    Raises:
        KeyError: If a band is not one of ``ECOSTRESS_BANDS``.
        ValueError: If the grid is not 1-D and strictly increasing.

    Examples:
        >>> srf = ecostress_srf()
        >>> [round(float(c) / 1e3, 2) for c in srf.band_centers_nm]  # µm
        [8.29, 8.79, 9.2, 10.49, 12.1]
    """
    table = _srf_table()
    names = ECOSTRESS_BANDS if bands is None else tuple(bands)
    missing = [b for b in names if b not in table]
    if missing or not names:
        raise MissingKeyError(
            f"ecostress_srf: no band(s) {missing}; bundled: {sorted(table)}."
        )
    samples = [table[b] for b in names]
    return _tabulated_srf("ecostress_srf", samples, names, wavelengths_hr_nm)

ecostress_instrument

ecostress_instrument(nu_cm: Float[ArrayLike, ' G'], *, bands: Sequence[str] | None = None, shift_cm: float = 0.0, per_um: bool = True) -> BandInstrument

ECOSTRESS bands as a :class:~radtrax.instrument.BandInstrument.

The SRFs (:func:ecostress_srf) projected onto the wavelengths of the model grid, λ_g = 10⁷ / ν_g, and turned onto the wavenumber grid by :meth:BandInstrument.from_srf <radtrax.instrument.BandInstrument.from_srf>. With per_um=True (the default) the instrument takes radiance per cm⁻¹ [W m⁻² sr⁻¹ (cm⁻¹)⁻¹] and returns the band average per µm [W m⁻² sr⁻¹ µm⁻¹], the unit of :func:ecostress_radiance:

\[ y_b = \sum_g w_{bg} \frac{\tilde\nu_g^2}{10^4} L_{\tilde\nu}(\tilde\nu_g + \delta). \]

Build it with the granule's bands (:func:ecostress_bands): the band tuple is static, so filter_jit compiles one model per mode, and the three-band instrument's rows are the five-band one's rows of B2, B4, B5.

Parameters:

Name	Type	Description	Default
nu_cm	Float[ArrayLike, ' G']	The model's wavenumber grid [cm⁻¹], strictly increasing, (G,); it should cover each band's response (about 770–1290 cm⁻¹ for all five).	required
bands	Sequence[str] | None	Band names, in order; all five by default.	None
shift_cm	float	Spectral shift δ [cm⁻¹], a differentiable leaf.	0.0
per_um	bool	Band radiance per µm from radiance per cm⁻¹; False gives the band average in the input's units.	True

Returns:

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

Warns:

Type	Description
UserWarning	If the grid misses more than 0.1 % of a band's response (from :func:ecostress_srf).

Raises:

Type	Description
KeyError	If a band is not one of ECOSTRESS_BANDS.
ValueError	If nu_cm is not 1-D and strictly increasing.

Examples:

>>> import numpy as np
>>> nu = np.arange(760.0, 1300.0, 0.5)  # 7.7–13.2 µm
>>> instrument = ecostress_instrument(nu, bands=ECOSTRESS_3BAND)
>>> instrument.band_names
('B2', 'B4', 'B5')
>>> [round(float(y), 9) for y in instrument(3.0 * 1e4 / nu**2)]  # 3 per µm
[3.0, 3.0, 3.0]

Source code in src/radtrax/geo/ecostress.py

def ecostress_instrument(
    nu_cm: Float[ArrayLike, " G"],
    *,
    bands: Sequence[str] | None = None,
    shift_cm: float = 0.0,
    per_um: bool = True,
) -> BandInstrument:
    """ECOSTRESS bands as a :class:`~radtrax.instrument.BandInstrument`.

    The SRFs (:func:`ecostress_srf`) projected onto the wavelengths of the
    model grid, ``λ_g = 10⁷ / ν_g``, and turned onto the wavenumber grid by
    :meth:`BandInstrument.from_srf
    <radtrax.instrument.BandInstrument.from_srf>`. With ``per_um=True`` (the
    default) the instrument takes radiance per cm⁻¹
    [W m⁻² sr⁻¹ (cm⁻¹)⁻¹] and returns the band average per µm
    [W m⁻² sr⁻¹ µm⁻¹], the unit of :func:`ecostress_radiance`:

    $$
    y_b = \\sum_g w_{bg} \\frac{\\tilde\\nu_g^2}{10^4} L_{\\tilde\\nu}(\\tilde\\nu_g + \\delta).
    $$

    Build it with the granule's bands (:func:`ecostress_bands`): the band
    tuple is static, so ``filter_jit`` compiles one model per mode, and the
    three-band instrument's rows are the five-band one's rows of B2, B4, B5.

    Args:
        nu_cm: The model's wavenumber grid [cm⁻¹], strictly increasing,
            ``(G,)``; it should cover each band's response (about
            770–1290 cm⁻¹ for all five).
        bands: Band names, in order; all five by default.
        shift_cm: Spectral shift ``δ`` [cm⁻¹], a differentiable leaf.
        per_um: Band radiance per µm from radiance per cm⁻¹; ``False``
            gives the band average in the input's units.

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

    Warns:
        UserWarning: If the grid misses more than ``0.1 %`` of a band's
            response (from :func:`ecostress_srf`).

    Raises:
        KeyError: If a band is not one of ``ECOSTRESS_BANDS``.
        ValueError: If ``nu_cm`` is not 1-D and strictly increasing.

    Examples:
        >>> import numpy as np
        >>> nu = np.arange(760.0, 1300.0, 0.5)  # 7.7–13.2 µm
        >>> instrument = ecostress_instrument(nu, bands=ECOSTRESS_3BAND)
        >>> instrument.band_names
        ('B2', 'B4', 'B5')
        >>> [round(float(y), 9) for y in instrument(3.0 * 1e4 / nu**2)]  # 3 per µm
        [3.0, 3.0, 3.0]
    """
    nu = np.asarray(nu_cm, dtype=float)
    if nu.ndim != 1 or nu.size < 2 or np.any(np.diff(nu) <= 0.0):
        raise GridError(
            "ecostress_instrument: `nu_cm` must be 1-D, strictly increasing, "
            "with at least two points."
        )
    srf = ecostress_srf(bands=bands, wavelengths_hr_nm=1e7 / nu[::-1])
    return BandInstrument.from_srf(srf, nu, shift_cm=shift_cm, per_um=per_um)

ecostress_noise

ecostress_noise(instrument: BandInstrument, *, netd_k: Mapping[str, float] | None = None, t_ref_k: Real[ArrayLike, ''] = 300.0) -> NEdTNoise

NEΔT as a noise model on instrument's bands.

\[ σ_b = \mathrm{NEΔT}_b \, \frac{\partial \bar B_b}{\partial T}\Big|_{T_\mathrm{ref}}, \qquad S_ε = \mathrm{diag}(σ^2), \]

with \(\bar B_b\) the Planck radiance integrated by instrument. The default NEΔT, ECOSTRESS_NETD_K, is the 0.1 K requirement at 300 K for every band: an upper bound, not an on-orbit measurement. A uniform NEΔT is not a uniform radiance noise: \(\partial \bar B / \partial T\) at 300 K is about 1.5 times larger at 8.3 µm than at 12.1 µm. \(σ\) is in the instrument's units (per µm with the default per_um=True), and is a larger NEΔT in a colder scene (:meth:NEdTNoise.netd_at <radtrax.instrument.NEdTNoise.netd_at>).

Parameters:

Name	Type	Description	Default
instrument	BandInstrument	The bands, from :func:ecostress_instrument; its band_names pick the NEΔT values.	required
netd_k	Mapping[str, float] | None	NEΔT per band [K] at t_ref_k; ECOSTRESS_NETD_K by default.	None
t_ref_k	Real[ArrayLike, '']	The temperature NEΔT is quoted at [K], a scalar.	300.0

Returns:

Name	Type	Description
An	NEdTNoise	class:~radtrax.instrument.NEdTNoise in the "reference" mode.

Raises:

Type	Description
KeyError	If a band has no NEΔT in netd_k.

Examples:

>>> import numpy as np
>>> nu = np.arange(760.0, 1300.0, 0.5)
>>> eco = ecostress_instrument(nu, bands=("B2", "B5"))
>>> sigma = ecostress_noise(eco).radiance_sigma(np.zeros(2))
>>> [f"{float(s):.2e}" for s in sigma]  # W m⁻² sr⁻¹ µm⁻¹
['1.78e-02', '1.20e-02']

Source code in src/radtrax/geo/ecostress.py

def ecostress_noise(
    instrument: BandInstrument,
    *,
    netd_k: Mapping[str, float] | None = None,
    t_ref_k: Real[ArrayLike, ""] = 300.0,
) -> NEdTNoise:
    """NEΔT as a noise model on ``instrument``'s bands.

    $$
    σ_b = \\mathrm{NEΔT}_b \\, \\frac{\\partial \\bar B_b}{\\partial T}\\Big|_{T_\\mathrm{ref}},
    \\qquad S_ε = \\mathrm{diag}(σ^2),
    $$

    with $\\bar B_b$ the Planck radiance integrated by ``instrument``. The
    default NEΔT, ``ECOSTRESS_NETD_K``, is the 0.1 K requirement at 300 K
    for every band: an upper bound, not an on-orbit measurement. A uniform
    NEΔT is not a uniform radiance noise: $\\partial \\bar B / \\partial T$
    at 300 K is about 1.5 times larger at 8.3 µm than at 12.1 µm. $σ$ is in
    the instrument's units (per µm with the default ``per_um=True``), and is
    a larger NEΔT in a colder scene (:meth:`NEdTNoise.netd_at
    <radtrax.instrument.NEdTNoise.netd_at>`).

    Args:
        instrument: The bands, from :func:`ecostress_instrument`; its
            ``band_names`` pick the NEΔT values.
        netd_k: NEΔT per band [K] at ``t_ref_k``; ``ECOSTRESS_NETD_K`` by
            default.
        t_ref_k: The temperature NEΔT is quoted at [K], a scalar.

    Returns:
        An :class:`~radtrax.instrument.NEdTNoise` in the ``"reference"`` mode.

    Raises:
        KeyError: If a band has no NEΔT in ``netd_k``.

    Examples:
        >>> import numpy as np
        >>> nu = np.arange(760.0, 1300.0, 0.5)
        >>> eco = ecostress_instrument(nu, bands=("B2", "B5"))
        >>> sigma = ecostress_noise(eco).radiance_sigma(np.zeros(2))
        >>> [f"{float(s):.2e}" for s in sigma]  # W m⁻² sr⁻¹ µm⁻¹
        ['1.78e-02', '1.20e-02']
    """
    table = ECOSTRESS_NETD_K if netd_k is None else netd_k
    missing = [b for b in instrument.band_names if b not in table]
    if missing:
        raise MissingKeyError(
            f"ecostress_noise: no NEΔT for band(s) {missing}; "
            f"tabulated: {sorted(table)}."
        )
    netd = np.array([table[b] for b in instrument.band_names], dtype=float)
    return NEdTNoise(netd, instrument, t_ref_k=t_ref_k)

ecostress_radiance

ecostress_radiance(radiance: Mapping[str, Real[ArrayLike, ' H W']], *, bands: Sequence[str] | None = None) -> Float[ndarray, ' B H W']

L1B radiance per band and pixel, fill values as NaN.

Stacks the ECO_L1B_RAD datasets radiance_1 to radiance_5 (group Radiance), already calibrated radiance per µm [W m⁻² sr⁻¹ µm⁻¹, ECOSTRESS L1 User Guide, table 3.1.1]; the special values -9997 (pixel not seen), -9998 (missing stripe, not filled in) and -9999 (missing or bad data) become NaN, as does every pixel of a band a three-band granule fills. Model it with ecostress_instrument(..., per_um=True).

Parameters:

Name	Type	Description	Default
radiance	Mapping[str, Real[ArrayLike, ' H W']]	The datasets by name ("radiance_1", ...): an h5py Radiance group or any mapping of name to array, (H, W) each.	required
bands	Sequence[str] | None	Band names, in order; all five by default. Pass :func:ecostress_bands of the granule to skip fill bands.	None

Returns:

Type	Description
Float[ndarray, ' B H W']	Radiance [W m⁻² sr⁻¹ µm⁻¹], float64, (B, H, W).

Raises:

Type	Description
KeyError	If a band is unknown or its dataset is missing.
ShapeError	If the datasets differ in shape.

Examples:

>>> import numpy as np
>>> rad = {
...     "radiance_2": np.array([[8.5, -9998.0]]),
...     "radiance_4": np.array([[9.25, 9.5]]),
... }
>>> ecostress_radiance(rad, bands=("B2", "B4")).tolist()
[[[8.5, nan]], [[9.25, 9.5]]]

Source code in src/radtrax/geo/ecostress.py

def ecostress_radiance(
    radiance: Mapping[str, Real[ArrayLike, " H W"]],
    *,
    bands: Sequence[str] | None = None,
) -> Float[np.ndarray, " B H W"]:
    """L1B radiance per band and pixel, fill values as NaN.

    Stacks the ECO_L1B_RAD datasets ``radiance_1`` to ``radiance_5`` (group
    ``Radiance``), already calibrated radiance per µm
    [W m⁻² sr⁻¹ µm⁻¹, ECOSTRESS L1 User Guide, table 3.1.1]; the special
    values -9997 (pixel not seen), -9998 (missing stripe, not filled in) and
    -9999 (missing or bad data) become NaN, as does every pixel of a band a
    three-band granule fills. Model it with
    ``ecostress_instrument(..., per_um=True)``.

    Args:
        radiance: The datasets by name (``"radiance_1"``, ...): an h5py
            ``Radiance`` group or any mapping of name to array, ``(H, W)``
            each.
        bands: Band names, in order; all five by default. Pass
            :func:`ecostress_bands` of the granule to skip fill bands.

    Returns:
        Radiance [W m⁻² sr⁻¹ µm⁻¹], float64, ``(B, H, W)``.

    Raises:
        KeyError: If a band is unknown or its dataset is missing.
        ShapeError: If the datasets differ in shape.

    Examples:
        >>> import numpy as np
        >>> rad = {
        ...     "radiance_2": np.array([[8.5, -9998.0]]),
        ...     "radiance_4": np.array([[9.25, 9.5]]),
        ... }
        >>> ecostress_radiance(rad, bands=("B2", "B4")).tolist()
        [[[8.5, nan]], [[9.25, 9.5]]]
    """
    names = ECOSTRESS_BANDS if bands is None else tuple(bands)
    unknown = [b for b in names if b not in ECOSTRESS_BANDS]
    if unknown:
        raise MissingKeyError(
            f"ecostress_radiance: unknown band(s) {unknown}; use {ECOSTRESS_BANDS}."
        )
    keys = [f"radiance_{b[1:]}" for b in names]
    absent = [k for k in keys if k not in radiance]
    if absent:
        raise MissingKeyError(f"ecostress_radiance: no dataset(s) {absent}.")
    layers = [np.asarray(radiance[k], dtype=float) for k in keys]
    if len({x.shape for x in layers}) != 1:
        raise ShapeError(
            "ecostress_radiance: the radiance datasets differ in shape "
            f"{[x.shape for x in layers]}."
        )
    stacked = np.stack(layers)
    return np.where(np.isin(stacked, _FILL_VALUES), np.nan, stacked)

ecostress_geometry

ecostress_geometry(solar_zenith: Real[ArrayLike, ' H W'], solar_azimuth: Real[ArrayLike, ' H W'], view_zenith: Real[ArrayLike, ' H W'], view_azimuth: Real[ArrayLike, ' H W']) -> ViewingGeometry

Per-pixel viewing geometry from the L1B geolocation angles.

ECO_L1B_GEO holds solar_zenith, solar_azimuth, view_zenith and view_azimuth [deg] for every radiance pixel (group Geolocation, ECOSTRESS L1 User Guide, table 4.1.1); the ISS's precessing orbit makes them vary from granule to granule, so they are taken per pixel. The relative azimuth

\[ \phi = \arccos\big(\cos(\varphi_v - \varphi_s)\big) \in [0, \pi], \]

sensor minus sun, as :func:~radtrax.geo.landsat_geometry forms it, does not depend on whether the azimuths are stored in [−180°, 180°] or [0°, 360°). The user guide does not say whether view_azimuth points toward or away from the sensor; it is read as toward the sensor, as EMIT's obs file and Landsat's VAA are.

Parameters:

Name	Type	Description	Default
solar_zenith	Real[ArrayLike, ' H W']	Solar zenith angle [deg], (H, W).	required
solar_azimuth	Real[ArrayLike, ' H W']	Solar azimuth angle [deg], (H, W).	required
view_zenith	Real[ArrayLike, ' H W']	Sensor zenith angle [deg], (H, W).	required
view_azimuth	Real[ArrayLike, ' H W']	Sensor azimuth angle [deg], (H, W).	required

Returns:

Type	Description
ViewingGeometry	A ViewingGeometry whose leaves are (H, W), built under
ViewingGeometry	eqx.filter_vmap.

Raises:

Type	Description
ValueError	If the four grids differ in shape.
RuntimeError	If a cosine is outside (0, 1] (a zenith angle at or beyond 90°, or NaN; equinox.EquinoxRuntimeError, from :class:~radtrax.atmosphere.ViewingGeometry).

Examples:

>>> import numpy as np
>>> g = ecostress_geometry(
...     np.full((1, 2), 60.0),
...     np.array([[100.0, -170.0]]),
...     np.full((1, 2), 20.0),
...     np.array([[-80.0, 120.0]]),
... )
>>> np.rad2deg(np.asarray(g.relative_azimuth_rad)).round(3).tolist()
[[180.0, 70.0]]

Source code in src/radtrax/geo/ecostress.py

def ecostress_geometry(
    solar_zenith: Real[ArrayLike, " H W"],
    solar_azimuth: Real[ArrayLike, " H W"],
    view_zenith: Real[ArrayLike, " H W"],
    view_azimuth: Real[ArrayLike, " H W"],
) -> ViewingGeometry:
    """Per-pixel viewing geometry from the L1B geolocation angles.

    ECO_L1B_GEO holds ``solar_zenith``, ``solar_azimuth``, ``view_zenith``
    and ``view_azimuth`` [deg] for every radiance pixel (group
    ``Geolocation``, ECOSTRESS L1 User Guide, table 4.1.1); the ISS's
    precessing orbit makes them vary from granule to granule, so they are
    taken per pixel. The relative azimuth

    $$
    \\phi = \\arccos\\big(\\cos(\\varphi_v - \\varphi_s)\\big) \\in [0, \\pi],
    $$

    sensor minus sun, as :func:`~radtrax.geo.landsat_geometry` forms it,
    does not depend on whether the azimuths are stored in ``[−180°, 180°]``
    or ``[0°, 360°)``. The user guide does not say whether ``view_azimuth``
    points toward or away from the sensor; it is read as toward the sensor,
    as EMIT's obs file and Landsat's ``VAA`` are.

    Args:
        solar_zenith: Solar zenith angle [deg], ``(H, W)``.
        solar_azimuth: Solar azimuth angle [deg], ``(H, W)``.
        view_zenith: Sensor zenith angle [deg], ``(H, W)``.
        view_azimuth: Sensor azimuth angle [deg], ``(H, W)``.

    Returns:
        A ``ViewingGeometry`` whose leaves are ``(H, W)``, built under
        ``eqx.filter_vmap``.

    Raises:
        ValueError: If the four grids differ in shape.
        RuntimeError: If a cosine is outside ``(0, 1]`` (a zenith angle at
            or beyond 90°, or NaN; ``equinox.EquinoxRuntimeError``, from
            :class:`~radtrax.atmosphere.ViewingGeometry`).

    Examples:
        >>> import numpy as np
        >>> g = ecostress_geometry(
        ...     np.full((1, 2), 60.0),
        ...     np.array([[100.0, -170.0]]),
        ...     np.full((1, 2), 20.0),
        ...     np.array([[-80.0, 120.0]]),
        ... )
        >>> np.rad2deg(np.asarray(g.relative_azimuth_rad)).round(3).tolist()
        [[180.0, 70.0]]
    """
    angles = [
        np.asarray(x, dtype=float)
        for x in (solar_zenith, solar_azimuth, view_zenith, view_azimuth)
    ]
    return _geometry_from_degrees("ecostress_geometry", *angles)

matched_filter

GeoTensor operators for the matched filter (behind the [geo] extra).

pipekit Operator wrappers around the one matched-filter core, :mod:radtrax.retrieval.matched_filter — ported from geotoolz's matched_filter operators (radtrax#29), which now delegate here. They keep geotoolz's carrier contract: a GeoTensor in gives a GeoTensor out on the same grid (a plain array gives a plain array), score maps are a new quantity holding NaN at nodata pixels with fill_value_default=NaN, and estimators return non-carrier results (vectors, covariance operators, backgrounds).

Layout: the band axis is -3 (georeader's (C, H, W) / (T, C, H, W)), so these operators move it to the core's band-last (H, W, B) and back. On a (T, C, H, W) stack the background is fitted once over the pixels of every frame and score maps come back (T, 1, H, W).

Nodata: a pixel that is non-finite or equals the carrier's fill_value_default in any band is left out of every statistic.

Precision: the core is JAX; the operators run it in float64 (jax.enable_x64 for the call only), matching geotoolz's NumPy results whatever the global JAX setting. Only the operators do: a covariance operator they return (EstimateCov*, MatchedFilter.cov_op_, StreamingBackgroundResult.cov_op) holds float64 arrays, and using it directly with the core under float32 JAX truncates it.

Conventions that differ from geotoolz's own implementation:

• covariances are the 1/n (maximum-likelihood) scatter, as scikit-learn and the rest of radtrax; geotoolz used 1/(n − 1);
• :class:AdaptiveWindowBackground gives the population variance, takes a guard_window and always reflects at the edges (no pad_mode);
• :class:StreamingBackground names its shrinkage "ledoit_wolf" / "oas" (geotoolz: "shrunk");
• :class:GMMClusterBackground fits scikit-learn's Gaussian mixture.

AdaptiveBackground

Bases: NamedTuple

Local window statistics, in the layout of the cube they came from.

Attributes:

Name	Type	Description
mean	Cube	Per-pixel local mean, the cube's shape.
variance	Cube	Per-pixel per-band local (population) variance, the cube's shape.

StreamingBackgroundResult

Bases: NamedTuple

Background from :class:StreamingBackground.

Attributes:

Name	Type	Description
mean	Spectrum	Mean spectrum (B,).
cov_op	AbstractLinearOperator	Covariance operator (B, B).

MatchedFilter

MatchedFilter(*, mean: SpectrumLike | None = None, cov_op: CovLike | None = None, target: SpectrumLike | None = None, fit_on_call: bool = False, mean_method: MeanMethod = 'median', cov_method: CovMethod = 'ledoit_wolf', axis: int = -3)

Bases: Operator

Per-pixel matched-filter score map.

y = (x − μ)ᵀ Σ⁻¹ t / (tᵀ Σ⁻¹ t) for every pixel, through :func:radtrax.retrieval.matched_filter.apply_image, as a single-band map on the input's carrier.

A fitted operator: :meth:fit estimates whichever of mean / cov_op the constructor left unset (with mean_method / cov_method) into mean_ / cov_op_; :meth:transform scores with that background; calling the operator fits the missing statistics once, on the first call, or — with fit_on_call=True — estimates them from every cube without storing anything. Fitted arrays never reach get_config().

Parameters:

Name	Type	Description	Default
mean	SpectrumLike | None	Background mean spectrum (B,); fitted when None.	None
cov_op	CovLike | None	Background covariance, a lineax operator or a (B, B) matrix (a nested list from a config reload is accepted); fitted when None.	None
target	SpectrumLike | None	Target signature (B,); required before scoring.	None
fit_on_call	bool	Estimate the background from each cube (per scene) instead of fitting the missing statistics once.	False
mean_method	MeanMethod	Mean estimator when fitting. Default "median".	'median'
cov_method	CovMethod	"empirical" (with a 1e-8 ridge), "ledoit_wolf" (default), "oas" or "lowrank".	'ledoit_wolf'
axis	int	Position of the band axis. Default -3.	-3

Attributes:

Name	Type	Description
mean_	Spectrum | None	Fitted mean, or None.
cov_op_	CovLike | None	Fitted covariance operator, or None.

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self,
    *,
    mean: SpectrumLike | None = None,
    cov_op: CovLike | None = None,
    target: SpectrumLike | None = None,
    fit_on_call: bool = False,
    mean_method: MeanMethod = "median",
    cov_method: CovMethod = "ledoit_wolf",
    axis: int = -3,
) -> None:
    self.mean = None if mean is None else np.asarray(mean, dtype=float)
    self.cov_op = _as_cov(cov_op)
    self.target = None if target is None else np.asarray(target, dtype=float)
    self.fit_on_call = fit_on_call
    self.mean_method = mean_method
    self.cov_method = cov_method
    self.axis = axis
    self.mean_: Spectrum | None = None
    self.cov_op_: CovLike | None = None

fit

fit(gt: Carrier) -> MatchedFilter

Estimate the unset background statistics from the valid pixels of gt.

Returns:

Type	Description
MatchedFilter	This operator, for chaining.

Source code in src/radtrax/geo/matched_filter.py

@_float64
def fit(self, gt: Carrier) -> MatchedFilter:
    """Estimate the unset background statistics from the valid pixels of ``gt``.

    Returns:
        This operator, for chaining.
    """
    samples, _, _ = _samples(gt, self.axis, type(self).__name__)
    mean, cov_op = self._estimate(samples, keep_supplied=True)
    self.cov_op_ = cov_op
    self.mean_ = mean  # published last: the fit-once check reads mean_
    return self

transform

transform(gt: Carrier) -> ScoreMap

Score gt against the fitted (or supplied) background.

Raises:

Type	Description
ValueError	If no background is available or target is unset.

Source code in src/radtrax/geo/matched_filter.py

@_float64
def transform(self, gt: Carrier) -> ScoreMap:
    """Score ``gt`` against the fitted (or supplied) background.

    Raises:
        ValueError: If no background is available or ``target`` is unset.
    """
    mean, cov_op = self._background()
    return self._score(gt, mean, cov_op)

MatchedFilterPixel

MatchedFilterPixel(*, mean: SpectrumLike, cov_op: CovLike, target: SpectrumLike)

Bases: Operator

Matched-filter score of one pixel spectrum, as a float.

Parameters:

Name	Type	Description	Default
mean	SpectrumLike	Background mean spectrum (B,).	required
cov_op	CovLike	Background covariance, a lineax operator or a (B, B) matrix.	required
target	SpectrumLike	Target signature (B,).	required

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self, *, mean: SpectrumLike, cov_op: CovLike, target: SpectrumLike
) -> None:
    self.mean = np.asarray(mean, dtype=float)
    self.cov_op = _as_cov(cov_op)
    self.target = np.asarray(target, dtype=float)

MatchedFilterSNR

MatchedFilterSNR(*, amplitude: float, cov_op: CovLike, target: SpectrumLike)

Bases: Operator

Theoretical detection SNR amplitude · √(tᵀ Σ⁻¹ t).

A graph source: computed from its configuration; an optional input is accepted and ignored, so it can close a Sequential.

Parameters:

Name	Type	Description	Default
amplitude	float	Target amplitude.	required
cov_op	CovLike	Background covariance, a lineax operator or a (B, B) matrix.	required
target	SpectrumLike	Target signature (B,).	required

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self, *, amplitude: float, cov_op: CovLike, target: SpectrumLike
) -> None:
    self.amplitude = amplitude
    self.cov_op = _as_cov(cov_op)
    self.target = np.asarray(target, dtype=float)

DetectionThreshold

DetectionThreshold(*, false_alarm_rate: float, cov_op: CovLike, target: SpectrumLike)

Bases: Operator

Score threshold with a given false-alarm rate under the Gaussian null.

A graph source, like :class:MatchedFilterSNR.

Parameters:

Name	Type	Description	Default
false_alarm_rate	float	Exceedance probability, in (0, 1).	required
cov_op	CovLike	Background covariance, a lineax operator or a (B, B) matrix.	required
target	SpectrumLike	Target signature (B,).	required

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self, *, false_alarm_rate: float, cov_op: CovLike, target: SpectrumLike
) -> None:
    self.false_alarm_rate = false_alarm_rate
    self.cov_op = _as_cov(cov_op)
    self.target = np.asarray(target, dtype=float)

ValidateMFInputs

ValidateMFInputs(*, cov_op: CovLike, target: SpectrumLike)

Bases: Operator

Pass-through guard: fail early on a zero target or a singular covariance.

Returns its input unchanged, whatever its type.

Parameters:

Name	Type	Description	Default
cov_op	CovLike	Background covariance, a lineax operator or a (B, B) matrix.	required
target	SpectrumLike	Target signature (B,).	required

Source code in src/radtrax/geo/matched_filter.py

def __init__(self, *, cov_op: CovLike, target: SpectrumLike) -> None:
    self.cov_op = _as_cov(cov_op)
    self.target = np.asarray(target, dtype=float)

EstimateMean

EstimateMean(*, method: MeanMethod = 'mean', trim_proportion: float = 0.1, huber_c: float = 1.345, axis: int = -3)

Bases: Operator

Background mean spectrum (B,) of a cube's valid pixels.

Parameters:

Name	Type	Description	Default
method	MeanMethod	"mean" (default), "median", "trimmed" or "huber".	'mean'
trim_proportion	float	Fraction cut from each tail for "trimmed".	0.1
huber_c	float	Huber threshold in MAD units for "huber".	1.345
axis	int	Position of the band axis. Default -3.	-3

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self,
    *,
    method: MeanMethod = "mean",
    trim_proportion: float = 0.1,
    huber_c: float = 1.345,
    axis: int = -3,
) -> None:
    self.method = method
    self.trim_proportion = trim_proportion
    self.huber_c = huber_c
    self.axis = axis

EstimateCovEmpirical

EstimateCovEmpirical(*, mean: SpectrumLike | None = None, ridge: float = 0.0, axis: int = -3)

Bases: Operator

Empirical (1/n) covariance operator of a cube's valid pixels.

Parameters:

Name	Type	Description	Default
mean	SpectrumLike | None	Mean spectrum (B,) to centre on; the sample mean when None.	None
ridge	float	Diagonal ridge added to the covariance.	0.0
axis	int	Position of the band axis. Default -3.	-3

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self, *, mean: SpectrumLike | None = None, ridge: float = 0.0, axis: int = -3
) -> None:
    self.mean = None if mean is None else np.asarray(mean, dtype=float)
    self.ridge = ridge
    self.axis = axis

EstimateCovShrunk

EstimateCovShrunk(*, mean: SpectrumLike | None = None, method: ShrinkMethod = 'ledoit_wolf', axis: int = -3)

Bases: Operator

Ledoit–Wolf or OAS covariance operator of a cube's valid pixels.

Parameters:

Name	Type	Description	Default
mean	SpectrumLike | None	Mean spectrum (B,) to centre on; the sample mean when None.	None
method	ShrinkMethod	"ledoit_wolf" (default) or "oas".	'ledoit_wolf'
axis	int	Position of the band axis. Default -3.	-3

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self,
    *,
    mean: SpectrumLike | None = None,
    method: ShrinkMethod = "ledoit_wolf",
    axis: int = -3,
) -> None:
    self.mean = None if mean is None else np.asarray(mean, dtype=float)
    self.method = method
    self.axis = axis

EstimateCovLowRank

EstimateCovLowRank(*, mean: SpectrumLike | None = None, rank: int = LOWRANK_RANK, tikhonov: float = LOWRANK_TIKHONOV, seed: int | None = 0, n_oversamples: int = 10, axis: int = -3)

Bases: Operator

Low-rank plus Tikhonov covariance λI + V D Vᵀ of a cube's valid pixels.

Parameters:

Name	Type	Description	Default
mean	SpectrumLike | None	Mean spectrum (B,) to centre on; the sample mean when None.	None
rank	int	Leading components kept.	LOWRANK_RANK
tikhonov	float	Diagonal floor λ > 0.	LOWRANK_TIKHONOV
seed	int | None	Seed of the randomised SVD.	0
n_oversamples	int	Extra random directions beyond rank.	10
axis	int	Position of the band axis. Default -3.	-3

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self,
    *,
    mean: SpectrumLike | None = None,
    rank: int = LOWRANK_RANK,
    tikhonov: float = LOWRANK_TIKHONOV,
    seed: int | None = 0,
    n_oversamples: int = 10,
    axis: int = -3,
) -> None:
    self.mean = None if mean is None else np.asarray(mean, dtype=float)
    self.rank = rank
    self.tikhonov = tikhonov
    self.seed = seed
    self.n_oversamples = n_oversamples
    self.axis = axis

GMMClusterBackground

GMMClusterBackground(*, n_clusters: int, cov_estimator: Literal['empirical', 'ledoit_wolf', 'oas'] = 'ledoit_wolf', seed: int | None = 0, bayesian: bool = False, axis: int = -3)

Bases: Operator

Per-cluster backgrounds from a Gaussian mixture over the valid pixels.

Returns a :class:~radtrax.retrieval.matched_filter.ClusterBackground whose labels have the cube's spatial shape (-1 at nodata pixels); wire it into :class:ApplyClusterMF.

Parameters:

Name	Type	Description	Default
n_clusters	int	Number of mixture components.	required
cov_estimator	Literal['empirical', 'ledoit_wolf', 'oas']	Per-cluster covariance, "empirical", "ledoit_wolf" (default) or "oas".	'ledoit_wolf'
seed	int | None	Seed of the mixture fit.	0
bayesian	bool	Use a Bayesian mixture, which prunes unused components.	False
axis	int	Position of the band axis. Default -3.	-3

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self,
    *,
    n_clusters: int,
    cov_estimator: Literal["empirical", "ledoit_wolf", "oas"] = "ledoit_wolf",
    seed: int | None = 0,
    bayesian: bool = False,
    axis: int = -3,
) -> None:
    self.n_clusters = n_clusters
    self.cov_estimator = cov_estimator
    self.seed = seed
    self.bayesian = bayesian
    self.axis = axis

ApplyClusterMF

ApplyClusterMF(*, target: SpectrumLike, axis: int = -3)

Bases: Operator

Matched filter with each pixel's own cluster background.

The background arrives at call time, so the config (target, axis) stays YAML-safe: ApplyClusterMF(target=t)(gt, background) with the background from :class:GMMClusterBackground.

Parameters:

Name	Type	Description	Default
target	SpectrumLike	Target signature (B,).	required
axis	int	Position of the band axis. Default -3.	-3

Source code in src/radtrax/geo/matched_filter.py

def __init__(self, *, target: SpectrumLike, axis: int = -3) -> None:
    self.target = np.asarray(target, dtype=float)
    self.axis = axis

AdaptiveWindowBackground

AdaptiveWindowBackground(*, window: int = 7, guard_window: int = 0, axis: int = -3)

Bases: Operator

Per-pixel local mean and variance over a square window.

Returns an :class:AdaptiveBackground of two cubes in the input's layout (not a carrier); nodata pixels carry no weight and are NaN in both.

Parameters:

Name	Type	Description	Default
window	int	Side length of the window, an odd integer ≥ 3. Default 7.	7
guard_window	int	Side length of a central guard region left out of each window (0, the default, keeps the pixel itself; 1 drops it, so a plume does not bias its own background).	0
axis	int	Position of the band axis. Default -3.	-3

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self, *, window: int = 7, guard_window: int = 0, axis: int = -3
) -> None:
    self.window = window
    self.guard_window = guard_window
    self.axis = axis

ApplyAdaptiveMF

ApplyAdaptiveMF(*, target: SpectrumLike, ridge: float = 0.0, axis: int = -3)

Bases: Operator

Matched filter against each pixel's local diagonal background.

α = tᵀ Λ⁻¹ (x − μ) / (tᵀ Λ⁻¹ t) with Λ = diag(σ² + ridge) from :class:AdaptiveWindowBackground, passed at call time: ApplyAdaptiveMF(target=t)(gt, background).

Parameters:

Name	Type	Description	Default
target	SpectrumLike	Target signature (B,).	required
ridge	float	Non-negative variance floor. Default 0.	0.0
axis	int	Position of the band axis. Default -3.	-3

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self, *, target: SpectrumLike, ridge: float = 0.0, axis: int = -3
) -> None:
    self.target = np.asarray(target, dtype=float)
    self.ridge = ridge
    self.axis = axis

StreamingBackground

StreamingBackground(*, cov_kind: Literal['empirical', 'ledoit_wolf', 'oas'] = 'ledoit_wolf', ridge: float = 1e-08, axis: int = -3)

Bases: Operator

Mean and covariance accumulated over an iterable of cubes.

Streams the valid pixels of every cube (GeoTensor or array) through :func:radtrax.retrieval.matched_filter.streaming_background; the cubes arrive at call time, so the config stays YAML-safe.

Parameters:

Name	Type	Description	Default
cov_kind	Literal['empirical', 'ledoit_wolf', 'oas']	"ledoit_wolf" (default), "oas" or "empirical"; the shrunk forms equal scikit-learn's on all the pixels at once.	'ledoit_wolf'
ridge	float	Diagonal ridge added to the covariance. Default 1e-8.	1e-08
axis	int	Position of the band axis in each cube. Default -3.	-3

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self,
    *,
    cov_kind: Literal["empirical", "ledoit_wolf", "oas"] = "ledoit_wolf",
    ridge: float = 1e-8,
    axis: int = -3,
) -> None:
    self.cov_kind = cov_kind
    self.ridge = ridge
    self.axis = axis

LinearTargetFromObs

LinearTargetFromObs(*, obs_model: Any, vmr_background: VMRField | Cube | None = None, pattern: str | VMRField | Cube = 'uniform', pixel: Pixel | None = None, axis: int = -3, fd_step: float | None = None)

Bases: Operator

Tangent-linear target signature (B,) from an observation model.

With a :class:~radtrax.retrieval.var3d.RadianceObservationModel the target is the exact JVP of its forward model (:func:radtrax.retrieval.matched_filter.linear_target_from_obs) about vmr_background (default: the model's reference VMR on the cube's grid; pass it on the high-resolution grid when the model has a GSD). Any other callable maps a state cube to a response with its band axis at axis (or to a (B,) vector), about vmr_background (default: the cube itself, its nodata pixels set to the band means of the valid ones); a cube response is averaged over the valid pixels. A JAX-traceable callable is linearised exactly with :func:jax.jvp; one that is not (NumPy code converting its input) by a central difference, which is only as good as its step: fd_step when given, else 1e-4 times the state's largest magnitude (1e-4 itself for an all-zero state, where the state carries no scale — set fd_step from the model's then). (geotoolz always took a one-sided step of 1e-4 in the state's own units, far outside the linear regime for a VMR of ~2e-6.)

Parameters:

Name	Type	Description	Default
obs_model	Any	An observation model, or a callable.	required
vmr_background	VMRField | Cube | None	Background state; see above.	None
pattern	str | VMRField | Cube	"uniform" (default), "impulse" or an explicit perturbation array.	'uniform'
pixel	Pixel | None	Impulse location on the state's grid; an observation model reads the target at the output pixel holding it (pixel // f for a GSD of factor f).	None
axis	int	Position of the band axis. Default -3.	-3
fd_step	float | None	Absolute central-difference step for a callable that JAX cannot trace; see above.	None

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self,
    *,
    obs_model: Any,
    vmr_background: VMRField | Cube | None = None,
    pattern: str | VMRField | Cube = "uniform",
    pixel: Pixel | None = None,
    axis: int = -3,
    fd_step: float | None = None,
) -> None:
    self.obs_model = obs_model
    self.vmr_background = vmr_background
    self.pattern = pattern
    self.pixel = pixel
    self.axis = axis
    self.fd_step = fd_step

NonlinearTargetFromObs

NonlinearTargetFromObs(*, obs_model: Any, vmr_background: VMRField | Cube | None = None, amplitude: float = 1.0, pattern: str | VMRField | Cube = 'uniform', pixel: Pixel | None = None, axis: int = -3)

Bases: Operator

Finite-amplitude target H(x_b + α δx) − H(x_b), (B,).

Keeps the forward model's nonlinearity at amplitude α. With a :class:~radtrax.retrieval.var3d.RadianceObservationModel it is :func:radtrax.retrieval.matched_filter.nonlinear_target_from_obs about vmr_background (default: the model's reference VMR on the cube's grid). Any other callable maps a state cube to a response with its band axis at axis (or to a (B,) vector), about vmr_background (default: the cube, nodata pixels set to the valid band means); a cube response is averaged over the valid pixels.

Parameters:

Name	Type	Description	Default
obs_model	Any	An observation model, or a callable.	required
vmr_background	VMRField | Cube | None	Background state; see above (the incoming cube for a callable).	None
amplitude	float	Perturbation amplitude α. Default 1.	1.0
pattern	str | VMRField | Cube	"uniform" (default), "impulse" or an explicit perturbation array.	'uniform'
pixel	Pixel | None	Impulse location for "impulse".	None
axis	int	Position of the band axis. Default -3.	-3

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self,
    *,
    obs_model: Any,
    vmr_background: VMRField | Cube | None = None,
    amplitude: float = 1.0,
    pattern: str | VMRField | Cube = "uniform",
    pixel: Pixel | None = None,
    axis: int = -3,
) -> None:
    self.obs_model = obs_model
    self.vmr_background = vmr_background
    self.amplitude = amplitude
    self.pattern = pattern
    self.pixel = pixel
    self.axis = axis

UnitAbsorptionTarget

UnitAbsorptionTarget(*, model: ForwardModel, atmosphere: AtmosphericState, surface: Surface, geometry: ViewingGeometry, gas: str = 'CH4', profile: Float[ndarray, ' L'] | None = None, log_cube: bool = False, mean_method: MeanMethod = 'median', axis: int = -3)

Bases: Operator

Target (B,) from v0's layered :class:~radtrax.ForwardModel.

The scene's unit absorption spectrum t = ∂ ln y / ∂ΔX at ΔX = 0 [per mol/m²], from :func:~radtrax.unit_absorption_target on the background atmosphere, surface and geometry (one geometry for the scene, e.g. its mean angles), turned into the target of the cube it is called on:

• a cube of radiance y (default): μ ⊙ t, with μ the cube's background mean (:class:EstimateMean, mean_method), since ∂y/∂ΔX = y t; the matched filter then scores ΔX in mol/m²;
• a cube of ln y (log_cube=True): t itself.

The single-slab :class:~radtrax.retrieval.var3d.RadianceObservationModel goes to :class:LinearTargetFromObs instead (radtrax#98).

Parameters:

Name	Type	Description	Default
model	ForwardModel	The forward model; its instrument's bands are the cube's bands, in order.	required
atmosphere	AtmosphericState	The background atmospheric state (state of :func:~radtrax.unit_absorption_target; pipekit's operators reserve state).	required
surface	Surface	The surface (a Lambertian albedo cancels from t).	required
geometry	ViewingGeometry	Sun and sensor directions, one for the scene.	required
gas	str	The plume's gas. Default "CH4".	'CH4'
profile	Float[ndarray, ' L'] | None	Relative weight of the plume per layer, (L,); the bottom layer when None.	None
log_cube	bool	Whether the cube holds ln y rather than y.	False
mean_method	MeanMethod	Mean estimator for μ. Default "median".	'median'
axis	int	Position of the band axis. Default -3.	-3

Raises:

Type	Description
ShapeError	(on call) If the model's band count differs from the cube's.

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self,
    *,
    model: ForwardModel,
    atmosphere: AtmosphericState,
    surface: Surface,
    geometry: ViewingGeometry,
    gas: str = "CH4",
    profile: Float[np.ndarray, " L"] | None = None,
    log_cube: bool = False,
    mean_method: MeanMethod = "median",
    axis: int = -3,
) -> None:
    self.model = model
    self.atmosphere = atmosphere
    self.surface = surface
    self.geometry = geometry
    self.gas = gas
    self.profile = None if profile is None else np.asarray(profile, dtype=float)
    self.log_cube = log_cube
    self.mean_method = mean_method
    self.axis = axis

ColumnEnhancement

ColumnEnhancement(*, gas: str = 'CH4', sensor: str = 'EMIT', obs_model: Any | None = None, mean_method: MeanMethod = 'median', cov_method: ShrinkMethod | Literal['lowrank'] = 'ledoit_wolf', target_pattern: str = 'uniform', vmr_background: VMRField | Cube | None = None, target_pixel: Pixel | None = None, fd_step: float | None = None, rank: int | None = None, tikhonov: float | None = None, axis: int = -3)

Bases: Operator

Trace-gas enhancement map: mean, covariance, target and matched filter.

Chains :class:EstimateMean, :class:EstimateCovShrunk (or :class:EstimateCovLowRank), :class:LinearTargetFromObs (a uniform unit target without obs_model) and :class:MatchedFilter.

Parameters:

Name	Type	Description	Default
gas	str	Gas label (metadata only). Default "CH4".	'CH4'
sensor	str	Sensor label (metadata only). Default "EMIT".	'EMIT'
obs_model	Any | None	Observation model or callable for the target.	None
mean_method	MeanMethod	Mean estimator. Default "median".	'median'
cov_method	ShrinkMethod | Literal['lowrank']	"ledoit_wolf" (default), "oas" or "lowrank".	'ledoit_wolf'
target_pattern	str	Perturbation pattern of the target.	'uniform'
vmr_background	VMRField | Cube | None	Background state of the target, forwarded to :class:LinearTargetFromObs; pass it on the model's high-resolution grid when the observation model has a GSD.	None
target_pixel	Pixel | None	The target's pixel: impulse location on the state's grid, read at the output pixel holding it.	None
fd_step	float | None	Central-difference step of the target for an obs_model JAX cannot trace, forwarded to :class:LinearTargetFromObs.	None
rank	int | None	Rank for "lowrank" (default 10 when None).	None
tikhonov	float | None	Floor for "lowrank" (default 1e-3 when None).	None
axis	int	Position of the band axis. Default -3.	-3

Source code in src/radtrax/geo/matched_filter.py

def __init__(
    self,
    *,
    gas: str = "CH4",
    sensor: str = "EMIT",
    obs_model: Any | None = None,
    mean_method: MeanMethod = "median",
    cov_method: ShrinkMethod | Literal["lowrank"] = "ledoit_wolf",
    target_pattern: str = "uniform",
    vmr_background: VMRField | Cube | None = None,
    target_pixel: Pixel | None = None,
    fd_step: float | None = None,
    rank: int | None = None,
    tikhonov: float | None = None,
    axis: int = -3,
) -> None:
    self.gas = gas
    self.sensor = sensor
    self.obs_model = reject_config_summary(obs_model, "obs_model")
    self.mean_method = mean_method
    self.cov_method = cov_method
    self.target_pattern = target_pattern
    self.vmr_background = (
        None if vmr_background is None else np.asarray(vmr_background, dtype=float)
    )
    self.target_pixel = None if target_pixel is None else tuple(target_pixel)
    self.fd_step = fd_step
    self.rank = rank
    self.tikhonov = tikhonov
    self.axis = axis

mbmp

The multi-pass operators on GeoTensors (behind the [geo] extra).

:class:MBMP is ported from geotoolz's plume.SBMP (radtrax#30), renamed after radtrax · 2: with a reference scene it computes the multi-band multi-pass log signal, :func:radtrax.retrieval.mbmp.mbmp_log_signal, or with scale_factors=True :func:radtrax.retrieval.mbmp.mbmp_signal; without one, the single-pass prefilter :func:radtrax.retrieval.mbmp.swir_contrast. :class:SBMP (one band, two passes) and :class:MBSP (two bands, one pass) complete radtrax · 2's family (radtrax#111). Their scale factors are fitted on every valid pixel of each pass (each time step of a stack); a plume in the fit biases them, so fit on a background mask with the functions in :mod:radtrax.retrieval.mbmp when a plume is large.

MBMP

MBMP(*, swir1: BandRef = 'B11', swir2: BandRef = 'B12', axis: int = -3, eps: float = 1e-10, scale_factors: bool = False)

Bases: Operator

Sentinel-2 multi-band multi-pass methane signal, with or without scale factors.

Called as MBMP()(scene, reference_scene): band 11 against band 12 on the observed and on a clean-air reference pass of the same pixel grid,

.. math:: s = \ln\frac{R_{11} + ε}{R_{12} + ε} - \ln\frac{R'{11} + ε}{R' ,} + ε

positive where a plume depresses band 12 (Varon et al., 2021). A scene-wide change between the passes stays in :math:s as an offset unless scale_factors subtracts ln(c / c′), each pass's factor fitted on its valid pixels (geotoolz called the signal SBMP; radtrax · 2 shows it is an MBMP). Without a reference scene the output is the single-pass prefilter (R₁₁ − R₁₂) / (R₁₁ + R₁₂ + ε), not a column retrieval. Negative values are clipped to zero first.

With a reference, a pixel whose band 11 or 12 is not positive on either pass (shadow, water, noise, or a -9999 sentinel in an unnamed reference) is NaN: its log would saturate at about ±ln(R/ε) ≈ ±22 and outrank every real plume (0.01–0.2). geotoolz's SBMP scored those pixels.

Bands resolve through the carrier's names (attrs band_names, descriptions, bands); only a plain array with exactly 12 bands falls back to the Sentinel-2 L2A order (B1–B8, B8A, B9, B11, B12), so "B11" never picks B10 on a 13-band L1C stack. A GeoTensor reference resolves its own names; a plain one shares the scene's band positions.

The output is single-band ((T, 1, H, W) for a stack), NaN where either scene has a nodata pixel, with fill_value_default=NaN. The band arithmetic runs in float64.

Parameters:

Name	Type	Description	Default
swir1	BandRef	Band 11 (weak absorption), an index or a name. Default "B11".	'B11'
swir2	BandRef	Band 12 (strong absorption), an index or a name. Default "B12".	'B12'
axis	int	Band axis. Default -3.	-3
eps	float	Guard against log and division by zero. Default 1e-10.	1e-10
scale_factors	bool	Fit and remove each pass's scale factor (radtrax#111). Default False, the signal as before.	False

Source code in src/radtrax/geo/mbmp.py

def __init__(
    self,
    *,
    swir1: BandRef = "B11",
    swir2: BandRef = "B12",
    axis: int = -3,
    eps: float = 1e-10,
    scale_factors: bool = False,
) -> None:
    self.swir1 = swir1
    self.swir2 = swir2
    self.axis = axis
    self.eps = eps
    self.scale_factors = scale_factors

SBMP

SBMP(*, swir2: BandRef = 'B12', axis: int = -3, eps: float = 1e-10)

Bases: Operator

Sentinel-2 single-band multi-pass signal, −ln(c_s R₁₂ / R′₁₂).

Called as SBMP()(scene, reference_scene): band 12 on the observed and on a reference pass of the same grid, with R′₁₂ ≈ c_s R₁₂ fitted on the valid pixels of each time step (:func:radtrax.retrieval.mbmp.sbmp_signal). Invert through τ₁₂. A scene-wide change between passes cancels; a local one (a wetted field) reads as methane (radtrax · 2). Not geotoolz's SBMP, which was the MBMP log signal (:class:MBMP).

Band resolution, NaN for nodata or non-positive pixels, the single-band output and float64 arithmetic are as in :class:MBMP.

Parameters:

Name	Type	Description	Default
swir2	BandRef	Band 12 (strong absorption), an index or a name. Default "B12".	'B12'
axis	int	Band axis. Default -3.	-3
eps	float	Guard against log of zero. Default 1e-10.	1e-10

Source code in src/radtrax/geo/mbmp.py

def __init__(
    self, *, swir2: BandRef = "B12", axis: int = -3, eps: float = 1e-10
) -> None:
    self.swir2 = swir2
    self.axis = axis
    self.eps = eps

MBSP

MBSP(*, swir1: BandRef = 'B11', swir2: BandRef = 'B12', axis: int = -3, eps: float = 1e-10)

Bases: Operator

Sentinel-2 multi-band single-pass signal, −ln(c R₁₂ / R₁₁).

Called as MBSP()(scene): band 11 against band 12 on one pass, with R₁₁ ≈ c R₁₂ fitted on the valid pixels of each time step (:func:radtrax.retrieval.mbmp.mbsp_signal). Invert through τ₁₂ − τ₁₁. A surface whose band ratio departs from the scene's reads as methane (radtrax · 2).

Band resolution, NaN for nodata or non-positive pixels, the single-band output and float64 arithmetic are as in :class:MBMP.

Parameters:

Name	Type	Description	Default
swir1	BandRef	Band 11 (weak absorption), an index or a name. Default "B11".	'B11'
swir2	BandRef	Band 12 (strong absorption), an index or a name. Default "B12".	'B12'
axis	int	Band axis. Default -3.	-3
eps	float	Guard against log of zero. Default 1e-10.	1e-10

Source code in src/radtrax/geo/mbmp.py

def __init__(
    self,
    *,
    swir1: BandRef = "B11",
    swir2: BandRef = "B12",
    axis: int = -3,
    eps: float = 1e-10,
) -> None:
    self.swir1 = swir1
    self.swir2 = swir2
    self.axis = axis
    self.eps = eps