Surface¶

Surfaces give the bidirectional reflectance factor R(μ₀, μ) a solver reflects, behind the Surface protocol. v0 has Lambertian, R = A; PolynomialAlbedo makes A a Legendre series in wavenumber; RossLi is the MODIS kernel BRDF (RossThick and LiSparse-Reciprocal). RPV and Cox–Munk come later in the roadmap.

surface ¶

Surface reflectance behind the Surface protocol: Lambertian in v0 (radtrax#35), a spectral (Legendre-series) albedo (radtrax#106) and the Ross–Li kernel BRDF (radtrax#105).

`lambertian`¶

The Surface protocol and the Lambertian surface.

A surface turns the viewing geometry into a bidirectional reflectance factor R (dimensionless): the reflected radiance relative to that of a white Lambertian surface under the same illumination, so the radiance leaving the surface is L = E₀ μ₀ R / π. A Lambertian surface reflects the same in every direction, R = A, its albedo (radtrax · 3). Kernel BRDFs (Ross–Li, RPV) and Cox–Munk arrive with v1 of the roadmap; they use the relative azimuth in :class:~radtrax.atmosphere.ViewingGeometry and, for a spectral albedo model, the wavenumber grid.

Shapes: G wavenumbers; *G means a scalar or (G,).

Surface ¶

Bases: Module

Protocol: the surface's bidirectional reflectance factor.

An abstract Equinox module, so filter_jit compiles each concrete surface and a solver can take any of them.

reflectance `abstractmethod` ¶

reflectance(geometry: ViewingGeometry, nu_cm: Float[Array, ' G']) -> Float[Array, '*G']

Bidirectional reflectance factor R (dimensionless).

Parameters:

Name	Type	Description	Default
`geometry`	`ViewingGeometry`	Sun and sensor directions.	required
`nu_cm`	`Float[Array, ' G']`	Wavenumber grid [cm⁻¹], `(G,)`.	required

Returns:

Type	Description
`Float[Array, '*G']`	`R`, a scalar or one value per wavenumber, `(G,)`.

Source code in src/radtrax/surface/lambertian.py

@abc.abstractmethod
def reflectance(
    self, geometry: ViewingGeometry, nu_cm: Float[Array, " G"]
) -> Float[Array, "*G"]:
    """Bidirectional reflectance factor ``R`` (dimensionless).

    Args:
        geometry: Sun and sensor directions.
        nu_cm: Wavenumber grid [cm⁻¹], ``(G,)``.

    Returns:
        ``R``, a scalar or one value per wavenumber, ``(G,)``.
    """

Lambertian ¶

Lambertian(albedo: Real[ArrayLike, '*G'])

Bases: Surface

A Lambertian surface: R = A in every direction.

Attributes:

Name	Type	Description
`albedo`	`Float[Array, '*G']`	Surface albedo `A` (dimensionless), a scalar or one value per wavenumber of the grid it is used on, `(G,)`.

Examples:

>>> import jax.numpy as jnp
>>> from radtrax.atmosphere import ViewingGeometry
>>> surface = Lambertian(0.25)
>>> geometry = ViewingGeometry.from_angles(35.0, 0.0)
>>> float(surface.reflectance(geometry, jnp.array([4300.0])))
0.25

Source code in src/radtrax/surface/lambertian.py

def __init__(self, albedo: Real[ArrayLike, "*G"]) -> None:
    self.albedo = jnp.asarray(albedo, dtype=float)

reflectance ¶

reflectance(geometry: ViewingGeometry, nu_cm: Float[Array, ' G']) -> Float[Array, '*G']

R = A, whatever the geometry and wavenumber.

Source code in src/radtrax/surface/lambertian.py

def reflectance(
    self, geometry: ViewingGeometry, nu_cm: Float[Array, " G"]
) -> Float[Array, "*G"]:
    """``R = A``, whatever the geometry and wavenumber."""
    del geometry, nu_cm
    return self.albedo

`polynomial`¶

A spectral albedo: a low-order Legendre series in wavenumber.

Real SWIR albedo has a slope across a fitting window, and a slope aliases into a methane band (plumax#98: a 10 % slope biases the column). Optimal estimation retrieves it as nuisance state, which needs a few well-conditioned coefficients rather than one value per wavenumber:

\[ A(\nu) = \sum_{k=0}^{K-1} a_k\, P_k(x),\qquad x = \frac{2(\nu - \nu_{\min})}{\nu_{\max} - \nu_{\min}} - 1, \]

with P_k the Legendre polynomials (geonnax.basis.legendre_polynomials) on the window [ν_min, ν_max], where they are orthogonal. The surface is Lambertian at each wavenumber: R(ν) = A(ν) in every direction.

Shapes: G wavenumbers; K coefficients (degree K − 1).

PolynomialAlbedo ¶

PolynomialAlbedo(coefficients: Real[ArrayLike, ' K'], nu_range_cm: Real[ArrayLike, ' 2'])

Bases: Surface

A Lambertian surface whose albedo is a Legendre series in wavenumber.

A(ν) = Σ_k a_k P_k(x) with x the wavenumber mapped onto [−1, 1] across nu_range_cm. With one coefficient it is Lambertian(a_0). Outside the window the series extrapolates; keep the model grid inside it. The window is checked to increase when it is concrete; a window traced under jit is not checked.

Attributes:

Name	Type	Description
`coefficients`	`Float[Array, ' K']`	Legendre coefficients `a_0 … a_{K−1}` (dimensionless), `(K,)`; `a_0` is the window's mean albedo.
`nu_range_cm`	`Float[Array, ' 2']`	The window `[ν_min, ν_max]` [cm⁻¹], `(2,)`, increasing.

Examples:

>>> import jax.numpy as jnp
>>> from radtrax.atmosphere import ViewingGeometry
>>> surface = PolynomialAlbedo(
...     jnp.array([0.25, 0.05]), nu_range_cm=jnp.array([4200.0, 4400.0])
... )
>>> geometry = ViewingGeometry.from_angles(35.0, 0.0)
>>> nu = jnp.array([4200.0, 4300.0, 4400.0])
>>> [round(v, 6) for v in surface.reflectance(geometry, nu).tolist()]
[0.2, 0.25, 0.3]

Source code in src/radtrax/surface/polynomial.py

def __init__(
    self, coefficients: Real[ArrayLike, " K"], nu_range_cm: Real[ArrayLike, " 2"]
) -> None:
    self.coefficients = jnp.asarray(coefficients, dtype=float)
    self.nu_range_cm = jnp.asarray(nu_range_cm, dtype=float)

albedo ¶

albedo(nu_cm: Float[Array, ' G']) -> Float[Array, ' G']

A(ν) on the grid (dimensionless), (G,).

Source code in src/radtrax/surface/polynomial.py

def albedo(self, nu_cm: Float[Array, " G"]) -> Float[Array, " G"]:
    """``A(ν)`` on the grid (dimensionless), ``(G,)``."""
    lo, hi = self.nu_range_cm[0], self.nu_range_cm[1]
    x = 2.0 * (jnp.asarray(nu_cm) - lo) / (hi - lo) - 1.0
    degree = jnp.shape(self.coefficients)[0] - 1
    return legendre_polynomials(x, degree) @ self.coefficients

reflectance ¶

reflectance(geometry: ViewingGeometry, nu_cm: Float[Array, ' G']) -> Float[Array, '*G']

R(ν) = A(ν), whatever the geometry.

Source code in src/radtrax/surface/polynomial.py

def reflectance(
    self, geometry: ViewingGeometry, nu_cm: Float[Array, " G"]
) -> Float[Array, "*G"]:
    """``R(ν) = A(ν)``, whatever the geometry."""
    del geometry
    return self.albedo(nu_cm)

`ross_li`¶

The Ross–Li kernel BRDF: RossThick and LiSparse-Reciprocal.

The MODIS kernel-driven model (Lucht et al. 2000) writes the bidirectional reflectance factor as a linear combination of an isotropic term and two kernels of the geometry:

\[ R(\theta_0, \theta, \phi) = f_{\mathrm{iso}} + f_{\mathrm{vol}} K_{\mathrm{vol}} + f_{\mathrm{geo}} K_{\mathrm{geo}}. \]

RossThick (volume scattering by a dense canopy), with the phase angle ξ, cos ξ = cos θ₀ cos θ + sin θ₀ sin θ cos φ: K_vol = ((π/2 − ξ) cos ξ + sin ξ) / (cos θ₀ + cos θ) − π/4.
LiSparse-Reciprocal (geometric-optical shadowing by sparse crowns) with MODIS's crown shape h/b = 2, b/r = 1 (so θ' = θ): K_geo = O − sec θ₀ − sec θ + ½ (1 + cos ξ) sec θ₀ sec θ, with the overlap O = (t − sin t cos t)(sec θ₀ + sec θ) / π and cos t = (h/b) √(D² + (tan θ₀ tan θ sin φ)²) / (sec θ₀ + sec θ), clipped to [−1, 1], D² = tan² θ₀ + tan² θ − 2 tan θ₀ tan θ cos φ.

φ is :class:~radtrax.atmosphere.ViewingGeometry's relative azimuth, sensor minus sun, so φ = 0 with θ₀ = θ is the hotspot (backscatter). Both kernels are 0 at nadir/nadir and symmetric in θ₀ ↔ θ.

First-order gradients are finite everywhere, nadir and the hotspot included (second derivatives are not: the Hessian at exact nadir is NaN):

(π/2 − ξ) cos ξ + sin ξ is c · arcsin c + √(1 − c²) in c = cos ξ, whose derivative is arcsin c; and t − sin t cos t is arccos u − u √(1 − u²) in u = cos t, whose derivative is −2 √(1 − u²). Both are written with that derivative (custom_jvp), where arccos differentiated naively gives 0 · ∞ at ±1.
√(D² + …) vanishes at nadir and at the hotspot, where the kernel has a kink (a distance); the double-where there gives the symmetric (sub)gradient 0, which is what a central difference measures.

Shapes: G wavenumbers; *G means a scalar or (G,).

RossLi ¶

RossLi(f_iso: Real[ArrayLike, '*G'], f_vol: Real[ArrayLike, '*G'], f_geo: Real[ArrayLike, '*G'])

Bases: Surface

The Ross–Li kernel BRDF: R = f_iso + f_vol K_vol + f_geo K_geo.

With f_vol = f_geo = 0 it is Lambertian(f_iso).

Attributes:

Name	Type	Description
`f_iso`	`Float[Array, '*G']`	Isotropic weight (dimensionless), a scalar or `(G,)`.
`f_vol`	`Float[Array, '*G']`	RossThick (volume) weight, a scalar or `(G,)`.
`f_geo`	`Float[Array, '*G']`	LiSparse-Reciprocal (geometric) weight, a scalar or `(G,)`.

Examples:

>>> import jax.numpy as jnp
>>> from radtrax.atmosphere import ViewingGeometry
>>> surface = RossLi(0.3, 0.1, 0.05)
>>> nadir = ViewingGeometry(1.0, 1.0)
>>> float(surface.reflectance(nadir, jnp.array([4300.0])))
0.3

Source code in src/radtrax/surface/ross_li.py

def __init__(
    self,
    f_iso: Real[ArrayLike, "*G"],
    f_vol: Real[ArrayLike, "*G"],
    f_geo: Real[ArrayLike, "*G"],
) -> None:
    self.f_iso = jnp.asarray(f_iso, dtype=float)
    self.f_vol = jnp.asarray(f_vol, dtype=float)
    self.f_geo = jnp.asarray(f_geo, dtype=float)

reflectance ¶

reflectance(geometry: ViewingGeometry, nu_cm: Float[Array, ' G']) -> Float[Array, '*G']

R = f_iso + f_vol K_vol + f_geo K_geo; the kernels do not depend on ν.

Source code in src/radtrax/surface/ross_li.py

def reflectance(
    self, geometry: ViewingGeometry, nu_cm: Float[Array, " G"]
) -> Float[Array, "*G"]:
    """``R = f_iso + f_vol K_vol + f_geo K_geo``; the kernels do not depend on ``ν``."""
    del nu_cm
    return (
        self.f_iso
        + self.f_vol * ross_thick(geometry)
        + self.f_geo * li_sparse_reciprocal(geometry)
    )

ross_thick ¶

ross_thick(geometry: ViewingGeometry) -> Float[Array, '']

The RossThick volume-scattering kernel K_vol (dimensionless).

Parameters:

Name	Type	Description	Default
`geometry`	`ViewingGeometry`	Sun and sensor directions.	required

Returns:

Type	Description
`Float[Array, '']`	`K_vol`, 0 at nadir/nadir.

Examples:

>>> from radtrax.atmosphere import ViewingGeometry
>>> round(float(ross_thick(ViewingGeometry(1.0, 1.0))), 12)
0.0

Source code in src/radtrax/surface/ross_li.py

def ross_thick(geometry: ViewingGeometry) -> Float[Array, ""]:
    """The RossThick volume-scattering kernel ``K_vol`` (dimensionless).

    Args:
        geometry: Sun and sensor directions.

    Returns:
        ``K_vol``, 0 at nadir/nadir.

    Examples:
        >>> from radtrax.atmosphere import ViewingGeometry
        >>> round(float(ross_thick(ViewingGeometry(1.0, 1.0))), 12)
        0.0
    """
    mu0, mu = geometry.mu0, geometry.mu
    cos_xi = _cos_phase(geometry)
    return _ross_numerator(cos_xi) / (mu0 + mu) - jnp.pi / 4

li_sparse_reciprocal ¶

li_sparse_reciprocal(geometry: ViewingGeometry) -> Float[Array, '']

The LiSparse-Reciprocal geometric kernel K_geo (dimensionless), MODIS's shape.

Parameters:

Name	Type	Description	Default
`geometry`	`ViewingGeometry`	Sun and sensor directions.	required

Returns:

Type	Description
`Float[Array, '']`	`K_geo`, 0 at nadir/nadir.

Examples:

>>> from radtrax.atmosphere import ViewingGeometry
>>> round(float(li_sparse_reciprocal(ViewingGeometry(1.0, 1.0))), 12)
0.0

Source code in src/radtrax/surface/ross_li.py

def li_sparse_reciprocal(geometry: ViewingGeometry) -> Float[Array, ""]:
    """The LiSparse-Reciprocal geometric kernel ``K_geo`` (dimensionless), MODIS's shape.

    Args:
        geometry: Sun and sensor directions.

    Returns:
        ``K_geo``, 0 at nadir/nadir.

    Examples:
        >>> from radtrax.atmosphere import ViewingGeometry
        >>> round(float(li_sparse_reciprocal(ViewingGeometry(1.0, 1.0))), 12)
        0.0
    """
    mu0, mu, phi = geometry.mu0, geometry.mu, geometry.relative_azimuth_rad
    sec0, sec1 = 1.0 / mu0, 1.0 / mu
    tan0, tan1 = _safe_sqrt(1.0 - mu0**2) / mu0, _safe_sqrt(1.0 - mu**2) / mu
    tan01 = tan0 * tan1
    # D² = tan²θ₀ + tan²θ − 2 tanθ₀ tanθ cos φ, rearranged with
    # 1 − cos φ = 2 sin²(φ/2) so it does not cancel near the hotspot.
    d_sq = (tan0 - tan1) ** 2 + 4.0 * tan01 * jnp.sin(phi / 2) ** 2
    root = _safe_sqrt(d_sq + (tan01 * jnp.sin(phi)) ** 2)
    cos_t = _H_OVER_B * root / (sec0 + sec1)
    overlap = _overlap_shape(cos_t) * (sec0 + sec1) / jnp.pi
    cos_xi = _cos_phase(geometry)
    return overlap - sec0 - sec1 + 0.5 * (1.0 + cos_xi) * sec0 * sec1

Surface¶

surface ¶

lambertian¶

Surface ¶

reflectance abstractmethod ¶

Lambertian ¶

reflectance ¶

polynomial¶

PolynomialAlbedo ¶

albedo ¶

reflectance ¶

ross_li¶

RossLi ¶

reflectance ¶

ross_thick ¶

li_sparse_reciprocal ¶

`lambertian`¶

reflectance `abstractmethod` ¶

`polynomial`¶

`ross_li`¶