Merging Tiled Lambda Outputs into a Cloud Optimized GeoTIFF

Merge the per-tile outputs by building a VRT over them with gdal.BuildVRT — referencing each tile in place through a /vsis3 path rather than downloading it — then translate that VRT to a Cloud Optimized GeoTIFF with gdal.Translate (COG driver) or rio-cogeo, writing internal tiling and overviews in one pass. The VRT is a kilobyte-scale XML index, so peak memory stays bounded to a few blocks and only the final product touches /tmp; guard the 10 GB ephemeral ceiling with an os.statvfs check before the translate and the whole fan-in runs comfortably inside a single Lambda invocation.

Merging tiled outputs into a COG The merge Lambda lists per-tile GeoTIFFs in S3, builds a VRT that references them in place via /vsis3, translates the VRT to a Cloud Optimized GeoTIFF with overviews while guarding /tmp, and uploads the product back to S3. tile 000_000.tif tile 000_512.tif tile N.tif S3 /vsis3 BuildVRT XML index, KB Translate → COG overviews + tiling block-by-block read product.tif valid COG → S3 os.statvfs guard: /tmp ≤ 10 GB Tiles referenced in place; only the final COG lands on /tmp — the VRT keeps memory and disk bounded.

Context

This is the fan-in stage of the 10 GB GeoTIFF recipe: after the distributed Map state has produced one GeoTIFF per window — each written to a deterministic key by the partitioning step described in partitioning a GeoTIFF into Map-state tiles — something has to stitch them back into one deliverable. The naive approach, gdal_merge.py, copies every pixel of the mosaic into a new raster before writing, which for a reassembled 10 GB scene overruns both the memory ceiling and the ephemeral /tmp limit.

A VRT sidesteps that entirely. gdal.BuildVRT writes a small XML file that references the tiles in place through /vsis3 paths; no pixels are copied at build time. gdal.Translate then reads through the VRT block by block to produce the final COG, so peak memory is a handful of blocks rather than the whole mosaic. The tradeoff is that the merge Lambda must read every tile once during translation — the same chunked, range-based reads the workers used — but it never materialises the mosaic in RAM. The one resource to watch is /tmp, where the VRT and the output COG live.

The overviews the COG driver builds are what make the merged product usable downstream: a map client zoomed out to a whole-country view reads a small overview level instead of decoding the full-resolution mosaic, and the same overview pyramid lets the next pipeline stage sample the product cheaply. The resampling method chosen for those overviews matters to correctness, not just speed. For continuous data — elevation, reflectance, NDVI — use average or bilinear so downsampled pixels represent the mean of their parents; for categorical rasters such as land-cover class codes, use nearest or mode, because averaging class integers produces meaningless intermediate values. The COG driver’s OVERVIEWS=AUTO picks a sensible default, but a categorical product almost always needs an explicit resampling override.

Prerequisites

  • Runtime: Python 3.11 on AWS Lambda with a GDAL 3.9+ layer exposing osgeo.gdal; rio-cogeo optional for validation. Build per Packaging & Dependency Management for Serverless GIS.
  • Ephemeral storage: raise the function’s /tmp allocation toward the 10 GB maximum if the product approaches multi-gigabyte size — see ephemeral storage limits in AWS Lambda.
  • Memory/timeout: a 2–3 GB tier gives enough vCPU for overview generation per the memory and CPU allocation model; allow up to the 15-minute ceiling for large mosaics.
  • IAM: s3:ListBucket and s3:GetObject on the tile prefix, s3:PutObject on the product prefix, scoped per IAM security boundaries for Cloud GIS.
  • Environment variables (applied before importing osgeo.gdal):
    code 
    GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR
CPL_VSIL_CURL_ALLOWED_EXTENSIONS=.tif,.tiff
GDAL_NUM_THREADS=ALL_CPUS
GDAL_DATA=/opt/share/gdal
PROJ_LIB=/opt/share/proj

Implementation

The handler enumerates tile outputs, builds a VRT over their /vsis3 paths, checks free /tmp, and translates to a COG with overviews. Tiles are never downloaded; only the VRT and the product COG occupy /tmp.

