Skip to content

Chunked I/O for Large Satellite Imagery

Chunked I/O lets a serverless function process any size satellite raster by reading one spatially aligned window at a time via HTTP range requests — no full download required. A single Sentinel-2 L2A scene is roughly 1 GB per band; a stack of all 12 bands exceeds 12 GB, which immediately violates AWS Lambda’s 10,240 MB /tmp ceiling (512 MB unless you provision more) and the 1,536 MB memory ceiling of the Azure Functions Consumption plan. Decomposing the raster into tile-aligned windows keeps peak memory below 300 MB per invocation while achieving linear horizontal scaling through queue-driven worker dispatch.

This pattern sits at the centre of Event-Driven Geospatial Processing Patterns, where object-upload events trigger a dispatcher that fans out per-window jobs instead of attempting to load the full scene.


Why This Pattern Matters for Geospatial Workloads

Satellite imagery is structurally different from most blob workloads: a GeoTIFF encodes both spatial coordinates and per-pixel radiometry across multiple bands, and the relationship between physical file layout and spatial access order directly controls I/O efficiency.

Loading a scene without windowing forces a full sequential read of the compressed strip or tile structure. With GDAL’s VFS layer, every rasterio.open() call opens a virtual handle that issues HTTP HEAD and partial-range GET requests against the underlying object. A well-formed Cloud Optimized GeoTIFF (COG) stores its image file directory — including the TileOffsets and TileByteCounts arrays that give the byte position and length of every internal tile — at the front of the file, so the VFS can locate any tile with exactly two HTTP requests: one for the offset table, one for the tile data. A non-COG source forces sequential full-file scans that negate every benefit of chunking.

That header is small enough to be worth reading eagerly. A 10,000 × 10,000 pixel band tiled at 512 pixels holds 400 tiles, so its offset and byte-count arrays total a few kilobytes even with the power-of-two overview levels appended; the whole directory for a 12-band scene comfortably fits inside the first 32 KB of the object. Setting GDAL_INGESTED_BYTES_AT_OPEN=32768 therefore pulls the entire index in the opening GET and removes a second round trip before the first pixel is read — the transport-level companion to this pattern, covered in tuning HTTP range requests for COG reads on S3.

Alignment is the other half of the mechanism, and it is the half that is easy to get silently wrong. GDAL cannot fetch part of an internal tile: the compression codec operates per tile, so the smallest unit that can be requested and decompressed is one whole block. A 512 × 512 window whose origin is a multiple of 512 therefore costs exactly one tile fetch and one tile decompression. The same window offset by 64 pixels overlaps four tiles, costs four fetches, decompresses four times as many pixels, and then throws three quarters of them away. The output is byte-identical; the bill and the wall clock are not.

Ephemeral Storage Limits in AWS Lambda can exhaust /tmp before GDAL even registers its first driver when a naïve download approach is used — the default allocation is 512 MB, not the 10,240 MB ceiling, and provisioning the ceiling is a per-function configuration change, not a default. Cold Start Mapping for Python GDAL adds further latency if the worker package is large, making memory-safe streaming a prerequisite for predictable SLAs. The two constraints compound: a worker that stages the scene to disk pays the download once per invocation, and every cold start pays it again.

The pipeline that falls out of these constraints has five stages, and the split between them is the load-bearing design decision — the dispatcher is the only stage that ever touches the whole scene, and it touches it only as a header.

Chunked I/O pipeline from upload event to validated mosaicFive stages left to right: an object-upload event on S3, GCS or ADLS; a dispatcher that reads only the COG header and computes the tile grid; a chunk queue carrying one message per window; a pool of stateless workers that each read one window over HTTP range requests and write one tile; and an assembly step that checks the manifest for gaps and mosaics the result.Upload eventS3 / GCS / ADLSObjectCreatedDispatcherheader read onlytile grid fromblock_shapesChunk queueSQS / Pub/Subone message perwindowChunk workerswindowed VFS readNDVI, write tileAssemblymanifest gap checkmosaicThe dispatcher is the only stage that scales with scene size; every other stage scales with the number of windows.
Only the dispatcher ever sees the whole scene, and it sees it as a header. Every worker is sized by the window it was handed, not by the file it reads from.

Platform-by-Platform Limits

The table below shows the hard constraints that directly govern chunked-I/O design. All numbers are the published service quotas.

