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.

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

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.

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.

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.