Atmospheric Correction and Surface Reflectance

A satellite does not measure the ground. It measures light that has passed through the atmosphere twice — once on the way down, once on the way back — and has been scattered, absorbed and added to along the way. The digital numbers in a raw product encode that total signal; surface reflectance is what remains once the atmosphere’s contribution has been estimated and removed. Whether you need to care about the difference depends entirely on what you are comparing, and this topic — part of Satellite Processing Workflows & Index Pipelines — is about knowing when, and how.

The short version is to use the provider’s surface reflectance product wherever one exists, apply your own correction only when one does not, and never mix top-of-atmosphere and surface reflectance in a single time series.

---

Prerequisites

pip install "rasterio>=1.3.0" "numpy>=1.23" "pystac-client>=0.7" "xarray>=2024.1"

The numerical ground rules — casting before scaling, masking on raw values, applying per-scene offsets — are set out in converting int16 reflectance to float safely, and they apply at every step below.

---

What the Atmosphere Does to the Signal

Path radiance is the dominant effect for most land applications. It is additive, strongest at short wavelengths, and varies with aerosol load from day to day — which is why an uncorrected blue band looks hazy and why an uncorrected time series shows changes in the atmosphere as if they were changes on the ground. Attenuation is multiplicative and matters most in bands near water vapour and oxygen absorption features. The adjacency effect is local and usually second-order, visible mainly along coastlines and snow edges.

---

Step-by-Step Workflow

Step 1 — Prefer the provider’s surface reflectance product

For Sentinel-2 that is L2A; for Landsat Collection 2 it is Level-2 Surface Reflectance. Both are produced with physical radiative-transfer models, aerosol retrievals and ancillary atmospheric data that are impractical to reproduce yourself.

import pystac_client

client = pystac_client.Client.open("https://earth-search.aws.element84.com/v1")
items = client.search(
    collections=["sentinel-2-l2a"],                 # surface reflectance, not L1C
    bbox=[34.6, 0.25, 34.9, 0.55],
    datetime="2026-06-01/2026-06-30",
).item_collection()

baselines = {it.properties.get("s2:processing_baseline") for it in items}
print(f"{len(items)} L2A items, processing baselines {sorted(baselines)}")

Printing the processing baselines up front is worth the line: a search spanning a baseline change returns scenes with different offsets, and knowing that before computing anything avoids the false step described in comparing L1C and L2A Sentinel-2 products.

Step 2 — Calibrate counts where only L1 exists

Older archives, some commercial sensors and a few specialised products are available only as calibrated digital numbers or top-of-atmosphere reflectance. Converting counts to top-of-atmosphere reflectance uses gains, offsets and the sun’s elevation from the scene metadata.

import numpy as np


def toa_reflectance(dn: np.ndarray, *, mult: float, add: float,
                    sun_elevation_deg: float) -> np.ndarray:
    """Landsat-style DN to TOA reflectance with sun-angle correction."""
    rho = dn.astype("float32") * np.float32(mult) + np.float32(add)
    return rho / np.float32(np.sin(np.deg2rad(sun_elevation_deg)))

The sun-angle term matters: the same surface is less brightly lit in winter or at high latitude, and dividing by the sine of the solar elevation normalises that out. The full treatment, including where the coefficients live in the metadata, is in converting Landsat DN to TOA reflectance.

Step 3 — Apply an image-based correction if nothing better exists

Where no surface reflectance product exists and a physical model is out of reach, dark object subtraction is the standard image-based approximation. It assumes the darkest pixels in each band — deep water, dense shadow — should have near-zero reflectance, and treats whatever they actually show as path radiance to subtract.

import numpy as np


def dark_object_subtraction(band: np.ndarray, percentile: float = 0.1) -> np.ndarray:
    valid = band[np.isfinite(band)]
    dark = np.float32(np.percentile(valid, percentile))
    return np.clip(band - dark, 0, None)

It is crude — it corrects only the additive term and assumes a spatially uniform atmosphere — but it removes most of the haze from the blue and green bands and makes a single scene’s colours plausible. Its limits and a more careful variant are in applying dark object subtraction in Python.

Step 4 — Harmonise before combining sensors

Surface reflectance from two sensors is not identical even after correction, because their spectral bands differ. Combining Landsat and Sentinel-2 in one time series without adjustment introduces a step at every sensor switch.

Linear per-band coefficients, fitted on near-simultaneous acquisitions or taken from published harmonisation studies, remove most of the difference. The procedure and its limits are covered in harmonising Landsat and Sentinel-2 reflectance, and the model-feature perspective on the same problem in harmonising features across sensors for transfer.

Step 5 — Validate against stable targets

A corrected product should give plausible reflectance over targets whose reflectance is known approximately: deep clear water near zero in the near infrared, dense vegetation between roughly 0.3 and 0.5 in the near infrared, bright desert or salt flats high across the visible.

import numpy as np

def check_targets(refl: dict[str, np.ndarray], masks: dict[str, np.ndarray]) -> dict:
    """Median reflectance per band over labelled stable targets."""
    return {target: {band: float(np.nanmedian(arr[m])) for band, arr in refl.items()}
            for target, m in masks.items()}

---

Top-of-Atmosphere or Surface: Which Analyses Need Which

The organising question is whether the analysis compares values across acquisitions. Within one scene the atmosphere is roughly constant, so relative contrasts survive it and a classifier trained on that scene learns around it. Across scenes it varies, so every comparison mixes changes on the ground with changes in the sky. Change detection is the most sensitive case, and the illumination half of the same problem is covered in separating real change from illumination differences.

---

Where Correction Sits in a Pipeline

Atmospheric correction is the second step of almost every optical pipeline, immediately after reading and before anything that compares values. Its position matters for two reasons.

