Tiling and Serving Orthomosaics as XYZ Tiles

A finished orthomosaic is routinely 20–60 GB. A client wants to look at it in a browser, on a laptop, over a hotel connection. The bridge between those two facts is the XYZ tile scheme: a fixed pyramid of 256-pixel PNG or JPEG tiles addressed by zoom, column and row, which every web map library speaks natively.

This page covers what that grid actually is, how to decide which zoom levels are worth producing, the reprojection into Web Mercator that every tile set requires, and the choice between pre-rendering a pyramid of files and serving tiles on demand from a Cloud-Optimized GeoTIFF. The second option is right more often than its lower profile suggests.

The source raster is the one produced in exporting Cloud-Optimized GeoTIFFs with rasterio, and the reprojection concerns are those handled in managing coordinate reference systems in GDAL.

Audience and prerequisites. Python 3.10+, GDAL ≥ 3.6, and an orthomosaic with a declared CRS.

Prerequisites

Library / tool Minimum version Install command Role
GDAL ≥ 3.6 conda install gdal gdal2tiles.py, gdalwarp
rasterio ≥ 1.3 pip install "rasterio>=1.3" Reading the source, checking the result
morecantile ≥ 5.0 pip install "morecantile>=5.0" Tile-grid arithmetic without hand-rolled maths
rio-tiler ≥ 6.4 pip install "rio-tiler>=6.4" On-demand tiles from a COG

Conceptual architecture

The XYZ grid is a quadtree over the Web Mercator plane. Zoom 0 is one tile covering the world; each level doubles the divisions, so zoom z has 4^z tiles. A tile is conventionally 256 × 256 pixels, which makes the ground resolution at zoom z roughly 156,543 / 2^z metres per pixel at the equator, shrinking with the cosine of latitude.

Two consequences follow immediately. The zoom level that matches a survey’s ground sample distance is fixed by that arithmetic rather than chosen — a 3 cm orthomosaic at 50° north corresponds to about zoom 22 — and the number of tiles quadruples per level, so the top two levels of any pyramid contain the overwhelming majority of the files.

Tile count and ground resolution by zoom level for one survey A chart over zoom levels sixteen to twenty-two for a one-square-kilometre survey. Tile count rises by a factor of four per level, from a handful at zoom sixteen to tens of thousands at zoom twenty-two, drawn on a logarithmic axis. Ground resolution halves per level, from about two and a half metres per pixel at zoom sixteen to under four centimetres at zoom twenty-two. A marker shows where the source orthomosaic's own ground sample distance of three centimetres sits, and a shaded region beyond it marks levels that add files without adding information because they are upsampling. upsampling 16 17 18 19 20 21 22 zoom level tiles (log scale) — ×4 per level ground resolution source GSD 3 cm ≈ zoom 21 Generating past the source resolution multiplies the file count without adding a single pixel of information.

Figure 1 — Where to stop. The last useful zoom is the one whose ground resolution matches the orthomosaic; every level beyond it quadruples the tile count while upsampling.

Step 1: Reproject into Web Mercator, once and deliberately

XYZ tiles are defined on EPSG:3857. A survey orthomosaic is almost never in it — it is in a UTM zone or a national grid — so tiling always involves a reprojection, and doing that reprojection as an explicit step rather than letting the tiler do it implicitly is worth the extra file.

import subprocess


def to_web_mercator(src: str, dst: str, resampling: str = "cubic") -> None:
    """Explicit reprojection to EPSG:3857, aligned to the tile grid.

    -tap snaps the output to whole pixels of the target resolution, which is
    what stops every tile boundary from falling mid-pixel and blurring.
    """
    subprocess.run(
        ["gdalwarp", "-t_srs", "EPSG:3857", "-r", resampling, "-tap",
         "-multi", "-wo", "NUM_THREADS=ALL_CPUS",
         "-co", "TILED=YES", "-co", "COMPRESS=DEFLATE",
         "-dstalpha", src, dst],
        check=True,
    )

-tap — target aligned pixels — is the flag whose absence produces the most complaints. Without it the warped raster’s origin is wherever the reprojected corner landed, which is generally not on a tile boundary, so every tile is resampled from a half-pixel offset and the whole set looks softer than the source. -dstalpha matters equally: an orthomosaic has an irregular footprint, and without an alpha band the area outside it becomes opaque black rather than transparent.

The resampling choice follows the raster type, exactly as for overviews: cubic for visual imagery, near for anything categorical.

Step 2: Choose between pre-rendering and serving on demand

Two architectures produce the same client experience, and they differ in almost every operational respect.

A pre-rendered pyramid is a directory of PNG or JPEG files produced by gdal2tiles.py. It needs no server beyond static hosting, it is trivially cacheable, and it is entirely predictable. It also takes hours to generate, produces hundreds of thousands of small files that are slow to copy and awkward to store, and must be regenerated in full whenever the orthomosaic changes.

