radiomax: design details

The reference behind radiomax: core types, the arrays callers pass in, the algorithms with their pitfalls, JAX notes, tests, dependencies, and a survey of the code radiomax takes over, read on 5 Nov 2026.

Core types¶

Every type is an Equinox module of arrays, so it traces, batches and differentiates. Sensor-specific coefficients are leaves; the sensor’s name is static. Angles are in radians and lengths in metres, as everywhere below geophysix.

# Shapes: B bands; Ny × Nx pixels; K coefficients of an empirical formula

import equinox as eqx
from jaxtyping import Array, Float


# proposed (radiomax)
class ViewGeometry(eqx.Module):
    """Solar and viewing angles per pixel, in radians."""

    solar_zenith: Float[Array, "y x"]
    solar_azimuth: Float[Array, "y x"]
    view_zenith: Float[Array, "y x"]
    view_azimuth: Float[Array, "y x"]


# proposed (radiomax)
class BandCalibration(eqx.Module):
    """Counts to radiance per band: L = g · DN + o."""

    gain: Float[Array, "band"]
    offset: Float[Array, "band"]


# proposed (radiomax)
class PlanckLinear(eqx.Module):
    """GOES-style band inverse: T = (fk₂ / ln(fk₁ / L + 1) − b₁) / b₂."""

    fk1: Float[Array, "band"]
    fk2: Float[Array, "band"]
    bc1: Float[Array, "band"]
    bc2: Float[Array, "band"]


# proposed (radiomax)
class PlanckQuadratic(eqx.Module):
    """Himawari-style: Planck at the central wavelength, then T = c₀ + c₁ Tₑ + c₂ Tₑ²."""

    wavelength: Float[Array, "band"]
    c0: Float[Array, "band"]
    c1: Float[Array, "band"]
    c2: Float[Array, "band"]


# proposed (radiomax)
class EmpiricalCoefficients(eqx.Module):
    """Fitted coefficients of one formula for one sensor (split-window, BRDF, …)."""

    coefficients: Float[Array, "coef"]
    formula: str = eqx.field(static=True)
    sensor: str = eqx.field(static=True)

Results are plain arrays: radiance, reflectance and brightness temperature come back with the shape they went in.

Types to adhere to¶

radiomax reads no files. Readers and Operators load sensor metadata and pass arrays:

# Shapes: B bands; W wavelengths; Ny × Nx pixels

# user: a reader turns file metadata into coefficient records
calibration = rmx.radiometry.BandCalibration(gain=gain, offset=offset)  # (B,), (B,)
# proposed (radiomax)
# JD from a UTC timestamp, on the host, once per scene
julian_date = rmx.geometry.julian_date("2026-11-05T10:30:00Z")  # → () float64

Units are SI at the boundary. Himawari’s radiances are per micrometre, so its reader converts (the factor of 10⁶ its brightness-temperature code applies today) before calling radiomax.

Algorithms¶

Solar position and Earth–Sun distance¶

The low-precision formulas of the Astronomical Almanac, about 0.01° between 1950 and 2050. With JD the Julian date and n = JD − 2451545.0 the days since J2000:

1. Mean longitude L and mean anomaly g from n (above).

2. Ecliptic longitude λ☉ = L + 1.915° sin g + 0.020° sin 2g, obliquity ε = 23.439° − 4 × 10⁻⁷° n.

3. Declination δ = asin(sin ε sin λ☉) and right ascension α = atan2(cos ε sin λ☉, cos λ☉).

4. Hour angle ω = θ_GMST(n) + λ − α, with λ the pixel’s longitude.

5. cos θₛ = sin φ sin δ + cos φ cos δ cos ω; the azimuth from atan2 of the east and north components.

6. Earth–Sun distance d = 1.00014 − 0.01671 cos g − 0.00014 cos 2g, in AU.

Cost: O(1) per pixel, about twenty transcendental calls. NREL’s SPA, behind an argument, reaches about 0.0003° at many times the cost (design decision 6).

