🔧 Fix Raster & rasterio Data cleaning Raster Analysis · step 12 of 14 Updated

Rasterio Returns the Wrong Values: NoData, Scaling and dtype Fixes

Problem statement

The elevation raster opens fine. The numbers are nonsense.

import rasterio

with rasterio.open("dem.tif") as src:
    band = src.read(1)

band.min(), band.max()      # (-32768, 1344)
band.mean()                 # -18422.7

Mean elevation of minus eighteen kilometres. The -32768 is not a measurement — it is the NoData marker, read as if it were data, and it has dragged every statistic with it.

The same class of problem has three other faces:

# temperatures that are 100× too big
band.max()                  # 3150   ← should be 31.5 °C, scale_factor is 0.01

# a difference that wraps around
band = src.read(1)          # dtype uint8
diff = band - 200           # 50 - 200 = 106, not -150

# reflectance values that are all zero
band.mean()                 # 0.0    ← integer division on an int16 band

All four are the same underlying mistake: treating the array rasterio hands you as the measurements, when it is actually the storage encoding. NoData, scale/offset and dtype are all part of the decoding, and rasterio deliberately does not apply them for you unless asked.

Quick answer

Read with masking, apply scale and offset, and convert to float before any arithmetic:

import rasterio
import numpy as np

with rasterio.open("dem.tif") as src:
    print("nodata     ", src.nodata)          # -32768.0
    print("dtype      ", src.dtypes[0])       # int16
    print("scales     ", src.scales)          # (1.0,)
    print("offsets    ", src.offsets)         # (0.0,)

    band = src.read(1, masked=True)           # a MaskedArray — nodata is masked out
    band = band.astype("float32")             # before arithmetic
    band = band * src.scales[0] + src.offsets[0]

print(band.min(), band.max(), band.mean())    # 0.4 1344.0 187.3
Symptom Cause Fix
huge negative min, absurd mean NoData read as data read(1, masked=True)
values 10×, 100× or 10000× out scale_factor / add_offset not applied multiply by src.scales, add src.offsets
subtraction wraps to a huge number unsigned dtype underflow .astype("float32") first
everything is 0 or 1 integer division cast to float before dividing
stats differ from QGIS QGIS applies nodata and scaling; rasterio does not do both explicitly
nodata is None the file never declared one find the sentinel yourself

What the array actually is

Flow from stored integers through nodata masking, scale and offset, to real-world values.
Three decoding steps between the stored bytes and a measurement. Rasterio applies none by default.

Step-by-step solution

Triage rows pairing each wrong-value symptom with its cause and fix.
Read the shape of the wrongness — it names the cause.

1. Read the metadata before the data 

with rasterio.open("dem.tif") as src:
    print(f"dtype    {src.dtypes}")
    print(f"nodata   {src.nodata}")
    print(f"scales   {src.scales}")
    print(f"offsets  {src.offsets}")
    print(f"units    {src.units}")
    print(f"crs      {src.crs}")
    print(f"tags     {src.tags(1)}")
dtype    ('int16',)
nodata   -32768.0
scales   (1.0,)
offsets  (0.0,)
units    (None,)
crs      EPSG:27700
tags     {'STATISTICS_MAXIMUM': '1344', 'STATISTICS_MINIMUM': '0.4'}

Two lines are doing the work. nodata tells you which value is not a measurement. tags(1) often carries the raster's own statistics — computed correctly by whoever wrote it, and therefore a free check on your decoding. If your min and max do not match STATISTICS_MINIMUM and STATISTICS_MAXIMUM, you have decoded it wrong.

2. Mask NoData rather than filtering it 

band = src.read(1, masked=True)          # numpy MaskedArray
type(band)                                # numpy.ma.core.MaskedArray
band.mean()                               # ignores masked cells automatically
band.count()                              # number of valid cells

A masked array is better than filtering because it keeps the raster's shape. Filtering with band[band != -32768] gives you a 1-D array, losing every spatial relationship — no neighbourhood operations, no writing it back out, no plotting.

If you prefer plain arrays, convert deliberately:

band = src.read(1).astype("float32")
band[band == src.nodata] = np.nan        # NaN needs a float dtype
print(np.nanmean(band))

Note the cast comes first. np.nan cannot be stored in an integer array — assigning it silently produces a garbage integer.

3. Apply scale and offset — they are not automatic

Many climate, satellite and sensor products store small integers with a scale factor to save space:

# stored: int16, actual: °C
real = band.astype("float32") * src.scales[0] + src.offsets[0]

src.scales and src.offsets are per-band tuples, defaulting to (1.0,) and (0.0,) when the file declares none — so applying them unconditionally is safe and costs nothing.

Watch for the same information hiding in tags rather than in the proper fields:

