Introduction

GeoPackage is an OGC standard for storing geospatial data in a single SQLite database file. A GeoPackage can hold vector features, attribute tables and raster tiles, together with the spatial reference systems and spatial indexes that go with them. QGIS, GDAL and most other desktop GIS software can read it, and because it is just a file, there is no server to run.

This book shows you how to read and write GeoPackage files from Rust with the geopackage crate, and how to use the gpkg command-line tool built on it.

The crates

The project is split into four crates:

• geopackage is the library most people want. It creates and opens files, reads and writes feature and attribute layers, manages spatial indexes, handles tile pyramids and tiled gridded coverages, and can move data in and out as Apache Arrow record batches.
• geopackage-core contains the parts of the format that don’t need a database, such as the geometry encoding and the table definitions. You only need it directly if you’re working with the format without SQLite.
• geopackage-cli provides gpkg, a command-line tool for inspecting, validating, indexing, repairing and copying files, and for reading tiles and coverages.
• geopackage-ffi is a C library over the same code, which passes data through the Arrow C Data Interface.

Installation

To use the library in a Rust project:

$ cargo add geopackage

By default, the library links against your system’s SQLite, which must include the RTree module (SQLITE_ENABLE_RTREE). The SQLite that ships with macOS and with Linux distributions has it. On Linux you also need the development package, such as libsqlite3-dev on Debian and Ubuntu.

If you’d rather not depend on the system SQLite, or you’re on Windows, which doesn’t have one, enable the bundled feature. It compiles its own copy of SQLite, so you’ll need a C compiler:

$ cargo add geopackage --features bundled

To install the command-line tool:

$ cargo install geopackage-cli

The tool uses the bundled SQLite by default. To link the system library instead, pass --no-default-features.

All the crates need Rust 1.95 or newer.

How this book is organised

The book has four parts:

• The tutorial walks you through creating and querying your first GeoPackage. Start there if you’re new to the library.
• The how-to guides show you how to do specific tasks, such as adding a spatial index or copying features between files.
• The explanation chapters describe how parts of the format work, and why the library behaves as it does.
• The reference page links to the API documentation on docs.rs, the C header, and a summary of the gpkg commands.

Still to come

These chapters haven’t been written yet:

• Migration notes for people coming from gdal, gpkg-rs or rusqlite-gpkg, showing the equivalent calls and where the models differ.
• A guide to the C library: building and linking geopackage-ffi, its conventions for handles and errors, and passing data through the Arrow C Data Interface. Until then, the five example programs listed on the reference page are the best place to start.
• More how-to guides, covering tile pyramids, attribute tables, metadata and the Related Tables Extension.

Your first GeoPackage

In this tutorial we start from an empty directory and build a GeoPackage file containing four cities. We query it by bounding box from Rust, then inspect it from the shell with the gpkg command-line tool.

You need a Rust toolchain, version 1.95 or newer, and an SQLite to link against. On Linux (install libsqlite3-dev or your distribution’s equivalent) and macOS that is the system library. On Windows, use the vendored copy via the bundled feature (see the introduction).

Create the project

First, we make a new binary project and change into it:

$ cargo new cities
$ cd cities

Now we add two dependencies. geopackage is the library, and geo-types provides our point geometries:

$ cargo add geopackage geo-types

The output ends with a summary of what was added:

      Adding geopackage v0.10.2 to dependencies
      Adding geo-types v0.7.20 to dependencies

The first build compiles every dependency, so it takes a minute or so. Later builds are quick.

Create the file

Open src/main.rs and replace its contents with this:

use geopackage::GeoPackage;

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let gpkg = GeoPackage::create("cities.gpkg")?;
    println!("created cities.gpkg");

    gpkg.close()?;
    Ok(())
}

GeoPackage::create fails rather than overwrite an existing file, so we delete the previous file before each run. That way every run starts from nothing:

$ rm -f cities.gpkg && cargo run

The output should end with:

created cities.gpkg

Let’s check the directory:

$ ls
Cargo.lock  Cargo.toml  cities.gpkg  src  target

cities.gpkg is a single file, with no -wal or -shm files beside it. This is the library’s default, so you can hand the file to someone else as it is.

Define a layer

The file has no layers yet, so it contains nothing. Next we declare a table called cities with a text column, an integer column and a point geometry column in EPSG:4326.

Add the two new use lines and the create_layer call, so that main.rs reads:

use geopackage::core::types::{ColumnType, GeometryType};
use geopackage::{ColumnSpec, GeoPackage, GeometrySpec, TableSchemaBuilder};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let gpkg = GeoPackage::create("cities.gpkg")?;
    println!("created cities.gpkg");

    gpkg.create_layer(
        &TableSchemaBuilder::new("cities")
            .column(ColumnSpec::new("name", ColumnType::Text(None)))
            .column(ColumnSpec::new("population", ColumnType::Integer))
            .geometry(GeometrySpec::new(GeometryType::Point, 4326)),
    )?;
    println!("created layer \"cities\"");

    gpkg.close()?;
    Ok(())
}

Run it again:

$ rm -f cities.gpkg && cargo run
created cities.gpkg
created layer "cities"

create_layer also built an empty RTree spatial index on the geometry column, which we will use shortly. The spatial index explains why the index is built now, while empty, rather than after the rows arrive.

Write four cities

Now we add some data. write_all takes an iterator of features and a batch size. A batch size of 0 writes all of them in one transaction.

Add the geo_types::Point import, extend the geopackage import, and append the write:

use geo_types::Point;
use geopackage::core::types::{ColumnType, GeometryType};
use geopackage::{ColumnSpec, GeoPackage, GeometrySpec, NewFeature, TableSchemaBuilder, Value};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let gpkg = GeoPackage::create("cities.gpkg")?;
    println!("created cities.gpkg");

    gpkg.create_layer(
        &TableSchemaBuilder::new("cities")
            .column(ColumnSpec::new("name", ColumnType::Text(None)))
            .column(ColumnSpec::new("population", ColumnType::Integer))
            .geometry(GeometrySpec::new(GeometryType::Point, 4326)),
    )?;
    println!("created layer \"cities\"");

    let layer = gpkg.layer("cities")?;
    layer.write_all(
        vec![
            NewFeature::new(
                Point::new(-6.26, 53.35),
                vec![Value::Text("Dublin".into()), Value::Integer(592_713)],
            ),
            NewFeature::new(
                Point::new(-0.13, 51.51),
                vec![Value::Text("London".into()), Value::Integer(8_866_180)],
            ),
            NewFeature::new(
                Point::new(2.35, 48.86),
                vec![Value::Text("Paris".into()), Value::Integer(2_048_472)],
            ),
            NewFeature::new(
                Point::new(13.40, 52.52),
                vec![Value::Text("Berlin".into()), Value::Integer(3_662_381)],
            ),
        ],
        0,
    )?;
    println!("wrote {} features", layer.features()?.len());

    gpkg.close()?;
    Ok(())
}

Run it:

$ rm -f cities.gpkg && cargo run
created cities.gpkg
created layer "cities"
wrote 4 features

The values in each NewFeature follow the order in which the columns were declared. The geometry and the feature id are not in that list. Each is handled separately.

Query by bounding box

Rather than read every city back, let’s ask a spatial question. features_in returns the rows whose geometry intersects a box given in the layer’s own coordinates.

Add BoundingBox to the import list, and add the loop after the write:

