Merge & stack #

Free functions for combining multiple rasters into one.

Hold "Ctrl" to enable pan & zoom

flowchart LR
    M["<b>pyramids.dataset.merge</b>"]
    M --> MR["<b>merge_rasters</b><br/>mosaic overlapping rasters<br/>into one (reproject on the fly)"]
    M --> SB["<b>stack_bands</b><br/>stack single-band rasters<br/>into one multi-band Dataset"]

merge_rasters — mosaic several (overlapping or adjacent) rasters into a single raster covering their union.
stack_bands — stack several single-band rasters into one multi-band raster.

See the Mosaic & merge notebook for runnable examples.

`pyramids.dataset.merge.merge_rasters(src, dst, no_data_value='0', init='nan', n='nan', method='last', dst_crs=None, resampling='nearest neighbor', signer=None)` #

Merge a group of rasters into one raster, resolving overlaps by method.

The overlap-resolution method selects how overlapping pixels are combined:

"last" (default) / "first" — z-order compositing: the last (or first) source covering a pixel wins. Implemented cheaply with :func:gdal.BuildVRT + :func:gdal.Translate.
"min" / "max" / "sum" — per-pixel reduction across every source overlapping that pixel, ignoring no-data. Each source is aligned onto the union grid and the bands are stacked and reduced with NaN-aware numpy.

Parameters:

Name	Type	Description	Default
`src`	`Sequence[str \| Path]`	Paths to all input rasters.	required
`dst`	`str \| Path`	Path to the output raster.	required
`no_data_value`	`float \| int \| str`	Stamped on the output bands as the nodata marker. For the reduction methods it also fills pixels with no source coverage.	`'0'`
`init`	`float \| int \| str`	Reported value for pixels with no source coverage in the VRT (z-order methods only). Maps to :func:`gdal.BuildVRTOptions` `VRTNodata`.	`'nan'`
`n`	`float \| int \| str`	Source pixels matching this value are ignored — both when building the VRT mosaic (z-order) and when reducing (the value is treated as source no-data).	`'nan'`
`method`	`str`	Overlap-resolution rule: one of `"first"`, `"last"` (default), `"min"`, `"max"`, `"sum"`.	`'last'`
`dst_crs`	`int \| str \| None`	Target CRS for the mosaic, as an EPSG code (`32632`) or any GDAL-parseable CRS string (`"EPSG:32632"`, a WKT, a PROJ string). Each source whose CRS differs from the target is reprojected onto it (via :func:`gdal.Warp`) before compositing, so tiles in different CRSs — e.g. a Sentinel-2 AOI straddling two UTM zones — mosaic correctly. `None` (default) keeps the previous behaviour: sources are assumed to share a CRS and are composited as-is, except when they are found to disagree, in which case they are reprojected onto the first source's CRS. Reprojection always happens before the `BuildVRT`/`Warp` compositing step because that step has no reprojection capability and assumes a single shared grid.	`None`
`resampling`	`str`	Resampling method used when a source is reprojected to `dst_crs` (or to the common CRS on auto-detect). Case-insensitive; any key of :data:`pyramids.base._utils.INTERPOLATION_METHODS`: `"nearest"` (alias `"nearest neighbor"`, the default), `"bilinear"`, `"cubic"`, `"cubic_spline"`, `"lanczos"`, `"average"`, `"mode"`, `"max"`, `"min"`, `"med"`, `"q1"`, `"q3"`, `"sum"`, and `"rms"`. Prefer `"bilinear"`/`"cubic"` for continuous data (reflectance, DEM) to avoid the blockiness nearest introduces across reprojection. Ignored when no source is reprojected.	`'nearest neighbor'`
`signer`	`Any`	Optional signer exposing `sign_href(str) -> str` and `gdal_env() -> dict[str, str]` (e.g. a :class:`pyramids.stac.signers.Signer`). When given, both hooks are applied — exactly as :func:`pyramids.stac.load_asset` does: `signer.sign_href` rewrites every source path first (e.g. grafting a SAS token onto a blob URL), then `signer.gdal_env()` is installed via :class:`~pyramids.base.remote.CloudConfig` for the duration of the merge. This means URL-signing signers (Planetary Computer SAS, whose credential rides the href and whose `gdal_env()` is empty) and env-based signers (Requester-Pays, bearer) both authenticate without wrapping the call in a `with CloudConfig(...)` block. `None` (default) leaves source hrefs untouched and installs no extra config.	`None`

