Coordinate Reference Systems & SRID Handling

← Back to Core Geospatial API Architecture

Every spatial API makes one irreversible decision on its first day: which coordinate reference system the geometry column is typed in. Get it right and reprojection is a cheap output-side concern; get it wrong and you inherit distorted distances, mixed-SRID errors that only fire in production, and sequential scans on tables that have a perfectly good index. This page covers the decision itself, the transformation rules that keep queries index-backed, and the validation layer that stops a client’s unlabelled coordinates from silently corrupting a dataset.

The problem is subtle because nothing crashes. A polygon stored in EPSG:3857 but labelled 4326 still draws on a map; it just draws in the wrong place, several hundred kilometres from where the surveyor put it. A distance filter written in degrees still returns rows; it just returns the wrong ones, and the error scales with latitude. These are data-integrity failures wearing the costume of a working feature, which is why the defences belong in the schema and in the Pydantic geometry validators at the edge rather than in a code review checklist.

# Prerequisites & Environment

The examples assume PostgreSQL 15 or 16 with PostGIS 3.3+, proj 9.x, FastAPI 0.110+, SQLAlchemy 2.0, and asyncpg 0.29. Confirm the PROJ database is present before relying on any transform — a container built without proj-data can resolve EPSG:4326 and EPSG:3857 from the built-in tables while failing on national grids such as EPSG:27700 or EPSG:2154:

-- Version and PROJ availability
SELECT postgis_full_version();

-- Does the target system resolve at all?
SELECT srid, auth_name, proj4text
FROM   spatial_ref_sys
WHERE  srid IN (4326, 3857, 27700, 2154);

-- Round-trip sanity check: transform out and back, expect sub-millimetre drift
SELECT ST_Distance(
         ST_SetSRID(ST_MakePoint(-0.1276, 51.5072), 4326)::geography,
         ST_Transform(
           ST_Transform(ST_SetSRID(ST_MakePoint(-0.1276, 51.5072), 4326), 27700),
           4326)::geography
       ) AS round_trip_error_m;

A round-trip error above a millimetre means the grid-shift files are missing and PROJ has fallen back to a coarse seven-parameter approximation. That is a deployment bug, not a data bug — pin the PostGIS image tag and the proj-data package together, as covered in Pinning PostGIS Versions in Production Images.

# Decision Matrix: which system to store in

The storage SRID is a schema decision with API-wide consequences. The table below compares the four candidates that come up in practice for a general-purpose feature API.

For most APIs the answer is geometry(…, 4326) with a cast to geography at the point of measurement. It keeps the stored value identical to what clients send and receive, and confines the projection question to the two or three endpoints that actually measure something.

The tree above resolves the storage question. Everything after it is mechanical: transform on the way out, measure in metres, and never let an unlabelled coordinate reach the table.

# Why degrees are not a unit of length

The single most common spatial bug in a young API is a distance filter written against a 4326 geometry. PostGIS answers the question you asked — the Cartesian distance between two points on a longitude/latitude plane — and that answer is in degrees. Because a degree of longitude shrinks with the cosine of latitude, the same numeric radius covers a wildly different ground distance depending on where the user is standing.

The practical rule: the moment an endpoint accepts a radius, a tolerance, an area threshold or a buffer, that number is in metres and the query must run in a metric system. Use geography for point-to-point work and a projected SRID for heavy polygon overlays. The K-Nearest Neighbor Routing Algorithms topic covers the ordering side of the same problem, where the <-> operator returns degrees for geometry and metres for geography.

# Step-by-Step Implementation

# 1. Pin the storage SRID in the column type

An untyped geometry column accepts anything: a 4326 point, a 3857 point and an SRID-0 point can coexist in the same table, and the failure only surfaces when a query compares two of them.

CREATE TABLE features (
    id          bigserial PRIMARY KEY,
    layer       text NOT NULL,
    -- Typed: PostgreSQL rejects any insert whose SRID is not 4326
    geom        geometry(MultiPolygon, 4326) NOT NULL,
    captured_at timestamptz NOT NULL DEFAULT now()
);

CREATE INDEX features_geom_gix ON features USING GIST (geom);
-- Second index for metric proximity queries against the same column
CREATE INDEX features_geog_gix ON features USING GIST ((geom::geography));