use geo_types::Point;
use geopackage::core::types::{ColumnType, GeometryType};
use geopackage::{
    BoundingBox, ColumnSpec, GeoPackage, GeometrySpec, NewFeature, TableSchemaBuilder, Value,
};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let gpkg = GeoPackage::create("cities.gpkg")?;
    println!("created cities.gpkg");

    gpkg.create_layer(
        &TableSchemaBuilder::new("cities")
            .column(ColumnSpec::new("name", ColumnType::Text(None)))
            .column(ColumnSpec::new("population", ColumnType::Integer))
            .geometry(GeometrySpec::new(GeometryType::Point, 4326)),
    )?;
    println!("created layer \"cities\"");

    let layer = gpkg.layer("cities")?;
    layer.write_all(
        vec![
            NewFeature::new(
                Point::new(-6.26, 53.35),
                vec![Value::Text("Dublin".into()), Value::Integer(592_713)],
            ),
            NewFeature::new(
                Point::new(-0.13, 51.51),
                vec![Value::Text("London".into()), Value::Integer(8_866_180)],
            ),
            NewFeature::new(
                Point::new(2.35, 48.86),
                vec![Value::Text("Paris".into()), Value::Integer(2_048_472)],
            ),
            NewFeature::new(
                Point::new(13.40, 52.52),
                vec![Value::Text("Berlin".into()), Value::Integer(3_662_381)],
            ),
        ],
        0,
    )?;
    println!("wrote {} features", layer.features()?.len());

    println!("features in the box:");
    for feature in layer.features_in(BoundingBox::new(-7.0, 53.0, -6.0, 54.0))? {
        let feature = feature?;
        let name = feature.value("name").and_then(|v| v.as_str()).unwrap_or("");
        let population = feature
            .value("population")
            .and_then(|v| v.as_i64())
            .unwrap_or(0);
        println!("  {name}: {population}");
    }

    gpkg.close()?;
    Ok(())
}

Run it again:

$ rm -f cities.gpkg && cargo run
created cities.gpkg
created layer "cities"
wrote 4 features
features in the box:
  Dublin: 592713

Only Dublin came back, although we wrote four cities. The box spans one degree of longitude west of Greenwich and one degree of latitude north of 53, and Dublin is the only one of our cities inside it. The spatial index that create_layer built answered the query. A layer without an index returns exactly the same rows from a full scan.

Read the geometries back

So far we have read back the name and population, but not the point itself. A feature returns its geometry lazily. feature.geometry() parses the stored blob into a wrapper over its bytes, and to_geo on that wrapper converts it to an owned geo_types::Geometry. This mirrors the write side, which accepted our geo_types::Point values directly. to_geo returns an Option because a GeoPackage can hold curve types that geo-types cannot represent. Our points are not among them.

Extend the loop body to fetch the location:

use geo_types::Point;
use geopackage::core::types::{ColumnType, GeometryType};
use geopackage::{
    BoundingBox, ColumnSpec, GeoPackage, GeometrySpec, NewFeature, TableSchemaBuilder, Value,
};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let gpkg = GeoPackage::create("cities.gpkg")?;
    println!("created cities.gpkg");

    gpkg.create_layer(
        &TableSchemaBuilder::new("cities")
            .column(ColumnSpec::new("name", ColumnType::Text(None)))
            .column(ColumnSpec::new("population", ColumnType::Integer))
            .geometry(GeometrySpec::new(GeometryType::Point, 4326)),
    )?;
    println!("created layer \"cities\"");

    let layer = gpkg.layer("cities")?;
    layer.write_all(
        vec![
            NewFeature::new(
                Point::new(-6.26, 53.35),
                vec![Value::Text("Dublin".into()), Value::Integer(592_713)],
            ),
            NewFeature::new(
                Point::new(-0.13, 51.51),
                vec![Value::Text("London".into()), Value::Integer(8_866_180)],
            ),
            NewFeature::new(
                Point::new(2.35, 48.86),
                vec![Value::Text("Paris".into()), Value::Integer(2_048_472)],
            ),
            NewFeature::new(
                Point::new(13.40, 52.52),
                vec![Value::Text("Berlin".into()), Value::Integer(3_662_381)],
            ),
        ],
        0,
    )?;
    println!("wrote {} features", layer.features()?.len());

    println!("features in the box:");
    for feature in layer.features_in(BoundingBox::new(-7.0, 53.0, -6.0, 54.0))? {
        let feature = feature?;
        let name = feature.value("name").and_then(|v| v.as_str()).unwrap_or("");
        let population = feature
            .value("population")
            .and_then(|v| v.as_i64())
            .unwrap_or(0);
        let location = match feature.geometry()?.and_then(|g| g.to_geo()) {
            Some(geo_types::Geometry::Point(p)) => format!("{}, {}", p.x(), p.y()),
            _ => "no location".into(),
        };
        println!("  {name}: {population} at {location}");
    }

    gpkg.close()?;
    Ok(())
}

Run it again:

$ rm -f cities.gpkg && cargo run
created cities.gpkg
created layer "cities"
wrote 4 features
features in the box:
  Dublin: 592713 at -6.26, 53.35

The point comes back exactly as we wrote it. You don’t have to convert, though. An algorithm written against the geo-traits interfaces can read coordinates directly from the wrapper, without ever creating an owned geometry. The geometry chapter covers this.

Look at the file from outside

So far only our own program has reported on the file. Let’s check it with a separate tool. Install the command-line tool:

$ cargo install geopackage-cli

Then ask it what our file contains:

$ gpkg info cities.gpkg
cities.gpkg
  version: 1.4

layer "cities" (features)
  rows:     4
  geometry: geom (POINT, srs_id 4326)
  srs:      WGS 84 geodetic (EPSG:4326)
  index:    current
  columns:
    fid                  INTEGER PRIMARY KEY
    name                 TEXT
    population           INTEGER
    geom                 POINT

extensions:
  gpkg_rtree_index                   cities.geom            implemented [write-only]

rows: 4 matches what our program reported, and index: current is the spatial index we never named: create_layer built it, and it uses the GeoPackage 1.4 trigger set. We never declared the fid and geom columns either. fid is the default primary key and geom the default geometry column name. You can change both with TableSchemaBuilder.

Let’s also check that the file is well formed:

$ gpkg validate cities.gpkg
cities.gpkg
  no findings (this crate's checks; the OGC ETS is the authority)

What we have built

We wrote a program that creates a GeoPackage, declares an indexed point layer, writes features into it, answers a spatial query and reads the geometries back as geo-types values. We then confirmed the result with a separate tool. cities.gpkg is an ordinary GeoPackage, so QGIS, GDAL and anything else that reads the format will open it.

From here:

• The how-to guides cover specific tasks: indexing an existing layer, repairing a file, copying features between files, reading a layer as Arrow batches, and validating a file.
• The explanation chapters cover why the extent, the transactions, the spatial index, the geometry encoding and tiled gridded coverages work the way they do.
• The reference page points at the full API documentation.

How to add a spatial index to an existing layer

A feature layer without an RTree index answers bounding-box queries by a full scan. This guide shows you how to build the index so those queries can use it.

From the shell

If the file is on disk and you don’t need to do anything else to it:

$ gpkg index roads.gpkg roads
roads: index built over 1204553 rows

gpkg index only builds an index where none exists. If the layer already has one, the command reports what it found and exits non-zero:

$ gpkg index roads.gpkg roads
gpkg: roads already has an index (legacy trigger set); run `gpkg repair`

In that case, see How to repair a file’s spatial indexes.

From Rust

Open the file read-write and take the layer handle:

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::GeoPackage;

let gpkg = GeoPackage::open("roads.gpkg")?;
let layer = gpkg.layer("roads")?;
Ok(()) }

Check the status first. create_spatial_index fails with Error::SpatialIndexExists if the rtree_<table>_<column> virtual table already exists, whatever state it is in:

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::{GeoPackage, SpatialIndexStatus};
let gpkg = GeoPackage::open("roads.gpkg")?;
let layer = gpkg.layer("roads")?;
match layer.spatial_index_status()? {
    SpatialIndexStatus::Absent => layer.create_spatial_index()?,
    SpatialIndexStatus::Current => {}
    // Present but not usable: repair it instead of building over it.
    SpatialIndexStatus::Legacy | SpatialIndexStatus::Stale => {
        layer.repair_spatial_index()?
    }
    _ => {}
}
Ok(()) }

create_spatial_index creates the virtual table, installs the GeoPackage 1.4 trigger set, populates the index from the existing rows, and registers gpkg_rtree_index in gpkg_extensions, all in one transaction. Like the triggers, it skips NULL and empty geometries.