Returns:

Type	Description
`None`	None

Note

The z-order methods ("first"/"last") preserve each source's data type via BuildVRT + Translate. The reduction methods ("min"/"max"/"sum") align every source onto the union grid with gdal.Warp (nearest resampling — exact for already-aligned tiles) and write a single-precision-safe Float64 output regardless of the source dtype, so they may differ in dtype from a z-order merge of the same integer inputs.

Raises:

Type	Description
`TypeError`	`resampling` is not a string.
`ValueError`	`method`/`resampling` is not a supported value, `dst_crs` cannot be parsed as a CRS, or a source carries no CRS.
`RuntimeError`	GDAL failed to open a source, reproject it, or build the source mosaic.

Examples:

Mosaic two tiles, keeping the larger value wherever they overlap:

>>> from pyramids.dataset.merge import merge_rasters
>>> merge_rasters(  # doctest: +SKIP
...     ["tile_a.tif", "tile_b.tif"],
...     "mosaic_max.tif",
...     no_data_value=-9999.0,
...     method="max",
... )

Default last-wins compositing (unchanged from the previous behaviour):

>>> merge_rasters(["tile_a.tif", "tile_b.tif"], "mosaic.tif")  # doctest: +SKIP

Mosaic tiles from two UTM zones into a single CRS:

>>> merge_rasters(  # doctest: +SKIP
...     ["utm32_tile.tif", "utm33_tile.tif"],
...     "mosaic_utm32.tif",
...     dst_crs=32632,
... )

Mosaic Requester-Pays S3 tiles by passing a signer (no with block):

>>> from pyramids.stac import AWSRequesterPaysSigner  # doctest: +SKIP
>>> merge_rasters(  # doctest: +SKIP
...     ["s3://bucket/a.tif", "s3://bucket/b.tif"],
...     "mosaic.tif",
...     signer=AWSRequesterPaysSigner(region="us-west-2"),
... )

Source code in src/pyramids/dataset/merge.py

