Transforming SRIDs in API Responses with ST_Transform

← Back to Coordinate Reference Systems & SRID Handling

This page shows how to let clients request geometry in a coordinate system other than the one you store in, without turning a 7 ms index scan into a 6-second sequential scan.

# Context & When to Use

Storage is one system; consumers want several. A web map wants EPSG:3857 so it can draw without reprojecting, a surveying client wants the national grid its instruments are calibrated to, and everything else wants plain longitude and latitude. Doing the conversion in PostGIS is almost always cheaper than doing it in Python, because the database already has the geometry in memory and PROJ is C.

The rule that makes it safe is positional: transform on output, filter on storage. A ST_Transform call in the SELECT list runs once per returned row, after the planner has already used the GiST index to cut the candidate set. The same call in the WHERE clause runs once per row examined and, worse, hides the indexed column inside an expression the planner cannot match to the index — the comparison laid out in Coordinate Reference Systems & SRID Handling.

Use this approach whenever the output system varies per request. If a single alternative projection accounts for essentially all traffic — a tile service reading 3857 exclusively — a functional index or a materialized projected column is worth measuring instead.

# Runnable Implementation

from typing import Annotated, Any

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

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

STORAGE_SRID = 4326
# Allow-list: every code here must resolve in spatial_ref_sys on this server
OUTPUT_SRIDS = {4326: 6, 3857: 2, 27700: 2, 2154: 2}   # srid -> decimal places

FEATURES_SQL = """
SELECT f.id,
       f.layer,
       ST_AsGeoJSON(
           -- Transform ONLY the rows that survived the filter…
           CASE WHEN $5::int = 4326 THEN f.geom ELSE ST_Transform(f.geom, $5::int) END,
           $6::int                                   -- …at an explicit precision
       )::json AS geometry
FROM   features f
-- …and keep the predicate on the RAW column so features_geom_gix is usable
WHERE  f.geom && ST_MakeEnvelope($1, $2, $3, $4, 4326)
ORDER  BY f.id
LIMIT  $7
"""


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 EPSG:4326")],
    out_srid: Annotated[int, Query(alias="crs")] = STORAGE_SRID,
    limit: Annotated[int, Query(ge=1, le=1000)] = 200,
    pool: asyncpg.Pool = Depends(get_pool),
) -> dict[str, Any]:
    if out_srid not in OUTPUT_SRIDS:
        raise HTTPException(
            422,
            detail={"error": "unsupported_crs", "received": out_srid,
                    "supported": sorted(OUTPUT_SRIDS)},
        )
    try:
        minx, miny, maxx, maxy = (float(v) for v in bbox.split(","))
    except ValueError:
        raise HTTPException(422, detail={"error": "bbox_must_be_four_numbers"})

    async with pool.acquire() as conn:
        rows = await conn.fetch(
            FEATURES_SQL, minx, miny, maxx, maxy,
            out_srid, OUTPUT_SRIDS[out_srid], limit,
        )

    return {
        "type": "FeatureCollection",
        # Tell the client what it received — never leave the system implicit
        "crs": {"type": "name", "properties": {"name": f"EPSG:{out_srid}"}},
        "features": [
            {"type": "Feature", "id": r["id"],
             "geometry": r["geometry"], "properties": {"layer": r["layer"]}}
            for r in rows
        ],
    }

# Key Parameters & Options

Short-circuiting the identity case matters more than it looks: on a service where 80 % of traffic wants 4326, the CASE removes four fifths of all PROJ work for one line of SQL.

The precision argument deserves the same scrutiny. It is the cheapest payload reduction available and the one most often left at the default.

# Gotchas & Failure Modes

Most transform failures are domain-of-validity problems: the geometry is fine, but it sits outside the area the target system was defined for.

• ERROR: transform: couldn't project point … latitude or longitude exceeded limits — a coordinate outside the target system’s domain of validity, typically a global dataset being pushed into a national grid. Clip to the projection’s bounds first, or reject the request for out-of-area features.
• Web Mercator above 85°. ST_Transform(geom, 3857) on polar data produces coordinates that grow without bound and renderers cannot draw. Filter latitude to ±85.06 when the requested output is 3857.
• Missing grid-shift files. Transforms to national grids silently degrade from centimetre to metre accuracy when proj-data is absent from the image. Verify at startup with a round-trip assertion, and pin the package alongside the PostGIS version as described in Pinning PostGIS Versions in Production Images.
• Cache keys that ignore the output system. A cached 3857 payload served to a client that asked for 4326 is a silent 20 000 km error. Include the EPSG code in every cache key — see Redis Caching for Spatial Queries.
• Cold PROJ pipeline after a pool restart. The first transform on each backend costs several milliseconds. Warm it in the connection-setup hook: SELECT ST_Transform(ST_SetSRID(ST_MakePoint(0,0),4326), 3857).

# Advertising the systems you support

Clients should not have to discover the allow-list by trial and error. Expose it, and expose which one is the default, so an integrator can negotiate rather than guess:

@router.get("/crs")
async def supported_crs() -> dict[str, object]:
    return {
        "default": STORAGE_SRID,
        "supported": [
            {"epsg": code, "decimals": dp,
             "uri": f"http://www.opengis.net/def/crs/EPSG/0/{code}"}
            for code, dp in sorted(OUTPUT_SRIDS.items())
        ],
    }

Using the OGC CRS URI form alongside the bare EPSG code costs nothing and makes the endpoint legible to OGC API Features clients, which expect that identifier shape. When a system is later added or withdrawn, this endpoint and the crs response member are the two places the change becomes visible to clients — treat a withdrawal as a breaking change and route it through the deprecation process in API Versioning for GIS Endpoints rather than removing the code silently.

# Verification Snippet

Confirm both that the index is still used and that the coordinates actually moved:

EXPLAIN (ANALYZE, BUFFERS)
SELECT ST_AsGeoJSON(ST_Transform(geom, 3857), 2)
FROM   features
WHERE  geom && ST_MakeEnvelope(-0.2, 51.4, 0.0, 51.6, 4326)
LIMIT  200;
-- Expect: Index Scan using features_geom_gix …  (NOT Seq Scan)

# Same feature, two systems: the numbers must differ by orders of magnitude
curl -s "localhost:8000/v1/features?bbox=-0.2,51.4,0,51.6&limit=1&crs=4326" | jq '.features[0].geometry.coordinates'
# [-0.127761, 51.507351]
curl -s "localhost:8000/v1/features?bbox=-0.2,51.4,0,51.6&limit=1&crs=3857" | jq '.features[0].geometry.coordinates'
# [-14222.34, 6711533.9]

# Related

• Coordinate Reference Systems & SRID Handling — choosing the storage system this page projects out of
• GeoJSON vs GeoParquet Serialization — how output precision interacts with format choice
• Reading EXPLAIN ANALYZE for Spatial Query Optimization — confirming the index survived

← Back to Coordinate Reference Systems & SRID Handling