The Raster Data Model Explained: Bands, dtype, NoData and the Transform #

Problem statement #

A raster looks simple. It is a grid of numbers. Then this happens:

import rasterio

with rasterio.open("elevation.tif") as src:
    data = src.read(1)

print(data.mean())

-9382.117

Britain is not 9 kilometres below sea level. Or this:

print(data.max())

255

on a dataset whose documentation says the peak is 1,345 m. Or the classic: you plot a raster and a vector layer together and the raster is in the Gulf of Guinea while the vector is over Manchester.

None of these raise an exception. Every one of them is the same underlying problem — a raster carries four separate pieces of state, and reading only the array throws three of them away.

Quick answer #

A raster is an array plus three things that give the array meaning:

with rasterio.open("elevation.tif") as src:
    data = src.read(1, masked=True)      # masked=True honours NoData
    print(src.profile)

{'driver': 'GTiff', 'dtype': 'int16', 'nodata': -9999.0, 'width': 4000,
 'height': 3000, 'count': 1, 'crs': CRS.from_epsg(27700),
 'transform': Affine(25.0, 0.0, 320000.0, 0.0, -25.0, 480000.0)}

print(data.mean())      # 214.83 — sensible, because -9999 is now masked out

masked=True is the single highest-value habit in this article. It turns the NoData value from a number your statistics silently swallow into a mask NumPy respects.

Step-by-step solution #

1. The array: bands, rows and columns #

src.read() returns a 3-D array in (band, row, col) order. src.read(1) returns a 2-D array for band 1 alone.

with rasterio.open("landsat.tif") as src:
    print(src.count, src.height, src.width)   # 4 3000 4000
    all_bands = src.read()                    # (4, 3000, 4000)
    red = src.read(3)                         # (3000, 4000)

Two details trip people up constantly:

• Bands are 1-indexed. src.read(1) is the first band. src.read(0) raises IndexError: band index 0 out of range.
• Row 0 is the top. Array index [0, 0] is the north-west corner, not the origin of a coordinate system that grows upward. This is why the transform's y-scale is negative.

print(red[0, 0])          # the top-left pixel
print(red[-1, -1])        # the bottom-right pixel

2. dtype: what the numbers can hold #

The dtype is not a formality. It is a hard ceiling on what the raster can represent.

with rasterio.open("elevation.tif") as src:
    print(src.dtypes)     # ('int16',)

The failure mode is writing values a dtype cannot hold. NumPy does not warn — it wraps:

import numpy as np

ndvi = np.array([[0.82, -0.14]], dtype="float32")
print(ndvi.astype("uint8"))

[[0 0]]

Both values became zero. The NDVI is gone, no exception raised. If you must store a float range in an integer dtype, scale it deliberately and record the scale factor:

scaled = np.round(ndvi * 10_000).astype("int16")   # −1400 … 8200, recoverable

3. NoData: the value that is not a value #

A raster grid is rectangular; the world is not. Every raster of anything other than a perfect rectangle has cells that mean "no measurement here". They are marked with a sentinel value.

with rasterio.open("elevation.tif") as src:
    print(src.nodata)     # -9999.0

That sentinel is an ordinary number in the array. NumPy has no idea it is special:

plain = src.read(1)
print(plain.min(), plain.mean())      # -9999 -9382.117   ← garbage

Three ways to handle it, in order of preference:

# 1. masked array — the mask travels with the data through every operation
data = src.read(1, masked=True)
print(data.mean())                     # 214.83

# 2. convert to float and use NaN — good for arithmetic, loses integer dtype
arr = src.read(1).astype("float32")
arr[arr == src.nodata] = np.nan
print(np.nanmean(arr))                 # 214.83

# 3. boolean mask — explicit, fine for one-off statistics
raw = src.read(1)
valid = raw != src.nodata
print(raw[valid].mean())               # 214.83

4. The transform: where the grid sits #

The affine transform is six numbers mapping array indices to coordinates.

print(src.transform)

| 25.00, 0.00, 320000.00|
| 0.00,-25.00, 480000.00|
| 0.00, 0.00, 1.00|

Read as Affine(a, b, c, d, e, f):

Converting between the two spaces:

x, y = src.xy(row=100, col=250)        # pixel centre → coordinates
print(x, y)                            # 326265.0 477485.0

row, col = src.index(326265.0, 477485.0)   # coordinates → pixel
print(row, col)                            # 100 250

