Writing a Tile Pyramid into a GeoPackage

Run gdal_translate -of GPKG to create the deepest level and gdaladdo to build the coarser ones. The two decisions that actually matter are the tile format — JPEG for imagery, PNG for anything with sharp edges or transparency — and the maximum zoom, which should stop at the source’s ground resolution because every level below it quadruples the container for no additional detail.

This page belongs to the Raster Tiles & Coverage in GeoPackage guide.

Why This Matters #

A pyramid is where a container’s size is decided. The deepest level holds roughly three quarters of the tiles, the two deepest hold about nine tenths, and adding one more below the source resolution multiplies the file by about four while showing upsampled pixels that carry no information. On a field device with fixed storage, that single choice frequently decides whether the basemap ships at all.

Format matters almost as much. The same aerial imagery stored as PNG rather than JPEG can produce a container several times larger, at a quality difference nobody can see on a phone screen — while the reverse choice, JPEG for a cartographic basemap, produces visible artefacts along every line.

Prerequisites #

• GDAL 3.4+ with the GPKG driver, including raster support
• A georeferenced source raster whose reference system and extent you know
• Enough temporary disk for the intermediate — a pyramid is larger than its source
• Familiarity with the registry model from GeoPackage Specification Deep Dive

Primary Method #

#!/usr/bin/env bash
# Build a GeoPackage tile pyramid from a georeferenced raster
set -euo pipefail

SRC=${1:?usage: build-pyramid.sh <source.tif> <out.gpkg>}
OUT=${2:?}

# The deepest level. TILING_SCHEME commits the pyramid to a grid and a CRS.
gdal_translate -of GPKG "$SRC" "$OUT" \
    -co RASTER_TABLE=basemap \
    -co TILE_FORMAT=JPEG \
    -co QUALITY=80 \
    -co TILING_SCHEME=GoogleMapsCompatible

# The coarser levels, averaged down from the one below
gdaladdo -r average "$OUT" 2 4 8 16 32

# Confirm what was produced before anything else touches it
gdalinfo "$OUT" | sed -n '1,12p'

TILING_SCHEME=GoogleMapsCompatible is the option worth choosing consciously. It commits the pyramid to Web Mercator and the standard web-tile grid, which is what mobile map libraries expect; omitting it produces a pyramid in the source raster’s own reference system, which is more faithful to the data and less portable.

Step-by-Step Walkthrough #

1. Work out the maximum zoom from the source resolution #

The right maximum is the level whose ground pixel size matches the source’s. Anything deeper upsamples.

# Ground sample distance -> the deepest useful zoom, for a Web Mercator pyramid
import math

EQUATOR_M = 40_075_016.686
TILE_PX = 256


def deepest_useful_zoom(gsd_m: float, latitude_deg: float = 0.0) -> int:
    """The first zoom whose pixel size is at least as fine as the source."""
    for z in range(0, 25):
        px = EQUATOR_M * math.cos(math.radians(latitude_deg)) / (TILE_PX * 2 ** z)
        if px <= gsd_m:
            return z
    return 24


print(deepest_useful_zoom(0.25, latitude_deg=52))   # 25 cm imagery at 52°N

Building one level past that figure roughly quadruples the container. Building two past it multiplies it by sixteen, for pixels that were interpolated from the same source data.

2. Choose the format from what the raster contains #

# Photographic imagery — JPEG, and the quality is worth tuning
gdal_translate -of GPKG ortho.tif basemap.gpkg \
    -co RASTER_TABLE=basemap -co TILE_FORMAT=JPEG -co QUALITY=80

# Cartography, hillshade with alpha, anything with hard edges — PNG
gdal_translate -of GPKG cartographic.tif basemap.gpkg \
    -co RASTER_TABLE=basemap -co TILE_FORMAT=PNG

# Mixed content — PNG8 keeps transparency at a fraction of PNG's size
gdal_translate -of GPKG overlay.tif overlay.gpkg \
    -co RASTER_TABLE=overlay -co TILE_FORMAT=PNG8

JPEG has no alpha channel, so a pyramid with transparent regions cannot use it. GDAL’s TILE_FORMAT=AUTO picks per tile — JPEG where a tile is fully opaque, PNG where it is not — which is often the best answer for imagery with a nodata mask.

3. Build the overviews #

# Average is right for imagery; nearest preserves categorical values
gdaladdo -r average basemap.gpkg 2 4 8 16 32

# For a classified raster, averaging invents classes that do not exist
gdaladdo -r nearest landcover.gpkg 2 4 8 16

The resampling choice is not cosmetic. Averaging a land-cover raster whose values are class codes produces intermediate numbers that correspond to no class, and the coarser levels of the pyramid then show categories that were never in the data.

