Purpose: correct-and-fast spatial SQL. The two recurring failure modes are\nsemantic (geometry vs geography, SRID mismatches → wrong answers) and\nperformance (missing index usage → hour-long joins); this skill guards\nboth.
\nMove from files/GeoPandas to PostGIS when any of: features > a few\nmillion, concurrent readers/writers, repeated ad-hoc querying, a serving\nAPI on top, or transactional integrity needs. For single-shot analytical\nscans over GeoParquet, DuckDB Spatial is often the fastest\nzero-install path — same SQL mindset, no server.
\nWhen requirements are incomplete, do not turn this heuristic into a final\nrecommendation. First obtain current and forecast data volume, concurrency,\ndelivery and mutation pattern, latency/SLA, serving needs, and operational\nownership (including backup and recovery). Define representative ingestion,\njoin, and read queries for both viable backends; compare runtime and resource\nuse only after row counts, join cardinality, SRID, geometry validity, and sample\noutputs agree. Include this benchmark and correctness plan in the current\nresponse; do not merely offer to draft it later.
\nThis runnable example assumes the data is contained in UTM zone 33N. Replace\nEPSG:32633 with a projected CRS verified for the actual area of interest.
\nCREATE TABLE parcels (\n id bigint GENERATED ALWAYS AS IDENTITY PRIMARY KEY,\n parcel_no text NOT NULL,\n landuse text,\n area_m2 double precision, -- unit in the name, always\n geom geometry(MultiPolygon, 32633) NOT NULL\n);\nCREATE INDEX parcels_geom_gix ON parcels USING gist (geom);\nANALYZE parcels;\nType the geometry column fully: geometry(MultiPolygon, SRID) — an\nuntyped geometry column happily accepts mixed garbage.
Promote to Multi* on load (ST_Multi) so Polygon/MultiPolygon mixing\nnever bites.
geometry vs geography: geometry in a projected SRID for regional\nanalysis (fast, full function set); geography (SRID 4326) when the\nextent is global/cross-zone and you want meters without picking a\nprojection (slower, smaller function set). Never store in 4326 geometry\nand call ST_Area expecting m² — that's square degrees.
Never use EPSG:3857/Web Mercator for area or length measurement. When the\nanalysis CRS is not yet known, either use 4326 geography for a geodesic\nresult or stop and select a verified local/equal-area CRS; do not present a\nknown-distorting CRS as a runnable measurement alternative.
\nAny stored geometry column you recommend must be typed with its SRID.\nAdvising a \"second projected geometry column\" for repeated measurement is\nincomplete until it is written as geometry(<Type>, <SRID>) with the index\nand the populating ST_Transform. An untyped column recommended as a fix\nreintroduces the mixed-SRID problem it was meant to solve:
ALTER TABLE parcels ADD COLUMN geom_32633 geometry(MultiPolygon, 32633);\nUPDATE parcels SET geom_32633 = ST_Transform(geom, 32633);\nCREATE INDEX parcels_geom_32633_gix ON parcels USING gist (geom_32633);\nGiST index on every geometry column, ANALYZE after bulk loads; BRIN\nonly for huge, spatially-ordered, append-only tables.
Load paths: ogr2ogr -f PostgreSQL, shp2pgsql, or GeoPandas\nto_postgis (small/medium). COPY beats INSERT by orders of magnitude.
ST_Intersects for \"touches at all\", ST_Contains/ST_Within for\ncontainment, ST_DWithin(a, b, dist) for proximity — never\nST_Distance(a,b) < dist (that form can't use the index).SELECT p.id, a.district\nFROM points p\nJOIN admin a ON ST_Intersects(a.geom, p.geom); -- GiST on both sides\nSELECT h.id, h.name\nFROM hospitals h\nORDER BY h.geom <-> (SELECT geom FROM incident WHERE id = 42)\nLIMIT 3;\n<-> gives true-distance ordering on modern PostGIS for geometry; wrap\nwith ST_DWithin to bound the search when tables are huge.
EXPLAIN (ANALYZE, BUFFERS) first — confirm the GiST index is used\n(look for \"Index Scan ... _gix\"); a Seq Scan on a big spatial join\nmeans a rewrite, not a bigger server.ST_Transform inside a\njoin predicate kills index use; store a transformed, indexed copy\ninstead.ST_Subdivide into a work table (typical 10-100× speedup on\njoins against them).The following example assumes countries(country_id, geom).
CREATE TABLE country_parts AS\nSELECT c.country_id, part.geom\nFROM countries AS c\nCROSS JOIN LATERAL ST_Subdivide(c.geom, 256) AS part(geom);\n\nCREATE INDEX country_parts_geom_gix ON country_parts USING gist (geom);\nANALYZE country_parts;\nST_Subdivide is a set-returning function; do not access its result as\n(ST_Subdivide(...)).geom.
ST_IsValid audit, ST_MakeValid repair, add a\nCHECK (ST_IsValid(geom)) if writers are untrusted.ST_SimplifyPreserveTopology copies or vector tiles\n(ST_AsMVT) for the web tier.VACUUM ANALYZE after churn.-- Area-weighted aggregation (e.g., population into custom zones)\nSELECT z.zone_id,\n SUM(b.pop * ST_Area(ST_Intersection(z.geom, b.geom)) / ST_Area(b.geom)) AS pop_est\nFROM zones z JOIN blocks b ON ST_Intersects(z.geom, b.geom)\nGROUP BY z.zone_id;\n\n-- Dissolve with attribute\nSELECT landuse, ST_Multi(ST_Union(geom))::geometry(MultiPolygon, 32633) AS geom\nFROM parcels GROUP BY landuse;\nArea-weighted interpolation assumes uniform density within source units —\nstate that assumption when reporting. Validity repair is ST_MakeValid,\nnever ST_Buffer(geom, 0).
INSTALL spatial; LOAD spatial;\nSELECT a.name, count(*)\nFROM 'admin.parquet' a, 'points.parquet' p\nWHERE ST_Intersects(a.geom, p.geom)\nGROUP BY a.name;\nReads GeoParquet/Shapefile/GPKG directly, parallel by default — ideal for\none-off large joins and pipeline steps without a server. No GiST; it plans\nits own joins — benchmark, don't assume.
\nSELECT DISTINCT ST_SRID(geom), GeometryType(geom) on every table\ntouched — one query kills two classic bug families.sql fence presented as runnable as a syntax and alias\nboundary: it must execute top-to-bottom after stated schema assumptions.\nNever put angle-bracket placeholders, ellipses, pseudocode, abandoned joins,\nor incomplete aliases inside it. If a schema value such as an SRID is\nunknown, ask for it or keep the template in a labeled text block.ST_Area/ST_Length on 4326 geometry (square degrees).ST_Distance < x instead of ST_DWithin (no index).ST_Transform in join predicates.ST_Subdivide.buffer(0) as validity repair (silent part loss) — ST_MakeValid.