How to Batch Process Rasters with Rasterio in Python

On this page

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.

Each raster should come out as a correctly georeferenced, compressed GeoTIFF, with one bad file reported rather than fatal and a summary of what succeeded.

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.

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}")

scene_01.tif                   ok
scene_02.tif                   ok
scene_03.tif                   ok
scene_04.tif                   ok

Four for four — including scene_04.tif, which is 99.9% nodata over Bristol: that September pass is listed at 0.12% cloud because cloud cover is a statistic of the tile's valid pixels, and almost none of the tile was imaged on this orbit. Writing it is not an error, so nothing fails; Example 5 checks the valid fraction and skips it.

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.

Open, read, update the profile, write

Most broken raster output is a profile key that was copied instead of updated.

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)

scene_01.tif
scene_02.tif
scene_03.tif
scene_04.tif

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])

Bands: 4
Size: 1709 x 1634
CRS: EPSG:32630
Dtype: uint16

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}")

OK   scene_01.tif
OK   scene_02.tif
OK   scene_03.tif
OK   scene_04.tif

Batch work on Sentinel-2 scenes of Bristol

Every batch reprojection below is this shape wrapped in a folder loop.

The data used on this page

Real imagery: four Sentinel-2 Level-2A scenes over Bristol from Earth Search (Copernicus Sentinel data 2025), cut to the city's ONS boundary box and stacked into 4-band GeoTIFFs in data/rasters/ — blue, green, red and near-infrared (B02, B03, B04, B08), 10 m, in UTM zone 30N. Three are clear summer and autumn scenes. The fourth, scene_04.tif, is a pass the catalogue lists at 0.12% cloud that holds almost no data over Bristol at all, because the city sits at the edge of that orbit's swath — batches meet files like this. Each band keeps the scale and offset the catalogue gives for turning its integers into reflectance. data/boundary.gpkg is the ONS boundary of Bristol.

from pathlib import Path

import geopandas as gpd
import numpy as np
import rasterio
import requests
from rasterio.windows import from_bounds

Path("data/rasters").mkdir(parents=True, exist_ok=True)
if not Path("data/boundary.gpkg").exists():
    gpd.read_file("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
                  "Local_Authority_Districts_May_2024_Boundaries_UK_BFC/FeatureServer/0/query"
                  "?where=LAD24CD%3D%27E06000023%27&outFields=LAD24CD,LAD24NM&outSR=27700&f=geojson"
                  ).to_file("data/boundary.gpkg")
boundary = gpd.read_file("data/boundary.gpkg")
SCENES = ["S2A_30UWC_20250711_0_L2A", "S2C_30UWC_20250712_0_L2A", "S2B_30UWC_20250922_2_L2A",
          "S2A_30UWC_20250922_1_L2A"]
for i, scene in enumerate(SCENES, start=1):
    path = Path(f"data/rasters/scene_{i:02d}.tif")
    if path.exists():
        continue
    item = requests.get(f"https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/items/{scene}",
                        timeout=60).json()
    bands = []
    for name in ("blue", "green", "red", "nir"):
        with rasterio.open(item["assets"][name]["href"]) as src:
            window = from_bounds(*boundary.to_crs(src.crs).total_bounds, src.transform).round_offsets().round_lengths()
            bands.append(src.read(1, window=window))
            transform, crs = src.window_transform(window), src.crs
    meta = item["assets"]["red"]["raster:bands"][0]
    stack = np.stack(bands)
    with rasterio.open(path, "w", driver="GTiff", width=stack.shape[2], height=stack.shape[1], count=4,
                       dtype="uint16", crs=crs, transform=transform, nodata=0, compress="deflate", tiled=True) as dst:
        dst.write(stack)
        dst.scales, dst.offsets = (meta["scale"],) * 4, (meta["offset"],) * 4
        dst.descriptions = ("B02 blue", "B03 green", "B04 red", "B08 nir")
        dst.update_tags(scene=scene, date=item["properties"]["datetime"][:10],
                        cloud_cover=item["properties"]["eo:cloud_cover"],
                        boa_offset_applied=item["properties"].get("earthsearch:boa_offset_applied"))

for path in sorted(Path("data/rasters").glob("*.tif")):
    with rasterio.open(path) as src:
        valid = (src.read_masks(1) > 0).mean()
        print(f"{path.name}  {src.tags()['date']}  {src.width}×{src.height}×{src.count}  "
              f"cloud {float(src.tags()['cloud_cover']):4.2f}%  valid pixels {valid:.1%}  "
              f"scale {src.scales[0]} offset {src.offsets[0]} (already applied: {src.tags()['boa_offset_applied']})")

scene_01.tif  2025-07-11  1709×1634×4  cloud 2.39%  valid pixels 100.0%  scale 0.0001 offset -0.1 (already applied: True)
scene_02.tif  2025-07-12  1709×1634×4  cloud 0.00%  valid pixels 100.0%  scale 0.0001 offset -0.1 (already applied: True)
scene_03.tif  2025-09-22  1709×1634×4  cloud 2.84%  valid pixels 100.0%  scale 0.0001 offset -0.1 (already applied: True)
scene_04.tif  2025-09-22  1709×1634×4  cloud 0.12%  valid pixels 0.1%  scale 0.0001 offset -0.1 (already applied: True)

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}")