Choose the build path explicitly

The index is populated by one of two paths, chosen by comparing the row count with BulkIndexOptions::bulk_threshold (10,000 rows by default). To choose the path yourself, use create_spatial_index_with:

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::{BulkIndexOptions, BulkVerification};

use geopackage::GeoPackage;
let gpkg = GeoPackage::open("roads.gpkg")?;
let layer = gpkg.layer("roads")?;
// Force the bulk build, and verify the result.
layer.create_spatial_index_with(
    BulkIndexOptions::always_bulk().with_verification(BulkVerification::Structure),
)?;
Ok(()) }

Set a verification level if a bad build must be caught at build time rather than found later:

• BulkVerification::None is the default and checks nothing.
• BulkVerification::Contents reads every entry back and checks it against the envelopes accumulated while writing.
• BulkVerification::Structure adds rtreecheck over the tree.
• BulkVerification::Database adds a whole-database PRAGMA integrity_check, which is O(database) rather than O(index).

Every level above None enables an automatic fallback to the per-row triggered build. A failed check makes the build slower, not wrong.

If you will append many rows to the layer after the build, lower the node fill so that the first append into a node doesn’t split it immediately:

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::{BulkIndexOptions, GeoPackage};
let gpkg = GeoPackage::open("roads.gpkg")?;
let layer = gpkg.layer("roads")?;
layer.create_spatial_index_with(BulkIndexOptions::default().with_fill_factor(0.7))?;
Ok(()) }

Index a layer as it is written

If the rows are not in the file yet, don’t build the index afterwards. Leave TableSchemaBuilder::spatial_index at its default of true, which creates an empty index, and let write_all fill it. The bulk path reuses the envelopes computed while encoding the geometries, which is between 20% and 47% cheaper than a separate build.

Then

• Confirm the result with gpkg info <file>, whose index: line reports the status, or with Layer::spatial_index_status.
• The spatial index: structure, contents, and the bulk build explains why checking an index’s structure and checking its contents differ in cost, and what the bulk build does.
• API reference: Layer::create_spatial_index, Layer::create_spatial_index_with, BulkIndexOptions, BulkVerification.

How to repair a file’s spatial indexes

An index in a file you receive may be maintained by a pre-1.4 trigger set, or may have lost its triggers entirely. This library trusts neither for queries, so bounding-box reads fall back to a full scan until the index is repaired. This guide shows you how to repair it.

Repair is never automatic. The library rewrites an existing trigger set only when you ask it to.

From the shell

To repair every layer in a file that needs it:

$ gpkg repair places.gpkg
points: legacy trigger set index repaired to the 1.4 trigger set

To repair one named layer, add its name:

$ gpkg repair places.gpkg points

If nothing needs repairing, the command says so and exits zero:

$ gpkg repair places.gpkg
nothing to repair

gpkg repair leaves a layer with no index alone, because an unindexed layer reads correctly. If you want an index built, see How to add a spatial index to an existing layer.

From Rust

Open the file read-write, preferably with the lenient open. A file that needs repair already has something wrong with it, and a strict open can fail on a legacy application_id before you reach the index.

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::{GeoPackage, SpatialIndexStatus};

let gpkg = GeoPackage::open_lenient("places.gpkg")?;

for layer in gpkg.layers()? {
    if layer.geometry_column().is_none() {
        continue;
    }
    match layer.spatial_index_status()? {
        SpatialIndexStatus::Legacy | SpatialIndexStatus::Stale => {
            layer.repair_spatial_index()?;
        }
        // Current needs nothing; Absent is a choice, not a defect.
        _ => {}
    }
}
Ok(()) }

repair_spatial_index drops every RTree trigger on the table, installs the 1.4 set, and rebuilds the index content, in one transaction. It returns without doing anything when the index is already current, and fails with Error::NoSpatialIndex when there is no index to repair.

If the structure is right but the contents are in doubt

spatial_index_status checks only the structure. It cannot tell whether the entries agree with the geometries. An index can be Current and still be wrong if rows were written while the triggers were absent, or if another tool populated it incompletely. To check the contents, audit the index:

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::GeoPackage;
let gpkg = GeoPackage::open_lenient("places.gpkg")?;
let layer = gpkg.layer("points")?;
let audit = layer.audit_spatial_index()?;
if !audit.is_consistent() {
    eprintln!(
        "{} missing, {} extra, {} not covering",
        audit.missing, audit.extra, audit.not_covering
    );
    layer.rebuild_spatial_index()?;
}
Ok(()) }

The audit reads every geometry in the layer, so run it as a deliberate check, not before each query. It writes nothing. To fix what it finds, call rebuild_spatial_index. It does unconditionally what the repair does only when the structure is wrong.

gpkg validate runs the same audit across the whole file. This is the cheaper way to find out whether any layer needs rebuilding. See How to validate a GeoPackage and act on the findings.

Then

• The spatial index: structure, contents, and the bulk build explains why repair is never automatic, and why a legacy trigger set is worth replacing.
• API reference: Layer::spatial_index_status, Layer::repair_spatial_index, Layer::audit_spatial_index, Layer::rebuild_spatial_index, SpatialIndexAudit.

How to copy features between files without decoding geometry

If you read a geometry into a geo-types value and write it back, every coordinate is parsed for nothing, and the non-linear curve types are lost entirely, because geo-traits has no representation for an arc. This guide copies rows between files as WKB bytes instead. Nothing is re-encoded, and a CIRCULARSTRING arrives unchanged.

From the shell

For a whole file, gpkg copy already works this way:

$ gpkg copy places.gpkg places-copy.gpkg

It copies feature and attribute layers, with their schemas, spatial reference systems and rows, and a spatial index wherever the source had one. It does not copy tile pyramids, tiled gridded coverages or the extension tables, and it names whatever it left behind at the end.

Mirror the source layer

In Rust, start with the source file open and the destination file created. The destination’s geometry column must match the source in three respects: the geometry type, the z and m flags, and the SRS. These are easy to leave at their defaults, and the copy then fails on the first row.

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::core::types::ColumnType;
use geopackage::{ColumnSpec, GeoPackage, GeometrySpec, TableSchemaBuilder};

let src = GeoPackage::open_read_only("places.gpkg")?;
let dst = GeoPackage::create("subset.gpkg")?;

let source = src.layer("points")?;
let geom = source.geometry_column().expect("a feature layer");

// The SRS must exist in the destination before a layer can refer to it.
// `add_srs` does nothing if it is already there.
if let Some(srs) = src.srs(geom.srs_id)? {
    dst.add_srs(&srs)?;
}

let mut builder = TableSchemaBuilder::new("points").geometry(
    GeometrySpec::new(geom.geometry_type, geom.srs_id)
        .column_name(geom.column_name.clone())
        .z(geom.z)
        .m(geom.m),
);
// Copy every column except the geometry and the primary key. Those two are
// declared elsewhere, and `Feature::values` skips them.
for column in &source.schema().columns {
    if column.is_primary_key() || column.name == geom.column_name {
        continue;
    }
    // A declared type outside the spec vocabulary parses to `None`. BLOB is
    // the typeless case under SQLite's affinity rules.
    let column_type = column
        .column_type
        .clone()
        .unwrap_or(ColumnType::Blob(None));
    builder = builder.column(ColumnSpec::new(column.name.clone(), column_type));
}
dst.create_layer(&builder)?;
Ok(()) }

If you leave z and m at their defaults, the dimension is declared prohibited, and the first geometry with a Z fails. If the source’s primary key is not named fid, also pass its name to TableSchemaBuilder::primary_key.

Copy the rows by their bytes

Feature::geometry_bytes returns the stored GPB blob, including its header. gpb::body_offset reads the header to find where the ISO WKB body starts. FeatureWriter::insert_wkb copies the body into the new blob without re-serialising it. The geometry is never parsed.

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::core::gpb;

use geopackage::GeoPackage;
let src = GeoPackage::open_read_only("places.gpkg")?;
let dst = GeoPackage::create("subset.gpkg")?;
let source = src.layer("points")?;
let target = dst.layer("points")?;
let mut writer = target.writer()?;
let mut cursor = source.cursor()?;

