Spatial Epi Public Health GIS

Remote Sensing Indices for Health Studies

This guide is part of the Environmental Exposure & Risk Surface Modeling section. Where air pollutants require a sparse network and a model, land-surface properties are measured everywhere by satellites — which moves the difficulty from prediction to interpretation.

Concept & Epidemiological Alignment

A remote sensing index is an arithmetic combination of reflectance or radiance bands, designed so that the result correlates with a physical property of the surface. The normalized difference vegetation index contrasts near-infrared with red reflectance, exploiting the fact that living vegetation reflects strongly in the near-infrared and absorbs red. Land surface temperature is derived from thermal infrared radiance with an emissivity correction. Imperviousness is usually a classification product rather than a band ratio.

The critical epidemiological point is that each index measures the surface, not the human experience of it. Land surface temperature is the radiometric temperature of roofs, roads and canopy tops as seen from above; it is correlated with, but systematically different from, the near-surface air temperature a person experiences, often by several degrees and with a bias that depends on land cover. A greenness index counts photosynthetic activity within a cell, whether or not anyone can reach it — a golf course, a median strip and a public park contribute alike.

Three assumptions therefore travel with every index used as an exposure. That the measured surface property is a usable proxy for the human-relevant quantity; that the scale of the product matches the scale of the exposure; and that the compositing window matches the health-relevant window. Each is testable, and none is automatic.

Where these indices are strong is in gradient and ranking. Comparing the greenness of neighbourhoods, or identifying the hottest blocks in a city for a heat plan, is work these products do reliably at fine resolution and low cost. Where they are weak is absolute interpretation: a coefficient per index unit is difficult to interpret physically, and a study that reports a risk per 0.1 greenness units has said something hard for a reader to act on.

Index Selection

Index Bands / product Health use Typical resolution Caveat
Vegetation difference index Near-infrared and red Greenness, mental health, birth outcomes 10–30 m Saturates over dense canopy; seasonal
Enhanced vegetation index Near-infrared, red, blue Greenness in dense canopy 10–250 m Less comparable with older literature
Land surface temperature Thermal infrared with emissivity Heat exposure, heat vulnerability 30–1000 m Surface not air temperature; clear-sky only
Imperviousness Classification product Urban form, runoff, heat 30 m Class definitions change between vintages
Tree canopy cover Classification product Shade, usable greenness 1–30 m Expensive; irregular update cycle
Night-time lights Low-light imaging Urbanicity, electrification 500 m Saturates in city cores; sensor changes
Two Gaps Between the Index and the Exposure Two panels. The left panel plots surface temperature against air temperature for three land-cover types: over asphalt the surface reads far hotter than the air, over grass the two are close, and under canopy the surface reads cooler than the air. The gap is therefore not a constant offset and cannot be corrected by subtracting one number. The right panel shows three cells with identical greenness values arising from a public park, a golf course and a highway median, illustrating that the index counts photosynthetic activity without regard to whether anyone can reach it. The index is the surface; the exposure is the person Surface temperature is not air temperature asphalt grass canopy air temperature → surface temperature Identical greenness, different exposure park golf course highway median index value: 0.62 in all three cells accessible, restricted and unreachable Pair greenness with an accessibility layer where the hypothesis is about use Neither gap is a defect of the sensor — both are questions the study design has to answer

Spatial Data Prerequisites

  • Imagery or derived product with documented sensor, processing level and calibration, covering the study period rather than a convenient year.
  • Coordinate reference system: rasters reprojected once to the analysis system with bilinear resampling for continuous indices and nearest neighbour for classifications.
  • Cloud and quality masks retained as their own band, never silently applied and discarded.
  • Scale match: the product resolution must be finer than the geographic variation the study intends to estimate. A 1 km thermal product cannot describe block-level heat.
  • Seasonal definition: which months compose the index, fixed before extraction and applied identically across years.
  • Sensor continuity: where a study spans a sensor change, a documented cross-calibration or an explicit break in the series.

Production Implementation

The routine below computes a vegetation index and a cloud-masked seasonal composite from a band stack, and extracts area summaries for a set of polygons — the common pattern behind a neighbourhood greenness exposure.

# Vegetation index, cloud-masked seasonal composite, and area summaries.
# Pinned: rasterio==1.3.10, numpy==1.26.4, geopandas==0.14.4, pandas==2.2.2
import logging

import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio
from rasterio.mask import mask as rio_mask

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

MIN_CLEAR_FRACTION = 0.5     # a composite cell needs this share of clear observations


