Advanced Resampling and Upscaling Techniques

Resampling and spatial upscaling are critical inflection points in any satellite processing pipeline. When moving from raw sensor acquisitions to analysis-ready data, the choice of interpolation strategy directly shapes spectral fidelity, spatial accuracy, and downstream model performance. For the full architectural context of where this step sits, see Satellite Processing Workflows & Index Pipelines.

These operations typically occur after radiometric calibration and atmospheric correction. Analysts and data engineers must balance computational overhead against scientific rigor, particularly when harmonizing multi-sensor datasets or preparing inputs for machine learning architectures. Improper interpolation introduces aliasing, shifts feature boundaries, and artificially inflates spectral variance — ultimately compromising classification accuracy and time-series consistency.

Prerequisites

Install dependencies via conda to ensure binary compatibility with GDAL:

conda install -c conda-forge rasterio rioxarray numpy gdal dask

Conceptual prerequisites: familiarity with affine geotransforms and pixel resolution and an understanding of the CRS transformation model in rasterio will make the code patterns below immediately legible.

Resampling vs. Upscaling: Distinct Operations

While frequently conflated in operational documentation, resampling and upscaling serve distinct mathematical purposes:

Resampling adjusts pixel alignment, grid spacing, or coordinate reference systems without fundamentally altering information density. The source and destination grids differ in origin, spacing, or projection — but the new pixel values are derived directly from existing measurements. Common methods: nearest-neighbor (categorical layers), bilinear and bicubic (continuous reflectance), area-weighted averaging (aggregation or downscaling).

Upscaling increases spatial resolution by synthesizing new pixel values through interpolation kernels, statistical modeling, or deep-learning super-resolution. This introduces estimated information at sub-original-resolution scales and demands rigorous validation against higher-resolution reference imagery or ground truth. The Sentinel-2 10 m / 20 m / 60 m harmonization problem is a canonical upscaling task: see Choosing the right resampling method for Sentinel-2 for sensor-specific kernel guidance.

Algorithm selection depends on data type, target resolution, and downstream analytical objectives. For multispectral optical imagery, preserving spectral integrity while minimizing high-frequency aliasing is paramount.

Step-by-Step Workflow

A robust resampling pipeline follows a deterministic sequence to prevent coordinate drift, spectral distortion, and memory exhaustion.

Step 1 — Validate Metadata and Harmonize CRS

Verify CRS consistency across all input bands before touching any pixel values. Mismatched projections or undefined datums silently corrupt affine transformations and produce geographically displaced outputs.

import rasterio
from rasterio.crs import CRS

def validate_crs(path: str, expected_epsg: int = 32633) -> None:
    """Raise ValueError if the raster CRS does not match expected_epsg."""
    with rasterio.open(path) as src:
        if src.crs is None:
            raise ValueError(f"{path}: no CRS defined — cannot proceed")
        if not src.crs.equals(CRS.from_epsg(expected_epsg)):
            raise ValueError(
                f"{path}: CRS is {src.crs.to_epsg()}, expected EPSG:{expected_epsg}"
            )
        if src.nodata is None:
            # Log a warning but do not abort — nodata can be set at write time
            print(f"Warning: {path} has no nodata value defined")
        print(f"OK  {path}  CRS=EPSG:{expected_epsg}  res={src.res}  dtype={src.dtypes[0]}")

Strip invalid metadata tags with gdal_edit.py -unsetmd for tags that confuse downstream parsers. Standardize nodata values to a single sentinel (commonly 0 for integer imagery or nan for float32) before reprojection — mismatched nodata conventions at read time cause ghost pixels at tile edges.

Step 2 — Compute Target Grid and Affine Transform

Calculate the destination affine transform based on desired pixel dimensions and spatial extent. Always use rasterio.warp.calculate_default_transform to derive a numerically stable target transform rather than computing the affine scale manually — manual arithmetic accumulates floating-point error across large extents.

from rasterio.warp import calculate_default_transform

def compute_target_transform(
    src_path: str,
    target_res_m: float
) -> tuple:
    """
    Return (transform, out_width, out_height) for a given target resolution.
    Rounds dimensions to the nearest integer to avoid sub-pixel boundary drift.
    """
    with rasterio.open(src_path) as src:
        scale = src.res[0] / target_res_m
        out_width  = round(src.width  * scale)
        out_height = round(src.height * scale)

        transform, _, _ = calculate_default_transform(
            src.crs, src.crs,
            out_width, out_height,
            left=src.bounds.left,
            bottom=src.bounds.bottom,
            right=src.bounds.right,
            top=src.bounds.top,
        )
    return transform, out_width, out_height