The typed column is the cheapest validation in the stack: it is enforced by the database, applies to every writer including migrations and manual psql sessions, and costs nothing at query time. Auditing an existing table for stragglers takes one query:

SELECT ST_SRID(geom) AS srid, count(*)
FROM   features
GROUP  BY 1
ORDER  BY 2 DESC;

# 2. Accept a declared input CRS and reject the rest

GeoJSON deliberately removed the crs member in RFC 7946 and mandates 4326, but real clients still post state-plane coordinates into a field labelled geometry. Accept an explicit query parameter, validate it against an allow-list, and fail loudly on anything else.

from typing import Annotated
from fastapi import Depends, HTTPException, Query

# Only systems the API is prepared to transform from
SUPPORTED_SRIDS: dict[int, str] = {
    4326: "WGS 84 lon/lat",
    3857: "Web Mercator",
    27700: "OSGB36 / British National Grid",
    2154: "RGF93 / Lambert-93",
}
STORAGE_SRID = 4326


def input_srid(
    crs: Annotated[int, Query(description="EPSG code of the posted geometry")] = STORAGE_SRID,
) -> int:
    if crs not in SUPPORTED_SRIDS:
        raise HTTPException(
            status_code=422,
            detail={
                "error": "unsupported_crs",
                "received": crs,
                "supported": sorted(SUPPORTED_SRIDS),
            },
        )
    return crs

Returning 422 with the supported list turns a silent 300-metre offset into a client-side fix on the first request. Pair it with the geometry-shape validation described in Validating WKT and GeoJSON with Pydantic v2, which rejects malformed rings before the coordinates ever reach PostGIS.

# 3. Transform inbound, once, at the boundary

Normalise on write. Every geometry that lands in the table is already in the storage SRID, so no read path ever has to think about it.

INSERT INTO features (layer, geom)
VALUES (
  $1,
  -- $2 is GeoJSON text, $3 the client's declared EPSG code
  ST_Transform(ST_SetSRID(ST_GeomFromGeoJSON($2), $3), 4326)
)
RETURNING id;

ST_SetSRID labels the coordinates; ST_Transform moves them. Calling ST_SetSRID alone is the classic corruption path — it relabels 3857 metres as degrees, producing a point somewhere past the edge of the map with no error raised.

# 4. Keep the predicate on the raw column

Transform placement decides whether the query uses the index. The predicate must compare the stored column against a constant that has already been transformed into storage space, not the other way round.

The measurements come from a 4.2 million row parcel table on a db.r6g.large-class instance; the ratio, not the absolute number, is the point. Reading these plans is covered in depth in Reading EXPLAIN ANALYZE for Spatial Query Optimization.

If a projected form really is queried on every request — a tile server reading 3857 exclusively — add a functional index instead of moving the transform:

CREATE INDEX features_geom_3857_gix
    ON features USING GIST (ST_Transform(geom, 3857));

That index is only usable when the query expression matches exactly, and it doubles the write cost of the table, so reach for it after measuring rather than before.

# 5. Project on output, next to the serializer

Output projection belongs in the SELECT list, after filtering and pagination have already cut the row count down. Combine it with the serialization decision described in GeoJSON vs GeoParquet Serialization so the coordinate precision matches the format.

SELECT id,
       layer,
       ST_AsGeoJSON(
         CASE WHEN $2::int = 4326 THEN geom
              ELSE ST_Transform(geom, $2::int) END,
         6   -- 6 decimal places ≈ 0.11 m; more is noise for most APIs
       )::json AS geometry
FROM   features
WHERE  ST_Intersects(geom, ST_MakeEnvelope($3, $4, $5, $6, 4326))
ORDER  BY id
LIMIT  $7;

# Production Code Example

A complete FastAPI route that accepts a client CRS for both input bounds and output geometry, keeps the predicate index-backed, and measures in metres.

import json
from typing import Annotated, Any

import asyncpg
from fastapi import APIRouter, Depends, HTTPException, Query

router = APIRouter(prefix="/v1/features", tags=["features"])

STORAGE_SRID = 4326
SUPPORTED_SRIDS = {4326, 3857, 27700, 2154}

