Radiometric Calibration in Python

import exiftool


REQUIRED = ["ExposureTime", "ISOSpeed", "BlackLevel", "BandName",
            "CentralWavelength", "Irradiance"]


def frame_metadata(path: str) -> dict:
    """Pull the calibration terms from a frame's metadata, failing loudly.

    A missing tag is not something to default around: a frame without an
    exposure time cannot be normalised, and silently assuming one produces a
    reflectance value that is wrong by whatever the real exposure was.
    """
    with exiftool.ExifToolHelper() as et:
        tags = et.get_metadata(path)[0]

    out = {}
    for key in REQUIRED:
        hit = next((v for k, v in tags.items() if k.endswith(key)), None)
        if hit is None:
            raise KeyError(f"{path}: metadata tag {key} is absent — "
                           "this frame cannot be calibrated")
        out[key] = hit
    out["path"] = path
    return out

import numpy as np


def normalise_frame(raw: np.ndarray, meta: dict) -> np.ndarray:
    """Remove the sensor offset and the camera's exposure choices.

    Clipping at zero after the dark subtraction is deliberate: genuinely dark
    pixels sit at or slightly below the black level through noise, and a
    small negative radiance is physically meaningless but numerically
    destructive once it reaches a ratio index.
    """
    dn = raw.astype(np.float64)
    dark = float(meta["BlackLevel"])
    exposure = float(meta["ExposureTime"])
    gain = float(meta["ISOSpeed"]) / 100.0

    normalised = (dn - dark) / (exposure * gain)
    return np.maximum(normalised, 0.0)

import numpy as np


def vignette_from_coefficients(shape: tuple[int, int], centre: tuple[float, float],
                               coeffs: list[float]) -> np.ndarray:
    """Evaluate a vendor's radial vignette polynomial over the frame.

    The polynomial gives the correction factor as a function of radius from
    the optical centre — which is not the image centre, and using the image
    centre instead leaves a visible gradient across the frame.
    """
    h, w = shape
    cy, cx = centre
    yy, xx = np.mgrid[0:h, 0:w]
    r = np.hypot(xx - cx, yy - cy)
    factor = np.zeros_like(r, dtype=np.float64)
    for power, c in enumerate(coeffs):
        factor += c * r ** power
    return 1.0 / np.maximum(factor, 1e-6)


def vignette_from_flat_field(flat: np.ndarray, *, order: int = 4) -> np.ndarray:
    """Fit a smooth radial model to a flat-field image.

    Fitting a smooth model rather than using the flat field directly is what
    keeps sensor noise and dust specks out of every calibrated frame.
    """
    h, w = flat.shape
    cy, cx = h / 2, w / 2
    yy, xx = np.mgrid[0:h, 0:w]
    r = np.hypot(xx - cx, yy - cy).ravel()
    values = flat.ravel().astype(np.float64)
    good = np.isfinite(values) & (values > 0)

    design = np.vander(r[good], order + 1, increasing=True)
    coeffs, *_ = np.linalg.lstsq(design, values[good], rcond=None)
    model = np.vander(r, order + 1, increasing=True) @ coeffs
    model = model.reshape(h, w)
    return model / np.max(model)

import numpy as np


def reflectance(radiance: np.ndarray, *, irradiance: float,
                panel_radiance: float, panel_reflectance: float,
                panel_irradiance: float) -> np.ndarray:
    """Convert corrected radiance to surface reflectance.

    The panel gives a radiance-to-reflectance ratio measured under the panel
    image's own irradiance. Scaling that ratio by the ratio of irradiances is
    what transfers it to a frame taken under different light — which is the
    entire purpose of the downwelling sensor.
    """
    if panel_radiance <= 0 or panel_irradiance <= 0:
        raise ValueError("panel radiance and irradiance must be positive")
    factor = panel_reflectance / panel_radiance
    illumination_ratio = panel_irradiance / max(irradiance, 1e-9)
    return radiance * factor * illumination_ratio

import numpy as np
import rasterio


