Reprojecting a Raster with rasterio.warp.reproject

Q: Which resampling method should I use when reprojecting a raster?

Match the method to the data type. Use Resampling.nearest for categorical rasters such as land-cover or soil-class maps, because it copies existing integer codes without inventing new ones. Use Resampling.bilinear or Resampling.cubic for continuous data such as elevation, reflectance, or NDVI, where interpolated values are physically meaningful. Never use bilinear on class labels — averaging class 3 and class 7 into 5 produces a category that does not exist.

Q: How do I keep nodata pixels from becoming valid data after reprojection?

Pass src_nodata and dst_nodata explicitly to reproject and set the same value in the output profile. rasterio then fills pixels outside the source footprint with the nodata value instead of zero, and prevents interpolation from blending nodata into neighbouring valid pixels. Write the nodata value into the profile with dst_profile['nodata'] so downstream readers and masks recognise it.

TL;DR Answer

Call calculate_default_transform to derive the target grid (transform, width, height) for the destination CRS, allocate an output array, then loop over bands calling rasterio.warp.reproject with an explicit Resampling method and matching src_nodata/dst_nodata. Update the source profile with the new CRS, transform, and dimensions before writing the output GeoTIFF. This keeps every band on one aligned grid — essential for stacking rasters into a clean feature array. It is a core step in CRS Alignment and Projection Handling for geospatial ML.

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

dst_crs = "EPSG:32633"  # UTM zone 33N

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

    with rasterio.open("output_utm.tif", "w", **profile) as dst:
        for i in range(1, src.count + 1):
            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.bilinear,
            )

Part of: CRS Alignment and Projection Handling

Why Reprojection Matters Before You Train a Model

A raster stored in geographic coordinates such as EPSG:4326 measures position in degrees, not metres. That is fine for display, but it breaks almost every spatial feature you might compute for a model. Pixel area varies with latitude, so a “per-pixel” NDVI statistic is not comparable between a scene in Norway and one in Kenya. Distance-based operations — buffers, focal windows, neighbourhood statistics — assume metric units and produce nonsense on degrees. Any model that consumes elevation slope, texture, or spatial lag features needs its inputs in an equal-area or metric projection like a UTM zone first.

Reprojecting a raster is not the same as reprojecting a vector layer. A vector reprojection moves a finite set of coordinates through a transform. A raster reprojection must resample a regular grid of pixels onto a new regular grid whose cell boundaries do not line up with the old ones. Every output pixel samples one or more input pixels, and the rule you choose for that sampling — nearest neighbour, bilinear, cubic — directly changes the values your model sees. Get it wrong and you either smear categorical labels into impossible classes or introduce blocky artefacts into continuous surfaces.

The corruption is quiet. A land-cover raster reprojected with bilinear resampling will still open, still render, and still stack — but its class codes will contain averaged values that map to the wrong legend entries. A drift-detection baseline or an NDVI and EVI calculation built on such a raster is silently wrong from the first epoch. This is why raster reprojection deserves the same care you give to projection mismatches in GeoDataFrames.

Core Principles for a Sound Raster Reprojection

Match the resampling method to the data type. Use Resampling.nearest for categorical rasters (land cover, soil class, flags) so no new codes are invented. Use Resampling.bilinear or Resampling.cubic for continuous surfaces (elevation, reflectance, temperature) where interpolation is physically meaningful.
Derive the grid once with calculate_default_transform. It returns the output transform, width, and height that best cover the source extent in the destination CRS. Do not hand-guess dimensions.
Preserve nodata explicitly. Pass src_nodata and dst_nodata to reproject and write the value into the output profile. Otherwise pixels outside the source footprint fill with zero and pollute your statistics.
Share one grid across all bands. When you plan to stack bands into a feature cube, reproject every input onto an identical target transform, width, and height. Independent per-band grids drift out of registration.
Update the full profile before writing. Copy the source profile and update crs, transform, width, height, and nodata. A stale transform or width silently truncates or offsets the output.
Verify after writing. Reopen the output and assert its CRS and transform equal what you intended. A five-line check catches a broken pipeline before it reaches training.

Production-Ready Code

