How to Export a Static Tile Pyramid for Offline Use

On this page

Problem statement

A static pyramid is tiles on disk, served by any file host or bundled into an application — no server process, no database, no runtime rendering.

The constraint is arithmetic. Tile count grows as $4^z$ , so each level costs roughly as much as all previous levels combined, three times over. For the Dutch coast from Katwijk to Rotterdam — the extent this page exports:

zoom       tiles
   8           2
  10          12
  12          96
  13         352
  14       1,302
total 8-14 1,796

Extending to zoom 15 would add 5,063 more. A full world pyramid to zoom 16 is 5,726,623,061 tiles — the sum of $4^z$ for every level up to 16 — which is why global static pyramids stop well below that.

Quick answer

Count before you generate:

import mercantile


def plan(bounds, zooms, bytes_per_tile=8_000):
    total = 0
    for zoom in zooms:
        count = len(list(mercantile.tiles(*bounds, zooms=[zoom])))
        total += count
        print(f"  z{zoom:<3} {count:9,} tiles  "
              f"{count * bytes_per_tile / 1e6:8.1f} MB")
    print(f"  total {total:,} tiles, {total * bytes_per_tile / 1e6:.0f} MB")
    return total


BOUNDS = (3.85, 51.95, 4.75, 52.35)          # Katwijk to Rotterdam, with sea to the west
plan(BOUNDS, range(8, 15))

  z8           2 tiles       0.0 MB
  z9           4 tiles       0.0 MB
  z10         12 tiles       0.1 MB
  z11         28 tiles       0.2 MB
  z12         96 tiles       0.8 MB
  z13        352 tiles       2.8 MB
  z14      1,302 tiles      10.4 MB
  total 1,796 tiles, 14 MB

Then decide between loose files and an archive. Loose files are simple and slow to copy — 1,796 files is fine, 1,796,000 is not. An MBTiles or PMTiles archive is one file and needs a reader.

Each level costs more than every level below it combined. The last zoom you add is most of the job.

A zoom range from the data, skip empties, deduplicate

1. Choose the zoom range from the data and the use

The maximum zoom should match the data's resolution. At 52° north, zoom 12 is 23.5 m per pixel and zoom 14 is 5.9 m. The elevation model here is 30 m, so zoom 14 already draws each source pixel about five times across; serving it higher only magnifies the same information.

The minimum zoom should be where the layer stops being useful — often the point at which everything has been generalised away.

2. Skip empty tiles

Over water, outside the data extent, or in sparse rural areas, most tiles have nothing in them. Writing a transparent PNG or an empty MVT for each multiplies the pyramid for no benefit.

Clients treat a 404 as "nothing here" correctly. Over the coastal extent here, 490 of the 1,796 tiles had no data at all and were never written (Example 1).

3. Deduplicate identical tiles

Where you do write empties — or where large areas are uniform — many tiles are byte-identical.

digest = hashlib.sha256(payload).hexdigest()
if digest in seen:
    os.link(seen[digest], path)          # hard link, not a copy

Hard links cost a directory entry. In an archive format, several tile ids can point at one blob. Here 529 tiles — open sea at a single elevation, rendered to identical bytes — became hard links (Example 1).

4. Choose loose files or an archive

	loose files	MBTiles	PMTiles
serving	any file host	needs a process	any host with ranges
copying	slow at scale	one file	one file
inspecting	trivial	SQL	a reader
updating	per tile	per row	rewrite

Loose files up to a few thousand; an archive beyond that. Copying a hundred thousand small files to a device takes far longer than copying one archive holding the same bytes.

5. Generate in parallel, by zoom or by tile range

Tiles are independent, so a process pool scales linearly. Splitting by zoom is simplest; splitting by tile range within a zoom balances better, because the top zoom has hundreds of times more tiles than the bottom.

Skipping empties is usually a larger saving than any compression setting.

Generating with a report

The data used on this page

The raster is one tile of the Copernicus 30 m digital elevation model, 52–53° N, 4–5° E: the Dutch coast from Katwijk to the Hook of Holland, the cities behind it, and a strip of the North Sea. The export extent reaches a little west and south of the tile, so some tiles have no data at all. Tiles are rendered by a small module — the overview-aware renderer from the COG tiles page, with a fixed stretch of −5 to 25 m — that returns PNG bytes, or None where a tile has no data.

import os

import requests

DEM = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N52_00_E004_00_DEM/"
       "Copernicus_DSM_COG_10_N52_00_E004_00_DEM.tif")
