Introduction to Rasterio: Reading Raster Data in Python

On this page

Problem statement

Before any raster analysis, open the file and read what it actually is. A GeoTIFF of elevation, land cover or imagery has to tell you:

the raster size
how many bands it has
its coordinate reference system
its extent
its nodata value
the actual pixel values

This is usually the first step in a larger workflow such as clipping, reprojection, raster calculations, or batch processing. If you do not inspect the raster first, it is easy to use the wrong band, miss a missing CRS, or read data into the wrong array shape.

This page is a practical Rasterio tutorial focused only on reading raster data and inspecting basic properties. It does not cover reprojection, masking, or raster algebra. Every example reads a real file: the Copernicus 30 m DEM, ESA WorldCover land cover and a four-band Sentinel-2 image, all over Rutland.

Quick answer

Bands are read with 1-based numbers; the NumPy array is 0-based with shape (bands, rows, columns).

Use Rasterio to open the raster with a context manager, inspect its metadata, and read one band or all bands into NumPy arrays.

import rasterio

raster_path = "data/dem.tif"

with rasterio.open(raster_path) as src:
    print("Width:", src.width)
    print("Height:", src.height)
    print("Band count:", src.count)
    print("CRS:", src.crs)
    print("Bounds:", src.bounds)
    print("Transform:", src.transform)
    print("Dtype:", src.dtypes[0])
    print("NoData:", src.nodata)
    print("Resolution:", src.res)

    band1 = src.read(1)
    print("Band 1 shape:", band1.shape)

Width: 1008
Height: 900
Band count: 1
CRS: EPSG:4326
Bounds: BoundingBox(left=-0.84, bottom=52.52, right=-0.42, top=52.77)
Transform: | 0.00, 0.00,-0.84|
| 0.00,-0.00, 52.77|
| 0.00, 0.00, 1.00|
Dtype: float32
NoData: None
Resolution: (0.00041666666666666664, 0.0002777777777777778)
Band 1 shape: (900, 1008)

This is the standard way to open a raster with Rasterio and inspect a GeoTIFF in Python. Note the resolution: the GLO-30 DEM's pixels are 1.5 arc-seconds wide and 1 tall at this latitude, which only src.res makes obvious.

Install, open, inspect metadata, read a band

Four pieces of metadata are what separate a raster from an ordinary image.

Install Rasterio

In a working GIS Python environment, Rasterio is commonly installed with pip or conda.

pip install rasterio

or:

conda install -c conda-forge rasterio

The pip wheels bundle their own GDAL, so the version Rasterio reports is not necessarily the GDAL on your system. If installation fails, handle that separately before continuing.

Open a raster file with Rasterio

Use rasterio.open() to open a GeoTIFF or another supported raster format. A context manager is the safest option because it closes the file automatically.

raster_path = "data/landcover.tif"

with rasterio.open(raster_path) as src:
    print(src)
print("closed after the with block:", src.closed)

<open DatasetReader name='data/landcover.tif' mode='r'>
closed after the with block: True

If the file opens successfully, src is a Rasterio dataset object.

Check raster metadata

Before reading cell values, inspect the core dataset properties.

with rasterio.open(raster_path) as src:
    print("Width:", src.width)
    print("Height:", src.height)
    print("Band count:", src.count)
    print("CRS:", src.crs)
    print("Bounds:", src.bounds)
    print("Transform:", src.transform)
    print("Dtypes:", src.dtypes)
    print("NoData:", src.nodata)
    print("Resolution:", src.res)
    print("Colour map on band 1:", len(src.colormap(1)), "entries")

Width: 2040
Height: 960
Band count: 1
CRS: EPSG:4326
Bounds: BoundingBox(left=-0.7199999999999998, bottom=52.620000000000005, right=-0.5499999999999998, top=52.7)
Transform: | 0.00, 0.00,-0.72|
| 0.00,-0.00, 52.70|
| 0.00, 0.00, 1.00|
Dtypes: ('uint8',)
NoData: 0.0
Resolution: (8.333333333333333e-05, 8.333333333333333e-05)
Colour map on band 1: 256 entries

These properties answer practical questions:

width, height: number of columns and rows
count: number of raster bands
crs: spatial reference system
bounds: raster extent as left, bottom, right, top coordinates
transform: maps pixel positions to real-world coordinates
dtypes: data type for each band
nodata: missing-data value

Two things to notice in that output. The printed transform rounds to two decimals, so a pixel a few hundred-thousandths of a degree wide shows as 0.00; read src.res for the pixel size. And the land-cover file also carries a colour map — WorldCover's class colours travel inside the GeoTIFF, which is why it displays in the right colours in QGIS with no style.

Read one raster band

For a single-band DEM or one specific spectral band, read only that band.

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