Pitfalls. A Julian date near 2.46 × 10⁶ has a float32 resolution of about a quarter of a day, so n is formed in float64 before any cast. Refraction is ignored, so angles within a degree of the horizon are off by up to half a degree. geoops’ distance today is d = 1 − 0.01673 cos(0.0172 (t − 4)) on the day of year t, which differs from the formula above by up to about 10⁻⁴ AU; tests state which one they compare to.

Viewing angles from the satellite’s position¶

With p the pixel and s the satellite in Earth-centred coordinates from geophysix, the line of sight is ℓ = s − p. At geodetic latitude φ and longitude λ the local frame is

1. ℓ = s − p and its norm.

2. cos θᵥ = (ℓ · û) / |ℓ|.

3. view azimuth = atan2(ℓ · ê, ℓ · n̂).

Cost: O(1) per pixel. Pitfalls: the normal must be the geodetic one, not the geocentric one, or θᵥ is off by up to 0.19° at mid-latitudes; ℓ · û < 0 means the satellite is below the horizon (off the limb) and the pixel is masked.

Band solar irradiance¶

1. Bring the solar spectrum onto the response’s wavelengths (a functax interpolation plan), or both onto a common grid.

2. Build an integrix rule on that grid (trapezoid by default).

3. Two rule applications per band and a division.

Cost: O(B W). Pitfalls: a response sampled more coarsely than the solar spectrum’s lines can bias E₀ where the solar spectrum has strong lines, so the common grid is the finer of the two; a response with negative tails (some published SRFs) is clipped at zero first.

Reflectance¶

Cost: O(B N). Pitfalls: cos θₛ → 0 at the terminator makes ρ blow up, so pixels past a threshold (85° by default) are masked rather than divided; geoops today evaluates θₛ once at the scene centre, which is wrong by more than a degree across a wide swath, and radiomax evaluates it per pixel.

Planck’s law and its band inverses¶

The two published band forms:

1. GOES (PlanckLinear): T = (fk₂ / ln(fk₁ / L + 1) − b₁) / b₂.

2. Himawari (PlanckQuadratic): Tₑ from the central wavelength, then T = c₀ + c₁ Tₑ + c₂ Tₑ².

3. From a full response: invert the band-integrated Planck function by a few Newton steps on T, starting from the central-wavelength Tₑ.

Cost: O(B N), plus a fixed number of Newton steps in form 3. Pitfalls: L ≤ 0 (noise in a cold, dark band) makes the logarithm undefined, so those pixels become NaN, as both geoproducts readers do today; radiance per micrometre against per metre changes the argument by 10⁶; dT/dL is finite for L > 0, so the inverse differentiates.

Geometric air-mass factor¶

Cost: O(N). Pitfalls: the plane-parallel form overestimates the path past about 75°, as radtrax’s docstring already warns; a spherical (Chapman) correction would still be geometry, but radtrax’s roadmap puts it with its v1 forward model, so the two must agree on where it lives before either writes it.

Parallax from geostationary orbit¶

A target at height H above a pixel at p is seen along the ray from the satellite s through the elevated point q = p + H û. Its apparent position is where that ray meets the Earth’s surface:

with |·|ₑ the ellipsoid’s norm (x²/a² + y²/a² + z²/b²), a quadratic in t.

1. q from p, H and the geodetic normal.

2. Solve the quadratic for the near root t₁; no real root means off the limb.

3. The apparent point s + t₁ (q − s), back to lat/lon.

4. Resample the field from the apparent to the true position (a functax plan).

Cost: O(N) plus the resampling. Pitfalls: geoops solves this on a sphere whose radius defaults to the equatorial 6,378,137 m, which shifts high-latitude pixels; radiomax takes the ellipsoid from geophysix.

Ross–Li BRDF kernels¶

With θₛ, θᵥ and the relative azimuth ϕ, the phase angle ξ satisfies cos ξ = cos θₛ cos θᵥ + sin θₛ sin θᵥ cos ϕ, and the volumetric (RossThick) kernel is

