Removing Seams in Multi-Scene Mosaics with Feathering

Feathering removes seams from multi-scene mosaics by computing a Euclidean distance transform on each tile’s valid-pixel mask, then blending overlapping regions as a distance-weighted average: V_final = Σ(V_i × α_i) / Σ(α_i). Pixels deep inside a tile receive high weights; edge pixels receive near-zero weights, dissolving radiometric discontinuities without manual seamline digitization. This page belongs to the Seamless Mosaicking and Edge Blending workflow within the broader Satellite Processing Workflows & Index Pipelines.

The minimal snippet — two spatially aligned tiles, one blended output:

import rasterio
import numpy as np
from scipy.ndimage import distance_transform_edt

def feather_two_tiles(path_a: str, path_b: str, out_path: str, nodata: float = -9999.0) -> None:
    num, den, profile = None, None, None
    for path in (path_a, path_b):
        with rasterio.open(path) as src:
            if profile is None:
                profile = src.profile.copy()
                profile.update(dtype="float32", nodata=np.nan)
                num = np.zeros((src.count, src.height, src.width), dtype=np.float32)
                den = np.zeros((src.height, src.width), dtype=np.float32)
            data = src.read().astype(np.float32)
            valid = np.all(data != nodata, axis=0)
            dist = distance_transform_edt(valid)
            alpha = dist / (dist.max() + 1e-6)
            num += data * alpha[np.newaxis]
            den += alpha
    with np.errstate(divide="ignore", invalid="ignore"):
        result = np.where(den[np.newaxis] > 0, num / den[np.newaxis], np.nan)
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(result.astype("float32"))

Why seams appear and where feathering fits

Satellite mosaics are built from scenes acquired at different times and angles. Even after orthorectification and atmospheric correction, adjacent tiles rarely share identical surface reflectance values at their shared boundary: atmospheric haze gradients, BRDF effects from differing sun-zenith angles, and sensor gain drift all leave radiometric offsets that manifest as visible seams.

Feathering is a terminal compositing step, not a radiometric correction. It assumes that upstream steps — cloud and shadow masking, atmospheric normalization, and precise geometric registration — have already been applied. What it solves is the blending geometry: given two overlapping, already-normalized tiles, how should pixel values be combined in the overlap zone to produce a visually continuous surface? The answer is to weight each pixel by its distance from the nearest invalid boundary, so only interior, well-observed pixels dominate the output.

The diagram below shows how weights shift from Scene A to Scene B across the overlap region:

Feathering weight transition across a mosaic overlap zone Two overlapping raster tiles. Scene A weight starts high on the left and falls to zero at the right edge. Scene B weight starts at zero on the left and rises to high on the right. The combined weight sum is constant across the overlap zone. Scene A Scene B overlap α_A = 1 α_B = 1 ← pixel position across boundary → α_A + α_B = const Weight derived from distance_transform_edt on valid-pixel mask

Environment and setup

Package Minimum version Role
rasterio 1.3 Reading/writing georeferenced rasters
numpy 1.24 Array accumulation and weighted average
scipy 1.10 distance_transform_edt for weight generation 
pip install rasterio>=1.3 numpy>=1.24 scipy>=1.10

All inputs must share the same CRS, spatial resolution, and pixel grid origin. Verify alignment with rasterio.transform before running any compositing operation — misaligned grids will cause ghosting that feathering cannot fix.

Complete working example

The function below handles any number of spatially aligned tiles, arbitrary band counts, and mixed nodata values. It initializes accumulator arrays on the first tile, then streams through the remaining tiles without loading more than one tile at a time.

import rasterio
import numpy as np
from scipy.ndimage import distance_transform_edt
from pathlib import Path


