Batch and Async Orthomosaic Processing
A single drone flight over one field is a laptop-scale problem. A cooperative running weekly flights over four hundred fields across a growing season is not — you are now processing thousands of multi-gigabyte orthomosaics, each needing stitching, cloud masking, and a stack of vegetation indices, with new imagery landing faster than any single process can consume it. This guide lays out the architecture for orchestrating multi-flight, multi-field orthomosaic and index jobs at scale in Python: a task queue that separates I/O-bound stitching from CPU-bound band math, asyncio for concurrent object-storage reads, Dask for chunked per-tile parallelism, and the idempotency, retry, and back-pressure machinery that keeps a run of ten thousand tiles from collapsing when the fiftieth upload times out. The concrete output is a pipeline that turns a directory of raw flight captures into a validated set of Cloud Optimized GeoTIFF index rasters, tracked by a job manifest that tells you exactly which tiles succeeded, which failed, and which are still in flight.
This page is part of the Drone Imagery Processing & Vegetation Index Workflows guide. See that page for the full context from ingest through temporal aggregation; here we focus purely on running those single-scene workflows at fleet scale.
Prerequisites
Python packages (exact versions tested):
celery==5.3.6redis==5.0.4dask==2024.5.0distributed==2024.5.0rasterio>=1.3.9aioboto3==12.4.0xarray==2024.5.0numpy>=1.24
pip install celery==5.3.6 redis==5.0.4 dask==2024.5.0 distributed==2024.5.0 \
rasterio==1.3.10 aioboto3==12.4.0 xarray==2024.5.0 numpy==1.26.4
Infrastructure assumptions:
- A reachable Redis instance (broker and result backend). A single node is fine for hundreds of concurrent tasks; use Redis Cluster or a managed queue beyond that.
- S3-compatible object storage holding raw captures and receiving COG outputs. Tiles are read as Cloud Optimized GeoTIFFs so that windowed HTTP range reads work without downloading the whole file.
Input data requirements:
- Per-flight orthomosaics (or the raw captures that produce them) in a projected CRS — always carry an explicit EPSG code (e.g.
EPSG:32615for UTM zone 15N). Mixing geographic and projected tiles in one batch is the most common cause of silent misalignment downstream. - Consistent band order across flights from a single sensor family — MicaSense RedEdge-MX (Blue 475, Green 560, Red 668, RedEdge 717, NIR 840 nm), DJI P4 Multispectral, or Parrot Sequoia. Record the sensor in the manifest so band indices are never guessed.
- Nodata value set in every source profile so masked pixels do not leak into index math.
1. Architecture & Concurrency Model
The defining property of a fleet-scale imagery pipeline is that its two dominant workloads are unlike each other. Stitching and object-storage transfer are I/O-bound: a worker spends most of its wall-clock time waiting on the network while a 4 GB orthomosaic uploads or a range read completes. Band math — NDVI, NDRE, SAVI across every tile — is CPU-bound and memory-heavy: the worker is pinned at 100% on the floating-point ratio while holding several bands of a tile in RAM. Running both on the same pool is the original sin of these systems. A prefork pool sized to physical cores will sit idle waiting on uploads; a large thread pool tuned for I/O will thrash the CPU under band math. The fix is to route the two workloads to separate named queues with separate worker pools, each sized and configured for its bottleneck.
Three layers of concurrency stack on top of that split:
Task-level (Celery + Redis). The unit of distribution is one task per flight for stitching, and one task per tile for index computation. Redis is the broker that holds the queues and applies back-pressure; a chord or group aggregates the per-tile results back into a per-field summary. This is the coarse-grained layer that spreads work across machines, and it is covered end-to-end in orchestrating multi-flight jobs with Celery and Redis.
Connection-level (asyncio). Within a single stitching or read task, dozens of COG windows must be pulled from object storage. Issuing those range reads serially wastes the network; asyncio with aioboto3 fires them concurrently and awaits the batch, cutting the read phase of a multi-tile mosaic from minutes to seconds. This is concurrency without parallelism — one thread, many in-flight requests — and it is exactly right for the I/O-bound half of the pipeline.
Array-level (Dask). Inside a band-math task, a single orthomosaic is still too large to hold in memory. Dask chunks the raster into blocks, builds a lazy task graph over them, and computes NDVI block by block under a hard memory cap, spilling to disk only when necessary. This fine-grained parallelism is the subject of parallel tile processing with Dask.
Why the layering matters agronomically: a batch that silently drops or corrupts tiles produces a vegetation-index mosaic with holes, and a hole in an NDVI raster is indistinguishable from genuinely bare or stressed ground once it reaches a management-zone map. The orchestration layer’s job is not speed for its own sake — it is completeness under failure, so that every field’s prescription is built from every tile that was flown.
2. Step-by-Step Implementation
Step 1 — Define the job manifest
The manifest is the source of truth for a run. It enumerates every unit of work, carries the CRS and sensor metadata that tasks must not guess, and records the state of each tile so a re-run resumes instead of restarting. Model it as a plain dataclass serialised to JSON (or a row per tile in a database for large fleets).
from dataclasses import dataclass, field, asdict
from enum import Enum
import json
class TileState(str, Enum):
PENDING = "pending"
RUNNING = "running"
DONE = "done"
FAILED = "failed"
@dataclass
class TileJob:
flight_id: str
field_id: str
tile_id: str
src_uri: str # s3://bucket/raw/flight/tile.tif
dst_uri: str # s3://bucket/ndvi/flight/tile.tif
epsg: int # explicit CRS, e.g. 32615 — never inferred
sensor: str # "micasense_rededge_mx"
state: TileState = TileState.PENDING
attempts: int = 0
@dataclass
class JobManifest:
run_id: str
tiles: list = field(default_factory=list)
def to_json(self, path: str) -> None:
with open(path, "w") as fh:
json.dump(
{"run_id": self.run_id, "tiles": [asdict(t) for t in self.tiles]},
fh, indent=2, default=str,
)
# Build a manifest for two flights, guarding against a missing CRS
tiles = [
TileJob("F001", "north-80", "t0001",
"s3://ag-raw/F001/t0001.tif", "s3://ag-ndvi/F001/t0001.tif",
epsg=32615, sensor="micasense_rededge_mx"),
TileJob("F002", "south-40", "t0001",
"s3://ag-raw/F002/t0001.tif", "s3://ag-ndvi/F002/t0001.tif",
epsg=32615, sensor="micasense_rededge_mx"),
]
for t in tiles:
assert t.epsg is not None, f"Tile {t.tile_id} has no EPSG — refuse to enqueue"
manifest = JobManifest(run_id="2026-07-12-weekly", tiles=tiles)
manifest.to_json("run_manifest.json")
print(f"Manifest: {len(manifest.tiles)} tiles enqueued")
Step 2 — Configure separate queues for I/O and CPU work
Two queues, two pools. The Celery app declares the routing; the workers are launched with different pool implementations and concurrency levels. Prefetch is set to 1 on the compute worker so a single greedy worker cannot hoard tiles it will process slowly.
from celery import Celery
app = Celery(
"ortho",
broker="redis://localhost:6379/0",
backend="redis://localhost:6379/1",
)
app.conf.update(
task_acks_late=True, # ack only after success — safe redelivery
task_reject_on_worker_lost=True, # requeue if a worker is killed mid-task
worker_prefetch_multiplier=1, # no hoarding of slow compute tasks
task_default_queue="cpu",
task_routes={
"ortho.stitch_flight": {"queue": "io"},
"ortho.compute_index_tile": {"queue": "cpu"},
},
broker_transport_options={"visibility_timeout": 3 * 3600}, # > longest task
)
Launch the two pools as separate processes (bash):
# I/O-bound stitching: high concurrency, green threads
celery -A ortho worker -Q io --pool gevent --concurrency 50 -n io@%h
# CPU-bound band math: one process per physical core
celery -A ortho worker -Q cpu --pool prefork --concurrency 8 -n cpu@%h
Step 3 — Read COG windows concurrently with asyncio
Before a tile can be stitched or index-computed, its windows must come off object storage. aioboto3 lets one task issue many range reads concurrently. The pattern below fetches a batch of byte ranges for the tiles of one flight and returns them in order, bounding concurrency with a semaphore so a single task cannot open a thousand sockets at once.
import asyncio
import aioboto3
async def _fetch_range(session, bucket, key, byte_range, sem):
async with sem:
async with session.client("s3") as s3:
resp = await s3.get_object(Bucket=bucket, Key=key, Range=byte_range)
return await resp["Body"].read()
async def fetch_windows(bucket, keys_ranges, max_concurrency=16):
"""Concurrently read (key, byte_range) pairs from S3-compatible storage."""
sem = asyncio.Semaphore(max_concurrency)
session = aioboto3.Session()
tasks = [
_fetch_range(session, bucket, key, rng, sem)
for key, rng in keys_ranges
]
return await asyncio.gather(*tasks)
# Example: pull the header ranges of 200 tiles at once
keys = [(f"raw/F001/t{ i:04d}.tif", "bytes=0-65535") for i in range(200)]
chunks = asyncio.run(fetch_windows("ag-raw", keys, max_concurrency=16))
assert len(chunks) == 200, "Lost a window read — check retry policy"
print(f"Fetched {len(chunks)} COG header windows concurrently")
Step 4 — Compute the tile index under bounded memory
The band-math task opens the tile, computes NDVI over Dask-chunked arrays, and writes to a temporary key. The full Dask mechanics — map_blocks, lazy compute, memory caps — are detailed in parallel tile processing with Dask; the excerpt below is the task body that the queue drives.
import numpy as np
import rasterio
from rasterio.windows import Window
RED_BAND, NIR_BAND = 3, 5 # MicaSense RedEdge-MX band order
def compute_ndvi_windowed(src_path, dst_path, epsg, chunk=1024):
with rasterio.open(src_path) as src:
assert src.crs is not None, "Source has no CRS"
assert src.crs.to_epsg() == epsg, (
f"Manifest EPSG {epsg} != raster EPSG {src.crs.to_epsg()}"
)
profile = src.profile.copy()
profile.update(count=1, dtype="float32", nodata=np.nan,
tiled=True, blockxsize=256, blockysize=256, compress="deflate")
with rasterio.open(dst_path, "w", **profile) as dst:
for row in range(0, src.height, chunk):
for col in range(0, src.width, chunk):
win = Window(col, row,
min(chunk, src.width - col),
min(chunk, src.height - row))
red = src.read(RED_BAND, window=win).astype("float32")
nir = src.read(NIR_BAND, window=win).astype("float32")
with np.errstate(invalid="ignore", divide="ignore"):
ndvi = (nir - red) / (nir + red)
dst.write(ndvi.astype("float32"), 1, window=win)
return dst_path
Step 5 — Make the task idempotent with atomic writes and retries
An idempotent task can be redelivered any number of times and converge on the same output. Write to a temporary URI keyed on the task, validate, then rename atomically. On redelivery, short-circuit if the validated output already exists.
import os
from celery import shared_task
from celery.exceptions import Retry
@shared_task(
bind=True, name="ortho.compute_index_tile",
autoretry_for=(IOError, OSError),
retry_backoff=True, retry_backoff_max=600,
retry_jitter=True, max_retries=5, acks_late=True,
)
def compute_index_tile(self, tile: dict):
dst = tile["dst_uri"]
tmp = f"{dst}.{self.request.id}.tmp"
# Idempotency guard: already done on a previous delivery
if _validated_output_exists(dst):
return {"tile_id": tile["tile_id"], "state": "done", "cached": True}
compute_ndvi_windowed(tile["src_uri"], tmp, tile["epsg"])
if not _validate_index(tmp, tile["epsg"]):
_remove(tmp)
raise self.retry(countdown=30) # transient corruption — try again
_atomic_rename(tmp, dst) # only now is the tile officially done
return {"tile_id": tile["tile_id"], "state": "done", "cached": False}
Step 6 — Fan out with a group and aggregate with a chord
One task per tile, aggregated into a per-field completion record. The chord’s callback runs only after every tile task in the group returns, giving a natural join point to update the manifest.
from celery import group, chord
def run_field(field_tiles: list):
header = group(
compute_index_tile.s(t) for t in field_tiles
)
callback = summarise_field.s(field_id=field_tiles[0]["field_id"])
return chord(header)(callback)
@shared_task(name="ortho.summarise_field")
def summarise_field(results, field_id):
done = sum(1 for r in results if r["state"] == "done")
return {"field_id": field_id, "tiles_done": done, "tiles_total": len(results)}
3. Key Parameters & Tuning
| Parameter | Type | Default | Agronomic Effect |
|---|---|---|---|
worker_prefetch_multiplier |
int | 1 | Keep at 1 for long band-math tasks so one worker cannot reserve dozens of tiles and stall a whole field’s completion. Raise to 4 only on the short I/O queue where tasks finish in under a second. |
visibility_timeout (Redis) |
int (s) | 10800 | Must exceed the longest task. A stitching task that outruns the timeout is redelivered and runs twice, doubling load and risking a duplicate write during a peak-flight backlog. |
max_concurrency (asyncio) |
int | 16 | Concurrent range reads per task. Too high saturates the storage endpoint and triggers throttling that manifests as failed tiles; 12–20 is safe for most S3-compatible backends. |
chunk / Dask block size |
int (px) | 1024 | Per-tile compute block. Smaller blocks cut peak RAM but add scheduler overhead; 512 for memory-tight nodes, 2048 for wide fields on high-memory workers. |
max_retries |
int | 5 | Retries for transient upload or read failures. Too low drops tiles on brief network blips, leaving NDVI holes; too high masks a genuinely bad source file for hours. |
retry_backoff_max |
int (s) | 600 | Caps exponential backoff so a struggling storage endpoint gets breathing room instead of a retry storm that keeps it down through the flight-processing window. |
--concurrency (prefork) |
int | cores | Band-math worker processes. Set to physical cores, not hyperthreads — index math is FLOP-bound and gains little from SMT while doubling memory pressure. |
4. Edge Cases & Failure Modes
Duplicate writes from redelivered tasks. With acks_late=True, a worker killed after finishing but before acknowledging causes redelivery. Without the atomic-rename idempotency guard from Step 5, the second run overwrites a valid tile mid-read by a downstream index job, producing a truncated NDVI raster that passes a naive open-and-read check. Always write to a temporary key and rename only after validation.
Mixed CRS across a batch. A cooperative flying fields on either side of a UTM zone boundary will deliver tiles in two EPSG codes. Enqueuing them into one mosaic job produces a seam or a silent reprojection. The manifest carries the EPSG per tile and Step 4 asserts the raster matches; resolve conflicts using the CRS-validation approach in understanding CRS in precision agriculture before the batch runs, not after.
The large-payload anti-pattern. Passing whole NumPy arrays or multi-megabyte GeoTIFF bytes as Celery task arguments serialises them through Redis, bloating the broker and slowing every worker. Pass URIs and small metadata only; let the worker read the pixels from object storage. A task argument should never exceed a few kilobytes.
Retry storms against a throttling backend. When object storage starts returning throttling errors, naive fixed-interval retries hammer it in lockstep and prolong the outage. Exponential backoff with jitter (Step 5) desynchronises the retries so the endpoint recovers. This matters most during the post-flight ingest burst when hundreds of tiles land at once.
Cloud-contaminated tiles poisoning the aggregate. Batch throughput tempts teams to skip masking. A single un-masked cloud pixel shifts NDVI by 0.15–0.40 at that location, and a temporal composite built from unmasked tiles carries that error into the seasonal mean. Run cloud masking for agricultural imagery as an upstream task in the same manifest, never as an afterthought.
Stragglers and back-pressure. One field with ten thousand tiny tiles can flood the CPU queue and starve every other field. Bound the enqueue rate — enqueue field by field, or cap in-flight tasks per field — so no single flight monopolises the workers. Redis queue length is the back-pressure signal to watch.
5. Verification & Output Validation
Completeness is the property that matters: every enqueued tile must reach a terminal state, and every done tile must have a readable, correctly-projected output. Reconcile the manifest against object storage after the run.
import json
import rasterio
def reconcile_run(manifest_path: str) -> dict:
with open(manifest_path) as fh:
manifest = json.load(fh)
counts = {"done": 0, "failed": 0, "pending": 0, "bad_output": 0}
for tile in manifest["tiles"]:
state = tile["state"]
counts[state] = counts.get(state, 0) + 1
if state == "done":
with rasterio.open(tile["dst_uri"]) as src:
# Output must exist, be single-band NDVI, and match manifest CRS
if src.count != 1 or src.crs.to_epsg() != tile["epsg"]:
counts["bad_output"] += 1
total = len(manifest["tiles"])
assert counts["pending"] == 0, (
f"{counts['pending']} tiles never reached a terminal state — run is incomplete"
)
assert counts["bad_output"] == 0, (
f"{counts['bad_output']} tiles marked done but failed output validation"
)
print(f"Reconciled {total} tiles: {counts}")
return counts
report = reconcile_run("run_manifest.json")
Two independent checks confirm a healthy run: the queue drains to zero (no task stuck in RUNNING past the visibility timeout), and the reconciliation above reports zero pending and zero bad outputs. If failed tiles remain, they are isolated in the manifest with their attempt count and error, so a targeted re-run touches only those tiles rather than reprocessing the whole fleet.
6. Integration with the Pipeline
This orchestration layer wraps the single-scene workflows described elsewhere in this section rather than replacing them.
Upstream. Raw captures become per-flight orthomosaics through the orthomosaic stitching workflows — the stitching task on the I/O queue is a batch wrapper around exactly that process, with each flight producing the COG tiles this pipeline then consumes.
Masking, then math. Each tile passes through cloud masking for agricultural imagery before the CPU queue computes indices, so contaminated pixels never enter the ratio. The index computation itself reuses the windowed patterns from band math and raster algebra in Python, lifted into a Dask-parallel, queue-driven task.
The two child guides. The coarse-grained task distribution — a runnable Celery app, chords and groups, retry and result verification — is built step by step in orchestrating multi-flight jobs with Celery and Redis. The fine-grained per-tile parallelism — LocalCluster, dask.array, map_blocks, and bounded-memory compute — is built in parallel tile processing with Dask.
Downstream. The validated index COGs feed straight into temporal aggregation of vegetation indices, where clean per-date tiles are composited across the season into the surfaces that drive prescription generation.
Frequently Asked Questions
Should stitching and band math share the same Celery queue?
No. Stitching is I/O-bound and network-heavy while band math is CPU-bound and memory-heavy, so they scale on different resources. Route them to separate named queues with dedicated worker pools, using a gevent or thread pool for the I/O workers and a prefork pool sized to physical cores for the compute workers. Sharing one queue lets slow uploads starve the CPU workers and inflates end-to-end latency.
How do I stop a retried task from writing a corrupt half-finished orthomosaic?
Make every task idempotent by writing to a temporary path keyed on a deterministic task id, then atomically renaming to the final output only after the write completes and a validation check passes. On retry the task first checks whether the validated output already exists and returns early if so. Never mutate the final artifact in place, because a worker killed mid-write leaves a truncated GeoTIFF that downstream index jobs will read as valid.
Why does my batch run exhaust memory even though each tile is small?
Unbounded concurrency is the usual cause: a prefetch multiplier that pulls hundreds of tile tasks into one worker, or a Dask graph that materialises every block before writing. Set the Celery prefetch multiplier to 1 for long compute tasks, cap the Dask cluster memory limit per worker, and process tiles in windows rather than loading the full orthomosaic. Back-pressure on the queue is what keeps peak RAM flat regardless of how many flights are enqueued.
Related
- Orchestrating Multi-Flight Jobs with Celery and Redis — a runnable Celery app with fan-out, chords, exponential backoff, and result verification
- Parallel Tile Processing with Dask — computing NDVI across the tiles of a large orthomosaic under bounded memory
- Orthomosaic Stitching Workflows — the per-flight stitching the I/O queue wraps at fleet scale
- Cloud Masking for Agricultural Imagery — the upstream masking task every tile passes through before index math
- Band Math & Raster Algebra in Python — the windowed NDVI, NDRE, and SAVI patterns each compute task runs