for feature in cursor.features()? {
    let feature = feature?;
    let values: Vec<_> = feature.values().collect();
    match feature.geometry_bytes() {
        Some(blob) => {
            let offset =
                gpb::body_offset(blob).map_err(|e| geopackage::Error::Core(e.into()))?;
            writer.insert_wkb(Some(feature.fid()), &blob[offset..], &values)?;
        }
        // A NULL geometry cell, or an attribute table.
        None => {
            writer.insert_row(Some(feature.fid()), &values)?;
        }
    }
}
writer.commit()?;
Ok(()) }

Feature::values yields borrowed ValueRefs, so text and blob cells bind into the destination statement without being copied.

To update rows in place instead of inserting them, use update_wkb. It gives the same guarantee.

For a long copy

A single writer is a single transaction, so the journal grows for the whole length of the copy. Commit periodically and open a new writer, as gpkg copy does every 10,000 rows. Add a counter to the loop above, and end each batch with:

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::GeoPackage;
let dst = GeoPackage::create("subset.gpkg")?;
let target = dst.layer("points")?;
let mut writer = target.writer()?;
let mut in_batch = 0usize;
const BATCH: usize = 10_000;
in_batch += 1;
if in_batch >= BATCH {
    writer.commit()?;
    writer = target.writer()?;
    in_batch = 0;
}
Ok(()) }

If the copy then fails part-way, only the current batch is rolled back, not the whole layer.

If you have opened your own transaction on GeoPackage::connection, these commits only stage the rows, and you must issue the durable commit yourself. See Transactions, and who commits.

Then

• Geometry storage: GPB, WKB, and curve types explains what the GPB header contains, why the WKB body passes through untouched, and what geo-traits cannot represent.
• API reference: Feature::geometry_bytes, FeatureWriter::insert_wkb, FeatureWriter::update_wkb, gpb::body_offset.

How to read a layer as Arrow record batches

The columnar path builds Arrow arrays directly from the statement’s column values, without creating a feature object for each row. This guide shows you how to read a layer as RecordBatches, filter the read, and pass the result to the rest of the Arrow ecosystem.

Enable the feature

The columnar paths are off by default, because they add the arrow-array and arrow-schema dependencies, which you don’t need if you only use the scalar API:

$ cargo add geopackage --features arrow

Read the whole layer

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::GeoPackage;
use geopackage::arrow::ArrowReadOptions;

let gpkg = GeoPackage::open_read_only("roads.gpkg")?;
let layer = gpkg.layer("roads")?;

for batch in layer.read_arrow(ArrowReadOptions::default())? {
    let batch = batch?;
    println!("{} rows", batch.num_rows());
}
Ok(()) }

Attribute columns follow the type mapping documented on the arrow module. The geometry column is WKB with the geoarrow.wkb extension name, and its metadata includes the CRS as PROJJSON. Batches arrive in primary-key order, whether or not the read is threaded.

Filter the read

Three variants take filters. All three are single-threaded:

• To read the rows intersecting a box, use read_arrow_in:

  fn main() -> Result<(), Box<dyn std::error::Error>> {
  use geopackage::{BoundingBox, GeoPackage};
  use geopackage::arrow::ArrowReadOptions;
  let gpkg = GeoPackage::open_read_only("roads.gpkg")?;
  let layer = gpkg.layer("roads")?;
  let bbox = BoundingBox::new(-7.0, 53.0, -6.0, 54.0);
  let batches = layer.read_arrow_in(bbox, ArrowReadOptions::default())?;
  Ok(()) }

  It uses the RTree index if the layer has one, and a full scan if not. Either way, it re-tests every candidate against its true f64 envelope. The candidate set is fixed when the read is opened, so a row inserted while you are still pulling batches is not returned.

• To filter on attributes, use read_arrow_where, which appends a raw SQL WHERE clause. The library trusts the clause and does not parse or sanitise it. Its placeholders are ?1 to ?N, bound in slice order:

  fn main() -> Result<(), Box<dyn std::error::Error>> {
  use geopackage::{GeoPackage, ValueRef};
  use geopackage::arrow::ArrowReadOptions;
  let gpkg = GeoPackage::open_read_only("roads.gpkg")?;
  let layer = gpkg.layer("roads")?;
  let batches = layer.read_arrow_where(
      "class = ?1",
      &[ValueRef::Text("motorway")],
      ArrowReadOptions::default(),
  )?;
  Ok(()) }

  You can also use it to read a single row in columnar form: pass fid = ?1 with the key as its parameter.

• To apply both at once, use read_arrow_in_where, which intersects the two.

Tune the read

ArrowReadOptions controls three things:

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::arrow::ArrowReadOptions;

let options = ArrowReadOptions::with_batch_size(16_384)
    .with_threads(1)
    .with_max_batch_bytes(256 * 1024 * 1024);
Ok(()) }

with_batch_size creates a new set of options from the defaults. The other two chain onto any options value.

• batch_size is rows per batch, 65,536 by default.
• threads defaults to min(4, available parallelism). Set it to 1 to keep the read on your own thread. The read uses threads only when the database is a file, the primary key is dense, there are at least two batches of rows, and the layer reads every column (no with_columns or without_geometry). Otherwise it runs single-threaded rather than failing.
• max_batch_bytes limits the geometry bytes in one batch, and defaults to min(INT32_MAX, RAM / 4). The geometry column’s Arrow offsets are 32-bit, so 2 GB is a hard ceiling. A batch that would cross the limit is emitted with fewer rows. Lower it for a layer of very large geometries, where the default multiplied by the thread count is more memory than you want in use at once.

Hand the batches to something else

ArrowBatches implements RecordBatchReader, so you can pass it to anything in the Arrow ecosystem that consumes one, such as a Parquet writer or the Arrow C Data Interface:

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::GeoPackage;
use geopackage::arrow::ArrowReadOptions;
use arrow_array::RecordBatchReader;

let gpkg = GeoPackage::open_read_only("roads.gpkg")?;
let layer = gpkg.layer("roads")?;
let batches = layer.read_arrow(ArrowReadOptions::default())?;
let schema = batches.schema();
Ok(()) }

Copy a layer through Arrow

Layer::arrow_schema and TableSchemaBuilder::from_arrow_schema apply the same type mapping in opposite directions, so you can copy a layer without restating its schema:

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::{GeoPackage, TableSchemaBuilder};
use geopackage::arrow::ArrowReadOptions;

let src = GeoPackage::open_read_only("roads.gpkg")?;
let roads = src.layer("roads")?;

let dst = GeoPackage::create("copy.gpkg")?;
let schema = roads.arrow_schema()?;
dst.create_layer(&TableSchemaBuilder::new("roads").from_arrow_schema(&schema)?)?;

let batches = roads.read_arrow(ArrowReadOptions::default())?;
dst.layer("roads")?.write_arrow(batches, 0)?;
Ok(()) }

Here the geometry column is declared with its source type. A layer built from an Arrow schema alone gets GEOMETRY, because the GeoArrow WKB encoding does not record a specific type.

Then

• Each batch is a separate query, paginated on the primary key, so a concurrent writer can change the table between batches. If you need one snapshot of the whole layer, open your own transaction on GeoPackage::connection. See Transactions, and who commits.
• API reference: Layer::read_arrow, Layer::read_arrow_in, Layer::read_arrow_where, ArrowReadOptions, and the arrow module for the type mapping.

How to validate a GeoPackage and act on the findings

GeoPackage::validate checks a file in one pass and returns what is wrong with it. It modifies nothing. This guide shows you how to run the check, read the result, and apply the repairs it names.

The report covers only what these checks can see. It does not establish conformance in every respect the spec defines. For that, the OGC executable test suite remains the authority.

In a script or a CI job

$ gpkg validate places.gpkg
places.gpkg
  warning: spatial index on "points" is maintained by a pre-1.4 or mixed trigger set
    repair: upgrade the trigger set with Layer::repair_spatial_index

  0 errors, 1 warning, 0 advisories

