Mosaicking Tiled GHI Rasters with rasterio.merge #

Satellite and reanalysis GHI archives ship as adjacent tiles — one GeoTIFF per NSRDB grid cell or per PVGIS download box — and turning them into a single seamless surface is the first thing that breaks before any of the alignment work in Solar Irradiance Raster Processing can proceed. The naive call merge([rasterio.open(p) for p in tile_paths]) returns something that looks like a mosaic but carries visible seams along tile boundaries, black rectangular borders where nodata was never declared, or it dies outright with MemoryError on a continental extent. Worse, the borders that bleed into the array are not obviously wrong: a 0 from an undeclared nodata region reads as “zero irradiance,” and a mean taken over that mosaic silently understates the resource a lender treats as ground truth. This page names the four failure signatures, shows a preflight that catches them before the merge runs, and gives a corrected merge plus a windowed VRT fallback for archives too large to hold in RAM.

The whole point of a mosaic is to compute honest area statistics over it. The area-weighted mean of a merged GHI surface is only meaningful once nodata is excluded from both the numerator and the weight:

where $A_i$ is the ground area of pixel $i$. If black-border pixels leak in with value 0 and area weight, the denominator inflates while the numerator does not, and $\overline{\mathrm{GHI}}$ drifts low in direct proportion to how much of the mosaic footprint is padding.

Root-cause analysis #

Four compounding causes account for nearly every broken GHI mosaic, and each maps to a specific fix below:

1. Tiles in different CRS. rasterio.merge does not reproject. It assumes every input shares one CRS and pastes pixels by their affine transforms. Feed it one tile in EPSG:4326 and its neighbour in EPSG:32610 and you get seams, gaps, or a mosaic where half the tiles land in the wrong place — no exception is raised because each file is individually valid. Uniform coordinate reference system alignment across the tile set is a precondition, not an afterthought.
2. nodata not set, so black borders bleed. When a tile has no nodata declared, merge treats its fill value (often 0) as valid data. Overlapping fill from one tile overwrites real irradiance in its neighbour, and the padded edges enter every downstream mean as spurious zeros — the black-border artefact.
3. dtype or resolution mismatch. merge requires a single output dtype and a single pixel size. A uint16 scaled-integer tile beside a float32 tile, or a 0.5° tile beside a 0.01° tile, triggers a dtype error or forces a silent, unintended resample that shifts pixel registration by a fraction of a cell.
4. Overlapping tiles, wrong merge method. Adjacent downloads usually overlap by a few pixels. The default method="first" keeps whichever tile was read first in the overlap zone; if that tile’s edge is cloud-contaminated or nodata-padded, its garbage wins. Choosing the method deliberately ("last", "min", "max", or a custom callable) is what controls the seam.

Beyond these, a whole-mosaic merge that materialises every tile in memory at once is the classic MemoryError on continental extents — addressed by the VRT fallback rather than by buying RAM.

Pre-flight validation #

Surface all four causes before merge runs. The validator opens each tile’s header only — never its full array — and raises on the exact defect so a CI/CD run fails fast with a precise message instead of writing a mosaic riddled with seams and bleed:

import rasterio
from rasterio.crs import CRS


def preflight_tile_set(tile_paths: list[str], target_epsg: int = 32610) -> None:
    """Raise on CRS, nodata, dtype, or resolution divergence before merge()."""
    target_crs = CRS.from_epsg(target_epsg)
    crs_seen, dtype_seen, res_seen, missing_nodata = set(), set(), set(), []

    for path in tile_paths:
        with rasterio.open(path) as tile:
            crs_seen.add(tile.crs.to_epsg() if tile.crs else None)
            dtype_seen.add(tile.dtypes[0])
            res_seen.add((round(tile.res[0], 6), round(tile.res[1], 6)))
            if tile.nodata is None:
                missing_nodata.append(path)

    # Cause 1: every tile must share the one metric CRS merge() will assume
    if crs_seen != {target_epsg}:
        raise ValueError(
            f"CRS divergence across tiles: {crs_seen}. "
            f"Reproject all tiles to EPSG:{target_epsg} before mosaicking."
        )
    # Cause 3: merge() cannot reconcile mixed dtypes or pixel sizes
    if len(dtype_seen) > 1:
        raise ValueError(f"dtype mismatch across tiles: {dtype_seen}. Cast to one dtype first.")
    if len(res_seen) > 1:
        raise ValueError(f"resolution mismatch across tiles: {res_seen}. Resample to one grid first.")
    # Cause 2: undeclared nodata is what bleeds black borders into the mean
    if missing_nodata:
        raise ValueError(
            f"{len(missing_nodata)} tile(s) have no nodata declared; "
            "fill pixels will bleed into the mosaic. Set nodata before merge()."
        )

Fix implementation #

The corrected function runs the preflight, then calls rasterio.merge with an explicit nodata and a deliberately chosen method so overlap zones resolve predictably. Parameter choices are justified for GHI use: nodata=np.nan with a float32 output keeps padded edges out of every downstream mean; method="max" favours the cloud-free reading in overlaps (clouds depress GHI, so the maximum of two co-located samples is the clearer-sky value); and resampling=Resampling.bilinear preserves radiometric continuity if merge must nudge a tile onto the common grid.

import numpy as np
import rasterio
from rasterio.merge import merge


