Building XYZ Tile Pyramids with gdal2tiles

gdal2tiles.py turns a raster into a directory of web tiles in a single invocation, and the defaults are wrong for a survey orthomosaic in several ways at once: the wrong tile scheme, a zoom range guessed rather than computed, no alpha handling, and a single process on a job that will take hours. Each is a flag, and each one silently produces output that looks finished.

This page covers the six decisions the command actually encodes, and how to verify the run completed rather than merely stopped.

What the defaults get wrong

The tile scheme. Without --xyz, tiles are written in the TMS convention, whose row index counts from the south. Every tile is correct and every tile is in the wrong place for a web map client, which renders the map mirrored about the equator.

The zoom range. --zoom defaults to a range derived from the raster, which is usually one or two levels beyond what the imagery supports. Each surplus level quadruples the tile count while upsampling.

Transparency. An orthomosaic has an irregular footprint. Without an alpha band in the source, the area outside it is written as opaque black, which sits over whatever basemap the client has beneath.

Parallelism. --processes defaults to one. On a survey-sized raster that is the difference between forty minutes and six hours.

Four defaults and what each costs Four rows comparing a default setting against the correct one for a survey orthomosaic. The tile scheme defaults to TMS, whose row index counts from the south, producing a map mirrored about the equator; the correct setting is XYZ. The zoom range defaults beyond the source resolution, quadrupling the tile count per surplus level; the correct range is computed from the ground sample distance. Transparency defaults to none, writing opaque black outside the footprint; the correct approach adds a destination alpha band during the warp. Parallelism defaults to a single process, taking hours instead of minutes. default correct for a survey tile scheme zoom range transparency parallelism TMS map is mirrored guessed from the raster ×4 tiles per surplus level none black outside the footprint 1 process hours instead of minutes --xyz computed from the GSD -dstalpha at the warp --processes = cores

Figure 1 — Four flags separate a working tile set from a plausible-looking one. Three of the four defaults produce output that renders without error.

Minimal reproducible solution

Compute the zoom range, warp with alpha and pixel alignment, then generate — as one function so the three steps cannot drift apart.

import math
import subprocess
from pathlib import Path

import rasterio
from rasterio.warp import transform_bounds


def build_pyramid(src: Path, out_dir: Path, processes: int = 8,
                  tile_format: str = "PNG") -> tuple[int, int]:
    """Warp to Web Mercator and generate an XYZ pyramid. Returns (min, max) zoom."""
    warped = out_dir.parent / f"{src.stem}_3857.tif"

    subprocess.run(
        ["gdalwarp", "-t_srs", "EPSG:3857", "-r", "cubic", "-tap",
         "-dstalpha", "-multi", "-wo", "NUM_THREADS=ALL_CPUS",
         "-co", "TILED=YES", "-co", "COMPRESS=DEFLATE",
         "-overwrite", str(src), str(warped)],
        check=True,
    )

    with rasterio.open(warped) as ds:
        west, south, east, north = transform_bounds(ds.crs, "EPSG:4326", *ds.bounds)
        gsd_m = abs(ds.transform.a)          # already metres in 3857

    lat = math.radians((north + south) / 2)
    max_zoom = max(0, math.ceil(math.log2(156_543.03 * math.cos(lat) / gsd_m)))
    span_deg = max(east - west, north - south)
    min_zoom = max(0, int(math.floor(math.log2(360.0 / max(span_deg, 1e-9)))))

    subprocess.run(
        ["gdal2tiles.py", "--xyz", f"--zoom={min_zoom}-{max_zoom}",
         f"--processes={processes}", "--resampling=cubic",
         f"--tiledriver={tile_format}", "--webviewer=none",
         "--exclude",                       # skip fully transparent tiles
         str(warped), str(out_dir)],
        check=True,
    )
    return min_zoom, max_zoom

--exclude is worth knowing about: it skips tiles that are entirely transparent, which on a survey with an irregular footprint removes a substantial fraction of the files and every one of them would have been a 300-byte blank PNG. The client’s map renders nothing there either way, and the deployment is meaningfully smaller.

Deriving max_zoom from the ground sample distance including the cosine-of-latitude term is what stops the run generating levels that upsample. At 60° north that term alone saves a full zoom level, which is four times the tiles of everything below it combined.

Edge-case matrix

Situation Symptom Handling
--xyz omitted Map mirrored vertically Regenerate; renaming does not work
No alpha in the source Black outside the footprint Warp with -dstalpha
JPEG tiles with an alpha band Alpha silently dropped Use PNG or WebP, or ship a mask
Zoom beyond the source GSD Huge tile count, no detail Compute the range
Single process Hours of wall clock Set --processes
Source not in Web Mercator gdal2tiles warps implicitly Warp explicitly, with -tap
Very large raster Memory pressure per process Fewer processes, or tile the source
Run interrupted Partial pyramid, looks complete Verify coverage, do not assume

The interrupted-run row is the one that reaches production. A pyramid that stopped at zoom 19 of a 14–21 range contains hundreds of thousands of correct tiles and is missing the two levels the client will actually use at close range; nothing about the directory says so.

Verification snippet

