⚙️ How-To Raster & rasterio CRS & projections Raster Analysis · step 4 of 14 Updated

How to Reproject a Raster in Python with Rasterio

Problem statement

Your vector layer is in British National Grid and your raster is in WGS 84. GeoPandas reprojects a vector in one line:

gdf = gdf.to_crs(27700)

Rasterio has no equivalent. There is no src.to_crs(), because reprojecting a raster is not a transformation of coordinates — it is the construction of an entirely new grid whose cells do not correspond one-to-one with the old ones.

Get it wrong and the failures are quiet:

# a common first attempt
profile = src.profile
profile.update(crs="EPSG:27700")          # ← relabels, does not reproject
with rasterio.open("out.tif", "w", **profile) as dst:
    dst.write(src.read())

That file now claims to be British National Grid while holding degrees. It opens. It plots. It is roughly 5,000 km from where it should be.

Quick answer 

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

DST_CRS = "EPSG:27700"

with rasterio.open("elevation_wgs84.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("elevation_bng.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,   # nearest for categorical data
            )
Vertical steps: compute the destination grid, update the profile, warp each band, verify bounds and values.
Four steps. Skipping step one is what produces a relabelled raster instead of a reprojected one.
Line Why it is there
calculate_default_transform works out the new grid's size, origin and pixel size
profile.update(...) width and height change — copying the old ones truncates the output
rasterio.band(src, i) streams band by band instead of holding both rasters in memory
resampling= always set explicitly; see resampling explained

Step-by-step solution

1. Establish what you actually have 

import rasterio

with rasterio.open("elevation_wgs84.tif") as src:
    print(src.crs)          # EPSG:4326
    print(src.res)          # (0.000277777, 0.000277777)  ← degrees
    print(src.bounds)       # (-2.7, 53.3, -1.9, 53.9)
    print(src.dtypes)       # ('int16',)
    print(src.nodata)       # -9999.0

Two things to note before going further. The resolution is in degrees, so the output resolution in metres is not something you can carry across. And src.crs is None means there is nothing to reproject from — see CRS not found; assigning a CRS is a different operation from converting one.

2. Compute the destination grid 

from rasterio.warp import calculate_default_transform

transform, width, height = calculate_default_transform(
    src.crs, "EPSG:27700", src.width, src.height, *src.bounds
)
print(transform)
print(width, height)
| 20.83, 0.00, 353000.65|
| 0.00,-20.83, 399000.19|
| 0.00, 0.00, 1.00|
2823 3210

calculate_default_transform projects the source corners, takes the bounding box of the result, and picks a pixel size that keeps roughly the same total cell count. The output is 2823 × 3210 where the input was 2880 × 2160 — the shape changed, which is why the profile must be updated rather than reused.

A projected raster's footprint is not a rectangle in the source CRS, so the output grid is the bounding box of a slightly curved quadrilateral. The corners of that box fall outside the source data and become NoData. This is normal and expected.

Two panels: a rectangular grid in the source CRS becomes a slightly curved quadrilateral inside a larger rectangular output grid.
The data is a quadrilateral; the grid must be a rectangle. The difference is NoData.

3. Choose the resolution deliberately

The default resolution is a guess. Override it when you have a reason:

transform, width, height = calculate_default_transform(
    src.crs, "EPSG:27700", src.width, src.height, *src.bounds,
    resolution=25,           # exactly 25 m cells
)
print(width, height)         # 2352 2674

Reasons to override:

  • Round numbers matter downstream. 20.83 m cells make every area calculation awkward; 25 m cells do not.
  • You are aligning to another raster. Then take the grid from that raster rather than computing one — see Example 2.
  • The default over-samples. Reprojecting from geographic to projected often produces slightly finer cells than the source really supports, inflating the file with no extra information.

4. Set the resampling method for what the data means 

from rasterio.warp import Resampling

Resampling.bilinear    # elevation, temperature, reflectance
Resampling.nearest     # land cover, soil class, any coded raster
Resampling.average     # continuous data being made much coarser
Resampling.mode        # categorical data being made much coarser

There is no safe default here. Interpolating class codes produces classes that do not exist, and the output rarely looks wrong enough to notice. The full decision is in raster resampling explained.

5. Handle NoData explicitly 

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,
    src_nodata=src.nodata,
    dst_nodata=src.nodata,
    resampling=Resampling.bilinear,
)

Without src_nodata, the sentinel value is treated as data and interpolated with its neighbours. A -9999 cell next to a 210 cell produces a smear of impossible intermediate elevations along every NoData boundary. Without dst_nodata, the corner fill defaults to 0, which on an elevation raster is indistinguishable from sea level.

6. Verify before trusting 

