Slab models

plumax's single-slab Beer–Lambert models (exact, Maclaurin, Taylor) with hand-written Jacobians, the normalised-brightness table, and the observation configuration. Kept as the regression oracle for the layered v0 forward model.

slab

plumax's single-slab Beer–Lambert models, kept as the regression oracle.

One homogeneous slab at one (p, T): forward holds the exact, Maclaurin and Taylor forward models with their hand-written Jacobians, nb_lut the normalised-brightness table for fast per-pixel plume injection, and config the observation geometry, instrument spec and unit constants. The layered v0 forward model (Wave 3) must reproduce these at L = 1 (gate V1).

Ported from plumax 0.1.2 at 28f257e (plumax.radtran, the forward half; its matched-filter half was retired in favour of radtrax.retrieval.matched_filter).

config

Observation-geometry and instrument configuration for radtran.

The forward models in :mod:radtrax.slab.forward need three groups of inputs: atmospheric state (T, P, background VMR), a gas-species LUT (σ(ν, T, P), from :mod:radtrax.spectroscopy.hitran), and the viewing geometry and instrument response.

This module factors the latter two into frozen dataclasses:

• :class:ObservationGeometry — solar/viewing zenith angles, vertical path length, optional AMF override.
• :class:InstrumentSpec — band centres and widths for a synthetic multispectral instrument. Real instruments should build a :class:~radtrax.instrument.srf.SpectralResponseFunction from measured SRFs; InstrumentSpec is a convenient constructor for the tutorial notebooks.