Serving from a COG keeps a single file and renders tiles on request. There is no generation step at all, updates are a file replacement, and storage is the size of the raster. The cost is a running service, a cold-cache latency of tens of milliseconds per tile, and a dependency that must be operated.

Pre-rendered pyramid against on-demand tiles from a COG Two architectures compared across five rows. The pre-rendered pyramid requires a multi-hour generation step, produces hundreds of thousands of files, needs only static hosting, has no per-request latency, and must be fully regenerated on any update. The on-demand path has no generation step, stores one file, requires a running tile service, adds tens of milliseconds of cold-cache latency per tile, and updates by replacing the file. A note observes that for a survey deliverable that is revised during review, the update column usually decides the choice. pre-rendered pyramid on-demand from a COG generation files stored hosting latency update 2–6 hours 100 000+ small files static only none regenerate everything none one raster a running service 20–80 ms cold replace the file For a deliverable revised during client review, the update row usually settles it before any of the others matter.

Figure 2 — The comparison in operational terms. The pyramid wins on serving simplicity; the COG wins on everything to do with change, which is what a deliverable under review is full of.

# Pre-rendered: one command, several hours.
subprocess.run(
    ["gdal2tiles.py", "--zoom=14-21", "--processes=8", "--xyz",
     "--resampling=cubic", "--webviewer=none",
     "ortho_3857.tif", "tiles/"],
    check=True,
)

The --xyz flag is not optional for web use. Without it gdal2tiles.py writes the TMS scheme, whose row index counts from the south rather than the north — the tiles are correct and every one of them appears in the wrong place, mirrored vertically about the equator.

Step 3: Decide the zoom range from the data, not from habit

The minimum zoom is where the whole survey fits in a few tiles; the maximum is where the tile resolution matches the source GSD. Both are computable rather than conventional.

import math

import morecantile
import rasterio
from rasterio.warp import transform_bounds

TMS = morecantile.tms.get("WebMercatorQuad")


def zoom_range(src_path: str, max_overview_tiles: int = 4) -> tuple[int, int]:
    """(min_zoom, max_zoom) from the raster's extent and resolution."""
    with rasterio.open(src_path) as ds:
        west, south, east, north = transform_bounds(ds.crs, "EPSG:4326", *ds.bounds)
        gsd_m = abs(ds.transform.a)
        if ds.crs.to_epsg() == 4326:
            gsd_m *= 111_320 * math.cos(math.radians((north + south) / 2))

    lat = math.radians((north + south) / 2)
    # Native zoom: where tile resolution first equals or exceeds the source GSD.
    max_zoom = max(0, math.ceil(
        math.log2(156_543.03 * math.cos(lat) / gsd_m)))

    min_zoom = max_zoom
    while min_zoom > 0:
        tiles = list(TMS.tiles(west, south, east, north, [min_zoom - 1]))
        if len(tiles) > max_overview_tiles:
            break
        min_zoom -= 1
    return min_zoom, max_zoom

The cosine-of-latitude term is the part most hand-rolled versions omit. At 60° north a Web Mercator pixel covers half the ground it does at the equator, so a survey there reaches its native resolution one zoom level earlier — and generating that extra level costs four times the tiles for nothing.

What the client actually receives

A tile set is not self-describing. A directory of PNGs carries no extent, no zoom range, no attribution and no indication of which survey or which date it represents, so a client handed one has no way to load it without being told the parameters separately — and no way to tell two revisions apart.

Three small files fix that, and all three are generated rather than written by hand. A TileJSON document states the URL template, the zoom range, the bounds and the centre, which is enough for every common web map library to load the set with no further configuration. A metadata file recording the source raster, its checksum, the run identifier and the generation date makes a revision identifiable months later. And a minimal viewer page — twenty lines against a mapping library — turns “here is a tile directory” into something a non-technical client can open, which materially reduces the number of support conversations a delivery generates.

None of this is required for the tiles to work, which is exactly why it is routinely skipped and routinely missed. The generation step already knows every one of these values; writing them out costs a few lines at the end of the pipeline and makes the difference between a deliverable and a directory.

Parameter deep-dive

Parameter Type Typical Effect
--zoom range 14–21 Levels generated; the top level dominates count and time
tile size px 256 512 halves the file count and suits high-density displays
--resampling enum cubic near for categorical rasters, average for continuous
--xyz flag on Without it the row index is TMS and the map is mirrored
-tap flag on Aligns the warp to whole pixels; prevents soft tiles
-dstalpha flag on Irregular footprints stay transparent instead of black
--processes int cores Parallel tile generation; I/O-bound past about eight
JPEG quality int 80–90 Below 80, compression artefacts are visible on imagery

Format and quality, per raster type

The tile format is a per-deliverable decision with real consequences. JPEG tiles are three to five times smaller than PNG for photographic imagery and cannot carry transparency, which makes them wrong at the survey boundary unless a separate mask layer is served alongside. PNG carries an alpha channel and is lossless, which matters for anything a client will measure from rather than look at. WebP does both — alpha and photographic compression — at roughly half the JPEG size, and is now supported widely enough to be a reasonable default for new work.

