Parametrizing Spatial DAGs by Tile Index

To fan a spatial DAG across a tile grid, derive the tile list from a raster’s extent and CRS, then map one idempotent task per tile so each z/x/y cell runs, retries, and materializes independently. Parametrizing spatial DAGs by tile index turns a single monolithic raster job into hundreds of small, parallel units of work, where a failure on one tile never forces a rerun of the whole scene. This guide shows how to compute the tile set from bounds, express the fan-out with Prefect .map() or Dagster dynamic outputs, and keep every per-tile task safe to replay. It builds directly on the DAG Design Principles for Spatial ETL section, sharpening the partitioning principle into a concrete, indexable parameter.

When to Use This Pattern#

Reach for tile-index parametrization when your workload has these traits:

• The input raster or mosaic is larger than a single worker’s memory, so it must be read in windows rather than loaded whole.
• The transform is embarrassingly parallel per tile — resampling, reprojection, NDVI, hillshade, or cloud masking that needs no cross-tile state.
• You want granular retries: one corrupt or timed-out tile should re-execute in isolation, not restart the run.
• Your output is itself a tiled product (XYZ tiles, a Cloud-Optimized GeoTIFF built from overviews, or per-tile Parquet) where each cell maps to a stable, addressable path.
• You need to backfill or reprocess a subset of the grid — a single UTM zone, a bounding box over one country — without touching neighboring tiles.

If tiles must share running state (a global histogram, a seam-aware blend), pure fan-out is the wrong tool; use it for the independent stage and reserve the reduction for a dedicated aggregation task downstream.

Choosing the zoom level is the one decision that governs everything else. Tile count roughly quadruples per level, so a raster that fits in nine tiles at zoom 8 explodes to well over a hundred at zoom 12. Match the zoom to the native ground sample distance of the source: rendering a 30-meter Landsat scene at zoom 15 (roughly 4.7 m/px at the equator) wastes compute upsampling data that carries no extra detail, while rendering sub-meter imagery at zoom 10 throws away resolution you paid to acquire. A useful rule is to pick the coarsest zoom whose pixel size is still finer than the source, then let overviews handle everything below it.

Complete Working Example#

The example below derives an XYZ tile list from a raster’s geographic bounds, then maps a per-tile reprojection-and-write task. The tile enumeration is a plain function so it is unit-testable in isolation; the per-tile task takes the full tuple it needs and writes to a deterministic path.

import math
from dataclasses import dataclass
from pathlib import Path

import mercantile
import rasterio
from pyproj import CRS
from rasterio.warp import transform_bounds


@dataclass(frozen=True)
class TileJob:
    z: int
    x: int
    y: int
    src_uri: str


def tiles_for_raster(src_uri: str, zoom: int) -> list[TileJob]:
    """Enumerate Web Mercator XYZ tiles covering a raster's extent."""
    with rasterio.open(src_uri) as src:
        src_crs: CRS = src.crs
        if src_crs is None:
            raise ValueError(f"{src_uri} has no CRS; cannot derive a tile index.")
        # Reproject the dataset bounds to EPSG:4326 for mercantile.
        west, south, east, north = transform_bounds(
            src_crs, CRS.from_epsg(4326), *src.bounds, densify_pts=21
        )
    tiles = mercantile.tiles(west, south, east, north, zooms=zoom)
    return [TileJob(z=t.z, x=t.x, y=t.y, src_uri=src_uri) for t in tiles]


def render_tile(job: TileJob, out_root: str, tile_px: int = 256) -> str:
    """Reproject one tile window to EPSG:3857 and write it idempotently."""
    from rasterio.warp import calculate_default_transform, reproject, Resampling

    out_path = Path(out_root) / str(job.z) / str(job.x) / f"{job.y}.tif"
    if out_path.exists():
        return str(out_path)  # already materialized; safe to skip

    bounds = mercantile.xy_bounds(job.x, job.y, job.z)  # meters, EPSG:3857
    dst_crs = CRS.from_epsg(3857)
    tmp_path = out_path.with_suffix(".tif.part")
    out_path.parent.mkdir(parents=True, exist_ok=True)

    with rasterio.open(job.src_uri) as src:
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, tile_px, tile_px,
            left=bounds.left, bottom=bounds.bottom,
            right=bounds.right, top=bounds.top,
        )
        profile = src.profile.copy()
        profile.update(
            crs=dst_crs, transform=transform, width=width, height=height,
            driver="GTiff", nodata=src.nodata if src.nodata is not None else 0,
        )
        with rasterio.open(tmp_path, "w", **profile) as dst:
            for band in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, band),
                    destination=rasterio.band(dst, band),
                    src_crs=src.crs, dst_crs=dst_crs,
                    dst_transform=transform, resampling=Resampling.bilinear,
                    src_nodata=src.nodata, dst_nodata=profile["nodata"],
                )
    tmp_path.replace(out_path)  # atomic rename: readers never see a partial tile
    return str(out_path)

Wire the fan-out with Prefect .map(), which submits one task run per element and returns a list of futures:

from prefect import flow, task
from prefect.logging import get_run_logger


@task(retries=3, retry_delay_seconds=[10, 30, 90], persist_result=True)
def render_tile_task(job: TileJob, out_root: str) -> str:
    return render_tile(job, out_root)


