GeoSilo

What you can do with one column type

The single most useful pattern: take an existing GEOMETRY column and re-declare it as GEOSILO with ZSTD compression. Every existing ST_* query still works — and the introspection / bounding-box / measurement calls get faster:

INSTALL spatial; LOAD spatial;
INSTALL geosilo FROM community; LOAD geosilo;

-- Compact geometry, transparent ST_* surface
CREATE TABLE parcels (
  id    INTEGER,
  name  VARCHAR,
  geom  GEOSILO('EPSG:4326') USING COMPRESSION zstd
);

INSERT INTO parcels SELECT id, name, geom FROM raw_parcels;

-- Bounding-box prefilter — entirely native, no WKB decode
SELECT id, ST_Area(geom)
FROM parcels
WHERE ST_XMin(geom) > -78
  AND ST_XMax(geom) < -76;

The GEOSILO type delta-encodes coordinates as scaled integers — each subsequent vertex of a ring is a 2-byte int16 offset from the previous one. On the TIGER/Line 2025 US Census dataset, GEOSILO + ZSTD is ~21% smaller than DuckDB’s columnar GEOMETRY storage.

Setup

GeoSilo overlays spatial’s ST_* namespace, so load both — spatial first:

INSTALL spatial; LOAD spatial;
INSTALL geosilo FROM community; LOAD geosilo;

The GEOSILO column type

Declare a column as GEOSILO, optionally tagging the CRS. With a CRS the integer scale is auto-detected (degrees → 1e7, meters → 100); without one, the default is 1e7:

-- With CRS — scale auto-detected
CREATE TABLE parcels (
  id    INTEGER,
  name  VARCHAR,
  geom  GEOSILO('EPSG:4326') USING COMPRESSION zstd
);

-- Without CRS — default scale 1e7
CREATE TABLE shapes (geom GEOSILO USING COMPRESSION zstd);

USING COMPRESSION zstd is the recommended setting — it’s where the storage win comes from. See CREATE TABLE for the full DDL.

The implicit GEOMETRY → GEOSILO cast auto-encodes on insert, so you can populate from any GEOMETRY source — Parquet, CSV-then-ST_GeomFromText, another DuckDB table:

INSERT INTO parcels SELECT id, geom FROM raw_data;

What’s native, what auto-decodes

Standard spatial functions all work on GEOSILO columns. The ones below have native GEOSILO overloads that read the silo header or walk the integer delta stream — no WKB round-trip:

Everything else (ST_Buffer, ST_Union, ST_Intersection, ST_AsText, ST_Centroid, etc.) auto-decodes via the implicit cast and runs through the GEOS-backed spatial implementation.

Encoding format, briefly

The format is verifiable and small enough to summarize:

1. Coordinates scale to integers (-75.509491 × 1e7 = -755094910).
2. The first vertex of each ring is a 4-byte int32 absolute coordinate.
3. Subsequent vertices are stored as 2-byte int16 deltas from the previous; ~90% fit. The rest use a 6-byte escape sequence.
4. Polygon rings omit the closing duplicate vertex — it’s reconstructed on decode.
5. Points use a 1-byte compact header.

The header (geometry_type, vertex_type, scale) is readable without decoding the body — see geosilo_metadata.

Scale and precision

Auto-applied when you declare GEOSILO('EPSG:...'). Override per-call with the two-argument form:

SELECT geosilo_encode(geom, 1000) FROM millimeter_data;

Arrow IPC transport

GeoSilo registers an Arrow extension type named queryfarm.geosilo. A producer tags an Arrow binary field with the right metadata and the DuckDB consumer (with GeoSilo loaded) auto-decodes silo blobs to GEOMETRY on read and re-encodes on write — CRS rides through the field’s metadata:

geom_field = pa.field("geom", pa.binary(), metadata={
  b"ARROW:extension:name":     b"queryfarm.geosilo",
  b"ARROW:extension:metadata": b'{"crs":"EPSG:4326"}',
})

Sibling extensions in the geo bucket

• a5 — pentagonal global geospatial index for spatial aggregation, joins, and equal-area binning. Pairs well: index with a5, store the underlying geometry with geosilo.
• lindel — Hilbert / Morton space-filling curves for ordering multi-dimensional data on disk. Sort a GEOSILO table by a lindel-encoded centroid to keep nearby geometries on the same Parquet row group / DuckDB block.

Cookbook setup

Load spatial first, then geosilo:

INSTALL spatial; LOAD spatial;
INSTALL geosilo FROM community; LOAD geosilo;

The GEOSILO column type and the native ST_* overloads only register once both are loaded — geosilo overlays spatial’s ST_* namespace.

Define a compact geometry table

The recommended shape — explicit CRS, ZSTD compression:

CREATE TABLE parcels (
  id    INTEGER,
  name  VARCHAR,
  geom  GEOSILO('EPSG:4326') USING COMPRESSION zstd
);

The CRS lets GeoSilo pick a sensible integer scale automatically (1e7 for degrees, 100 for projected meters). USING COMPRESSION zstd is where the storage win comes from — without it, raw GEOSILO blobs are larger than DuckDB’s columnar GEOMETRY storage. See CREATE TABLE for the DDL surface.

