FlatGeobuf and GeoPackage Ingestion in DuckDB

FlatGeobuf (.fgb) and GeoPackage (.gpkg) are the two container formats you reach for when a dataset must stay a single portable file yet still answer bounding-box reads without a full scan, and DuckDB ingests both through the GDAL-backed st_read table function. This guide sits under the DuckDB Spatial architecture and fundamentals reference and covers the concrete ingestion workflow: how each container lays out geometry and its index, when to prefer one over the other or over columnar Parquet, how to select layers and push a spatial pre-filter into the driver’s own index, and how to validate that the scan is doing what you asked. The recurring nuance is that st_read is an opaque OGR feature reader, not DuckDB’s native parallel columnar scanner — so the levers that make a Parquet read fast are not the levers that make a FlatGeobuf or GeoPackage read fast, and confusing the two is the most common source of a slow ingest.

Runtime Configuration & Memory Guardrails

Ingestion through st_read streams features one batch at a time out of GDAL and materializes them into DuckDB GEOMETRY vectors. It does not parallelize a single layer scan the way read_parquet does, so throughput is bounded by the OGR driver, not by thread count. Configure the session for the materialization and any downstream transform rather than for the scan itself.

INSTALL spatial;          -- one-time fetch of the extension (bundles GDAL/OGR + GEOS)
LOAD spatial;             -- per-session; st_read, ST_Read_Meta and COPY ... FORMAT GDAL live here

-- Threads help the CTAS materialization, ST_MakeValid, and any reprojection that
-- follows the scan; the single-layer OGR read itself stays serial, so do not expect
-- more threads to speed the raw st_read.
SET threads = 8;

-- Cap the working set below total RAM. A GeoPackage attribute join or an ST_Union
-- after ingest can multiply the footprint well past the file size on disk.
SET memory_limit = '10GB';

-- Drop row-order guarantees so the downstream CREATE TABLE ... AS pipelines
-- in parallel. Feature order from the FlatGeobuf/GPKG file is not preserved anyway.
SET preserve_insertion_order = false;

-- Spill target for the materialization stage. Point at local NVMe; a network mount
-- turns a spilling reprojection into an I/O cliff.
SET temp_directory = '/var/lib/duckdb/spill';

Before running an ingest into a production table, confirm the following hold:

Trade-off: because the layer scan is serial, the fastest way to shorten a large ingest is to read fewer features, not to add cores. That is exactly what the format-level spatial index buys you, and pushing a bounding box into the driver — covered next — is the single highest-leverage move for both formats.

How the Three Containers Differ

FlatGeobuf and GeoPackage solve the “one portable file” problem in opposite ways, and understanding the layout tells you which read patterns are cheap. A FlatGeobuf file is a magic header, a schema-and-CRS header, an optional packed Hilbert R-tree, then the features themselves as length-delimited FlatBuffers. The features are written in Hilbert-curve order, so spatially near features are byte-near on disk; the packed R-tree lets a reader seek directly to the byte ranges that overlap a query box without decoding anything in between. It is append-friendly and streamable — a consumer can begin reading before the whole file arrives — but it is write-once: there is no in-place update.

A GeoPackage is an ordinary SQLite database. Each vector layer is a table whose geometry column holds GeoPackageBinary (a small envelope header wrapping standard WKB), and its spatial index is a companion SQLite R*Tree virtual table (rtree_<table>_<geom>). Because it is SQLite, one .gpkg can hold many feature layers, attribute-only tables, and even raster tiles, and it supports random reads and in-place edits. The cost is that random access needs a seekable, low-latency file handle — SQLite issues many small reads — which makes a large .gpkg a poor fit for high-latency object storage.

GeoParquet is the columnar contrast: geometry is a WKB column, CRS lives in the file’s geo key-value metadata as PROJJSON, and pruning comes from per-row-group bounding-box statistics rather than a tree. DuckDB reads GeoParquet through its native vectorized read_parquet scanner — parallel and column-projecting — whereas it reads FlatGeobuf and GeoPackage through the serial OGR path. That difference dominates the decision, and the head-to-head numbers against legacy formats are quantified in the GeoParquet vs Shapefile performance reference.

