Converting Landsat DN to TOA Reflectance

Read the per-band coefficients and sun elevation from the scene’s MTL file, rescale, then correct for sun angle:

import numpy as np

rho = dn.astype("float32") * mult + add                 # REFLECTANCE_MULT/ADD_BAND_n
rho /= np.sin(np.deg2rad(sun_elevation))                # SUN_ELEVATION from the MTL
rho[dn == 0] = np.nan                                   # 0 is fill in Landsat L1

Three numbers per band, one per scene, and a fill value — that is the whole conversion, provided they come from the right file. This page belongs to atmospheric correction and surface reflectance in Satellite Processing Workflows & Index Pipelines.

Why the Sun-Angle Term Matters

Same surface, different illumination A field is imaged in June with the sun at 65 degrees and in December with the sun at 30 degrees. The raw rescaled value in December is roughly half the June value because the surface receives about half the light per unit area. Dividing each by the sine of its sun elevation brings both to the same reflectance. June and December over one field sun at 65° raw 0.27 → /sin 65° → 0.30 sun at 30° raw 0.15 → /sin 30° → 0.30 Without the correction, the field appears to lose half its reflectance in winter.

The correction assumes flat terrain; on slopes the local illumination angle differs from the solar elevation, which is a separate correction covered in applying a terrain illumination correction.

Environment & Setup

Package Version pin Used for
rasterio >=1.3.0 Reading the DN bands
numpy >=1.23 The rescaling arithmetic
re stdlib Parsing the plain-text MTL file 
pip install "rasterio>=1.3.0" "numpy>=1.23"

Complete Working Example

import re
from pathlib import Path

import numpy as np
import rasterio


def parse_mtl(path: str) -> dict[str, float]:
    """Pull every NAME = number pair out of a Landsat MTL text file."""
    text = Path(path).read_text()
    pairs = re.findall(r"^\s*(\w+)\s*=\s*\"?([-\d.Ee+]+)\"?\s*$", text, re.M)
    return {k: float(v) for k, v in pairs}


def toa_band(dn_path: str, mtl: dict[str, float], band: int) -> np.ndarray:
    mult = mtl[f"REFLECTANCE_MULT_BAND_{band}"]
    add = mtl[f"REFLECTANCE_ADD_BAND_{band}"]
    sun_el = mtl["SUN_ELEVATION"]

    with rasterio.open(dn_path) as src:
        dn = src.read(1)
        fill = src.nodata if src.nodata is not None else 0

    rho = dn.astype("float32") * np.float32(mult) + np.float32(add)
    rho /= np.float32(np.sin(np.deg2rad(sun_el)))
    rho[dn == fill] = np.nan                              # mask on raw DN
    return rho


def write_toa(scene_dir: str, bands=(2, 3, 4, 5, 6, 7)) -> None:
    scene = Path(scene_dir)
    mtl = parse_mtl(next(scene.glob("*_MTL.txt")))
    for b in bands:
        src_path = next(scene.glob(f"*_B{b}.TIF"))
        rho = toa_band(str(src_path), mtl, b)
        with rasterio.open(src_path) as src:
            profile = src.profile | {"dtype": "float32", "nodata": np.nan,
                                     "compress": "zstd", "predictor": 3, "tiled": True}
        with rasterio.open(scene / f"TOA_B{b}.tif", "w", **profile) as dst:
            dst.write(rho, 1)
            dst.set_band_description(1, f"B{b}_TOA")
            dst.update_tags(processing_level="TOA", sun_elevation=str(mtl["SUN_ELEVATION"]))

Recording processing_level="TOA" in the output tags is the habit that keeps top-of-atmosphere data from being mistaken for surface reflectance later. The two look identical as arrays; only metadata distinguishes them, and mixing them in one time series produces exactly the kind of false change discussed on the parent page.

Where the Coefficients Come From

The three MTL fields that matter An MTL file contains many groups of metadata. Three kinds of field are needed for the conversion: REFLECTANCE_MULT_BAND_n and REFLECTANCE_ADD_BAND_n in the radiometric rescaling group, one pair per band, and SUN_ELEVATION in the image attributes group, one per scene. From the scene's *_MTL.txt GROUP = IMAGE_ATTRIBUTES SUN_ELEVATION = 58.21834 GROUP = LEVEL1_RADIOMETRIC_RESCALING REFLECTANCE_MULT_BAND_4 = 2.0000E-05 REFLECTANCE_ADD_BAND_4 = -0.100000 … one MULT and one ADD per reflective band END_GROUP Read them per scene; never copy them from one scene into another's processing.