def vegetation_index(nir: np.ndarray, red: np.ndarray,
                     eps: float = 1e-6) -> np.ndarray:
    """Normalized difference of near-infrared and red reflectance.

    Inputs must be surface reflectance with scale factors ALREADY applied; using
    raw digital numbers silently changes the index because the ratio is not
    invariant to an additive offset.
    """
    nir = nir.astype("float32")
    red = red.astype("float32")
    denom = nir + red
    out = np.where(np.abs(denom) < eps, np.nan, (nir - red) / np.where(denom == 0, eps, denom))
    # Physically valid range; anything outside is a processing artefact.
    out[(out < -1.0) | (out > 1.0)] = np.nan
    return out


def seasonal_composite(scene_index: list[np.ndarray],
                       scene_clear: list[np.ndarray]) -> tuple[np.ndarray, np.ndarray]:
    """Per-cell median over clear observations, plus the clear-observation count."""
    stack = np.stack([np.where(c, v, np.nan) for v, c in zip(scene_index, scene_clear)])
    clear_n = np.sum(~np.isnan(stack), axis=0)
    with np.errstate(all="ignore"):
        comp = np.nanmedian(stack, axis=0)
    enough = clear_n >= max(1, int(MIN_CLEAR_FRACTION * len(scene_index)))
    comp = np.where(enough, comp, np.nan)
    log.info("composite from %d scenes: median clear obs=%d, cells below threshold=%d",
             len(scene_index), int(np.median(clear_n)), int((~enough).sum()))
    return comp.astype("float32"), clear_n.astype("int16")


def area_summaries(raster_path: str, areas: gpd.GeoDataFrame,
                   id_col: str = "geoid", band: int = 1) -> pd.DataFrame:
    """Mean, median and valid-cell count of the index within each area."""
    rows = []
    with rasterio.open(raster_path) as src:
        ar = areas.to_crs(src.crs) if areas.crs != src.crs else areas
        for rec in ar.sort_values(id_col).itertuples():
            try:
                arr, _ = rio_mask(src, [rec.geometry.__geo_interface__],
                                  crop=True, filled=False, indexes=[band])
            except ValueError:            # polygon does not overlap the raster
                rows.append({id_col: getattr(rec, id_col), "mean": np.nan,
                             "median": np.nan, "n_cells": 0})
                continue
            layer = arr[0]
            valid = layer.compressed()
            valid = valid[np.isfinite(valid)]
            rows.append({id_col: getattr(rec, id_col),
                         "mean": float(valid.mean()) if valid.size else np.nan,
                         "median": float(np.median(valid)) if valid.size else np.nan,
                         "n_cells": int(valid.size)})
    df = pd.DataFrame(rows)
    log.info("area summaries: %d areas, %d with no valid cells",
             len(df), int((df["n_cells"] == 0).sum()))
    return df

Two details in that code are where remote-sensing exposures usually go wrong. The index is computed on surface reflectance with scale factors applied — a normalized difference is invariant to a multiplicative scale but not to an additive offset, so computing it from raw digital numbers gives a number that looks like an index and is not one. And the composite keeps a clear-observation count per cell, which is the only way to know later whether a low greenness value describes a genuinely sparse area or a cell that was cloudy all summer.

The area summary returns both mean and median. For a skewed index over a heterogeneous area — a tract containing a lake, or a large parking structure — the two differ materially, and which one the health model should use is a study design decision rather than a default.

Parameter Selection & Tuning

Compositing window. A summer median is the usual choice for greenness because it captures peak canopy and avoids leaf-off confusion, but it is the wrong window for a study of winter exposure. Fix the window from the hypothesis and apply it identically in every year.

Composite statistic. A median is robust to residual cloud that survived the mask; a maximum-value composite is the traditional choice for vegetation indices and is biased upward by exactly that residual cloud edge brightness. Prefer the median unless comparability with older work forces otherwise.

Buffer radius for a point exposure. Greenness around a residence is usually summarised within 100 m to 1 km, and the choice changes the exposure substantially. It is a radius decision with the same status as any other, and belongs to the reasoning in Choosing a Buffer Radius for Proximity Exposure.

Clear-observation threshold. Requiring half the scenes to be clear is a reasonable default in temperate regions and impossible in persistently cloudy ones. Set it from the actual distribution of clear counts and report the fraction of cells it nulls.

Emissivity for land surface temperature. Thermal products depend on an emissivity estimate that is usually derived from land cover. Using a constant emissivity across a city introduces a systematic difference between vegetated and built areas — which is exactly the contrast a heat study is measuring.

Edge Cases & Failure Modes