Use src.xy and src.index. Doing the arithmetic by hand is where the half-pixel offsets come from — c and f describe the pixel edge, while xy() returns the centre, and forgetting that shifts everything by half a cell.

5. The CRS: which earth those coordinates are on #

print(src.crs)                    # EPSG:27700
print(src.crs.is_projected)       # True

A transform without a CRS is ambiguous — 320000, 480000 could be British National Grid metres, UTM metres, or nothing at all. src.crs is None means the file never recorded it, and every downstream overlay is a guess. See CRS not found error in GeoPandas for the vector version of the same problem.

Code examples #

Example 1: a five-line audit that catches most raster surprises #

import rasterio
import numpy as np

def audit(path):
    with rasterio.open(path) as src:
        data = src.read(masked=True)
        print(f"{path}")
        print(f"  size      {src.width} × {src.height} × {src.count} band(s)")
        print(f"  dtype     {src.dtypes[0]}")
        print(f"  nodata    {src.nodata}")
        print(f"  crs       {src.crs}")
        print(f"  pixel     {src.res[0]} × {src.res[1]} {'m' if src.crs and src.crs.is_projected else 'deg'}")
        print(f"  bounds    {tuple(round(v, 1) for v in src.bounds)}")
        covered = 100 * (1 - data.mask.mean()) if np.ma.is_masked(data) else 100.0
        print(f"  valid     {covered:.1f}% of cells")
        print(f"  range     {data.min()} … {data.max()}")

audit("elevation.tif")

elevation.tif
  size      4000 × 3000 × 1 band(s)
  dtype     int16
  nodata    -9999.0
  crs       EPSG:27700
  pixel     25.0 × 25.0 m
  bounds    (320000.0, 405000.0, 420000.0, 480000.0)
  valid     61.4% of cells
  range     -3 … 1344

Two lines here do real work. valid 61.4% tells you the grid is mostly NoData — normal for a coastal tile, alarming for an inland one. range -3 … 1344 matches the documented peak, so the dtype is not clipping.

Example 2: writing a raster that other software reads correctly #

The reliable pattern is to copy the source profile and change only what you mean to change.

with rasterio.open("elevation.tif") as src:
    profile = src.profile.copy()
    data = src.read(1, masked=True)

slope_proxy = np.abs(np.gradient(data.astype("float32"))[0])

profile.update(
    dtype="float32",
    nodata=np.nan,
    count=1,
    compress="deflate",
    tiled=True,
    blockxsize=256,
    blockysize=256,
)

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

profile already holds width, height, transform and CRS. Building a profile from scratch is where people forget the transform and produce a raster that opens at the origin of the coordinate system.

compress="deflate" and tiled=True are close to free: lossless, universally readable, and tiling means later window reads touch only the blocks they need. See how to reduce GIS file size for what compression costs.

Example 3: bands that mean different things #

A four-band Landsat composite is not four rasters of the same quantity. Record what each band is:

with rasterio.open("landsat.tif") as src:
    print(src.descriptions)      # (None, None, None, None) ← nobody wrote them

with rasterio.open("landsat.tif") as src:
    profile = src.profile
    stack = src.read()

with rasterio.open("landsat_labelled.tif", "w", **profile) as dst:
    dst.write(stack)
    for i, name in enumerate(["blue", "green", "red", "nir"], start=1):
        dst.set_band_description(i, name)

Then band selection stops being a guessing game:

with rasterio.open("landsat_labelled.tif") as src:
    idx = {name: i + 1 for i, name in enumerate(src.descriptions)}
    red = src.read(idx["red"], masked=True).astype("float32")
    nir = src.read(idx["nir"], masked=True).astype("float32")

ndvi = (nir - red) / (nir + red)

Computing NDVI from bands 3 and 4 when the file is ordered nir, red, green, blue produces a plausible-looking array of entirely wrong numbers. Band descriptions are the cheapest defence against that.

Explanation #

The raster data model is old, and its awkward corners are all inherited from the constraint that a grid must be a rectangle while the thing being measured is not.

NoData exists because of that mismatch. There is no room in an array of int16 for a separate "missing" flag, so one value is sacrificed to mean it. This is exactly the design problem that gave databases NULL and floating point NaN, and the raster answer is the weakest of the three: the sentinel is indistinguishable from data unless something outside the array remembers which value is special. That "something" is a tag in the file header — one that survives a GeoTIFF round-trip but is routinely lost by a conversion to a plain array, a numpy.save, or a naïve copy. Hence masked=True: it moves the information out of the value space and into a mask that NumPy operations respect.