def merge_rasters(
    src: Sequence[str | Path],
    dst: str | Path,
    no_data_value: float | int | str = "0",
    init: float | int | str = "nan",
    n: float | int | str = "nan",
    method: str = "last",
    dst_crs: int | str | None = None,
    resampling: str = "nearest neighbor",
    signer: Any = None,
) -> None:
    """Merge a group of rasters into one raster, resolving overlaps by ``method``.

    The overlap-resolution ``method`` selects how overlapping pixels are
    combined:

    * ``"last"`` (default) / ``"first"`` — z-order compositing: the last (or
      first) source covering a pixel wins. Implemented cheaply with
      :func:`gdal.BuildVRT` + :func:`gdal.Translate`.
    * ``"min"`` / ``"max"`` / ``"sum"`` — per-pixel reduction across every
      source overlapping that pixel, ignoring no-data. Each source is aligned
      onto the union grid and the bands are stacked and reduced with NaN-aware
      numpy.

    Args:
        src (Sequence[str | Path]):
            Paths to all input rasters.
        dst (str | Path):
            Path to the output raster.
        no_data_value (float | int | str):
            Stamped on the output bands as the nodata marker. For the reduction
            methods it also fills pixels with no source coverage.
        init (float | int | str):
            Reported value for pixels with no source coverage in the VRT (z-order
            methods only). Maps to :func:`gdal.BuildVRTOptions` ``VRTNodata``.
        n (float | int | str):
            Source pixels matching this value are ignored — both when building
            the VRT mosaic (z-order) and when reducing (the value is treated as
            source no-data).
        method (str):
            Overlap-resolution rule: one of ``"first"``, ``"last"`` (default),
            ``"min"``, ``"max"``, ``"sum"``.
        dst_crs (int | str | None):
            Target CRS for the mosaic, as an EPSG code (``32632``) or any
            GDAL-parseable CRS string (``"EPSG:32632"``, a WKT, a PROJ string).
            Each source whose CRS differs from the target is reprojected onto it
            (via :func:`gdal.Warp`) **before** compositing, so tiles in different
            CRSs — e.g. a Sentinel-2 AOI straddling two UTM zones — mosaic
            correctly. ``None`` (default) keeps the previous behaviour: sources
            are assumed to share a CRS and are composited as-is, *except* when
            they are found to disagree, in which case they are reprojected onto
            the first source's CRS. Reprojection always happens before the
            ``BuildVRT``/``Warp`` compositing step because that step has no
            reprojection capability and assumes a single shared grid.
        resampling (str):
            Resampling method used when a source is reprojected to ``dst_crs``
            (or to the common CRS on auto-detect). Case-insensitive; any key
            of :data:`pyramids.base._utils.INTERPOLATION_METHODS`: ``"nearest"``
            (alias ``"nearest neighbor"``, the default), ``"bilinear"``,
            ``"cubic"``, ``"cubic_spline"``, ``"lanczos"``, ``"average"``,
            ``"mode"``, ``"max"``, ``"min"``, ``"med"``, ``"q1"``, ``"q3"``,
            ``"sum"``, and ``"rms"``.
            Prefer ``"bilinear"``/``"cubic"`` for continuous data (reflectance,
            DEM) to avoid the blockiness nearest introduces across reprojection.
            Ignored when no source is reprojected.
        signer (Any):
            Optional signer exposing ``sign_href(str) -> str`` and
            ``gdal_env() -> dict[str, str]`` (e.g. a
            :class:`pyramids.stac.signers.Signer`). When given, **both** hooks
            are applied — exactly as :func:`pyramids.stac.load_asset` does:
            ``signer.sign_href`` rewrites every source path first (e.g. grafting
            a SAS token onto a blob URL), then ``signer.gdal_env()`` is installed
            via :class:`~pyramids.base.remote.CloudConfig` for the duration of
            the merge. This means URL-signing signers (Planetary Computer SAS,
            whose credential rides the href and whose ``gdal_env()`` is empty)
            and env-based signers (Requester-Pays, bearer) both authenticate
            without wrapping the call in a ``with CloudConfig(...)`` block.
            ``None`` (default) leaves source hrefs untouched and installs no
            extra config.

    Returns:
        None

    Note:
        The z-order methods (``"first"``/``"last"``) preserve each source's data
        type via ``BuildVRT`` + ``Translate``. The reduction methods
        (``"min"``/``"max"``/``"sum"``) align every source onto the union grid
        with ``gdal.Warp`` (nearest resampling — exact for already-aligned tiles)
        and write a single-precision-safe **Float64** output regardless of the
        source dtype, so they may differ in dtype from a z-order merge of the
        same integer inputs.

    Raises:
        TypeError: ``resampling`` is not a string.
        ValueError: ``method``/``resampling`` is not a supported value,
            ``dst_crs`` cannot be parsed as a CRS, or a source carries no CRS.
        RuntimeError: GDAL failed to open a source, reproject it, or build the
            source mosaic.

    Examples:
        - Mosaic two tiles, keeping the larger value wherever they overlap:
            ```python
            >>> from pyramids.dataset.merge import merge_rasters
            >>> merge_rasters(  # doctest: +SKIP
            ...     ["tile_a.tif", "tile_b.tif"],
            ...     "mosaic_max.tif",
            ...     no_data_value=-9999.0,
            ...     method="max",
            ... )

            ```
        - Default last-wins compositing (unchanged from the previous behaviour):
            ```python
            >>> merge_rasters(["tile_a.tif", "tile_b.tif"], "mosaic.tif")  # doctest: +SKIP

            ```
        - Mosaic tiles from two UTM zones into a single CRS:
            ```python
            >>> merge_rasters(  # doctest: +SKIP
            ...     ["utm32_tile.tif", "utm33_tile.tif"],
            ...     "mosaic_utm32.tif",
            ...     dst_crs=32632,
            ... )

            ```
        - Mosaic Requester-Pays S3 tiles by passing a signer (no ``with`` block):
            ```python
            >>> from pyramids.stac import AWSRequesterPaysSigner  # doctest: +SKIP
            >>> merge_rasters(  # doctest: +SKIP
            ...     ["s3://bucket/a.tif", "s3://bucket/b.tif"],
            ...     "mosaic.tif",
            ...     signer=AWSRequesterPaysSigner(region="us-west-2"),
            ... )

            ```
    """
    if method not in _MERGE_METHODS:
        raise ValueError(
            f"method must be one of {list(_MERGE_METHODS)}, got {method!r}."
        )

    # SMELL: `init` and `n` default to the string `"nan"`, which
    # round-trips through GDAL as float NaN. For integer-typed
    # rasters (e.g. UInt16) GDAL emits a warning per band:
    # `Band data type of <T> cannot represent the specified NoData
    # value of nan`. The defaults are kept for backwards-compat
    # with the previous gdal_merge.main-based signature; callers
    # that hit integer rasters should pass an explicit numeric
    # value instead of relying on the default.
    src_paths = [str(p) for p in src]
    if signer is not None:
        # Apply the signer's href rewrite to every source (e.g. graft a SAS
        # token onto a blob URL) so URL-signing signers authenticate. A no-op
        # for signers that authenticate via gdal_env() only — the base
        # sign_href returns the href unchanged. This mirrors load_asset, which
        # applies BOTH signer hooks (sign_href + gdal_env); applying only the
        # env half here would silently read URL-signed sources unauthenticated.
        src_paths = [signer.sign_href(p) for p in src_paths]

    # All GDAL reads/writes run under the signer's cloud config (a no-op when
    # signer is None) so authenticated remote sources open with the right
    # credentials for the whole merge.
    with _cloud_config(signer):
        # Put every source on one CRS before compositing. The BuildVRT/Warp
        # mosaic below cannot reproject — it stitches pixel grids assuming a
        # shared CRS — so mismatched sources must be warped first or they would
        # mis-align silently. `_keepalive` holds the in-memory warped VRTs so
        # GDAL does not free them while the mosaic is built.
        sources, _keepalive = _prepare_sources(src_paths, dst_crs, resampling)

        if method in _REDUCE_METHODS:
            _merge_reduce(sources, str(dst), method, no_data_value, n)
            return

        # z-order: "last" keeps natural order (last source wins); "first"
        # reverses so the original first source is placed last in the VRT and
        # therefore wins.
        ordered = list(reversed(sources)) if method == "first" else sources
        vrt_opts = gdal.BuildVRTOptions(
            srcNodata=str(n),
            VRTNodata=str(init),
        )
        vrt_ds = gdal.BuildVRT("", ordered, options=vrt_opts)
        if vrt_ds is None:
            raise RuntimeError(
                f"gdal.BuildVRT returned None for sources {src_paths!r}; "
                "check that all paths are readable rasters with consistent "
                "band counts and CRS."
            )
        translate_opts = gdal.TranslateOptions(
            creationOptions=["COMPRESS=LZW"],
            noData=str(no_data_value),
        )
        out_ds = gdal.Translate(str(dst), vrt_ds, options=translate_opts)
        out_ds.FlushCache()
        out_ds = None
        vrt_ds = None

`pyramids.dataset.merge.stack_bands(files, *, band_names=None, align=False, no_data_value=_INHERIT_NO_DATA, path=None, signer=None)` #

Stack N single-band rasters into one multi-band :class:Dataset.

Free-function alias for :meth:pyramids.dataset.Dataset.from_band_files — see that method for the full contract, edge cases, and examples.