if not os.path.exists("dem_nl.tif"):
    with requests.get(DEM, stream=True, timeout=300) as reply:
        reply.raise_for_status()
        with open("dem_nl.tif", "wb") as out:
            for chunk in reply.iter_content(1 << 20):
                out.write(chunk)
print(f"dem_nl.tif {os.path.getsize('dem_nl.tif') / 1e6:.1f} MB")

dem_nl.tif 17.0 MB

%%writefile tilesrc.py
import io

import mercantile
import numpy as np
import rasterio
from matplotlib import colormaps
from PIL import Image
from rasterio.warp import Resampling, reproject

VMIN, VMAX = -5.0, 25.0                       # metres; one stretch for every tile


def overview_level(src, tile, size=256):
    west, _, east, _ = mercantile.bounds(tile)
    levels = [i for i, f in enumerate(src.overviews(1)) if src.res[0] * f <= (east - west) / size]
    return levels[-1] if levels else -1


def render_tile_from(source, tile, size=256):
    """PNG bytes for one XYZ tile, or None if the tile holds no data."""
    with rasterio.open(source) as src:
        level = overview_level(src, tile, size)
    with rasterio.open(source, **({"overview_level": level} if level >= 0 else {})) as src:
        data = np.full((size, size), np.nan, dtype="float32")
        reproject(rasterio.band(src, 1), data,
                  dst_transform=rasterio.transform.from_bounds(*mercantile.xy_bounds(tile), size, size),
                  dst_crs="EPSG:3857", resampling=Resampling.bilinear, dst_nodata=np.nan)
    valid = np.isfinite(data)
    if not valid.any():
        return None
    scaled = np.clip((np.nan_to_num(data) - VMIN) / (VMAX - VMIN), 0, 1)
    rgba = (colormaps["terrain"](scaled) * 255).astype("uint8")
    rgba[..., 3] = np.where(valid, 255, 0)
    buffer = io.BytesIO()
    Image.fromarray(rgba).save(buffer, format="PNG", optimize=True)
    return buffer.getvalue()

Writing tilesrc.py

Example 1 — generating with skipping, deduplication and a report

import gzip
import hashlib
import os
import mercantile


def export_pyramid(render, bounds, out_dir, zooms=range(10, 17),
                   extension="pbf", compress=True, deduplicate=True):
    """render(tile) -> bytes or None. Writes z/x/y.ext, skipping empties."""
    os.makedirs(out_dir, exist_ok=True)
    seen, written, skipped, linked, total_bytes = {}, 0, 0, 0, 0

    for zoom in zooms:
        tiles = list(mercantile.tiles(*bounds, zooms=[zoom]))
        for tile in tiles:
            payload = render(tile)
            if payload is None:
                skipped += 1
                continue
            if compress and extension == "pbf":
                payload = gzip.compress(payload, 6)

            directory = os.path.join(out_dir, str(tile.z), str(tile.x))
            os.makedirs(directory, exist_ok=True)
            path = os.path.join(directory, f"{tile.y}.{extension}")

            if deduplicate:
                digest = hashlib.sha256(payload).hexdigest()
                if digest in seen:
                    os.link(seen[digest], path)
                    linked += 1
                    continue
                seen[digest] = path

            with open(path, "wb") as handle:
                handle.write(payload)
            written += 1
            total_bytes += len(payload)

        print(f"  z{zoom}: {len(tiles):,} tiles considered")

    print(f"  {written:,} written, {linked:,} hard-linked duplicates, "
          f"{skipped:,} empty and skipped")
    print(f"  {total_bytes / 1e6:.1f} MB on disk, "
          f"mean {total_bytes / max(written, 1) / 1024:.1f} kB per tile")
    return written

import time

from tilesrc import render_tile_from

started = time.perf_counter()
export_pyramid(lambda tile: render_tile_from("dem_nl.tif", tile), BOUNDS, "pyramid",
               zooms=range(8, 15), extension="png", compress=False)
print(f"  {time.perf_counter() - started:.0f} s")

  z8: 2 tiles considered
  z9: 4 tiles considered
  z10: 12 tiles considered
  z11: 28 tiles considered
  z12: 96 tiles considered
  z13: 352 tiles considered
  z14: 1,302 tiles considered
  777 written, 529 hard-linked duplicates, 490 empty and skipped
  35.3 MB on disk, mean 44.3 kB per tile
  52 s