FEATURES_SQL = """
SELECT f.id,
       f.layer,
       ST_AsGeoJSON(
         CASE WHEN $6::int = 4326 THEN f.geom ELSE ST_Transform(f.geom, $6::int) END,
         6
       )::json                                        AS geometry,
       -- Metric distance from the viewport centre: cast, never subtract degrees
       ROUND(ST_Distance(
         f.geom::geography,
         ST_Centroid(ST_MakeEnvelope($1, $2, $3, $4, 4326))::geography
       )::numeric, 1)                                 AS distance_m
FROM   features f
-- Predicate compares the RAW column, so features_geom_gix is usable
WHERE  ST_Intersects(f.geom, ST_MakeEnvelope($1, $2, $3, $4, 4326))
ORDER  BY f.id
LIMIT  $5
"""


def validated_srid(code: int, field: str) -> int:
    if code not in SUPPORTED_SRIDS:
        raise HTTPException(
            status_code=422,
            detail={"error": "unsupported_crs", "field": field, "received": code,
                    "supported": sorted(SUPPORTED_SRIDS)},
        )
    return code


async def get_pool() -> asyncpg.Pool:      # wired at app startup
    raise NotImplementedError


@router.get("")
async def list_features(
    bbox: Annotated[str, Query(description="minx,miny,maxx,maxy in bbox_crs")],
    bbox_crs: Annotated[int, Query(description="EPSG code of the bbox")] = 4326,
    out_crs: Annotated[int, Query(description="EPSG code for returned geometry")] = 4326,
    limit: Annotated[int, Query(ge=1, le=1000)] = 200,
    pool: asyncpg.Pool = Depends(get_pool),
) -> dict[str, Any]:
    validated_srid(bbox_crs, "bbox_crs")
    validated_srid(out_crs, "out_crs")

    try:
        minx, miny, maxx, maxy = (float(v) for v in bbox.split(","))
    except ValueError:
        raise HTTPException(422, detail={"error": "bbox_must_be_four_numbers"})
    if minx >= maxx or miny >= maxy:
        raise HTTPException(422, detail={"error": "bbox_min_must_precede_max"})

    async with pool.acquire() as conn:
        if bbox_crs != STORAGE_SRID:
            # Transform the CONSTANT into storage space — one call, not one per row
            minx, miny, maxx, maxy = await conn.fetchrow(
                """
                SELECT ST_XMin(e), ST_YMin(e), ST_XMax(e), ST_YMax(e)
                FROM  (SELECT ST_Transform(
                                ST_MakeEnvelope($1, $2, $3, $4, $5), 4326) AS e) t
                """,
                minx, miny, maxx, maxy, bbox_crs,
            )

        rows = await conn.fetch(FEATURES_SQL, minx, miny, maxx, maxy, limit, out_crs)

    return {
        "type": "FeatureCollection",
        "crs": {"type": "name", "properties": {"name": f"EPSG:{out_crs}"}},
        "features": [
            {
                "type": "Feature",
                "id": r["id"],
                "geometry": json.loads(r["geometry"]),
                "properties": {"layer": r["layer"], "distance_m": float(r["distance_m"])},
            }
            for r in rows
        ],
    }

Two details carry most of the value. The bounding box is transformed once, in a single round trip, before the main query runs — so the predicate still matches the index. And distance_m casts to geography rather than subtracting degrees, so the number means the same thing in Oslo as it does in Nairobi.

# Verification & Testing

Reprojection bugs are invisible to eyeball testing, so test them numerically. Assert against known control points with published coordinates in both systems.

import pytest
from httpx import ASGITransport, AsyncClient

from app.main import app

# Nelson's Column, London — published 4326 and OSGB36 grid coordinates
LON, LAT = -0.12776, 51.50735
EASTING, NORTHING = 530034.0, 180381.0


@pytest.mark.asyncio
async def test_bng_bbox_matches_wgs84_bbox(seeded_db):
    """The same viewport expressed in two systems returns the same features."""
    transport = ASGITransport(app=app)
    async with AsyncClient(transport=transport, base_url="http://t") as client:
        wgs = await client.get("/v1/features", params={
            "bbox": f"{LON - 0.01},{LAT - 0.01},{LON + 0.01},{LAT + 0.01}"})
        bng = await client.get("/v1/features", params={
            "bbox": f"{EASTING - 700},{NORTHING - 1100},{EASTING + 700},{NORTHING + 1100}",
            "bbox_crs": 27700})

    assert wgs.status_code == bng.status_code == 200
    wgs_ids = {f["id"] for f in wgs.json()["features"]}
    bng_ids = {f["id"] for f in bng.json()["features"]}
    # Allow a small edge difference from the non-identical footprints
    assert len(wgs_ids ^ bng_ids) <= 2