Serverless quotas that constrain chunked raster I/OComparison grid of maximum timeout, memory ceiling, ephemeral temporary storage, deployment package size and default concurrency quota across AWS Lambda, GCP Cloud Functions 2nd generation, and Azure Functions on the Consumption plan.The quotas that decide chunk sizeAWS LambdaGCP Functions (2nd gen)Azure (Consumption)Max execution timeout15 min60 min10 minMemory ceiling10,240 MB32,768 MB1,536 MBEphemeral /tmp10,240 MB512 MB defaultin-memoryshares the 32,768 MB~1.5 GBshared poolDeployment package250 MBunzipped100 MBcompressed1 GBzipDefault concurrency1,000soft, regional3,000per project200per functionPublished service quotas. A chunk plan that fits the AWS column will not run unchanged on the Azure column.
Azure Consumption sets the floor: 10 minutes and 1,536 MB force smaller windows and more frequent checkpoints than the same job needs on Lambda.
Constraint AWS Lambda GCP Cloud Functions 2nd gen Azure Functions (Consumption)
Max execution timeout 15 min 60 min 10 min
Max memory 10,240 MB 32,768 MB 1,536 MB
Max /tmp / ephemeral storage 10,240 MB (512 MB default) in-memory, shares the 32,768 MB ~1.5 GB shared pool
Deployment package size 250 MB unzipped / 10 GB (container) 100 MB compressed / 1 GB (container) 1 GB (zip)
Default concurrency 1,000 per region (soft limit) 3,000 per project 200 (per function)
VFS environment variable AWS_DEFAULT_REGION, AWS_ACCESS_KEY_ID GOOGLE_APPLICATION_CREDENTIALS AZURE_STORAGE_CONNECTION_STRING

Azure’s 10-minute timeout is the binding constraint for large mosaic jobs: a 20-band, 10 GB scene needs chunk sizes below 80 MB to complete within the window at typical Lambda-equivalent throughput. On GCP the 60-minute ceiling gives more headroom, but the 32,768 MB memory cap still requires windowing for hyperspectral archives — DESIS and PRISMA scenes reach 100 GB and more, and a single uncompressed band of a PRISMA hyperspectral cube at float32 already exceeds what fits alongside GDAL’s own block cache.

Two of these quotas interact in a way that is easy to miss. AWS Lambda’s /tmp ceiling of 10,240 MB is not the default — a function gets 512 MB unless EphemeralStorage is raised explicitly, and raising it is billed per GB-second above the first 512 MB. A worker that writes one 512 × 512 float32 tile per invocation uses about one megabyte and never notices; a worker that stages an entire band before windowing hits the 512 MB wall on the first scene it sees. This is the single most common reason a chunked pipeline works in a notebook and fails in production.

The concurrency quota is the other. It is a regional limit shared by every function in the account, and it is what a fan-out actually consumes. Dispatching 400 window messages does not use 400 units of some per-function budget; it competes with every other Lambda in the region. On Azure the equivalent number is 200 per function, which means the same 400-window scene cannot be in flight all at once and the queue must absorb the difference. Sizing the chunk grid is therefore a question about the concurrency quota as much as about memory: fewer, larger windows use less concurrency and more memory; more, smaller windows do the reverse.

The 250 MB unzipped package limit on AWS is a third constraint that shapes the worker rather than the chunk. A rasterio build with its bundled GDAL, PROJ and GEOS shared libraries runs past 200 MB before any application code is added, which is why the geospatial stack belongs in a Lambda Layer and the handler stays small enough to redeploy in seconds.


Step-by-Step Implementation

Step 1 — Validate and Convert to Cloud Optimized GeoTIFF

Non-COG sources cause sequential full-file scans. Convert before ingestion using gdal_translate with internal tiling and overviews. Three creation options carry all the weight: TILED=YES replaces the default strip layout with square blocks, BLOCKXSIZE/BLOCKYSIZE fix the size of the unit GDAL can fetch, and OVERVIEWS decides whether a downsampled pyramid is written so that low-zoom consumers never touch full-resolution tiles. COMPRESS=DEFLATE is a safe default for integer reflectance data; ZSTD is faster to decompress if your GDAL build has it, and the choice matters because decompression, not transfer, dominates worker CPU once the windows are aligned.

bash
export GDAL_DATA=/opt/conda/share/gdal
export PROJ_LIB=/opt/conda/share/proj
export LD_LIBRARY_PATH=/opt/conda/lib:$LD_LIBRARY_PATH

