Versioning FlatGeobuf Streams for Web Delivery

FlatGeobuf packs a spatial index directly into the file and lays features out for progressive streaming, which makes it both fast to serve to a web map and — once its output is made deterministic — cheap to version. This page covers producing byte-stable FlatGeobuf for clean diffs and versioning per-tile files for delivery, as part of Cloud-Native Spatial Formats for Versioned Pipelines.

Concept & Context

FlatGeobuf is a single-file binary vector format built around a static, packed Hilbert R-tree written at the head of the file, ahead of the feature records. Because the index lives in the file, a client that supports HTTP range requests can read the header, binary-search the tree for a bounding box, and fetch only the byte ranges holding the matching features — no sidecar .qix or database connection required. The feature section is a flat, sequential run of length-prefixed records, so a reader can also stream the whole layer front-to-back as it downloads. That combination of in-file index plus streamable layout is what makes it a natural web-delivery format.

The delivery advantages compound in a versioned pipeline. Because a tile is one self-contained file, a content delivery layer can cache it by its content hash and invalidate only the tiles that changed on a release. A map front-end never re-downloads unchanged geography, and a reviewer inspecting a commit sees precisely which tiles moved — a per-tile hash manifest reads like a spatial changelog. Contrast this with a shapefile set, where a single edit rewrites the .shp, .dbf, and .shx together and the sidecar .qix index must be regenerated, so the whole bundle re-hashes on any change.

For versioning, the streamable layout is a double-edged sword. FlatGeobuf writes one contiguous binary; there is no internal chunk grid the way Zarr has, so a naive rewrite produces a brand-new file and a full re-hash on every edit. Two disciplines recover cheap versioning. First, make the output deterministic so an unchanged layer serialises to identical bytes and a content hash is meaningful. Second, tile the layer on a fixed grid so each tile is a small self-contained FlatGeobuf that only re-hashes when its own features change — the vector analogue of the chunk manifest used for rasters, and a close cousin of the file-level pointer synchronization for raster datasets pattern.

Core Pipeline

  1. Normalize inputs. Reproject to one delivery CRS (commonly EPSG:4326 for web maps), round coordinates to a fixed decimal precision, and enforce a stable attribute field order and types.
  2. Sort features by a stable spatial key. Order by a Hilbert or Morton curve over the geometry centroid so feature order is reproducible and spatially coherent — this also improves the packed index quality.
  3. Write FlatGeobuf with the index enabled. Use the GDAL FlatGeobuf driver so the packed R-tree is embedded, producing a self-contained streamable file.
  4. Hash and diff. Take a content hash of each output file and compare against the prior version’s manifest to find which tiles changed.
  5. Publish per tile. Upload only the changed tile files to object storage and commit the updated manifest, so a web map fetches the current tile set and the version history stays small.

Working Implementation

The script reprojects and normalizes a GeoDataFrame, sorts features by a Hilbert key for reproducibility, tiles them on a fixed grid, and writes one deterministic FlatGeobuf per non-empty tile using GDAL via geopandas/pyogrio. It reports a content hash per tile for the version manifest.

#!/usr/bin/env python3
"""
tile_flatgeobuf.py

Produce deterministic, per-tile FlatGeobuf files for web delivery and
report a content hash per tile for incremental versioning.

Requirements: geopandas>=0.14, pyogrio, shapely>=2.0, GDAL FlatGeobuf driver
"""

import hashlib
import math
from pathlib import Path

import geopandas as gpd
from shapely import STRtree  # noqa: F401  (ensures shapely 2.x)


DELIVERY_CRS = "EPSG:4326"
COORD_PRECISION = 7          # ~1 cm at the equator; pin for stable bytes
TILE_DEG = 1.0               # fixed 1-degree delivery grid
OUT_DIR = Path("tiles")


