Pre-Seeding a WMTS Tile Cache with gdal2tiles

TL;DR: Warp your source raster to EPSG:3857, then run gdal2tiles with --profile mercator to pre-render a {z}/{x}/{y}.png pyramid whose grid is identical to the GoogleMapsCompatible WMTS tile matrix set. Drive it from Python with subprocess or the osgeo_utils.gdal2tiles module, use --processes to parallelise across cores, and choose --xyz or the default TMS row ordering to match how your tile server addresses TileRow.

The Core Challenge

A Web Map Tile Service answers requests by looking up a pre-computed tile at a fixed (TileMatrix, TileRow, TileCol) address — it does no rendering at request time. That is the whole point of the tiled model covered in WMTS Tile Matrix Sets Explained: the expensive rasterisation happens once, ahead of time, and every subsequent request is a cheap file read. But that only works if the tiles you generate land on exactly the grid the WMTS TileMatrixSet describes.

The trap is coordinate alignment. A WMTS TileMatrixSet pins down four things for every zoom level: the coordinate reference system, the TopLeftCorner origin, the tile pixel size (almost always 256×256), and a ScaleDenominator per TileMatrix. The near-universal GoogleMapsCompatible set defines these over EPSG:3857 (Web Mercator) with the origin at the top-left of the world extent. If gdal2tiles emits tiles on a differently-aligned grid — a different projection, a different origin, or the opposite row numbering — a WMTS client will request address 12/2048/1362 and get a tile that covers the wrong ground. The map does not error; it simply shows the wrong place, or shows it upside down.

gdal2tiles produces exactly the GoogleMapsCompatible grid when you invoke it with --profile mercator. The remaining decisions — zoom range, resampling, parallelism, and tile-row ordering — determine how fast the seed runs and whether the output registers with your service.

Production-Ready Code

The script below warps a source raster to Web Mercator, then seeds a tile pyramid. It shows both entry points: the osgeo_utils.gdal2tiles module API (preferred — no shelling out, and it raises Python exceptions on failure) and a subprocess fallback for environments where only the gdal2tiles.py command-line utility is on the PATH. The only dependency is GDAL 3.4+ with its Python bindings (pip install gdal==$(gdal-config --version)).

"""
seed_wmts_cache.py — pre-render a WMTS/XYZ tile pyramid from a source raster.

Dependencies:
    GDAL >= 3.4 with Python bindings (osgeo, osgeo_utils)

The output directory follows the {z}/{x}/{y}.png convention and, with
--profile mercator, aligns tile-for-tile with the GoogleMapsCompatible
WMTS TileMatrixSet (EPSG:3857, 256 px tiles, top-left world origin).
"""
from __future__ import annotations

import os
import shutil
import subprocess
from pathlib import Path

from osgeo import gdal

gdal.UseExceptions()

WEB_MERCATOR = "EPSG:3857"


def warp_to_web_mercator(src: Path, dst: Path, resampling: str = "bilinear") -> Path:
    """
    Reproject the source raster to EPSG:3857 so tile boundaries align cleanly
    to the Web Mercator grid. Pre-warping (rather than letting gdal2tiles warp
    on the fly at every zoom) gives explicit control over resampling and nodata.
    """
    gdal.Warp(
        str(dst),
        str(src),
        dstSRS=WEB_MERCATOR,
        resampleAlg=resampling,          # bilinear/cubic for imagery, near for categorical
        dstAlpha=True,                   # write an alpha band so edges are transparent
        multithread=True,
        warpOptions=["NUM_THREADS=ALL_CPUS"],
    )
    return dst


def seed_via_module(
    src: Path,
    out_dir: Path,
    zoom: str = "0-14",
    processes: int = os.cpu_count() or 4,
    xyz: bool = True,
    resampling: str = "average",
) -> Path:
    """
    Seed the pyramid using the osgeo_utils.gdal2tiles module API.

    --profile mercator  -> GoogleMapsCompatible grid (EPSG:3857)
    --xyz               -> row 0 at the north edge (matches WMTS TileRow)
    --processes         -> fan tile generation across CPU cores
    --resampling average-> good default for down-sampled overview levels
    """
    from osgeo_utils import gdal2tiles

    argv = [
        "gdal2tiles.py",
        "--profile", "mercator",
        "--zoom", zoom,
        "--processes", str(processes),
        "--resampling", resampling,
        "--webviewer", "none",           # skip the HTML/leaflet viewer output
        "--tilesize", "256",
    ]
    if xyz:
        argv.append("--xyz")             # omit for TMS (bottom-origin) row order
    argv += [str(src), str(out_dir)]

    gdal2tiles.main(argv)
    return out_dir


def seed_via_subprocess(
    src: Path,
    out_dir: Path,
    zoom: str = "0-14",
    processes: int = os.cpu_count() or 4,
    xyz: bool = True,
    resampling: str = "average",
) -> Path:
    """Fallback: invoke the gdal2tiles.py CLI when only the binary is available."""
    exe = shutil.which("gdal2tiles.py") or shutil.which("gdal2tiles")
    if exe is None:
        raise FileNotFoundError("gdal2tiles.py not found on PATH")

    cmd = [
        exe,
        "--profile", "mercator",
        "--zoom", zoom,
        "--processes", str(processes),
        "--resampling", resampling,
        "--webviewer", "none",
        "--tilesize", "256",
    ]
    if xyz:
        cmd.append("--xyz")
    cmd += [str(src), str(out_dir)]

    # check=True raises CalledProcessError on a non-zero exit code.
    subprocess.run(cmd, check=True, text=True, capture_output=True)
    return out_dir