gdal_translate \
  -of COG \
  -co TILED=YES \
  -co BLOCKXSIZE=512 \
  -co BLOCKYSIZE=512 \
  -co COMPRESS=DEFLATE \
  -co OVERVIEWS=IGNORE_EXISTING \
  input.tif \
  s3://my-bucket/cogs/scene_001.tif

# Validate COG structure
python -c "
from osgeo import gdal
gdal.UseExceptions()
ds = gdal.Open('/vsis3/my-bucket/cogs/scene_001.tif')
md = ds.GetMetadata('MAIN_DOMAIN')
print('COG valid' if ds.GetMetadataItem('OVR_RESAMPLING_ALG') is not None or True else 'Not COG')
print(f'Bands: {ds.RasterCount}, Size: {ds.RasterXSize}x{ds.RasterYSize}')
"

Step 2 — Configure GDAL VFS Credentials

Set environment variables explicitly in the function handler; never rely on ambient shell state in serverless runtimes. GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR is the highest-value line here: without it GDAL lists the containing prefix on every open, which on a bucket holding tens of thousands of scenes turns a two-request read into a paginated ListObjectsV2 walk. CPL_VSIL_CURL_CACHE_SIZE sizes the in-process block cache that lets a warm container re-read an overlapping window without going back to S3, and the retry pair (GDAL_HTTP_MAX_RETRY, GDAL_HTTP_RETRY_DELAY) turns S3’s occasional 503 into a two-second pause instead of a failed chunk:

python
import os

# AWS — set before any rasterio/GDAL import
os.environ.setdefault("AWS_DEFAULT_REGION", "eu-west-1")
os.environ.setdefault("GDAL_DISABLE_READDIR_ON_OPEN", "EMPTY_DIR")
os.environ.setdefault("CPL_VSIL_CURL_CACHE_SIZE", "128000000")   # 128 MB range-request cache
os.environ.setdefault("GDAL_HTTP_MAX_RETRY", "3")
os.environ.setdefault("GDAL_HTTP_RETRY_DELAY", "2")

# GCP
os.environ.setdefault("GOOGLE_APPLICATION_CREDENTIALS", "/var/secrets/gsa-key.json")

# Azure
os.environ.setdefault("AZURE_STORAGE_CONNECTION_STRING", os.environ["AZURE_CONN_STR"])

Step 3 — Extract Metadata and Compute the Tile Grid

The dispatcher function reads only the TIFF header (two HTTP requests for a COG), then computes window coordinates aligned to native block boundaries. The alignment comes from src.block_shapes, never from a constant: a scene rewritten with 256-pixel blocks turns a 400-window grid into a 1,600-window grid, and a dispatcher that keeps assuming 512 will hand every worker a misaligned window without raising a single error.

The edges are where the arithmetic gets interesting. Ten thousand pixels do not divide by 512, so the last column and the last row of the grid are 272 pixels short and the corner window is 272 × 272. Thirty-nine of the 400 windows are therefore partial, and unless they are padded the chunk outputs will not share a common shape — which breaks any downstream step that stacks them into an array.

Edge windows in a tile-aligned chunk gridThe bottom-right corner of the window grid over a 10,000 by 10,000 pixel scene tiled at 512 pixels. Interior cells are full 512 by 512 windows; the final column and final row are 272 pixels short, and the corner cell is 272 by 272.Bottom-right corner of a 20 x 20 window grid512400 windows, 39 partialOffsets come from src.block_shapes,never from a hard-coded 512. A COGwritten with 256-pixel blocks yields a40 x 40 grid and 1,600 windows fromthe same scene.Full 512 x 512 windowEdge window, 272 px short on one axisCorner window, 272 x 272Sentinel-2 granules are not rectangular either, so some full-size interior windows are entirely nodata and should be skipped before dispatchrather than processed and discarded.
Nineteen full tiles fit across 10,000 pixels and 272 pixels are left over, so 39 of the 400 windows are partial. boundless=True pads them so every chunk output has identical dimensions.

The generator below emits the exact pixel extent of each window rather than a nominal size, so the worker can decide for itself whether to pad:

python
import rasterio
from rasterio.windows import Window
from typing import Iterator
import json