Load from any GEOMETRY source

The implicit GEOMETRY → GEOSILO cast auto-encodes on insert:

-- From another DuckDB table
INSERT INTO parcels
SELECT id, name, geom FROM raw_parcels;

-- From WKT in CSV
INSERT INTO parcels
SELECT id, name, ST_GeomFromText(wkt) FROM read_csv('parcels.csv');

-- From Parquet with a WKB column
INSERT INTO parcels
SELECT id, name, ST_GeomFromWKB(geom_wkb)
FROM read_parquet('parcels.parquet');

Sources can be any expression returning GEOMETRY — see spatial’s constructor functions for the full list.

Native bounding-box prefilter

The four ST_XMin / ST_XMax / ST_YMin / ST_YMax calls have GEOSILO-native implementations — they walk the integer delta stream without rebuilding a GEOMETRY. Use them in WHERE clauses to prefilter before any GEOS-backed work:

-- Everything in a longitude band — entirely native
SELECT id, name
FROM parcels
WHERE ST_XMin(geom) > -78
  AND ST_XMax(geom) < -76;

-- Hand off to GEOS only for the survivors of the prefilter
SELECT id, ST_AsText(ST_Centroid(geom))
FROM parcels
WHERE ST_XMin(geom) > -78
  AND ST_XMax(geom) < -76
  AND ST_YMin(geom) >  39
  AND ST_YMax(geom) <  41;

See ST_XMin / ST_XMax / ST_YMin / ST_YMax.

Native area, length, perimeter

ST_Area on a GEOSILO column runs the shoelace formula on int64 coordinates rebuilt from the delta stream — no WKB decode:

-- Per-feature area
SELECT id, ST_Area(geom) AS area
FROM parcels
ORDER BY area DESC
LIMIT 10;

-- Total length of a road network
SELECT SUM(ST_Length(geom)) AS total_road_length
FROM roads;

-- Perimeter histogram
SELECT FLOOR(LOG(ST_Perimeter(geom))) AS bucket, COUNT(*) AS n
FROM parcels
GROUP BY bucket
ORDER BY bucket;

See ST_Area, ST_Length, ST_Perimeter.

Per-point coordinate access

ST_X / ST_Y on point geometries decode straight from the silo header:

SELECT id,
       ST_X(geom) AS lon,
       ST_Y(geom) AS lat
FROM stations;

For non-point geometries, use ST_Centroid (which auto-decodes to GEOMETRY) and then ST_X / ST_Y on the resulting point.

Filter by geometry type without decoding

geosilo_metadata reads the silo header directly. Combine with ST_GeometryType (which is also native) to filter cheaply by type:

-- Counts by geometry type — no WKB decode per row
SELECT ST_GeometryType(geom) AS gtype, COUNT(*) AS n
FROM mixed_geom_table
GROUP BY gtype
ORDER BY n DESC;

-- Inspect the integer scale of stored blobs
SELECT geosilo_metadata(silo_blob) AS meta
FROM compact_table
LIMIT 5;
-- {'geometry_type':'POLYGON','vertex_type':'XY','scale':10000000}

Manual encode / decode for BLOB pipelines

When you control the BLOB pipeline directly — exporting to a file format, shipping over a non-Arrow protocol, staging — use geosilo_encode and geosilo_decode:

-- Encode for export
COPY (
  SELECT id, geosilo_encode(geom) AS silo
  FROM parcels_geometry
) TO 'parcels_silo.parquet' (FORMAT 'parquet');

-- Decode on read
SELECT id, ST_AsText(geosilo_decode(silo)) AS wkt
FROM read_parquet('parcels_silo.parquet')
LIMIT 5;

For a custom integer scale (e.g. millimeter-precision projected data), pass the second argument:

SELECT geosilo_encode(geom, 1000) FROM mm_precision_data;

Mix native fast paths with GEOS-backed operations

Anything outside the native overload set auto-decodes through the implicit GEOSILO → GEOMETRY cast. You don’t need to think about it — but if you’re scanning many rows, push the native operations first to shrink the candidate set:

-- ST_Buffer auto-decodes; ST_XMin/XMax/YMin/YMax run native first
SELECT id, ST_AsText(ST_Buffer(geom, 0.001)) AS buffered_wkt
FROM parcels
WHERE ST_XMin(geom) > -78
  AND ST_XMax(geom) < -76
LIMIT 100;

See the spatial extension for the GEOS-backed surface that auto-decodes.

Pair with sibling geo extensions

GeoSilo handles compact storage; pair it with the other geo bucket extensions for indexing and ordering:

-- a5: pentagonal geospatial index for aggregation / joins
INSTALL a5 FROM community; LOAD a5;

-- Bin parcels by A5 cell at resolution 12
SELECT a5_lonlat_to_cell(ST_X(ST_Centroid(geom)),
                         ST_Y(ST_Centroid(geom)),
                         12) AS cell,
       COUNT(*) AS parcels_in_cell
FROM parcels
GROUP BY cell
ORDER BY parcels_in_cell DESC;

See a5 for the indexing surface, and lindel for space-filling curve ordering — sort a GEOSILO table by a Hilbert-encoded centroid to keep nearby geometries on the same DuckDB block / Parquet row group.