---

name: postgis-spatial-sql description: >- Invoke whenever spatial SQL or its execution backend is the decision: PostGIS, DuckDB Spatial, SpatiaLite, ST_* functions, recurring spatial joins, concurrent/growing workloads, or large GeoParquet queries. Covers backend selection, schemas, GiST/BRIN indexes, KNN, geometry versus geography, correctness benchmarks, and EXPLAIN optimization. Use PostGIS for managed concurrent services and embedded engines for bounded local analytics when evidence supports that choice. Use geo-data-engineering for acquisition, conversion, and file-based ETL without spatial SQL. license: MIT metadata: author: Muhammed Enes Duran

PostGIS & Spatial SQL

Purpose: correct-and-fast spatial SQL. The two recurring failure modes are semantic (geometry vs geography, SRID mismatches → wrong answers) and performance (missing index usage → hour-long joins); this skill guards both.

When the database is the right tool

Move from files/GeoPandas to PostGIS when any of: features > a few million, concurrent readers/writers, repeated ad-hoc querying, a serving API on top, or transactional integrity needs. For single-shot analytical scans over GeoParquet, DuckDB Spatial is often the fastest zero-install path — same SQL mindset, no server.

When requirements are incomplete, do not turn this heuristic into a final recommendation. First obtain current and forecast data volume, concurrency, delivery and mutation pattern, latency/SLA, serving needs, and operational ownership (including backup and recovery). Define representative ingestion, join, and read queries for both viable backends; compare runtime and resource use only after row counts, join cardinality, SRID, geometry validity, and sample outputs agree. Include this benchmark and correctness plan in the current response; do not merely offer to draft it later.

Schema fundamentals

This runnable example assumes the data is contained in UTM zone 33N. Replace EPSG:32633 with a projected CRS verified for the actual area of interest.

CREATE TABLE parcels (
  id          bigint GENERATED ALWAYS AS IDENTITY PRIMARY KEY,
  parcel_no   text NOT NULL,
  landuse     text,
  area_m2     double precision,          -- unit in the name, always
  geom        geometry(MultiPolygon, 32633) NOT NULL
);
CREATE INDEX parcels_geom_gix ON parcels USING gist (geom);
ANALYZE parcels;

• Type the geometry column fully: geometry(MultiPolygon, SRID) — an untyped geometry column happily accepts mixed garbage.

• Promote to Multi* on load (ST_Multi) so Polygon/MultiPolygon mixing never bites.

• geometry vs geography: geometry in a projected SRID for regional analysis (fast, full function set); geography (SRID 4326) when the extent is global/cross-zone and you want meters without picking a projection (slower, smaller function set). Never store in 4326 geometry and call ST_Area expecting m² — that's square degrees.

• Never use EPSG:3857/Web Mercator for area or length measurement. When the analysis CRS is not yet known, either use 4326 geography for a geodesic result or stop and select a verified local/equal-area CRS; do not present a known-distorting CRS as a runnable measurement alternative.

• Any stored geometry column you recommend must be typed with its SRID. Advising a "second projected geometry column" for repeated measurement is incomplete until it is written as geometry(<Type>, <SRID>) with the index and the populating ST_Transform. An untyped column recommended as a fix reintroduces the mixed-SRID problem it was meant to solve:

  ALTER TABLE parcels ADD COLUMN geom_32633 geometry(MultiPolygon, 32633);
  UPDATE parcels SET geom_32633 = ST_Transform(geom, 32633);
  CREATE INDEX parcels_geom_32633_gix ON parcels USING gist (geom_32633);
  

• GiST index on every geometry column, ANALYZE after bulk loads; BRIN only for huge, spatially-ordered, append-only tables.

• Load paths: ogr2ogr -f PostgreSQL, shp2pgsql, or GeoPandas to_postgis (small/medium). COPY beats INSERT by orders of magnitude.

Correct spatial predicates