Exit codes:

• Zero when there are no findings, or only warnings and advisories.
• Non-zero when any finding is an error. An error means a reader can get a wrong answer from the file.

To fail on warnings too, add --strict:

$ gpkg validate --strict places.gpkg
$ echo $?
1

The command opens the file read-only and leniently, so a file with problems is reported on rather than rejected.

In Rust

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::{GeoPackage, Severity};

let gpkg = GeoPackage::open_read_only_lenient("places.gpkg")?;
let findings = gpkg.validate()?;

for finding in &findings {
    println!("{}: {finding}", finding.severity());
    if let Some(repair) = finding.repair() {
        println!("  repair: {repair}");
    }
}

let failed = findings.iter().any(|f| f.severity() == Severity::Error);
Ok(()) }

Findings are sorted most severe first, in a stable order within each severity. If you diff the output of two runs, any differences are in the findings, not in their order.

Read the severities

Severity reflects the consequence of a problem, not how hard it is to fix:

• Severity::Error: a reader can get a wrong answer. For example, a gpkg_contents row naming a missing table, a spatial index that no longer describes its rows, or a metadata or relation row that points at something that no longer exists.
• Severity::Warning: the file reads correctly but does not follow the current spec. For example, a pre-1.2 application_id, a pre-1.4 index trigger set, an extension these checks cannot identify, or a tile pyramid that breaks the matrix rules.
• Severity::Advisory: a remark, not a defect. The advisory case is a feature table with no spatial index. It reads correctly, and someone may have chosen not to index the layer.

Apply the repairs

Finding::repair returns the advice as text, naming the method that performs it. None means this library cannot make the fix. It needs the writer that produced the file, or a decision about the data that the library should not make for you. To act on the findings instead of printing them, match on the variant:

fn main() -> Result<(), Box<dyn std::error::Error>> {
use geopackage::{Finding, GeoPackage};

let gpkg = GeoPackage::open_lenient("places.gpkg")?;
for finding in gpkg.validate()? {
    match finding {
        // Contents disagree with the geometries: rebuild.
        Finding::SpatialIndexOutOfStep { table_name, .. } => {
            gpkg.layer(&table_name)?.rebuild_spatial_index()?;
        }
        // Structure is a pre-1.4 or mixed trigger set: repair.
        Finding::LegacySpatialIndexTriggers { table_name } => {
            gpkg.layer(&table_name)?.repair_spatial_index()?;
        }
        // An advisory, so build one only if you want the queries indexed.
        Finding::NoSpatialIndex { table_name } => {
            gpkg.layer(&table_name)?.create_spatial_index()?;
        }
        _ => {}
    }
}
Ok(()) }

The handle must be open read-write for any of these repairs. Finding is #[non_exhaustive], so a variant added in a later version falls through to the catch-all arm instead of breaking the build.

To repair indexes without validating the whole file first, see How to repair a file’s spatial indexes.

Then

• The spatial index: structure, contents, and the bulk build explains what an audit checks that a status cannot.
• API reference: GeoPackage::validate, Finding, Severity.

The layer extent, and why it cannot be trusted

Every row of gpkg_contents has four bounding-box columns: min_x, min_y, max_x and max_y. They look like the layer’s extent, and a lot of software treats them as the extent. They are the least reliable values in a GeoPackage, by the standard’s design rather than through neglect.

What the spec says

The relevant text is prose, not a numbered requirement, and it is identical in every version from 1.2.0 to 1.4.0. The bounding box “provides an informative bounding box of the content”. Applications “may use this bounding box as the extents of a default view”. And “there are no requirements that this bounding box be exact or represent the minimum bounding box of the content”. The columns are nullable, and no requirement governs their values.

This is intentional. The GeoPackage SWG considered adding a test that features fall inside the recorded extent, and declined in 2018. Maintaining that invariant on every insert and delete is too expensive to demand of a writer.

So if you read files you did not produce, a file can record an extent that is stale, inflated or simply wrong, and still be conformant.

No query uses the recorded extent

No read path in this library consults the recorded extent. features_in and cursor_in go through the RTree index, or through a full scan when there is none, and re-filter every candidate against its exact f64 envelope. The recorded bounds never short-circuit anything.

An obvious optimisation would be to skip the scan and return nothing when the query box misses the recorded extent. The library does not do this. It would be correct for every file this library wrote and wrong for files it did not. On those files, a query would silently return no rows even though matching rows exist. A read path’s correctness should not depend on which tool wrote the file.

Why a wrong extent is worse than an absent one

Readers treat a NULL extent and a wrong extent very differently, so the nullable columns matter.

GDAL validates only that the four values are present and that min <= max per axis. It returns a well-ordered but wrong extent verbatim and never recomputes it, even when the caller passes bForce. A NULL extent, by contrast, makes GDAL compute the true one. QGIS behaves the same way. Since 3.34 its “Update Extents” action is a no-op for a local GeoPackage, so a user who receives a file with a wrong extent cannot correct it without a full edit session.

An absent extent costs a reader one scan, and then the reader fixes it. Every reader believes a wrong extent indefinitely, and the users most likely to notice have no simple way to fix it.

So this library follows one rule: never record an extent that is not guaranteed to cover the layer. NULL is spec-legal and accurate, and readers repair it. A wrong box is none of those things. Four things follow from the rule.

The writer grows the recorded box to cover what it writes, and never shrinks it. The result is exact or an over-estimate, which the spec expressly permits. A delete leaves the box alone and does not trigger a rescan.

A writer that cannot rely on its starting point leaves the extent alone. It does not replace it with a box covering only the rows it wrote. This is the case where the stored extent is absent, NULL or inverted while the table already contains rows. A fold over what this writer wrote is a lower bound, not the true extent. Recording it would replace an accurate “unknown” with a box that excludes every pre-existing row.

A layer with no geometry to measure gets NULL, not an invented value. GDAL does the same, through UpdateContentsToNullExtent.

A measurement is recorded only when the whole layer could be measured inside one transaction. If it loses that race to another writer, the measurement describes a layer that changed underneath it, and the rule forbids recording that. So lock contention is not an error here. The file keeps what it had, and the caller still gets the measurement.

Reading an extent can write one

Layer::extent records what it had to measure. If a file’s recorded extent is unusable, reading the extent modifies the file. Both its content and its modification time change.

This is deliberate, and it is the only place where this library writes during what looks like a read. Later readers, in this library or elsewhere, find a usable box and do not have to scan for it. GDAL does the same, persisting a computed extent through SaveExtent on any dataset open for update.

The cost falls on some callers. For example, a pipeline that checksums its inputs will see them change. There are two cheap ways to avoid this. Open the file read-only, which suppresses the write and still returns the measurement. Or read GeoPackage::contents, which reports the recorded values without measuring or writing anything.

Compared with index repair

Spatial-index repair, by contrast, never happens unless you request it. The two can look like the same question, but they differ in cost and in what happens when they are missing.

Rebuilding an index reads every geometry in the layer and rewrites the tree. That is too expensive to run as the side effect of a question. And a stale or older-generation index is not silently believed: has_spatial_index reports it as unusable, and queries fall back to a correct full scan. An extent is the opposite on both counts. Measuring one is comparatively cheap, and every reader that opens the file believes a wrong one indefinitely. See The spatial index for the other half of that comparison.

One deviation from GDAL

GDAL prefers the RTree when it has to compute an extent. That is faster, but it yields a box rounded outward to f32, and it trusts the index about what the layer contains.

This library always measures the geometries. The value will be written into the file and believed indefinitely, so it is worth measuring exactly. Measuring also removes any dependence on the index being current, and many files arrive with indexes that are not.

Transactions, and who commits

SQLite has no nested transactions. A BEGIN issued while one is already open fails with “cannot start a transaction within a transaction”. This shapes how every write path in this library starts its work. You should understand it before you open your own transaction on the underlying connection.

The original problem

Originally, every write path here issued a BEGIN unconditionally. That works while the library is the only thing driving the connection. It fails as soon as a caller uses the escape hatch.

