SpatiaLite & GeoPackage Automation with Python

How to Read a Raster Tile from a GeoPackage in Python

Read the pyramid's extent and grid dimensions from gpkgtilematrixset and gpkgtilematrix, convert your ground coordinate into a (tilecolumn, tilerow) pair…

Read the pyramid’s extent and grid dimensions from gpkg_tile_matrix_set and gpkg_tile_matrix, convert your ground coordinate into a (tile_column, tile_row) pair — remembering that rows count downward from the top of the extent — and SELECT tile_data from the tile table. The BLOB is a PNG or JPEG, so any image decoder opens it and no GDAL is involved anywhere.

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

Why This Matters

A field application that renders a basemap needs one thing from the container: the image covering the area on screen. Going through GDAL for that pulls in a large dependency to do work that is two queries and a decoder call — and on a mobile or embedded target, the dependency is frequently the reason the basemap ends up in a separate file instead of the GeoPackage where it belongs.

The arithmetic is also worth understanding rather than delegating. The one detail that catches everyone is row direction: GeoPackage numbers tile_row from the top of the extent downward, matching the common web-tile convention and inverting the TMS one. Getting it wrong produces a map that is mirrored vertically at every zoom level, which reads as a projection problem and is not.

Prerequisites

  • Python 3.9+ with sqlite3 — nothing else is required to fetch the BLOB
  • Pillow 10+, or any PNG/JPEG decoder, to turn the bytes into pixels
  • A GeoPackage containing a tile table, per Writing a Tile Pyramid into a GeoPackage
  • The tile table’s name, or the enumeration query that finds it

Primary Method

python
# Ground coordinate -> tile address -> decoded image, with no GDAL
import io
import sqlite3
from dataclasses import dataclass
from PIL import Image


@dataclass(frozen=True)
class Matrix:
    min_x: float; min_y: float; max_x: float; max_y: float
    srs_id: int
    width: int; height: int          # grid dimensions, in tiles
    tile_w: int; tile_h: int         # tile dimensions, in pixels