def normalize(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Reproject, round coordinates, and impose a stable column order."""
    gdf = gdf.to_crs(DELIVERY_CRS)
    gdf["geometry"] = gdf.geometry.set_precision(10 ** -COORD_PRECISION)
    cols = ["geometry"] + sorted(c for c in gdf.columns if c != "geometry")
    return gdf[cols]


def hilbert_key(gdf: gpd.GeoDataFrame, bits: int = 16) -> gpd.GeoDataFrame:
    """Sort features by a Hilbert index of their centroid for reproducibility."""
    c = gdf.geometry.centroid
    order = gpd.GeoSeries(c).hilbert_distance(level=bits)
    return gdf.assign(_h=order).sort_values("_h").drop(columns="_h")


def tile_index(geom) -> tuple[int, int]:
    """Return the fixed-grid (col, row) tile a geometry's centroid falls in."""
    c = geom.centroid
    return (math.floor(c.x / TILE_DEG), math.floor(c.y / TILE_DEG))


def sha256_of(path: Path, chunk: int = 1 << 20) -> str:
    """Content hash of an output file for the version manifest."""
    h = hashlib.sha256()
    with path.open("rb") as fh:
        for block in iter(lambda: fh.read(chunk), b""):
            h.update(block)
    return h.hexdigest()


def write_tiles(gdf: gpd.GeoDataFrame) -> dict[str, str]:
    """Write one deterministic FlatGeobuf per non-empty tile; return {name: hash}."""
    OUT_DIR.mkdir(exist_ok=True)
    gdf = hilbert_key(normalize(gdf))
    gdf = gdf.assign(_tile=gdf.geometry.map(tile_index))

    manifest: dict[str, str] = {}
    for (col, row), part in gdf.groupby("_tile"):
        part = part.drop(columns="_tile")
        name = f"tile_{col}_{row}.fgb"
        out = OUT_DIR / name
        # SPATIAL_INDEX=YES embeds the packed R-tree for range reads.
        part.to_file(out, driver="FlatGeobuf",
                     layer_options=["SPATIAL_INDEX=YES"])
        manifest[name] = sha256_of(out)
    return manifest


def main() -> None:
    gdf = gpd.read_file("input_features.gpkg")
    manifest = write_tiles(gdf)
    for name, digest in sorted(manifest.items()):
        print(f"{digest[:12]}  {name}")
    # Diff step: compare `manifest` to the committed one; upload changed tiles only.


if __name__ == "__main__":
    main()

The determinism comes from three pinned choices: a single CRS, a fixed coordinate precision via set_precision, and a reproducible feature order from the Hilbert sort. With those fixed, an unchanged tile serialises to identical bytes, so its sha256 is stable and the diff step only flags tiles whose features genuinely changed.

Validation & Output Verification

Confirm the index is embedded, the file is valid, and hashes are reproducible:

# Confirm the FlatGeobuf carries a spatial index and correct CRS
ogrinfo -so tiles/tile_5_51.fgb tile_5_51 | grep -E "FID Column|Extent|SRS"

# Re-running the writer on unchanged input must reproduce identical hashes
python tile_flatgeobuf.py > run1.txt
python tile_flatgeobuf.py > run2.txt
diff run1.txt run2.txt && echo "DETERMINISTIC: hashes stable"

Assert feature counts survive the round trip and that a bounding-box read returns a subset rather than the whole file:

import geopandas as gpd
tile = gpd.read_file("tiles/tile_5_51.fgb")
assert len(tile) > 0, "empty tile written"
# A bbox read should return fewer features than the full tile
subset = gpd.read_file("tiles/tile_5_51.fgb", bbox=(5.1, 51.1, 5.3, 51.3))
assert len(subset) <= len(tile)
print(f"tile features: {len(tile)}, bbox subset: {len(subset)}")

Failure Modes

  • Symptom: identical input produces different file hashes across runs. Root cause: unpinned coordinate precision or non-deterministic feature order. Fix: apply set_precision and sort by a Hilbert key before writing, as in the script.
  • Symptom: a web client downloads the entire file instead of a bounding-box subset. Root cause: the FlatGeobuf was written without SPATIAL_INDEX=YES, so no packed R-tree exists. Fix: re-export with the spatial index option enabled.
  • Symptom: a small attribute edit re-hashes every tile. Root cause: editing before normalizing, so attribute reordering or reprojection perturbs all outputs. Fix: normalize first, then edit and re-tile, so only affected tiles change.
  • Symptom: features vanish near tile edges. Root cause: centroid-based tile assignment drops large geometries whose centroid sits outside their extent. Fix: assign by bounding-box overlap and duplicate boundary-spanning features across the tiles they touch.

Back to Cloud-Native Spatial Formats for Versioned Pipelines