⚙️ How-To

How to Batch Process Rasters with Rasterio in Python

You have a folder of GeoTIFFs — a season of satellite scenes, a directory of elevation tiles, or a stack of drone orthophotos — and you need to do the same thing to every one of them: reproject to a common CRS, clip to a study area, or compute an index. Doing it by hand in QGIS forty times is not a plan. This guide shows you how to march through that folder with rasterio, apply one operation per file, and — the part beginners get wrong — copy and update the raster profile correctly so the output GeoTIFFs are valid, compressed, and georeferenced. One bad file will not kill the run.

Problem statement

You have a directory of raster files and you need to apply a uniform operation to all of them. The usual pain points:

  • Repetitive manual work — opening each .tif in a GUI, running the same tool, and exporting is slow and error-prone across dozens of files.
  • Broken output metadata — you write a processed array back to disk but forget to update the profile, so the new file has the wrong width, height, transform, or crs and no longer lines up in space.
  • Memory blow-ups on large rasters — calling src.read() on a multi-gigabyte scene loads the whole thing into RAM and the process is killed.
  • Bloated file sizes — output written with no compression is several times larger than it needs to be.
  • One corrupt file stops everything — an unreadable or truncated raster throws halfway through and you lose the whole batch.

The goal: one repeatable script that discovers every raster, opens it safely, applies your operation, writes a correctly georeferenced and compressed GeoTIFF, and hands you a summary of what succeeded and what failed.

Quick answer

Loop over the folder with pathlib, open each raster with rasterio.open() as a context manager, read the band(s) you need, run your operation, then copy the source profile and update() only the keys that changed before writing. Wrap every file in its own try/except so one bad raster does not abort the batch, and turn on compress="deflate" on output.

A folder of GeoTIFF rasters passing through open, read, process, and write, with the profile copy-and-update step highlighted.
Open each raster, read it, process the array, then copy the profile and update only the keys that changed before writing.
import rasterio
from pathlib import Path

src_dir = Path("data/rasters")
out_dir = Path("data/processed")
out_dir.mkdir(parents=True, exist_ok=True)

summary = []
for path in sorted(src_dir.glob("*.tif")):
    try:
        with rasterio.open(path) as src:
            band = src.read(1)                 # read the data
            result = band * 2                  # your per-raster operation

            profile = src.profile.copy()       # copy the metadata
            profile.update(dtype=result.dtype, compress="deflate")

        out_path = out_dir / f"{path.stem}_out.tif"
        with rasterio.open(out_path, "w", **profile) as dst:
            dst.write(result, 1)
        summary.append((path.name, "ok"))
    except Exception as exc:
        summary.append((path.name, f"error: {exc}"))

for name, status in summary:
    print(f"{name:30} {status}")

That is the whole pattern. The rest of this article swaps in real operations — reproject, clip, band math, windowed reads — and explains the profile idiom so you can adapt it without breaking the georeferencing.

Step-by-step solution

Discover the raster files

Use pathlib.Path and glob to find every GeoTIFF. Globbing keeps the file discovery readable and lets you control exactly which extensions you accept. Use rglob instead of glob if your tiles are nested in subfolders.

from pathlib import Path

src_dir = Path("data/rasters")
PATTERNS = ("*.tif", "*.tiff")

files = sorted(p for pat in PATTERNS for p in src_dir.glob(pat))
for f in files:
    print(f.name)

GeoTIFFs come with either .tif or .tiff extensions, so glob both. If your imagery uses another driver-supported format (.img, .jp2), add its pattern too — rasterio reads anything GDAL can.

Open each raster with a context manager

Open every file with rasterio.open() inside a with block. The context manager closes the dataset and releases the underlying GDAL handle automatically, even if an error is raised — essential when you are opening hundreds of files in a loop.

import rasterio

with rasterio.open("data/rasters/scene_01.tif") as src:
    print("Bands:", src.count)
    print("Size:", src.width, "x", src.height)
    print("CRS:", src.crs)
    print("Dtype:", src.dtypes[0])

If you need a refresher on what these properties mean, see the introduction to rasterio. Inspecting one file before you loop is a good habit — it tells you the band count and dtype your loop has to handle.

Read the bands or a window

Read the data you need and nothing more. For a single-band DEM read band 1; for imagery read all bands or just the ones your operation touches. On large rasters, read a Window instead of the whole array to keep memory flat.

import rasterio
from rasterio.windows import Window

with rasterio.open("data/rasters/scene_01.tif") as src:
    full = src.read(1)                          # whole band into memory
    tile = src.read(1, window=Window(0, 0, 512, 512))   # a 512x512 block

