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:

import rasterio

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

with rasterio.open("out.tif") as bad:
    print(bad.crs, tuple(round(v, 2) for v in bad.bounds))

EPSG:27700 (-2.0, 53.6, -1.6, 53.85)

That file now claims to be British National Grid while holding degrees: a raster about half a metre wide, sitting a couple of metres from the grid's false origin in the sea south-west of the Isles of Scilly — nearly 600 km from Halifax, where the data belongs. It opens. It plots.

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
            )
with rasterio.open("elevation_bng.tif") as out:
    print(out.crs, out.shape, tuple(round(v) for v in out.bounds))

EPSG:27700 (957, 908) (400084, 411511, 426542, 439397)

What you have, the destination grid, a deliberate resolution

1. Establish what you actually have

with rasterio.open("elevation_wgs84.tif") as src:
    print(src.crs)
    print(src.res, "← degrees, and not square")
    print(tuple(round(v, 4) for v in src.bounds))
    print(src.dtypes, src.nodata)

EPSG:4326
(0.00041666666666666664, 0.0002777777777777778) ← degrees, and not square
(-2.0002, 53.6001, -1.6002, 53.8501)
('float32',) None

Two things to note before going further. The resolution is in degrees — and for Copernicus data north of 50° the cells are 1.5 arc-seconds wide and 1 tall — 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 cannot transform naive geometries; 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(f"pixel {transform.a:.2f} × {-transform.e:.2f} m; {width} × {height} cells, from {src.width} × {src.height}")

pixel 29.14 × 29.14 m; 908 × 957 cells, from 960 × 900

calculate_default_transform projects the source corners, takes the bounding box of the result, and picks a square pixel size that keeps roughly the same total cell count. 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.

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=30,           # exactly 30 m cells, the data's nominal resolution
)
print(width, height)

882 930

Reasons to override:

• Round numbers matter downstream. Cells of 23.something metres make every area calculation awkward; 30 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

import numpy as np

with rasterio.open("elevation_wgs84.tif") as src:
    transform, width, height = calculate_default_transform(src.crs, DST_CRS, src.width, src.height,
                                                           *src.bounds, resolution=30)
    for label, nodata in [("no dst_nodata", None), ("dst_nodata=-9999", -9999)]:
        out = np.full((height, width), -9999 if nodata else 0, dtype="float32")
        reproject(source=rasterio.band(src, 1), destination=out,
                  src_transform=src.transform, src_crs=src.crs, src_nodata=src.nodata,
                  dst_transform=transform, dst_crs=DST_CRS, dst_nodata=nodata,
                  resampling=Resampling.bilinear)
        print(f"{label:17} cells at exactly 0 m: {(out == 0).sum():>6,}; at -9999: {(out == -9999).sum():>6,}")

no dst_nodata     cells at exactly 0 m:  4,432; at -9999:      0
dst_nodata=-9999  cells at exactly 0 m:      0; at -9999:  4,432

Without src_nodata, a sentinel value in the source 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. (This Copernicus window has no NoData cells, so there is nothing to smear here.) Without dst_nodata, the corners of the rotated footprint are filled with 0, which on an elevation raster is indistinguishable from sea level — and the Pennines have no sea.

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

with rasterio.open("elevation_wgs84.tif") as src:
    transform, width, height = calculate_default_transform(src.crs, DST_CRS, src.width, src.height,
                                                           *src.bounds, resolution=30)
    profile = src.profile | {"crs": DST_CRS, "transform": transform, "width": width, "height": height,
                             "nodata": -9999}
    with rasterio.open("elevation_bng.tif", "w", **profile) as dst:
        reproject(rasterio.band(src, 1), rasterio.band(dst, 1), dst_nodata=-9999, resampling=Resampling.bilinear)

check("elevation_wgs84.tif", "elevation_bng.tif")

crs        EPSG:4326 → EPSG:27700
shape      (900, 960) → (930, 882)
res        (0.00041666666666666664, 0.0002777777777777778) → (30.0, 30.0)
value min  29.0 → 29.022972106933594
value max  480.2169494628906 → 480.13629150390625
mean       189.96 → 189.98
valid %    100.0 → 99.5

Read that output carefully — it is the whole test:

• min and max barely moved — bilinear interpolation averages, so it can pull the extremes in slightly but never beyond them.

• mean barely moved — the grid changed, the data did not.

• valid dropped by half a per cent — the NoData corners of the rotated footprint. Expected. A drop to 60% would mean something else went wrong.

A reusable function with the right defaults

The data used on this page

Two real rasters for the same patch of the South Pennines around Halifax and Huddersfield (2.0–1.6°W, 53.6–53.85°N), both in WGS84 as their publishers distribute them:

• elevation_wgs84.tif — a window of the Copernicus 30 m DEM (GLO-30), float32, read from its public bucket.
• landcover_wgs84.tif — the same window of ESA WorldCover 2021, 10 m, uint8 class codes.

raw/ holds copies of both plus scan_1987.tif, a block of pixels written with no georeferencing at all, standing in for a scanned map sheet.