4. Add the feature layers into the same container #

The point of tiles in a GeoPackage is that they share a file with the vector data. Appending features is an ordinary -update write.

# One artefact: the basemap pyramid plus the survey layers over it
ogr2ogr -f GPKG -update basemap.gpkg parcels.shp -nln parcels \
    -t_srs EPSG:3857 -lco SPATIAL_INDEX=YES

ogr2ogr -f GPKG -update basemap.gpkg observations.shp -nln observations \
    -t_srs EPSG:3857 -lco SPATIAL_INDEX=YES

Reprojecting the features to the pyramid’s system is not required — each layer records its own — but it saves the client transforming every geometry on every frame, which on a field device is measurable.

5. Checkpoint before shipping #

A container written in WAL mode is three files, and a transfer that ships only the main one loses the newest commits. Fold the log in and switch the journal mode before the artefact leaves the build.

sqlite3 basemap.gpkg "PRAGMA wal_checkpoint(TRUNCATE); PRAGMA journal_mode=DELETE;"

Verification #

-- Every tile table must have a matrix set row and at least one matrix level
SELECT c.table_name,
       (SELECT count(*) FROM gpkg_tile_matrix_set s
         WHERE s.table_name = c.table_name) AS has_set,
       (SELECT count(*) FROM gpkg_tile_matrix m
         WHERE m.table_name = c.table_name) AS levels
FROM gpkg_contents c WHERE c.data_type = 'tiles';

-- No tile may sit outside the grid its level declares
SELECT t.zoom_level, count(*) AS out_of_range
FROM basemap t
JOIN gpkg_tile_matrix m
  ON m.table_name = 'basemap' AND m.zoom_level = t.zoom_level
WHERE t.tile_column >= m.matrix_width OR t.tile_row >= m.matrix_height
GROUP BY t.zoom_level;

The second query is the one that catches a row-convention mismatch before anyone opens the file. Then confirm the tiles themselves decode, since a BLOB that is not an image satisfies every structural check:

import io, sqlite3
from PIL import Image, UnidentifiedImageError

conn = sqlite3.connect("file:basemap.gpkg?mode=ro", uri=True)
for (z, tw, th) in conn.execute(
    "SELECT zoom_level, tile_width, tile_height FROM gpkg_tile_matrix "
    "WHERE table_name='basemap' ORDER BY zoom_level"
):
    for (blob,) in conn.execute(
        "SELECT tile_data FROM basemap WHERE zoom_level=? LIMIT 5", (z,)
    ):
        try:
            img = Image.open(io.BytesIO(blob))
        except UnidentifiedImageError:
            raise AssertionError(f"undecodable tile at zoom {z}")
        assert img.size == (tw, th), f"zoom {z}: {img.size} != ({tw}, {th})"

Alternative Approaches or Edge Cases #

Building the registry by hand. Possible, and rarely worth it. The matrix rows must be internally consistent — grid dimensions, tile size and pixel size all describing the same geometry — and getting one wrong produces a pyramid that renders offset. Where a custom scheme is genuinely required, generate the rows from a single source of truth rather than typing them.

Updating part of a pyramid. Tiles are rows, so replacing the tiles covering a changed area is an ordinary transaction. The work is propagating the change upward: every coarser tile above the changed ones is a downsampled composite that is now stale, and recomputing that chain is what keeps the pyramid consistent.

Very large sources. gdal_translate builds the deepest level in one pass and needs temporary space proportional to the output. For a source too large for that, build in tiles with gdal_retile and merge, or build directly at a coarser maximum zoom and accept the resolution.

Troubleshooting #

ERROR 1: Only single band, or 3 or 4 band on translate #

Cause: The source has a band count the GeoPackage tile format cannot represent — a multispectral raster, for instance. Fix: Select bands explicitly with -b 1 -b 2 -b 3, and produce a separate pyramid for any additional band you need.

The container is far larger than expected #

Cause: PNG on photographic imagery, or a maximum zoom past the source resolution. Fix: Check both. Rebuilding with JPEG and one level shallower routinely reduces a pyramid by an order of magnitude with no visible difference.

Coarse levels look blocky or show wrong colours #

Cause: -r nearest on continuous imagery, or -r average on categorical data. Fix: Match the resampling to the content, as in step 3, and rebuild the overviews — gdaladdo -clean removes them first.

Frequently Asked Questions #

Related #

• Raster Tiles & Coverage in GeoPackage — parent guide: the pyramid model and its registry
• How to Read a Raster Tile from a GeoPackage in Python — the read path for what this builds
• GeoPackage Tile Matrix Sets Explained — the registry rows the build writes
• How to Convert Shapefiles with ogr2ogr on the Command Line — the vector half of the same container