def read_matrix(conn, table: str, zoom: int) -> Matrix:
    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
        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 = ?
    """, (table, zoom)).fetchone()
    if row is None:
        raise LookupError(f"no tile matrix for {table!r} at zoom {zoom}")
    return Matrix(*row)


def tile_address(m: Matrix, x: float, y: float) -> tuple[int, int]:
    """Ground coordinate -> (tile_column, tile_row). Rows count DOWN from max_y."""
    span_x = (m.max_x - m.min_x) / m.width
    span_y = (m.max_y - m.min_y) / m.height

    col = int((x - m.min_x) // span_x)
    row = int((m.max_y - y) // span_y)      # the line people get wrong

    if not (0 <= col < m.width and 0 <= row < m.height):
        raise ValueError(f"({x}, {y}) is outside the pyramid extent")
    return col, row


def read_tile(conn, table: str, zoom: int, x: float, y: float):
    m = read_matrix(conn, table, zoom)
    col, row = tile_address(m, x, y)

    blob = conn.execute(
        f'SELECT tile_data FROM "{table}" '
        "WHERE zoom_level = ? AND tile_column = ? AND tile_row = ?",
        (zoom, col, row),
    ).fetchone()

    if blob is None:
        return None, (col, row)          # sparse pyramids are normal
    return Image.open(io.BytesIO(blob[0])), (col, row)

Returning None rather than raising for a missing tile is deliberate. Pyramids are commonly sparse — a container covering a coastal survey stores nothing for open water — and a reader that treats absence as an error fails constantly on perfectly good data.

Ground coordinate to tile addressThe pyramid extent is divided into a grid whose dimensions come from the tile matrix for that zoom level. The column is the horizontal distance from the extent's minimum x, divided by the tile span. The row is the distance downward from the extent's maximum y, not upward from the minimum, which is the difference between the GeoPackage convention and the TMS one.extent + grid sizefrom the two matrix tablesper zoom levelcolumn(x − min_x) ÷ span_xleft to right, as expectedSELECTby the three keysa primary-key hitrow counts DOWN from max_y(max_y − y) ÷ span_yusing min_y instead mirrors the map at every level
Two of the three steps are unsurprising. The amber box is the entire difficulty of tile addressing.

Step-by-Step Walkthrough

1. Find the tile tables

sql
-- Tile layers, with the zoom range each covers
SELECT c.table_name,
       min(m.zoom_level) AS zmin,
       max(m.zoom_level) AS zmax,
       count(*)          AS levels
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;

2. Open read-only

Nothing here writes, and opening read-only means the reader cannot be the cause of a lock or a stray modification to a container a field application is also using.

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

3. Fetch a whole window rather than a tile at a time

Rendering a viewport needs several tiles, and one query for the range beats one query per tile. The addresses are contiguous, so the range is expressible directly.

python
def read_window(conn, table, zoom, x0, y0, x1, y1):
    """Every tile covering a ground rectangle, in one query."""
    m = read_matrix(conn, table, zoom)
    c0, r1 = tile_address(m, min(x0, x1), min(y0, y1))   # bottom-left -> max row
    c1, r0 = tile_address(m, max(x0, x1), max(y0, y1))   # top-right   -> min row

    rows = conn.execute(
        f'SELECT tile_column, tile_row, tile_data FROM "{table}" '
        "WHERE zoom_level = ? AND tile_column BETWEEN ? AND ? "
        "  AND tile_row BETWEEN ? AND ?",
        (zoom, c0, c1, r0, r1),
    ).fetchall()
    return {(c, r): blob for c, r, blob in rows}

The row bounds swap relative to the ground coordinates, because the ground y increases upward and tile_row increases downward. Deriving both from tile_address rather than computing them separately keeps that inversion in one place.

4. Decode, and handle the absent tile

python
tiles = read_window(conn, "basemap", 14, 400_000, 300_000, 402_000, 302_000)

for (col, row) in expected_addresses:
    blob = tiles.get((col, row))
    if blob is None:
        canvas.paste(TRANSPARENT, position_of(col, row))   # normal, not an error
        continue
    canvas.paste(Image.open(io.BytesIO(blob)), position_of(col, row))

5. Read the pixel size when you need scale

gpkg_tile_matrix records pixel_x_size and pixel_y_size per level, which is what tells you the ground resolution — and therefore which zoom level to request for a given screen scale.

python
def best_zoom(conn, table: str, target_m_per_px: float) -> int:
    """The deepest level whose pixels are no coarser than the target."""
    rows = conn.execute(
        "SELECT zoom_level, pixel_x_size FROM gpkg_tile_matrix "
        "WHERE table_name = ? ORDER BY zoom_level", (table,)
    ).fetchall()
    candidates = [z for z, px in rows if px <= target_m_per_px]
    return candidates[0] if candidates else rows[-1][0]
One range query instead of one query per tileRendering a viewport at zoom fourteen might need twelve tiles. Fetching them individually costs twelve primary-key lookups and twelve round trips through the driver. Computing the address range from the two corner coordinates and fetching with a single BETWEEN query costs one, and returns the same rows. The saving grows with viewport size and matters most on a device where every call is expensive.one query per tiletwelve lookups for a viewporttwelve driver round tripssimple to writecosts scale with the viewportone range querytwo BETWEEN clausesone round tripsame rows returnedconstant cost per framethe tile table's primary key already orders by zoom, column and row
The range query is not a clever optimisation — it is what the primary key was designed for.

Verification

Confirm the address arithmetic against the extent, and confirm the tiles decode at the size the matrix declares.

python
# Corners of the extent must map to the corner tiles
m = read_matrix(conn, "basemap", 14)

assert tile_address(m, m.min_x + 0.001, m.max_y - 0.001) == (0, 0)
assert tile_address(m, m.max_x - 0.001, m.min_y + 0.001) == (m.width - 1, m.height - 1)
python
# Sampled tiles must decode, at the declared dimensions
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 3", (z,)
    ):
        img = Image.open(io.BytesIO(blob))
        assert img.size == (tw, th), f"zoom {z}: {img.size} != declared ({tw}, {th})"

The first assertion is the one that catches an inverted row calculation, because a mirrored implementation maps the top-left corner to (0, height - 1) instead of (0, 0).

The corner assertion that catches an inverted row calculationUnder the GeoPackage convention, the top-left corner of the extent maps to tile column zero, row zero, and the bottom-right corner maps to the last column and last row. An implementation that counts rows upward from the minimum y instead maps the top-left corner to column zero and the last row, which the corner assertion catches immediately without needing a rendered image.correcttop-left of extent → (0, 0)bottom-right → (w−1, h−1)rows count down from max_ythe GeoPackage conventioninvertedtop-left of extent → (0, h−1)bottom-right → (w−1, 0)rows count up from min_ythe TMS conventiontwo assertions on the corners settle it without rendering anything
Both conventions are used in the wild; only one is what a GeoPackage means.

Alternative Approaches or Edge Cases

Reading through GDAL. gdal.Open("GPKG:basemap.gpkg:basemap") gives a raster dataset with the usual windowed-read interface, handles the addressing, and is the right choice when you need resampling, reprojection or a mosaic. For “give me the tile at this address”, it is a large dependency for two queries.

Tiles in a different reference system than the features. Perfectly legal — each registers its own — and it means the client transforms geometry to match the basemap on every frame. Reprojecting the feature layers at build time to match the pyramid removes that, at the cost of storing them in a system that may be poor for measurement.

Web Mercator pyramids and the standard scheme. A pyramid built with TILING_SCHEME=GoogleMapsCompatible uses the grid every web map library expects, so the addressing above reduces to the familiar z/x/y arithmetic. Pyramids in a national grid do not, which is exactly why the matrix tables exist rather than a formula.

Troubleshooting

The map is mirrored vertically

Cause: The row calculation counts up from min_y instead of down from max_y. Fix: Use (max_y - y) // span_y. The corner assertion in the verification section catches this immediately.

Every tile query returns nothing

Cause: The wrong zoom level, or the coordinate is outside the pyramid extent. Fix: Print the matrix bounds and the computed address; an address outside the grid raises from tile_address, but a valid address at a zoom level with no stored tiles returns None silently.

UnidentifiedImageError on a tile

Cause: The BLOB is not a decodable image — a truncated write, or a container built by a tool that stored something else. Fix: Sample tiles across every level as in the verification step; an undecodable tile passes every structural check the format defines.

Frequently Asked Questions

Do I need the spatial extension to read tiles?

No. Tile tables hold ordinary integers and BLOBs, the matrix tables are ordinary metadata, and nothing in the read path calls a spatial function. This is one of the few parts of a GeoPackage that is fully accessible from a plain sqlite3 connection, which makes it well suited to restricted environments.

Is there an index on the tile table?

Yes — the primary key on (zoom_level, tile_column, tile_row), which is what makes both the single-tile lookup and the range query efficient. No R-tree is involved, and none would help: the client computes the addresses it wants arithmetically rather than asking which tiles overlap a rectangle.

How should a missing tile be rendered?

As transparent, in almost every case. A sparse pyramid is a design decision — the producer stored nothing where there was nothing to show — so a gap is information rather than an error. Rendering a placeholder or an error tile turns a correct container into one that looks broken.