import shutil
from pathlib import Path

import numpy as np
import rasterio
from rasterio.windows import from_bounds

BOUNDS = (-2.0, 53.6, -1.6, 53.85)
SOURCES = {
    "elevation_wgs84.tif": "https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_W002_00_DEM/"
                           "Copernicus_DSM_COG_10_N53_00_W002_00_DEM.tif",
    "landcover_wgs84.tif": "https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/"
                           "ESA_WorldCover_10m_2021_v200_N51W003_Map.tif",
}
for name, url in SOURCES.items():
    with rasterio.open(url) as src:
        window = from_bounds(*BOUNDS, transform=src.transform).round_offsets().round_lengths()
        profile = src.profile | {"driver": "GTiff", "width": int(window.width), "height": int(window.height),
                                 "transform": src.window_transform(window), "compress": "deflate"}
        with rasterio.open(name, "w", **profile) as dst:
            dst.write(src.read(window=window))

Path("raw").mkdir(exist_ok=True)
shutil.copy("elevation_wgs84.tif", "raw/elevation.tif")
shutil.copy("landcover_wgs84.tif", "raw/landcover.tif")
with rasterio.open("elevation_wgs84.tif") as src:
    block = src.read(1)[:500, :500]
with rasterio.open("raw/scan_1987.tif", "w", driver="GTiff", width=500, height=500, count=1, dtype="float32") as dst:
    dst.write(block, 1)
for name in SOURCES:
    with rasterio.open(name) as src:
        print(f"{name:20} {src.width:>5} × {src.height:<5} {src.dtypes[0]:8} nodata {src.nodata}")

elevation_wgs84.tif    960 × 900   float32  nodata None
landcover_wgs84.tif   4800 × 3000  uint8    nodata 0.0

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


reproject_raster("landcover_wgs84.tif", "landcover_bng.tif", "EPSG:27700", resolution=10)
with rasterio.open("landcover_wgs84.tif") as a, rasterio.open("landcover_bng.tif") as b:
    print(f"classes before: {sorted(np.unique(a.read(1)).tolist())}")

    print(f"classes after:  {sorted(np.unique(b.read(1)).tolist())}   (0 is WorldCover's NoData, filling the new corners)")

classes before: [10, 20, 30, 40, 50, 60, 80, 90]
classes after:  [0, 10, 20, 30, 40, 50, 60, 80, 90]   (0 is WorldCover's NoData, filling the new corners)

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("elevation_wgs84.tif", "landcover_bng.tif", "elevation_on_landcover_grid.tif", Resampling.bilinear)

NAMES = {10: "tree cover", 20: "shrubland", 30: "grassland", 40: "cropland", 50: "built-up",
         60: "bare", 80: "water", 90: "wetland", 100: "moss and lichen"}
with rasterio.open("elevation_on_landcover_grid.tif") as e, rasterio.open("landcover_bng.tif") as c:
    assert e.transform == c.transform and e.shape == c.shape
    elev, cover = e.read(1), c.read(1)
for code in sorted(np.unique(cover[cover > 0]).tolist()):
    sel = cover == code
    print(f"{NAMES.get(code, code):>11}: {sel.sum():>9,} cells, mean elevation {elev[sel].mean():5.0f} m")

 tree cover: 2,051,980 cells, mean elevation   152 m
  shrubland:       552 cells, mean elevation   368 m
  grassland: 3,427,350 cells, mean elevation   232 m
   cropland:   200,069 cells, mean elevation   156 m
   built-up: 1,617,459 cells, mean elevation   149 m
       bare:     6,716 cells, mean elevation   166 m
      water:    36,710 cells, mean elevation   248 m
    wetland:       660 cells, mean elevation   422 m

Afterwards the two rasters share transform, CRS and shape exactly, so elev[cover == 50].mean() — the mean elevation of built-up land — 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

results = reproject_folder("raw/", "bng/", "EPSG:27700", resolution=30)
print(next(r["error"] for r in results if r["status"] == "failed"))

  ✓ elevation.tif                bilinear
  ✓ landcover.tif                nearest
  ✗ scan_1987.tif                bilinear
raw/scan_1987.tif has no CRS — assign one before reprojecting

scan_1987.tif is the setup's stand-in for a scanned map sheet — pixels with no georeferencing — and the function refuses it with a reason instead of inventing a location. 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.

Why rasters cannot simply change CRS

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.

Relabelling without warping and the antimeridian

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

Internal links

• Raster resampling explained — the choice this operation forces on you
• The raster data model explained — transform, CRS, dtype and NoData
• How to reproject spatial data in Python (GeoPandas) — the vector equivalent, and how it differs
• Raster and vector do not line up in Python — diagnosing the result when this goes wrong
• How to choose the right projected CRS for your study area — deciding what to reproject to
• How to standardise and repair CRS across a folder of files — the vector-side folder sweep
• How to merge and mosaic rasters in Python — the next step after aligning several tiles
• How to reproject between datums correctly — when metres of accuracy matter

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.

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.