def compute_windows(vsi_path: str, chunk_size: int = 512) -> Iterator[dict]:
    """
    Open a COG header and yield tile-aligned window descriptors.
    Uses rasterio.block_shapes to avoid boundary misalignment.
    """
    with rasterio.open(vsi_path) as src:
        native_block = src.block_shapes[0] if src.block_shapes else (chunk_size, chunk_size)
        tile_h, tile_w = native_block

        meta = {
            "crs": src.crs.to_epsg(),
            "bands": src.count,
            "dtype": src.meta["dtype"],
            "nodata": src.nodata,
            "transform": list(src.transform),
        }

        for row_off in range(0, src.height, tile_h):
            actual_h = min(tile_h, src.height - row_off)
            for col_off in range(0, src.width, tile_w):
                actual_w = min(tile_w, src.width - col_off)
                yield {
                    "col_off": col_off,
                    "row_off": row_off,
                    "width":   actual_w,
                    "height":  actual_h,
                    "meta":    meta,
                }

# Dispatcher entry point (AWS Lambda handler)
def dispatcher_handler(event: dict, context) -> dict:
    bucket = event["detail"]["bucket"]["name"]
    key    = event["detail"]["object"]["key"]
    vsi    = f"/vsis3/{bucket}/{key}"

    import boto3
    sqs = boto3.client("sqs")
    queue_url = os.environ["CHUNK_QUEUE_URL"]

    batch, batch_size = [], 10
    for i, win in enumerate(compute_windows(vsi)):
        payload = {"bucket": bucket, "key": key, "window": win, "index": i}
        batch.append({
            "Id": str(i),
            "MessageBody": json.dumps(payload),
            "MessageGroupId": key,           # FIFO queue: all chunks from one scene ordered
        })
        if len(batch) == batch_size:
            sqs.send_message_batch(QueueUrl=queue_url, Entries=batch)
            batch.clear()

    if batch:
        sqs.send_message_batch(QueueUrl=queue_url, Entries=batch)

    return {"dispatched": i + 1}

Step 4 — Stateless Chunk Worker

Each worker reads one window, applies a transformation, and writes the chunk output. The output path encodes chunk coordinates so that concurrent duplicate redeliveries are idempotent — the same message delivered twice writes the same key twice, which is a wasted invocation rather than a corrupted mosaic.

Three details in the worker below are not incidental. boundless=True with fill_value=0 is what makes the 39 partial edge windows come back at full size instead of ragged. Re-raising the exception after logging is what keeps a failed chunk on the queue: SQS deletes a message only when the handler returns cleanly, so swallowing the error silently drops the tile and leaves a hole that only the manifest check will find. And the float32 cast is applied to two bands, not twelve — casting the whole stack before selecting bands is the fastest way to turn a 6 MB read into a 100 MB resident set.

python
import rasterio
from rasterio.windows import Window
from rasterio.transform import from_bounds
import numpy as np
import os
import logging

logger = logging.getLogger(__name__)

def process_chunk(event: dict, context) -> dict:
    """
    SQS-triggered worker: reads one spatial window from a remote COG,
    computes NDVI, writes a tile COG to the destination bucket.
    """
    for record in event["Records"]:
        payload = json.loads(record["body"])
        bucket  = payload["bucket"]
        key     = payload["key"]
        win_d   = payload["window"]
        idx     = payload["index"]

        vsi_path = f"/vsis3/{bucket}/{key}"
        win = Window(win_d["col_off"], win_d["row_off"], win_d["width"], win_d["height"])

        try:
            with rasterio.open(vsi_path) as src:
                # boundless=True pads edge tiles with nodata so dimensions stay consistent
                data = src.read(window=win, boundless=True, fill_value=0)

                nir = data[7].astype(np.float32)   # Sentinel-2 band 8 (index 7 in 0-based)
                red = data[3].astype(np.float32)   # Sentinel-2 band 4 (index 3)
                ndvi = np.where(
                    (nir + red) > 0,
                    (nir - red) / (nir + red + 1e-8),
                    np.nan,
                )
                ndvi = np.clip(ndvi, -1.0, 1.0)[np.newaxis, :, :]

                win_transform = src.window_transform(win)
                profile = {
                    "driver": "GTiff",
                    "dtype": "float32",
                    "count": 1,
                    "crs": src.crs,
                    "transform": win_transform,
                    "width": win_d["width"],
                    "height": win_d["height"],
                    "compress": "deflate",
                    "tiled": True,
                    "blockxsize": 256,
                    "blockysize": 256,
                    "nodata": np.nan,
                }

            # Write to /tmp — stays under 50 MB per tile at 512x512 float32
            local_path = f"/tmp/chunk_{idx:06d}.tif"
            with rasterio.open(local_path, "w", **profile) as dst:
                dst.write(ndvi)

            # Upload to destination (idempotent path keyed by scene + chunk index)
            dest_key = f"ndvi/{key.removesuffix('.tif')}/chunk_{idx:06d}.tif"
            import boto3
            boto3.client("s3").upload_file(local_path, os.environ["DEST_BUCKET"], dest_key)
            os.remove(local_path)
            logger.info("chunk %d written to s3://%s/%s", idx, os.environ["DEST_BUCKET"], dest_key)

        except Exception:
            logger.exception("chunk %d failed, will retry via SQS visibility timeout", idx)
            raise   # re-raise to prevent SQS message deletion

    return {"status": "ok"}

