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

Problem statement

A raster looks simple. It is a grid of numbers. Here is the Environment Agency's lidar terrain model of Robin Hood's Bay on the North Yorkshire coast, and its mean height:

import rasterio

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

print(data.mean(), data.min())

-inf -3.4028235e+38

The Yorkshire coast is not infinitely deep. Or this — the same file saved with a uint8 profile to make it smaller. Rasterio refuses the first attempt, because the NoData value does not fit in eight bits, so NoData becomes 0 as well:

with rasterio.open("elevation.tif") as src:
    try:
        rasterio.open("elevation_uint8.tif", "w", **(src.profile | {"dtype": "uint8"}))
    except ValueError as err:
        print(f"ValueError: {err}")
    small = src.profile | {"dtype": "uint8", "nodata": 0}
    with rasterio.open("elevation_uint8.tif", "w", **small) as dst:
        dst.write(src.read(1), 1)
with rasterio.open("elevation_uint8.tif") as src:
    print("max:", src.read(1).max())

ValueError: Given nodata value, -3.4028234663852886e+38, is beyond the valid range of its data type, uint8.
max: 255

The terrain here tops out at 266 m on the moor above Ravenscar, so a maximum of 255 looks like clipping. It is worse than that, as step 2 shows. Or the classic: an array saved without its georeferencing opens with no CRS and a grid whose corner sits at (0, 0) — which, drawn on a web map, is the Gulf of Guinea:

with rasterio.open("elevation.tif") as src:
    bare = {k: v for k, v in src.profile.items() if k not in ("crs", "transform")}
    with rasterio.open("elevation_bare.tif", "w", **bare) as dst:
        dst.write(src.read(1), 1)
with rasterio.open("elevation_bare.tif") as src:
    print(src.crs, src.bounds)

None BoundingBox(left=0.0, bottom=1600.0, right=2000.0, top=0.0)

None of these raises 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({k: src.profile[k] for k in ("driver", "dtype", "nodata", "width", "height", "count")})
    print(src.crs, src.transform[:6])
print(f"mean {data.mean():.2f} m")

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.4028234663852886e+38, 'width': 2000, 'height': 1600, 'count': 1}
EPSG:27700 (5.0, 0.0, 490000.0, 0.0, -5.0, 508000.0)
mean 91.34 m

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.

The array, the dtype, NoData and the transform

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. A four-band Sentinel-2 stack of the same bay:

with rasterio.open("sentinel.tif") as src:
    print(src.count, src.height, src.width)
    all_bands = src.read()
    red = src.read(3)
    print(all_bands.shape, red.shape)
    try:
        src.read(0)
    except IndexError as err:
        print(f"IndexError: {err}")

4 814 1011
(4, 814, 1011) (814, 1011)
IndexError: band index 0 out of range (not in (1, 2, 3, 4))

Two details trip people up constantly:

• Bands are 1-indexed. src.read(1) is the first band, and src.read(0) raises.
• 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.

with rasterio.open("sentinel.tif") as src:

    for label, (row, col) in [("top-left", (0, 0)), ("bottom-right", (src.height - 1, src.width - 1))]:
        x, y = src.xy(row, col)
        print(f"{label:12} row {row:>3}, col {col:>4} → pixel centre {x:.0f} E, {y:.0f} N")

top-left     row   0, col    0 → pixel centre 654715 E, 6037405 N
bottom-right row 813, col 1010 → pixel centre 664815 E, 6029275 N

The northing falls as the row number rises.

2. dtype: what the numbers can hold

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

The failure mode is writing values a dtype cannot hold, and neither NumPy nor Rasterio stops you. Back to the uint8 copy of the terrain:

import numpy as np

with rasterio.open("elevation.tif") as src:
    heights = src.read(1, masked=True)
with rasterio.open("elevation_uint8.tif") as src:
    stored = src.read(1)

for label, pick in [("highest point", heights.argmax()), ("lowest point", heights.argmin())]:
    r, c = np.unravel_index(pick, heights.shape)
    print(f"{label:14} {heights[r, c]:7.2f} m  stored as {stored[r, c]}")

