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

Problem statement

The raster opens fine. The numbers are nonsense. Here is NOAA's daily sea-surface temperature analysis for 1 June 2026, as a GeoTIFF converted from NOAA's NetCDF (the setup is under The data used on this page):

import rasterio

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

print(f"min {band.min()}, max {band.max()}, mean {band.mean():.1f}")

min -999, max 3282, mean 603.3

A sea surface at minus nine hundred and ninety-nine degrees, a maximum over three thousand, and a mean that is neither. The -999 is not a measurement — it is the NoData marker for land, read as if it were data, and it has dragged every statistic with it. The 3282 is a measurement, stored a hundred times too big.

The same class of problem has two other faces:

import numpy as np

with rasterio.open("worldcover.tif") as lc:
    classes = lc.read(1)                  # uint8 land-cover codes: 10 tree cover … 80 water
print(f"{classes.dtype}: 50 - 200 = {(np.array([50], dtype=classes.dtype) - 200)[0]}")

try:
    band /= 100                           # int16, rescaled in place
except TypeError as err:
    print(f"{type(err).__name__}: {str(err)[:70]}…")

uint8: 50 - 200 = 106
UFuncTypeError: Cannot cast ufunc 'divide' output from dtype('float64') to dtype('int1…

All of them 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:

with rasterio.open("sst.tif") as src:
    print("nodata     ", src.nodata)
    print("dtype      ", src.dtypes[0])
    print("scales     ", src.scales)
    print("offsets    ", src.offsets)

    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(f"{band.min():.2f} {band.max():.2f} {band.mean():.2f}   °C")

nodata      -999.0
dtype       int16
scales      (0.009999999776482582,)
offsets     (0.0,)
-1.80 32.82 14.05   °C

What the array actually is

Metadata first, mask NoData, apply scale and offset

1. Read the metadata before the data

with rasterio.open("sst.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"bounds   {tuple(src.bounds)}")
    print(f"tags     {({k: v for k, v in src.tags(1).items() if k.startswith(('STAT', 'scale', 'add', 'units'))})}")

dtype    ('int16',)
nodata   -999.0
scales   (0.009999999776482582,)
offsets  (0.0,)
units    ('Celsius',)
crs      None
bounds   (0.0, -90.0, 360.0, 90.0)
tags     {'add_offset': '0', 'scale_factor': '0.0099999998', 'STATISTICS_MAXIMUM': '3282', 'STATISTICS_MEAN': '1404.6062837301', 'STATISTICS_MINIMUM': '-180', 'STATISTICS_STDDEV': '1188.7154610972', 'STATISTICS_VALID_PERCENT': '66.66', 'units': 'Celsius'}

Most of the lines are doing work. nodata tells you which value is not a measurement; scales says the stored integers are hundredths. Two more are worth noticing: the conversion lost the CRS, and the grid runs from 0° to 360° of longitude rather than −180° to 180° — both of which matter the moment you sample it at a place (example 3). And tags(1) carries the raster's own statistics, which are a free check on your decoding — as long as you read them in the right units (example 2).

2. Mask NoData rather than filtering it

with rasterio.open("sst.tif") as src:
    band = src.read(1, masked=True)       # numpy MaskedArray
print(type(band).__name__, f"{band.count():,} valid of {band.size:,} cells ({band.count() / band.size:.1%} ocean)")

MaskedArray 691,150 valid of 1,036,800 cells (66.7% ocean)

A masked array is better than filtering because it keeps the raster's shape. Filtering with band[band != -999] 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:

with rasterio.open("sst.tif") as src:
    plain = src.read(1).astype("float32")
    plain[plain == src.nodata] = np.nan        # NaN needs a float dtype
print(f"nanmean of the stored values: {np.nanmean(plain):.1f}")

stored = np.array([1404], dtype="int16")
try:
    stored[0] = np.nan
except ValueError as err:
    print(f"assigning NaN into int16: ValueError: {err}")

nanmean of the stored values: 1404.6
assigning NaN into int16: ValueError: cannot convert float NaN to integer

Note the cast comes first. np.nan cannot be stored in an integer array — NumPy refuses outright, which is kinder than the silent garbage some older code paths produced.

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; this one stores hundredths of a degree in an int16:

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. This file has both — GDAL copied NetCDF's scale_factor and add_offset attributes into the band's tags and set the proper fields — but a file written by other tools may carry only the tags:

with rasterio.open("sst.tif") as src:
    tags = src.tags(1)
    scale = float(tags.get("scale_factor", src.scales[0]))
    offset = float(tags.get("add_offset", src.offsets[0]))
print(f"from the tags: scale {scale}, offset {offset}; from the fields: {src.scales[0]}, {src.offsets[0]}")

from the tags: scale 0.0099999998, offset 0.0; from the fields: 0.009999999776482582, 0.0

4. Cast before arithmetic, every time

with rasterio.open("worldcover.tif") as lc:
    classes = lc.read(1)
print(f"{classes.dtype}; codes present: {sorted(np.unique(classes).tolist())}")
print("uint8 - 200:         ", (classes[:1, :3] - 200).tolist())
print("float32 - 200:       ", (classes[:1, :3].astype("float32") - 200).tolist())

with rasterio.open("sst.tif") as src:
    raw = src.read(1)
ocean = raw != -999
print("int16 / 100 gives   ", (raw[ocean][:3] / 100).tolist(), "(a float array)")
print("int16 // 100 gives  ", (raw[ocean][:3] // 100).tolist(), "(floored integers)")

uint8; codes present: [10, 30, 40, 50, 60, 80, 90]
uint8 - 200:          [[86, 86, 106]]
float32 - 200:        [[-170.0, -170.0, -150.0]]
int16 / 100 gives    [-1.8, -1.8, -1.8] (a float array)
int16 // 100 gives   [-2, -2, -2] (floored integers)

Unsigned integer underflow produces plausible-looking positive numbers rather than an error, which is what makes it hard to spot. Division is the other half of the trap: in modern NumPy / returns floats, but //, in-place /=, and any result written back into an integer array floor or refuse. 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 — carefully

The Copernicus 30 m elevation tile for 54–55°N, 3–2°W declares no NoData value:

with rasterio.open("dem_tile.tif") as dem:
    print("nodata", dem.nodata)
    z = dem.read(1)
vals, counts = np.unique(z, return_counts=True)
for v, c in sorted(zip(vals, counts), key=lambda t: -t[1])[:3]:
    print(f"{v:>8}  {c:>10,}  {c / z.size:6.1%}")

nodata None
     0.0     215,975    2.5%
    38.5      17,748    0.2%
   144.5      10,714    0.1%

The usual advice is that one value far more common than any other — and -9999, -32768, 0, 255 and 65535 are the usual suspects — is a NoData marker. Here the most common value is 0, twelve times as common as the next, and it is not a marker at all: it is the sea, and sea level is a real elevation. Masking it would punch holes in every coastline. Look at where the suspect value is before deciding:

rows, cols = np.nonzero(z == 0)
xs, ys = rasterio.transform.xy(dem.transform, rows[::500], cols[::500])
print(f"zero cells lie between {min(xs):.2f} and {max(xs):.2f}°E, {min(ys):.2f} and {max(ys):.2f}°N")

zero cells lie between -3.00 and -2.80°E, 54.00 and 54.20°N

All in the tile's south-west corner: Morecambe Bay and the Irish Sea off Barrow. When you have established that a value really is a fill, declare it, 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

with rasterio.open("sst.tif") as src:
    profile = src.profile.copy()
    decoded = src.read(1, masked=True).astype("float32") * src.scales[0] + src.offsets[0]
profile.update(dtype="float32", nodata=np.nan, count=1)

with rasterio.open("sst_celsius.tif", "w", **profile) as dst:
    dst.write(decoded.filled(np.nan), 1)
    dst.scales, dst.offsets = (1.0,), (0.0,)            # the values are now degrees
with rasterio.open("sst_celsius.tif") as check:
    back = check.read(1, masked=True)
    print(f"{check.dtypes[0]}, nodata {check.nodata}, scale {check.scales[0]}: {back.min():.2f} to {back.max():.2f} °C")

float32, nodata nan, scale 1.0: -1.80 to 32.82 °C

Three 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 -999 on a float band that now uses NaN is a contradiction the next reader will trip over.
• Keeping the old scale. The copied profile says nothing about scale, but writers that copy band metadata can carry scale_factor=0.01 onto values that are already degrees, and the next reader divides them by a hundred again. Set it back to 1 explicitly.

A decode function for all four cases

The data used on this page

• sst.tif — NOAA's OISST v2.1 daily sea-surface temperature for 1 June 2026, downloaded as NetCDF from NCEI and converted to GeoTIFF with gdal_translate: a 0.25° global grid, int16 hundredths of a degree Celsius, -999 over land. gdalinfo -stats has been run on it, as a GIS would, so it carries its own statistics.
• worldcover.tif — a 1,000 × 1,000 cell window of ESA WorldCover 2021 at 10 m, uint8 class codes, read straight from its public bucket.
• dem_tile.tif — the Copernicus 30 m elevation tile for 54–55°N, 3–2°W, which declares no NoData value.

import subprocess
from pathlib import Path

import numpy as np
import rasterio
import requests
from rasterio.windows import Window

OISST = ("https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/"
         "202606/oisst-avhrr-v02r01.20260601.nc")
Path("oisst.nc").write_bytes(requests.get(OISST, timeout=300).content)
subprocess.run(["gdal_translate", "-q", "-of", "GTiff", "netcdf:oisst.nc:sst", "sst.tif"], check=True)
subprocess.run(["gdalinfo", "-stats", "sst.tif"], check=True, capture_output=True)      # writes sst.tif.aux.xml

WC = "https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/ESA_WorldCover_10m_2021_v200_N51W003_Map.tif"
with rasterio.open(WC) as src:
    window = Window(20000, 12000, 1000, 1000)
    profile = src.profile | {"width": 1000, "height": 1000, "transform": src.window_transform(window)}
    with rasterio.open("worldcover.tif", "w", **profile) as dst:
        dst.write(src.read(1, window=window), 1)

DEM = "https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N54_00_W003_00_DEM/Copernicus_DSM_COG_10_N54_00_W003_00_DEM.tif"
Path("dem_tile.tif").write_bytes(requests.get(DEM, timeout=600).content)
print({p: f"{Path(p).stat().st_size / 1e6:.1f} MB" for p in ("sst.tif", "worldcover.tif", "dem_tile.tif")})

{'sst.tif': '2.1 MB', 'worldcover.tif': '0.1 MB', 'dem_tile.tif': '29.8 MB'}

Example 1: a decode function that handles all four cases

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("sst.tif")
print(f"{meta['valid_cells']:,}/{meta['total_cells']:,} valid ({np.nanmin(arr):.2f} to {np.nanmax(arr):.2f} °C)")

691,150/1,036,800 valid (-1.80 to 32.82 °C)

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

def verify_decoding(path, band=1):
    """Compare decoded min/max/mean with the statistics stored in the file."""
    with rasterio.open(path) as src:
        tags = src.tags(band)
        scale, offset = src.scales[band - 1], src.offsets[band - 1]
        arr = src.read(band, masked=True).astype("float32") * scale + offset

    checks = {}
    for tag, ours in [("STATISTICS_MINIMUM", float(arr.min())),
                      ("STATISTICS_MAXIMUM", float(arr.max())),
                      ("STATISTICS_MEAN", float(arr.mean()))]:
        if tag in tags:
            stored = float(tags[tag])
            checks[tag] = {"stored": stored, "stored × scale": round(stored * scale + offset, 3),
                           "ours": round(ours, 3),
                           "match": abs(stored * scale + offset - ours) < max(0.01, abs(ours) * 0.001)}
    return checks


for tag, c in verify_decoding("sst.tif").items():
    print(f"{tag:20} {c}")

STATISTICS_MINIMUM   {'stored': -180.0, 'stored × scale': -1.8, 'ours': -1.8, 'match': True}
STATISTICS_MAXIMUM   {'stored': 3282.0, 'stored × scale': 32.82, 'ours': 32.82, 'match': True}
STATISTICS_MEAN      {'stored': 1404.6062837301, 'stored × scale': 14.046, 'ours': 14.046, 'match': True}

This is the single most useful check on this page — with one catch this file demonstrates. The statistics GDAL stores are computed on the stored values, NoData excluded but scale not applied: a minimum of -180 means −1.80 °C. Compare your decoded numbers with the stored statistics directly and every row fails; decode the stored statistics the same way and they match. When a check like this does fail, your decoding is wrong and the file is telling you so — no external reference needed.

Example 3: sampling points without the NoData trap

import geopandas as gpd
from shapely import Point


def sample_points(path, gdf, band=1, *, assume_crs=None):
    """Extract decoded raster values at point locations, nodata as NaN."""
    with rasterio.open(path) as src:
        crs = src.crs or assume_crs
        if crs is None:
            raise ValueError("the raster has no CRS — pass assume_crs from the product documentation")
        pts = gdf.to_crs(crs)
        xs = pts.geometry.x.to_numpy()
        if src.bounds.right > 180:                     # a 0–360 longitude grid
            xs = np.where(xs < 0, xs + 360, xs)
        scale, offset, nodata = src.scales[band - 1], src.offsets[band - 1], src.nodata
        out = [np.nan if (nodata is not None and v[0] == nodata) else v[0] * scale + offset
               for v in src.sample(zip(xs, pts.geometry.y), indexes=band)]
    return gdf.assign(sst=out)


places = gpd.GeoDataFrame({"name": ["North Sea off Whitby", "Bay of Biscay", "Birmingham (land)"]}, crs=4326,
                          geometry=[Point(0.2, 54.6), Point(-5.0, 45.0), Point(-1.9, 52.5)])
print(sample_points("sst.tif", places, assume_crs="EPSG:4326").drop(columns="geometry").round(2).to_string(index=False))

                name   sst
North Sea off Whitby 12.88
       Bay of Biscay 18.51
   Birmingham (land)   NaN

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. The two other lines this file forced are just as easy to miss: the CRS has to come from the product documentation (OISST is a regular latitude-longitude grid, so EPSG:4326), and the Bay of Biscay's longitude of −5° has to become 355° to land on a grid that runs from 0 to 360. See how to extract raster values at point locations.

read() returns what is stored, not what it means

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 $\text{stored} \times \text{scale} + \text{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.

Masked arrays and NaN in integers

• 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; current NumPy raises ValueError: cannot convert float NaN to integer (step 2), and code paths that bypass the check store garbage.

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

• Tools disagree about the default. rioxarray.open_rasterio leaves values raw unless you pass mask_and_scale=True (the default was False in rioxarray 0.23), while xarray.open_dataset on the original NetCDF masks and scales by default — so the same data can give different numbers depending on which door it came through.

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

Internal links

• Introduction to Rasterio: reading raster data in Python — the basics this page assumes
• How to extract raster values at point locations with Rasterio — where the sampling trap in example 3 lives
• RasterioIOError: not recognized as a supported file format — when the file will not open at all
• How to batch process rasters with Rasterio in Python — applying this decode across a folder
• Vector vs raster data in Python GIS — what a band actually is
• How to convert a raster to a vector in Python — where NoData decides the polygon boundaries
• How to handle missing and null values in spatial datasets — the vector-side equivalent
• Null, empty, missing and invalid: four kinds of broken geometry — the same "absent means several things" problem
• A hazard raster's nodata is read as zero — the same bug on data that is mostly nodata
• ppb and µg/m³ mixed up — a unit error that passes every range check
• LST values on the wrong scale — scale, offset and fill in surface temperature

FAQ

Why does my DEM have a minimum of -32768?

That is the int16 minimum being used as a NoData marker, read as data — the same thing the -999 in this page's sea-surface temperatures does. Use src.read(1, masked=True) and it disappears from your statistics. If the file declares no NoData value, find the sentinel as in step 5, and check where it is before masking it.

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.