python 
#!/usr/bin/env python3
"""merge_tiles_to_cog.py — Lambda: fan-in per-tile outputs into one COG.

Builds a VRT that references tiles in place via /vsis3, then translates it to a
Cloud Optimized GeoTIFF with overviews, guarding the 10 GB /tmp ceiling.
"""
import os

import boto3
from osgeo import gdal

gdal.UseExceptions()
os.environ.setdefault("GDAL_DISABLE_READDIR_ON_OPEN", "EMPTY_DIR")
os.environ.setdefault("CPL_VSIL_CURL_ALLOWED_EXTENSIONS", ".tif,.tiff")
os.environ.setdefault("GDAL_NUM_THREADS", "ALL_CPUS")

s3 = boto3.client("s3")
VRT_PATH = "/tmp/mosaic.vrt"
COG_PATH = "/tmp/product.tif"


def _free_tmp_bytes(path: str = "/tmp") -> int:
    """Bytes currently free on the ephemeral mount."""
    st = os.statvfs(path)
    return st.f_bavail * st.f_frsize


def _list_tiles(bucket: str, prefix: str) -> list[str]:
    """Return /vsis3 paths for every tile under the run's output prefix."""
    paths = []
    for page in s3.get_paginator("list_objects_v2").paginate(Bucket=bucket, Prefix=prefix):
        for obj in page.get("Contents", []):
            if obj["Key"].endswith(".tif"):
                paths.append(f"/vsis3/{bucket}/{obj['Key']}")
    return sorted(paths)


def handler(meta, context):
    bucket, prefix, run_id = meta["output_bucket"], meta["out_prefix"], meta["run_id"]

    vsis3_tiles = _list_tiles(bucket, prefix)
    if not vsis3_tiles:
        raise RuntimeError(f"No tiles found under s3://{bucket}/{prefix}")

    # 1) Virtual mosaic — references tiles in place, costs kilobytes.
    gdal.BuildVRT(VRT_PATH, vsis3_tiles, options=gdal.BuildVRTOptions(resolution="highest"))

    # 2) Guard /tmp before materialising the product. Estimate the COG at the
    #    uncompressed extent of the VRT and require 20% headroom for overviews.
    vrt = gdal.Open(VRT_PATH)
    bytes_per_px = {"Byte": 1, "UInt16": 2, "Int16": 2, "Float32": 4}.get(
        gdal.GetDataTypeName(vrt.GetRasterBand(1).DataType), 4)
    est_bytes = vrt.RasterXSize * vrt.RasterYSize * vrt.RasterCount * bytes_per_px
    if est_bytes * 1.2 > _free_tmp_bytes():
        # Product would not fit /tmp — write straight to S3 via /vsis3 instead.
        out_target = f"/vsis3/{bucket}/runs/{run_id}/product.tif"
    else:
        out_target = COG_PATH

    # 3) Translate through the VRT block by block into a COG with overviews.
    gdal.Translate(
        out_target, VRT_PATH,
        options=gdal.TranslateOptions(
            format="COG",
            creationOptions=[
                "COMPRESS=DEFLATE", "PREDICTOR=2",
                "BLOCKSIZE=512", "OVERVIEWS=AUTO",
                "NUM_THREADS=ALL_CPUS",
            ],
        ),
    )

    out_key = f"runs/{run_id}/product.tif"
    if out_target == COG_PATH:  # local write → upload, then release /tmp
        s3.upload_file(COG_PATH, bucket, out_key)
        os.remove(COG_PATH)
    os.remove(VRT_PATH)

    return {
        "run_id": run_id,
        "product_bucket": bucket,
        "product_key": out_key,
        "tiles_merged": len(vsis3_tiles),
        "crs": meta["profile"]["crs"],
    }

To use rio-cogeo instead of the raw Translate call — its default profile writes web-optimised overviews and validates the result — swap step 3 for:

python 
from rio_cogeo.cogeo import cog_translate
from rio_cogeo.profiles import cog_profiles