reprojected scene_01.tif -> EPSG:4326
reprojected scene_02.tif -> EPSG:4326
reprojected scene_03.tif -> EPSG:4326
reprojected scene_04.tif -> EPSG:4326

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}")

clipped scene_01.tif
clipped scene_02.tif
clipped scene_03.tif
clipped scene_04.tif

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}")

NDVI scene_01.tif -> scene_01_ndvi.tif
NDVI scene_02.tif -> scene_02_ndvi.tif
NDVI scene_03.tif -> scene_03_ndvi.tif
NDVI scene_04.tif -> scene_04_ndvi.tif

Those integers are reflectance × 10,000, which cancels out of a ratio like NDVI. The catalogue also lists an offset of −0.1 for every band: the +1,000 that ESA has added to Level-2A values since processing baseline 04.00 in January 2022. These files come from Earth Search, though, which has already taken it off — its items say earthsearch:boa_offset_applied: true — and still list the offset. Apply it anyway and see what happens:

with rasterio.open("data/rasters/scene_01.tif") as src:
    red, nir = (src.read(b, masked=True).astype("float32") for b in (RED, NIR))
    scale, offset = src.scales[RED - 1], src.offsets[RED - 1]

ndvi = (nir - red) / (nir + red)
red_r, nir_r = red * scale + offset, nir * scale + offset          # the listed offset, applied again
ndvi_twice = (nir_r - red_r) / (nir_r + red_r)
print(f"median NDVI {np.ma.median(ndvi):.3f}; with the offset applied a second time "
      f"{np.ma.median(ndvi_twice):.3f}, and {(ndvi_twice > 1).mean():.0%} of pixels above 1")

median NDVI 0.437; with the offset applied a second time 1.079, and 55% of pixels above 1

An NDVI above 1 cannot exist. Scale before any calculation that is not a ratio, and apply an offset only when the data still carries it — Earth Search records that in earthsearch:boa_offset_applied; files straight from ESA carry it from baseline 04.00 on, and the product metadata's BOA_ADD_OFFSET says by how much.

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.

The catalogue’s offset had already been applied to these files; applying it again pushes NDVI above 1 for 55% of pixels.

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:
        scale = src.scales[0]                                   # from the file, not from memory
        offset = 0.0 if src.tags().get("boa_offset_applied") == "True" else src.offsets[0]
        profile = src.profile.copy()
        profile.update(count=1, dtype="float32", nodata=float("nan"), 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(3, window=window, masked=True).astype("float32")
                dst.write((block * scale + offset).filled(float("nan")), 1, window=window)
    print(f"processed {path.name} in {len(list(src.block_windows(1)))} blocks")

with rasterio.open("data/rasters/scene_01.tif") as src:
    red = src.read(3, masked=True).astype("float32")
print(f"median red reflectance over Bristol: {np.ma.median(red * 0.0001):.4f}; "
      f"with the listed offset applied again: {np.ma.median(red * 0.0001 - 0.1):.4f}")

processed scene_01.tif in 49 blocks
processed scene_02.tif in 49 blocks
processed scene_03.tif in 49 blocks
processed scene_04.tif in 49 blocks
median red reflectance over Bristol: 0.0950; with the listed offset applied again: -0.0050

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()
            valid = (src.read_masks(1) > 0).mean()
            profile = src.profile.copy()
            profile.update(compress="deflate", predictor=2)
        if valid < 0.5:                                    # an empty pass is not a success
            summary.append((path.name, "skipped: no data", profile["count"], f"{valid:.1%} valid"))
            continue

        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}")

3/4 rasters processed
scene_01.tif                 ok               4 bands  16.1 MB
scene_02.tif                 ok               4 bands  16.0 MB
scene_03.tif                 ok               4 bands  16.4 MB
scene_04.tif                 skipped: no data 4 bands  0.1% valid

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.

The profile as a contract

The profile dictionary is the contract

The loop looks the same from a distance; everything inside it is different.

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.

Nodata, dtypes and very large rasters

The second one is silent — writing float data into an integer profile truncates it.

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.

Internal links

How to Batch Process GIS Files in Python — the batch-processing cornerstone that this raster workflow specialises.
Introduction to Rasterio: Reading Raster Data in Python — the fundamentals of opening and reading the files you are looping over.
Extract Raster Values at Point Locations with Rasterio — a common per-raster operation to run inside the same loop.
How to Convert a Raster to a Vector in Python — polygonise each processed raster as a follow-on batch step.
How to Batch Process Multiple Shapefiles in Python — the vector equivalent of this pattern.
Speed Up Batch Jobs with Parallel Processing — run the per-file work across multiple cores once the folder is large.

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