Spatial Epi Public Health GIS

Deriving Land Surface Temperature for Heat Vulnerability

This guide, part of Remote Sensing Indices for Health Studies, walks the derivation of a surface temperature exposure for a heat vulnerability assessment, and is blunt about the three places where the resulting number stops describing what a heat plan needs.

Problem Context & Constraints

A thermal sensor measures radiance. Converting that to a temperature requires an emissivity, because a surface at a given temperature emits differently depending on what it is made of. Asphalt, dry soil, water and canopy have materially different emissivities, and the common shortcut of assuming a single value across a city introduces a bias that is correlated with land cover — which means it is correlated with exactly the urban–vegetated contrast a heat vulnerability study is measuring. The bias does not average out; it inflates or deflates the specific comparison the study exists to make.

Two further gaps separate the derived number from heat risk. Overpass time. Polar-orbiting sensors observe at a fixed local time, commonly late morning or early afternoon for the daytime pass and after midnight for the nighttime one. Heat mortality is driven substantially by night-time temperatures, when the body cannot recover, and the urban heat island is strongest at night. A daytime scene describes the hottest surfaces; a night scene describes the ones that fail to cool. They rank neighbourhoods differently, and the night ranking is closer to the health-relevant one.

Clear-sky retrieval. Thermal retrieval requires an unobstructed view, so a composite of available scenes is a composite of clear days. Clear days in summer are hot days. The resulting statistic is not the summer mean surface temperature; it is the mean over the clearest and generally hottest subset, and describing it as a seasonal average is wrong in a direction that matters.

Day and Night Rank Neighbourhoods Differently Three neighbourhood types plotted twice, once for the daytime overpass and once for the night overpass. In the daytime the industrial area with large flat roofs is by far the hottest, an older dense residential area is intermediate, and a leafy suburb is coolest. At night the ranking changes: the dense residential area retains the most heat because of its thermal mass and lack of sky view, the industrial roofs cool quickly, and the leafy suburb is coolest again. A note records that night-time retention is closer to the quantity that drives heat mortality. Which neighbourhood is hottest depends on which overpass you used Daytime overpass industrial dense residential leafy suburb Night overpass industrial dense residential leafy suburb Night-time heat retention is closer to what drives mortality a daytime-only vulnerability map will prioritise the wrong blocks

Prerequisites

  • python 3.11
  • numpy 1.26.4, rasterio 1.3.10, pandas 2.2.2, geopandas 0.14.4
  • Thermal band radiance or brightness temperature with the sensor’s calibration constants
  • A vegetation index from the same scene, used to estimate emissivity
  • Cloud and quality masks for each scene
  • Scene metadata including acquisition local time and processing collection

Step-by-Step Solution

Estimate emissivity per cell from the vegetation fraction, apply the single-channel correction to brightness temperature, and composite separately for day and night.

# Land surface temperature from brightness temperature with per-cell emissivity.
# Pinned: numpy==1.26.4, rasterio==1.3.10, pandas==2.2.2
import logging

import numpy as np

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("exposure.lst")

# Representative emissivities. Replace with values for your sensor's band and
# your region's surfaces — these are defaults, not constants of nature.
EMIS_VEG = 0.985
EMIS_SOIL = 0.960
EMIS_WATER = 0.991
LAMBDA_M = 10.9e-6          # effective wavelength of the thermal band, metres
RHO = 1.438e-2              # h*c/sigma, metres kelvin


def vegetation_fraction(ndvi: np.ndarray, ndvi_soil: float = 0.2,
                        ndvi_veg: float = 0.6) -> np.ndarray:
    """Fractional vegetation cover from a vegetation index, clipped to [0, 1]."""
    fv = ((ndvi - ndvi_soil) / (ndvi_veg - ndvi_soil)) ** 2
    return np.clip(fv, 0.0, 1.0)


def emissivity(ndvi: np.ndarray, water_mask: np.ndarray | None = None) -> np.ndarray:
    """Per-cell emissivity from vegetation fraction, with water handled separately."""
    fv = vegetation_fraction(ndvi)
    e = EMIS_VEG * fv + EMIS_SOIL * (1.0 - fv)
    if water_mask is not None:
        e = np.where(water_mask, EMIS_WATER, e)
    return e.astype("float32")


def land_surface_temperature(bt_kelvin: np.ndarray, ndvi: np.ndarray,
                             water_mask: np.ndarray | None = None) -> np.ndarray:
    """Single-channel emissivity correction of brightness temperature.

    LST = BT / (1 + (lambda * BT / rho) * ln(emissivity))
    Returns degrees Celsius. Inputs must be from the SAME scene and grid.
    """
    if bt_kelvin.shape != ndvi.shape:
        raise ValueError("brightness temperature and vegetation index grids differ")
    e = emissivity(ndvi, water_mask)
    with np.errstate(invalid="ignore", divide="ignore"):
        lst_k = bt_kelvin / (1.0 + (LAMBDA_M * bt_kelvin / RHO) * np.log(e))
    lst_c = lst_k - 273.15
    # Physically implausible values are retrieval failures, not cold neighbourhoods.
    lst_c[(lst_c < -40.0) | (lst_c > 75.0)] = np.nan
    const_e = float(np.nanmean(e))
    log.info("LST derived: mean emissivity=%.4f range=%.4f–%.4f, valid cells=%d",
             const_e, float(np.nanmin(e)), float(np.nanmax(e)),
             int(np.isfinite(lst_c).sum()))
    return lst_c.astype("float32")