@flow(name="tile-index-mosaic")
def mosaic_flow(src_uri: str, out_root: str, zoom: int = 10) -> int:
    logger = get_run_logger()
    jobs = tiles_for_raster(src_uri, zoom)
    logger.info("Derived %d tiles at zoom %d", len(jobs), zoom)
    futures = render_tile_task.map(job=jobs, out_root=out_root)
    written = [f.result() for f in futures]
    return len(written)

The Dagster equivalent expresses the same fan-out with dynamic outputs. DynamicOut yields one DynamicOutput per tile, and .map() on the resulting handle runs the per-tile op in parallel across the executor:

from dagster import op, job, DynamicOut, DynamicOutput, OpExecutionContext, RetryPolicy


@op(out=DynamicOut(TileJob))
def enumerate_tiles(context: OpExecutionContext, src_uri: str, zoom: int):
    for t in tiles_for_raster(src_uri, zoom):
        # mapping_key must be a stable, unique string per tile.
        yield DynamicOutput(t, mapping_key=f"z{t.z}_x{t.x}_y{t.y}")


@op(retry_policy=RetryPolicy(max_retries=3, delay=10))
def render_tile_op(job: TileJob, out_root: str) -> str:
    return render_tile(job, out_root)


@job
def tile_mosaic_job():
    enumerate_tiles().map(lambda j: render_tile_op(j))

Both versions share the same tiles_for_raster and render_tile core, so the orchestration choice stays orthogonal to the spatial logic — the split emphasized throughout Prefect vs Dagster for GIS Workloads.

The fan-out mechanics differ in one operationally important way. Prefect’s .map() submits every tile as its own task run and returns a list of futures immediately; the task runner (thread, process, or Dask/Ray) decides how many execute concurrently, and you cap that with a task-runner concurrency limit or a work-pool queue so a thousand tiles do not open a thousand simultaneous GDAL handles. Dagster’s dynamic .map() behaves similarly but is bounded by the executor configuration on the job — the multiprocess executor’s max_concurrent setting, or per-op concurrency tags. In both engines, unbounded fan-out is the classic way to exhaust file descriptors or saturate S3 request quotas, so treat the concurrency cap as a required parameter, not an afterthought. Because each tile is a discrete run, the orchestrator’s UI also gives you per-tile retry visibility: you can re-run the seventeen tiles that hit a transient read timeout without disturbing the thousands that already succeeded.

Parameter & Option Reference#

Tile enumeration follows the OGC/Web Mercator XYZ convention documented by mercantile, which keeps the index compatible with slippy-map tile servers and Cloud-Optimized GeoTIFF overviews.

Verification & Testing#

Test the enumeration and the per-tile output separately. The tile list is deterministic, so assert its size and membership against a known extent:

def test_tiles_cover_and_are_unique():
    jobs = tiles_for_raster("tests/data/scene_utm.tif", zoom=8)
    keys = {(j.z, j.x, j.y) for j in jobs}
    assert len(keys) == len(jobs)          # no duplicate tile indices
    assert all(j.z == 8 for j in jobs)     # every tile at requested zoom
    assert len(jobs) > 0                   # extent actually intersects the grid

Confirm that mapping twice does not corrupt or duplicate outputs — the out_path.exists() guard plus the atomic .part rename makes render_tile safe to replay, which is the core of tile-level idempotency covered in generating idempotent keys for shapefile uploads. Verify a written tile with the GDAL CLI:

gdalinfo build/10/163/395.tif | grep -E "Size is|EPSG|NoData"
# Size is 256, 256
# ... ID["EPSG",3857]]
# NoData Value=0

A quick end-to-end assertion checks that the tile count returned by the flow matches the enumeration, catching silently dropped map elements:

def test_flow_writes_every_tile(tmp_path):
    src = "tests/data/scene_utm.tif"
    expected = len(tiles_for_raster(src, zoom=8))
    written = mosaic_flow(src, str(tmp_path), zoom=8)
    assert written == expected

Common Pitfalls#

• Enumerating tiles in the source CRS. mercantile.tiles() expects geographic degrees. Passing raw UTM meters yields absurd tile indices. Always reproject the dataset bounds to EPSG:4326 first with transform_bounds and densify_pts so curved edges are not clipped.
• Non-deterministic mapping keys. In Dagster, a mapping_key built from a timestamp or uuid4() breaks retries and backfills because the key must be stable across runs. Encode the tile index itself, exactly as the example does.
• Writing directly to the final path. A worker killed mid-write leaves a truncated GeoTIFF that later reads as valid-looking garbage. Write to a .part file and atomically rename — the same durability contract used when checkpointing large raster mosaics.
• Ignoring edge tiles with partial coverage. Tiles straddling the raster boundary are mostly nodata. Dropping nodata from the destination profile turns those margins into opaque black borders in the final mosaic.

Related: Ground the fan-out in DAG Design Principles for Spatial ETL, pair it with checkpointing large raster mosaics for resumable builds, offload the heavy per-tile compute with running async geopandas tasks safely, and avoid recomputing unchanged cells using caching reprojected rasters with content hashing.

← Back to DAG Design Principles for Spatial ETL