GeoPackage::connection returns the underlying rusqlite connection, because SQLite is the query engine and anything the API does not cover should be a query away. A caller who began a transaction on that connection and then used the ordinary API got the nesting error instead of their writes. The C ABI had a worse version of the problem. It could not offer gpkg_begin and gpkg_commit at all, because the sequence a C consumer would reasonably write could not work.

Inheriting rather than nesting

Now a write path checks whether a transaction is already open. If one is, it runs its statements inside it.

Atomicity is unchanged. The statements are still grouped, they still commit together, and a failure still rolls them back together. What changes is who decides when that happens: whoever issued the BEGIN.

Every write path in the library behaves this way, from a single FeatureWriter insert to a bulk write_all and the extent recording described in the extent chapter.

Why not savepoints

Savepoints are the other way to get nesting out of SQLite. They were considered and rejected.

rusqlite::Connection::savepoint needs &mut Connection, and this library has only &Connection, because the read path shares it. Savepoints would therefore have to be driven as raw SAVEPOINT and RELEASE text. That is more machinery, for a finer rollback granularity that nobody has asked for. The question is not closed. If a need arises, the transaction type is the one place that would have to change.

What inheritance costs the caller

These consequences apply to every write path, not only the feature writer.

Commit does not commit when the transaction was inherited. It still does everything else a commit does, including flushing the gpkg_contents last_change and bounding box, but the caller issues the durable commit.

Dropping a writer does not roll back when the transaction was inherited. An error part-way through leaves partial work staged in the caller’s transaction, and the caller has to discard it.

write_all’s batch_size stops bounding transactions, because every batch belongs to the caller’s transaction. It bounds nothing else either: the rows are written in the same order, through the same statements.

A writer does not expose whether its transaction was inherited, and this is deliberate. A caller who opened a transaction knows they did. A caller who did not cannot get this behaviour at all.

Why the writer owns a transaction rather than taking one

The writer could instead accept a rusqlite::Transaction from the caller. That design was rejected for two reasons.

The first is API surface. rusqlite types do not appear in the public API here. Geometry is impl geo_traits::GeometryTrait<T = f64> and non-geometry values are this library’s own value types, so a major-version bump in rusqlite does not cascade into every caller. Accepting a transaction object would put that type in the signature of the most-used write path in the library.

The second is bookkeeping. The writer maintains a running bounding-box fold and a last_change timestamp, and flushes both at one point. Owning the transaction makes that point well defined.

You can still drive your own transaction through the raw connection. Inheritance is what makes that work.

Beginning without locking

Layer::writer opens a transaction without writing anything, and does not take a lock. SQLite’s BEGIN DEFERRED acquires nothing until the first statement needs it. So the first failure a caller sees comes from the first row, not from the call that returned the writer.

This matters most in two cases. On a read-only connection, opening a writer succeeds and writing the first row fails. Under contention, the wait for another connection’s lock happens at the first write, not when the writer opens. So a successful call to open a writer does not show that the file is writable.

Reading and writing on one connection at once

A writer and a feature cursor share the connection, so a scan can drive its own updates: read a row, recompute a column, write it back.

SQLite does not define what such a scan sees. Its isolation documentation says explicitly that a SELECT on one connection has no isolation from writes on that same connection, and that an application “can UPDATE the current row or any prior row, though doing so might cause that row to reappear in a subsequent sqlite3_step()”. SQLite promises only that the file will not be harmed. It does not promise that the result set is stable, complete, or free of repeats.

In practice the scan is stable when all three of these are true:

• the cursor is a plain table scan, not one driven by an index;
• the columns written are not ones the scan’s index reads;
• the primary key is not written, since moving a row’s id moves it within a rowid scan.

Avoid writing a geometry during a bounding-box cursor. That cursor is driven by a join against the RTree, and writing a geometry moves the row inside that index through the triggers. This is the pattern that makes a scan return rows it has already returned. To recompute geometries, use two passes: collect the feature ids, finish the scan, then write.

None of this is specific to this library, and none of it is a defect that could be fixed here. It is SQLite’s stated contract for one connection reading and writing at once. The alternative, a second connection, has its own locking to manage.

Getting the work done

• How to copy features between files without decoding geometry shows a long write committed in batches, and what changes when the caller owns the transaction.
• GeoPackage::connection is the escape hatch this chapter is about.

The spatial index: structure, contents, and the bulk build

A GeoPackage indexes a feature layer with an SQLite RTree virtual table named rtree_<table>_<column>. A set of triggers, written out in Annex F.3 of the spec, keeps it in step with the rows. The table and the triggers can each go wrong, in different ways, and checking each costs a different amount.

Structure and contents

The first check is structural: does the virtual table exist, and which generation of trigger set maintains it? The SQLite catalogue answers this cheaply, and repair acts on the answer. It cannot tell you whether the entries in the index agree with the geometries in the table.

The second check is about contents: does the index describe exactly the rows it should, and cover each one? Answering it means reading every geometry in the layer and comparing. It reports the indexable rows, the entries, the rows with no entry, the entries whose row is gone, and the entries that fail to cover their geometry’s true envelope.

An index can pass the first check and fail the second. This happens in a file where rows were written while the triggers were absent, or that another tool populated incompletely. This library’s own writes aren’t expected to produce it, but with bulk verification off, which is the default, nothing checks each bulk build for it. Because an index can pass one check and fail the other, the two have separate remedies. Repair puts the structure right and rebuilds the content, and does nothing when the structure is already current. Rebuild does the content work unconditionally.

A file can arrive in one of four structural states:

• Absent: neither the virtual table nor any trigger is there.
• Current: the virtual table exists with the GeoPackage 1.4 trigger set.
• Legacy: a pre-1.4 or mixed trigger set.
• Stale: the two halves disagree, with a virtual table but no triggers, or triggers but no table.

A legacy trigger set is a real problem, not just untidy. The pre-1.4 update1 trigger corrupts the index under UPSERT. That is why 1.4 renamed the fixed triggers, and why the repaired state can be detected by name. QGIS tracked the same failure in QGIS#36935. A mixed set, with triggers from more than one generation, is the worst case. You get one when two tools have edited a file in turn.

This library does not trust a Legacy or Stale index. has_spatial_index reports both as unusable, so a bounding-box query falls back to a correct full scan instead of returning whatever the suspect index says.

Why repair is never automatic, when extent recording is

This library writes during a read in one place: reading a layer’s extent records the measurement when the recorded bounds are unusable. It does not do the same for an index, and the difference is deliberate.

The two operations differ in cost, and in what happens without them. Rebuilding an index reads every geometry in the layer and rewrites the tree. That is too expensive to run as the side effect of a question, and a caller who asked about an index’s state has not consented to it. Without a usable index, queries fall back to a scan and still return the correct rows, more slowly.

An extent is the opposite on both counts. Measuring it is comparatively cheap, and every reader that opens the file believes a wrong one indefinitely. See The layer extent for that side of the comparison.

There is a second reason to keep repair explicit. Repairing a legacy trigger set rewrites triggers in a file the caller may only have opened to look at. Doing that unasked would change files on disk just because someone inspected them.

Populating an index

Once the virtual table and the triggers exist, the index has to be filled from the rows already in the table. There are two ways to do it.

Below a threshold, currently 10,000 candidate rows, population is a single INSERT INTO rtree SELECT over the existing rows. It uses the registered ST_* functions and skips empty and NULL geometries exactly as the triggers do. It is simple, it goes through SQLite’s own RTree module, and for a small table it is cheaper.

At or above that threshold, the library uses the bulk build.

The bulk build

SQLite’s RTree module has no bulk-load entry point. Its only write path is the per-row virtual-table update, so every insert descends the tree and may split nodes. Populating an index one row at a time pays that cost for every row. GDAL reached the same conclusion in gdal#7614.

The obvious alternative is to build the RTree in an ATTACHed scratch database, through the RTree module, and copy its shadow tables into the target. Measurements show that this build has the same dominant cost, because it still inserts every entry through the module one row at a time. The work moves but does not go away.