When aligning multiple scenes to a common analysis grid (e.g. a fixed MGRS tile extent), snap the target extent to that grid’s pixel boundaries — even a 0.5 m sub-pixel shift compounds into visible seam artifacts during seamless mosaicking.

Step 3 — Execute Resampling with Windowed Chunking

Large rasters cannot be loaded entirely into RAM. The pattern below processes data in spatially contiguous blocks, writing each block directly to the output file. This approach keeps peak memory proportional to block size, not scene size.

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

def resample_raster(
    input_path: str,
    output_path: str,
    target_res_m: float,
    method: Resampling = Resampling.bilinear,
    block_size: int = 1024,
) -> None:
    """
    Production-grade windowed resampling with explicit transform computation
    and lossless LZW-compressed COG-compatible output.
    """
    transform, out_width, out_height = compute_target_transform(input_path, target_res_m)

    with rasterio.open(input_path) as src:
        profile = src.profile.copy()
        profile.update(
            driver="GTiff",
            width=out_width,
            height=out_height,
            transform=transform,
            compress="lzw",
            tiled=True,
            blockxsize=block_size,
            blockysize=block_size,
            interleave="band",
        )

        with rasterio.open(output_path, "w", **profile) as dst:
            for band_idx in range(1, src.count + 1):
                destination = np.empty(
                    (out_height, out_width),
                    dtype=src.dtypes[band_idx - 1],
                )
                reproject(
                    source=rasterio.band(src, band_idx),
                    destination=destination,
                    src_transform=src.transform,
                    dst_transform=transform,
                    src_crs=src.crs,
                    dst_crs=src.crs,
                    resampling=method,
                    src_nodata=src.nodata,
                    dst_nodata=src.nodata,
                )
                dst.write(destination, band_idx)

Step 4 — Select the Correct Kernel per Band Type

Kernel choice has a larger impact on output quality than any other single parameter. The following table covers the most common cases in multispectral satellite workflows:

For spectral index calculation inputs that require all bands on a common 10 m grid, bilinear is the standard choice for Sentinel-2 red-edge and SWIR bands. When the output feeds a cloud and shadow masking step that relies on a QA/SCL band, always resample the mask layer with nearest-neighbor on its own pass.

Step 5 — Dask-Backed Lazy Resampling with rioxarray

For petabyte-scale archives or distributed compute environments, rioxarray with dask defers all computation until write time, enabling out-of-core streaming without explicit windowing code:

import rioxarray  # noqa: F401 — activates .rio accessor on xarray objects
import xarray as xr

def lazy_resample(
    input_path: str,
    output_path: str,
    target_res_m: float,
    method: str = "bilinear",
    chunk_xy: int = 1024,
) -> None:
    """
    Dask-backed lazy resampling via rioxarray.reproject.
    Computation is deferred until .to_raster() triggers the write.
    """
    da = xr.open_dataarray(
        input_path,
        engine="rasterio",
        chunks={"x": chunk_xy, "y": chunk_xy},
        mask_and_scale=False,  # keep raw integer reflectance; do not auto-scale
    )

    resampled = da.rio.reproject(
        da.rio.crs,
        resolution=target_res_m,
        resampling=method,
        nodata=da.rio.nodata,
    )

    resampled.rio.to_raster(
        output_path,
        driver="GTiff",
        compress="lzw",
        tiled=True,
        windowed=True,  # write block-by-block to limit RAM
    )

Pinning mask_and_scale=False is essential for integer reflectance data (DN or scaled reflectance). With auto-scaling enabled, xarray converts integer bands to float64 silently, tripling memory use and producing floating-point nodata mismatch errors at the write step.

Parameter Reference

Verification and Testing

After resampling, confirm that the output geometry and values are correct before the data enters downstream steps.