@pytest.mark.asyncio
async def test_unknown_crs_is_rejected(seeded_db):
    transport = ASGITransport(app=app)
    async with AsyncClient(transport=transport, base_url="http://t") as client:
        r = await client.get("/v1/features", params={"bbox": "0,0,1,1", "bbox_crs": 99999})
    assert r.status_code == 422
    assert r.json()["detail"]["error"] == "unsupported_crs"

The database side deserves an assertion too. This query fails loudly if any row has drifted out of the declared system:

-- Should return zero rows on a healthy table
SELECT id, ST_SRID(geom)
FROM   features
WHERE  ST_SRID(geom) <> 4326
LIMIT  10;

Wire both checks into the pipeline described in GitHub Actions Integration Tests with a PostGIS Service Container so a reprojection regression fails the build rather than the map.

# Failure Modes & Edge Cases

1. ERROR: Operation on mixed SRID geometries — two operands carry different SRIDs, usually because a constant was built with ST_GeomFromText without the SRID argument. Always pass it: ST_GeomFromText('POINT(0 0)', 4326).
2. ST_SetSRID used where ST_Transform was meant. The coordinates do not move, only the label. Symptom: features land in the Gulf of Guinea (0, 0) or several hundred kilometres off. There is no error and no way to recover the original values once the raw import is gone.
3. Axis order confusion. EPSG:4326 formally defines latitude first, but GeoJSON, PostGIS and every web map use longitude first. A WMS or WFS client following the strict definition sends the pair reversed, producing points in the wrong hemisphere. Validate that latitude is within ±90 and reject the request rather than clamping.
4. Antimeridian-crossing bounding boxes. ST_MakeEnvelope(179, -1, -179, 1, 4326) produces an envelope that wraps the entire globe the wrong way. Split the request into two envelopes at ±180 and union the result sets.
5. geography distance on huge polygons. ST_Distance on geography uses geodesic maths and is markedly slower for complex polygons. For polygon-to-polygon work at national scale, transform both operands to a local equal-area projection instead.
6. Missing grid-shift files. Transformations to national grids silently degrade from centimetre to metre accuracy when proj-data is absent. The round-trip query in the prerequisites section is the canary — run it in a startup health check.
7. Precision inflation. ST_AsGeoJSON(geom) defaults to 15 significant digits, which triples payload size for no benefit. Six decimal places is roughly 11 cm; specify it explicitly.

# Performance Notes

ST_Transform costs roughly 1–3 µs per simple geometry once PROJ has cached the transformation pipeline, and the first call per connection pays an extra 2–10 ms while that pipeline is built. On a pooled connection this warm-up is amortised away, but it does show up as a latency spike after a PgBouncer restart in transaction pooling mode, where every backend is effectively cold.

Casting to geography is free in storage terms — it is the same coordinates with different semantics — but it needs its own GiST index, since a geometry index cannot serve a geography predicate. Budget roughly 15–20 % of the table’s size for the second index and measure whether the proximity endpoint justifies it.

For output projection at volume, transforming 10 000 features in the SELECT list adds around 20–30 ms. That is usually cheaper than transforming client-side, but if the same viewport is requested repeatedly it is cheaper still to cache the projected payload, which is exactly what Redis Caching for Spatial Queries is for — include the output EPSG code in the cache key so a 3857 response never gets served to a 4326 client.

# Related

• Spatial Resource Modeling Patterns — where the geometry column sits in the wider resource design
• GeoJSON vs GeoParquet Serialization — coordinate precision and format selection on output
• Strict Pydantic Validation for Geometry — rejecting malformed and out-of-range coordinates at the edge
• Bounding Box & Spatial Index Queries — the index behaviour that transform placement protects
• API Versioning for GIS Endpoints — how to change a default output system without breaking clients

← Back to Core Geospatial API Architecture