print("Array shape:", elevation.shape, elevation.dtype)
print("Top-left value:", elevation[0, 0])

Array shape: (900, 1008) float32
Top-left value: 97.08997

This is the usual way to read a GeoTIFF in Python when only one band is needed.

Read all raster bands

For multiband imagery, read the full stack.

with rasterio.open("data/sentinel_subset.tif") as src:
    all_bands = src.read()
    print("Band descriptions:", src.descriptions)

print("Array shape:", all_bands.shape)

Band descriptions: ('blue', 'green', 'red', 'nir')
Array shape: (4, 800, 1200)

The returned shape is (bands, rows, columns) — four Sentinel-2 bands here, not the (rows, columns, bands) that image libraries and matplotlib's imshow expect.

The files this page reads, as rasterio reported them: each has its own CRS, dtype and idea of missing data.

Understand Rasterio band indexing

Rasterio uses 1-based band indexing when reading from the dataset, and NumPy uses 0-based indexing once the data is loaded:

with rasterio.open("data/sentinel_subset.tif") as src:
    first_band = src.read(1)                    # first band: 1
    print("src.read(1) is all_bands[0]:", bool((first_band == all_bands[0]).all()))
    try:
        src.read(0)
    except Exception as exc:
        print(type(exc).__name__ + ":", exc)

src.read(1) is all_bands[0]: True
IndexError: band index 0 out of range (not in (1, 2, 3, 4))

This difference causes mistakes in many beginner workflows; at least read(0) fails loudly rather than returning the wrong band.

Print basic raster statistics or sample values

Basic validation helps confirm that the raster loaded correctly — and shows what nodata does to it. dem_clipped.tif is the same DEM with everything outside Rutland set to its nodata value, -9999:

import numpy as np

with rasterio.open("data/dem_clipped.tif") as src:
    band1 = src.read(1)
    masked = src.read(1, masked=True)
    print("NoData:", src.nodata)

print("Plain read:  min", np.min(band1), " max", round(float(np.max(band1)), 1))
print("masked=True: min", round(float(masked.min()), 1), " max", round(float(masked.max()), 1),
      f"| {int(masked.mask.sum()):,} of {masked.size:,} cells masked")
print("Sample values:", band1[0:2, 0:3].round(1).tolist())

NoData: -9999.0
Plain read:  min -9999.0  max 197.6
masked=True: min 14.0  max 197.6 | 344,768 of 796,835 cells masked
Sample values: [[-9999.0, -9999.0, -9999.0], [-9999.0, -9999.0, -9999.0]]

The plain read reported a minimum elevation of -9999 m for a county whose lowest ground is well above sea level. masked=True returns a NumPy masked array with the nodata cells hidden, and its statistics are the real ones.

Core properties of a Copernicus DEM

The data used on this page

The Copernicus GLO-30 DEM (© DLR e.V. 2010–2014 and © Airbus Defence and Space GmbH 2014–2018, provided under COPERNICUS by the European Union and ESA), ESA WorldCover 10 m 2021 (© ESA WorldCover project 2021, CC BY 4.0) and one Sentinel-2 L2A scene of 19 June 2025 (contains modified Copernicus Sentinel data 2025), all read as windows over Rutland from their cloud-optimised GeoTIFFs on AWS. The county boundary is the ONS December 2024 district boundary.

from pathlib import Path

import numpy as np
import rasterio
import requests
from osgeo import gdal
from rasterio.windows import from_bounds

gdal.UseExceptions()
Path("data").mkdir(exist_ok=True)
RUTLAND = (-0.84, 52.52, -0.42, 52.77)


def window_copy(url, path, bounds, descriptions=None):
    with rasterio.open(url) as src:
        window = from_bounds(*bounds, src.transform)
        data = src.read(window=window)
        profile = src.profile | {"driver": "GTiff", "height": data.shape[1], "width": data.shape[2],
                                 "transform": src.window_transform(window), "tiled": False, "compress": "deflate"}
        profile.pop("blockxsize", None), profile.pop("blockysize", None)
        colormap = src.colormap(1) if src.count == 1 and src.dtypes[0] == "uint8" else None
    with rasterio.open(path, "w", **profile) as dst:
        dst.write(data)
        if colormap:
            dst.write_colormap(1, colormap)