Reading a window pulls only the requested block off disk, which is how you process rasters too big to load whole. More on that in the windowed-processing example below.

Copy the profile and update only what changed

This is the step that separates working batch scripts from broken ones. The profile is a dictionary of everything needed to write a valid GeoTIFF: driver, dtype, nodata, width, height, count, crs, and transform. Start from the source profile, then update() only the keys your operation actually changed.

import rasterio

with rasterio.open("data/rasters/scene_01.tif") as src:
    band = src.read(1).astype("float32")
    profile = src.profile.copy()        # inherit crs, transform, size, nodata

# We changed the data type; everything else is unchanged.
profile.update(dtype="float32", compress="deflate")

If you reproject, you update crs, transform, width, and height. If you go from three bands to one, you update count. If you change the number format, you update dtype. Everything you do not touch — extent, pixel size, nodata — is inherited correctly from the source. Reconstructing the profile from scratch is where beginners lose the georeferencing.

Write the output and handle errors per file

Write with rasterio.open(path, "w", **profile) inside its own context manager, and wrap the whole read-process-write cycle for each file in try/except. A truncated tile or an unreadable file then produces a logged error rather than a stack trace that stops the batch.

import rasterio
from pathlib import Path

out_dir = Path("data/processed")
out_dir.mkdir(parents=True, exist_ok=True)

for path in sorted(Path("data/rasters").glob("*.tif")):
    try:
        with rasterio.open(path) as src:
            data = src.read(1)
            profile = src.profile.copy()
            profile.update(compress="deflate")
        with rasterio.open(out_dir / path.name, "w", **profile) as dst:
            dst.write(data, 1)
        print(f"OK   {path.name}")
    except Exception as exc:
        print(f"FAIL {path.name}: {exc}")

Code examples

Example 1: Batch reproject a folder to a target CRS

Reprojection is the classic batch job. For each file, calculate_default_transform() works out the output grid, then you update the profile with the new crs, transform, width, and height, and reproject() each band into the destination dataset.

import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from pathlib import Path

src_dir = Path("data/rasters")
out_dir = Path("data/reprojected")
out_dir.mkdir(parents=True, exist_ok=True)

DST_CRS = "EPSG:4326"

for path in sorted(src_dir.glob("*.tif")):
    with rasterio.open(path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, DST_CRS, src.width, src.height, *src.bounds
        )
        profile = src.profile.copy()
        profile.update(
            crs=DST_CRS, transform=transform,
            width=width, height=height, compress="deflate",
        )

        with rasterio.open(out_dir / path.name, "w", **profile) as dst:
            for i in range(1, src.count + 1):        # reproject every band
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform, src_crs=src.crs,
                    dst_transform=transform, dst_crs=DST_CRS,
                    resampling=Resampling.nearest,
                )
    print(f"reprojected {path.name} -> {DST_CRS}")

Use Resampling.nearest for categorical data such as land cover, and Resampling.bilinear or Resampling.cubic for continuous data such as elevation or reflectance.

Example 2: Batch clip every raster to a boundary

To clip a folder to one study-area polygon, read the boundary once with GeoPandas, then rasterio.mask.mask() each raster. mask returns the clipped array and the new transform, and you update the profile with the new height, width, and transform.

import rasterio
import rasterio.mask
import geopandas as gpd
from pathlib import Path

boundary = gpd.read_file("data/boundary.gpkg")
out_dir = Path("data/clipped")
out_dir.mkdir(parents=True, exist_ok=True)

for path in sorted(Path("data/rasters").glob("*.tif")):
    with rasterio.open(path) as src:
        shapes = boundary.to_crs(src.crs).geometry     # match the raster CRS
        clipped, clip_transform = rasterio.mask.mask(src, shapes, crop=True)

        profile = src.profile.copy()
        profile.update(
            height=clipped.shape[1], width=clipped.shape[2],
            transform=clip_transform, compress="deflate",
        )

    with rasterio.open(out_dir / path.name, "w", **profile) as dst:
        dst.write(clipped)
    print(f"clipped {path.name}")

Reprojecting the boundary to src.crs inside the loop means it lines up with each raster even if the folder contains mixed CRSs. The mask array already has shape (bands, rows, cols), so dst.write(clipped) writes every band at once.

Example 3: Batch NDVI (band math) to a new single-band raster

Band math collapses several input bands into one output, so the key profile changes are count=1 and a float dtype. This example computes NDVI from red and near-infrared bands and writes a single-band GeoTIFF.

import rasterio
import numpy as np
from pathlib import Path

