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.

A doubled offset, integer wraparound and other causes

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

dn = rasterio.open("B04.tif").read(1)
print(dn.dtype, dn.min())
print(dn.min() - 1000)          # wraps: uint16 has no negative numbers

uint16 1
64537

uint16 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.

dn = dn.astype("float32")       # then subtract

Cast 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 DN

A 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}"

Diagnosing which cause it is

The data used on this page

One Sentinel-2 L2A scene from Earth Search, S2A_30UVD_20250922_0_L2A, cut to 4.16–4.00° W, 53.03–53.12° N around Snowdon and written as B02.tif, B03.tif, B04.tif, B08.tif and SCL.tif on the 10 m grid (the classification resampled with nearest neighbour). red and nir are those bands read the way a pipeline that trusts the declared STAC metadata reads them: DN × 0.0001 − 0.1.

import os

import numpy as np
import pystac_client
import rasterio
from rasterio.enums import Resampling
from rasterio.warp import transform_bounds

AOI = (-4.16, 53.03, -4.00, 53.12)
item = pystac_client.Client.open("https://earth-search.aws.element84.com/v1").get_collection(
    "sentinel-2-l2a").get_item("S2A_30UVD_20250922_0_L2A")
files = {"blue": "B02.tif", "green": "B03.tif", "red": "B04.tif", "nir": "B08.tif", "scl": "SCL.tif"}
if not all(os.path.exists(f) for f in files.values()):
    with rasterio.open(item.assets["red"].href) as ds:
        win = rasterio.windows.from_bounds(*transform_bounds("EPSG:4326", ds.crs, *AOI),
                                           transform=ds.transform).round_offsets().round_lengths()
        shape, crs, transform = (int(win.height), int(win.width)), ds.crs, ds.window_transform(win)
    for name, path in files.items():
        with rasterio.open(item.assets[name].href) as ds:
            w = rasterio.windows.from_bounds(*transform_bounds("EPSG:4326", ds.crs, *AOI),
                                             transform=ds.transform).round_offsets().round_lengths()
            data = ds.read(1, window=w, out_shape=shape, resampling=Resampling.nearest)
        with rasterio.open(path, "w", driver="GTiff", width=shape[1], height=shape[0], count=1,
                           dtype=data.dtype, crs=crs, transform=transform) as out:
            out.write(data, 1)

declared = item.assets["red"].extra_fields["raster:bands"][0]
red = rasterio.open("B04.tif").read(1).astype("float32") * declared["scale"] + declared["offset"]
nir = rasterio.open("B08.tif").read(1).astype("float32") * declared["scale"] + declared["offset"]
print(f"{red.shape}; declared scale {declared['scale']}, offset {declared['offset']}")

(1017, 1087); declared scale 0.0001, offset -0.1

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")


diagnose_index("B04.tif", "B08.tif", offset=-0.1)

  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 True


paths = {"blue": "B02.tif", "red": "B04.tif", "nir": "B08.tif"}
for offset in (0.0, -0.1):
    try:
        print(f"offset {offset}: passes", verify_scaling(paths, "SCL.tif", 1e-4, offset))
    except ValueError as exc:
        print(f"offset {offset}: {exc}")

offset 0.0: passes True
offset -0.1: scale=0.0001, offset=-0.1 is wrong: class 6 blue=-0.085 outside 0.0-0.1; class 6 nir=-0.094 outside 0.0-0.06; class 4 red=-0.044 outside 0.01-0.12

Run 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")


dn_red = rasterio.open("B04.tif").read(1).astype("float32")
dn_nir = rasterio.open("B08.tif").read(1).astype("float32")
good = normalised_difference(dn_nir * 1e-4, dn_red * 1e-4, "NDVI")
print(f"NDVI {np.nanmin(good):.3f} to {np.nanmax(good):.3f}")
try:
    normalised_difference(nir, red, "NDVI with the declared offset")
except ValueError as exc:
    print(exc)

NDVI -1.000 to 0.942
NDVI with the declared offset: 17.7% of denominators are negative — the inputs are not both non-negative, so the result is unbounded

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.

Example 4 — the numbers at the top of this page

scl = rasterio.open("SCL.tif").read(1)
veg = scl == 4
for label, a, b in (("computed correctly", dn_nir * 1e-4, dn_red * 1e-4),
                    ("with the offset applied", nir, red)):
    with np.errstate(divide="ignore", invalid="ignore"):
        v = ((a - b) / (a + b))[veg]
    print(f"median NDVI over {int(veg.sum()):,} vegetation pixels, {label}: {np.nanmedian(v):.3f}")

with np.errstate(divide="ignore", invalid="ignore"):
    print(f"median NDVI over water, computed correctly: "
          f"{np.nanmedian(((dn_nir - dn_red) / (dn_nir + dn_red))[scl == 6]):.3f}")
print("water median reflectance     blue    green     red      nir")
water = scl == 6
med = {b: np.median(rasterio.open(f).read(1)[water].astype("float32")) * 1e-4
       for b, f in (("blue", "B02.tif"), ("green", "B03.tif"), ("red", "B04.tif"), ("nir", "B08.tif"))}
print("  DN * 1e-4               " + "  ".join(f"{v:6.3f}" for v in med.values()))
print("  DN * 1e-4 - 0.1         " + "  ".join(f"{v - 0.1:6.3f}" for v in med.values()))

median NDVI over 466,954 vegetation pixels, computed correctly: 0.730
median NDVI over 466,954 vegetation pixels, with the offset applied: 1.376
median NDVI over water, computed correctly: -0.096
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

Why the bound exists and how it breaks

Why the bound exists and how it breaks

For non-negative $a$ and $b$, $|a - b| \le 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 $(\text{nir} - \text{red})/(\text{nir} + \text{red}) \approx 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.

Clipping, negative reflectance and infinities

• 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/0 is NaN, $x/0$ is ±inf. Guard on abs(total) < eps, not on total == 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.
• nodata is 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.