SpatiaLite & GeoPackage Automation with Python

How to Read GeoPackage Geometry Blobs in Python

Select the geometry column with plain sqlite3, read the two-byte GP magic and the flags byte of the GeoPackage Binary header to compute its length, slice…

Select the geometry column with plain sqlite3, read the two-byte GP magic and the flags byte of the GeoPackage Binary header to compute its length, slice the header off, and hand the trailing standard WKB to shapely.from_wkb — no GDAL required.

This page is part of the Spatial Data Serialization Patterns guide. The other pages in that section cover the write direction — encoding geometries to Well-Known Binary for insertion — while this walkthrough covers decoding: reading a GeoPackage geometry BLOB straight out of the SQLite file and turning it back into a Shapely object.

Why This Matters

A GeoPackage stores geometry not as raw WKB but as a GeoPackage Binary (GPB) BLOB: a small binary header defined by the OGC GeoPackage specification followed by standard WKB. If you SELECT geom FROM my_layer with sqlite3 and pass the bytes straight to shapely.from_wkb, the parse fails, because Shapely sees the GP magic bytes instead of a WKB byte-order marker. You must strip the envelope first.

Doing this by hand — instead of reaching for GDAL or GeoPandas — matters on locked-down field devices and in minimal containers where the GDAL stack is not installed. sqlite3 and shapely are lightweight, pure-enough dependencies, and understanding the GPB header is also what lets you read the bounding-box envelope for a fast pre-filter without decoding the full geometry.

The whole decoding problem is one of knowing where the header stops. Everything before that offset is GeoPackage bookkeeping; everything after it is a byte-for-byte standard WKB payload that Shapely already knows how to read. The header itself is eight fixed bytes followed by an optional envelope whose size the flags byte announces:

Anatomy of a GeoPackage Binary geometry BLOBA byte strip left to right. The first four fields form the eight-byte GeoPackageBinaryHeader: the two-byte GP magic, a one-byte version, a one-byte flags field, and a four-byte srs_id integer. After them comes a variable-length envelope of zero, thirty-two, forty-eight or sixty-four bytes, then the standard OGC WKB payload. A connector drops from the flags field to a second strip expanding its bits: bit 0 is the byte order, bits 1 to 3 are the envelope code, bit 4 marks an empty geometry, and bits 5 to 7 are reserved.8-byte GeoPackageBinaryHeadervariable-length payloadGPmagic · 2 B0x00version · 1 Bflags1 Bsrs_idint32 · 4 Benvelope0 / 32 / 48 / 64 BWKB payloadwhat from_wkb() readsbit 0byte orderbits 1–3envelope code 0–4bit 4empty geometrybits 5–7reservedheader_len = 8 + ENVELOPE_LEN[envelope code]
The flags byte is the key: bits 1–3 give the envelope code, which fixes where the WKB payload starts.

Prerequisites

  • Python 3.9+ (the sqlite3 module is in the standard library)
  • shapely 2.0+ (pip install shapely)
  • A GeoPackage file with at least one feature layer
  • Basic familiarity with struct and Python bytes slicing
  • Optional: a reading of the GeoPackage specification deep dive for the full header layout

Primary Method

python
# GeoPackage Binary (GPB) BLOB -> Shapely geometry, using stdlib sqlite3 + shapely 2.0+
import sqlite3
import struct
from shapely import from_wkb


def gpb_to_wkb(blob: bytes) -> bytes:
    """
    Strip the GeoPackage Binary header from a geometry BLOB and return
    the trailing standard WKB.

    Header layout (GeoPackage spec, clause on GeoPackageBinaryHeader):
      bytes 0-1 : magic  -> must be b'GP'
      byte  2   : version
      byte  3   : flags  -> bit 0 = byte order, bits 1-3 = envelope code
      bytes 4-7 : srs_id (int32)
      then      : envelope (0/32/48/64 bytes) then the WKB payload
    """
    if blob[:2] != b"GP":
        raise ValueError(f"Not a GeoPackage geometry BLOB: magic={blob[:2]!r}")

    flags = blob[3]
    little_endian = flags & 0x01           # bit 0: 1 = little-endian header ints
    envelope_code = (flags >> 1) & 0x07    # bits 1-3: envelope contents indicator

    # Envelope byte length keyed by the 3-bit indicator (each ordinate is a float64)
    envelope_bytes = {0: 0, 1: 32, 2: 48, 3: 48, 4: 64}[envelope_code]

    # 8-byte fixed header (magic, version, flags, srs_id) + optional envelope
    header_len = 8 + envelope_bytes
    return blob[header_len:]