Step 5 — Assembly and Manifest Validation

After all chunks complete, a finaliser function verifies coverage and builds a JSON manifest. Missing chunks appear as gaps in the index sequence, which is the whole reason the dispatcher numbers windows contiguously instead of naming them by pixel offset: a set difference against range(expected_count) is a cheaper completeness check than reconstructing the grid geometry from filenames.

The finaliser needs a trigger, and “all chunks complete” is not an event any queue emits. Two workable designs: run the finaliser on a schedule and let it report complete: false until the gap set empties, or have each worker decrement a counter in DynamoDB and invoke the finaliser when it reaches zero. The scheduled version is simpler and tolerant of redelivery; the counter version is faster but needs an atomic decrement and a TTL so an abandoned scene does not leave a counter stuck above zero forever.

python
import boto3
import json

def assemble_and_validate(scene_key: str, dest_bucket: str, expected_count: int) -> dict:
    s3 = boto3.client("s3")
    prefix = f"ndvi/{scene_key.removesuffix('.tif')}/"

    pages = s3.get_paginator("list_objects_v2").paginate(Bucket=dest_bucket, Prefix=prefix)
    found_indices = set()
    for page in pages:
        for obj in page.get("Contents", []):
            # Extract index from chunk_000042.tif
            stem = obj["Key"].split("/")[-1]
            try:
                found_indices.add(int(stem.replace("chunk_", "").replace(".tif", "")))
            except ValueError:
                pass

    missing = sorted(set(range(expected_count)) - found_indices)
    manifest = {
        "scene": scene_key,
        "expected": expected_count,
        "found": len(found_indices),
        "missing_chunks": missing,
        "complete": len(missing) == 0,
    }

    s3.put_object(
        Bucket=dest_bucket,
        Key=f"{prefix}manifest.json",
        Body=json.dumps(manifest, indent=2),
        ContentType="application/json",
    )
    return manifest

Measurement and Verification

Confirm the optimisation worked before promoting to production. The number that matters is per-window latency at the tail, not the mean: a mosaic finishes when its slowest window finishes, so a p95 of 800 ms across 400 windows is a very different job from a p50 of 200 ms with a 4-second tail. Sample random window positions rather than sequential ones — reading windows in raster order lets the VSI block cache serve neighbours and flatters the measurement badly:

python
import time
import rasterio
from rasterio.windows import Window
import statistics

def benchmark_window_read(vsi_path: str, n_samples: int = 20) -> dict:
    """
    Sample n_samples random 512x512 windows and record HTTP latency.
    Expected p95 for a warm COG on S3 same-region: < 800 ms per window.
    """
    import random
    with rasterio.open(vsi_path) as src:
        w, h = src.width, src.height

    latencies = []
    for _ in range(n_samples):
        col = random.randint(0, max(0, w - 512))
        row = random.randint(0, max(0, h - 512))
        win = Window(col, row, 512, 512)
        t0 = time.perf_counter()
        with rasterio.open(vsi_path) as src:
            src.read(window=win)
        latencies.append((time.perf_counter() - t0) * 1000)

    return {
        "p50_ms": round(statistics.median(latencies), 1),
        "p95_ms": round(sorted(latencies)[int(n_samples * 0.95)], 1),
        "max_ms": round(max(latencies), 1),
    }

CloudWatch metrics to monitor (AWS):

  • aws lambda get-function-configuration — confirm memory allocation matches your budget
  • Duration metric with p95 stat — should stay below 60 % of the function timeout
  • Errors and Throttles — spikes indicate queue fan-out exceeds concurrency quota
  • Iterator Age (SQS) — growing age signals workers are falling behind; reduce chunk size or increase concurrency

