Applying Dark Object Subtraction in Python

Mask clouds and shadows, take a low percentile of each band as its haze, subtract it:

import numpy as np

def dos(band: np.ndarray, clear: np.ndarray, pct: float = 0.1) -> np.ndarray:
    dark = np.float32(np.percentile(band[clear & np.isfinite(band)], pct))
    return np.clip(band - dark, 0, None)                    # remove additive path radiance

Crude, fast, needing nothing beyond the scene itself — and good enough to remove most visible haze. This page belongs to atmospheric correction and surface reflectance in Satellite Processing Workflows & Index Pipelines.

What Is Being Subtracted

Shifting the histogram back to zero Before correction the blue band's histogram starts at about 0.08, because even the darkest pixels in the scene carry the atmosphere's scattered light. Dark object subtraction takes that floor as the path radiance and shifts the whole histogram left so the darkest pixels sit near zero, as deep clear water physically should. Blue band, before and after dark floor ≈ 0.08 before after 0 0.3 The shift is largest in blue and smallest in the shortwave infrared, following the physics of scattering.

The size of the shift per band is itself a sanity check. Rayleigh and aerosol scattering fall off steeply with wavelength, so the dark value should be largest in the blue, smaller in the green and red, and close to zero in the near and shortwave infrared. A dark value that is higher in the near infrared than in the blue means the “dark object” was not dark in the near infrared — often vegetation rather than water.

Environment & Setup

Package Version pin Used for
numpy >=1.23 Percentiles and subtraction
rasterio >=1.3.0 Reading bands and writing the result
scipy >=1.11 Optional dilation of the cloud mask 
pip install "numpy>=1.23" "rasterio>=1.3.0" "scipy>=1.11"

Complete Working Example

import numpy as np
from scipy import ndimage


def clear_mask(cloud: np.ndarray, shadow: np.ndarray, *, buffer_px: int = 3) -> np.ndarray:
    """Pixels safely away from clouds and cloud shadows."""
    bad = ndimage.binary_dilation(cloud | shadow, iterations=buffer_px)
    return ~bad


def dark_values(bands: dict[str, np.ndarray], clear: np.ndarray,
                pct: float = 0.1, min_pixels: int = 5_000) -> dict[str, float]:
    out = {}
    for name, arr in bands.items():
        sel = clear & np.isfinite(arr)
        if sel.sum() < min_pixels:
            raise ValueError(f"too few clear pixels ({sel.sum()}) to estimate haze in {name}")
        out[name] = float(np.percentile(arr[sel], pct))
    return out


def apply_dos(bands: dict[str, np.ndarray], dark: dict[str, float]) -> dict[str, np.ndarray]:
    return {name: np.clip(arr - np.float32(dark[name]), 0, None).astype("float32")
            for name, arr in bands.items()}


if __name__ == "__main__":
    # bands: TOA reflectance arrays keyed by name; cloud/shadow: boolean masks
    clear = clear_mask(cloud, shadow)
    dark = dark_values(bands, clear)
    print({k: round(v, 4) for k, v in dark.items()})       # should fall with wavelength
    corrected = apply_dos(bands, dark)

Requiring a minimum number of clear pixels is the guard against the method’s worst failure: a mostly cloudy scene where the “darkest” clear pixels are just the edges of shadows the mask missed. Refusing to correct such a scene is better than silently over-darkening it, and the refusal points at a masking problem worth fixing — covered in dilating cloud masks to catch thin cirrus.

Choosing the Percentile

How the dark value moves with the percentile At the absolute minimum the dark value is dominated by a handful of noisy or defective pixels and is too low. Between roughly the 0.05th and the 1st percentile it is stable, reflecting genuinely dark surfaces. Above a few per cent it rises into ordinary land surfaces and the correction becomes too aggressive. Blue-band dark value by percentile stable plateau noisy minimum into land cover min 0.1 1 10 Percentile (log scale). Pick a value inside the plateau; 0.1 is a safe default.