The function below wraps the whole operation with input validation, explicit nodata handling, typed signatures, and logging. It reprojects every band of a source GeoTIFF onto a single shared grid, so the result is immediately safe to stack with other rasters reprojected the same way.

import logging
from pathlib import Path

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

logger = logging.getLogger(__name__)


def reproject_raster(
    src_path: str | Path,
    dst_path: str | Path,
    dst_crs: str,
    resampling: Resampling = Resampling.bilinear,
    dst_nodata: float | None = None,
) -> Path:
    """Reproject a GeoTIFF to a new CRS onto a single aligned grid.

    Parameters
    ----------
    src_path : str or Path
        Path to the source raster.
    dst_path : str or Path
        Path to write the reprojected GeoTIFF.
    dst_crs : str
        Target CRS, e.g. "EPSG:32633" for UTM zone 33N.
    resampling : Resampling, default=Resampling.bilinear
        Use Resampling.nearest for categorical data (land cover, class codes)
        and Resampling.bilinear or Resampling.cubic for continuous surfaces.
    dst_nodata : float or None, default=None
        Nodata value for the output. Falls back to the source nodata.

    Returns
    -------
    Path
        The path to the written raster.
    """
    src_path, dst_path = Path(src_path), Path(dst_path)
    if not src_path.exists():
        raise FileNotFoundError(f"Source raster not found: {src_path}")

    with rasterio.open(src_path) as src:
        if src.crs is None:
            raise ValueError(
                f"{src_path} has no CRS. Assign one before reprojecting."
            )

        nodata = dst_nodata if dst_nodata is not None else src.nodata
        if nodata is None and resampling is not Resampling.nearest:
            logger.warning(
                "No nodata value set; edge pixels will fill with 0. "
                "Pass dst_nodata to avoid contaminating band statistics."
            )

        # Derive the target grid once — reused for every band.
        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,
            nodata=nodata,
        )

        logger.info(
            "Reprojecting %s (%s) -> %s onto %dx%d grid using %s",
            src_path.name, src.crs, dst_crs, width, height, resampling.name,
        )

        with rasterio.open(dst_path, "w", **profile) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    src_nodata=src.nodata,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    dst_nodata=nodata,
                    resampling=resampling,
                )

    logger.info("Wrote %s", dst_path)
    return dst_path

The dst_nodata argument threads all the way through to the output profile and to reproject itself, so pixels beyond the source footprint are labelled instead of quietly set to zero. The warning branch flags the one combination that most often produces contaminated features: a continuous raster with no nodata value.

Step-by-Step Walkthrough

1. Inspect the source and choose the destination CRS

import rasterio

with rasterio.open("landsat_ndvi_4326.tif") as src:
    print(f"CRS: {src.crs}")            # EPSG:4326 — geographic degrees
    print(f"Size: {src.width}x{src.height}")
    print(f"Bands: {src.count}, dtype: {src.dtypes[0]}")
    print(f"Bounds: {src.bounds}")
    print(f"Nodata: {src.nodata}")

Pick a metric CRS that fits your area of interest. For a single Landsat scene, the covering UTM zone (here EPSG:32633) minimises area and shape distortion far better than a global projection like EPSG:3857.

2. Compute the target grid

from rasterio.warp import calculate_default_transform

dst_crs = "EPSG:32633"

with rasterio.open("landsat_ndvi_4326.tif") as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds
    )

print(f"Target size: {width}x{height}")
print(f"Target transform:\n{transform}")

calculate_default_transform returns an affine transform describing pixel size and origin in the new CRS. Passing *src.bounds unpacks the source extent so the output grid covers exactly the same ground area.

3. Reproject with the right resampling method

For a continuous NDVI raster, bilinear resampling is appropriate:

from reproject_utils import reproject_raster
from rasterio.warp import Resampling

reproject_raster(
    "landsat_ndvi_4326.tif",
    "landsat_ndvi_utm.tif",
    dst_crs="EPSG:32633",
    resampling=Resampling.bilinear,
    dst_nodata=-9999.0,
)

For a categorical land-cover raster, switch a single argument — nearest neighbour preserves the class codes exactly:

reproject_raster(
    "landcover_4326.tif",
    "landcover_utm.tif",
    dst_crs="EPSG:32633",
    resampling=Resampling.nearest,
    dst_nodata=255,
)