def feather_mosaic(
    input_paths: list[str],
    output_path: str,
    nodata: float = -9999.0,
) -> None:
    """
    Create a distance-weighted feathered mosaic from spatially aligned rasters.

    All inputs must share identical CRS, resolution, and grid alignment.
    Inputs may have non-overlapping or partially overlapping valid-data extents.

    Args:
        input_paths: Ordered list of raster file paths (GeoTIFF or any GDAL-readable format).
        output_path: Destination path for the blended output raster.
        nodata: Sentinel value indicating invalid or missing pixels in the inputs.
    """
    if len(input_paths) < 2:
        raise ValueError("At least two overlapping rasters are required.")

    # Accumulator arrays; initialised lazily on the first tile
    mosaic_num: np.ndarray | None = None   # weighted sum of pixel values
    mosaic_den: np.ndarray | None = None   # sum of weights
    profile: dict | None = None

    for path in input_paths:
        with rasterio.open(path) as src:
            # Initialise once, using the first tile's spatial metadata
            if profile is None:
                profile = src.profile.copy()
                profile.update(dtype="float32", nodata=np.nan, count=src.count)
                mosaic_num = np.zeros(
                    (src.count, src.height, src.width), dtype=np.float32
                )
                mosaic_den = np.zeros((src.height, src.width), dtype=np.float32)

            data = src.read().astype(np.float32)

            # 2D validity mask: True where every band has a valid pixel.
            # Using np.all(..., axis=0) prevents spectral leakage — a pixel
            # masked in *any* band is excluded from blending entirely.
            valid_mask = np.all(data != nodata, axis=0)

            # Skip tiles that are entirely nodata (e.g. cloud-only scenes)
            if not np.any(valid_mask):
                continue

            # distance_transform_edt assigns each True (valid) pixel its
            # Euclidean distance to the nearest False (invalid/edge) pixel.
            # Interior pixels receive large values; edge pixels receive ~0.
            dist = distance_transform_edt(valid_mask)

            # Normalise to [0, 1]; add epsilon to guard against single-pixel tiles
            max_dist = float(dist.max())
            alpha = dist / (max_dist + 1e-6)

            # Broadcast alpha across all bands for element-wise multiplication
            alpha_3d = alpha[np.newaxis, :, :]  # shape: (1, H, W) → broadcast to (C, H, W)

            # Accumulate weighted values and the weights themselves
            mosaic_num += data * alpha_3d
            mosaic_den += alpha

    assert profile is not None and mosaic_num is not None and mosaic_den is not None

    # Divide accumulated weighted values by total weight; NaN where no tile covered
    with np.errstate(divide="ignore", invalid="ignore"):
        mosaic_final = np.where(
            mosaic_den[np.newaxis, :, :] > 0,
            mosaic_num / mosaic_den[np.newaxis, :, :],
            np.nan,
        )

    with rasterio.open(output_path, "w", **profile) as dst:
        dst.write(mosaic_final.astype("float32"))

Key implementation notes

  • Distance transform direction. distance_transform_edt(valid_mask) measures distance from each valid pixel to the nearest invalid pixel (not the other way). Interior valid pixels get large distances — high weights. Edge pixels approach zero. This is the correct direction: interior observations dominate the blend.
  • Spectral consistency. np.all(data != nodata, axis=0) masks a pixel if any band is invalid. If you used np.any(...) instead, you would blend a fully valid pixel in Band 1 against a nodata pixel in Band 2, producing physically meaningless spectral vectors.
  • Numerical stability. The 1e-6 epsilon in the alpha normalisation prevents zero-division for pathological single-pixel valid regions. The np.errstate guard suppresses runtime warnings from the final denominator division.
  • Profile propagation. Copying src.profile from the first tile preserves CRS, affine transform, and driver metadata. profile.update(dtype="float32", nodata=np.nan) ensures the output uses NaN as the nodata sentinel, consistent with float32 conventions.

Variant patterns

Variant 1: Normalised weight (prevent oversaturation on high-count stacks)

When compositing many tiles — for example, a full Sentinel-2 annual stack — the raw dist.max() normalisation can produce very similar alpha values across tiles, reducing the discriminating power of the distance weighting. Normalise each tile’s alpha to sum to 1 before accumulation to keep total weight contributions balanced:

# Replace the normalisation block with:
alpha_sum = float(dist.sum())
alpha = dist / (alpha_sum + 1e-6)  # sum-normalised instead of max-normalised

This variant is useful when overlap zones span hundreds of pixels and you want the interior of each tile to contribute proportionally rather than uniformly.

Variant 2: Dask-enabled out-of-core blending

Full-orbit Sentinel-2 or Landsat scenes can exceed available RAM. Replace numpy with dask.array to compute distance transforms and accumulation lazily, executing only when the output array is written:

import dask.array as da
from scipy.ndimage import distance_transform_edt
import rasterio
import numpy as np


