Spatial Epi Public Health GIS

Cloud Masking and Compositing for Exposure Rasters

This guide, part of Remote Sensing Indices for Health Studies, covers the step every satellite-derived exposure depends on and nobody writes up: deciding which observations are usable and how to reduce the survivors to one number per cell.

Problem Context & Constraints

Cloud does not remove data at random. It removes it in coherent patches, correlated with terrain, coastlines and season. That has two consequences for an exposure raster. Cells that survive masking in a cloudy region describe a biased subset of the period; and cells that fail to survive cluster geographically, so the missing data have a map.

Residual contamination — cloud that the mask failed to catch — biases indices in a direction that depends on the index. For a vegetation index, cloud is bright in red and relatively less bright in near-infrared, so a contaminated cell reads low, and a maximum-value composite systematically selects the clearest observations, which is why that composite became the tradition. For a thermal product, cloud top is cold, so contamination reads low as well, and a mean composite drags the whole city cooler in cloudier areas. For an aerosol product, thin cloud reads as high aerosol.

Three constraints shape the fix. Quality bits must be decoded, not thresholded. Modern products ship a packed bitmask where individual bits mean cloud, shadow, cirrus, snow and confidence level; comparing the packed integer against a number is a bug that happens to work on some scenes. The mask must be dilated. Cloud edges and shadow are under-detected, and a one-to-three cell dilation removes most residual contamination at a modest cost in coverage. And the clear count must survive into the output, because a composite value without it cannot be distinguished from a guess.

Why the Mask Needs Dilating Two panels of the same scene. On the left the cloud mask covers the cloud core only, leaving a ring of partially contaminated cells at the edge and an undetected shadow offset to one side, both of which pass into the composite carrying biased values. On the right the mask has been dilated by two cells and a shadow mask applied, removing the contaminated ring and the shadow at the cost of a modest reduction in usable cells. A note records that the discarded cells are few but carry the largest errors in the scene. The worst values in a scene sit just outside the mask Mask as delivered cloud contaminated ring survives shadow Dilated, with shadow masked cloud + 2 cells masked a few per cent fewer cells, far less bias Dilation trades coverage for correctness, and the trade is almost always worth taking provided the resulting clear-observation count is carried forward

Prerequisites

  • python 3.11
  • numpy 1.26.4, scipy 1.13.1, rasterio 1.3.10
  • Scene-level index arrays and their quality assurance bands, on one grid
  • The product’s quality-band bit definitions from its user guide — not from memory
  • A declared compositing window and clear-observation threshold

Step-by-Step Solution

Decode the quality bits into a boolean clear mask, dilate it, then composite with a statistic chosen for the index and keep the count.

# Quality-bit decoding, mask dilation and compositing for an exposure raster.
# Pinned: numpy==1.26.4, scipy==1.13.1, rasterio==1.3.10
import logging

import numpy as np
from scipy import ndimage

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

# Bit positions from the product user guide. These differ between products and
# between collections of the SAME product — read the guide for your collection.
QA_BITS = {"fill": 0, "dilated_cloud": 1, "cirrus": 2, "cloud": 3,
           "cloud_shadow": 4, "snow": 5, "water": 7}
REJECT = ("fill", "dilated_cloud", "cirrus", "cloud", "cloud_shadow", "snow")


def clear_mask(qa: np.ndarray, reject=REJECT, bits=QA_BITS,
               dilate_cells: int = 2) -> np.ndarray:
    """Boolean clear mask from a packed quality band, dilated to catch edges."""
    bad = np.zeros(qa.shape, dtype=bool)
    for name in reject:
        if name not in bits:
            raise KeyError(f"no bit position defined for {name!r}")
        # Test the individual BIT, never the packed integer against a threshold.
        bad |= ((qa >> bits[name]) & 1).astype(bool)
    if dilate_cells > 0:
        k = 2 * dilate_cells + 1
        bad = ndimage.binary_dilation(bad, structure=np.ones((k, k), dtype=bool))
    clear = ~bad
    log.info("clear mask: %.1f%% clear after %d-cell dilation",
             100.0 * clear.mean(), dilate_cells)
    return clear


def composite(scenes: list[np.ndarray], clears: list[np.ndarray],
              statistic: str = "median", min_clear: int = 4) -> dict:
    """Reduce a scene stack to one layer, keeping the clear-observation count."""
    if len(scenes) != len(clears):
        raise ValueError("one clear mask per scene is required")
    stack = np.stack([np.where(c, s.astype("float32"), np.nan)
                      for s, c in zip(scenes, clears)])
    n_clear = np.sum(~np.isnan(stack), axis=0).astype("int16")

    with np.errstate(all="ignore"):
        if statistic == "median":
            out = np.nanmedian(stack, axis=0)
        elif statistic == "mean":
            out = np.nanmean(stack, axis=0)
        elif statistic == "max":
            # Traditional for vegetation indices; biased upward by residual haze.
            out = np.nanmax(stack, axis=0)
        elif statistic == "p95":
            out = np.nanpercentile(stack, 95, axis=0)
        else:
            raise ValueError(f"unknown statistic: {statistic}")

    out = np.where(n_clear >= min_clear, out, np.nan).astype("float32")
    nulled = int((n_clear < min_clear).sum())
    log.info("composite statistic=%s scenes=%d min_clear=%d nulled_cells=%d (%.2f%%)",
             statistic, len(scenes), min_clear, nulled,
             100.0 * nulled / n_clear.size)
    return {"value": out, "n_clear": n_clear, "statistic": statistic,
            "n_scenes": len(scenes), "min_clear": min_clear}


