DuckDB Spatial Analytical SQL for GIS

Reading GeoPackage Layers and Attribute Tables

A GeoPackage is an SQLite database wearing a GIS format, and reading one well means treating it as a database — enumerate the catalogue, name the layer, and expect a schema — rather than as a file you point a reader at. This walkthrough, part of the FlatGeobuf and GeoPackage ingestion reference, covers listing what a package actually contains, reading a named layer and its attribute tables, and the two failures specific to a container whose spatial index is maintained by triggers.

Root-Cause Analysis: what goes wrong reading a GeoPackage

  • The wrong layer is read. A package can hold many layers, and a read that does not name one takes the first — which is correct by accident and changes when the producer reorders them.
  • Attribute-only tables are invisible. A GeoPackage may contain tables with no geometry at all: lookup tables, metadata, relationship tables. A spatial reader does not surface them, and the attributes the analysis needs are in one of them.
  • The spatial index disagrees with the data. The R-tree lives in a virtual table kept in step by triggers. A package written by a tool that bypassed those triggers has an index that misses rows, and the symptom is a bounding-box query returning too few results rather than an error.
  • Column projection does not happen. GeoPackage is row-major, so selecting three of forty columns still reads all forty off the disk. The projection happens after the read rather than instead of it.
  • The declared CRS is dropped on import. The package records its frame in gpkg_spatial_ref_sys, and DuckDB has nowhere to keep it. Read it deliberately or it is gone.

The distinguishing question is whether you know what the package contains. Almost every failure above follows from reading a layer without having listed the catalogue first.

What is inside a GeoPackage A catalogue table, a CRS table, feature layers, attribute-only tables, and trigger-maintained R-tree virtual tables. COMPONENT WHAT IT HOLDS WHY IT MATTERS gpkg_contents every layer, type and extent read it first, always gpkg_spatial_ref_sys the CRS definitions DuckDB will not carry them one table per feature layer attributes plus a geom blob this is what a reader shows you attribute-only tables lookups and relationships invisible to a spatial reader rtree_<layer>_geom the spatial index maintained by triggers, can drift

A spatial reader shows you the third row. The other four are where the surprises live.

Deterministic Configuration 

INSTALL spatial; LOAD spatial;
INSTALL sqlite;  LOAD sqlite;      -- for the attribute-only tables

SET memory_limit = '4GB';
SET threads = 8;

Optimized Execution Pattern

The pattern is two reads rather than one: the catalogue first, then the named layer — with the SQLite reader used for anything the spatial reader does not surface.

-- ANTI-PATTERN: takes whichever layer happens to be first, and silently
-- changes meaning when the producer adds one.
CREATE TABLE features AS SELECT * FROM st_read('survey.gpkg');
-- PATTERN: list, then name. The catalogue is one query and it removes the
-- entire class of "which layer did I get" question.
SELECT table_name, data_type, srs_id, min_x, min_y, max_x, max_y
FROM sqlite_scan('survey.gpkg', 'gpkg_contents');

CREATE OR REPLACE TABLE parcels AS
SELECT * FROM st_read('survey.gpkg', layer = 'parcels_2024');

-- And the attribute-only tables a spatial reader will not show you.
CREATE OR REPLACE TABLE land_use_lookup AS
SELECT * FROM sqlite_scan('survey.gpkg', 'land_use_codes');

sqlite_scan is the part most readers of GeoPackage documentation never reach, and it is what makes the format usable as what it is — a database with a schema, rather than a bag of features with a lookup table stranded inside it.

Two readers, two purposes The spatial reader parses geometry and uses the package index; the SQLite reader reaches every table including the non-spatial ones. READER GIVES YOU USE IT FOR st_read native GEOMETRY, index-backed bbox the feature layers sqlite_scan every table, geometry as a blob catalogue, lookups, metadata

Using only the first leaves half the package unreachable.

Reading the CRS before it disappears

A GeoPackage records its coordinate reference system properly, in gpkg_spatial_ref_sys, keyed by the srs_id each layer declares in the catalogue. DuckDB has nowhere to put that, so the information is available at exactly one moment: while you are reading the package. Capturing it then is two lines; recovering it afterwards means going back to the source.

-- The layer's frame, joined out of the catalogue. Record this alongside the
-- table you create, because nothing downstream will carry it for you.
SELECT c.table_name, c.srs_id, s.organization, s.organization_coordsys_id, s.definition
FROM sqlite_scan('survey.gpkg', 'gpkg_contents')       c
JOIN sqlite_scan('survey.gpkg', 'gpkg_spatial_ref_sys') s USING (srs_id)
WHERE c.table_name = 'parcels_2024';