How GeoParquet, FlatGeobuf and GeoPackage compare across layout, index, streaming and best use A four-row by three-column matrix. GeoParquet: columnar WKB column, row-group bbox statistics, parallel column scan, best for analytics on object storage. FlatGeobuf: row FlatBuffers in one file, packed Hilbert R-tree, seekable stream, best for portable indexed reads. GeoPackage: SQLite tables, SQLite R-tree tables, random serial access, best for multi-layer editable GIS. ONE FILE, THREE STRATEGIES — PICK BY READ PATTERN GeoParquet FlatGeobuf GeoPackage Layout Index Streaming Best use Columnar (WKB column) Row FlatBuffers, one file SQLite tables Row-group bbox stats Packed Hilbert R-tree SQLite R-tree tables Parallel column scan Seekable stream Random access, serial Analytics on object storage Portable indexed reads Multi-layer editable GIS

The practical decision rule follows from the matrix. Reach for GeoParquet as the default analytical store, because only it gets DuckDB’s parallel columnar scan and column projection. Reach for FlatGeobuf when you need one portable file that still supports index-accelerated bounding-box reads over HTTP or object storage — its seekable layout tolerates high latency far better than SQLite’s many small reads. Reach for GeoPackage when the data is genuinely multi-layer, needs in-place editing, or must interoperate with desktop GIS tooling that expects it. Whichever you land on, the ingestion path into DuckDB is the same st_read call.

Primary Execution Patterns

The canonical FlatGeobuf ingest is a CREATE TABLE ... AS SELECT over st_read, materializing the OGR feature stream into native GEOMETRY. Everything else — layer selection, spatial pre-filtering, CRS handling — is expressed as options on that single call.

-- Baseline: stream every feature of a FlatGeobuf into a native GEOMETRY table.
CREATE OR REPLACE TABLE roads AS
SELECT * FROM st_read('data/roads.fgb');

-- Inspect layers, CRS, and feature counts WITHOUT reading geometry.
-- ST_Read_Meta parses only the header, so it is cheap even over the network.
SELECT * FROM ST_Read_Meta('data/roads.fgb');

The index pays off when you narrow the read to a region. FlatGeobuf’s packed R-tree and GeoPackage’s SQLite R*Tree can both skip features that do not overlap a query box — but only if you pass the box as the spatial_filter_box option, so st_read hands it to the OGR driver. A plain WHERE clause cannot do this: st_read is opaque to the optimizer, so any predicate outside the option list runs after every feature has already been decoded.

-- Index-accelerated read: the box is pushed into the FlatGeobuf packed R-tree,
-- so only overlapping features are ever decoded off disk.
-- spatial_filter_box takes a BOX_2D; ST_Extent of an envelope yields one.
CREATE OR REPLACE TABLE roads_bbox AS
SELECT *
FROM st_read(
        'data/roads.fgb',
        spatial_filter_box = ST_Extent(ST_MakeEnvelope(-2.30, 51.35, -2.20, 51.42))
     );
The st_read ingestion path: spatial_filter_box pushes a box into the container's own index Left to right: a FlatGeobuf or GeoPackage file feeds the OGR driver inside st_read; a spatial_filter_box is pushed into the format index so only overlapping features decode; survivors materialize as native GEOMETRY vectors into a DuckDB table. A WHERE clause instead of the option decodes everything first. .fgb / .gpkg header + index + features st_read GDAL / OGR driver spatial_filter_box ↓ format index packed R-tree / SQLite R-tree skip non-overlap DuckDB table native GEOMETRY 2,048-row vectors A WHERE clause runs AFTER decode — only spatial_filter_box reaches the format index.

GeoPackage ingestion is identical except that a multi-layer container forces you to name the layer. Skipping layer= reads the first layer OGR reports, which is rarely the one you meant.

-- List the layers first; a GeoPackage routinely holds several feature tables.
SELECT layers FROM ST_Read_Meta('data/city.gpkg');

-- Pin the exact layer and push a bounding box into the SQLite R-tree.
CREATE OR REPLACE TABLE parcels AS
SELECT *
FROM st_read(
        'data/city.gpkg',
        layer = 'parcels',
        spatial_filter_box = ST_Extent(ST_MakeEnvelope(-2.30, 51.35, -2.20, 51.42))
     );

