Radiometric Levels Explained: DN, TOA and Surface Reflectance
Problem statement
Two people download the same Sentinel-2 scene, compute NDVI over the same field, and get 0.73 and 1.38. Neither made an arithmetic mistake. They disagreed about what the integers in the file mean.
Satellite imagery comes in processing levels, and each level stores a different physical quantity in the same uint16 array:
- DN — a raw detector count. Comparable to nothing.
- L1C / TOA reflectance — what the sensor saw, atmosphere included.
- L2A / surface reflectance — an estimate of what the ground would look like without an atmosphere.
The dangerous part is that the conversion between the stored integer and the physical value is described in metadata, the metadata is sometimes stale, and applying a correction twice produces numbers that are wrong in a way no exception will catch.
Quick answer
Never trust a declared scale and offset without checking it against a target whose reflectance you already know:
import numpy as np
import rasterio
from rasterio.enums import Resampling
with rasterio.open("clear_scl.tif") as ds:
scl = ds.read(1) # scene classification, 20 m
with rasterio.open("clear_blue.tif") as ds:
blue = ds.read(1, out_shape=scl.shape, resampling=Resampling.nearest)
water = blue[scl == 6].astype("float32") # class 6 = water
print("water DN ", np.median(water))
print("DN * 1e-4 ", np.median(water) * 1e-4) # 0.015 plausible
print("DN * 1e-4 - 0.1", np.median(water) * 1e-4 - 0.1) # -0.085 impossibleWater is dark but not negative. The second conversion is wrong for this product even though the STAC metadata declares exactly that scale and offset.
Step-by-step solution
1. Identify the level you actually have
The product name is the fastest signal. L1C is top-of-atmosphere; L2A is surface reflectance; a Landsat _SR_ product is surface reflectance and _TOA_ is not.
item = next(client.search(ids=["S2A_30UVD_20250922_0_L2A"],
collections=["sentinel-2-l2a"]).items())
print(item.properties["s2:processing_baseline"]) # 05.11
print(item.assets["red"].extra_fields["raster:bands"][0]){'nodata': 0, 'data_type': 'uint16', 'spatial_resolution': 10,
'scale': 0.0001, 'offset': -0.1}
Baseline 05.11 is after the January 2022 change that introduced a −1000 DN offset, so the declared offset: -0.1 is the correct description of raw baseline-04.00+ data.
2. Test the conversion against known targets
Metadata describes the product family. It does not always describe the copy you have, because archives are often harmonised — reprocessed so that scenes either side of a baseline change line up, which means the offset has already been subtracted from the stored integers.
Median reflectance over the Snowdonia window under each interpretation:
blue green red nir
vegetation no offset 0.039 0.067 0.056 0.361
+offset -0.061 -0.033 -0.044 0.261
water no offset 0.015 0.013 0.008 0.006
+offset -0.085 -0.087 -0.092 -0.094
cloud no offset 0.527 0.507 0.504 0.599
+offset 0.427 0.407 0.404 0.499
The offset interpretation gives water a reflectance of −0.09 in every band and vegetation a negative blue reflectance. Neither is physically possible. The no-offset interpretation gives a textbook water spectrum, a textbook vegetation spectrum and a bright, flat cloud.
The offset is already applied in this archive's COGs. Applying it again subtracts 0.1 from every band a second time.
3. Watch what a double correction does to an index
NDVI over pixels the scene classification calls vegetation
no offset 0.730
with offset 1.376
NDVI is bounded by −1 and 1 whenever both inputs are positive. A value of 1.376 is not a strange vegetation reading; it is proof that the denominator went negative. See NDVI values are above 1 or below −1.
4. Cast before you subtract
red = rasterio.open("clear_red.tif").read(1)
print(red.dtype, red.min()) # uint16 1
print(red.min() - 1000) # 64537Unsigned integers wrap. A dark pixel with DN 1 becomes 64,537, which after scaling is a reflectance of 6.45. This single line has produced more "impossible reflectance" bug reports than atmospheric correction ever has.
red = red.astype("float32") # then, and only then, apply the offset
5. Record the level in the output
Every derived product should carry the level it was computed from. A time series that silently mixes L1C and L2A has a step change in it that looks like a real event.
Code examples
Example 1 — a conversion check you can run on any product
import numpy as np
import rasterio
from rasterio.enums import Resampling
# expected surface reflectance ranges, from physics rather than from a header
EXPECTED = {
"water": {"blue": (0.00, 0.10), "nir": (0.00, 0.06)},
"vegetation": {"red": (0.01, 0.12), "nir": (0.15, 0.60)},
}
SCL_CODE = {"water": 6, "vegetation": 4}
def check_scaling(band_paths, scl_path, scale, offset):
"""Does this scale/offset put known surfaces where physics says they belong?"""
with rasterio.open(scl_path) as ds:
scl = ds.read(1)
problems = []
for surface, bands in EXPECTED.items():
mask = scl == SCL_CODE[surface]
if mask.sum() < 500:
continue
for band, (lo, hi) in bands.items():
with rasterio.open(band_paths[band]) as ds:
dn = ds.read(1, out_shape=scl.shape, resampling=Resampling.nearest)
value = float(np.median(dn[mask].astype("float32")) * scale + offset)
flag = "ok" if lo <= value <= hi else "OUT OF RANGE"
print(f" {surface:11} {band:5} {value:7.3f} expected {lo}-{hi} {flag}")
if not lo <= value <= hi:
problems.append(f"{surface}/{band} = {value:.3f}")
if problems:
raise ValueError(
f"scale={scale} offset={offset} is wrong for this product: "
+ "; ".join(problems)
)
print(f" scale={scale} offset={offset} is consistent with known surfaces") water blue 0.015 expected 0.0-0.1 ok
water nir 0.006 expected 0.0-0.06 ok
vegetation red 0.056 expected 0.01-0.12 ok
vegetation nir 0.361 expected 0.15-0.6 ok
scale=0.0001 offset=0 is consistent with known surfaces
And with the declared offset:
water blue -0.085 expected 0.0-0.1 OUT OF RANGE
water nir -0.094 expected 0.0-0.06 OUT OF RANGE
vegetation red -0.044 expected 0.01-0.12 OUT OF RANGE
vegetation nir 0.261 expected 0.15-0.6 ok
ValueError: scale=0.0001 offset=-0.1 is wrong for this product:
water/blue = -0.085; water/nir = -0.094; vegetation/red = -0.044
Three of four checks fail, and the one that passes is the one with the widest tolerance. That is the shape of a systematic offset error: it moves everything by the same amount, so the tightest constraint fails first.
Example 2 — reading a band at a stated level, safely
import numpy as np
import rasterio
def read_level(path, level, scale=1e-4, offset=0.0, nodata=0):
"""Read a band and label it with the processing level it came from."""
if level not in {"dn", "toa", "surface"}:
raise ValueError(f"unknown level {level!r}")
with rasterio.open(path) as ds:
dn = ds.read(1)
tags = ds.tags()
values = dn.astype("float32") # cast first: uint16 arithmetic wraps
values[dn == nodata] = np.nan
if level == "dn":
return values, {"level": "dn", "units": "detector counts"}
reflectance = values * scale + offset
return reflectance, {
"level": level,
"units": "reflectance",
"scale": scale,
"offset": offset,
"source": tags.get("TIFFTAG_IMAGEDESCRIPTION", path),
}Carrying the metadata dictionary alongside the array is the cheapest way to stop a downstream function guessing. When a composite is built from ten scenes, one function can assert that all ten metadata dictionaries agree.
Example 3 — detecting a level change inside a time series
import numpy as np
def find_level_breaks(dates, medians, threshold=0.05):
"""A step in the series baseline that no season would produce."""
medians = np.asarray(medians, dtype=float)
steps = np.diff(medians)
breaks = []
for i, step in enumerate(steps):
if abs(step) > threshold:
breaks.append({"between": (dates[i], dates[i + 1]),
"step": round(float(step), 3)})
for b in breaks:
print(f" {b['between'][0]} -> {b['between'][1]}: "
f"median moves {b['step']:+.3f}")
if not breaks:
print(" no step changes above threshold")
return breaksA processing-level change shows up as a step of roughly 0.1 in reflectance — exactly the size of the offset — applied to every pixel on the same day. Vegetation phenology does not do that. If you see a uniform step on a known reprocessing date, the level changed, not the ground.
Explanation
What each level actually removes
DN is proportional to photons hitting a detector, scaled by whatever gain that detector was set to. Two DN values from different sensors, or the same sensor with different gain settings, are not comparable.
TOA reflectance divides the measured radiance by the incoming solar irradiance, corrected for sun angle and earth-sun distance. It removes "the sun was lower in December" and produces a ratio that is comparable across dates and sensors — but it still contains the atmosphere.
Surface reflectance additionally models and removes atmospheric scattering and absorption. Over a hazy scene the atmosphere contributes a large fraction of the blue signal, so the correction is biggest in blue and smallest in the near-infrared.
Each rung removes one source of variation that is not the ground. Each rung is also an additional model that can be wrong — which is why L2A products carry quality bands and why some analyses deliberately stay at L1C.
Why the baseline change matters
Before January 2022, Sentinel-2 L2A stored reflectance as DN / 10000, and dark pixels clipped at zero. Since processing baseline 04.00, an offset of −1000 DN is added before storage, so genuinely dark and slightly negative reflectances survive rather than being clipped.
That means the same physical scene stored under the two baselines differs by 1000 DN. A time series that crosses the change without handling it has a 0.1 reflectance step in it.
Archives solve this in one of two ways: keep the raw DN and declare the offset in metadata, or harmonise — subtract the offset at ingest so every scene in the archive is on the old scale. Both are defensible. The trap is metadata that describes the first while the files were produced by the second, which is what the measurements above found.
Why a physical check beats reading the documentation
Documentation describes intent; the file in front of you is the fact. A check against known surfaces takes four lines, runs in a second, and fails loudly on exactly the error class that is otherwise silent.
Water is the best target because it is dark in every band and its reflectance falls with wavelength, so a constant offset shows up as a sign error rather than as a plausible-looking number. Deep, calm water is best; shallow or turbid water is brighter and less diagnostic.
Cloud is the second-best target because it is bright and nearly flat across the visible bands. If cloud comes out above 1.0, the scale is too large.
Why this is not pedantry
The measurement above is the difference between an NDVI of 0.73 and an NDVI of 1.38 on the same pixels. Downstream, that is the difference between "healthy vegetation" and a number outside the range of the index, which will be silently clipped, plotted, averaged into a regional statistic, and reported.
Nothing raises an exception. The array has the right shape, the right dtype and the right CRS, and every value in it is wrong by the same amount.
Edge cases or notes
- Harmonisation is per-archive, not per-satellite. The same scene from two providers can need different handling. Check each source once and record the answer.
offsetin STACraster:bandsis in reflectance units (−0.1), not DN (−1000). Mixing the two gives an error of 1000×.- Never subtract on an unsigned array. Cast to float first, or
1 - 1000becomes 64,537. nodatais 0. Mask it before computing statistics, or the fill area drags every median down.- Slightly negative surface reflectance is legitimate over very dark targets — around −0.01, from an imperfect atmospheric model. Around −0.09 across all bands is a double correction.
- L1C and L2A are not interchangeable for indices. NDVI from TOA is systematically lower than from surface reflectance because atmospheric scattering brightens red more than near-infrared.
- Thermal bands are not reflectance at all. They convert to brightness temperature with a completely different formula.
Internal links
- Spectral bands explained: what satellite imagery actually measures — the three things a band is
- NDVI values are above 1 or below −1: how to fix it — the most common symptom of this problem
- How to calculate NDVI from Sentinel-2 in Python — the conversion done correctly
- How to load Sentinel-2 bands into Python as an analysis-ready array — where the scaling belongs in a pipeline
- Rasterio returns the wrong values: NoData, scaling and dtype fixes — the same class of problem in general raster work
- STAC explained: how satellite imagery catalogues work — where scale and offset metadata come from
- How to extract a vegetation index time series for a polygon — where a level change shows up as a fake event
- Cloud masking explained — the other band you need before trusting a pixel
FAQ
What is the difference between L1C and L2A?
L1C is top-of-atmosphere reflectance — what the sensor saw, including haze and scattering. L2A is surface reflectance, with an atmospheric model removed. Use L2A for anything comparing dates or places.
What is a digital number?
The raw integer a detector produces, proportional to photons collected and to the sensor gain. It is not comparable between sensors, or sometimes even between scenes.
Should I apply the offset in the STAC metadata?
Check first. Compute median reflectance over water with and without it; water should be low and positive in every band. On the Earth Search Sentinel-2 L2A COGs measured here, applying it gives −0.09 and is wrong.
Why is my NDVI above 1?
Almost always a double-applied offset making the denominator negative, or unsigned-integer arithmetic wrapping. See NDVI values are above 1 or below −1.
Can surface reflectance be negative?
Slightly, over very dark targets, because atmospheric correction is a model. Values around −0.01 are normal; −0.09 across every band is a processing error.
Does the processing level matter for NDVI?
Yes. Atmospheric scattering brightens red more than near-infrared, so TOA NDVI is systematically lower than surface NDVI. A time series must not mix them.
How do I make a time series that crosses the 2022 baseline change safe?
Use a harmonised archive, or apply the offset only to scenes with processing baseline 04.00 and later. Then plot the scene medians and check for a 0.1 step on the changeover date.