Of 1,796 tiles, 490 were empty and 529 duplicated another tile, so only 777 were written. The plan's guess of 8 kB a tile was five times too low for these terrain PNGs, which averaged 44.3 kB; plan with a measured tile, not a default.

Hard links only work within one filesystem and are invisible to most copying tools, which will expand them. For distribution, an archive format handles duplicates properly.

Measured in Example 1.

Example 2 — parallel generation with progress

The worker function goes into the module, so the process pool can import it:

%%writefile -a tilesrc.py
import os


def render_chunk(args):
    tiles, source, out_dir, extension = args
    written = 0
    for tile in tiles:
        payload = render_tile_from(source, tile)
        if payload is None:
            continue
        directory = os.path.join(out_dir, str(tile.z), str(tile.x))
        os.makedirs(directory, exist_ok=True)
        with open(os.path.join(directory, f"{tile.y}.{extension}"), "wb") as f:
            f.write(payload)
        written += 1
    return written

from concurrent.futures import ProcessPoolExecutor, as_completed

import tilesrc


def export_parallel(source, bounds, out_dir, zooms=range(10, 17),
                    workers=8, chunk=200, extension="png"):
    """Split every zoom into chunks so the work is balanced across workers."""
    tasks = []
    for zoom in zooms:
        tiles = list(mercantile.tiles(*bounds, zooms=[zoom]))
        for i in range(0, len(tiles), chunk):
            tasks.append((tiles[i:i + chunk], source, out_dir, extension))

    print(f"  {sum(len(t[0]) for t in tasks):,} tiles in {len(tasks)} chunks")
    total, done = 0, 0
    with ProcessPoolExecutor(workers) as pool:
        futures = [pool.submit(tilesrc.render_chunk, task) for task in tasks]
        for future in as_completed(futures):
            total += future.result()
            done += 1
            if done % 20 == 0:
                print(f"    {done}/{len(tasks)} chunks, {total:,} tiles written")
    print(f"  {total:,} tiles written")
    return total


import importlib

importlib.reload(tilesrc)
started = time.perf_counter()
export_parallel("dem_nl.tif", BOUNDS, "pyramid_parallel", zooms=range(8, 15), workers=8, chunk=100)
print(f"  {time.perf_counter() - started:.0f} s with 8 workers")

  1,796 tiles in 23 chunks
    20/23 chunks, 1,042 tiles written
  1,306 tiles written
  9 s with 8 workers

Chunking across zooms rather than by zoom balances the work: zoom 14 has 1,302 tiles here and zoom 8 has two, so one process per zoom would leave most workers idle almost immediately. Chunked, eight workers took 9 s for what one process did in 52.

Example 3 — bundling for offline use

import json
import os
import zipfile


def bundle(tile_dir, out_path, bounds, zooms, name, attribution="",
           tile_format="pbf"):
    """A zip with the tiles and a manifest, for shipping to a device."""
    manifest = {
        "name": name,
        "format": tile_format,
        "bounds": list(bounds),
        "minzoom": min(zooms),
        "maxzoom": max(zooms),
        "attribution": attribution,
        "scheme": "xyz",
    }

    count, total_bytes = 0, 0
    with zipfile.ZipFile(out_path, "w", zipfile.ZIP_STORED) as archive:
        archive.writestr("manifest.json", json.dumps(manifest, indent=2))
        for root, _, files in os.walk(tile_dir):
            for name_ in files:
                full = os.path.join(root, name_)
                archive.write(full, os.path.relpath(full, tile_dir))
                count += 1
                total_bytes += os.path.getsize(full)

    print(f"  {count:,} tiles, {total_bytes / 1e6:.1f} MB raw, "
          f"{os.path.getsize(out_path) / 1e6:.1f} MB zipped")
    return out_path

bundle("pyramid", "coast.zip", BOUNDS, range(8, 15), "Dutch coast elevation",
       attribution="Copernicus DEM GLO-30, © DLR e.V. 2010-2014 and © Airbus Defence and Space GmbH 2014-2018",
       tile_format="png")

  1,306 tiles, 35.6 MB raw, 35.7 MB zipped

ZIP_STORED rather than ZIP_DEFLATED is deliberate: the tiles are already compressed, so re-compressing costs time and saves nothing. It also lets a reader seek to a tile without decompressing.

The 529 hard-linked sea tiles went into the zip as 529 separate copies — cheap here, because a uniform tile is tiny, but it is the expansion the note on hard links warns about. Recording "scheme": "xyz" in the manifest is the line that prevents the row-numbering confusion on the consuming side.