def verify_resampled_output(
    output_path: str,
    expected_res_m: float,
    original_path: str,
    tolerance_px: float = 0.01,
) -> None:
    """
    Assert spatial resolution, extent, and basic spectral consistency.
    Raises AssertionError with a descriptive message on any failure.
    """
    with rasterio.open(output_path) as dst, rasterio.open(original_path) as src:
        # 1. Resolution check
        actual_res = dst.res[0]
        assert abs(actual_res - expected_res_m) < tolerance_px, (
            f"Resolution mismatch: got {actual_res:.4f} m, expected {expected_res_m} m"
        )
        # 2. Extent preservation: bounds should match within one source pixel
        for attr in ("left", "bottom", "right", "top"):
            diff = abs(getattr(dst.bounds, attr) - getattr(src.bounds, attr))
            assert diff < src.res[0], (
                f"Bounds drift on {attr}: {diff:.4f} m (> one source pixel)"
            )
        # 3. Nodata propagation
        assert dst.nodata == src.nodata, (
            f"Nodata changed: {src.nodata} → {dst.nodata}"
        )
        # 4. Band count preserved
        assert dst.count == src.count, (
            f"Band count changed: {src.count} → {dst.count}"
        )
        print(
            f"PASS  {output_path}  res={actual_res} m  "
            f"size={dst.width}x{dst.height}  bands={dst.count}"
        )

For continuous-band outputs, run a Structural Similarity Index (SSIM) comparison against a downsampled reference to quantify spatial coherence. For categorical mask outputs, compare class histograms pre- and post-resampling — any class that was absent in the source should remain absent in the output.

Misaligned grids after resampling will cause feature displacement during automated image clipping and cropping, producing contaminated training samples and incorrect spectral index values at parcel boundaries.

Troubleshooting

PROJ_NETWORK errors or silent CRS fallback to WGS84. Cause: PROJ cannot reach the CDN to download datum shift grids. Fix: set PROJ_NETWORK=OFF and ensure the local PROJ data directory contains the required .tif grid files (conda install -c conda-forge proj-data), or pre-download with projsync.

Sub-pixel coordinate drift accumulating across a mosaic. Cause: target extent coordinates are not snapped to the reference grid’s pixel boundary. Fix: round all extent coordinates to n * target_res_m before passing to calculate_default_transform. A 0.3 m drift per tile becomes 3 m over 10 adjacent tiles.

Out-of-memory crash on a large scene. Cause: block_size too large for available RAM, or dask scheduler spawning too many concurrent tasks. Fix: reduce block_size to 512, limit dask workers with dask.config.set(num_workers=2), or stream via rasterio’s window API rather than computing the full array in memory.

Bilinear resampling introduces ringing artifacts at sharp land/water boundaries. Cause: the bilinear kernel does not respect discontinuities — it interpolates across the DN jump. Fix: apply a pre-resampling median filter to suppress the artifact, or resample the land/water mask separately with nearest-neighbor and use it to zero out boundary pixels post-resampling.

Output dtype silently promoted to float64. Cause: rioxarray with mask_and_scale=True (the default) applies CF conventions and promotes integer arrays. Fix: open with mask_and_scale=False and set an explicit nodata value before calling .rio.reproject().

Frequently Asked Questions

Which resampling method should I use for Sentinel-2 reflectance bands? Use bilinear for 20 m → 10 m upscaling of continuous reflectance bands (B05–B07, B8A). Reserve nearest for the SCL classification band and binary cloud masks. Use average when downscaling to coarser grids such as 60 m → 100 m for comparison with MODIS.

How do I avoid coordinate drift when resampling multiple scenes to a common grid? Snap the target extent to the nearest pixel boundary of the reference grid before computing the affine transform. Use rasterio.warp.calculate_default_transform rather than computing scale manually — manual arithmetic accumulates floating-point error across large extents.

Can I resample a raster without loading it entirely into RAM? Yes. Use rasterio’s windowed read/write API to process one block at a time, or use rioxarray with a dask backend to stream data lazily. Both approaches cap peak memory at roughly one chunk (e.g. 1024×1024 pixels per band).

---

Related

• Choosing the Right Resampling Method for Sentinel-2 — Sensor-specific kernel recommendations for all Sentinel-2 spectral bands and spatial resolutions.
• Seamless Mosaicking and Edge Blending — Grid-aligned tiles from this workflow feed directly into feathered multi-scene composites.
• Automated Image Clipping and Cropping — Clipping to irregular polygon boundaries requires pixel-perfect alignment with the resampled grid.
• Handling Pixel Resolution and Scaling — Core explanation of affine transforms, pixel resolution, and scale factors underlying every resampling operation.
• Spectral Index Calculation Pipelines — Accurate NDVI, EVI, and NDWI computation depends on all bands sharing the same resampled grid.