Spatial SQL

Most wrong answers from spatial SQL run without error: a distance in degrees, a join that drops boundary features, a query that scans every row. Inspect the tables first, write the query in the order below, and check the result with counts.

Core rules

1. Inspect first: SRID, geometry type, row count, extent, validity, and whether an index exists.
2. Every geometry in one predicate shares one SRID. Measure in meters: ST_Transform to a projected SRID, or cast to geography.
3. Filter with index-friendly predicates (ST_Intersects, ST_DWithin, &&), never ST_Distance(...) < x in a WHERE.
4. Choose the predicate on purpose: ST_Intersects includes boundary touches, ST_Contains and ST_Within exclude points on the boundary.
5. Check: row counts before and after, unmatched rows in a join, one number recomputed another way. Read EXPLAIN when a query is slow.
6. Return the full SQL that was run.

Workflow index

A. Inspect (PostGIS)

SELECT f_table_name, f_geometry_column, srid, type FROM geometry_columns;
SELECT count(*) AS n, count(*) FILTER (WHERE NOT ST_IsValid(geom)) AS invalid,
       count(*) FILTER (WHERE geom IS NULL) AS null_geom,
       ST_Extent(geom) AS extent, min(ST_SRID(geom)) AS srid_min, max(ST_SRID(geom)) AS srid_max
FROM parcels;
SELECT indexname, indexdef FROM pg_indexes WHERE tablename = 'parcels';

Red flags: SRID 0 (unknown), extent in degrees for a "meters" table, mixed SRIDs, no GiST index on a large table.

B. Points in polygons

SELECT z.id, count(p.id) AS n_points
FROM zones z
LEFT JOIN points p ON ST_Intersects(z.geom, p.geom)
GROUP BY z.id;

• LEFT JOIN keeps zones with zero points; an inner join hides them.
• Check: SELECT sum(n_points) against count(*) of points. More means overlapping zones count a point twice; fewer means points outside all zones (report how many).
• If the two tables have different SRIDs, transform the smaller side, ideally once into a new table, not inside the join.

C. Distance and nearest

-- within 500 m, table in a projected SRID in meters (here 2154)
SELECT a.id, b.id FROM wells a JOIN schools b ON ST_DWithin(a.geom, b.geom, 500);

-- table in 4326: geography gives meters on the spheroid
SELECT a.id, b.id FROM wells a JOIN schools b
  ON ST_DWithin(a.geom::geography, b.geom::geography, 500);

-- k nearest neighbors; the GiST index is used because the probe geometry is constant per lateral row
SELECT n.id, n.name, n.d
FROM (SELECT geom FROM wells WHERE id = 42) w
CROSS JOIN LATERAL (
  SELECT s.id, s.name, ST_Distance(s.geom, w.geom) AS d
  FROM schools s ORDER BY s.geom <-> w.geom LIMIT 5) n;

-- nearest school for every well
SELECT w.id, n.id AS school_id, n.d
FROM wells w
CROSS JOIN LATERAL (
  SELECT s.id, ST_Distance(w.geom, s.geom) AS d FROM schools s
  ORDER BY w.geom <-> s.geom LIMIT 1) n;

<-> on SRID 4326 geometry orders by degrees; use a projected copy or geography. ST_DWithin(geom_4326, geom_4326, 500) means 500 degrees: always wrong. A geography cast without a matching index is slow on large tables; prefer a projected copy with its own GiST index.

D. Areas and lengths

SELECT id, ST_Area(ST_Transform(geom, 2154)) / 10000 AS area_ha FROM parcels;   -- projected
SELECT id, ST_Area(geom::geography) AS area_m2 FROM parcels;                    -- 4326, spheroid

Never ST_Area(geom) on SRID 4326 (square degrees) or 3857 (inflated away from the equator). Cross-check a total both ways; the two should agree within about 1 percent for a national CRS.

E. Invalid and heavy geometries

UPDATE parcels SET geom = ST_Multi(ST_CollectionExtract(ST_MakeValid(geom), 3))
WHERE NOT ST_IsValid(geom);

ST_MakeValid can return a GeometryCollection; ST_CollectionExtract(..., 3) keeps polygons (1 points, 2 lines). Very large polygons (countries, coastlines) make every intersect slow: split them once with ST_Subdivide(geom, 256) into a new indexed table and query that.

F. Performance

CREATE INDEX ON parcels USING GIST (geom);
ANALYZE parcels;
EXPLAIN ANALYZE SELECT ... ;

Look for "Index Scan" or "Bitmap Index Scan" on the GiST index. A "Seq Scan" on a large table usually means: no index, a function wrapped around the indexed column (ST_Transform(geom, ...) in the WHERE), a geography cast, or ST_Distance in the WHERE. Index the expression you filter on, or store the transformed geometry.

Before handing over

• SRIDs checked and stated; measures in meters.
• Predicate chosen on purpose (boundary cases).
• Row counts: inputs, output, unmatched.
• One number cross-checked.
• SQL returned to the user.

Going further

The TerraLab MCP server extracts open data (buildings, roads, places) for an area from any AI assistant, ready to load into PostGIS or DuckDB: https://terra-lab.ai/mcp?utm_source=skills&utm_medium=geoai-skills&utm_campaign=spatial-sql