def feather_mosaic_dask(
    input_paths: list[str],
    output_path: str,
    nodata: float = -9999.0,
    chunk_size: int = 2048,
) -> None:
    """Dask-backed feathering for scenes larger than available RAM."""
    mosaic_num: da.Array | None = None
    mosaic_den: da.Array | None = None
    profile: dict | None = None

    for path in input_paths:
        with rasterio.open(path) as src:
            if profile is None:
                profile = src.profile.copy()
                profile.update(dtype="float32", nodata=np.nan)

            # Read into Dask array via numpy; for COG sources use rioxarray instead
            data_np = src.read().astype(np.float32)
            data = da.from_array(data_np, chunks=(src.count, chunk_size, chunk_size))

            valid_np = np.all(data_np != nodata, axis=0)
            if not valid_np.any():
                continue

            # distance_transform_edt is not natively lazy; compute on the mask
            # (mask is small compared to the full float array)
            dist_np = distance_transform_edt(valid_np).astype(np.float32)
            alpha = da.from_array(
                dist_np / (float(dist_np.max()) + 1e-6),
                chunks=(chunk_size, chunk_size),
            )

            if mosaic_num is None:
                mosaic_num = data * alpha[np.newaxis]
                mosaic_den = alpha
            else:
                mosaic_num = mosaic_num + data * alpha[np.newaxis]
                mosaic_den = mosaic_den + alpha

    assert profile is not None and mosaic_num is not None and mosaic_den is not None

    result = da.where(mosaic_den[np.newaxis] > 0, mosaic_num / mosaic_den[np.newaxis], np.nan)
    result_np = result.compute().astype("float32")

    with rasterio.open(output_path, "w", **profile) as dst:
        dst.write(result_np)

Variant 3: Windowed I/O for tile-by-tile flushing

For pipelines integrated with automated image clipping and cropping where individual window reads are already established, accumulate feathering per block and flush numerator/denominator arrays to disk between windows. This avoids holding the full mosaic in memory and is compatible with rasterio’s Window API for processing tiles of 1024×1024 pixels.

from rasterio.windows import Window
import rasterio
import numpy as np
from scipy.ndimage import distance_transform_edt


def feather_mosaic_windowed(
    input_paths: list[str],
    output_path: str,
    nodata: float = -9999.0,
    block_size: int = 1024,
) -> None:
    """Window-by-window feathering: constant memory footprint regardless of scene size."""
    # Read profile from first source
    with rasterio.open(input_paths[0]) as src0:
        profile = src0.profile.copy()
        profile.update(dtype="float32", nodata=np.nan)
        height, width, n_bands = src0.height, src0.width, src0.count

    with rasterio.open(output_path, "w", **profile) as dst:
        for row_off in range(0, height, block_size):
            for col_off in range(0, width, block_size):
                win = Window(
                    col_off, row_off,
                    min(block_size, width - col_off),
                    min(block_size, height - row_off),
                )
                h, w = win.height, win.width
                num = np.zeros((n_bands, h, w), dtype=np.float32)
                den = np.zeros((h, w), dtype=np.float32)

                for path in input_paths:
                    with rasterio.open(path) as src:
                        data = src.read(window=win).astype(np.float32)
                    valid = np.all(data != nodata, axis=0)
                    if not np.any(valid):
                        continue
                    dist = distance_transform_edt(valid)
                    alpha = dist / (dist.max() + 1e-6)
                    num += data * alpha[np.newaxis]
                    den += alpha

                with np.errstate(divide="ignore", invalid="ignore"):
                    block = np.where(den[np.newaxis] > 0, num / den[np.newaxis], np.nan)
                dst.write(block.astype("float32"), window=win)

Note: windowed distance transforms operate only on the block interior, so edge weights at block boundaries will underestimate true interior distances. For critical blends, increase block_size or add a halo region that is cropped on write.

Common errors

ValueError: operands could not be broadcast together The input tiles have different band counts or array shapes. Confirm that all tiles share the same height, width, and band count before passing them to the mosaic function. Use rasterio.open(path).count and .shape to verify.

Ghosting stripes at tile edges in the output Inputs are not on the same pixel grid. Even a sub-pixel offset causes each tile’s pixels to fall at different spatial positions, and feathering cannot compensate. Reproject all tiles to a common CRS and snapped grid using rasterio.warp.reproject or gdalwarp -tap before compositing. The Advanced Resampling and Upscaling Techniques cluster covers choosing the right resampling kernel for this step.

Feathered output is uniformly NaN or all-nodata The nodata argument does not match the actual nodata value encoded in the files. Read the stored nodata with src.nodata and pass it explicitly: feather_mosaic(paths, out, nodata=src.nodata). If inputs use different nodata values, normalise them to a single sentinel before blending.