out_dir = Path("data/ndvi")
out_dir.mkdir(parents=True, exist_ok=True)

RED, NIR = 3, 4        # band indices for this sensor (1-based)

for path in sorted(Path("data/rasters").glob("*.tif")):
    with rasterio.open(path) as src:
        red = src.read(RED).astype("float32")
        nir = src.read(NIR).astype("float32")

        denom = nir + red
        ndvi = np.where(denom == 0, np.nan, (nir - red) / denom)

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

    out_path = out_dir / f"{path.stem}_ndvi.tif"
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(ndvi.astype("float32"), 1)
    print(f"NDVI {path.name} -> {out_path.name}")

Casting to float32 before dividing avoids integer truncation, and np.where guards against divide-by-zero where red and NIR are both zero. Setting count=1 is what makes the output single-band.

Example 4: Windowed (block) processing for large rasters

When a raster is too big to load whole, process it block by block. src.block_windows() yields the file's native tiling, so you read, process, and write one window at a time and memory stays flat regardless of the total size.

import rasterio
from pathlib import Path

out_dir = Path("data/scaled")
out_dir.mkdir(parents=True, exist_ok=True)

for path in sorted(Path("data/rasters").glob("*.tif")):
    with rasterio.open(path) as src:
        profile = src.profile.copy()
        profile.update(dtype="float32", compress="deflate", tiled=True)

        with rasterio.open(out_dir / path.name, "w", **profile) as dst:
            for _, window in src.block_windows(1):     # iterate native blocks
                block = src.read(1, window=window).astype("float32")
                dst.write(block * 0.0001, 1, window=window)
    print(f"processed {path.name} in blocks")

Reading and writing by window means only one block is in memory at a time. Setting tiled=True on the output writes an internally tiled GeoTIFF, which makes later windowed reads of your result fast too. For genuinely huge jobs, see speeding up batch jobs with parallel processing.

Example 5: Robust batch loop with a summary and compression

Put it together with per-file error handling and an audit trail. Collect one record per file — status, band count, output size — and print a summary at the end so you know exactly what happened without scrolling through logs.

import rasterio
from pathlib import Path

src_dir = Path("data/rasters")
out_dir = Path("data/processed")
out_dir.mkdir(parents=True, exist_ok=True)

summary = []
for path in sorted(src_dir.glob("*.tif")):
    try:
        with rasterio.open(path) as src:
            data = src.read()
            profile = src.profile.copy()
            profile.update(compress="deflate", predictor=2)

        out_path = out_dir / path.name
        with rasterio.open(out_path, "w", **profile) as dst:
            dst.write(data)

        size_mb = out_path.stat().st_size / 1e6
        summary.append((path.name, "ok", profile["count"], f"{size_mb:.1f} MB"))
    except Exception as exc:
        summary.append((path.name, f"error: {exc}", 0, "-"))

ok = sum(1 for row in summary if row[1] == "ok")
print(f"{ok}/{len(summary)} rasters processed")
for name, status, bands, size in summary:
    print(f"{name:28} {status:16} {bands} bands  {size}")

compress="deflate" with predictor=2 (horizontal differencing) often shrinks continuous rasters substantially with no loss of data. The final tally tells you at a glance whether the batch was clean.

Explanation

The profile dictionary is the contract

Every rasterio dataset carries a profile: a dictionary describing how the pixels map to the world and disk. The keys that matter for writing are:

  • driver — the output format, usually "GTiff".
  • dtype — the pixel data type ("uint8", "int16", "float32", …). It must match the array you write.
  • count — the number of bands.
  • width, height — columns and rows.
  • crs — the coordinate reference system.
  • transform — an affine transform mapping pixel (row, col) to world coordinates, and encoding pixel size.
  • nodata — the value that marks missing pixels.

When you write output, **profile unpacks this dictionary into rasterio.open(..., "w"). The whole batch pattern rests on one rule: inherit the source profile, then override only what your operation changed. That way the output stays correctly georeferenced by default.

Which keys to update for which operation

The changed keys follow directly from what the operation does to the array:

  • Reproject: new grid and CRS, so update crs, transform, width, height.
  • Clip/mask: the array shrinks and shifts, so update width, height, transform.
  • Band math (e.g. NDVI): fewer bands and a new number format, so update count and dtype.
  • Rescale/convert: only the number format changes, so update dtype.

If the array shape you write disagrees with width, height, or count in the profile, rasterio raises an error — a useful guardrail that catches a forgotten update().

Transforms and windows