So the library builds the tree directly. It lays out the rtree_%_node, _rowid and _parent contents in memory from the accumulated entry set, and writes them as ordinary rows. No RTree module logic, no ST_* calls, no triggers, and no scratch database. This works because the GeoPackage RTree is always rtree(id, minx, maxx, miny, maxy), two-dimensional with no auxiliary columns. The three shadow tables therefore have a fixed shape whatever the index is named.

Avoiding the scratch database also matters for a reason more important than speed. ATTACH requires autocommit, which would force the build out of the caller’s transaction and leave a window in which a crash commits the rows but not the index. Built directly, the whole build is one transaction: dropping the triggers, every row insert, the gpkg_contents flush, the index work and reinstalling the triggers all commit together. The rows can never be committed against an index that was not brought up to date with them. The exception is a write inside a transaction you opened yourself. A failed call leaves its work staged there, and committing it anyway leaves a stale index.

The entry set comes either from a scan of the table, which parses each geometry once in Rust, or, when the caller can account for every indexable row, from the envelopes computed while encoding the geometries. Because of that second source, building the index during a write costs between 20% and 47% less than building it afterwards. It is also why creating a layer leaves an empty index in place for the first bulk write to fill.

The verification gate, and its cost

A hand-written RTree uses an on-disk format that SQLite does not document as an interface, so the library can check the result before anything relies on it.

The check has three levels. The first reads the index back and requires exactly the accumulated (fid, envelope) set, by row count and by a per-row containment check. The second adds rtreecheck on the index, and the third adds a whole-database PRAGMA integrity_check. Any anomaly discards the built index and falls back to the triggered population, so a failed check costs time, not correctness.

The cost is substantial. Profiled over one million points, the second level is about 45% of the build: roughly 745 ms of a 1593 ms build. The time is split about evenly between the bijection scan, which reads every entry back out of the index, and rtreecheck, which walks the whole tree. GDAL’s builder runs no equivalent check. With that check on, this library is level with GDAL, not ahead of it.

Verification is therefore opt-in, and off by default. With verification off, nothing is read back, so a build cannot fail its check and cannot fall back to the triggered path. The automatic fallback exists only where a check does. If you are writing files that matter, or bisecting a suspected index problem, turn verification on to get both.

Node fill

By default the bulk build packs each RTree node full. This gives the smallest tree, the shallowest descent and the best queries, and suits a bulk load.

A full node has no room for a later insert, so the first append into it splits immediately. Appends after a bulk load go through the triggers, which is the per-row path the bulk build exists to avoid. The default therefore favours the load and the queries. A lower fill factor suits a freshly built index that you are about to append to heavily.

Getting the work done

• How to add a spatial index to an existing layer
• How to repair a file’s spatial indexes

Geometry storage: GPB, WKB, and curve types

A GeoPackage stores each geometry as a blob in an ordinary SQLite column, in an encoding the spec calls GeoPackage Binary, or GPB. A GPB blob is a small header followed by a body in ISO Well-Known Binary (WKB). Most of how this library reads and writes geometry follows from the layout of that header, and from one thing WKB can express that the Rust geometry traits cannot.

The blob

The header is eight bytes: the magic GP, a version byte, a flags byte, and a 32-bit srs_id. The flags byte records the byte order of the header and envelope, an envelope indicator, an empty-geometry flag, and an extended-GPB flag. After the header comes an optional envelope of four, six or eight doubles depending on the indicator, and after that the WKB body, which declares its own byte order and includes its own type code.

The envelope is optional in the format, but this library always writes one: XY, or XYZ when the geometry has a Z. It does this for write throughput. The spec’s RTree triggers call the ST_MinX family of functions four times per row. These functions can answer from a header envelope in constant time, and otherwise have to traverse the WKB body. Writing 32 extra bytes per row is cheap compared with traversing every point of every geometry on every insert.

For points the trade looks unfavourable, because a point’s envelope is larger than the point. This library writes one for points anyway, because uniformity is worth more than the saving. GDAL omits envelopes for points, and that is the practical reason a traversal fallback has to exist at all. Blobs without envelopes are common in files this library did not write, and they still need bounds computed for them.

Reading without materialising

A geometry read from a row is a lazy wrapper, which pairs the parsed header with the WKB body that follows it. The body is read through the georust wkb crate, and the wrapper implements geo_traits::GeometryTrait by delegating to it.

This means you can traverse coordinates without ever creating an owned geometry object. An algorithm written against geo_traits reads coordinates directly from the stored bytes. The write path accepts any impl GeometryTrait<T = f64>, so you can stream a geometry out of one file, measure it, and write it into another without converting it to a geo-types value in either direction. Conversion to geo-types is available behind a default feature if you want it. You choose when to convert, and the read path does not do it for you.

Two things still allocate: each row’s blob, which is copied out of SQLite, and the new blob the writer serialises. An algorithm that produces new geometry also allocates its output. What you avoid is the intermediate object, not every allocation.

Parsing arbitrary bytes never panics. A malformed header, a truncated body, or a type the reader cannot handle all produce a typed error. There is one known limitation close to that guarantee. The wkb reader pre-allocates from element counts read out of the blob without checking them against the buffer size, so a malformed geometry declaring a 0xFFFFFFFF-member collection causes a very large allocation. The fix belongs upstream. Until it lands, be careful when parsing GeoPackage files from untrusted sources.

Bytes in, bytes out

Writes normally go through the trait, but you can also pass WKB bytes to the writer directly. They are copied into the new blob rather than re-serialised, so a body read out of one GeoPackage passes into another unchanged, with only the header rebuilt around it. The rest of this chapter explains why this path exists.

What the traits cannot represent

The GeoPackage geometry vocabulary in Annex G includes five non-linear types: CIRCULARSTRING, COMPOUNDCURVE, CURVEPOLYGON, MULTICURVE and MULTISURFACE. Their defining feature is the circular arc, which is given by three points. The first and last are the arc’s ends, and the middle one fixes the circle the arc lies on.

geo-traits, the interface this library uses to read and write geometry, has no representation for an arc, and neither does geo-types. This is not an oversight in either. An arc is not a sequence of coordinates, and a library built on sequences of coordinates has nowhere to put one. Densifying an arc into a linestring to approximate it would be a lossy transformation that no caller requested.

This library therefore handles the non-linear types as bytes. You can write them, index them and query them by extent, and their WKB is returned as stored. Only the call that would return a geometry object fails on them, because there is no object to return. For this reason, anything that moves geometry between files should use the byte path. Reading into a geo-types value and writing it back loses the curve types entirely, and also parses every coordinate for no purpose.

Why an arc’s envelope has to be exact

Indexing and querying a curve by extent requires an envelope, and that envelope has to be computed from the arc rather than from the three points that define it.

A circular arc bulges away from the chord joining its endpoints, and may bulge past its middle control point too. The bounding box of the three control points can therefore be strictly smaller than the arc’s own. An envelope that is too small is a correctness problem, not a tuning or quality question. The GPB header envelope and the RTree entry are both derived from it, and a reader that trusts either will silently drop features that do intersect the query window.

The exact box can be computed in closed form. The arc lies on a circle through the three points, so its extremes in x and y are either the endpoints or the points of that circle furthest along each axis. The chord between the endpoints decides which of those four candidates the arc actually reaches. This is the approach PostGIS takes in lw_arc_calculate_gbox_cartesian_2d. GDAL computes the same box by a different route, from swept angles, in OGRCircularString::ExtendEnvelopeWithCircular.

A degenerate case explains why the arithmetic uses a tolerance. As the three points approach collinearity, cancellation error dominates the circumcentre and the radius diverges. In that limit the arc is treated as its chord. This is correct, not a fudge. As the points straighten, the bulge past them falls to zero, so the error the tolerance admits shrinks with it.

Arcs are planar by definition

An arc is a circle in the coordinate space of the layer’s CRS, including when that CRS is geographic. Three longitude and latitude points define a circle in degree space, not a small circle on the ellipsoid.