land = ~heights.mask
off = np.abs(stored[land].astype("float32") - heights.compressed()) >= 1
print(f"surveyed cells stored more than 1 m from their height: {off.mean():.1%}")
print(f"sea cells (NoData) stored as 0: {(stored[~land] == 0).mean():.0%}")

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

highest point   265.99 m  stored as 9
lowest point     -3.10 m  stored as 253
surveyed cells stored more than 1 m from their height: 28.1%
sea cells (NoData) stored as 0: 100%
NDVI as uint8: [[0 0]]

This is not clipping at 255. The float heights were cast to 8 bits, and anything outside 0–255 wrapped around: the moor top came out as 9 m, and the foreshore below Ordnance Datum — a quarter of the surveyed cells sit around −2 m at the time of the flight — came out as 253 to 255, the highest values in the file. More than a quarter of the surveyed cells are now wrong by over a metre, and the sea, NoData at −3.4 × 10³⁸, became 0. The NDVI fared no better: both values truncated towards zero. Converting an out-of-range float to an integer is not even defined behaviour in C, so the exact numbers can differ by platform. What does not differ is that nothing warns you. If you must store a float range in an integer dtype, scale it deliberately and record the scale factor:

scaled = np.round(heights.compressed() * 10).astype("int16")    # decimetres, recoverable
print(scaled.min(), scaled.max())

-31 2660

3. NoData: the value that is not a value

A raster grid is rectangular; the world is not. Much of this rectangle is North Sea, beyond the edge of the lidar survey. Those cells are marked with a sentinel value:

with rasterio.open("elevation.tif") as src:
    print(src.nodata, "=", np.finfo("float32").min)
    plain = src.read(1)
print(f"{(plain == src.nodata).mean():.1%} of cells are NoData")

-3.4028234663852886e+38 = -3.4028235e+38
39.4% of cells are NoData

The Environment Agency uses the most negative number a float32 can hold, and it is an ordinary number in the array. NumPy has no idea it is special, which is how the mean became -inf: two cells in five at −3.4 × 10³⁸ overflow the sum. Three ways to handle it, in order of preference:

# 1. masked array — the mask travels with the data through every operation
with rasterio.open("elevation.tif") as src:
    data = src.read(1, masked=True)
    raw = src.read(1)
    nodata = src.nodata
print(f"masked array   {data.mean():.2f}")

# 2. convert to float and use NaN — good for arithmetic
arr = raw.astype("float32")
arr[arr == nodata] = np.nan
print(f"NaN            {np.nanmean(arr):.2f}")

# 3. boolean mask — explicit, fine for one-off statistics
valid = raw != nodata
print(f"boolean mask   {raw[valid].mean():.2f}")

masked array   91.34
NaN            91.34
boolean mask   91.34

4. The transform: where the grid sits

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

with rasterio.open("elevation.tif") as src:
    print(src.transform)

| 5.00, 0.00, 490000.00|
| 0.00,-5.00, 508000.00|
| 0.00, 0.00, 1.00|

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

Converting between the two spaces:

with rasterio.open("elevation.tif") as src:
    x, y = src.xy(row=100, col=250)                 # pixel centre → coordinates
    print(x, y)
    print(src.index(x, y))                          # coordinates → pixel
    print(src.transform * (250, 100))               # the corner of that pixel, not its centre

491252.5 507497.5
(100, 250)
(491250.0, 507500.0)

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, and multiplying the transform by (col, row) gives an edge too, while xy() returns the centre. Forgetting that shifts everything by half a cell, 2.5 m here.

5. The CRS: which earth those coordinates are on

with rasterio.open("elevation.tif") as src:
    print(src.crs, "| projected:", src.crs.is_projected, "| units:", src.crs.linear_units)

with rasterio.open("sentinel.tif") as src:
    print(src.crs, "| projected:", src.crs.is_projected, "| units:", src.crs.linear_units)