First, masking and correction interact. The provider’s L2A scene classification band is produced during correction and is usually the best cloud mask available, so a pipeline that uses surface reflectance gets its mask for free. A pipeline that corrects its own data must mask separately, and should do so before computing the dark-object statistics — otherwise a cloud shadow can become the “dark object” and the whole scene is over-corrected. The masking strategies are covered in masking clouds with the Sentinel-2 SCL band.

Second, every later step inherits whatever level the data is at. Indices computed on top-of-atmosphere reflectance are top-of-atmosphere indices; composites built from them mix the atmospheres of every contributing date; change maps built from those composites contain atmospheric change. None of that is visible in the output format, which is why the processing level should be recorded in the tags of every derived product, alongside the other provenance described in reading and writing GDAL tags and band descriptions.

A practical rule follows: fix the processing level at ingest, check it in every downstream step, and refuse to combine inputs at different levels. A single assertion comparing a processing_level tag across inputs prevents an entire category of quiet error.

---

Choosing a Correction When No Product Exists

When the provider offers only top-of-atmosphere data, three options remain, in increasing order of effort and accuracy.

Dark object subtraction is image-based and needs nothing beyond the scene itself. It removes most of the additive haze and suits visual products and single-scene work, but it assumes a uniform atmosphere and ignores attenuation, so its reflectance values are approximate.

Empirical line calibration uses targets of known reflectance within the scene — field-measured panels, or stable surfaces such as deep water and bright sand with well-established values — to fit a linear relationship per band. It is more accurate than dark object subtraction when good targets exist and unusable when they do not.

Physical radiative-transfer modelling, using a model such as 6S driven by aerosol and water vapour estimates, is the most accurate and by far the most involved. It is what the providers’ L2 products use, and reproducing it is rarely justified outside a research context. If an analysis genuinely needs this level of rigour and no L2 product exists, the effort is better spent finding an alternative sensor that does have one.

---

Parameter Reference

---

Cost and Storage Implications

Surface reflectance products are larger than top-of-atmosphere ones in practice, mostly because they carry extra layers — scene classification, aerosol optical thickness, water vapour — and because some providers distribute them at more than one resolution. For a pipeline that reads only four bands, that makes no difference: the extra assets are never fetched. For a pipeline that mirrors whole products locally, it can double the storage bill for layers nobody uses.

The efficient pattern is to read surface reflectance directly from the provider’s cloud-optimized assets, selecting only the bands and the classification layer the analysis needs, and to avoid mirroring at all unless the same data will be read many times. When mirroring is justified, mirror selectively — a band list, not a product — and record which bands were kept so later users do not assume the rest exist. The read-side mechanics are the same as in reading a COG over S3 without downloading, and the cost side in reducing S3 egress costs in raster pipelines.

---

Verification & Testing

import numpy as np

nir_water = np.nanmedian(nir[water_mask])
nir_forest = np.nanmedian(nir[forest_mask])
blue_water = np.nanmedian(blue[water_mask])

assert nir_water < 0.05, f"water too bright in NIR ({nir_water:.3f}) — haze or wrong level"
assert 0.2 < nir_forest < 0.6, f"forest NIR implausible ({nir_forest:.3f})"
assert blue_water < 0.08, f"blue over water {blue_water:.3f} suggests residual path radiance"

Clear water is the most useful single target. Its near-infrared reflectance is close to zero in reality, so any substantial value there is either residual path radiance, sun glint, or a scale error. Checking it takes one mask and three lines and catches most correction problems before they propagate.

---

Troubleshooting

The blue band looks milky and washed out

Path radiance has not been removed — the data is top-of-atmosphere. Use the L2 product or apply dark object subtraction.

A time series steps up at one date

Either a processing baseline changed or a sensor switched. Apply per-scene offsets and harmonise between sensors.

Water has negative reflectance

A small negative value over dark water is normal noise after correction. A large one means the correction overshot — often dark object subtraction with too high a percentile.

Indices differ from published values for the same cover

The published thresholds assume surface reflectance and you are using top-of-atmosphere, or vice versa. Match the level before comparing.

Correction brightens shadows unrealistically

The image-based method assumes a uniform atmosphere and cannot handle terrain shadow. Terrain effects need their own correction, covered in applying a terrain illumination correction.

---

Frequently Asked Questions

Q: Do I need atmospheric correction for NDVI? For a single date and a relative comparison within one scene, top-of-atmosphere NDVI is often adequate. For anything across dates, across scenes or against thresholds from the literature, use surface reflectance, because the atmosphere changes between acquisitions and biases the index differently each time.

Q: Should I run my own atmospheric correction? Rarely. Providers’ surface reflectance products for Sentinel-2 and Landsat use physical models, aerosol retrievals and ancillary data you would struggle to reproduce. Run your own only when no L2 product exists, and then prefer a documented method over an improvised one.

Q: What is the difference between L1C and L2A? L1C is top-of-atmosphere reflectance: calibrated, orthorectified, but still including the atmosphere’s contribution. L2A is bottom-of-atmosphere, or surface, reflectance with aerosols and water vapour corrected, plus a scene classification band useful for masking.

Q: Does atmospheric correction remove clouds? No. It corrects the clear-sky atmosphere; clouds and their shadows remain and must be masked separately, as covered in cloud and shadow masking strategies. Thin cirrus is partly corrected in some products and partly not, which is why masks often need dilation.

---

Related

• Converting Landsat DN to TOA Reflectance — radiometric calibration from metadata.
• Applying Dark Object Subtraction in Python — the image-based correction and its limits.
• Comparing L1C and L2A Sentinel-2 Products — what changes between levels, and baselines.
• Harmonising Landsat and Sentinel-2 Reflectance — combining sensors in one series.
• Spectral Index Calculation Pipelines — the indices this reflectance feeds.