This guide is part of Raster Mosaicking & Tiling for Pipeline Outputs, within the broader Automated Vector & Raster Cleaning Workflows reference.

Generating XYZ Tiles from a Raster

Web map tiles are a different product from an analytical raster: a fixed quadtree in Web Mercator, 256 pixels square, addressed by zoom, column and row. Producing them from a pipeline output is mostly a question of deciding how many to produce and when.

Why Tile Generation Goes Wrong

• Rendering every zoom level. Zoom 18 over a country is billions of tiles, almost none of which will ever be requested.
• Rendering empty tiles. Without a footprint intersection, most tiles in a bounding box contain no data and cost storage anyway.
• Implicit contrast stretching. A per-tile stretch makes neighbouring tiles disagree, producing visible checkerboarding in the map.
• Ignoring the projection change. Data in a national grid must be reprojected to Web Mercator, and doing that per tile without a common grid reintroduces the seam problem.

Version and Environment Compatibility

pip install "rio-tiler>=6.4" "morecantile>=5.0" "rasterio>=1.3"

Choosing the Zoom Range

The tiles_for_raster Recipe

from __future__ import annotations

import logging
import math
from pathlib import Path

import morecantile
import rasterio
from rasterio.warp import transform_bounds
from rio_tiler.io import Reader

logger = logging.getLogger(__name__)
TMS = morecantile.tms.get("WebMercatorQuad")


def max_useful_zoom(resolution_m: float, latitude: float = 0.0) -> int:
    """Highest zoom whose ground resolution is still finer than the source."""
    equator_circumference = 40_075_016.686
    scale = math.cos(math.radians(latitude))
    for zoom in range(0, 24):
        tile_resolution = equator_circumference * scale / (256 * 2 ** zoom)
        if tile_resolution <= resolution_m:
            return zoom
    return 23


def tiles_for_raster(cog_path: Path, min_zoom: int | None = None) -> list[morecantile.Tile]:
    """Enumerate only the tiles the raster's footprint actually covers."""
    with rasterio.open(cog_path) as src:
        resolution = abs(src.transform.a)
        west, south, east, north = transform_bounds(src.crs, "EPSG:4326", *src.bounds, densify_pts=21)

    centre_latitude = (south + north) / 2
    max_zoom = max_useful_zoom(resolution, centre_latitude)
    start = min_zoom if min_zoom is not None else max(0, max_zoom - 6)

    tiles = [
        tile
        for zoom in range(start, max_zoom + 1)
        for tile in TMS.tiles(west, south, east, north, [zoom])
    ]
    logger.info("%s: zooms %d–%d, %d tiles", cog_path.name, start, max_zoom, len(tiles))
    return tiles


def render_tile(cog_path: Path, tile: morecantile.Tile, out_root: Path,
                rescale: tuple[float, float], colormap: dict) -> Path | None:
    """Render one tile, skipping any whose window contains no data at all."""
    with Reader(str(cog_path)) as reader:
        image = reader.tile(tile.x, tile.y, tile.z, tilesize=256)

    if image.mask.max() == 0:
        return None   # entirely nodata — do not write an empty object

    image.rescale(in_range=(rescale,))
    content = image.render(img_format="PNG", colormap=colormap)

    destination = out_root / str(tile.z) / str(tile.x) / f"{tile.y}.png"
    destination.parent.mkdir(parents=True, exist_ok=True)
    destination.write_bytes(content)
    return destination

Key Implementation Notes

• The maximum zoom comes from the data, not from the map. Rendering past the source resolution multiplies storage by four per level while adding nothing a client could not interpolate itself.
• The rescale range is passed in, not computed per tile. A per-tile stretch is the cause of visible checkerboarding, because each tile then maps a different value range onto the same colours.
• Empty tiles are skipped, not written. A transparent PNG costs storage, a request and a cache entry; a 404 costs none of those and most clients handle it correctly.
• The colormap is an argument. It belongs in versioned configuration next to the product, because a changed ramp changes every published tile and needs the same version bump as changed data.
• Latitude enters the zoom calculation. Web Mercator resolution scales with the cosine of latitude, so a Scandinavian product reaches its natural zoom a level earlier than an equatorial one.
• Tiles are written under a zoom/x/y path. That is the convention every client expects, and it makes a whole zoom level deletable as a prefix.

Cache Headers and Versioned Paths

A tile that is served without a cache policy is fetched again on every pan, which wastes bandwidth on both sides and makes the map feel slow. A tile that is cached forever cannot be corrected. The way out is to make the path change when the content changes.

Put a version segment in the tile path — tiles/ndvi/v7/{z}/{x}/{y}.png — and set a long immutable cache lifetime on the objects. A new render publishes under v8, the client is pointed at the new prefix, and nothing needs to be invalidated because nothing was overwritten. The old version stays available until it is deliberately expired, which also makes a rollback a one-line change.

The alternative, overwriting tiles in place and relying on cache invalidation, works until a CDN somewhere holds a stale copy for a week and nobody can explain why one region of the map disagrees with the data.

Troubleshooting Tile Rendering

Integration Note

Tile generation is a publication concern, not a processing one: it consumes the validated COGs from writing cloud-optimized GeoTIFFs with rio-cogeo and produces a presentation artefact. Run it as a separate downstream job so that a slow render cannot delay the data pipeline, and key it on the same tile identifiers so that only the map tiles overlapping a rebuilt data tile are regenerated.

Serving Values, Not Only Colours

A tiled map answers “what does this look like”; analysts also ask “what is the value here”. Three additions cover that without a second product.

A point-query endpoint reads the single pixel under a coordinate straight from the COG. It is the same windowed read the tiler already performs, bounded to one pixel, and it returns the actual value rather than a colour.

Numeric tiles — sometimes called terrain-RGB — encode values into the colour channels so a client can decode them. Useful for elevation and continuous indices where a browser needs the numbers, at the cost of an encoding scheme every consumer must implement identically.

A legend document published alongside the tiles states the ramp, its breaks and the units. Without it a colour map is uninterpretable, and the legend is invariably the thing that goes stale first — which is another reason to version it in the same path segment as the tiles it describes.

Related

• Raster Mosaicking & Tiling for Pipeline Outputs
• Writing Cloud-Optimized GeoTIFFs with rio-cogeo
• Building Virtual Raster Mosaics with GDAL VRT
• Reprojecting Rasters to a Common Grid with rasterio.warp