This is not peculiar to this library. PostGIS and GDAL both read arcs the same way, and the GeoPackage RTree indexes the same degree space. A geodesic reading would produce envelopes that disagreed with the index they were stored in. The envelope computation therefore uses no geodesic arithmetic. Nor does it involve a coordinate transformation, which this library never performs.

Tight envelopes, and where slack is added

The computed box is the minimum bounding box, with no outward margin. Slack is added later, in the index, where it can be reasoned about. RTree columns are 32-bit floats, so each bound is rounded outward when it is narrowed. A stored cell therefore always contains the f64 box it was derived from. Annex F.3 of the spec also tells clients to expand query windows to absorb exactly that rounding.

Keeping the geometry-level envelope tight and the index-level entry loose puts the widening in one place, where it is a property of the index rather than of the data.

Two envelope walkers

As a result, the codebase has two envelope computations rather than one. The ordinary path traverses a geometry through the wkb reader, which cannot read the non-linear types. The curve path walks the WKB structure itself, byte by byte. It handles every Annex G type and every XY, XYZ, XYM and XYZM variant, so it yields an envelope for a curve body without waiting on upstream support.

The curve walker also has some defensive properties. It skips non-finite coordinates, so an empty point, which by convention has NaN coordinates, yields no envelope rather than a NaN-valued box. It caps nesting at 32 containers deep, so a hostile body cannot drive the recursion into a stack overflow. It allocates nothing, so a bad element count in a curve body fails immediately instead of reserving memory for it. The reader used for the linear types does not behave this way.

The split creates one error case of its own: a body whose own type is a core type but which contains a non-linear member, such as a GEOMETRYCOLLECTION containing a CIRCULARSTRING. The body’s own type code decides which walker it gets, so such a container reaches the reader that cannot read the member. This is reported as a distinct error rather than as a malformed body, because the body is not malformed, and otherwise no caller could tell the two apart.

Getting the work done

• How to copy features between files without decoding geometry

Tiled gridded coverages, and why they are not tile pyramids

A GeoPackage can store two kinds of tiled raster. One is a basemap: PNG or JPEG tiles that a renderer draws. The other is a tiled gridded coverage: a grid of measurements, usually elevation, stored tile by tile in the same pyramid structure. The two look similar on disk, but their uses are different, so this library keeps them separate.

What the extension adds

Coverages are an extension to the core standard, published separately as OGC 17-066r2 and registered in a file as gpkg_2d_gridded_coverage. The extension adds three things to a tile pyramid:

• a different gpkg_contents.data_type: 2d-gridded-coverage, not tiles;
• gpkg_2d_gridded_coverage_ancillary, one row for each coverage, which states what the samples mean: integer or float, the scale and offset that convert a stored sample to a value, the value for “no data”, and the unit of measure;
• gpkg_2d_gridded_tile_ancillary, one row for each tile, with the scale and offset of that tile and, optionally, its minimum, maximum, mean and standard deviation.

The payload of a tile is a TIFF or a 16-bit greyscale PNG, and the extension constrains it closely: one sample per grid cell, 32-bit float or 8, 16 or 32-bit integer samples, no compression or LZW compression, one image per file, and no internal tiles.

A coverage payload is not a GeoTIFF

The TIFF payloads of a coverage are often assumed to be GeoTIFFs. They are not.

The georeferencing is in gpkg_tile_matrix_set, the same extent and grid that a basemap pyramid uses, so the payload does not need any. The extension does not use GeoKeyDirectory, ModelPixelScale or ModelTiepoint. It specifies a small profile of baseline TIFF, and each rule in the profile is a statement about a tag, not about a pixel.

So this library can check a coverage payload completely without an image codec. geopackage_core::coverage::coverage_tiff walks the first IFD, reads seven tags, and checks the encoding requirements against them. coverage_png does the same for IHDR. Neither function decodes a sample.

The limits of a header check

This library cannot check one requirement. Requirement 21 says that every pixel in a tile must have a valid value, and that NaN and Inf must not occur. That is a statement about samples, and a check of the samples needs a decoder. A payload that this library accepts conforms in structure only, and the check says nothing about its contents.

The same limit applies to the write path. The four per-tile statistics come from the samples, so they are None unless you supply them:

use geopackage::core::tiles::TileCoord;
use geopackage::{GeoPackage, TileAncillary};
fn main() -> Result<(), Box<dyn std::error::Error>> {
let gpkg = GeoPackage::open("dem.gpkg")?;
let coverage = gpkg.coverage("elevation")?;
let payload: Vec<u8> = Vec::new();
let mut writer = coverage.writer()?;
writer.put_with_ancillary(
    TileCoord::new(0, 0, 0),
    &payload,
    // Values from the code that decoded the samples. Without them, the row has
    // nulls. That conforms, but the row has less data than a row from GDAL.
    &TileAncillary::defaults().with_statistics(0.5, 96.5, 48.53, 27.99),
)?;
writer.commit()?;
Ok(()) }

Why a coverage has a separate handle

GeoPackage::tiles rejects a coverage, and GeoPackage::coverage rejects a basemap. There are three reasons.

Every tile needs a second row. Requirement 10 specifies a matching gpkg_2d_gridded_tile_ancillary row for each tile in the user table. The normative table definition has no ON DELETE CASCADE, so the database does not keep the rows paired. CoverageWriter writes both rows and deletes both rows. A TilePyramid writer does not know about the second row, and a coverage written through it would collect tiles without rows.

The payload rules are opposites. A tile pyramid can store PNG, JPEG or WebP, and must not store TIFF. A coverage must store TIFF or 16-bit PNG, and nothing else. Both checks are strict, and neither is a special case of the other.

The result of a read is different. From a basemap, you need bytes for a decoder. From a coverage, you need a value: the stored sample, converted by two scale and offset pairs, and compared with the null:

use geopackage::core::tiles::TileCoord;
use geopackage::GeoPackage;
fn main() -> Result<(), Box<dyn std::error::Error>> {
let gpkg = GeoPackage::open_read_only("dem.gpkg")?;
let coverage = gpkg.coverage("elevation")?;
let tile = coverage.tile_ancillary(TileCoord::new(0, 0, 0))?.unwrap();

// `sample` comes from a decode of the payload, which your code does.
let sample = 1234.0;
if coverage.is_null(sample) {
    println!("no data here");
} else {
    println!("{} {}", coverage.value(&tile, sample),
             coverage.ancillary().uom.as_deref().unwrap_or(""));
}
Ok(()) }

In this example, the null is compared with the stored sample, before the scale and offset apply, because the extension says that the scale and offset do not apply to the null. If you apply them first, you lose your no-data cells, and nothing reports an error.

Files from other writers

As elsewhere in this library, a read accepts more than a write. A coverage with a datatype value that the requirement does not permit still opens, and Coverage::datatype returns None. A read accepts a payload with Deflate compression, although Requirement 15 specifies baseline TIFF and Deflate is not part of it. A write rejects the same payload.

GeoPackage::validate reports the defects in a file, requirement by requirement: ancillary rows that point to no tile, tiles without a row, a float coverage that scales its samples, and a payload that contradicts the datatype of its coverage. The last is an error, not a warning, because a reader that uses the ancillary row gets incorrect values.

Before you call validate on a large file, note that the coverage check is the only check with work that scales with the size of the file, not with the number of tables. The encoding requirements are statements about payloads, so the check must read every payload.

Reference

This page points to the reference documentation, which is generated from the source.

Crates

The C API

The header, geopackage-ffi/include/geopackage.h, is generated with cbindgen and checked in. cargo cinstall -p geopackage-ffi installs it with a versioned library and a pkg-config file named geopackage.

geopackage-ffi/examples/ has five C programs. CI compiles each against the committed header and runs it.

gpkg commands

Every command opens the file leniently. gpkg validate exits non-zero when a finding is an error. Every other command exits non-zero only when it fails.

Format and toolchain

• Specification: OGC GeoPackage 1.4.
• Minimum supported Rust version: 1.95.