Raster Tiles & Coverage in GeoPackage

Most work with GeoPackage stops at feature tables, and that leaves half the format unused. The same container can hold a raster basemap as a pyramid of image tiles — which is what makes a single .gpkg a complete offline map rather than a layer that needs a basemap from somewhere else. On a field device with no network, that difference is the whole product.

This guide is part of the Python Integration & Database Workflows section. It covers the tile side of the specification: how a pyramid is addressed, which registry tables describe it, and how to read and write tiles from Python with and without GDAL.

Prerequisites #

• Python 3.9+ with sqlite3 — tiles are BLOBs in an ordinary table, so no extension is required to read them
• Pillow 10+ for decoding tile images, or any PNG/JPEG decoder you prefer
• GDAL/OGR 3.4+ if you intend to build a pyramid rather than only read one
• Familiarity with the registry model from GeoPackage Specification Deep Dive
• A clear target reference system — a tile pyramid commits to one, and changing it later means rebuilding

Concept & Specification Reference #

A tile pyramid is a stack of grids. Each zoom level covers the same ground with more, smaller tiles, and every tile is addressed by three numbers: its zoom level, its column, and its row. Four registry tables describe the arrangement.

The addressing detail that catches everyone is row direction. GeoPackage numbers tile_row from the top of the extent downward, which matches the common web-tile convention but is the opposite of the TMS convention some tools emit. A pyramid built with the wrong convention renders as a vertically mirrored map, level by level — visually obvious once you know to look for it, and easy to mistake for a projection problem if you do not.

Step-by-Step Implementation #

1. Discover what the container actually holds #

A container can hold feature layers, tile layers, or both. gpkg_contents distinguishes them by data_type, and reading it first avoids opening a tile table as though it were features.

# List feature and tile layers, and the zoom range of each pyramid
import sqlite3

conn = sqlite3.connect("file:basemap.gpkg?mode=ro", uri=True)

for name, dtype in conn.execute(
    "SELECT table_name, data_type FROM gpkg_contents ORDER BY data_type, table_name"
):
    print(f"{dtype:>8}  {name}")

for name, zmin, zmax, n in conn.execute("""
    SELECT c.table_name,
           min(m.zoom_level), max(m.zoom_level), count(*)
    FROM gpkg_contents c
    JOIN gpkg_tile_matrix m ON m.table_name = c.table_name
    WHERE c.data_type = 'tiles'
    GROUP BY c.table_name
"""):
    print(f"{name}: zoom {zmin}–{zmax} across {n} matrix levels")

2. Read the tile matrix for a zoom level #

Everything needed to locate a tile geographically lives in gpkg_tile_matrix_set and gpkg_tile_matrix. The set gives the pyramid’s extent and reference system; the matrix gives the grid at each level.

# Everything needed to map a ground coordinate to a tile address
row = conn.execute("""
    SELECT s.min_x, s.min_y, s.max_x, s.max_y, s.srs_id,
           m.matrix_width, m.matrix_height,
           m.tile_width, m.tile_height,
           m.pixel_x_size, m.pixel_y_size
    FROM gpkg_tile_matrix_set s
    JOIN gpkg_tile_matrix m ON m.table_name = s.table_name
    WHERE s.table_name = ? AND m.zoom_level = ?
""", ("basemap", 12)).fetchone()

min_x, min_y, max_x, max_y, srs_id, mw, mh, tw, th, px, py = row
print(f"zoom 12: {mw}×{mh} tiles of {tw}×{th} px at {px} units/px in EPSG:{srs_id}")

3. Convert a ground coordinate to a tile address #

With the matrix in hand, the mapping is arithmetic — and the only subtlety is the downward row direction.

