Reprojecting a Raster from UTM to WGS84

To reproject a UTM GeoTIFF into WGS84 (EPSG:4326), compute the destination grid with calculate_default_transform, then warp the array with rasterio.warp.reproject — this preserves ground resolution and alignment instead of naively relabelling the CRS:

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

with rasterio.open("utm_scene.tif") as src:
    transform, width, height = calculate_default_transform(
        src.crs, "EPSG:4326", src.width, src.height, *src.bounds
    )
    profile = src.profile | {"crs": "EPSG:4326", "transform": transform,
                             "width": width, "height": height}
    with rasterio.open("wgs84_scene.tif", "w", **profile) as dst:
        reproject(rasterio.band(src, 1), rasterio.band(dst, 1),
                  resampling=Resampling.bilinear)

Getting the destination transform right is the crux of every CRS change, which is why it anchors Mastering CRS Transformations in rasterio.

Why This Arises in Remote Sensing Workflows

Sentinel-2, Landsat, and most analysis-ready satellite products ship in UTM because a metre-based projected grid keeps pixels square and distances undistorted within a zone. That is ideal for area and distance calculations, but it fragments a wide-area dataset across multiple UTM zones and makes overlaying data in a web map (which expects geographic or Web Mercator coordinates) awkward. Reprojecting to WGS84 gives a single global degree-based grid that every mapping library and STAC consumer understands.

The mistake to avoid is treating reprojection as a metadata edit. Overwriting src.crs alone leaves the pixel grid and transform describing UTM metres while claiming degrees, which throws the raster thousands of kilometres off position. A correct reprojection resamples the pixels onto a new grid: calculate_default_transform derives that grid — its origin, degree-based resolution, width, and height — and reproject moves every pixel value onto it. If your CRS is simply mislabelled rather than genuinely different, that is a separate problem covered in Fixing EPSG Mismatches in rasterio.open.

This topic sits inside Mastering CRS Transformations in rasterio and, more broadly, Core Raster Fundamentals & STAC Mapping.

What calculate_default_transform Actually Computes

The diagram shows why a UTM square becomes a subtly skewed quadrilateral in WGS84, and how the default transform bounds it into a clean axis-aligned grid.

calculate_default_transform inspects the source corners, projects them into the target CRS, takes the bounding box, and divides it into a grid whose cell size matches the original ground sampling distance. You get back a destination transform, width, and height that you drop straight into the output profile.

Environment & Setup

pip install "rasterio>=1.3.0" "numpy>=1.23"

An optional but recommended addition is rioxarray, used in the first variant below for a one-line labelled-array reprojection:

pip install rioxarray

Complete Working Example

This function reprojects every band of a UTM raster to WGS84, preserving nodata and choosing a resampling method per data class. It writes a tiled output so the result is itself efficient to read later.

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


def reproject_to_wgs84(
    src_path: str,
    dst_path: str,
    *,
    resampling: Resampling = Resampling.bilinear,
    dst_crs: str = "EPSG:4326",
) -> None:
    """
    Reproject a raster (e.g. UTM) to WGS84, preserving alignment and nodata.

    Parameters
    ----------
    src_path : str
        Source GeoTIFF in a projected CRS such as UTM.
    dst_path : str
        Output path for the reprojected GeoTIFF.
    resampling : Resampling
        Use nearest for categorical rasters, bilinear/cubic for continuous data.
    dst_crs : str
        Target CRS. Defaults to geographic WGS84.
    """
    with rasterio.open(src_path) as src:
        # Derive the destination grid: transform + dimensions in the target CRS
        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,
            driver="GTiff",
            tiled=True,          # write a tiled output for efficient later reads
            blockxsize=512,
            blockysize=512,
            compress="deflate",
        )

        with rasterio.open(dst_path, "w", **profile) as dst:
            # Warp each band independently onto the new grid
            for band in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, band),
                    destination=rasterio.band(dst, band),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    src_nodata=src.nodata,   # exclude fill pixels from interpolation
                    dst_nodata=src.nodata,   # fill uncovered output cells
                    resampling=resampling,
                    num_threads=4,           # parallelise the warp across cores
                )