if not Path("data/sentinel_subset.tif").exists():
    window_copy("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N52_00_W001_00_DEM/"
                "Copernicus_DSM_COG_10_N52_00_W001_00_DEM.tif", "data/dem.tif", RUTLAND)
    window_copy("https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/"
                "ESA_WorldCover_10m_2021_v200_N51W003_Map.tif", "data/landcover.tif", (-0.72, 52.62, -0.55, 52.70))
    gdal.Warp("data/elevation_utm.tif", "data/dem.tif", dstSRS="EPSG:32630", xRes=30, yRes=30, resampleAlg="bilinear")
    Path("data/rutland.geojson").write_bytes(requests.get(
        "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
        "Local_Authority_Districts_December_2024_Boundaries_UK_BFC/FeatureServer/0/query",
        params={"where": "LAD24NM='Rutland'", "outFields": "LAD24NM", "outSR": 4326, "f": "geojson"}, timeout=120).content)
    gdal.Warp("data/dem_clipped.tif", "data/dem.tif", cutlineDSName="data/rutland.geojson", cropToCutline=True,
              dstNodata=-9999)
    scene = "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/30/U/XD/2025/6/S2C_30UXD_20250619_0_L2A/"
    bands = {"blue": "B02", "green": "B03", "red": "B04", "nir": "B08"}
    arrays = []
    for name, code in bands.items():
        with rasterio.open(scene + code + ".tif") as src:
            window = from_bounds(654000, 5831000, 666000, 5839000, src.transform)     # Rutland Water, UTM 30N
            arrays.append(src.read(1, window=window))
            profile = src.profile | {"driver": "GTiff", "count": 4, "height": arrays[-1].shape[0],
                                     "width": arrays[-1].shape[1], "transform": src.window_transform(window),
                                     "tiled": False, "compress": "deflate"}
    profile.pop("blockxsize", None), profile.pop("blockysize", None)
    with rasterio.open("data/sentinel_subset.tif", "w", **profile) as dst:
        dst.write(np.stack(arrays))
        dst.descriptions = tuple(bands)
print("rasterio", rasterio.__version__, "with GDAL", rasterio.__gdal_version__, "|",
      sorted(p.name for p in Path("data").glob("*.tif")))

rasterio 1.5.1 with GDAL 3.12.4 | ['dem.tif', 'dem_clipped.tif', 'elevation_utm.tif', 'landcover.tif', 'sentinel_subset.tif']

Example 1: Open a GeoTIFF and print core dataset properties

raster_path = "data/elevation_utm.tif"

with rasterio.open(raster_path) as src:
    print("CRS:", src.crs)
    print("Bounds:", tuple(round(v) for v in src.bounds))
    print("Width:", src.width)
    print("Height:", src.height)
    print("Band count:", src.count)
    print("Dtype:", src.dtypes[0])
    print("Resolution:", src.res, "(metres, because the CRS is UTM)")

CRS: EPSG:32630
Bounds: (645718, 5821066, 675028, 5849806)
Width: 977
Height: 958
Band count: 1
Dtype: float32
Resolution: (30.0, 30.0) (metres, because the CRS is UTM)

Example 2: Read the first raster band into a NumPy array

with rasterio.open(raster_path) as src:
    band1 = src.read(1)
    nodata = src.nodata

print("Band 1 shape:", band1.shape)
print("Cell [0, 0]:", band1[0, 0], "| NoData:", nodata)
print("First 5 values in first row:", band1[0, :5].round(1).tolist())

Band 1 shape: (958, 977)
Cell [0, 0]: 0.0 | NoData: None
First 5 values in first row: [0.0, 0.0, 0.0, 0.0, 0.0]

Reprojecting the DEM to UTM left a margin of empty cells around the rotated grid, and because the source declared no nodata value, the warp filled that margin with 0.0 and declared none either. The first row reads as elevations of zero — valid-looking numbers that would drag down any mean. It is a good reason to look at a few values before trusting a raster, and to pass dstNodata when you warp.

Example 3: Read all bands from a multiband raster

with rasterio.open("data/sentinel_subset.tif") as src:
    data = src.read()
    names = src.descriptions

print("Full array shape:", data.shape)
print("First band shape:", data[0].shape)
red, nir = data[names.index("red")].astype("float32"), data[names.index("nir")].astype("float32")
ndvi = (nir - red) / (nir + red)
print(f"Median NDVI over the window: {np.nanmedian(ndvi):.2f}")

Full array shape: (4, 800, 1200)
First band shape: (800, 1200)
Median NDVI over the window: 0.71

Looking bands up by their description rather than by position is what keeps the NDVI right if the file's band order ever changes.

Example 4: Check nodata and metadata before processing

with rasterio.open("data/dem_clipped.tif") as src:
    print("NoData:", src.nodata)
    print("Transform:", src.transform)
    print("Profile:")
    for key, value in src.profile.items():
        print(f"  {key}: {value}")

NoData: -9999.0
Transform: | 0.00, 0.00,-0.82|
| 0.00,-0.00, 52.76|
| 0.00, 0.00, 1.00|
Profile:
  driver: GTiff
  dtype: float32
  nodata: -9999.0
  width: 943
  height: 845
  count: 1
  crs: EPSG:4326
  transform: | 0.00, 0.00,-0.82|