def read_geometries(gpkg_path: str, layer: str, geom_col: str = "geom"):
    conn = sqlite3.connect(gpkg_path)
    try:
        cur = conn.execute(f'SELECT "{geom_col}" FROM "{layer}"')
        for (blob,) in cur:
            if blob is None:
                continue
            yield from_wkb(gpb_to_wkb(blob))
    finally:
        conn.close()

Every geometry the generator yields is a native Shapely object you can buffer, measure, or reproject. The whole path uses only the standard library plus Shapely.

Step-by-Step Walkthrough

1. Read a single BLOB with sqlite3

python
import sqlite3

conn = sqlite3.connect("field.gpkg")
row = conn.execute('SELECT "geom" FROM "survey_points" LIMIT 1').fetchone()
blob = row[0]
print(type(blob), len(blob), blob[:8])   # <class 'bytes'> ...

2. Confirm the magic and read the flags byte

The first two bytes must be b"GP". Byte 3 (index 3) is the flags byte that tells you the header’s byte order and which envelope, if any, is present:

python
assert blob[:2] == b"GP", "not a GeoPackage geometry"
version = blob[2]
flags = blob[3]
little_endian = bool(flags & 0x01)
envelope_code = (flags >> 1) & 0x07
print(version, little_endian, envelope_code)

3. Compute the envelope length

The 3-bit envelope indicator maps to a fixed byte count. Code 0 means no envelope; 1 is a 2D [minx, maxx, miny, maxy] envelope (four float64s = 32 bytes); 2 and 3 add a Z or M range (48 bytes); 4 carries both Z and M (64 bytes):

python
ENVELOPE_LEN = {0: 0, 1: 32, 2: 48, 3: 48, 4: 64}
envelope_bytes = ENVELOPE_LEN[envelope_code]
header_len = 8 + envelope_bytes

Each envelope ordinate is a float64, so the envelope grows in 16-byte steps as dimensions are added — one pair of min/max doubles per axis. Plotting the resulting total header length makes the cost of the richer envelopes obvious, and explains why writers that do not need Z or M ranges should leave them out:

Total GeoPackage Binary header length by envelope codeA horizontal bar chart with five rows, one per envelope code. Each bar is split into a fixed eight-byte segment and an envelope segment. Code 0 with no envelope totals 8 bytes; code 1 carrying an XY envelope totals 40 bytes; codes 2 and 3, carrying XYZ and XYM, total 56 bytes each; code 4 carrying XYZM totals 72 bytes.Bytes before the WKB payload startscode 0 — no envelope8 Bcode 1 — XY40 Bcode 2 — XY + Z range56 Bcode 3 — XY + M range56 Bcode 4 — XY + Z + M72 Bfixed 8 Benvelope (16 B per axis pair)
Codes 2 and 3 are the same length — one carries a Z range, the other an M range — so you cannot infer dimensionality from the byte count alone.

4. Slice off the header and decode the WKB

Everything after the header is ordinary WKB that Shapely reads directly:

python
from shapely import from_wkb

wkb = blob[header_len:]
geom = from_wkb(wkb)
print(geom.geom_type, geom.wkt[:60])

5. Optionally read the envelope for a cheap bounding box

When an envelope is present you can grab the bounds without decoding the geometry — useful for a fast spatial pre-filter:

python
import struct

if envelope_code == 1:
    endian = "<" if little_endian else ">"
    minx, maxx, miny, maxy = struct.unpack_from(endian + "4d", blob, 8)
    print(f"bbox = ({minx}, {miny}) .. ({maxx}, {maxy})")

6. Read the SRS id from the header

Bytes 4 through 7 hold the spatial reference id as a 32-bit integer, so you can group geometries by coordinate reference system without a join:

python
endian = "<" if little_endian else ">"
srs_id = struct.unpack_from(endian + "i", blob, 4)[0]
print("srs_id:", srs_id)

Verification

Cross-check the decoded geometry count and type against what the GeoPackage reports natively. If Fiona or GeoPandas is available, compare a couple of geometries; otherwise confirm the WKT round-trips:

python
from shapely import from_wkb, to_wkb

geom = from_wkb(gpb_to_wkb(blob))
assert from_wkb(to_wkb(geom)).equals(geom)   # decode/encode is stable
print("round-trip OK:", geom.geom_type)

You can also confirm the layer registration and geometry column name against the container’s registry:

sql
-- GeoPackage: find the geometry column name for a feature table
SELECT table_name, column_name, geometry_type_name, srs_id
FROM gpkg_geometry_columns;

Alternative Approaches or Edge Cases

GPB versus SpatiaLite extended WKB. A native SpatiaLite geometry BLOB is not GPB. It uses a different envelope: a leading 0x00 start byte, a trailing 0xFE end marker, and an interior format SpatiaLite calls extended WKB. The b"GP" magic check cleanly distinguishes the two — if the first byte is 0x00, you are looking at a SpatiaLite BLOB and need SpatiaLite’s AsBinary()/GeomFromWKB() functions (or the mod_spatialite reader) instead of this header math. Reading both formats correctly is the reason the magic check comes first.

GeoPackage Binary compared with a native SpatiaLite geometry BLOBTwo byte strips stacked for comparison. The upper strip is a GeoPackage Binary BLOB: GP magic, version, flags, srs_id, optional envelope, then standard WKB, with no terminator. The lower strip is a SpatiaLite BLOB: a 0x00 start byte, an endian byte, a srid integer, a minimum bounding rectangle of four doubles, a geometry class code, extended WKB, and a trailing 0xFE end marker. The first byte differs, which is what lets a reader tell them apart before parsing.GeoPackage Binary — first byte 0x47 (“G”)47 50GP magicver1 Bflags1 Bsrs_idint32envelopeoptionalstandard WKBno terminatorSpatiaLite BLOB — first byte 0x0000startendian1 Bsridint32MBR4 × float647Cclassextended WKBSpatiaLite dialectFEendblob[:2] == b'GP' decides which parser to use
Both containers wrap WKB, but only the GeoPackage layout is safe to slice at a computed offset — the SpatiaLite dialect embeds its own class codes.

The practical consequence is that a reader which has to cope with both containers should branch on the first byte and never on the file extension. A .sqlite file produced by ogr2ogr with -f GPKG holds GeoPackage geometries despite the extension, and a .gpkg file that someone loaded mod_spatialite into can carry a SpatiaLite-format column added after the fact. The magic bytes are the only reliable signal, and checking them costs two bytes per row.

Empty geometries. The flags byte also has an “empty geometry” bit (bit 4). When set, the WKB payload still parses but represents an empty geometry; guard downstream code that assumes non-empty coordinates. This matters most in aggregation: shapely.union_all over a list containing empty geometries succeeds silently, but an .area or .centroid call on the result of a fully-empty group raises or returns nan, and the traceback points at the aggregation rather than at the row that was empty on ingest.

Big-endian headers in the wild. Bit 0 of the flags byte can legitimately be 0, meaning the header integers are big-endian. Almost every writer emits little-endian, so code that hard-codes <i for the srs_id unpack works until it meets a file from a writer that did not. The envelope and the srs_id follow the header’s byte order; the WKB payload carries its own byte-order marker in its first byte and is independent of it. That split is the single most common source of “the geometry decodes but the SRS id is nonsense” reports, and it is why the code above derives endian from the flags byte rather than assuming.

Troubleshooting

shapely.errors.GEOSException: ParseException: Unknown WKB type or a garbled geometry

Cause: You passed the full GPB BLOB (including the GP header) to from_wkb, or you miscomputed the envelope length so the slice starts inside the header. Fix: Verify blob[:2] == b"GP" and recompute header_len = 8 + ENVELOPE_LEN[envelope_code] from the flags byte, matching the technique in reading a GeoPackage file header in Python.

ValueError: Not a GeoPackage geometry BLOB: magic=b'\x00\x01'

Cause: The layer stores SpatiaLite geometries, not GeoPackage Binary. Fix: Read those with SpatiaLite functions via mod_spatialite — see using sqlite3 with SpatiaLite functions — rather than the GPB header parser.

struct.error: unpack_from requires a buffer of at least N bytes

Cause: The BLOB is shorter than the header you tried to read, usually a truncated or NULL geometry. Fix: Skip None values and check len(blob) >= header_len before unpacking the envelope or SRS id.