def day_night_composites(scenes: list[dict]) -> dict:
    """Separate composites by overpass. scenes: [{'lst':arr,'clear':arr,'pass':'day'|'night'}]"""
    out = {}
    for which in ("day", "night"):
        sel = [s for s in scenes if s["pass"] == which]
        if not sel:
            log.warning("no %s scenes available — the composite will be one-sided", which)
            continue
        stack = np.stack([np.where(s["clear"], s["lst"], np.nan) for s in sel])
        n_clear = np.sum(~np.isnan(stack), axis=0)
        with np.errstate(all="ignore"):
            out[which] = {"mean": np.nanmean(stack, axis=0).astype("float32"),
                          "p95": np.nanpercentile(stack, 95, axis=0).astype("float32"),
                          "n_clear": n_clear.astype("int16"),
                          "n_scenes": len(sel)}
        log.info("%s composite from %d scenes, median clear obs=%d",
                 which, len(sel), int(np.median(n_clear)))
    return out

Deriving emissivity from the scene’s own vegetation index, rather than from a constant, is what removes the land-cover-correlated bias. It is not exact — the vegetation-fraction method is an approximation and water and bare rock need explicit handling — but it is far better than a constant, and the log line reporting the emissivity range makes the correction’s magnitude visible rather than hidden.

The day and night composites are kept separate on purpose. Merging them produces an average of two different physical phenomena and is defensible only if the study genuinely has no timing hypothesis, which for heat vulnerability it almost never does. Both a mean and a high percentile are returned, because a heat plan is usually driven by the extreme rather than the average.

What a Constant Emissivity Costs the Contrast Two bars showing the derived temperature difference between a built area and a vegetated area. Using a constant emissivity across the city, the derived contrast is smaller than the truth because the built surface's lower emissivity is not accounted for. Using a per-cell emissivity derived from the scene's own vegetation fraction, the contrast is closer to the truth. Because the error follows land cover, it does not cancel across a city and it shrinks precisely the comparison a heat vulnerability study reports. A constant emissivity biases the built-to-vegetated contrast, not the mean constant emissivity 4.1 °C per-cell emissivity 6.3 °C contrast measured by ground sensors Derived built-minus-vegetated surface temperature difference in one city, one summer

Validation & Edge Cases

Water bodies. High emissivity and high thermal inertia make water look anomalously cool by day and warm at night. Mask it or handle it explicitly; leaving it in a neighbourhood mean makes waterfront blocks look cooler than they are for residents.

Cloud shadow. Shadowed land is genuinely cooler and is not a retrieval failure, but it is also not representative. Quality masks that flag shadow should be applied alongside the cloud mask.

Gaps in the thermal band. Some sensors have banding or striping artefacts in their thermal channels. These appear as regular linear features in the composite and must be caught visually; no statistical check reliably finds them.

Downscaled products. Thermal bands are frequently coarser than the optical bands and are downscaled using the vegetation index. That makes the temperature partly a function of greenness by construction, which is circular if the study also uses greenness as a covariate. Know whether your product is downscaled and how.

Comparing across years. Surface temperature depends on the weather on the day of the overpass. A year-to-year comparison of absolute values is largely a comparison of weather; a comparison of relative patterns within a city is far more stable, which is why heat vulnerability work usually standardises within scene.

INFO exposure.lst LST derived: mean emissivity=0.9714 range=0.9600–0.9910, valid cells=18442310
INFO exposure.lst day composite from 14 scenes, median clear obs=11
INFO exposure.lst night composite from 9 scenes, median clear obs=7

The standardisation point at the end of that list deserves a picture, because it is what makes a multi-year heat map interpretable at all:

Weather Moves the Level; the Ranking Stays Three neighbourhoods' composite surface temperatures in three different summers. The absolute values rise and fall together by several degrees from year to year because of weather on the overpass days. The order of the three neighbourhoods is the same every year: the dense residential area is always hottest and the leafy suburb always coolest. Standardising within each year's scene removes the weather and leaves the stable spatial pattern that a heat vulnerability index needs. Compare within a scene, not across years 2023 2024 2025 dense residential mixed leafy suburb composite surface temperature; the level moves with weather, the order does not

Compliance Notes

Record the sensor, thermal band and calibration constants used, the emissivity method with its parameters and the realised emissivity range, the scene list with acquisition local times and day or night designation, the cloud and shadow masking, the clear-observation counts, and whether the product was downscaled and by what auxiliary variable.

Describe the product accurately wherever it appears. The honest label is “clear-sky land surface temperature at the sensor’s overpass time, composited over N scenes” — not “summer temperature” and not “air temperature”. A heat vulnerability index built on this layer inherits that description, and a plan that allocates cooling resources on a daytime-only surface metric should say so, because the blocks it prioritises differ from the ones a night-time metric would choose.