def tile_address(x: float, y: float, zoom: int, matrix) -> tuple[int, int]:
    """Map a ground coordinate to (tile_column, tile_row) at a zoom level."""
    min_x, min_y, max_x, max_y, mw, mh = matrix

    span_x = (max_x - min_x) / mw
    span_y = (max_y - min_y) / mh

    col = int((x - min_x) // span_x)
    # tile_row counts DOWN from max_y — this is the line people get wrong
    row = int((max_y - y) // span_y)

    if not (0 <= col < mw and 0 <= row < mh):
        raise ValueError(f"({x}, {y}) is outside the pyramid extent at zoom {zoom}")
    return col, row

4. Fetch and decode a tile #

Tiles are ordinary BLOBs holding PNG or JPEG data, so reading one needs nothing beyond the standard library plus an image decoder.

# Fetch one tile and decode it
import io
from PIL import Image

blob = conn.execute("""
    SELECT tile_data FROM basemap
    WHERE zoom_level = ? AND tile_column = ? AND tile_row = ?
""", (12, col, row)).fetchone()

if blob is None:
    print("no tile stored at that address — the pyramid may be sparse")
else:
    img = Image.open(io.BytesIO(blob[0]))
    print(img.format, img.size, img.mode)

A missing tile is normal rather than exceptional. Pyramids are commonly sparse: a container covering a coastal survey area stores no tiles for the open sea, and a reader that treats None as an error will fail constantly on perfectly good data.

5. Write a pyramid #

Building a pyramid by hand is possible and rarely worth it. gdal_translate and gdaladdo produce a conforming container in two commands, and getting the matrix rows right by hand is fiddly enough that the driver route is the sensible default.

# Build a GeoPackage tile pyramid from a georeferenced raster
gdal_translate -of GPKG orthophoto.tif basemap.gpkg \
    -co RASTER_TABLE=basemap \
    -co TILE_FORMAT=JPEG \
    -co QUALITY=80 \
    -co TILING_SCHEME=GoogleMapsCompatible

# Add the coarser levels of the pyramid
gdaladdo -r average basemap.gpkg 2 4 8 16 32

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

6. Append feature layers into the same container #

The point of tiles in GeoPackage is that they share a container with features. Adding a feature layer to a container that already holds a pyramid is an ordinary append, following the pattern in How to Append Layers to an Existing GeoPackage with Fiona.

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

Reprojecting the features to match the pyramid’s reference system is not required by the format — each layer records its own — but it saves the client transforming every geometry at render time, which on a field device is the difference between a map that pans smoothly and one that does not.

Validation & Verification #

A pyramid can be structurally valid and still unusable, so check both the registry and the tiles.

-- 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 matrix 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;

A count above zero in the second query means the pyramid was written with a different row convention or a mismatched extent, and it is the check that catches the mirrored-map problem before anyone opens the file.

Sampling the tiles themselves is worth doing too, because a BLOB that is not a decodable image passes every structural check:

# Sample a handful of tiles per level and confirm they decode at the declared size
import io, random
from PIL import Image, UnidentifiedImageError

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"
):
    sample = conn.execute(
        "SELECT tile_data FROM basemap WHERE zoom_level=? LIMIT 5", (z,)
    ).fetchall()
    for (blob,) in sample:
        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} != declared ({tw}, {th})"

Common Failure Modes & Fixes #

The map renders upside down, level by level #

Diagnosis: The pyramid was written with TMS row numbering, which counts upward from the bottom, while GeoPackage counts downward from the top. Fix: Recompute tile_row as matrix_height - 1 - tile_row for every tile at every level. This is a pure key rewrite — no image data changes — but it must be applied consistently across all levels, since matrix_height differs per level.

Tiles are missing over part of the extent #

Diagnosis: Usually not a fault. Pyramids are sparse, and a container built from a raster with nodata regions stores no tiles there. Fix: Confirm against the source raster’s coverage before treating it as a defect; render a missing tile as transparent rather than failing.

The container is enormous relative to the source raster #

Diagnosis: The tile format is PNG where JPEG would do. Lossless PNG on aerial imagery commonly produces containers several times larger than JPEG at a quality nobody can distinguish on a phone screen. Fix: Rebuild with TILE_FORMAT=JPEG for imagery, keeping PNG for anything with sharp edges or transparency — cartographic basemaps, overlays, hillshades with an alpha channel.

A conforming client sees no tile layer at all #

Diagnosis: The tile table exists but gpkg_contents has no row with data_type = 'tiles', or the matrix set row is missing. This is the tile equivalent of an unregistered feature table. Fix: Insert the missing registry rows; the table itself is fine, and no image data needs rewriting.

Zoom levels beyond a point show no more detail #

Diagnosis: The pyramid was built with overviews below the source resolution, so the deepest levels are upsampled from the same pixels. Fix: Stop the pyramid at the level whose pixel_x_size matches the source ground sample distance. Levels beyond it quadruple storage per step and add no information — commonly the single largest avoidable component of a container’s size.

Performance Notes #

Storage in a pyramid is dominated by its deepest level. Because each level holds four times the tiles of the one above, the last level alone is about three quarters of the total and the last two are roughly nine tenths. That makes the maximum zoom the most consequential decision in the whole build: dropping one level from a pyramid removes most of its size, and adding one below the source resolution adds size for no detail.

Reading is fast for a different reason. A tile lookup is a primary-key hit on (zoom_level, tile_column, tile_row), which is an ordinary B-tree index — no spatial index is involved, and none is needed. What does hurt is a query pattern that scans: fetching “all tiles at zoom 15” to find the ones in view reads the whole level, where computing the address range and fetching by key reads only what is drawn.

Finally, tile blobs are already compressed. Storing them in a container that is itself compressed, or copying them through a compressing transport, spends CPU for a percent or two. Keep the compression in the image format where it belongs.

Child Pages #

Pages in this section go deeper on individual tile tasks:

• How to Read a Raster Tile from a GeoPackage in Python — address a tile and decode it with no GDAL dependency
• Writing a Tile Pyramid into a GeoPackage — building the pyramid and its registry rows
• GeoPackage Tile Matrix Sets Explained — the two matrix tables in detail
• How to Serve GeoPackage Tiles Offline — a minimal local tile endpoint over a container
• Comparing MBTiles and GeoPackage for Basemaps — two SQLite tile formats, and when each wins

Frequently Asked Questions #

Related #

• Python Integration & Database Workflows — parent section: the full Python spatial SQLite stack
• GeoPackage Specification Deep Dive — the registry model tile tables share with feature tables
• How to Append Layers to an Existing GeoPackage with Fiona — adding feature layers beside a pyramid
• Offline-First Sync Strategies — distributing the resulting container to field devices