The absolute minimum is a poor choice because a few defective or saturated-low pixels set it. A percentile in the plateau — 0.1 is a good default — reflects genuinely dark surfaces while ignoring outliers. Scenes with no dark surfaces at all, such as a dry landscape without water or deep shadow, have no plateau; there the method has nothing to anchor on and should not be applied.

Where the Method Breaks Down

Dark object subtraction removes an additive constant per band and nothing else. It does not correct attenuation, which scales the whole signal down in hazy conditions. It does not handle aerosol variation across the scene, which matters for large scenes near pollution sources or dust plumes. It does not handle adjacency effects near bright boundaries. And it cannot handle terrain: shadowed slopes are dark because of geometry, not atmosphere.

For those reasons its output is best treated as “haze-reduced reflectance” rather than surface reflectance. It is excellent for visual products and good for single-scene classification; for quantitative time series it narrows the gap to a physical correction without closing it, and a provider’s surface reflectance product should be preferred whenever one exists, as the comparison in comparing L1C and L2A Sentinel-2 products shows.

Verification

Dark values should fall with wavelength Estimated dark values for blue, green, red, near-infrared and shortwave infrared fall steadily from about 0.08 to near zero, matching the wavelength dependence of atmospheric scattering. A band that breaks the pattern indicates contaminated or missing dark objects. Expected pattern of per-band dark values blue green red NIR SWIR If NIR stands taller than red, the "dark object" was vegetation — do not apply the correction. 
dark = dark_values(bands, clear)
order = ["blue", "green", "red", "nir", "swir1"]
vals = [dark[b] for b in order if b in dark]
assert all(a >= b - 0.005 for a, b in zip(vals, vals[1:])), \
    f"dark values do not fall with wavelength: {dict(zip(order, vals))}"
assert dark["blue"] < 0.2, "blue dark value implausibly high — cloud leaking into the sample"

The monotonic check encodes the physics directly, and it costs nothing to run on every scene. It fails for exactly the scenes where the method is unreliable — no true dark objects, or cloud contamination in the sample — and turns a silent over-correction into an explicit error.

Common Errors

The corrected scene is too dark everywhere

A cloud shadow set the dark value. Mask and dilate clouds and shadows before estimating.

Near-infrared dark value is higher than blue

The darkest pixels are vegetation, not water or shadow. The scene lacks true dark objects; do not apply the method.

Many pixels clip to exactly zero

The percentile is too high and has reached into ordinary surfaces. Lower it into the plateau.

The correction is applied twice in a pipeline

A surface reflectance product was fed through the same step as top-of-atmosphere data. Check the processing-level tag before applying, and skip anything already corrected.

Results vary a lot between adjacent scenes

Each scene’s dark value is estimated independently, and aerosol load genuinely differs. That variation is the point, but it limits cross-scene consistency.

Frequently Asked Questions

Q: What does dark object subtraction assume? That somewhere in the scene there are pixels whose true reflectance is essentially zero in each band — deep clear water, dense shadow — so whatever they record is atmospheric path radiance. It also assumes that path radiance is the same across the scene.

Q: Why must clouds be masked first? Cloud shadows are among the darkest pixels in a scene, and they are dark for a reason unrelated to the clear-sky atmosphere. If a shadow sets the dark value, the whole scene is over-corrected and everything becomes too dark.

Q: Is it good enough for time series? It is much better than nothing and worse than a physical correction. It removes most day-to-day haze variation but not attenuation or spatial variation in aerosols, so residual differences between dates remain.

Q: Should it be applied to surface reflectance products? No. They are already corrected, and subtracting a dark value again removes real signal from dark surfaces. Apply it only to top-of-atmosphere data.

Q: Can the dark values be reused for the next scene? Only for scenes acquired on the same pass under the same conditions. Aerosol load changes from day to day, and sometimes within hours, so reusing another date’s dark values defeats the purpose of estimating them. Estimate per scene and record the values in the output’s tags for later comparison.