The constants, the number density and the air-mass factor have one home, :mod:radtrax.atmosphere (radtrax#27); this module re-exports the constants and wraps the number density with scalar validation.

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

ObservationGeometry dataclass

ObservationGeometry(sza_deg: float = 30.0, vza_deg: float = 0.0, path_length_cm: float = 840000.0, amf: float | None = None)

Viewing geometry for a single pixel.

Attributes

sza_deg : float Solar zenith angle [degrees], 0 (overhead sun) to ~90 (grazing). vza_deg : float Viewing zenith angle [degrees], 0 (nadir) to ~90 (limb). path_length_cm : float Vertical atmospheric path [cm] traversed by the line of sight before AMF correction. For a single-slab column over the full atmosphere, a common choice is L = p_surf / (ρ · g) ≈ 8.4 km = 8.4e5 cm (scale height of dry air at 1 atm, 288 K). amf : float, optional Air-mass factor. If None the plane-parallel two-way form 1/cos(SZA) + 1/cos(VZA) is used.

Notes

The plane-parallel AMF breaks down past ~75° zenith; beyond that a spherical correction (Chapman function) is needed. For the tutorial regime (SZA, VZA < 60°) the plane-parallel form is within a few percent of the exact value.

air_mass_factor property

air_mass_factor: float

Plane-parallel two-way AMF, or the explicit override if provided.

InstrumentSpec dataclass

InstrumentSpec(name: str = 'synthetic', band_centers_nm: ndarray = (lambda: array([1610.0, 2190.0], dtype=float))(), band_widths_nm: ndarray = (lambda: array([90.0, 180.0], dtype=float))(), band_names: tuple[str, ...] = ('B11', 'B12'), srf_type: str = 'gaussian')

Multispectral instrument description.

Attributes

name : str Free-text instrument identifier for dataset metadata. band_centers_nm : np.ndarray Band centre wavelengths [nm], shape (n_bands,). band_widths_nm : np.ndarray Full-width at half-maximum of each band [nm], shape (n_bands,). band_names : tuple of str Short names for each band (e.g. ('B11', 'B12')). Length matches band_centers_nm. srf_type : str Shape of the synthetic SRF built by :meth:InstrumentSpec.make_srf — one of 'gaussian', 'rectangular', 'triangular'.

Examples

Sentinel-2 B11/B12-like synthetic instrument::

InstrumentSpec(
    name="S2A-like",
    band_centers_nm=np.array([1610.0, 2190.0]),
    band_widths_nm=np.array([90.0, 180.0]),
    band_names=("B11", "B12"),
)

n_bands property

n_bands: int

Number of bands.

make_srf

make_srf(wavelengths_nm: Float[ArrayLike, ' Ghr']) -> SpectralResponseFunction

Construct a :class:SpectralResponseFunction over wavelengths_nm.

Imported locally to avoid a circular import at module load.

Source code in src/radtrax/slab/config.py

def make_srf(
    self, wavelengths_nm: Float[ArrayLike, " Ghr"]
) -> SpectralResponseFunction:
    """Construct a :class:`SpectralResponseFunction` over ``wavelengths_nm``.

    Imported locally to avoid a circular import at module load.
    """
    from radtrax.instrument.srf import SpectralResponseFunction

    return SpectralResponseFunction(
        wavelengths_hr_nm=np.asarray(wavelengths_nm, dtype=float),
        band_centers_nm=np.asarray(self.band_centers_nm, dtype=float),
        band_widths_nm=np.asarray(self.band_widths_nm, dtype=float),
        band_names=tuple(self.band_names),
        srf_type=self.srf_type,
    )

number_density_cm3

number_density_cm3(p_atm: float, T_K: float) -> float

Ideal-gas total number density [molecules / cm³].

Source code in src/radtrax/slab/config.py

def number_density_cm3(p_atm: float, T_K: float) -> float:
    """Ideal-gas total number density [molecules / cm³]."""
    if not (T_K > 0.0 and p_atm > 0.0):
        raise ParameterError(
            "number_density_cm3: `T_K` and `p_atm` must be > 0 "
            f"(got T_K={T_K!r}, p_atm={p_atm!r})"
        )
    return float(number_density(p_atm, T_K))

forward

Beer-Lambert forward models for atmospheric retrieval.

Three flavours, each returning radiance + Jacobian + transmittance:

• :func:forward_nonlinear — exact Beer-Lambert L = L₀ · exp(-τ).
• :func:forward_maclaurin — polynomial in VMR via Maclaurin of exp(-τ) around VMR = 0. Order 1 is the classical "linear" retrieval.
• :func:forward_taylor — Taylor-linearised around a background state VMR_bg (the state used in 3D-/4D-Var inner loops).

Plus normalised variants that divide by the background radiance:

• :func:forward_nonlinear_normalized — L_norm = exp(-Δτ); cancels surface reflectance, solar irradiance, and common aerosol slope.
• :func:forward_maclaurin_normalized / :func:forward_taylor_normalized.

All functions take an xarray.Dataset LUT carrying an absorption_cross_section variable with dims (wavenumber, temperature, pressure) — the output of :func:radtrax.spectroscopy.hitran.build_lut_dataset.

Ported and adapted from jej_vc_snippets/methane_retrieval/lut_model_beers.py — the Jacobian is now always returned alongside the radiance so the same functions can drive variational retrievals (H = dL/dVMR) and matched-filter target-spectrum construction (t = H · ΔVMR).

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

ForwardResult dataclass

ForwardResult(radiance: ndarray, jacobian: ndarray, transmittance: ndarray)

Bundle of (radiance, Jacobian, transmittance) returned by the forward models.

Using a frozen dataclass rather than a plain tuple gives the caller named access — result.jacobian — which matters because target spectra for the matched filter (:mod:radtrax.retrieval.matched_filter) consume only the Jacobian or only the transmittance.

Attributes

radiance : np.ndarray Simulated radiance at each wavenumber, shape (n_nu,). jacobian : np.ndarray dL/dVMR evaluated at the supplied state, shape (n_nu,). transmittance : np.ndarray Atmospheric transmittance exp(-τ) (non-normalised variants) or exp(-Δτ) (normalised variants), shape (n_nu,).

forward_nonlinear

forward_nonlinear(ds: Dataset, nu_obs: Float[ArrayLike, ' G'], *, T_K: float, p_atm: float, vmr: float, path_length_cm: float, amf: float, surface_reflectance: float = 1.0, solar_irradiance: float = 1.0, var: str = 'absorption_cross_section') -> ForwardResult

Exact Beer-Lambert forward model.

Radiance: L(ν) = (F₀ R / π) · exp(-τ(ν, VMR)) Optical τ: τ = σ · N_total · VMR · L · AMF Jacobian: dL/dVMR = -L · (dτ/dVMR), dτ/dVMR = σ · N_total · L · AMF

Defaults surface_reflectance = solar_irradiance = 1 collapse the prefactor to 1/π, which is convenient for toy tests where the absolute radiance scale does not matter.

Source code in src/radtrax/slab/forward.py

def forward_nonlinear(
    ds: xr.Dataset,
    nu_obs: Float[ArrayLike, " G"],
    *,
    T_K: float,
    p_atm: float,
    vmr: float,
    path_length_cm: float,
    amf: float,
    surface_reflectance: float = 1.0,
    solar_irradiance: float = 1.0,
    var: str = "absorption_cross_section",
) -> ForwardResult:
    """Exact Beer-Lambert forward model.

    Radiance:    ``L(ν) = (F₀ R / π) · exp(-τ(ν, VMR))``
    Optical τ:   ``τ = σ · N_total · VMR · L · AMF``
    Jacobian:    ``dL/dVMR = -L · (dτ/dVMR)``, ``dτ/dVMR = σ · N_total · L · AMF``

    Defaults ``surface_reflectance = solar_irradiance = 1`` collapse the
    prefactor to ``1/π``, which is convenient for toy tests where the
    absolute radiance scale does not matter.
    """
    sigma = _interp_sigma(ds, nu_obs, T_K, p_atm, var)
    N_total = number_density_cm3(p_atm, T_K)
    tau = sigma * N_total * vmr * path_length_cm * amf
    transmittance = np.exp(-tau)
    L0 = solar_irradiance * surface_reflectance / np.pi
    radiance = L0 * transmittance
    dtau_dvmr = sigma * N_total * path_length_cm * amf
    jacobian = -radiance * dtau_dvmr
    return ForwardResult(
        radiance=radiance, jacobian=jacobian, transmittance=transmittance
    )

forward_nonlinear_normalized

forward_nonlinear_normalized(ds: Dataset, nu_obs: Float[ArrayLike, ' G'], *, T_K: float, p_atm: float, vmr_background: float, delta_vmr: float, path_length_cm: float, amf: float, var: str = 'absorption_cross_section') -> ForwardResult

Exact Beer-Lambert normalised by the background radiance.

L_norm = L(VMR_bg + ΔVMR) / L(VMR_bg) = exp(-Δτ), which cancels any multiplicative surface/solar/aerosol factors that appear equally in the plume and background pixels.

Returns

ForwardResult radiance is the normalised transmittance exp(-Δτ); transmittance is the same quantity (kept for API symmetry with :func:forward_nonlinear); jacobian is d(L_norm)/d(ΔVMR) = -exp(-Δτ) · (σ · N · L · AMF).

Source code in src/radtrax/slab/forward.py

def forward_nonlinear_normalized(
    ds: xr.Dataset,
    nu_obs: Float[ArrayLike, " G"],
    *,
    T_K: float,
    p_atm: float,
    vmr_background: float,
    delta_vmr: float,
    path_length_cm: float,
    amf: float,
    var: str = "absorption_cross_section",
) -> ForwardResult:
    """Exact Beer-Lambert *normalised* by the background radiance.

    ``L_norm = L(VMR_bg + ΔVMR) / L(VMR_bg) = exp(-Δτ)``, which cancels any
    multiplicative surface/solar/aerosol factors that appear equally in the
    plume and background pixels.

    Returns
    -------
    ForwardResult
        ``radiance`` is the normalised transmittance ``exp(-Δτ)``;
        ``transmittance`` is the same quantity (kept for API symmetry with
        :func:`forward_nonlinear`); ``jacobian`` is
        ``d(L_norm)/d(ΔVMR) = -exp(-Δτ) · (σ · N · L · AMF)``.
    """
    sigma = _interp_sigma(ds, nu_obs, T_K, p_atm, var)
    N_total = number_density_cm3(p_atm, T_K)
    dtau_d_dvmr = sigma * N_total * path_length_cm * amf
    delta_tau = dtau_d_dvmr * delta_vmr
    L_norm = np.exp(-delta_tau)
    jacobian = -L_norm * dtau_d_dvmr
    return ForwardResult(radiance=L_norm, jacobian=jacobian, transmittance=L_norm)

forward_maclaurin

forward_maclaurin(ds: Dataset, nu_obs: Float[ArrayLike, ' G'], *, T_K: float, p_atm: float, vmr: float, path_length_cm: float, amf: float, surface_reflectance: float = 1.0, solar_irradiance: float = 1.0, order: int = 1, var: str = 'absorption_cross_section') -> ForwardResult

Maclaurin-series forward model: expand exp(-τ(VMR)) around VMR = 0.

With a = σ · N · L · AMF so τ = a · VMR:

• order 1: T ≈ 1 − a·VMR (linear in VMR).
• order 2: T ≈ 1 − a·VMR + ½ (a·VMR)².
• order 3: T ≈ 1 − a·VMR + ½ (a·VMR)² − ⅙ (a·VMR)³.

Accurate when the total optical depth a·VMR is ≪ 1 — the regime of classical linear retrievals for thin absorbers.

Source code in src/radtrax/slab/forward.py

def forward_maclaurin(
    ds: xr.Dataset,
    nu_obs: Float[ArrayLike, " G"],
    *,
    T_K: float,
    p_atm: float,
    vmr: float,
    path_length_cm: float,
    amf: float,
    surface_reflectance: float = 1.0,
    solar_irradiance: float = 1.0,
    order: int = 1,
    var: str = "absorption_cross_section",
) -> ForwardResult:
    """Maclaurin-series forward model: expand ``exp(-τ(VMR))`` around VMR = 0.

    With ``a = σ · N · L · AMF`` so ``τ = a · VMR``:

    - order 1: ``T ≈ 1 − a·VMR`` (linear in VMR).
    - order 2: ``T ≈ 1 − a·VMR + ½ (a·VMR)²``.
    - order 3: ``T ≈ 1 − a·VMR + ½ (a·VMR)² − ⅙ (a·VMR)³``.

    Accurate when the total optical depth ``a·VMR`` is ≪ 1 — the regime of
    classical linear retrievals for thin absorbers.
    """
    if order not in (1, 2, 3):
        raise ParameterError(
            f"forward_maclaurin: `order` must be 1, 2 or 3 (got {order!r})"
        )
    sigma = _interp_sigma(ds, nu_obs, T_K, p_atm, var)
    N_total = number_density_cm3(p_atm, T_K)
    a = sigma * N_total * path_length_cm * amf
    a_vmr = a * vmr

    if order == 1:
        transmittance = 1.0 - a_vmr
        dtrans_dvmr = -a
    elif order == 2:
        transmittance = 1.0 - a_vmr + 0.5 * a_vmr**2
        dtrans_dvmr = -a + (a**2) * vmr
    else:  # order == 3
        transmittance = 1.0 - a_vmr + 0.5 * a_vmr**2 - (1.0 / 6.0) * a_vmr**3
        dtrans_dvmr = -a + (a**2) * vmr - 0.5 * (a**3) * vmr**2

    L0 = solar_irradiance * surface_reflectance / np.pi
    radiance = L0 * transmittance
    jacobian = L0 * dtrans_dvmr
    return ForwardResult(
        radiance=radiance, jacobian=jacobian, transmittance=transmittance
    )

forward_maclaurin_normalized

forward_maclaurin_normalized(ds: Dataset, nu_obs: Float[ArrayLike, ' G'], *, T_K: float, p_atm: float, delta_vmr: float, path_length_cm: float, amf: float, order: int = 1, var: str = 'absorption_cross_section') -> ForwardResult

Maclaurin expansion of exp(-Δτ) around ΔVMR = 0.

At order 1 this reduces to L_norm ≈ 1 − σ · N · ΔVMR · L · AMF — the canonical linearised retrieval signal.

Source code in src/radtrax/slab/forward.py

def forward_maclaurin_normalized(
    ds: xr.Dataset,
    nu_obs: Float[ArrayLike, " G"],
    *,
    T_K: float,
    p_atm: float,
    delta_vmr: float,
    path_length_cm: float,
    amf: float,
    order: int = 1,
    var: str = "absorption_cross_section",
) -> ForwardResult:
    """Maclaurin expansion of ``exp(-Δτ)`` around ``ΔVMR = 0``.

    At order 1 this reduces to ``L_norm ≈ 1 − σ · N · ΔVMR · L · AMF`` — the
    canonical linearised retrieval signal.
    """
    if order not in (1, 2, 3):
        raise ParameterError("forward_maclaurin_normalized: `order` must be 1, 2 or 3.")
    sigma = _interp_sigma(ds, nu_obs, T_K, p_atm, var)
    N_total = number_density_cm3(p_atm, T_K)
    a = sigma * N_total * path_length_cm * amf
    a_dvmr = a * delta_vmr
    if order == 1:
        L_norm = 1.0 - a_dvmr
        jac = -a
    elif order == 2:
        L_norm = 1.0 - a_dvmr + 0.5 * a_dvmr**2
        jac = -a + (a**2) * delta_vmr
    else:
        L_norm = 1.0 - a_dvmr + 0.5 * a_dvmr**2 - (1.0 / 6.0) * a_dvmr**3
        jac = -a + (a**2) * delta_vmr - 0.5 * (a**3) * delta_vmr**2
    return ForwardResult(radiance=L_norm, jacobian=jac, transmittance=L_norm)

forward_taylor

forward_taylor(ds: Dataset, nu_obs: Float[ArrayLike, ' G'], *, T_K: float, p_atm: float, vmr: float, vmr_background: float, path_length_cm: float, amf: float, surface_reflectance: float = 1.0, solar_irradiance: float = 1.0, var: str = 'absorption_cross_section') -> ForwardResult

First-order Taylor expansion around vmr_background.

L(VMR) ≈ L_b + H · (VMR − VMR_b) with H = dL/dVMR |_{VMR_b}. This is the linearisation used in the inner loops of 3D-/4D-Var.

Source code in src/radtrax/slab/forward.py

def forward_taylor(
    ds: xr.Dataset,
    nu_obs: Float[ArrayLike, " G"],
    *,
    T_K: float,
    p_atm: float,
    vmr: float,
    vmr_background: float,
    path_length_cm: float,
    amf: float,
    surface_reflectance: float = 1.0,
    solar_irradiance: float = 1.0,
    var: str = "absorption_cross_section",
) -> ForwardResult:
    """First-order Taylor expansion around ``vmr_background``.

    ``L(VMR) ≈ L_b + H · (VMR − VMR_b)`` with ``H = dL/dVMR |_{VMR_b}``. This
    is the linearisation used in the inner loops of 3D-/4D-Var.
    """
    sigma = _interp_sigma(ds, nu_obs, T_K, p_atm, var)
    N_total = number_density_cm3(p_atm, T_K)
    dtau_dvmr = sigma * N_total * path_length_cm * amf
    tau_b = dtau_dvmr * vmr_background
    L0 = solar_irradiance * surface_reflectance / np.pi
    L_b = L0 * np.exp(-tau_b)
    H = -L_b * dtau_dvmr
    radiance = L_b + H * (vmr - vmr_background)
    transmittance = np.exp(-tau_b) + (-dtau_dvmr * np.exp(-tau_b)) * (
        vmr - vmr_background
    )
    return ForwardResult(radiance=radiance, jacobian=H, transmittance=transmittance)

forward_taylor_normalized

forward_taylor_normalized(ds: Dataset, nu_obs: Float[ArrayLike, ' G'], *, T_K: float, p_atm: float, delta_vmr: float, path_length_cm: float, amf: float, var: str = 'absorption_cross_section') -> ForwardResult

Taylor expansion of L_norm around ΔVMR = 0.

Since L_norm(0) = 1 and dL_norm/dΔVMR(0) = -σ · N · L · AMF, the first-order Taylor expansion coincides with the order-1 Maclaurin expansion. Kept as a named function for API symmetry with :func:forward_taylor.

Source code in src/radtrax/slab/forward.py

def forward_taylor_normalized(
    ds: xr.Dataset,
    nu_obs: Float[ArrayLike, " G"],
    *,
    T_K: float,
    p_atm: float,
    delta_vmr: float,
    path_length_cm: float,
    amf: float,
    var: str = "absorption_cross_section",
) -> ForwardResult:
    """Taylor expansion of ``L_norm`` around ``ΔVMR = 0``.

    Since ``L_norm(0) = 1`` and ``dL_norm/dΔVMR(0) = -σ · N · L · AMF``,
    the first-order Taylor expansion coincides with the order-1 Maclaurin
    expansion. Kept as a named function for API symmetry with
    :func:`forward_taylor`.
    """
    return forward_maclaurin_normalized(
        ds,
        nu_obs,
        T_K=T_K,
        p_atm=p_atm,
        delta_vmr=delta_vmr,
        path_length_cm=path_length_cm,
        amf=amf,
        order=1,
        var=var,
    )

nb_lut

Normalised-brightness LUT for fast per-pixel plume injection.

The band-integrated normalised brightness for a pixel with methane column enhancement ΔX over a background column X_bg is

nB_b(ΔX) = T_b(X_bg + ΔX) / T_b(X_bg),
T_b(X)   = ∫ f_b(λ) · exp(-AMF · σ(λ) · X) dλ  /  ∫ f_b(λ) dλ

where f_b is the instrument SRF for band b and σ is the target gas cross-section in cm²/molecule. This integral is the natural scalar "plume throughput" for band b once surface reflectance and solar irradiance have been divided out (differential Beer-Lambert).

The background matters: it has already darkened the strongest line cores, so a plume over it absorbs less than in clean air. plumax built the table with X_bg = 0, which overstates Sentinel-2 band 12's plume optical depth by 1.1 % (radtrax · 2, the GeoStack design note); pass the scene's background column as background_column (radtrax#36). For methane at 1.9 ppm over a 1013 hPa column it is 0.678 mol/m² (0.0109 kg/m²). From a profile, AtmosphericState.column(gas) is in molecules/cm²: divide by AVOGADRO_PER_MOL and multiply by CM2_PER_M2 for mol/m².

Evaluating the integral per pixel is expensive, so we follow the Eucalyptus radtran() trick: pre-tabulate nB over a 1-D grid of ΔX values, then look up per-pixel via piecewise-linear interpolation (np.interp). The grid is quadratic, ΔX_k = ΔX_max (k / (n − 1))²: the interpolation error of a small plume is first order in the step near 0, where the band optical depth curves most, so plumax's uniform grid of 4001 points over [0, 200] mol/m² missed the exact band integral by up to 1.4 % at small ΔX on radtrax · 3's synthetic table, and the quadratic one by 4e-6 (gate V5, radtrax#36). Plume injection is a single np.interp call per band.

Public surface

• :func:build_nb_lut — build the (n_bands, n_delta) LUT.
• :func:lookup_nb — look up nB per band for a 2-D plume map.
• :func:inject_plume — apply nB to a clean multispectral scene.

Ported from plumax 0.1.2 at 28f257e (src/plumax/radtran/nb_lut.py); the background column is radtrax's (radtrax#36).

NBLookup dataclass

NBLookup(delta_column: ndarray, nB: ndarray, band_names: tuple[str, ...], background_column: float = 0.0)

Pre-computed normalised-brightness table, indexed by band and ΔX.

Attributes

delta_column : np.ndarray Column-enhancement grid [mol/m²], shape (n_delta,). Strictly increasing, starts at 0 (quadratic for :func:build_nb_lut). nB : np.ndarray Normalised brightness per band at each grid point, shape (n_bands, n_delta). nB[:, 0] == 1 always. band_names : tuple of str Band labels, length n_bands. background_column : float Background column X_bg the table is relative to [mol/m²]; 0 for plumax's clean-air table.

to_dataset

to_dataset() -> Dataset

Wrap as an :class:xarray.Dataset for easy serialisation.

Source code in src/radtrax/slab/nb_lut.py

def to_dataset(self) -> xr.Dataset:
    """Wrap as an :class:`xarray.Dataset` for easy serialisation."""
    return xr.Dataset(
        data_vars={
            "nB": (["band", "delta_column"], self.nB),
        },
        coords={
            "band": list(self.band_names),
            "delta_column": (
                "delta_column",
                self.delta_column,
                {"units": "mol/m^2"},
            ),
        },
        attrs={
            "long_name": "Band-integrated normalised brightness",
            "description": (
                "nB_b(ΔX) = T_b(X_bg + ΔX) / T_b(X_bg), "
                "T_b(X) = ∫ f_b(λ) exp(-AMF σ(λ) X) dλ / ∫ f_b(λ) dλ"
            ),
            "background_column_mol_m2": self.background_column,
        },
    )

build_nb_lut

build_nb_lut(ds: Dataset, srf: SpectralResponseFunction, *, T_K: float, p_atm: float, amf: float, max_delta_column: float = 200.0, n_grid: int = 4001, var: str = 'absorption_cross_section', background_column: float = 0.0) -> NBLookup

Pre-tabulate band-integrated nB over a column-enhancement grid.

Parameters

ds : xr.Dataset HAPI absorption-cross-section LUT, output of :func:radtrax.spectroscopy.hitran.build_lut_dataset. srf : SpectralResponseFunction Instrument SRF mapping the LUT wavelength grid to bands. T_K, p_atm : float Slab atmospheric temperature [K] and pressure [atm], at which the cross-section σ(λ; T, P) is interpolated. ΔX is a column already, so the slab path length does not enter the integrand. amf : float Two-way air-mass factor multiplying the slant column. max_delta_column : float Upper bound of the ΔX grid [mol/m²]. Default 200. n_grid : int Number of grid points on the quadratic grid (default 4001: steps from 1.2e-5 mol/m² at 0 to 0.1 mol/m² at 200). var : str Cross-section variable name (default 'absorption_cross_section', matching :func:radtrax.spectroscopy.hitran.build_lut_dataset). background_column : float Background column X_bg [mol/m², not molecules/cm²] the plume sits on; nB is relative to it. The default, 0, is plumax's clean-air table, which overstates a strong band's plume optical depth (module docstring).

Returns

NBLookup nB[b, k] is the band-integrated normalised brightness at delta_column[k] for band b.

Source code in src/radtrax/slab/nb_lut.py

def build_nb_lut(
    ds: xr.Dataset,
    srf: SpectralResponseFunction,
    *,
    T_K: float,
    p_atm: float,
    amf: float,
    max_delta_column: float = 200.0,
    n_grid: int = 4001,
    var: str = "absorption_cross_section",
    background_column: float = 0.0,
) -> NBLookup:
    """Pre-tabulate band-integrated nB over a column-enhancement grid.

    Parameters
    ----------
    ds : xr.Dataset
        HAPI absorption-cross-section LUT, output of
        :func:`radtrax.spectroscopy.hitran.build_lut_dataset`.
    srf : SpectralResponseFunction
        Instrument SRF mapping the LUT wavelength grid to bands.
    T_K, p_atm : float
        Slab atmospheric temperature [K] and pressure [atm], at which the
        cross-section ``σ(λ; T, P)`` is interpolated. ΔX is a column
        already, so the slab path length does not enter the integrand.
    amf : float
        Two-way air-mass factor multiplying the slant column.
    max_delta_column : float
        Upper bound of the ΔX grid [mol/m²]. Default 200.
    n_grid : int
        Number of grid points on the quadratic grid (default 4001: steps
        from 1.2e-5 mol/m² at 0 to 0.1 mol/m² at 200).
    var : str
        Cross-section variable name (default ``'absorption_cross_section'``,
        matching :func:`radtrax.spectroscopy.hitran.build_lut_dataset`).
    background_column : float
        Background column ``X_bg`` [mol/m², not molecules/cm²] the plume
        sits on; ``nB`` is relative to it. The default, 0, is plumax's
        clean-air table, which overstates a strong band's plume optical
        depth (module docstring).

    Returns
    -------
    NBLookup
        ``nB[b, k]`` is the band-integrated normalised brightness at
        ``delta_column[k]`` for band ``b``.
    """
    if var not in ds:
        raise MissingKeyError(f"build_nb_lut: {var!r} not in dataset.")
    if n_grid < 2:
        raise ParameterError(f"build_nb_lut: `n_grid` must be ≥ 2 (got {n_grid!r})")
    if not (0.0 <= background_column <= _MAX_COLUMN_MOL_PER_M2):
        raise ParameterError(
            f"build_nb_lut: `background_column` must be in [0, "
            f"{_MAX_COLUMN_MOL_PER_M2:g}] mol/m² (got {background_column!r}); a "
            "column in molecules/cm² is 6.02e19 times larger than in mol/m²."
        )
    if not (max_delta_column > 0.0):
        raise ParameterError(
            f"build_nb_lut: `max_delta_column` must be > 0 (got {max_delta_column!r})"
        )

    # Cross-section on the LUT wavenumber grid, at (T, P).
    sigma_nu = ds[var].interp(temperature=T_K, pressure=p_atm, method="linear")

    # Convert to wavelength (nm) to match the SRF. Drop the unused coord to
    # avoid `interp` warning about non-dim-coord changes.
    sigma_lambda = sigma_nu.assign_coords(
        lam_nm=("wavenumber", np.asarray(1e7 / ds["wavenumber"].values, dtype=float))
    )
    # Interpolate σ onto the SRF's HR wavelength grid (decreasing in lam is fine
    # for np.interp after flipping).
    lam_src = sigma_lambda["lam_nm"].values
    sort = np.argsort(lam_src)
    # Flag bands whose SRF response reaches past the LUT wavelength coverage;
    # np.interp zero-fills σ there, producing a silently-attenuated nB (≡ 1).
    check_srf_lut_coverage(srf, lam_src, context="build_nb_lut")
    sigma_hr = np.interp(
        srf.wavelengths_hr_nm,
        lam_src[sort],
        sigma_lambda.values[sort],
        left=0.0,
        right=0.0,
    )

    # Convert ΔX [mol/m²] to the per-wavelength optical depth prefactor.
    #   τ_λ(ΔX) = σ_λ [cm²/molec] · ΔX [mol/m²] · N_A [molec/mol] · 1e-4 [m²/cm²] · AMF
    # No path_length_cm appears here — ΔX already integrates the column.
    optical_coef = sigma_hr * gas.AVOGADRO_PER_MOL / gas.CM2_PER_M2 * amf  # (n_lambda,)

    # ΔX_k = ΔX_max (k / (n − 1))²: dense near 0, where τ_b(ΔX) curves most.
    delta_column = float(max_delta_column) * np.linspace(0.0, 1.0, int(n_grid)) ** 2
    # Band transmittance over the background plus the plume:
    #   T_b(X_bg + ΔX) = Σ_λ f_bλ exp(−c_λ (X_bg + ΔX))       (n_bands, n_delta)
    # Chunked over ΔX, so memory stays at _CHUNK × n_lambda whatever n_grid.
    columns = float(background_column) + delta_column  # (n_delta,)
    band_transmittance = np.concatenate(
        [
            srf.apply(np.exp(-optical_coef[None, :] * chunk[:, None]))
            for chunk in np.array_split(columns, -(-columns.size // _CHUNK))
        ]
    ).T  # (n_bands, n_delta)
    # nB_b(ΔX) = T_b(X_bg + ΔX) / T_b(X_bg): exactly 1 at ΔX = 0.
    nB = band_transmittance / band_transmittance[:, :1]
    return NBLookup(
        delta_column=delta_column,
        nB=nB,
        band_names=srf.band_names,
        background_column=float(background_column),
    )

lookup_nb

lookup_nb(plume_column: Float[ArrayLike, '*shape'], nb_lookup: NBLookup) -> Float[ndarray, ' B *shape']

Look up per-band nB for every pixel of a 2-D column map.

Values outside the LUT range are clipped. Clipping at the top end is physically conservative (saturates the absorption); negative inputs are clipped to zero since they are not physical.

Parameters

plume_column : np.ndarray Column enhancement [mol/m²], shape (ny, nx) or any broadcastable shape. NaNs are preserved as NaNs in the output. nb_lookup : NBLookup Pre-built LUT.

Returns

nb_map : np.ndarray Shape (n_bands, *plume_column.shape). Dimensionless multiplier in [0, 1] (1 at ΔX=0).

Source code in src/radtrax/slab/nb_lut.py

def lookup_nb(
    plume_column: Float[ArrayLike, "*shape"],
    nb_lookup: NBLookup,
) -> Float[np.ndarray, " B *shape"]:
    """Look up per-band nB for every pixel of a 2-D column map.

    Values outside the LUT range are clipped. Clipping at the top end is
    physically conservative (saturates the absorption); negative inputs are
    clipped to zero since they are not physical.

    Parameters
    ----------
    plume_column : np.ndarray
        Column enhancement [mol/m²], shape ``(ny, nx)`` or any broadcastable
        shape. NaNs are preserved as NaNs in the output.
    nb_lookup : NBLookup
        Pre-built LUT.

    Returns
    -------
    nb_map : np.ndarray
        Shape ``(n_bands, *plume_column.shape)``. Dimensionless multiplier
        in ``[0, 1]`` (1 at ΔX=0).
    """
    x = np.asarray(plume_column, dtype=float)
    grid = nb_lookup.delta_column
    nan_mask = np.isnan(x)
    x_clipped = np.clip(np.where(nan_mask, 0.0, x), grid[0], grid[-1])

    # Piecewise-linear interp per band (avoids nearest-neighbour stair-step
    # artefacts that the Eucalyptus `lookup_nB` shows on smooth plumes).
    out_shape = (nb_lookup.n_bands, *x.shape)
    out = np.empty(out_shape, dtype=float)
    for b in range(nb_lookup.n_bands):
        out[b] = np.interp(x_clipped, grid, nb_lookup.nB[b])
    if nan_mask.any():
        out[:, nan_mask] = np.nan
    return out

inject_plume

inject_plume(scene: Float[ArrayLike, '*scene'], plume_column: Float[ArrayLike, 'H W'], nb_lookup: NBLookup, *, band_axis: int = 0) -> Float[ndarray, '*scene']

Multiplicative plume injection: scene_plume = scene · nB(ΔX).

Parameters

scene : np.ndarray Clean multispectral scene, shape (n_bands, ny, nx) by default (band_axis=0). plume_column : np.ndarray ΔX column map [mol/m²], shape (ny, nx). nb_lookup : NBLookup Pre-built LUT. nb_lookup.n_bands must match scene.shape[band_axis]. band_axis : int Axis of scene corresponding to bands. Default 0.

Returns

scene_plume : np.ndarray Same shape as scene.

Source code in src/radtrax/slab/nb_lut.py

def inject_plume(
    scene: Float[ArrayLike, "*scene"],
    plume_column: Float[ArrayLike, "H W"],
    nb_lookup: NBLookup,
    *,
    band_axis: int = 0,
) -> Float[np.ndarray, "*scene"]:
    """Multiplicative plume injection: ``scene_plume = scene · nB(ΔX)``.

    Parameters
    ----------
    scene : np.ndarray
        Clean multispectral scene, shape ``(n_bands, ny, nx)`` by default
        (``band_axis=0``).
    plume_column : np.ndarray
        ΔX column map [mol/m²], shape ``(ny, nx)``.
    nb_lookup : NBLookup
        Pre-built LUT. ``nb_lookup.n_bands`` must match ``scene.shape[band_axis]``.
    band_axis : int
        Axis of ``scene`` corresponding to bands. Default 0.

    Returns
    -------
    scene_plume : np.ndarray
        Same shape as ``scene``.
    """
    scene_arr = np.asarray(scene)
    if scene_arr.shape[band_axis] != nb_lookup.n_bands:
        raise ShapeError(
            f"inject_plume: scene has {scene_arr.shape[band_axis]} bands along "
            f"axis {band_axis}, but LUT has {nb_lookup.n_bands}."
        )
    nb_map = lookup_nb(plume_column, nb_lookup)  # (n_bands, ny, nx)
    # Move `scene`'s band axis to the front for broadcasting, then move back.
    scene_bf = np.moveaxis(scene_arr, band_axis, 0)
    if scene_bf.shape[1:] != nb_map.shape[1:]:
        raise ShapeError(
            f"inject_plume: spatial shape mismatch — scene {scene_bf.shape[1:]} "
            f"vs. plume {nb_map.shape[1:]}"
        )
    return np.moveaxis(scene_bf * nb_map, 0, band_axis)