Converting Radiometric JPEGs to Temperature Rasters

import io
import struct

import numpy as np
import exiftool
from PIL import Image


def extract_raw_thermal(path: str) -> dict:
    """Pull the raw thermal array and calibration constants from the file.

    The raw stream is usually a 16-bit PNG or a flat TIFF depending on the
    vendor, so it is decoded rather than assumed. Byte order varies too,
    which is why the sanity check on the value range is worth having.
    """
    with exiftool.ExifToolHelper() as et:
        tags = et.get_metadata(path)[0]
        blob = et.execute("-b", "-RawThermalImage", path, raw_bytes=True)

    if not blob:
        raise ValueError(f"{path} has no embedded raw thermal stream — "
                         "this is a colourised image, not a radiometric one")

    raw = np.array(Image.open(io.BytesIO(blob)))
    if raw.dtype != np.uint16:
        raw = raw.astype(np.uint16)

    # Vendors differ in byte order; a plausible sensor range is 5 000–40 000.
    if np.median(raw) > 45000 or np.median(raw) < 2000:
        raw = raw.byteswap()

    def tag(name, default=None):
        hit = next((v for k, v in tags.items() if k.endswith(name)), default)
        if hit is None:
            raise KeyError(f"{path}: thermal parameter {name} is absent")
        return hit

    return {
        "raw": raw,
        "planck": {"r1": float(tag("PlanckR1")), "b": float(tag("PlanckB")),
                   "f": float(tag("PlanckF")), "o": float(tag("PlanckO")),
                   "r2": float(tag("PlanckR2"))},
        "emissivity": float(tag("Emissivity", 0.95)),
        "reflected_temp_c": float(tag("ReflectedApparentTemperature", 20.0)),
        "datetime": tag("DateTimeOriginal"),
    }

import numpy as np


def raw_to_celsius(raw: np.ndarray, planck: dict, *, emissivity: float,
                   reflected_temp_c: float) -> np.ndarray:
    """Invert the camera's Planck calibration, with the reflected term removed.

    Removing the reflected component before inversion is what makes the
    emissivity correction physically meaningful: the sensor measured the sum
    of what the surface emitted and what it reflected, and only the first is
    a property of the surface.
    """
    r1, b, f, o, r2 = (planck[k] for k in ("r1", "b", "f", "o", "r2"))

    def radiance(temp_c: float) -> float:
        return r1 / (r2 * (np.exp(b / (temp_c + 273.15)) - f)) - o

    measured = raw.astype(np.float64)
    surface = (measured - (1.0 - emissivity) * radiance(reflected_temp_c)) / emissivity

    with np.errstate(invalid="ignore", divide="ignore"):
        inner = r1 / (r2 * (surface + o)) + f
        kelvin = b / np.log(np.maximum(inner, 1e-12))

    celsius = (kelvin - 273.15).astype("float32")
    return np.where(np.isfinite(celsius) & (celsius > -80) & (celsius < 600),
                    celsius, np.nan)

import numpy as np
import rasterio


def write_temperature_raster(celsius: np.ndarray, out_path: str, *,
                             transform, crs, emissivity: float,
                             reflected_temp_c: float, source: str) -> None:
    """A temperature raster that states its units and its assumptions."""
    profile = {"driver": "GTiff", "height": celsius.shape[0],
               "width": celsius.shape[1], "count": 1, "dtype": "float32",
               "crs": crs, "transform": transform, "nodata": np.nan,
               "compress": "deflate", "predictor": 3, "tiled": True}

    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(celsius.astype("float32"), 1)
        dst.set_band_description(1, "surface_temperature")
        dst.update_tags(1, UNITS="degrees_celsius",
                        EMISSIVITY=f"{emissivity:.3f}",
                        REFLECTED_TEMP_C=f"{reflected_temp_c:.1f}",
                        SOURCE=source)

import numpy as np


def sanity_check_temperatures(celsius: np.ndarray, *, air_temp_c: float,
                              tolerance_c: float = 30.0) -> dict:
    """Are these temperatures physically plausible for this scene?"""
    v = celsius[np.isfinite(celsius)]
    if v.size == 0:
        return {"problems": ["no finite temperatures"]}

    problems = []
    median = float(np.median(v))
    if abs(median - air_temp_c) > tolerance_c:
        problems.append(f"scene median {median:.1f} °C is {abs(median - air_temp_c):.0f} °C "
                        f"from the air temperature — check byte order and constants")
    if float(v.max() - v.min()) > 120:
        problems.append("range exceeds 120 °C — likely a conversion fault")
    return {"median_c": median, "p05_c": float(np.percentile(v, 5)),
            "p95_c": float(np.percentile(v, 95)), "problems": problems}