tags = src.tags(1)
scale = float(tags.get("scale_factor", src.scales[0]))
offset = float(tags.get("add_offset", src.offsets[0]))

NetCDF and HDF products converted to GeoTIFF frequently keep scale_factor and add_offset as tags, where rasterio will not apply them even if you ask.

4. Cast before arithmetic, every time 

band = src.read(1)                 # uint8
band.dtype                          # dtype('uint8')

band - 200                          # wraps: 50 - 200 → 106
band.astype("float32") - 200        # -150.0, correct

# integer division is the other half of the same trap
ndvi = (nir - red) / (nir + red)               # int16 → all zeros
ndvi = (nir.astype("f4") - red) / (nir.astype("f4") + red)   # correct

Unsigned integer underflow produces plausible-looking positive numbers rather than an error, which is what makes it hard to spot. The rule is simple: cast to float immediately after reading, unless you have a specific memory reason not to.

5. When nodata is None, find the sentinel yourself 

print(src.nodata)      # None

The file declares no NoData value, but the data almost certainly has one:

band = src.read(1)
import numpy as np
vals, counts = np.unique(band, return_counts=True)
for v, c in sorted(zip(vals, counts), key=lambda t: -t[1])[:5]:
    print(f"{v:>8}  {c:>10,}  {c / band.size:6.1%}")
  -32768   4,102,331   41.0%      ← 41% of the raster is one value
       0     182,004    1.8%
     412       9,221    0.1%

Forty-one percent of cells sharing one extreme value is a NoData marker, not terrain. The usual suspects are -9999, -32768, 0, 255 and 65535.

NODATA = -32768
band = np.where(band == NODATA, np.nan, band.astype("float32"))

Then fix the file so nobody has to work it out again:

with rasterio.open("dem.tif", "r+") as src:
    src.nodata = -32768

6. Preserve the decoding when you write 

profile = src.profile.copy()
profile.update(dtype="float32", nodata=np.nan, count=1)

with rasterio.open("dem_metres.tif", "w", **profile) as dst:
    dst.write(band.filled(np.nan).astype("float32"), 1)

Two mistakes to avoid here:

  • Writing a masked array directly writes the underlying data including the masked cells' raw values. Use .filled(nodata).
  • Keeping the old nodata in the profile after changing dtype. A nodata of -32768 on a float band that now uses NaN is a contradiction the next reader will trip over.

Code examples

Example 1: a decode function that handles all four cases 

import rasterio
import numpy as np

def read_decoded(path, band=1, *, sentinel=None):
    """Return real-world values as float32 with NaN for nodata."""
    with rasterio.open(path) as src:
        arr = src.read(band, masked=True).astype("float32")

        # scale/offset from the proper fields, falling back to tags
        tags = src.tags(band)
        scale = float(tags.get("scale_factor", src.scales[band - 1]))
        offset = float(tags.get("add_offset", src.offsets[band - 1]))
        if scale != 1.0 or offset != 0.0:
            arr = arr * scale + offset

        # an undeclared sentinel, if the caller found one
        if sentinel is not None:
            arr = np.ma.masked_equal(arr, sentinel * scale + offset)

        meta = {
            "crs": src.crs, "transform": src.transform,
            "nodata_declared": src.nodata, "scale": scale, "offset": offset,
            "valid_cells": int(arr.count()), "total_cells": int(arr.size),
        }
    return arr.filled(np.nan), meta

arr, meta = read_decoded("dem.tif")
print(f"{meta['valid_cells']:,}/{meta['total_cells']:,} valid "
      f"({np.nanmin(arr):.1f} to {np.nanmax(arr):.1f})")

Example 2: checking your decoding against the file's own statistics 

def verify_decoding(path, band=1):
    """The raster usually carries statistics computed by whoever wrote it."""
    with rasterio.open(path) as src:
        tags = src.tags(band)
        arr = src.read(band, masked=True).astype("float32")
        arr = arr * src.scales[band - 1] + src.offsets[band - 1]

    checks = {}
    for tag, ours in [("STATISTICS_MINIMUM", float(np.nanmin(arr))),
                      ("STATISTICS_MAXIMUM", float(np.nanmax(arr))),
                      ("STATISTICS_MEAN", float(np.nanmean(arr)))]:
        if tag in tags:
            theirs = float(tags[tag])
            checks[tag] = {"file": theirs, "ours": round(ours, 3),
                           "match": abs(theirs - ours) < max(0.01, abs(theirs) * 0.001)}
    return checks

print(verify_decoding("dem.tif"))
# {'STATISTICS_MINIMUM': {'file': 0.4, 'ours': 0.4, 'match': True},
#  'STATISTICS_MAXIMUM': {'file': 1344.0, 'ours': 1344.0, 'match': True}}

This is the single most useful check on this page. When match is False, your decoding is wrong and the file is telling you so — no external reference needed.