def check(src_path, dst_path):
    with rasterio.open(src_path) as s, rasterio.open(dst_path) as d:
        a, b = s.read(1, masked=True), d.read(1, masked=True)
        print(f"crs        {s.crs} → {d.crs}")
        print(f"shape      {s.shape} → {d.shape}")
        print(f"res        {s.res} → {d.res}")
        print(f"value min  {a.min()} → {b.min()}")
        print(f"value max  {a.max()} → {b.max()}")
        print(f"mean       {a.mean():.2f} → {b.mean():.2f}")
        print(f"valid %    {100*(1-a.mask.mean()):.1f} → {100*(1-b.mask.mean()):.1f}")

check("elevation_wgs84.tif", "elevation_bng.tif")
crs        EPSG:4326 → EPSG:27700
shape      (2160, 2880) → (2674, 2352)
res        (0.000278, 0.000278) → (25.0, 25.0)
value min  -3 → -3
value max  1344 → 1344
mean       214.83 → 214.81
valid %    100.0 → 92.6

Read that output carefully — it is the whole test:

  • min and max unchanged — the interpolation did not overshoot or clip.
  • mean barely moved — the grid changed, the data did not.
  • valid dropped to 92.6% — the NoData corners of the rotated footprint. Expected. A drop to 60% would mean something else went wrong.

Code examples

Example 1: a reusable function with the right defaults 

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

def reproject_raster(src_path, dst_path, dst_crs, *,
                     resampling=Resampling.nearest, resolution=None,
                     compress="deflate"):
    """Reproject a raster. `resampling` defaults to nearest — the only safe default."""
    src_path, dst_path = Path(src_path), Path(dst_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")

        kw = {"resolution": resolution} if resolution else {}
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds, **kw
        )
        profile = src.profile.copy()
        profile.update(crs=dst_crs, transform=transform,
                       width=width, height=height,
                       compress=compress, tiled=True,
                       blockxsize=256, blockysize=256)

        tmp = dst_path.with_suffix(".tmp.tif")
        with rasterio.open(tmp, "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=src.nodata,
                    resampling=resampling,
                    num_threads=4,
                )
        tmp.replace(dst_path)          # atomic: no half-written raster on failure
    return dst_path

Three details worth keeping. resampling defaults to nearest so a forgotten argument cannot corrupt categorical data. The write goes to a temporary path and is renamed only on success, so an interrupted run leaves no plausible-looking half-file — the same reasoning as in building a resumable batch job. And num_threads=4 is close to free: GDAL's warper parallelises well.

Example 2: reprojecting onto an existing raster's exact grid

Reprojecting two rasters to the same CRS does not make them cell-aligned. Each gets its own default grid, offset by an arbitrary fraction of a pixel. To compare them cell for cell, warp one onto the other's grid:

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

def match_grid(src_path, template_path, dst_path, resampling=Resampling.nearest):
    with rasterio.open(template_path) as tmpl:
        dst_transform = tmpl.transform
        dst_crs = tmpl.crs
        dst_shape = (tmpl.height, tmpl.width)

    with rasterio.open(src_path) as src:
        profile = src.profile.copy()
        profile.update(crs=dst_crs, transform=dst_transform,
                       width=dst_shape[1], height=dst_shape[0])
        out = np.empty((src.count, *dst_shape), dtype=src.dtypes[0])
        reproject(
            source=src.read(), destination=out,
            src_transform=src.transform, src_crs=src.crs, src_nodata=src.nodata,
            dst_transform=dst_transform, dst_crs=dst_crs, dst_nodata=src.nodata,
            resampling=resampling,
        )

    with rasterio.open(dst_path, "w", **profile) as dst:
        dst.write(out)

match_grid("slope_wgs84.tif", "landcover_bng.tif", "slope_aligned.tif",
           Resampling.bilinear)

Afterwards slope_aligned.tif and landcover_bng.tif share transform, CRS and shape exactly, so slope[cover == 3].mean() is meaningful. This is the standard preparation step before any raster algebra, and the reason raster and vector do not line up so often.

Example 3: reprojecting a folder, with per-file method selection 

from pathlib import Path
import numpy as np
import rasterio
from rasterio.warp import Resampling

CATEGORICAL = {"landcover", "soil", "zones", "mask"}

def method_for(path):
    """Categorical by name, or by the shape of the value distribution."""
    if any(word in path.stem.lower() for word in CATEGORICAL):
        return Resampling.nearest
    with rasterio.open(path) as src:
        if src.dtypes[0].startswith("float"):
            return Resampling.bilinear
        sample = src.read(1, out_shape=(256, 256), masked=True).compressed()
        return Resampling.nearest if len(np.unique(sample)) < 32 else Resampling.bilinear

def reproject_folder(src_dir, dst_dir, dst_crs, resolution=None):
    src_dir, dst_dir = Path(src_dir), Path(dst_dir)
    dst_dir.mkdir(parents=True, exist_ok=True)
    results = []
    for path in sorted(src_dir.glob("*.tif")):
        method = method_for(path)
        try:
            reproject_raster(path, dst_dir / path.name, dst_crs,
                             resampling=method, resolution=resolution)
            results.append({"file": path.name, "method": method.name, "status": "ok"})
        except Exception as exc:
            results.append({"file": path.name, "method": method.name,
                            "status": "failed", "error": str(exc)})
    for r in results:
        mark = "✓" if r["status"] == "ok" else "✗"
        print(f"  {mark} {r['file']:<28} {r['method']}")
    return results