For an orthomosaic that is purely visual, JPEG at quality 85 with a separate boundary polygon is the smallest useful answer. For a classification raster or an index product where pixel values carry meaning, PNG is the only correct choice, and the resampling must be nearest-neighbour all the way up the pyramid for the reasons set out under overview resampling.

The one combination to avoid is JPEG tiles of a categorical raster. The compression is lossy in a way that changes pixel values, so class codes come back subtly wrong — and unlike a blurred photograph, a wrong class code does not look wrong.

Caching, and why the pyramid’s shape decides your bandwidth bill

Tile requests are not uniform. A client opening a survey looks at the whole site once, then spends the rest of the session at high zoom over two or three areas of interest. That means the low zoom levels are requested by every visitor and the high ones by almost none — and it is the low levels that are cheap to cache, since there are only a handful of them.

The practical consequence is that a long cache lifetime on the low levels and a short one on the high levels serves both goals: repeat visitors never re-fetch the overview, while a revised orthomosaic propagates to detail views quickly. Where the tiles are versioned by a path segment carrying the run identifier, the lifetime can be effectively unlimited at every level, because a new revision is a new URL and nothing has to expire at all.

That last arrangement also removes the most common support question a tile deliverable generates, which is a client seeing a stale image after a revision and concluding the update was not delivered. A versioned path makes “which revision am I looking at?” answerable from the URL.

Verification and output inspection

Two failures account for most tile-set complaints, and both are checkable without a browser: tiles in the wrong place, and tiles that are blank.

import math
from pathlib import Path

import rasterio
from rasterio.warp import transform_bounds


def assert_tiles_cover_raster(tile_dir: Path, src_path: str, zoom: int) -> None:
    """Every tile the raster's extent implies must exist, and no others."""
    with rasterio.open(src_path) as ds:
        west, south, east, north = transform_bounds(ds.crs, "EPSG:4326", *ds.bounds)

    def deg_to_tile(lon: float, lat: float, z: int) -> tuple[int, int]:
        n = 2 ** z
        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 = deg_to_tile(west, north, zoom)      # north is the SMALL y in XYZ
    x1, y1 = deg_to_tile(east, south, zoom)

    missing = [(x, y)
               for x in range(x0, x1 + 1)
               for y in range(y0, y1 + 1)
               if not (tile_dir / str(zoom) / str(x) / f"{y}.png").exists()]
    assert not missing, f"{len(missing)} expected tiles missing at zoom {zoom}"

    present = {int(p.parent.name) for p in (tile_dir / str(zoom)).glob("*/*.png")}
    stray = present - set(range(x0, x1 + 1))
    assert not stray, f"tiles outside the raster extent at columns {sorted(stray)[:5]}"

Computing the expected tile range independently, from the raster’s own bounds, is what makes this a real check: it compares the tiler’s output against the geometry rather than against itself. The north → small y mapping in the third line is the XYZ convention, and getting it backwards is precisely the TMS mirroring the --xyz flag exists to prevent — so this check also catches that.

Troubleshooting

The map is mirrored vertically about the equator. TMS row indexing instead of XYZ. Regenerate with --xyz, or configure the client for TMS; do not attempt to rename the files, because the mapping is per zoom level.

Everything outside the survey footprint is black. No alpha band. Warp with -dstalpha and generate from the alpha-bearing raster; the tiler cannot invent transparency the source does not carry.

Tiles look softer than the orthomosaic. The warp was not pixel-aligned, so every tile samples the source at a half-pixel offset. Add -tap and regenerate.

Generation takes all night and fills the disk. The top zoom is one or two levels beyond the source resolution. Compute the native zoom rather than picking a round number; each surplus level is four times the tiles of the one before.

Tiles at the edges have a hairline seam. The source raster’s own edge pixels are partly transparent and the resampling is spreading them. Buffer the source slightly before warping, or use near at the edges.

A COG-backed tile server is slow on the first view and fast afterwards. That is the expected shape: the first request per region reads the relevant overview from storage. Pre-warm the low zoom levels for the survey extent after each deploy and the client never sees the cold path.

← DEM/DSM Generation & Raster Export Automation

XYZ and TMS number rows in opposite directions Two four-by-four tile grids side by side at the same zoom level. In the XYZ scheme the row index increases downward, so row zero is the northernmost. In the TMS scheme the row index increases upward, so row zero is the southernmost. A survey occupying the top-left tile is shown addressed as row zero under XYZ and as row three under TMS. A note explains that a client configured for one scheme and given tiles in the other renders every tile in the correct column and the mirrored row, which looks like a broken map rather than a naming problem. XYZ — row 0 at the north y = 0 1 2 3 survey TMS — row 0 at the south y = 3 2 1 0 survey Same tile, same column, different row number — which is why a mismatch mirrors the map rather than losing it.

Figure 3 — The one-flag difference that produces the most confusing tile bug. Both schemes are correct; only one of them matches what a web map client expects by default.