NDVI Values Are Above 1 or Below −1: How to Fix It
Problem statement
NDVI is (nir − red) / (nir + red). When both bands are non-negative that expression is mathematically bounded by −1 and 1. There is no vegetation, no sensor and no season that produces 1.376.
So an out-of-range NDVI is never a data story. It is always one of four bugs, and all four are in the code between reading the file and doing the division:
NDVI over pixels the scene classification calls vegetation
computed correctly 0.730
with the offset applied 1.376
Both were computed from the same two GeoTIFFs, on the same 466,954 pixels, five lines apart.
Quick answer
Check the denominator, not the numerator:
import numpy as np
total = nir + red
print(f"denominator: {np.nanmin(total):.4f} to {np.nanmax(total):.4f}")
print(f"negative denominators: {(total < 0).sum():,} of {total.size:,}")denominator: -0.1998 to 2.6600
negative denominators: 195,656 of 1,105,479
A negative denominator flips the sign of the whole fraction and removes the bound. Its cause is almost always a radiometric offset subtracted from bands that already had it subtracted.
Step-by-step solution
1. Cause one: the offset was applied twice
Sentinel-2 L2A products from processing baseline 04.00 onwards carry a −1000 DN offset, and STAC metadata declares it:
raster:bands -> {"scale": 0.0001, "offset": -0.1, "nodata": 0}
Many archives, however, harmonise their COGs — they subtract the offset at ingest so the whole archive is on one scale — and leave the metadata describing the unharmonised product. Applying the declared offset then subtracts it a second time.
Test it against water, which must be dark and positive in every band:
water median reflectance blue green red nir
DN * 1e-4 0.015 0.013 0.008 0.006
DN * 1e-4 - 0.1 -0.085 -0.087 -0.092 -0.094
Water at −0.09 reflectance in all four bands is not a dark lake. Drop the offset.
2. Cause two: unsigned integer wraparound
red = rasterio.open("B04.tif").read(1)
print(red.dtype, red.min()) # uint16 1
print(red.min() - 1000) # 64537uint16 has no negative numbers, so 1 - 1000 wraps to 64,537. After scaling that is a reflectance of 6.45, and any index built on it is nonsense — but a plausible-looking, finite nonsense.
red = red.astype("float32") # then subtractCast before any subtraction, always. This includes differences between dates, dark-object corrections and offsets.
3. Cause three: the two bands are on different scales
red = red_ds.read(1).astype("float32") * 1e-4 # reflectance
nir = nir_ds.read(1).astype("float32") # still DNA mismatch of 10,000× drives NDVI to almost exactly 1 everywhere, which looks like an implausibly healthy landscape rather than an error. If your NDVI is above 0.98 across the whole scene, check that both bands went through the same conversion.
The safest structure is one function that reads and scales, used for every band, so the two cannot drift apart.
4. Cause four: an unhandled zero denominator
Where both bands are zero — no-data fill, scene edge — the division is 0/0, which is NaN with a warning. Where only one is zero it is finite but extreme.
with np.errstate(divide="ignore", invalid="ignore"):
total = nir + red
ndvi = np.where(np.abs(total) < 1e-6, np.nan, (nir - red) / total)Do not substitute 0 for the undefined case; 0 is a real NDVI value meaning "equal red and near-infrared".
5. Add the assertion that would have caught it
finite = ndvi[np.isfinite(ndvi)]
assert finite.min() >= -1.001 and finite.max() <= 1.001, \
f"NDVI out of range: {finite.min():.3f} to {finite.max():.3f}"
Code examples
Example 1 — diagnose which of the four it is
import numpy as np
import rasterio
def diagnose_index(red_path, nir_path, scale=1e-4, offset=0.0):
"""Work out why a normalised difference left [-1, 1]."""
with rasterio.open(red_path) as ds:
red_raw = ds.read(1)
with rasterio.open(nir_path) as ds:
nir_raw = ds.read(1)
print(f" dtypes red {red_raw.dtype}, nir {nir_raw.dtype}")
print(f" shapes red {red_raw.shape}, nir {nir_raw.shape}")
print(f" raw ranges red {red_raw.min()}-{red_raw.max()}, "
f"nir {nir_raw.min()}-{nir_raw.max()}")
if red_raw.shape != nir_raw.shape:
print(" -> shapes differ: resample onto one grid first")
return
red = red_raw.astype("float32") * scale + offset
nir = nir_raw.astype("float32") * scale + offset
total = nir + red
negatives = int((total < 0).sum())
zeros = int((np.abs(total) < 1e-6).sum())
print(f" denominator {np.nanmin(total):.4f} to {np.nanmax(total):.4f}")
print(f" negative {negatives:,} ({negatives / total.size:.1%})")
print(f" ~zero {zeros:,} ({zeros / total.size:.1%})")
if negatives:
print(" -> negative denominator: the offset is probably applied twice")
if zeros:
print(" -> zero denominator: mask nodata and guard the division")
ratio = float(np.median(nir[nir > 0])) / max(float(np.median(red[red > 0])), 1e-9)
if ratio > 100 or ratio < 0.01:
print(f" -> band medians differ by {ratio:.0f}x: mixed scaling") dtypes red uint16, nir uint16
shapes red (1017, 1087), nir (1017, 1087)
raw ranges red 1-13184, nir 0-15528
denominator -0.1998 to 2.6600
negative 195,656 (17.7%)
~zero 138 (0.0%)
-> negative denominator: the offset is probably applied twice
-> zero denominator: mask nodata and guard the division
17.7% of denominators negative is decisive. A genuine dark-target artefact affects a fraction of a percent of pixels — here it is nearly a fifth of the scene, and it is exactly the fifth that is water, shadow and dark rock.
The 138 near-zero denominators are separate and real: pixels where both bands read zero at the scene edge. They need the guard whether or not the offset is right.
Example 2 — a scale-and-offset check that runs before the index
import numpy as np
import rasterio
from rasterio.enums import Resampling
# what physics says these surfaces must look like in surface reflectance
BOUNDS = {6: {"blue": (0.0, 0.10), "nir": (0.0, 0.06)}, # water
4: {"red": (0.01, 0.12), "nir": (0.15, 0.60)}} # vegetation
def verify_scaling(band_paths, scl_path, scale, offset):
with rasterio.open(scl_path) as ds:
scl = ds.read(1)
failures = []
for code, expectations in BOUNDS.items():
mask = scl == code
if mask.sum() < 500:
continue
for band, (lo, hi) in expectations.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)
if not lo <= value <= hi:
failures.append(f"class {code} {band}={value:.3f} "
f"outside {lo}-{hi}")
if failures:
raise ValueError(f"scale={scale}, offset={offset} is wrong: "
+ "; ".join(failures))
return TrueRun it once per data source, not once per scene. The answer is a property of the archive, and once you know it you can hard-code it with a comment saying how you found out.
Example 3 — a safe normalised difference for any pair of bands
import numpy as np
def normalised_difference(a, b, name="index", tolerance=1e-3):
"""(a - b) / (a + b) with the four failure modes handled."""
a = np.asarray(a, dtype="float64") # never integer arithmetic
b = np.asarray(b, dtype="float64")
if a.shape != b.shape:
raise ValueError(f"{name}: shapes {a.shape} and {b.shape} differ")
total = a + b
if np.nanmin(total) < 0:
share = float((total < 0).mean())
raise ValueError(
f"{name}: {share:.1%} of denominators are negative — "
"the inputs are not both non-negative, so the result is unbounded"
)
with np.errstate(divide="ignore", invalid="ignore"):
out = np.where(np.abs(total) < 1e-12, np.nan, (a - b) / total)
finite = out[np.isfinite(out)]
if finite.size and (finite.min() < -1 - tolerance or
finite.max() > 1 + tolerance):
raise ValueError(f"{name}: range {finite.min():.3f} to "
f"{finite.max():.3f} escapes [-1, 1]")
return out.astype("float32")Raising on a negative denominator rather than silently producing an unbounded result is the whole point. The failure is upstream and the message should say so.
Explanation
Why the bound exists and how it breaks
For non-negative a and b, |a − b| ≤ a + b, so the ratio lies in [−1, 1]. That is the entire argument, and it depends on one thing: the denominator being positive.
Make one input negative and the inequality fails. If nir = 0.261 and red = −0.044, the numerator is 0.305 and the denominator is 0.217 — a ratio of 1.41. The formula is doing exactly what it was asked to do; the inputs stopped satisfying its precondition.
That is why "clip the result to [−1, 1]" is not a fix. It hides an input error behind a plausible output.
Why the offset trap is so persistent
The metadata is not lying. scale: 0.0001, offset: -0.1 correctly describes Sentinel-2 L2A products at baseline 04.00 and later, as delivered by the mission.
What it does not describe is what a particular archive did on ingest. Harmonising — subtracting the offset so scenes either side of the 2022 baseline change are comparable — is a sensible thing for an archive to do, and it changes the stored integers without necessarily changing the item-level metadata.
The consequence is that you cannot answer "should I apply the offset?" from documentation. You answer it from a dark target, once per archive.
Why NDVI near 1 everywhere is the same class of bug
Mixed scaling is the quiet version. If nir is in DN (thousands) and red is in reflectance (hundredths), then (nir − red)/(nir + red) ≈ 1 for every pixel, because red is negligible.
The output is in range, contains no NaN, and maps as a uniformly bright, healthy landscape. Only the histogram gives it away: real NDVI has a spread, and a scene at 0.999 ± 0.001 has none.
Why the check belongs in code, not in your memory
Every one of these produces a finite array of the right dtype, shape and CRS. There is no exception, no warning, and no visual signature — an out-of-range NDVI plotted with a default colour ramp looks like a map.
Two comparisons at the end of the function catch three of the four causes, and the fourth (reading the wrong band) is caught by the class-median check in Example 2. Both are cheaper than the meeting where someone asks why the NDVI is above 1.
Edge cases or notes
- Clipping is not fixing.
np.clip(ndvi, -1, 1)converts a loud error into a silent bias. - Surface reflectance can be slightly negative over very dark targets — around −0.01 from an imperfect atmospheric model. Two thirds of a scene at −0.09 is a double correction.
0/0is NaN,x/0is ±inf. Guard onabs(total) < eps, not ontotal == 0.- Integer inputs truncate. In older NumPy
(a - b) / (a + b)on integers gave integer division; cast explicitly rather than relying on the version. - The scale cancels, the offset does not. For a ratio index you can skip the multiplicative scale entirely — but never the offset question.
- Other normalised indices have the same bound and the same failure: NDWI, NDBI, NDSI, NBR.
- A whole-scene NDVI of 0.999 is mixed scaling, not a rainforest.
nodatais 0 and 0 is a legal reflectance. Mask from the classification band.
Internal links
- Radiometric levels explained — why the declared offset can be wrong for your copy
- How to calculate NDVI from Sentinel-2 in Python — the correct order of operations
- Spectral indices explained — the other indices with the same bound
- Spectral bands explained — the reflectance values to check against
- Rasterio returns the wrong values: NoData, scaling and dtype fixes — the general form of this problem
- Satellite bands have different shapes or do not align — the shape half of the same class of bug
- How to load Sentinel-2 bands into Python as an analysis-ready array — a loader that cannot produce mixed scaling
- Cloud masking explained — the other reason an NDVI value is not trustworthy
FAQ
Why is my NDVI greater than 1?
The denominator nir + red went negative. Almost always a radiometric offset applied to bands that already had it applied, which makes one or both reflectances negative.
Why is my NDVI exactly 1 everywhere?
The two bands are on different scales — one in reflectance and one in raw digital numbers, differing by about 10,000×, which makes the smaller term negligible.
Can NDVI legitimately be negative?
Yes, down to about −0.1 over open water, and more negative over snow and cloud. Measured over water in this scene: −0.096.
Should I clip NDVI to [−1, 1]?
No. Clipping hides an input error. Fix the scaling and the range takes care of itself.
Why does 1 − 1000 give 64537?
Because the array is uint16, which has no negative values, so the subtraction wraps. Cast to float before any subtraction on a raw band.
How do I know whether to apply the STAC offset?
Compute the median reflectance over water both ways. Water must be low and positive in every band. If applying the offset makes it around −0.09, it was already applied.
Does this affect other indices?
Yes — every normalised difference index shares the bound and the failure. NDWI, NDBI, NDSI and NBR all escape [−1, 1] when a band goes negative.