Validating Geometries in DuckDB Spatial

DuckDB removed the most tedious step in file-based spatial QC: the import. With the spatial extension loaded it reads GeoParquet, GeoPackage, FlatGeobuf and shapefiles directly, exposes a GEOS-backed predicate set, and runs the whole check as SQL over the files as they sit. For a lakehouse-shaped pipeline that means a nightly validation job with no server, no loading window and no second copy of the data. This guide runs that job end to end — per-geometry checks, prefiltered topology joins, CRS assertions and a findings table written back as Parquet — implementing the in-place approach set out in Warehouse and Lakehouse Spatial Validation.

Prerequisites

• DuckDB 0.10+ with the spatial extension. The CLI, the Python API (duckdb 0.10+) or any client works; the SQL below is identical in all of them.
• GeoParquet or another readable spatial format on local disk or object storage. For S3 or GCS, the httpfs extension and credentials in the session.
• A declared canonical CRS and the expectation that the data is already in it. DuckDB does not reproject as part of a predicate.
• Roughly 4 GB of memory for comfortable work on a few million features; set memory_limit explicitly rather than relying on the default in a container.

Step-by-Step Procedure

Step 1 — Install the extension and pin what you are running

-- setup.sql
INSTALL spatial;
LOAD spatial;

-- Record what the checks were run against; GEOS version changes validity behaviour.
SELECT version() AS duckdb_version,
       (SELECT extension_version FROM duckdb_extensions() WHERE extension_name = 'spatial')
         AS spatial_version;

SET memory_limit = '6GB';
SET threads = 8;
SET preserve_insertion_order = false;   -- lets the reader parallelise freely

Verification: store both version strings with the run output. When a geometry validates in DuckDB and fails in PostGIS, the first question is always which GEOS each one used, and recording it turns a day of confusion into a one-line answer.

Step 2 — Read the dataset in place

-- Direct read: no import, no copy, predicates pushed into the Parquet scan.
CREATE OR REPLACE VIEW parcels AS
SELECT *
FROM read_parquet('s3://lake/parcels/load_date=2026-08-06/*.parquet');

-- For formats other than Parquet, ST_Read wraps the GDAL drivers.
CREATE OR REPLACE VIEW zoning AS
SELECT * FROM ST_Read('/data/zoning.gpkg', layer = 'zoning_2026');

-- Confirm what the reader thinks the geometry column is.
DESCRIBE parcels;

Verification: DESCRIBE should show the geometry column typed as GEOMETRY, not BLOB. A BLOB means the file is plain Parquet with WKB rather than GeoParquet, and every geometry function will need an explicit ST_GeomFromWKB wrapper — worth knowing before writing twenty rules.

Step 3 — Per-geometry checks in one pass

-- Every per-row rule evaluated in a single scan, emitting the shared contract.
CREATE OR REPLACE TABLE findings_geometry AS
WITH checked AS (
  SELECT
    parcel_id,
    geom,
    ST_IsValid(geom)                       AS is_valid,
    ST_IsValidReason(geom)                 AS reason,
    ST_GeometryType(geom)                  AS gtype,
    ST_NPoints(geom)                       AS vertices,
    ST_Area(geom)                          AS area_m2,
    ST_XMin(geom) AS xmin, ST_XMax(geom) AS xmax,
    ST_YMin(geom) AS ymin, ST_YMax(geom) AS ymax
  FROM parcels
)
SELECT parcel_id AS feature_id, 'GEOM_VALID_001' AS rule_id,
       'blocker' AS severity, reason AS message
FROM checked WHERE NOT is_valid
UNION ALL
SELECT parcel_id, 'GEOM_TYPE_001', 'blocker',
       'expected POLYGON or MULTIPOLYGON, found ' || gtype
FROM checked WHERE gtype NOT IN ('POLYGON', 'MULTIPOLYGON')
UNION ALL
SELECT parcel_id, 'GEOM_EMPTY_001', 'blocker', 'zero-area geometry'
FROM checked WHERE area_m2 <= 0
UNION ALL
SELECT parcel_id, 'GEOM_BOUNDS_001', 'blocker',
       'outside the jurisdiction envelope'
FROM checked
WHERE xmin < 400000 OR xmax > 560000 OR ymin < 100000 OR ymax > 220000
UNION ALL
SELECT parcel_id, 'GEOM_VERTEX_001', 'informational',
       'unusually dense geometry: ' || vertices || ' vertices'
FROM checked WHERE vertices > 10000;

Verification: SELECT rule_id, count(*) FROM findings_geometry GROUP BY 1 should return a plausible mix. A rule with zero hits on a real dataset is either genuinely clean or wrong — check it against a known-bad feature before believing it.

Step 4 — Pairwise topology with a bounding-box prefilter