The organization and organization_coordsys_id columns together give the EPSG code in the usual case, and the definition column gives the full WKT for the cases where they do not.

Diagnostic Queries & Plan Validation

The failure specific to this format is an index that disagrees with the data, and it is checkable in one query because both are ordinary tables.

-- Feature count against index count. A gap means the R-tree was not
-- maintained, so a bounding-box read will silently miss rows.
SELECT
    (SELECT count(*) FROM sqlite_scan('survey.gpkg', 'parcels_2024'))          AS features,
    (SELECT count(*) FROM sqlite_scan('survey.gpkg', 'rtree_parcels_2024_geom')) AS indexed;

If the two disagree, do not use the package index: read the layer in full and build a DuckDB R-tree over it instead. The drift cannot be repaired from outside the producing tool, and a partially populated index is worse than none because it returns a plausible subset.

Three symptoms and what each points at Too few rows means a drifted index; slow reads mean row-major storage; missing attributes mean a table the spatial reader did not show. SYMPTOM POINTS AT FIX bbox query returns too few rows a drifted trigger-maintained index read in full; build your own slower than the same data as Parquet row-major storage, no projection convert once attributes appear to be missing an attribute-only table list the catalogue

Only the first is a defect. The second is the format and the third is a habit.

Converting on arrival

GeoPackage is an excellent interchange format and a poor analytical one, for the same reason in both cases: it is a row-major database with an editable schema. That makes it ideal for handing a colleague a file they can open in desktop GIS and edit, and it means every analytical query reads every column of every row it touches.

The habit worth forming is to convert on arrival — read the catalogue, read each layer you need, join in the lookups, and write GeoParquet. That takes one pass, preserves the attributes and the CRS you deliberately captured, and turns every subsequent query into a partial read. Keeping the original package as the archived source is sensible; querying it repeatedly is not.

Geometry Validation & Fallback Routing

Where the package has to be read repeatedly, at least avoid re-reading the parts that do not change.

-- Convert once, with the lookups joined in and the frame recorded in the
-- table name, so the analytical copy is self-describing.
COPY (
    SELECT p.*, l.description AS land_use_description
    FROM st_read('survey.gpkg', layer = 'parcels_2024') p
    LEFT JOIN sqlite_scan('survey.gpkg', 'land_use_codes') l
           ON p.land_use_code = l.code
) TO 'parcels_2024_epsg27700.parquet' (FORMAT PARQUET, ROW_GROUP_SIZE 30000);

Why a GeoPackage arrives at all

It is worth understanding why this format keeps turning up, because the answer explains what to do with it. GeoPackage is the format desktop GIS produces when someone is asked for “the data”: it holds several layers in one file, it survives being emailed, it opens in every tool your colleagues use, and it can be edited in place. Every one of those properties is about handing data between people, and none of them is about querying it.

That is the whole disposition. Treat a GeoPackage as an arrival: list it, read what you need, capture the frame, join in the lookups, and write the analytical copy. What comes back out at the end of a project — the result someone else has to open — is very reasonably a GeoPackage again, written from the analytical copy rather than from the original.

Frequently Asked Questions

How do I see what layers a GeoPackage contains?

Read gpkg_contents with sqlite_scan. It names every layer, its type and its extent, and it is one query. Reading a package without listing it first is how a pipeline ends up processing whichever layer the producer happened to write first.

Why can I not see one of the tables?

Because it has no geometry, and a spatial reader only surfaces feature layers. Lookup tables, metadata and relationship tables are ordinary SQLite tables and need sqlite_scan. This is the most common reason an attribute “is not in the file” when it demonstrably is.

My bounding-box query returns too few rows. Why?

Almost certainly a drifted spatial index. The R-tree is a virtual table kept in step by triggers, so a package written by a tool that bypassed them has an index missing rows — and a bbox read served by that index returns a plausible subset with no error. Compare the feature count against the index count.

Does DuckDB use the package’s own index?

Through GDAL, a bounding-box restriction can be pushed down to it, which is what makes a spatial window over a GeoPackage fast. What cannot be pushed down is an attribute filter or a column projection, because the format is row-major.

Should I keep querying the GeoPackage or convert it?

Convert, if it will be queried more than a few times. GeoPackage is an interchange format: editable, portable, and read in full on every query. One conversion pass turns every later query into a partial read, and the original stays as the archived source.

What happens to the CRS?

It is in the package and it does not survive the import, because DuckDB geometry carries no frame. Join gpkg_contents to gpkg_spatial_ref_sys while you are reading and record the result, because that is the only moment the information is in front of you.

Up: FlatGeobuf and GeoPackage Ingestion in DuckDB