The transform is an affine matrix. Its first and fifth coefficients (transform.a and transform.e) are the pixel width and height; the third and sixth (transform.c, transform.f) are the coordinates of the top-left corner. Reprojecting or clipping produces a new top-left corner and possibly a new pixel size, which is why those operations hand you a fresh transform to put in the profile.

A Window is a (col_off, row_off, width, height) view into the pixel grid. Reading and writing by window lets you touch a slice of a raster without loading the rest, and rasterio keeps the georeferencing straight because the window is expressed in the dataset's own pixel space.

Edge cases or notes

Preserve or set the nodata value

If your operation introduces missing data — clipping leaves pixels outside the polygon, band math produces NaN — make sure the output nodata value is set in the profile and matches the dtype. NaN only works with a float dtype; for integer rasters pick a sentinel like -9999 and set profile.update(nodata=-9999). A wrong or missing nodata value is why downstream statistics come out skewed.

Match dtype to your data

The output dtype must hold every value you write. Reflectance scaled to 0..1 needs float32, not the source uint16; a classified result fits in uint8. Writing float32 data into a uint8 profile silently truncates and corrupts values. Cast the array and update dtype together so they never disagree.

Big files: window, do not read whole

Never call src.read() on a multi-gigabyte raster in a loop — one file can exhaust RAM. Iterate src.block_windows() (Example 4) so memory use is bounded by one block, not the file. Writing a tiled=True GeoTIFF also means your output supports fast windowed access later.

Mixed CRS or resolution in one folder

A folder can contain rasters in different CRSs or at different pixel sizes. That is fine for a per-file operation like reprojection to a common target — each file is handled on its own terms. But if you plan to stack, mosaic, or compare the outputs, reproject and resample them onto one common grid first, or they will not align pixel-to-pixel. Inspect src.crs and src.res up front to know what you are dealing with.

Compression is free size savings

Add compress="deflate" to the output profile on every write. It is lossless, costs a little CPU, and often halves file size or better. For continuous data, predictor=2 improves the ratio further; for floating-point rasters use predictor=3. There is rarely a reason to write an uncompressed GeoTIFF.

Do not write into the folder you are reading

Write outputs to a separate directory. If out_dir equals src_dir and you reuse filenames, a later iteration can glob and reprocess a file you just wrote, or overwrite an input mid-run. Keep inputs and outputs apart, as every example here does.

FAQ

How do I batch process a folder of GeoTIFFs with rasterio?

Loop over the folder with pathlib and glob("*.tif"), open each file with rasterio.open() inside a with block, read the bands, apply your operation, then copy the source profile, update() the keys that changed, and write with a second rasterio.open(path, "w", **profile). Wrap each file in try/except so one failure does not stop the batch.

Why do my output rasters lose their georeferencing?

Almost always because the profile was rebuilt from scratch or the wrong keys were updated. Always start from src.profile.copy() and change only what your operation touched — crs and transform for reprojection, width/height/transform for clipping, count and dtype for band math. Everything you leave alone is inherited correctly.

How do I process a raster that is too large to fit in memory?

Do not call src.read() on the whole file. Iterate src.block_windows(1) and read, process, and write one window at a time with src.read(1, window=window) and dst.write(block, 1, window=window). Memory use stays bounded by a single block regardless of the total file size.

Which resampling method should I use when reprojecting?

Use Resampling.nearest for categorical rasters such as land cover or classified data, because it never invents new class values. Use Resampling.bilinear or Resampling.cubic for continuous data such as elevation, temperature, or reflectance, where smooth interpolation is appropriate.

How do I keep one corrupt file from stopping the whole batch?

Wrap the read-process-write cycle for each file in its own try/except block, record the exception in a summary list, and continue the loop. The bad file ends up flagged as an error in your report while every other raster is still processed. Never let a single truncated tile abort an overnight run.

How do I compress the output GeoTIFFs?

Add compress="deflate" to the profile before writing: profile.update(compress="deflate"). Deflate is lossless. For continuous data add predictor=2, and for floating-point data add predictor=3, to improve the compression ratio further with no data loss.

How do I change the number of bands in the output?

Set count in the profile to the number of bands you write, and make sure the array shape agrees. Band math that produces a single index sets count=1 and writes with dst.write(array, 1); stacking three inputs into one output sets count=3 and writes an array of shape (3, rows, cols). A mismatch between count and the array raises an error.

Can I clip every raster in a folder to the same boundary?

Yes. Read the boundary once with GeoPandas, then inside the loop reproject it to each raster's CRS with boundary.to_crs(src.crs) and call rasterio.mask.mask(src, shapes, crop=True). Update the profile with the returned transform and the clipped array's new height and width before writing.