-- keep_wkb returns geometry as WKB_BLOB instead of GEOMETRY, skipping the decode
-- when you only need to pass bytes straight through to another writer.
CREATE OR REPLACE TABLE parcels_raw AS
SELECT * FROM st_read('data/city.gpkg', layer = 'parcels', keep_wkb = true);

For source formats that carry attributes you do not need, project the columns you keep in the SELECT list so the materialization is narrower — the OGR read still touches every field, but the table you persist stays lean. When the same region is queried repeatedly, materialize once into a DuckDB table or convert the layer to GeoParquet so subsequent reads use the parallel scanner instead of paying the OGR cost each time. The converse direction — turning a legacy Shapefile into a FlatGeobuf with an embedded index — uses the GDAL COPY writer and is worked through in its own walkthrough.

Execution Plan Validation

Because st_read is a table function, its plan node is a single TABLE_FUNCTION (labelled ST_READ) — you will not see the parallel PARQUET_SCAN nodes a GeoParquet read produces. What EXPLAIN ANALYZE does tell you is whether the spatial_filter_box actually pruned, by reporting how many rows the function emitted.

-- Measure the pushed-down read: the row count out of the function should equal the
-- features overlapping the box, NOT the file's total feature count.
EXPLAIN ANALYZE
SELECT count(*)
FROM st_read(
        'data/roads.fgb',
        spatial_filter_box = ST_Extent(ST_MakeEnvelope(-2.30, 51.35, -2.20, 51.42))
     );

What to assert in the output:

  • Rows emitted by ST_READ. Compare this to the total feature count from ST_Read_Meta. If the box covers a tenth of the extent and the function still emits every feature, the filter did not reach the index — check that you passed spatial_filter_box (a BOX_2D) and not a WHERE predicate, and that the file actually carries an index.
  • No downstream FILTER doing the geographic work. A FILTER node holding your bounding-box logic above the ST_READ means the box was applied after decode — every feature was read off disk first. Move the box into the option.
  • Single-threaded function scan is expected. Unlike read_parquet, one st_read call does not fan out across threads. If a large ingest is CPU-idle and I/O-idle, the bottleneck is OGR feature decode, and the remedy is to read fewer features (spatial filter, or a narrower layer), not more threads.

A FlatGeobuf without an index, or a query box that covers the whole extent, will read the entire file regardless of the option — the pushdown skips features, it does not conjure selectivity. Confirm the file is indexed via ST_Read_Meta before assuming the box will help.

Performance Trade-offs

Quantify these against your own data and access pattern, but the qualitative ordering is stable:

Dimension GeoParquet FlatGeobuf GeoPackage
DuckDB read path Native read_parquet (parallel) Serial OGR st_read Serial OGR st_read
Column projection Yes — skips column chunks No — every field decoded No — every field decoded
Bbox pushdown Row-group stats + optional bbox column Packed R-tree via spatial_filter_box SQLite R-tree via spatial_filter_box
Object-storage reads Excellent (range reads on row groups) Good (seekable, few large reads) Poor (SQLite needs many small reads)
Random single-feature access Weak Weak Strong (SQLite lookup)
In-place edit / append No Append-only Full read/write
Multiple layers in one file No No Yes
Typical best fit Analytical warehouse, lakehouse Portable indexed exchange Editable multi-layer GIS

The highest-leverage decision is the read path. If a dataset is queried analytically more than a handful of times, the one-time cost of converting it to GeoParquet is almost always repaid by the parallel, projectable scan — the GeoParquet parsing reference covers that ingest, and the storage-and-scan comparison against row formats is in GeoParquet vs Shapefile performance. FlatGeobuf and GeoPackage earn their place when the file must stay a single self-describing artifact — for handoff, for editing, or for indexed reads over the network — rather than a partitioned analytical store.

Edge Cases & Anti-Patterns

Multi-layer GeoPackage read pinned to the wrong layer (silent wrong data). A .gpkg with several feature tables returns the first one when layer= is omitted, so the query succeeds while reading data you never intended.

-- Enumerate layers, then always pin the one you want.
SELECT UNNEST(layers).name AS layer_name FROM ST_Read_Meta('data/city.gpkg');