The transform is affine because that is all a north-up grid needs. Six coefficients express scale, rotation and translation; for the overwhelming majority of rasters the two rotation terms are zero, and the interesting content is four numbers — pixel size in x and y, and the coordinates of one corner. The negative y-scale is not a quirk but a consequence of two conventions colliding: image rows increase downward, and projected northings increase upward. Something must be negative, and by convention it is e.

dtype is a storage decision that leaks into correctness. A uint8 categorical raster with 300 classes is not a slightly lossy version of the truth; classes 256 to 300 become 0 to 44, silently colliding with real classes. Because the wrap-around is arithmetic rather than an error, the only defence is checking the range against the dtype before writing — which the audit function above does.

The consequence for daily work: a raster is not portable as an array. Handing someone a .npy file, or a NumPy array across a function boundary, discards the georeferencing, the NoData marker and often the dtype intent. Keep the array with its profile, or pass the open dataset. Everything else in raster processing — reprojection, clipping, zonal statistics — is an operation on all four parts at once, not just on the numbers.

Edge cases or notes #

• src.nodata can be None even when the raster clearly has empty regions. Some producers use an internal mask (src.read_masks(1)) or an alpha band instead. Check src.mask_flag_enums when the array looks wrong but the NoData tag is empty.
• nodata=np.nan only works for float dtypes. Setting it on an integer raster is silently ineffective, because no integer equals NaN.
• A NoData value inside the real data range is a trap. 0 as NoData on an elevation raster erases sea level. Choose a sentinel outside the plausible range.
• Rotated rasters exist. If transform.b or transform.d is non-zero, src.res and simple bounds arithmetic are misleading. Reproject to a north-up grid first.
• src.bounds describes outer edges, src.xy() returns centres. Mixing them shifts results half a pixel — the single most common cause of raster and vector not lining up.
• Scale and offset tags (src.scales, src.offsets) mean the stored integers are not the physical values. rasterio does not apply them on read. See Rasterio returns the wrong values.
• count is bands, not features. A 4,000 × 3,000 four-band uint16 raster is 96 MB in memory regardless of how much of it is NoData.
• gdalinfo -stats file.tif prints the same audit from the command line if you would rather not open Python.

Internal links #

• Introduction to Rasterio: reading raster data in Python — the practical starting point
• Vector vs raster data in Python GIS — when a grid is the wrong model entirely
• Rasterio returns the wrong values: NoData, scaling and dtype fixes — the symptoms of getting this wrong
• Raster resampling explained — what happens to those values when the grid changes
• How to reproject a raster in Python with Rasterio — changing the transform and the CRS together
• Raster and vector do not line up in Python — transform and CRS problems, diagnosed
• How to extract raster values at point locations — the transform doing its job
• Coordinate reference systems explained for Python GIS — the fourth part, in depth

FAQ #

Why is the y-scale in the transform negative? #

Because array rows increase downward while northings increase upward. The negative sign reconciles the two. A positive e describes a south-up raster, which almost always means the transform was built by hand and built wrong.

What is the difference between src.read(1) and src.read(1, masked=True)? #

read(1) gives a plain array where the NoData sentinel is an ordinary number that statistics will include. masked=True gives a numpy.ma.MaskedArray where those cells are excluded from mean, min, max and arithmetic.

Which dtype should I use? #

The smallest one that holds your full value range with the precision you need. int16 for elevation in metres, uint16 for satellite reflectance, float32 for indices and model output. float64 is almost never justified and doubles the file size.

How do I find the coordinates of a pixel? #

src.xy(row, col) returns the coordinates of the pixel centre. src.index(x, y) goes the other way. Do not do the arithmetic manually — the transform's origin is a pixel edge, not a centre.

My raster has no CRS. Can I just set one? #

Only if you know what it is. Assigning a CRS relabels the coordinates without moving them, so guessing produces data that is confidently in the wrong place. Compare against a layer whose CRS you trust before deciding.

Why did my values change when I wrote the file? #

Almost always a dtype narrower than the data. Floats written to an integer dtype are truncated; values beyond the dtype's range wrap around. Check data.min() and data.max() against the target dtype before writing.

Do I need to close the dataset? #

Use with rasterio.open(...) as src:. On write in particular, the header and any internal overviews are finalised on close — a dataset that is never closed can produce a file other software refuses to open.