reproject_folder("raw/", "bng/", "EPSG:27700", resolution=25)
  ✓ elevation.tif               bilinear
  ✓ landcover.tif               nearest
  ✓ ndvi_2026.tif               bilinear
  ✗ scan_1987.tif               nearest

The heuristic in method_for is a convenience, not an authority: it samples a 256 × 256 overview and calls anything with fewer than 32 distinct values categorical. Log the method it picked — as the output above does — so a wrong guess is visible rather than silent. Per-file error handling follows the pattern in batch script stops at the first bad file.

Explanation

Panels contrasting vector reprojection, which moves vertices, with raster reprojection, which builds a new grid and resamples.
A vector keeps its structure and moves. A raster gets a whole new structure.

The reason there is no raster.to_crs() is that reprojecting a vector and reprojecting a raster are different kinds of operation.

A vector geometry is a list of coordinates. Reprojecting it transforms each coordinate and leaves the structure untouched: same number of vertices, same rings, same topology. The output is the same object in a different reference frame.

A raster is a grid whose cells are defined by the coordinate system. Project the four corners of a rectangular grid into a different CRS and you get a quadrilateral with slightly curved edges — not a rectangle, and therefore not a grid. Something has to give, and what gives is the cells: the output is a new rectangular grid, and every one of its cells must be assigned a value derived from source cells that do not line up with it. That derivation is resampling, and it is why resampling= is a required part of the operation rather than an optimisation.

This has three consequences that explain most raster reprojection problems.

The output is bigger than the input. The bounding box of a curved quadrilateral is larger than the quadrilateral, so the corners fill with NoData. Reprojecting back and forth grows the footprint each time.

Values change. Except with nearest, output values are computed from several source values. Reprojecting an elevation raster twice and comparing to the original will not give identical numbers, and the difference is real information loss, not rounding.

Alignment is not implied by a shared CRS. Two rasters reprojected independently to EPSG:27700 have grids offset by an arbitrary fraction of a cell, because each one's origin came from its own bounding box. Cell-for-cell comparison requires warping onto a shared transform, as in Example 2. This is the single most common reason raster algebra produces subtly wrong answers — the arrays have the same shape, numpy happily multiplies them, and each cell is comparing slightly different ground.

Finally, note what calculate_default_transform is not doing: it is not preserving area, cell count or resolution exactly. It aims for a comparable total pixel count and a north-up grid. When the numbers matter — a resolution you will publish, a grid you will align to — supply them rather than accepting the default.

Edge cases or notes

  • Updating crs in the profile without warping relabels the file. It is the raster equivalent of set_crs versus to_crs — see set_crs vs to_crs in GeoPandas.
  • Reprojecting across the antimeridian produces a grid spanning the whole globe. Split the raster at 180° first, or use a CRS centred on the data.
  • Datum shifts need grid files. OSGB36 ↔ WGS 84 without the OSTN15 grid is out by up to 5 m. pyproj downloads grids on demand if PROJ_NETWORK=ON; see reprojecting between datums correctly.
  • num_threads speeds up the warp but not the read or the write. On a compressed source, decompression often dominates.
  • Very large rasters should be warped with rio warp or gdalwarp, which stream, rather than by reading whole bands into memory.
  • dst_nodata defaults to 0 if you do not set it and the profile has no nodata, which silently turns the fill corners into valid data.
  • profile.copy() matters. src.profile is a live view; mutating it without copying can affect the open dataset.
  • Overviews are not reprojected. Rebuild them on the output with rio overview --build 2,4,8,16.

FAQ

Why is there no to_crs() for rasters?

Because reprojecting a raster builds a new grid rather than moving existing coordinates. Every output cell needs a value derived from source cells that do not align with it, so the operation needs a destination grid and a resampling method — neither of which can be inferred.

Can I just change the CRS in the profile?

No. That relabels the coordinates without moving them, putting the raster in the wrong place while claiming otherwise. Use reproject.

Why is my output larger than my input?

The projected footprint is not a rectangle, and the output grid is its bounding box. The corners are NoData. This is expected; it grows each time you reproject.

How do I make two rasters line up exactly?

Warp both onto the same explicit transform, or warp one onto the other's grid as in Example 2. Reprojecting both to the same CRS independently does not align them.

Which resampling method should I use?

nearest for categorical rasters, bilinear for continuous ones, average or mode when making the raster much coarser. Never rely on the default.

My reprojected raster is full of −9999 smears along the edges.

src_nodata was not passed, so the sentinel was interpolated with real values. Pass both src_nodata and dst_nodata to reproject.

How do I reproject a raster too large to fit in memory?

Use rio warp or gdalwarp on the command line — both stream through the file. Alternatively warp window by window, but the command-line tools already do this well.