cog_translate(
    VRT_PATH, COG_PATH, cog_profiles.get("deflate"),
    overview_level=5, overview_resampling="average",
    web_optimized=False, in_memory=False,  # in_memory=False keeps the mosaic off RAM
)

Verification

Confirm the product is a valid COG with correctly written overviews and internal tiling, and that the merge stayed within /tmp.

bash 
# Validate COG structure (overviews + internal tiling present)
rio cogeo validate "s3://$OUTPUT_BUCKET/runs/$RUN_ID/product.tif"

# Confirm overviews and block size via gdalinfo
gdalinfo "/vsis3/$OUTPUT_BUCKET/runs/$RUN_ID/product.tif" \
  | grep -E "Block=|Overviews"

Expected output for a merged 40000×40000 mosaic:

code 
/vsis3/geo-ingest/runs/scene/product.tif is a valid cloud optimized GeoTIFF
Band 1 Block=512x512 Type=Float32
  Overviews: 20000x20000, 10000x10000, 5000x5000, 2500x2500, 1250x1250

In the Lambda CloudWatch log, confirm Max Memory Used stayed well below the allocation (the VRT keeps it bounded) and that no No space left on device error appeared — the os.statvfs guard should have diverted an oversized product to /vsis3 output.

Gotchas and Edge Cases

  • Overlapping or gapped tiles produce seams. BuildVRT honours the tiles’ geotransforms; if the partitioning step emitted misaligned windows, the VRT will show black gaps or double-counted overlaps. Block-aligned windows from the partitioning stage prevent this; verify with gdalinfo that the VRT extent equals the source extent.

  • GDAL_NUM_THREADS=ALL_CPUS only helps at higher memory tiers. Overview generation is CPU-bound, and vCPU scales with the memory allocation. At 512 MB there is effectively one core, so multi-threaded translate gives no benefit — raise the tier if overview build dominates the runtime.

  • Writing a COG directly to /vsis3 needs a two-pass driver. The COG driver builds overviews after the main image, which requires seekable output; /vsis3 streaming writes can fail for very large products. Prefer local /tmp with the os.statvfs guard, and only fall back to /vsis3 output when the product genuinely will not fit.

  • Nodata mismatches merge incorrectly. If tile workers wrote differing nodata values, pass -srcnodata/-vrtnodata to BuildVRT so the mosaic treats them consistently; otherwise edge pixels bleed into the overviews.

  • A single merge Lambda can still overrun the 15-minute ceiling. Overview generation for a very large mosaic is CPU-bound and reads every tile once; on a 100 GB reassembled product it can exceed the timeout even at maximum vCPU. When it does, split the merge itself: build intermediate COGs per row-band of tiles, then a final VRT-over-COGs translate, or move the merge to a Fargate task via the same circuit-breaker fallback the recipe uses for oversized tiles.

  • List consistency across a large fan-out. list_objects_v2 is strongly consistent on S3, but if the merge starts before the last tile worker’s write completes it will silently omit tiles. Gate the merge on the Map state’s completion — which Step Functions guarantees — rather than polling the prefix on a timer.

Frequently Asked Questions

Do I need to download the tiles to merge them?

No. gdal.BuildVRT references the per-tile GeoTIFFs in place through /vsis3 paths, so the VRT is a small XML index rather than a copy of the pixels. Only the final COG is materialised on /tmp, and even that can be streamed to /vsis3 output if it approaches the 10 GB ephemeral ceiling.

Why build a VRT instead of merging with gdal_merge.py?

gdal_merge.py copies every pixel into a new raster before writing, which can exceed both the memory and /tmp ceilings for a large mosaic. A VRT is a virtual index that costs kilobytes; gdal.Translate then reads through it block by block, keeping peak memory bounded to a few blocks.

How do I keep the merge within the 10 GB /tmp limit?

Reference tiles via /vsis3 rather than downloading them, write only the VRT (kilobytes) and output COG to /tmp, and check free space with os.statvfs before translating. If the product would exceed /tmp, write it directly to /vsis3 output or raise the function’s ephemeral storage allocation toward 10 GB.

Back to Process a 10 GB GeoTIFF with AWS Step Functions and Lambda