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.

Reprojecting a UTM grid to a WGS84 grid A source raster in UTM metres has a rectangular footprint. Under reprojection its corners map to a curved quadrilateral in WGS84 degrees. calculate_default_transform computes an axis-aligned bounding grid in degrees with a resolution that preserves the original ground sampling distance, and reproject resamples every pixel onto it. Source: UTM (metres) square pixels, EPSG:326xx reproject Warped footprint (degrees) dashed = default_transform grid Output: EPSG:4326 resampled to degree grid GSD preserved: 10 m ≈ 0.00009° — the numeric resolution changes, the ground footprint does not.

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

Package Minimum version Why required
rasterio 1.3.0 Provides warp.calculate_default_transform, warp.reproject, Resampling
GDAL (C library) 3.4.0 Runs the underlying warp and CRS operations
numpy 1.23 Backs the source and destination pixel arrays 
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.

What happens to the fill border during a warp A bilinear kernel averages a two by two neighbourhood. Without src_nodata, the fill value is treated as data, so a ring of blended garbage appears along every edge of the reprojected scene. Passing src_nodata excludes those cells from the average and dst_nodata fills the uncovered corners cleanly. no nodata passed src_nodata + dst_nodata passed valid pixels a 1–2 px halo of averaged fill along every edge valid pixels edge stays sharp; corners hold the declared fill value The halo survives every downstream step: it biases indices, statistics and any mosaic that blends these edges. 
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

Choosing the destination grid calculate_default_transform derives a grid that preserves ground sampling distance, which suits one-off conversions. Copying a reference dataset's transform guarantees pixel-for-pixel alignment for stacking or differencing. Specifying an explicit resolution produces a predictable tiling scheme, at the cost of resampling everything onto it. Three ways to decide where the output pixels land default transform grid derived from the source corners and its GSD best for: one-off conversion, publishing a single scene alignment with others: none match a reference copy crs, transform, width and height from another file best for: stacking, differencing, time series, band math alignment: exact, by construction fixed resolution pass resolution= to get a round degree or metre step best for: tile pyramids and fixed-grid archives alignment: exact within the scheme If two rasters will ever be subtracted, use the middle column — everything else invites a half-pixel argument later. Each warp resamples; chaining conversions compounds the smoothing, so reproject once from the original.

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.

Frequently Asked Questions

Q: Why does my WGS84 output have a different pixel size than the UTM input? WGS84 measures distance in degrees, not metres, so a 10 m UTM pixel becomes roughly 0.00009 degrees. calculate_default_transform derives an appropriate degree-based resolution automatically; the ground sampling distance is preserved even though the numeric value changes.

Q: Which resampling method should I use when reprojecting? Use nearest for categorical data such as land-cover classes or masks to avoid inventing new values. Use bilinear or cubic for continuous data such as reflectance or elevation. Never use bilinear on a class map.

Q: How do I stop nodata pixels from bleeding into valid data? Pass src_nodata and dst_nodata to reproject so the warper excludes fill pixels from interpolation and fills uncovered output cells with the nodata value. Also set nodata in the destination profile so downstream tools honour it.