4. Reproject a whole stack onto one shared grid

When bands come from separate files, compute the grid from a reference raster and force every band onto it. This is the step that guarantees pixel-for-pixel registration for a model-ready feature cube.

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

band_paths = ["blue.tif", "green.tif", "red.tif", "nir.tif"]
dst_crs = "EPSG:32633"

# Derive one grid from the first band and reuse it for all others.
with rasterio.open(band_paths[0]) as ref:
    transform, width, height = calculate_default_transform(
        ref.crs, dst_crs, ref.width, ref.height, *ref.bounds
    )
    profile = ref.profile.copy()

profile.update(crs=dst_crs, transform=transform, width=width,
               height=height, count=len(band_paths))

with rasterio.open("stack_utm.tif", "w", **profile) as dst:
    for idx, path in enumerate(band_paths, start=1):
        with rasterio.open(path) as src:
            reproject(
                source=rasterio.band(src, 1),
                destination=rasterio.band(dst, idx),
                src_transform=src.transform,
                src_crs=src.crs,
                src_nodata=src.nodata,
                dst_transform=transform,
                dst_crs=dst_crs,
                dst_nodata=src.nodata,
                resampling=Resampling.bilinear,
            )

Every band now shares the same transform, width, and height, so dst.read() returns a clean (bands, rows, cols) array with no offset between layers.

How the Reprojection Pipeline Fits Together

The diagram traces a source raster in EPSG:4326 through the two functions that do the real work and out to an aligned GeoTIFF. calculate_default_transform sizes the target grid; reproject samples every source pixel onto it using the resampling rule you pick for the data type.

Verification

Never trust a reprojection you have not reopened. The check below asserts the output CRS and grid match the target and confirms nodata survived the round trip.

import rasterio
from rasterio.warp import calculate_default_transform


def verify_reprojection(
    src_path: str, dst_path: str, expected_crs: str
) -> None:
    """Assert the output raster matches the intended CRS, grid, and nodata."""
    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, expected_crs, src.width, src.height, *src.bounds
        )
        src_nodata = src.nodata

    with rasterio.open(dst_path) as dst:
        assert dst.crs == expected_crs, f"CRS mismatch: {dst.crs}"
        assert (dst.width, dst.height) == (width, height), "Grid size mismatch"
        assert dst.transform == transform, "Transform mismatch"
        assert dst.nodata == src_nodata, "Nodata not preserved"

    print(f"Verified: {dst_path} is {expected_crs}, {width}x{height}, nodata OK")


verify_reprojection(
    "landsat_ndvi_4326.tif", "landsat_ndvi_utm.tif", "EPSG:32633"
)

If the transform assertion fails, a band was reprojected against a stale grid — the exact failure that misaligns stacked features. If the CRS assertion fails, check that you passed a valid authority code; malformed strings surface as CRSError: Input is not a CRS from pyproj. Once every band passes, the raster is ready to feed the wider feature engineering workflow.

FAQ

Which resampling method should I use when reprojecting a raster?

Match the method to the data type. Use Resampling.nearest for categorical rasters such as land-cover or soil-class maps, because it copies existing integer codes without inventing new ones. Use Resampling.bilinear or Resampling.cubic for continuous data such as elevation, reflectance, or NDVI, where interpolated values are physically meaningful. Never use bilinear on class labels — averaging class 3 and class 7 into 5 produces a category that does not exist in your legend.

How do I keep nodata pixels from becoming valid data after reprojection?

Pass src_nodata and dst_nodata explicitly to reproject and set the same value in the output profile. rasterio then fills pixels outside the source footprint with the nodata value instead of zero, and stops interpolation from blending nodata into neighbouring valid pixels. Write the value into the profile with profile.update(nodata=...) so downstream readers and masks recognise it.

Why do my reprojected bands end up misaligned when I stack them for ML?

Misalignment happens when each band is reprojected independently and calculate_default_transform picks a slightly different width, height, or transform per file. Compute one target transform, width, and height once, then reproject every band onto that identical grid. Sharing a single destination transform guarantees pixel-for-pixel registration so the stacked array is a valid multi-band feature cube.

Part of: CRS Alignment and Projection Handling — Spatial Feature Engineering for Machine Learning