The geometric kernel (LiSparse-Reciprocal) uses the crown shape ratios h/b = 2 and b/r = 1, as MODIS does. The reflectance is R = f_iso + f_vol K_vol + f_geo K_geo, with f fitted per band.

Cost: O(N) per band. Pitfalls: the LiSparse overlap term takes an arccos of a quantity that leaves [−1, 1] by round-off near the hot spot, so it is clipped; NBAR evaluates the same kernels at nadir view and a reference sun angle.

Lee despeckling¶

With μ and σ² the local mean and variance over a w × w window, C_u the speckle’s coefficient of variation (about 1 / √L for L looks) and Cᵢ = σ / μ the local one:

1. Box sums of x and x² by cumulative sums, so the cost does not grow with w.

2. μ, σ² and Cᵢ per pixel.

3. The weight and the update.

Cost: O(N) per image. Pitfalls: Cᵢ → 0 in flat regions, so the weight is computed as max(0, 1 − C_u² μ² / σ²) with σ² floored; edges use reflected padding. geoops’ default C_u = 0.523 corresponds to about 3.7 looks. Refined Lee and Frost follow the same pattern with directional sub-windows and an exponential kernel.

JAX compatibility¶

Static: window sizes, formula and sensor names, quadrature grids. Traced: every angle, radiance and coefficient. Every function keeps the input’s dtype, except that the Julian date and n are formed in float64.

Testing¶

Every public forward function gets the four tests of the house style.

External dependencies¶

radiomax imports no SRF tables, no file readers and no pysolar: geoops’ solar functions reach pysolar through georeader today, and that path goes.

Where it stands¶

Read on 5 Nov 2026 at geotoolz d2f53f4, radtrax ed04971 and plumax 28f257e.

geoops radiometry/_src/array.py. dn_to_radiance, radiance_to_dn and dn_to_reflectance (a linear decode with scale and offset per band); bt_from_radiance, the two-coefficient form T = K₂ / ln(K₁ / L + 1) used for Landsat, without GOES’s b₁ and b₂; dos1, a dark-object subtraction by np.nanpercentile. All NumPy. They move; the display helpers in the same file (min_max_normalize, percentile_clip, gamma_correct) are not radiometry and stay in geoops.

geoops radiometry/_src/solar.py. earth_sun_distance_correction_factor, compute_sza and observation_date_correction_factor are thin wrappers over georeader.reflectance, which uses pysolar. They take a datetime and return one float per scene, evaluated at the scene centre. radiomax replaces them with per-pixel, traceable versions.

geoops radiometry/_src/operators.py. The Operators, including RadianceToReflectance over georeader.reflectance.radiance_to_reflectance, stay in geoops and call radiomax.

geoops geom/_src/array.py and operators.py. bowtie_detector_positions and scan_pixel_growth give the MODIS and VIIRS scan geometry; they move. GeostationaryParallaxCorrect does its ray intersection inside _apply, on NumPy only, on a sphere of 6,378,137 m; its maths moves and the Operator stays.

geoops restore/_src/array.py. despeckle_lee (window 7, C_u = 0.523), despeckle_refined_lee and despeckle_frost (damping 2); they move.

geoproducts. himawari/reader.py _brightness_temperature applies Planck at the central wavelength and then the quadratic c₀ + c₁ T + c₂ T² from the file header; goes/l1b.py applies (fk₂ / ln(fk₁ / L + 1) − b₁) / b₂. Both move as PlanckQuadratic and PlanckLinear; the readers keep reading the coefficients. In _src/geostationary.py, scan_angle_transform and on_earth are navigation maths and move; geos_crs and FixedGrid build CRSs and stay.

radtrax and plumax. radtrax’s atmosphere/geometry.py already consolidated the air-mass factor (radtrax#27), and slab/config.py’s property calls it. plumax still has two copies, hapi_lut/beers.py and radtran/config.py. So there are three implementations, not the four the Design tab counts.