• ST_Intersects for "touches at all", ST_Contains/ST_Within for containment, ST_DWithin(a, b, dist) for proximity — never ST_Distance(a,b) < dist (that form can't use the index).
• The classic point-in-polygon join:

SELECT p.id, a.district
FROM points p
JOIN admin a ON ST_Intersects(a.geom, p.geom);   -- GiST on both sides

• KNN nearest-neighbor with the distance operator (index-assisted):

SELECT h.id, h.name
FROM hospitals h
ORDER BY h.geom <-> (SELECT geom FROM incident WHERE id = 42)
LIMIT 3;

<-> gives true-distance ordering on modern PostGIS for geometry; wrap with ST_DWithin to bound the search when tables are huge.

Performance playbook

1. EXPLAIN (ANALYZE, BUFFERS) first — confirm the GiST index is used (look for "Index Scan ... _gix"); a Seq Scan on a big spatial join means a rewrite, not a bigger server.
2. Same SRID on both sides of every predicate — ST_Transform inside a join predicate kills index use; store a transformed, indexed copy instead.
3. Big-polygon problem: country/basin-sized geometries make index bboxes useless → ST_Subdivide into a work table (typical 10-100× speedup on joins against them).

The following example assumes countries(country_id, geom).

CREATE TABLE country_parts AS
SELECT c.country_id, part.geom
FROM countries AS c
CROSS JOIN LATERAL ST_Subdivide(c.geom, 256) AS part(geom);

CREATE INDEX country_parts_geom_gix ON country_parts USING gist (geom);
ANALYZE country_parts;

ST_Subdivide is a set-returning function; do not access its result as (ST_Subdivide(...)).geom.

4. Validity in-database: ST_IsValid audit, ST_MakeValid repair, add a CHECK (ST_IsValid(geom)) if writers are untrusted.
5. Simplify for serving, not for analysis: keep full-resolution geometry; generate ST_SimplifyPreserveTopology copies or vector tiles (ST_AsMVT) for the web tier.
6. Batch updates in transactions; VACUUM ANALYZE after churn.

Common analytical patterns

-- Area-weighted aggregation (e.g., population into custom zones)
SELECT z.zone_id,
       SUM(b.pop * ST_Area(ST_Intersection(z.geom, b.geom)) / ST_Area(b.geom)) AS pop_est
FROM zones z JOIN blocks b ON ST_Intersects(z.geom, b.geom)
GROUP BY z.zone_id;

-- Dissolve with attribute
SELECT landuse, ST_Multi(ST_Union(geom))::geometry(MultiPolygon, 32633) AS geom
FROM parcels GROUP BY landuse;

Area-weighted interpolation assumes uniform density within source units — state that assumption when reporting. Validity repair is ST_MakeValid, never ST_Buffer(geom, 0).

DuckDB Spatial quick path

INSTALL spatial; LOAD spatial;
SELECT a.name, count(*)
FROM 'admin.parquet' a, 'points.parquet' p
WHERE ST_Intersects(a.geom, p.geom)
GROUP BY a.name;

Reads GeoParquet/Shapefile/GPKG directly, parallel by default — ideal for one-off large joins and pipeline steps without a server. No GiST; it plans its own joins — benchmark, don't assume.

Verification protocol

1. Row-count accounting query after each join/overlay CTE.
2. SELECT DISTINCT ST_SRID(geom), GeometryType(geom) on every table touched — one query kills two classic bug families.
3. Sample 5 output features rendered over a basemap (QGIS connects directly) — numbers can pass while geometries are garbage.
4. Treat every sql fence presented as runnable as a syntax and alias boundary: it must execute top-to-bottom after stated schema assumptions. Never put angle-bracket placeholders, ellipses, pseudocode, abandoned joins, or incomplete aliases inside it. If a schema value such as an SRID is unknown, ask for it or keep the template in a labeled text block.

Pitfalls checklist

• ST_Area/ST_Length on 4326 geometry (square degrees).
• EPSG:3857/Web Mercator for area or length measurement (systematic distortion).
• ST_Distance < x instead of ST_DWithin (no index).
• ST_Transform in join predicates.
• Untyped geometry columns with mixed SRIDs.
• Country-sized polygons joined without ST_Subdivide.
• buffer(0) as validity repair (silent part loss) — ST_MakeValid.
• Serving full-resolution geometries to web clients.

Execution contract

• Workflow: inspect schema, SRID, geometry type, size, and query goal; choose predicates and indexes; write auditable CTEs; inspect the plan; reconcile results; operationalize safely.
• Decision rules: use PostGIS for concurrent, repeated, or transactional spatial workloads; use file pipelines or DuckDB Spatial for bounded one-off transformations when a server adds no value.
• Verification protocol: assert SRID and geometry invariants, account for rows at each join, compare indexed plans and timings, sample geometries on a map, and test boundary semantics.
• Failure modes: block release for mixed SRIDs, accidental many-to-many explosion, invalid geometries, non-indexable predicates, geography/geometry unit confusion, or unexplained plan regressions.
• Deliverables: self-contained parameterized SQL or migration with consistent CTE/table aliases, indexes and rationale, query plan evidence, row accounting, sample validation, expected schema, performance notes, and rollback guidance.
• Source freshness: consult the authoritative source registry for the deployed database and extension versions before selecting functions or plans.