def calibrate_frame(raw_path: str, out_path: str, *, meta: dict,
                    vignette: np.ndarray, panel: dict) -> dict:
    """Full calibration of one band frame, written as float reflectance."""
    with rasterio.open(raw_path) as src:
        raw = src.read(1)
        profile = src.profile

    radiance = normalise_frame(raw, meta) / np.maximum(vignette, 1e-6)
    refl = reflectance(radiance,
                       irradiance=float(meta["Irradiance"]),
                       panel_radiance=panel["radiance"],
                       panel_reflectance=panel["reflectance"],
                       panel_irradiance=panel["irradiance"])

    # Physically, reflectance is in [0, 1]; values above about 1.1 indicate a
    # specular return or a calibration error, and are worth counting.
    over = float(np.count_nonzero(refl > 1.1) / refl.size)

    profile.update(dtype="float32", count=1, nodata=np.nan, compress="deflate")
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(refl.astype("float32"), 1)
        dst.set_band_description(1, str(meta["BandName"]))
        dst.update_tags(1, WAVELENGTH_NM=str(meta["CentralWavelength"]),
                        PANEL_REFLECTANCE=str(panel["reflectance"]))

    return {"band": meta["BandName"], "over_unity_fraction": over,
            "mean_reflectance": float(np.nanmean(refl))}

import numpy as np


def panel_factor_series(pre: dict, post: dict, frame_times: np.ndarray,
                        irradiance: np.ndarray | None = None) -> np.ndarray:
    """Per-frame panel factor, interpolated between the two panel captures.

    With an irradiance record, the interpolation follows the measured light
    rather than assuming a linear drift — which matters on a flight where the
    cloud cleared halfway through.
    """
    t0, t1 = pre["time_s"], post["time_s"]
    f0 = pre["reflectance"] / pre["radiance"]
    f1 = post["reflectance"] / post["radiance"]

    if irradiance is None or t1 <= t0:
        weight = np.clip((frame_times - t0) / max(t1 - t0, 1e-9), 0.0, 1.0)
    else:
        # Follow the irradiance record's own trajectory between the endpoints.
        i0, i1 = irradiance[0], irradiance[-1]
        weight = np.clip((irradiance - i0) / max(i1 - i0, 1e-9), 0.0, 1.0)
    return f0 * (1 - weight) + f1 * weight

def audit_flight_metadata(frames: list[dict]) -> dict:
    """Consistency checks across a whole flight's frames before calibrating."""
    exposures = [float(f["ExposureTime"]) for f in frames]
    bands = {f["BandName"] for f in frames}
    problems = []

    spread = max(exposures) / max(min(exposures), 1e-9)
    if spread > 8:
        problems.append(f"exposure varies by {spread:.0f}× across the flight — "
                        "auto-exposure was active on a radiometric capture")
    if len(bands) != 1:
        problems.append(f"frames report {len(bands)} band names: {sorted(bands)}")
    if any(float(f["Irradiance"]) <= 0 for f in frames):
        problems.append("some frames report zero irradiance — the downwelling "
                        "sensor was obstructed or not recording")
    return {"frames": len(frames), "exposure_spread": spread,
            "problems": problems}

import numpy as np


def calibration_report(reflectance_by_band: dict[str, np.ndarray],
                       known_target_mask: np.ndarray,
                       known_reflectance: dict[str, float],
                       tol: float = 0.03) -> dict:
    """Compare calibrated values over a known target against its true values."""
    findings = []
    measured = {}
    for band, arr in reflectance_by_band.items():
        vals = arr[known_target_mask & np.isfinite(arr)]
        if vals.size < 100:
            findings.append(f"{band}: too few target pixels to judge")
            continue
        got = float(np.median(vals))
        measured[band] = got
        expected = known_reflectance.get(band)
        if expected is not None and abs(got - expected) > tol:
            findings.append(
                f"{band}: measured {got:.3f} against a true {expected:.3f} — "
                f"off by {got - expected:+.3f}")
    return {"measured": measured, "findings": findings,
            "ok": not findings}