Check coverage per zoom level against what the raster’s bounds imply, rather than checking that the command exited zero.

import math
from pathlib import Path

import rasterio
from rasterio.warp import transform_bounds


def _tile_range(west, south, east, north, z):
    n = 2 ** z
    def xy(lon, lat):
        x = int((lon + 180.0) / 360.0 * n)
        y = int((1.0 - math.asinh(math.tan(math.radians(lat))) / math.pi) / 2.0 * n)
        return x, y
    x0, y0 = xy(west, north)
    x1, y1 = xy(east, south)
    return x0, y0, x1, y1


def assert_pyramid_complete(tile_dir: Path, warped_src: Path,
                            min_zoom: int, max_zoom: int,
                            min_present_frac: float = 0.6) -> None:
    """Every zoom level must exist and be substantially populated."""
    with rasterio.open(warped_src) as ds:
        west, south, east, north = transform_bounds(ds.crs, "EPSG:4326", *ds.bounds)

    for z in range(min_zoom, max_zoom + 1):
        level = tile_dir / str(z)
        assert level.is_dir(), f"zoom {z} missing entirely — the run stopped early"

        x0, y0, x1, y1 = _tile_range(west, south, east, north, z)
        expected = (x1 - x0 + 1) * (y1 - y0 + 1)
        present = sum(1 for _ in level.glob("*/*"))
        # --exclude removes transparent tiles, so require a fraction not all.
        frac = present / max(expected, 1)
        assert frac >= min_present_frac, (
            f"zoom {z}: {present} of ~{expected} tiles present ({frac:.0%}) — "
            "either the run was interrupted or the footprint is unexpectedly sparse")

The fractional threshold rather than an exact count is what makes this compatible with --exclude: a survey whose footprint is a diagonal strip legitimately fills about half its bounding box, and requiring every tile would fail on correct output. Sixty percent is a reasonable floor for a survey block and should be lowered deliberately for genuinely thin corridors.

Tile count per zoom level, complete and interrupted A bar chart of tiles present per zoom level from fourteen to twenty-one, for two runs. In the complete run, the count quadruples at each level and the top level holds the overwhelming majority of the files. In the interrupted run, the counts are identical up to zoom nineteen and then stop: zoom twenty and twenty-one have no tiles at all. A note observes that the interrupted directory contains over ninety percent of the correct tile names by count of levels but under a tenth of the files, and that nothing in the directory signals the omission. 14 15 16 17 18 19 20 21 zoom level complete interrupted run stops here Six of eight levels are present and under a tenth of the files; the directory looks populated and the client sees nothing at close zoom.

Figure 2 — Why per-level verification matters. Tile counts are so lopsided toward the top levels that an interrupted run still produces a convincing directory.

Making the pyramid self-describing

A tile directory carries no metadata, so the last step of the generation is to write the three small files that turn it into something a client can open without instructions.

A TileJSON document at the root states the URL template, the zoom range, the bounds and a sensible centre. Every common web-map library reads it directly, which removes the round of correspondence that otherwise begins with “what zoom levels are in this?”.

A provenance record — the source raster’s name and checksum, the run identifier, the generation date, the resampling and the tile format — is what makes two revisions distinguishable six months later. Without it, two directories differing by one reprocessing are indistinguishable except by comparing pixels.

A minimal viewer page, twenty lines against a mapping library, is what makes the deliverable openable by a non-technical recipient. It is the single highest-leverage file in the set for reducing support traffic, and it is generated rather than written.

All three values are already known at the moment the pyramid is generated, so writing them costs a few lines at the end of the function that produced the tiles.

When to escalate

  • Generation is I/O-bound rather than CPU-bound. Past about eight processes the bottleneck is writing hundreds of thousands of small files. A filesystem tuned for large files handles this badly; writing to a local disk and syncing afterwards is usually faster than generating directly onto network storage.
  • The pyramid is correct and the client still sees blank areas. The tiles exist and the server is not serving them, most often because the deployment stopped at a directory-count limit. Compare the file count on disk against the count at the destination.
  • The raster is large enough that generation is impractical. This is the case the on-demand serving path exists for: a COG plus a tile service has no generation step at all, and for a deliverable still under review it is usually the better answer.

← Tiling and Serving Orthomosaics as XYZ Tiles

Tile format against total pyramid size and transparency support Three tile formats compared for a survey orthomosaic pyramid. PNG is lossless and supports an alpha channel, producing the largest pyramid at roughly four point two gigabytes. JPEG is lossy and has no alpha, producing the smallest at about one point one gigabytes but rendering the area outside the footprint as opaque black. WebP is lossy and supports alpha, at about zero point six gigabytes, and is the smallest option that still handles the irregular footprint. A note observes that only PNG is appropriate where pixel values carry meaning rather than appearance. PNG JPEG WebP 4.2 GB 1.1 GB 0.6 GB no alpha — black outside the footprint alpha supported alpha supported, lossless For a visual orthomosaic, WebP is the smallest option that still handles the footprint. For any raster whose pixel values carry meaning, PNG is the only correct choice regardless of size.

Figure 3 — The format decision in the terms that matter: size, transparency, and whether the pixels are a picture or a measurement.