Example 3: sampling points without the NoData trap 

def sample_points(path, gdf, band=1):
    """Extract raster values at point locations, nodata as NaN."""
    with rasterio.open(path) as src:
        if gdf.crs != src.crs:
            gdf = gdf.to_crs(src.crs)
        coords = [(p.x, p.y) for p in gdf.geometry]
        scale, offset = src.scales[band - 1], src.offsets[band - 1]
        nodata = src.nodata

        out = []
        for value in src.sample(coords, indexes=band):
            v = float(value[0])
            out.append(np.nan if (nodata is not None and v == nodata)
                       else v * scale + offset)
    return gdf.assign(value=out)

src.sample has no masked= option, so the NoData check has to be explicit — and it is the most common place this bug survives after the main read has been fixed. See how to extract raster values at point locations.

Explanation

Two panels contrasting statistics computed with nodata included against masked correctly.
The same raster. One column is a measurement, the other is an artefact of the storage format.

Rasterio's design decision is deliberate: read() gives you what is stored, not what it means. That is the right default for a low-level library — it makes reads fast, lossless and predictable, and it means writing a file back out unchanged is trivial. But it puts the decoding on you, and the decoding has three parts.

NoData exists because a raster is a rectangle and the world is not. Everything outside the survey area, behind a cloud, or off the coast has to hold something, and that something is a sentinel value chosen to be outside the plausible range. -32768 is the minimum of int16; -9999 is a convention older than GeoTIFF. Statistics computed over those cells are not slightly wrong, they are meaningless — and because the sentinel is usually extreme, it dominates the mean.

Scale and offset exist to save space. Storing temperature to 0.01 °C precision as a float32 costs four bytes per cell; storing it as an int16 with scale_factor = 0.01 costs two, halving a 4 GB file. The decoded value is stored × scale + offset, and the metadata carries both. Whoever wrote the file expects readers to apply them — GDAL's gdalinfo, QGIS and rioxarray all do by default, which is precisely why rasterio's numbers can disagree with QGIS's for the same file.

dtype matters because NumPy arithmetic follows the array's type, not your intent. A uint8 array cannot represent -150, so 50 - 200 wraps to 106 — correct modular arithmetic, wrong answer. Integer division truncates, so an NDVI computed on int16 bands is all zeros. Neither raises.

The unifying rule: the array is an encoding, and you decode it before you reason about it. masked=True, astype("float32"), then scale and offset — in that order, every time. And when the file carries STATISTICS_* tags, use them; they are the author's own answer to whether you decoded it correctly.

Edge cases or notes

  • masked=True returns a MaskedArray, which most NumPy functions respect but some third-party code does not. .filled(np.nan) when passing it onward.
  • NaN cannot be stored in an integer array. Cast to float before assigning it, or the value becomes garbage silently.
  • Multi-band rasters have per-band nodata, scale and offset. src.nodatavals, src.scales and src.offsets are tuples — index by band - 1.
  • nodata=0 is ambiguous when zero is also a legitimate measurement, which it often is for reflectance. Prefer an out-of-range sentinel when writing.
  • NaN as a nodata value works only for float bands, and nan == nan is False, so comparisons must use np.isnan.
  • src.tags() with no argument gives dataset-level tags; src.tags(1) gives band-level. Scale factors usually live at band level.
  • rioxarray applies scale, offset and masking by default (mask_and_scale=True), which is why the same file gives different numbers through the two libraries.
  • Overviews may have their own statistics. Reading at a lower resolution with out_shape resamples, and the default resampling is nearest — fine for categorical data, wrong for continuous.

FAQ

Why does my DEM have a minimum of -32768?

That is the int16 minimum being used as a NoData marker, read as data. Use src.read(1, masked=True) and it disappears from your statistics.

Why do my values differ from QGIS?

QGIS applies NoData masking and scale/offset by default; rasterio's read() does not. Apply both explicitly and the numbers will match.

What if src.nodata is None?

The file declares none. Look for a value that occupies an implausible share of the raster with np.unique(..., return_counts=True) — 40% of cells sharing one extreme value is a sentinel.

Why did subtracting give me a huge positive number?

Unsigned integer underflow. uint8 cannot hold a negative, so it wraps. Cast to float before any arithmetic.

Should I use masked=True or NaN?

masked=True for analysis, since NumPy's masked functions ignore masked cells automatically. Convert to NaN with .filled(np.nan) when handing the array to code that does not understand masks.

Does scale_factor get applied automatically anywhere?

In rioxarray with mask_and_scale=True (the default) and in GDAL's own tools. Never in plain rasterio.read().

How do I know my decoding is right?

Compare against src.tags(1) — most rasters carry STATISTICS_MINIMUM, STATISTICS_MAXIMUM and STATISTICS_MEAN computed by the writer. If yours match, you decoded correctly.