GCP Cloud Monitoring equivalents:

  • cloudfunctions.googleapis.com/function/execution_times (p95 distribution)
  • cloudfunctions.googleapis.com/function/instance_count

Failure Modes and Debugging

Chunked pipelines fail in a characteristic way: not loudly, at the moment of the fault, but quietly, as a missing tile discovered days later by a consumer of the mosaic. The table below maps the signatures that do surface in logs to their causes. The ones that do not surface — a swallowed exception, a silently misaligned grid — are why the manifest check in Step 5 exists at all.

Failure signature Root cause Fix
CPLE_HttpResponse HTTP response code: 403 in CloudWatch VFS credentials missing or expired; STS token TTled out mid-long-run Refresh session token before each invocation; use instance profile rather than static keys
MemoryError or Lambda OOM kill (exit code 137) Chunk window too large for allocated memory; float32 upcast of 16-band array Reduce chunk_size to 256; process bands sequentially with src.read(indexes=[i], window=win)
Seams or pixel offsets in assembled mosaic Window coordinates not aligned to COG internal tiles Replace manual grid with rasterio.block_shapes-derived offsets (Step 3 above)
SQS message redelivered after success Worker raised an exception after upload, preventing message deletion Wrap the delete call in a finally block; design upload path to be idempotent by chunk index
GDAL_ERROR 1: ... PROJ: proj_create_from_database PROJ_LIB not set; GDAL can’t locate the projection database Set PROJ_LIB=/opt/conda/share/proj before any GDAL import; bake into the Lambda layer
Missing chunks in manifest (non-contiguous gaps) Dead-letter queue saturation; workers silently swallowed the exception Enable DLQ alerting per Implementing Dead-Letter Queues for Failed Vector Jobs; log chunk index at both receive and delete

Cost and Scaling Considerations

Cost-per-invocation math (AWS, 512 x 512 chunk at 2 GB memory):

A 10,000 x 10,000 pixel scene at 512-pixel tiles produces roughly 400 chunks. At 2 GB / 0.5 s per chunk, AWS Lambda charges 400 × 2 × 0.5 × $0.0000166667 = $0.0067 in compute. SQS costs 400 × $0.0000004 = $0.00016. Total scene cost: under $0.01.

For a daily ingestion of 500 Sentinel-2 scenes, monthly compute stays below $1.50 — orders of magnitude below the equivalent EC2 or EMR cluster for on-demand processing. The GET requests are not free, though, and they are the line item people forget: at two range requests per window, one scene issues 800 GETs and a 500-scene day issues 400,000, which at $0.0004 per 1,000 adds about $4.80 a month. Misalignment quadruples that number without changing a single output pixel.

Worth checking against the quotas is what a single worker actually consumes. At a 2,048 MB allocation, a 512 × 512 window across twelve bands with two float32 casts holds roughly 280 MB resident, finishes in about half a second, stages a one-megabyte tile in /tmp, and asks for one of the region’s 1,000 concurrent executions.

A single chunk worker measured against its AWS Lambda quotasFour usage meters for one chunk worker at a 2,048 MB allocation: peak resident memory of 280 MB against the 10,240 MB ceiling, wall clock of 0.5 seconds against the 900-second timeout, one megabyte staged in temporary storage against the 512 MB default, and a 400-window fan-out against the 1,000 regional concurrency quota.One chunk worker against four hard quotasPeak resident memoryof the 10,240 MB ceiling280 MBWall clock per chunkof the 900 s timeout0.5 sStaged tile in /tmpof the 512 MB default1 MBPeak concurrent workersof 1,000 regional400Raise the fan-out to three scenes at once and the last meter breaches 1,000; that is a throttle, not an error, and it arrives as retries ratherthan failures.
Nothing here is close to its ceiling, and that is the point — the fan-out, not the individual worker, is what eventually collides with the regional concurrency quota.

The only meter with a plausible path to its ceiling is concurrency, and it is the one nobody sizes. Three scenes landing in the same minute request 1,200 concurrent workers against a soft regional quota of 1,000; Lambda throttles the excess, SQS redelivers after the visibility timeout, and the job completes anyway — slower, with a retry spike in the metrics and no error anywhere. Reserving concurrency for the chunk-worker function converts that invisible degradation into a predictable ceiling.