if __name__ == "__main__":
    source = Path("input/orthophoto.tif")
    warped = Path("build/orthophoto_3857.tif")
    tiles = Path("build/tiles")

    warp_to_web_mercator(source, warped)
    seed_via_module(warped, tiles, zoom="0-16", xyz=True)
    print(f"Pyramid written to {tiles}/{{z}}/{{x}}/{{y}}.png")

Step-by-Step Walkthrough

Pre-warp to Web Mercator (warp_to_web_mercator). gdal.Warp reprojects the source to EPSG:3857 before any tiling happens. Doing this once, up front, is better than letting gdal2tiles warp on the fly: you pick the resampleAlg explicitly (bilinear or cubic for continuous imagery, near for categorical land-cover so class codes are not blended), and dstAlpha=True adds an alpha band so tiles outside the data footprint come out transparent rather than black. Because the WMTS GoogleMapsCompatible matrix set is defined in EPSG:3857, warping to that CRS is what makes the eventual tile boundaries fall on the matrix-set grid lines. Reprojection between EPSG:4326 and EPSG:3857 has its own axis-order subtleties covered in SRS and Coordinate Reference System Handling.

Choosing the profile (--profile mercator). This is the single most important flag. mercator selects the EPSG:3857 grid with a 256-pixel tile size and the world origin at the top-left corner — byte-for-byte the geometry of the GoogleMapsCompatible tile matrix set. The alternative --profile geodetic builds an EPSG:4326 grid that lines up with the WorldCRS84Quad matrix set instead, and --profile raster builds a non-geographic pyramid that no WMTS client can address. Match the profile to the TileMatrixSet your service advertises; for aligning custom or global grids, see Configuring Tile Matrix Sets for Global WMTS Deployments.

The zoom range (--zoom). --zoom 0-16 generates every TileMatrix level from 0 (one tile for the whole world) through 16. Each level down quadruples the tile count, so the deepest level dominates both render time and disk usage. Seed only the levels your service will serve directly; deeper levels are usually cheaper to render lazily on first request.

Parallelism (--processes). gdal2tiles splits the base-zoom tile generation across a process pool. Setting --processes to os.cpu_count() typically gives near-linear speedup on the base level because each tile is independent. The overview levels (built by down-sampling the level below) parallelise less well, but the base level is where the bulk of the work is.

Tile-row ordering (--xyz). With --xyz, row 0 is the northernmost tile — the same numbering WMTS uses for TileRow and the same as the XYZ/Google/OSM convention. Omit the flag and gdal2tiles writes the classic TMS layout, where row 0 is the southernmost tile. The two differ by y_tms = 2^z - 1 - y_xyz. Pick the one your tile server maps to TileRow without flipping; getting it wrong renders every map vertically mirrored.

Resampling (--resampling). This controls how overview (lower-zoom) tiles are built from the base level. average is a solid default for imagery; bilinear/cubic are smoother; mode or near preserves discrete class values for categorical rasters. It is independent of the resampleAlg used during the initial warp.

Verification

Confirm the pyramid exists and that a known tile coordinate resolves to the right ground. The tile count per zoom level should match the source extent, and the top-level tile should always be 0/0/0.png.

# 1. Count tiles per zoom level (directories are the z values)
for z in build/tiles/[0-9]*; do
  printf "zoom %s: %s tiles\n" "$(basename "$z")" "$(find "$z" -name '*.png' | wc -l)"
done

# 2. Inspect one tile's georeferencing — gdalinfo reads the embedded/implied extent
gdalinfo build/tiles/0/0/0.png

Expected output shape:

zoom 0: 1 tiles
zoom 1: 4 tiles
zoom 2: 16 tiles
...
zoom 16: 24187 tiles

Driver: PNG/Portable Network Graphics
Size is 256, 256

A quick Python check confirms an XYZ tile address maps back to the expected Web Mercator extent, so you can spot an off-by-one or a flipped row before wiring the cache into a service:

import math


def xyz_tile_bounds_3857(z: int, x: int, y: int) -> tuple[float, float, float, float]:
    """Return (minx, miny, maxx, maxy) in EPSG:3857 for an XYZ tile (y from north)."""
    world = 20037508.342789244            # half the Web Mercator world span (metres)
    span = (2 * world) / (2 ** z)         # tile side length at this zoom
    minx = -world + x * span
    maxy = world - y * span               # y counts down from the north edge
    return (minx, maxy - span, minx + span, maxy)


print(xyz_tile_bounds_3857(2, 2, 1))
# (0.0, 5009377.085697312, 5009377.085697311, 10018754.171394624)

Gotchas & Edge Cases

Back to WMTS Tile Matrix Sets Explained

Related

• WMTS Tile Matrix Sets Explained — how TileMatrixSet, ScaleDenominator, and TopLeftCorner define the grid your tiles must land on
• Configuring Tile Matrix Sets for Global WMTS Deployments — aligning custom and global grids beyond the default GoogleMapsCompatible set
• Understanding OGC Web Map Service Specifications — when to render dynamically with WMS instead of pre-seeding a static tile cache