def missingness_map(n_clear: np.ndarray, min_clear: int) -> dict:
    """Summarise where the composite is uninformed — the map of the gaps."""
    short = n_clear < min_clear
    lab, n = ndimage.label(short)
    sizes = ndimage.sum(short, lab, index=range(1, n + 1)) if n else np.array([])
    return {"nulled_fraction": float(short.mean()),
            "n_patches": int(n),
            "largest_patch_cells": int(sizes.max()) if sizes.size else 0}

Decoding individual bits rather than thresholding the packed value is the correctness point. A quality band is an integer whose bits carry independent flags, and a comparison like qa < 300 happens to exclude the right cells on some scenes and the wrong ones on others, with no error either way.

The choice of composite statistic should follow the contamination direction of the index rather than tradition. A median is robust to a minority of contaminated observations in both directions and is the right default. A maximum-value composite is appropriate only for a vegetation index where contamination is known to push values down, and it is actively harmful for a thermal product where the extreme is what the study wants but contamination pushes in the same direction.

Matching the Composite Statistic to the Contamination Direction A table-like layout relating three index families to the direction residual cloud pushes their values and the composite statistic that suits each. For a vegetation index, contamination reads low, so a maximum or high-percentile composite rejects it and a median is safe. For a thermal product, contamination reads low as well, but the study usually wants the extreme heat, so a maximum would select clear hot days correctly while a mean drags cooler in cloudy areas; a median with an explicit clear-count threshold is preferred. For an aerosol product, thin cloud reads high, so a maximum composite selects contamination and a median or low percentile is required. Residual cloud pushes different indices in different directions Vegetation index contamination reads low max composite rejects it median is safe and less sensitive to a bright edge use: median Thermal product cloud top is cold mean cools cloudy areas the study usually wants the hot extreme anyway use: median plus a percentile Aerosol product thin cloud reads high max composite selects it the traditional choice is exactly wrong here use: median Choose the statistic from the physics of the contamination, not from what the literature used and record the choice, because it changes every downstream exposure value

Validation & Edge Cases

Bit definitions that changed between collections. A reprocessing can renumber quality bits. Code that hard-codes positions will silently mask the wrong things. Load the bit map from configuration, version it, and assert the expected clear fraction is plausible.

Snow and ice flagged as cloud. In northern study areas a large share of winter observations are removed by a snow flag, which is correct for a vegetation index and wrong for a thermal study of cold exposure. Choose the reject list per application.

Water flagged but wanted. Water is often a quality flag rather than a defect. Whether to reject it depends on the index; for greenness it should be masked, for a heat study it may be informative.

Persistent-cloud regions. Where no cell reaches the threshold across a whole region, the honest outcome is a hole in the exposure raster. Filling it from a neighbouring region imports a value from a different climate.

Dilation on a coarse grid. A two-cell dilation at 30 m removes 60 m; at 1 km it removes 2 km and can erase a small study area. Scale the dilation to the cell size.

INFO exposure.mask clear mask: 61.4% clear after 2-cell dilation
INFO exposure.mask composite statistic=median scenes=18 min_clear=4 nulled_cells=1842190 (4.71%)

Dilation radius should scale with the grid, and the arithmetic is worth doing before it silently removes a small study area:

A Two-Cell Dilation Means Different Things at Different Resolutions Two grids with the same cloud and a two-cell dilation applied. On a thirty-metre grid the dilation adds a sixty-metre margin around the cloud, a small loss that removes edge contamination. On a one-kilometre grid the same two cells add a two-kilometre margin, which can remove an entire small town from the composite. The dilation should be specified in metres and converted to cells for each product rather than fixed in cells. Specify dilation in metres, then convert to cells per product 30 m grid 2 cells = 60 m margin 1 km grid town 2 cells = 2 km margin the same configuration removes a thin edge on one product and a whole town on the other

Compliance Notes

Record the product and collection identifier, the quality bit map with its source document and version, the reject list and the reason for each entry, the dilation radius in cells and in metres, the scene list with acquisition dates, the composite statistic and the reason it was chosen, the minimum clear-observation threshold, and the nulled fraction with its patch structure.

Ship the clear-observation count as a band of the composite, not as a separate file that gets lost. Every downstream consumer of the exposure raster needs it to distinguish an informed cell from an uninformed one, and separating them is how a value derived from two observations ends up being used as though it were derived from eighteen.