The last level is most of the work

Why the last zoom level is most of the work

At each level the tile count quadruples, so level z has as many tiles as all levels below it combined, plus one.

For the measured extent: zooms 8 through 13 total 494 tiles, and zoom 14 alone is 1,302. Adding zoom 15 would add 5,063 more.

That makes the maximum zoom the single most consequential parameter. Choosing one level lower cuts the pyramid by roughly 75%.

Why to skip rather than write empty tiles

import io

import mapbox_vector_tile
import numpy as np
from PIL import Image

empty_png = io.BytesIO()
Image.fromarray(np.zeros((256, 256, 4), dtype="uint8")).save(empty_png, format="PNG", optimize=True)
empty_mvt = mapbox_vector_tile.encode([{"name": "layer", "features": []}])
print(f"transparent PNG {len(empty_png.getvalue())} bytes, empty MVT {len(empty_mvt)} bytes")

transparent PNG 334 bytes, empty MVT 14 bytes

Neither is large, and there can be very many.

More importantly, they carry a cost at every later stage: copying, uploading, listing, and — on object storage — a request charge each.

Clients handle 404 correctly: MapLibre and Leaflet treat a missing tile as an empty one. The only case for writing empties is a host that returns an HTML error page for 404, which some CDNs do by default.

Why deduplication beats compression on sparse data

A pyramid over a coastal extent contains many identical tiles: the same empty tile, or the same uniform sea colour — 529 of them in the export here.

Compression works within a tile and cannot see across them, so each identical tile is stored in full.

Deduplication — by hash, then hard links or archive-level pointers — removes them entirely. On sparse data that is a much larger saving than any compression setting, and the two are complementary.

Why an archive beats loose files at scale

Filesystems handle a few thousand small files well and hundreds of thousands badly. Directory listings slow down, copying takes minutes, and every file carries block-size overhead — a 200-byte tile on a 4 kB block leaves 95% of the block unused.

An archive is one file: fast to copy, exact in size, and with an internal index that is faster than a directory walk.

The threshold is roughly ten thousand tiles. Below that, loose files are simpler and easier to debug; above it, the archive pays for the extra reader.

Below a few thousand tiles, loose files are simpler. Above that, the archive pays for its reader.

Growth by four and hard links

Tile count grows as $4^z$ . The maximum zoom dominates everything.
Skip empty tiles; clients treat 404 as empty.
Deduplicate identical tiles — usually a bigger saving than compression.
Hard links do not survive copying. Use an archive for distribution.
ZIP_STORED for already-compressed tiles.
Chunk across zooms when parallelising, or most workers idle.
Record the scheme (xyz or tms) in a manifest.
Match the maximum zoom to the data's resolution, not to what looks smooth.

Internal links

Web map tiles explained: XYZ, zoom levels and the tile pyramid — where the counts come from
PMTiles and MBTiles explained — archive formats
How to build vector tiles from a GeoDataFrame in Python — generating vector tiles
How to render raster tiles from a COG in Python — generating raster tiles
My tiles are offset or in the wrong place — validating the output
How to speed up batch GIS jobs with parallel processing in Python — the parallel pattern
How to package GIS deliverables into zipped bundles with Python — bundling for delivery
Generalisation for zoom levels explained — what goes into each level

FAQ

How many tiles will my pyramid have?

Count them with mercantile.tiles. For the Dutch coast across zooms 8–14 it was 1,796, of which 1,302 are at zoom 14 alone.

Should I write empty tiles?

No. Clients treat a 404 as an empty tile, and empties cost copying, uploading and per-request charges.

How do I reduce the pyramid size?

Lower the maximum zoom — each level is most of the total — then skip empties and deduplicate identical tiles.

Loose files or an archive?

Loose files up to a few thousand tiles; an archive beyond that. Filesystems handle hundreds of thousands of small files badly.

How do I parallelise generation?

Chunk tiles across all zooms and hand chunks to a process pool. One process per zoom leaves most workers idle, because the top zoom has hundreds of times more tiles.

Should I zip the bundle with compression?

Use ZIP_STORED. The tiles are already compressed, so deflating again costs time and saves nothing.

What maximum zoom should I choose?

The one whose ground resolution matches the data. At 52° north zoom 12 is 23.5 m per pixel; serving 30 m data much beyond zoom 13 gives pixels finer than the information.