from pyproj import Transformer
for crs in ("EPSG:27700", "EPSG:32630"):
    lon, lat = Transformer.from_crs(crs, "EPSG:4326", always_xy=True).transform(490000, 508000)
    print(f"(490000, 508000) in {crs}: {lat:.2f}°N, {lon:.2f}°E")

EPSG:27700 | projected: True | units: metre
EPSG:32630 | projected: True | units: metre
(490000, 508000) in EPSG:27700: 54.46°N, -0.61°E
(490000, 508000) in EPSG:32630: 4.60°N, -3.09°E

Both are in metres, and they are not the same metres: 490000, 508000 is the north-west corner of the terrain model in British National Grid, and a point 4½ degrees north of the equator in UTM zone 30N, just off the coast of West Africa. A transform without a CRS is ambiguous, and src.crs is None — as in the bare file above — means every downstream overlay is a guess. See cannot transform naive geometries for the vector version of the same problem.

A short audit that catches most surprises

The data used on this page

Two open sources over the same 10 × 8 km of coast, from Robin Hood's Bay north towards Whitby. The Environment Agency's LIDAR Composite DTM is served by its Web Coverage Service; asking for it at SCALEFACTOR=0.2 resamples the 1 m model to 5 m on the server, and the file comes back exactly as the service writes it. The four 10 m bands of the clear Sentinel-2 L2A scene of 12 July 2026 come from Element 84's Earth Search catalogue, cut to the same box.

import numpy as np
import rasterio
import requests
from pystac_client import Client
from rasterio.warp import transform_bounds
from rasterio.windows import from_bounds

BOX = (490000, 500000, 500000, 508000)        # British National Grid, metres
wcs = requests.get(
    "https://environment.data.gov.uk/spatialdata/lidar-composite-digital-terrain-model-dtm-1m/wcs",
    params=[("service", "WCS"), ("version", "2.0.1"), ("request", "GetCoverage"),
            ("CoverageId", "13787b9a-26a4-4775-8523-806d13af58fc__Lidar_Composite_Elevation_DTM_1m"),
            ("format", "image/tiff"), ("subset", f"E({BOX[0]},{BOX[2]})"),
            ("subset", f"N({BOX[1]},{BOX[3]})"), ("SCALEFACTOR", "0.2")],
    timeout=300)
wcs.raise_for_status()
with open("elevation.tif", "wb") as f:
    f.write(wcs.content)

item = Client.open("https://earth-search.aws.element84.com/v1") \
    .get_collection("sentinel-2-l2a").get_item("S2B_30UXF_20260712_0_L2A")
bands = []
for asset in ("blue", "green", "red", "nir"):
    with rasterio.open(item.assets[asset].href) as src:
        window = from_bounds(*transform_bounds("EPSG:27700", src.crs, *BOX), src.transform)
        window = window.round_offsets().round_lengths()
        bands.append(src.read(1, window=window))
        profile = src.profile | {"driver": "GTiff", "count": 4, "width": bands[0].shape[1],
                                 "height": bands[0].shape[0], "transform": src.window_transform(window)}
with rasterio.open("sentinel.tif", "w", **profile) as dst:
    dst.write(np.stack(bands))

for name in ("elevation.tif", "sentinel.tif"):
    with rasterio.open(name) as src:
        print(f"{name:14} {src.width} × {src.height} × {src.count}  {src.dtypes[0]}  {src.res[0]:.0f} m")
print(f"scene cloud cover {item.properties['eo:cloud_cover']:.2f}%")

elevation.tif  2000 × 1600 × 1  float32  5 m
sentinel.tif   1011 × 814 × 4  uint16  10 m
scene cloud cover 0.02%

Example 1: a short audit that catches most raster surprises

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 '(no units)'}")
        print(f"  bounds    {tuple(round(v, 1) for v in src.bounds)}")
        covered = 100 * (1 - np.ma.getmaskarray(data).mean())
        print(f"  valid     {covered:.1f}% of cells")
        print(f"  range     {data.min():.2f} … {data.max():.2f}")

audit("elevation.tif")
audit("elevation_uint8.tif")