Cloud masks that under-mask. Thin cirrus and cloud shadow frequently survive automatic masks and depress a vegetation index or a temperature. A composite median absorbs occasional contamination; a single-date extraction does not.

Saturation. Vegetation difference indices saturate over dense canopy, so a forest and a dense park both read near the maximum and the contrast between them vanishes. Where the study population lives in high-canopy areas, an enhanced index or canopy cover product discriminates better.

Sensor and processing-collection changes. A reprocessing collection can shift values across an entire archive. A study spanning such a change must cross-calibrate on overlapping dates or treat it as a break.

Mixed pixels at 30 m. A single cell containing a house, a garden and a road returns one blended value, which is the correct measurement and a poor description of any of the three. For within-parcel questions, the resolution is simply insufficient.

Areas with no valid cells. Small urban tracts after cloud masking can end up with none. They must be reported as missing, not filled with a city mean.

Thermal products and clear-sky bias. Land surface temperature is only retrievable under clear sky, so a summer composite systematically describes the hottest, clearest days. That is often the exposure of interest, but it is not the summer mean, and calling it one is wrong.

Clear Observations Are Not Evenly Distributed A schematic map of a study area divided into bands with different clear-observation counts over one summer. A coastal band has few clear scenes because of persistent marine layer, a central band has many, and a mountain band has few because of orographic cloud. A threshold requiring half the scenes to be clear nulls the coastal and mountain bands, which means the greenness exposure is missing in two specific and non-random parts of the study area. A note states that filling those areas with a regional mean would place the missing data exactly where the population differs most. Cloud is geographic, so missing composite cells are too coastal 4 clear of 18 interior valley 15 clear of 18 mountain 6 clear of 18 nulled retained nulled Report the clear-observation count per area; a filled value hides which places were never seen

Saturation is the edge case most likely to mislead a greenness study in a well-treed city, and the two index families diverge exactly where it matters:

Where the Common Index Stops Discriminating Index value plotted against true leaf area for two vegetation indices. The normalized difference index rises quickly at low vegetation and flattens near its maximum as canopy becomes dense, so a dense park and a forest read almost identically. The enhanced index rises more slowly and keeps increasing through the dense range. A shaded band marks the dense-canopy range typical of leafy residential neighbourhoods, where only the enhanced index still separates them. In leafy neighbourhoods the common index has already flattened out dense canopy range normalized difference enhanced index true leaf area index value Choose the index from where your study population actually lives on this curve

Compliance & Audit Controls

Remote-sensing inputs are public, so the audit burden is about reproducibility rather than disclosure until the index meets an address. Store the sensor and collection identifier, the scene list with acquisition dates, the processing level and whether scale factors were applied, the cloud-mask algorithm and its version, the compositing window and statistic, the clear-observation threshold and the fraction of cells it nulled, the reprojection and resampling method, and the emissivity treatment for any thermal product.

Publish the derived index raster where policy allows, with the clear-observation count as a band. It is a public-good artifact that other agencies can validate and reuse, and reuse is cheaper than re-derivation for everyone. The protected step remains the join to residences, under the controls in Assigning Exposure to Residential Locations.

Production Implementation Checklist

Frequently Asked Questions

Is a vegetation index a good measure of greenspace exposure? It is a good measure of photosynthetic activity within a cell, which is a useful proxy for neighbourhood greenness and a poor one for usable greenspace. If the hypothesis concerns physical activity or restorative use, pair the index with a park accessibility layer, because a golf course and a public park have the same index value and very different meanings.

Can I use land surface temperature as a heat exposure? Yes, with the right framing. It is the clear-sky radiometric temperature of the surface, which correlates strongly with heat vulnerability at neighbourhood scale and differs from near-surface air temperature by an amount that varies with land cover. Report it as what it is rather than converting it to air temperature with a single offset.

What resolution do I need? Finer than the variation you are testing. Block-level heat analysis needs a 30 m thermal product or a downscaled one; a county-level comparison is well served by 1 km. Using a coarse product for a fine question produces a null result that says nothing about the hypothesis.

How do I handle a study period that spans a sensor change? Cross-calibrate on dates where both sensors observed the same area, and report the adjustment. Where no overlap exists, treat the change as a break in the series and do not estimate a trend across it — a step change in the sensor is indistinguishable from a step change in the world.

Should I use a single date or a composite? A composite, in almost every case. A single date carries whatever cloud, haze and phenological state happened to exist that day, and those vary systematically across a large study area. The exception is an acute event — a heat wave, a fire — where the date is the exposure.