if __name__ == "__main__":
    reproject_to_wgs84(
        "S2A_36NYF_20230615_B04.tif",
        "S2A_36NYF_20230615_B04_wgs84.tif",
        resampling=Resampling.bilinear,
    )
    with rasterio.open("S2A_36NYF_20230615_B04_wgs84.tif") as check:
        print("output CRS :", check.crs)
        print("output res :", check.res)      # now in degrees
        print("output size:", check.width, "x", check.height)

Two details matter for correctness. Passing rasterio.band(src, band) hands reproject a band handle so GDAL streams the warp rather than loading whole arrays, and passing both src_nodata and dst_nodata keeps fill values from smearing into valid pixels along the reprojected edges.

Variant Patterns

1. One-line reprojection with rioxarray

If you already work with labelled arrays, rioxarray wraps the whole calculate_default_transform + reproject dance in a single call and carries the CRS and nodata through automatically.

import rioxarray

# open_rasterio with masked=True promotes nodata to NaN
da = rioxarray.open_rasterio("utm_scene.tif", masked=True)

# Reproject to WGS84; resampling defaults to nearest, override for continuous data
from rasterio.enums import Resampling
da_wgs84 = da.rio.reproject("EPSG:4326", resampling=Resampling.bilinear)

da_wgs84.rio.to_raster("wgs84_scene.tif", compress="deflate", tiled=True)
print(da_wgs84.rio.crs, da_wgs84.rio.resolution())

This is convenient inside array-native workflows such as Band Math Operations with xarray, where you want the reprojected result to stay a labelled DataArray.

2. Reprojecting a categorical mask correctly

For a land-cover or cloud mask, interpolation would fabricate class values that never existed. Force nearest-neighbour resampling and keep the integer dtype.

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

with rasterio.open("landcover_utm.tif") as src:
    transform, width, height = calculate_default_transform(
        src.crs, "EPSG:4326", src.width, src.height, *src.bounds
    )
    profile = src.profile | {
        "crs": "EPSG:4326", "transform": transform,
        "width": width, "height": height,
    }
    with rasterio.open("landcover_wgs84.tif", "w", **profile) as dst:
        reproject(
            rasterio.band(src, 1), rasterio.band(dst, 1),
            src_crs=src.crs, dst_crs="EPSG:4326",
            src_transform=src.transform, dst_transform=transform,
            resampling=Resampling.nearest,   # never invent class values
            src_nodata=src.nodata, dst_nodata=src.nodata,
        )

3. Matching an existing target grid

When the output must align pixel-for-pixel with another dataset (for stacking or differencing), skip the default transform and supply the target grid explicitly.

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

with rasterio.open("reference_wgs84.tif") as ref:
    dst_crs, dst_transform = ref.crs, ref.transform
    dst_w, dst_h = ref.width, ref.height

with rasterio.open("utm_scene.tif") as src:
    dst_arr = np.full((dst_h, dst_w), src.nodata or 0, dtype=src.dtypes[0])
    reproject(
        source=rasterio.band(src, 1),
        destination=dst_arr,
        src_crs=src.crs, dst_crs=dst_crs,
        src_transform=src.transform, dst_transform=dst_transform,
        resampling=Resampling.bilinear,
        src_nodata=src.nodata, dst_nodata=src.nodata,
    )
    print("aligned to reference grid:", dst_arr.shape)

Common Errors

CRSError: Invalid projection: EPSG:4326

The installed PROJ database cannot resolve the code, usually because PROJ_LIB/PROJ_DATA points at a stale or missing directory after a manual GDAL install. Reinstall rasterio from wheels (pip install --force-reinstall rasterio) so it bundles a matching PROJ, or fix the PROJ_DATA environment variable.

Output is shifted or scaled by a huge factor

You reassigned crs without recomputing transform, so the pixels are still on the UTM grid but labelled as degrees. Always pair a CRS change with calculate_default_transform and a fresh transform, width, and height.

Nodata edges bleed into valid pixels

src_nodata/dst_nodata were not passed to reproject, so the resampler averaged fill values into real data along the warped border. Supply both, and set nodata in the destination profile so readers mask those cells.

Related

• Mastering CRS Transformations in rasterio — the parent section covering CRS objects, transforms, and coordinate operations end to end.
• Fixing EPSG Mismatches in rasterio.open — diagnose whether you need a reprojection or just a corrected CRS label.
• Band Math Operations with xarray — where the rioxarray reprojection variant fits into labelled-array pipelines.
• Core Raster Fundamentals & STAC Mapping — the broader context for grids, projections, and catalog alignment.