def mosaic_ghi_tiles(
    tile_paths: list[str],
    dst_path: str,
    target_epsg: int = 32610,
    merge_method: str = "max",   # clouds depress GHI; max favours clear-sky overlap
) -> dict:
    """Mosaic aligned GHI tiles with explicit nodata and overlap handling."""
    preflight_tile_set(tile_paths, target_epsg)  # fail fast on the four root causes

    srcs = [rasterio.open(p) for p in tile_paths]
    try:
        ghi_array, out_transform = merge(
            srcs,
            nodata=np.nan,          # padded edges excluded from data, not read as 0
            method=merge_method,    # deterministic winner in overlap zones
            resampling=rasterio.enums.Resampling.bilinear,
            dtype="float32",
        )
        profile = srcs[0].profile | {
            "driver": "GTiff",
            "height": ghi_array.shape[1],
            "width": ghi_array.shape[2],
            "count": ghi_array.shape[0],
            "dtype": "float32",
            "nodata": np.nan,
            "crs": rasterio.crs.CRS.from_epsg(target_epsg),
            "transform": out_transform,
            "tiled": True, "blockxsize": 512, "blockysize": 512,
            "compress": "lzw",      # GHI surfaces compress well; keeps archives small
        }
        with rasterio.open(dst_path, "w", **profile) as dst:
            dst.write(ghi_array)
            dst.update_tags(SOURCE="NSRDB tiles", MERGE_METHOD=merge_method, CRS_EPSG=str(target_epsg))
    finally:
        for src in srcs:
            src.close()

    valid = np.isfinite(ghi_array)
    return {
        "shape": tuple(ghi_array.shape),
        "valid_frac": float(valid.mean()),
        "mean_ghi": float(np.nanmean(ghi_array)),
    }

Passing nodata=np.nan is the single detail that eliminates the black-border artefact: merge writes NaN into every gap and every masked overlap, so np.nanmean and the area-weighted mean above see only real irradiance. Choosing method="max" over the default "first" removes the “whichever tile loaded first wins” nondeterminism that makes mosaics irreproducible between runs.

Fallback routing & performance tuning #

When the mosaic will not fit in RAM — a continental NSRDB stack at native resolution routinely exceeds a workstation’s memory — layer these strategies on top of the core function:

• Build a VRT instead of a materialised mosaic. gdal.BuildVRT("stack.vrt", tile_paths) (or gdalbuildvrt on the CLI) creates a virtual mosaic that references the tiles on disk with zero pixel copies. Read and write it in windows so peak memory stays proportional to one block, not the whole extent.
• Window the write. Iterate dst.block_windows() and merge(srcs, bounds=window_bounds, ...) per block, writing each tile of output as it is computed. This is the memory-safe equivalent of the in-RAM path and is what continental runs should default to.
• Cast to float32, compress with LZW. float32 halves the footprint of a float64 intermediate at no cost to GHI precision, and LZW compresses smooth irradiance fields well — both shrink the archive and the I/O the merge must stream.
• Cap GDAL_CACHEMAX. Set it (e.g. 512 MB) so GDAL’s block cache, not the data, stops driving peak memory during a batch mosaic.
• Match block size to the source tiling. Aligning blockxsize/blockysize with the tiles’ internal geometry (often 256 or 512) avoids re-blocking overhead on every window read. For the kernel trade-offs when a resample is unavoidable, see the resampling and raster kernel quick reference.

Downstream validation #

Before a mosaic feeds a yield model or a spatial join, gate it with an assertion suitable for a CI/CD pipeline. This catches nodata bleed, an unexpected footprint, and any band-count or dtype regression introduced upstream:

import numpy as np
import rasterio


def assert_mosaic_integrity(
    dst_path: str,
    expected_epsg: int = 32610,
    expected_bounds: tuple | None = None,
    max_zero_frac: float = 0.001,
) -> None:
    """CI/CD gate: fail the build if the GHI mosaic is not assessment-grade."""
    with rasterio.open(dst_path) as mosaic:
        assert mosaic.crs.to_epsg() == expected_epsg, "mosaic lost its target CRS"
        assert mosaic.count == 1, f"expected 1 GHI band, got {mosaic.count}"
        assert mosaic.dtypes[0] == "float32", f"expected float32, got {mosaic.dtypes[0]}"
        assert np.isnan(mosaic.nodata), "nodata must be NaN so borders stay excluded"

        ghi_array = mosaic.read(1)
        # nodata bleed check: near-zero pixels signal undeclared fill leaking in as 0
        finite = ghi_array[np.isfinite(ghi_array)]
        zero_frac = float((finite == 0).mean()) if finite.size else 1.0
        assert zero_frac <= max_zero_frac, (
            f"{zero_frac:.2%} of pixels are exactly 0 — black-border bleed suspected"
        )
        assert np.nanmax(ghi_array) <= 13.0, "GHI exceeds physical ceiling (kWh/m²/day)"

        if expected_bounds is not None:
            assert np.allclose(mosaic.bounds, expected_bounds, atol=mosaic.res[0]), (
                f"mosaic footprint {tuple(mosaic.bounds)} != expected {expected_bounds}"
            )

The zero-fraction test is the specific guard against Cause 2: a correctly masked mosaic has essentially no exactly-zero pixels, so a spike in that fraction is the fingerprint of undeclared nodata bleeding in as 0. Logging valid_frac and the merge method as provenance keeps the mosaic auditable — an independent reviewer can see how much of the footprint was real data versus padding, mirroring the alignment discipline the aligned surfaces carry into Stacking NASA POWER and PVGIS rasters in rasterio. Pin rasterio and GDAL versions in pyproject.toml so a default-method change in merge cannot silently shift the mosaic between runs.

Related #

• Solar Irradiance Raster Processing — the parent workflow this mosaic feeds, with the full reproject-and-validate stage.
• Stacking NASA POWER and PVGIS rasters in rasterio — resolving overlap and registration when combining multi-source surfaces.
• Resampling and Raster Kernel Quick Reference — choosing bilinear vs average vs nearest when a merge must resample.
• Coordinate Reference Systems for Energy Projects — the uniform-CRS enforcement mosaicking depends on.