elevation.tif
  size      2000 × 1600 × 1 band(s)
  dtype     float32
  nodata    -3.4028234663852886e+38
  crs       EPSG:27700
  pixel     5.0 × 5.0 m
  bounds    (490000.0, 500000.0, 500000.0, 508000.0)
  valid     60.6% of cells
  range     -3.10 … 265.99
elevation_uint8.tif
  size      2000 × 1600 × 1 band(s)
  dtype     uint8
  nodata    0.0
  crs       EPSG:27700
  pixel     5.0 × 5.0 m
  bounds    (490000.0, 500000.0, 500000.0, 508000.0)
  valid     59.5% of cells
  range     1.00 … 255.00

Two lines do the real work. valid says how much of the grid is NoData — 61% valid is normal for a coastal tile and would be alarming for an inland one. And the uint8 copy has quietly lost a point of it: every height between 0 and 1 m truncated to 0, which is now the NoData value. range is where the rest of the damage shows — 1 … 255 on terrain that runs from −3 m to 266 m says the values have been through a cast.

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. Here, slope in degrees from the terrain model:

with rasterio.open("elevation.tif") as src:
    profile = src.profile.copy()
    heights = src.read(1, masked=True).astype("float32").filled(np.nan)
    xres, yres = src.res

dz_dy, dz_dx = np.gradient(heights, yres, xres)
slope = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy))).astype("float32")

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, 1)

with rasterio.open("slope.tif") as src:
    s = src.read(1, masked=True)
    print(src.crs, src.transform[:6], src.nodata)
    print(f"slope {s.min():.1f} … {s.max():.1f}°, {(s > 45).sum():,} cells steeper than 45°")

EPSG:27700 (5.0, 0.0, 490000.0, 0.0, -5.0, 508000.0) nan
slope 0.0 … 81.6°, 14,669 cells steeper than 45°

profile already held width, height, transform and CRS, so the output lands exactly on the input. Building a profile from scratch is where people forget the transform and produce a raster that opens at the origin of the coordinate system — the bare file in the problem statement. Passing the cell size to np.gradient is what makes the angles real: without yres, xres it assumes 1 m spacing and every slope on this 5 m grid comes out five times too steep in its tangent. Under 1% of the surveyed ground steeper than 45°, on a coast of cliffs and steep-sided becks, is the plausibility check.

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 Sentinel-2 stack is not four rasters of the same quantity, and the file does not say which is which:

with rasterio.open("sentinel.tif") as src:
    print(src.descriptions)

(None, None, None, None)

Record what each band is when you write it:

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

with rasterio.open("sentinel_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)

with rasterio.open("sentinel_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")
    guessed_red = src.read(1, masked=True).astype("float32")      # "the first band is red"

print(f"NDVI, bands by name:    median {np.ma.median((nir - red) / (nir + red)):.2f}")
print(f"NDVI, red guessed as 1: median {np.ma.median((nir - guessed_red) / (nir + guessed_red)):.2f}")

NDVI, bands by name:    median 0.50
NDVI, red guessed as 1: median 0.66

The wrong band gives a plausible-looking array of wrong numbers: a median of 0.66 instead of 0.50, still between −1 and 1 as any normalised difference is, just not NDVI — band 1 is blue. Band descriptions are the cheapest defence against that.

A rectangle measuring things that are not

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 numbers 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 lost by a conversion to a plain array, a numpy.save, or a naïve copy, and made meaningless by a dtype that cannot hold it, as the uint8 copy showed. 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: cast from a wider integer, class 300 becomes 44, silently colliding with a real class. Because the wrap-around is arithmetic rather than an error, the only defence is checking the range against the dtype before writing.

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.

Missing NoData tags and NaN in integer types

• 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 needs a float dtype. Rasterio refuses it on an integer raster: ValueError: Given nodata value, nan, is beyond the valid range of its data type, int16.

• A NoData tag the dtype cannot hold is refused, as the uint8 copy showed. The usual workaround — nodata=0 — then collides with any real value that truncates to 0.

• 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 — a 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 whole metres, float32 for centimetre-precision lidar heights like the ones on this page, 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.

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.