| 0.00,-0.00, 52.76|
| 0.00, 0.00, 1.00|
  blockxsize: 943
  blockysize: 2
  tiled: False
  interleave: band

The profile dictionary is useful when you later write output rasters with matching metadata.

Example 5: Pixel coordinates and windows

from rasterio.windows import Window

with rasterio.open("data/elevation_utm.tif") as src:
    row, col = src.height // 2, src.width // 2
    x, y = src.xy(row, col)
    print(f"centre pixel ({row}, {col}) -> x={x:.0f}, y={y:.0f}; back to {src.index(x, y)}")
    block = src.read(1, window=Window(col - 2, row - 2, 5, 5))
    print("5 x 5 window around it:", block.shape, "min", round(float(block.min()), 1), "max", round(float(block.max()), 1))

centre pixel (479, 488) -> x=660373, y=5835421; back to (479, 488)
5 x 5 window around it: (5, 5) min 81.5 max 81.5

src.xy() gives the coordinates of a pixel's centre, src.index() goes back, and a window reads only the block you ask for — the way to work with rasters too large to load whole. The centre of this DEM happens to fall on Rutland Water, which is why all 25 cells of the window hold the same elevation: the DEM records the reservoir's flat surface.

What Rasterio gives you

Rasterio is a Python library for reading and writing raster datasets such as GeoTIFFs. In practical GIS work, it gives you direct access to both metadata and raster values.

This is why you keep the dataset open rather than passing the bare array around.

The three rasters on this page cover the usual shapes:

a single-band DEM with elevation values
multiband satellite imagery with multiple spectral bands
a classified raster where each cell stores a class code

Metadata tells you how the raster fits into space. That includes the CRS, bounds, transform, dimensions, band count, and nodata value. Raster values are the actual cell values stored in each band.

CRS and transform are especially important. A raster may open without errors, but if the CRS is missing or wrong, any later overlay with vectors or other rasters may fail. The transform tells you where each pixel sits in real-world coordinates and what the pixel size is.

If your goal is automation, reading metadata first is a reliable habit. It helps catch problems before you batch-process many files.

Unexpected values, CRS issues and memory

Raster file opens but data looks unexpected

Usually one of these:

you read the wrong band
nodata values are affecting min and max statistics, as the -9999 above did
you expected (rows, columns, bands) but Rasterio returned (bands, rows, columns)

CRS issues

Some rasters have missing, incomplete, or incorrect CRS information. Always inspect src.crs before spatial analysis. If the CRS is unclear, confirm it from the source data before continuing.

Large rasters and memory use

Reading an entire raster into memory may not be practical for large files. The setup of this page never downloaded a whole tile: each file is a window read straight from a cloud-optimised GeoTIFF, the same window= mechanism as example 5.

Not all rasters are GeoTIFFs

GeoTIFF is the most common starting format, but Rasterio supports other raster formats through GDAL drivers — see not recognised as a supported file format when one will not open.

Internal links

FAQ

How do I read a GeoTIFF in Python with Rasterio?

Open it with rasterio.open() and read a band with src.read(1).

import rasterio

with rasterio.open("data/dem.tif") as src:
    band1 = src.read(1)

Why does Rasterio use band numbers starting at 1?

Rasterio follows GIS raster conventions where bands are numbered starting from 1 at the dataset level. NumPy arrays still use 0-based indexing after the data is loaded.

What is the difference between raster metadata and raster values?

Metadata describes the raster structure and spatial reference, such as CRS, bounds, transform, and nodata. Raster values are the actual cell values stored in each band.

Should I read the whole raster at once or one band at a time?

Read one band if that is all you need. Read the full raster only when you actually need all bands and the dataset fits comfortably in memory.

How do I get the real-world coordinates of a specific pixel in Rasterio?

Use the dataset's transform via src.xy(row, col), which returns the coordinate of that pixel's center in the raster CRS. The inverse, src.index(x, y), converts a coordinate back to a (row, col) pair.

How do I read just a small region of a large raster without loading it all?

Read a window with src.read(1, window=Window(col_off, row_off, width, height)) after from rasterio.windows import Window. Windowed reads pull only the requested block into memory, which is the standard approach for rasters too large to load whole.

How do I find the pixel size (resolution) of a raster in Rasterio?

Use src.res, which returns (x_size, y_size) in the raster's CRS units, or read the pixel width and height from src.transform.a and src.transform.e. For a raster in a UTM CRS such as EPSG:32633 these are metres; for EPSG:4326 they are degrees.

Why are the min and max values of my band wrong or extreme?

Nodata cells are included in plain src.read(1) statistics, so a sentinel value like -9999 distorts the minimum and maximum. Read with src.read(1, masked=True) so nodata is excluded, then compute band.min() and band.max() on the masked array.