-- sequential_layer_scan forces a full driver scan to enumerate hidden/virtual layers
-- some GeoPackages expose; slower, but exhaustive when a layer name is missing.
SELECT * FROM st_read('data/city.gpkg', layer = 'transport', sequential_layer_scan = true);

Assuming CRS survives the read (silent unit mismatch). The container records its CRS — GeoPackage in gpkg_spatial_ref_sys, FlatGeobuf in the header — but DuckDB GEOMETRY has no SRID slot, so after st_read the coordinates are just numbers. Read the source CRS from metadata, record it, and reproject before joining against any other layer.

-- Recover the source CRS from the container header before trusting coordinates.
SELECT UNNEST(layers).name AS layer, UNNEST(layers).geometry_fields AS gf
FROM ST_Read_Meta('data/city.gpkg');
-- Then reproject explicitly; mixing a degree layer with a metric one returns wrong,
-- not empty, results. See the CRS reference for ST_Transform semantics.

Feeding a mixed-frame set into a spatial join produces meaningless overlaps rather than an error; the guardrails and ST_Transform mechanics are in how DuckDB Spatial handles coordinate systems.

Large GeoPackage over a network mount (latency cliff). SQLite issues many small random reads to walk its R*Tree and page the geometry table. Over /vsicurl, S3, or an NFS mount, each read is a round trip and the ingest slows by orders of magnitude. For remote data, prefer FlatGeobuf (few large seekable reads) or convert to partitioned GeoParquet; reserve GeoPackage for local disk. The resident-versus-spilled behaviour of the materialization once bytes arrive is governed by the in-memory vs disk storage boundaries.

Expecting a WHERE box to accelerate the read (no pushdown). As shown above, only spatial_filter_box reaches the index. WHERE ST_Intersects(geom, ...) after st_read decodes the whole file first — correct answer, full-scan cost.

Reading a FlatGeobuf that was written without a spatial index. GDAL can write an unindexed .fgb; spatial_filter_box then falls back to a linear feature scan. Verify the index exists before relying on selectivity, and when you control the write, keep the index on so downstream indexed reads and the R-tree index internals apply.

Ingest Verification & Regression Analysis

An ingest is only trustworthy if the feature count and CRS survive intact and the spatial pre-filter keeps pruning as data versions change. Capture both the metadata-reported count and the materialized count, assert they agree (for an unfiltered load), and confirm a boxed read stays selective. This harness slots into the orchestration covered in the Python and DuckDB integration workflows reference.

import duckdb

con = duckdb.connect(":memory:")
con.execute("INSTALL spatial; LOAD spatial;")
con.execute("SET threads = 8; SET memory_limit = '10GB';")

SRC = "data/roads.fgb"

# 1. Header-only metadata: cheap, reads no geometry.
meta = con.execute(
    "SELECT UNNEST(layers).feature_count AS fc FROM ST_Read_Meta(?)", [SRC]
).fetchone()
declared = meta[0]

# 2. Full materialization count must match the header's declared feature count.
loaded = con.execute(f"SELECT count(*) FROM st_read('{SRC}')").fetchone()[0]
assert loaded == declared, f"Feature-count drift: header={declared}, loaded={loaded}"

# 3. A boxed read must emit far fewer rows than the full file, proving the
#    packed R-tree pruned. A ratio near 1.0 means the pushdown was lost.
boxed = con.execute(f"""
    SELECT count(*) FROM st_read(
        '{SRC}',
        spatial_filter_box = ST_Extent(ST_MakeEnvelope(-2.30, 51.35, -2.20, 51.42))
    )
""").fetchone()[0]
assert boxed < declared * 0.5, f"spatial_filter_box did not prune: {boxed}/{declared}"
print(f"OK: {loaded} features, boxed read pruned to {boxed}")

Three signals are worth tracking across ingests: the header-versus-materialized count (a mismatch means the driver skipped or duplicated features, often an encoding or invalid-geometry issue), the boxed-read selectivity ratio (a drift toward 1.0 means the index was dropped on a re-export), and the source CRS string from ST_Read_Meta (a change means an upstream re-projection you must mirror downstream). For the GeoPandas handoff that consumes these validated tables without a WKT round-trip, see DuckDB-to-GeoPandas sync.

See also

Up: DuckDB Spatial Architecture & Fundamentals