The thermal bands use a different pair — RADIANCE_MULT and RADIANCE_ADD, followed by a brightness-temperature conversion with the K1 and K2 constants — and must not be passed through the reflectance function. Selecting coefficient names by band number, as the example does, naturally excludes them as long as the band list contains only reflective bands.

When to Stop at TOA and When Not To

For Landsat Collection 2, the Level-2 surface reflectance product exists for almost every scene and removes the atmosphere as well as the illumination. That makes the TOA conversion a fallback rather than a default: it is right for Level-1-only scenes, for quick looks, and for single-scene analyses where the classifier is trained on the same scene. For time series, change detection and anything compared with published values, the Level-2 product is the better starting point, and it has its own scale and offset described in converting int16 reflectance to float safely.

If only Level-1 exists and cross-date comparison is needed, TOA followed by dark object subtraction is the pragmatic next step, as described in applying dark object subtraction in Python.

Verification

Reading the median as a diagnostic Run it on the near-infrared band first, because vegetation makes that band's expected range tight and well known, and a problem with the coefficients or the sun term shows up there most clearly. A median near-infrared TOA reflectance between about 0.15 and 0.4 is typical over mixed land. Values in the thousands mean raw DNs were never rescaled. Values below about 0.05 over vegetated land usually mean the sun-angle term was applied twice or the wrong band's coefficients were used. Median NIR TOA reflectance over land below 0.05 sun term doubled or wrong band coefficients 0.15 – 0.40 typical mixed land conversion looks right in the thousands raw DN, never rescaled apply MULT and ADD Snow, desert and open water shift the median legitimately — check the scene before blaming the code. 
import numpy as np

rho = toa_band("LC09_..._B5.TIF", mtl, 5)                   # near infrared
valid = rho[np.isfinite(rho)]
p1, p50, p99 = np.percentile(valid, [1, 50, 99])
print(f"NIR TOA: p1 {p1:.3f}  median {p50:.3f}  p99 {p99:.3f}")
assert 0.0 < p50 < 0.6, "median NIR implausible — coefficients or fill wrong"
assert p99 < 1.2, "values above 1.2 suggest the sun-angle term was skipped or doubled"

A median near-infrared TOA reflectance between about 0.15 and 0.4 is typical over mixed land. A median in the thousands means the rescaling was skipped; one below 0.05 over a vegetated scene usually means the sun-angle correction was applied twice or the wrong band’s coefficients were used.

Common Errors

KeyError: 'REFLECTANCE_MULT_BAND_10'

A thermal band was passed to the reflectance function. Thermal bands use radiance coefficients and a temperature conversion instead.

Values around 20,000

The rescaling was not applied; these are raw DNs. Apply the multiplicative and additive coefficients, which for Collection 2 map the 16-bit range onto roughly −0.1 to 1.2 reflectance.

Winter scenes look much darker than summer

The sun-elevation correction was skipped. Divide by the sine of the scene’s own sun elevation.

A border of zero reflectance around the scene

Fill pixels were rescaled as data. Mask on the raw DN fill value before rescaling.

Two scenes on the same day disagree

They are adjacent path/row scenes with slightly different sun elevations in their metadata. Each must use its own value; the overlap will then agree to within the atmosphere’s variation.

Frequently Asked Questions

Q: Why divide by the sine of the sun elevation? Because the same surface receives less light when the sun is low. The raw rescaled value is reflectance for an overhead sun; dividing by the sine of the elevation corrects it for the actual illumination geometry, so winter and summer scenes of the same surface agree.

Q: Are the coefficients the same for every scene? For a given collection and band they are usually constant, but they are defined per scene in the metadata and have changed between collections. Always read them from the scene’s own MTL rather than hard-coding them.

Q: Is TOA reflectance good enough for my analysis? For single-scene work, often. For anything comparing dates or scenes, use Landsat Collection 2 Level-2 surface reflectance instead, which removes the atmosphere as well as the illumination difference.

Q: Does the same approach work for Sentinel-2 L1C? Sentinel-2 L1C is already delivered as scaled top-of-atmosphere reflectance, so there is no DN-to-reflectance step — only the quantification scale and, for newer baselines, an offset. The sun-angle normalisation has already been applied.

Q: Should the output be stored as float32 or packed back to integers? Packed integers are fine for storage once the conversion is done, using a fixed scale and offset recorded in the file, exactly as the provider does for its own products. Keep float32 only for intermediate files that will be consumed immediately; for an archive, the packing described in packing float rasters into int16 with scale and offset halves the size without losing anything measurable.