-- Overlap detection: prefilter on bounding boxes, then apply the exact predicate.
CREATE OR REPLACE TABLE findings_topology AS
WITH boxed AS (
  SELECT parcel_id, geom,
         ST_XMin(geom) AS xmin, ST_XMax(geom) AS xmax,
         ST_YMin(geom) AS ymin, ST_YMax(geom) AS ymax
  FROM parcels
  WHERE ST_IsValid(geom)                 -- never run topology over invalid input
),
candidates AS (
  SELECT a.parcel_id AS a_id, b.parcel_id AS b_id, a.geom AS a_geom, b.geom AS b_geom
  FROM boxed a
  JOIN boxed b
    ON a.parcel_id < b.parcel_id          -- each unordered pair exactly once
   AND a.xmax >= b.xmin AND b.xmax >= a.xmin
   AND a.ymax >= b.ymin AND b.ymax >= a.ymin
)
SELECT a_id AS feature_id, 'TOPO_OVERLAP_001' AS rule_id, 'blocker' AS severity,
       'overlaps ' || b_id || ' by ' ||
         printf('%.2f', ST_Area(ST_Intersection(a_geom, b_geom))) || ' m2' AS message
FROM candidates
WHERE ST_Intersects(a_geom, b_geom)
  AND NOT ST_Touches(a_geom, b_geom)
  AND ST_Area(ST_Intersection(a_geom, b_geom)) > 0.5;

Verification: run EXPLAIN ANALYZE on the candidate CTE and check the row count reaching the exact predicate. On a well-behaved parcel layer it should be a small multiple of the feature count, not its square. If it is not, the bounding-box join is not being applied — usually because the box columns were computed inside the join rather than materialised first.

Step 5 — Assert the CRS and write the findings back

-- The CRS is metadata, so check it explicitly: DuckDB will not do it for you.
CREATE OR REPLACE TABLE findings_crs AS
SELECT 'LAYER' AS feature_id, 'CRS_001' AS rule_id, 'blocker' AS severity,
       'GeoParquet CRS is ' || COALESCE(crs_authority, 'undeclared') ||
       ', expected EPSG:27700' AS message
FROM (
  SELECT json_extract_string(
           parquet_kv_metadata('s3://lake/parcels/load_date=2026-08-06/part-0.parquet'),
           '$.geo.columns.geometry.crs.id.code') AS crs_authority
)
WHERE crs_authority IS DISTINCT FROM '27700';

-- One findings table, one contract, written back beside the data.
COPY (
  SELECT * FROM findings_geometry
  UNION ALL SELECT * FROM findings_topology
  UNION ALL SELECT * FROM findings_crs
) TO 's3://lake/qa/findings/load_date=2026-08-06/findings.parquet'
  (FORMAT PARQUET, COMPRESSION ZSTD);

Verification: read the findings file back and count by severity. Writing the output as Parquet beside the data keeps the whole loop inside the lake, which means the reporting layer described in Building a Spatial Data Quality Scorecard can read it with the same engine.

Interpreting Results

The empty-glob case deserves emphasis. read_parquet over a path that matches nothing raises an error, but a path that matches an empty partition returns zero rows perfectly happily, and a validation run over zero rows reports zero defects. Assert the input row count as the first rule of any scheduled run.

Gotchas & Edge Cases

DuckDB treats geometry as planar, always. There is no geography type, so a dataset in EPSG:4326 will produce areas in square degrees and distances in degrees. That is not a bug, but it makes every metric rule meaningless unless the data is projected first — reproject on write, or convert inside the query with ST_Transform and pay the cost per row.

ST_Read goes through GDAL, so driver quirks come with it. Shapefile field truncation, GeoPackage layer selection and encoding surprises behave exactly as they do elsewhere, and the notes in Shapefile vs GeoPackage Schema Enforcement apply unchanged.

Self-joins need an explicit ordering predicate. a.parcel_id < b.parcel_id is what makes each pair appear once. Without it every overlap is reported twice, and the doubled count will be noticed by whoever reads the dashboard rather than by whoever wrote the SQL.

Spatial indexes exist but are not automatic. DuckDB’s R-tree index helps repeated point queries; for a single-pass batch job the bounding-box prefilter is usually faster than building an index first. Measure rather than assume.

Extension versions move quickly. Functions have been added, renamed and had their semantics tightened between minor releases. Pin the DuckDB version in the job image, and treat an extension upgrade as a change that requires re-running the test suite from Testing Spatial Validation Code.

Writing back to the same prefix you read from is a footgun. A findings file written into the data partition will be picked up by the next read_parquet glob. Keep findings in a separate prefix, as in Step 5.

When to Escalate

• Memory exhaustion on a topology join that is already prefiltered means the layer has genuine spatial clustering — thousands of features in one small area. Partition the job by tile and run it per tile, as described in Batch Processing Large Spatial Datasets.
• Persistent disagreement with PostGIS on validity after version alignment is worth reporting upstream to GEOS with the offending WKT; it is occasionally a real library bug and always worth documenting internally.
• Repair requirements — DuckDB can detect far more than it can safely repair in place over immutable files. Route repairs to the system that owns the writable copy.
• Concurrent writers or transactional needs are outside DuckDB’s model for this use case; that is the point at which a spatial database earns its operational cost.

Frequently Asked Questions

Related

• Warehouse and Lakehouse Spatial Validation — the in-place validation pattern and its cost model
• PostGIS vs GeoPandas for Validation at Scale — the engine comparison DuckDB now sits inside
• Validating GeoParquet Schemas with PyArrow — the file contract these queries depend on

Back to Warehouse and Lakehouse Spatial Validation