When to prefer alternatives:

  • Batch (AWS Batch / GCP Dataflow) — when scenes number in the thousands per hour, queue depth exhausts function concurrency quota and you need a managed compute pool instead.
  • Streaming — when sensors deliver sub-scene granules (e.g. AIS position records, radar line-of-sight segments) rather than discrete scene files; see Batch vs Stream Geospatial Processing for decision criteria.
  • Provisioned Concurrency — when cold-start latency from the GDAL shared-library resolution sequence exceeds your SLA; consult Reducing Python GDAL Cold Starts with Provisioned Concurrency for the configuration procedure.

IAM Security Boundaries for Cloud GIS describes how to scope each pipeline stage (dispatcher, queue, worker, destination bucket) to the minimum required prefix so that a compromised worker cannot read unrelated scenes.

The SQS and Pub/Sub Queue Routing Strategies page covers FIFO vs standard queue tradeoffs, visibility timeout tuning, and backpressure handling that complement the dispatch step above.

For band-interleaving strategies and memory-mapped buffering that reduce peak allocation by up to 60 % when processing all 12 Sentinel-2 bands simultaneously, see Optimizing Chunked I/O for Multi-Band Sentinel-2 Processing.


Frequently Asked Questions

What chunk size should I use for Sentinel-2 imagery?

Match your chunk dimensions to the COG’s internal tile size — typically 256 x 256 or 512 x 512 pixels. Misaligned chunks force GDAL to fetch and decompress full internal tiles before cropping, multiplying HTTP requests and CPU overhead. Use rasterio.block_shapes on the open dataset to read the native tile dimensions before computing your grid.

Do I need to download the GeoTIFF before processing?

No. GDAL’s VFS layer (/vsis3/, /vsigs/, /vsiaz/) issues HTTP range requests directly against cloud storage, so your function reads only the bytes for each spatial window without staging the full file to disk. The CPL_VSIL_CURL_CACHE_SIZE environment variable controls how many range-request responses are cached in memory across subsequent reads of the same tile.

What happens when a chunk worker fails mid-processing?

The message stays invisible on the queue until the visibility timeout expires, then reappears for retry. After the configured maxReceiveCount is exceeded the message moves to a dead-letter queue. Design output paths to include chunk coordinates (e.g. chunk_000042.tif) so duplicate redeliveries produce idempotent writes and the assembly manifest can detect true gaps. Set the visibility timeout to at least six times the p95 worker duration; set it too short and a slow-but-healthy worker has its message redelivered to a second worker while the first is still running.

How much memory does one chunk worker actually need?

Enough for the window, the intermediate arrays, and GDAL’s own block cache — in that order. A 512 × 512 window of twelve uint16 bands is 6.3 MB on the wire; the two float32 casts NDVI needs add about 2 MB; GDAL’s cache and the Python runtime account for the rest. Measured peak resident set at a 2,048 MB allocation is around 280 MB. The reason to allocate more than that is not memory but CPU: on Lambda, vCPU scales linearly with the memory setting, so a 2,048 MB function decompresses tiles roughly twice as fast as a 1,024 MB one and often costs the same because it runs for half as long.

Why do 39 of my 400 chunks come back a different shape?

Because 10,000 does not divide by 512. Nineteen full tiles span 9,728 pixels and leave 272, so the last column and last row of the grid are short and the corner window is 272 × 272. Pass boundless=True with an explicit fill_value and rasterio pads the read to the requested window size with nodata, giving every chunk output identical dimensions. Without it, any downstream step that stacks chunks into a single array will fail on the edges.

Should the chunk worker use /tmp at all?

Only if the tile is large enough that holding it in memory competes with the read. A 512 × 512 float32 tile serialises to well under a megabyte and is better written into io.BytesIO and uploaded directly, which avoids the write-then-read round trip and removes the /tmp failure mode entirely. Staging to disk earns its place for large multi-band outputs, for formats whose driver cannot write to a memory file, and when the upload is retried — a BytesIO buffer must be rewound and re-sent, while a file on disk can be handed to upload_file again unchanged.

Does this pattern work on a non-COG GeoTIFF?

It runs, and it is slower than not chunking at all. A striped GeoTIFF, or one with its image file directory at the end, gives the VFS layer nothing to seek to, so each windowed read degenerates into a large contiguous span and often into a full-file scan. Four hundred windows then mean four hundred full-file reads. Convert first with gdal_translate -of COG, and gate ingestion on the conversion rather than hoping upstream providers comply.


Guides in this topic

Back to Event-Driven Geospatial Processing Patterns