Measuring Distance and Area in Metres with Geography

← Back to Coordinate Reference Systems & SRID Handling

This page shows how to make every distance, radius and area an endpoint returns a real measurement in metres, instead of a degree value that happens to look like one.

# Context & When to Use

Any endpoint that accepts a radius, returns a distance, sorts by proximity or reports an area is making a measurement claim. If the underlying column is geometry(…, 4326) and nothing casts it, that claim is expressed in degrees — a unit whose ground length varies from 111 km to nothing depending on where you are and which axis you are moving along. Nothing errors; the API just returns a number that means something different for each user, as set out in Coordinate Reference Systems & SRID Handling.

There are two correct answers, and the choice is about workload rather than accuracy. Casting to geography gives geodesic maths on the ellipsoid: correct everywhere on the globe, no projection to choose, slightly slower per call. Transforming to a projected system gives planar maths in metres: faster for heavy polygon work, accurate only inside that projection’s area of use.

Use geography for point-based proximity — “stores within 5 km”, “nearest ten vehicles”, geofence membership. Use a projected system for bulk area and length computation over complex polygons, especially where the result feeds an aggregation rather than a single response.

# Runnable Implementation

-- One index per access pattern: the geometry index cannot serve a geography predicate
CREATE INDEX features_geom_gix ON features USING GIST (geom);
CREATE INDEX features_geog_gix ON features USING GIST ((geom::geography));

-- Proximity search that returns real metres, index-assisted on both sides
PREPARE nearby (float8, float8, float8, int) AS
SELECT f.id,
       f.name,
       ROUND(ST_Distance(
           f.geom::geography,
           ST_SetSRID(ST_MakePoint($1, $2), 4326)::geography
       )::numeric, 1) AS distance_m
FROM   features f
WHERE  ST_DWithin(                       -- $3 is METRES, because both sides are geography
           f.geom::geography,
           ST_SetSRID(ST_MakePoint($1, $2), 4326)::geography,
           $3
       )
ORDER  BY f.geom::geography <-> ST_SetSRID(ST_MakePoint($1, $2), 4326)::geography
LIMIT  $4;

EXECUTE nearby(-0.1276, 51.5072, 5000, 10);   -- within 5 km of Trafalgar Square

The FastAPI side keeps the unit contract explicit — the parameter is named for its unit, validated, and bounded:

from typing import Annotated, Any

import asyncpg
from fastapi import APIRouter, Depends, Query

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

NEARBY_SQL = """
SELECT f.id, f.name,
       ROUND(ST_Distance(f.geom::geography, p.pt)::numeric, 1) AS distance_m
FROM   features f,
       LATERAL (SELECT ST_SetSRID(ST_MakePoint($1, $2), 4326)::geography AS pt) p
WHERE  ST_DWithin(f.geom::geography, p.pt, $3)
ORDER  BY f.geom::geography <-> p.pt
LIMIT  $4
"""


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


@router.get("")
async def nearby(
    lon: Annotated[float, Query(ge=-180, le=180)],
    lat: Annotated[float, Query(ge=-90, le=90)],
    # Named for the unit: no caller can mistake this for degrees
    radius_m: Annotated[float, Query(gt=0, le=50_000)] = 1_000,
    limit: Annotated[int, Query(ge=1, le=200)] = 20,
    pool: asyncpg.Pool = Depends(get_pool),
) -> dict[str, Any]:
    async with pool.acquire() as conn:
        rows = await conn.fetch(NEARBY_SQL, lon, lat, radius_m, limit)
    return {
        "unit": "metre",
        "radius_m": radius_m,
        "results": [dict(r) for r in rows],
    }

# Key Parameters & Options

ST_DWithin rather than ST_Distance(...) < 5000 is the single most valuable substitution on this page: the former is index-assisted, the latter forces a distance computation for every row in the table.

# Cost of each measurement path

# Choosing between geography and a projection for area

Distance is settled — cast to geography and move on. Area is the case where the projected route genuinely competes, because geodesic area on a complex polygon is expensive and the accuracy advantage only matters over long distances. The deciding factors are how big the polygons are, how many of them there are per request, and whether they all fall inside one projection’s area of use.

A rule that holds up in practice: below a few thousand polygons per request, use geography and stop thinking about it. Above that, or where the polygons carry tens of thousands of vertices each, transform once into a local equal-area system and measure there — the accuracy difference across a single country is fractions of a percent, and the speed difference is an order of magnitude.

# Gotchas & Failure Modes

• A geography predicate with only a geometry index. The plan silently becomes a sequential scan. Confirm with EXPLAIN: the Index Cond should name the _geog_gix index, not appear as a Filter.
• Casting only one side. ST_DWithin(geom, point::geography, 5000) raises function st_dwithin(geometry, geography, integer) does not exist, or worse, silently resolves to the geometry overload if both are castable. Cast both operands explicitly.
• ST_Area(geography) on national-scale polygons. Geodesic area on a 200 000-vertex boundary can take seconds. Transform to an equal-area projection for bulk work and accept the area-of-use limits.
• ST_Buffer on geography. It converts to geometry internally, buffers, and converts back — accurate near the buffer’s centre and increasingly wrong at its edges for large radii. Above about 100 km, buffer in a local projection instead.
• Storing the cast instead of casting. A generated geography column doubles storage for data you already have. The functional index gives the same query performance at index cost only.

# Verification Snippet

-- The plan must show the GEOGRAPHY index, not a filter
EXPLAIN (ANALYZE, BUFFERS)
SELECT id FROM features
WHERE ST_DWithin(geom::geography,
                 ST_SetSRID(ST_MakePoint(-0.1276, 51.5072), 4326)::geography,
                 5000);
-- Index Scan using features_geog_gix on features
--   Index Cond: ((geom)::geography && _st_expand(...))

-- A known distance: Trafalgar Square to St Paul's is ~2.06 km
SELECT ROUND(ST_Distance(
         ST_SetSRID(ST_MakePoint(-0.12776, 51.50735), 4326)::geography,
         ST_SetSRID(ST_MakePoint(-0.09831, 51.51385), 4326)::geography
       )::numeric, 0) AS metres;
--  metres
-- --------
--    2131

curl -s "localhost:8000/v1/nearby?lon=-0.1276&lat=51.5072&radius_m=2000&limit=3" | jq
# {"unit":"metre","radius_m":2000,"results":[{"id":41,"name":"…","distance_m":184.6}, …]}

# Related

• Coordinate Reference Systems & SRID Handling — the storage decision that makes the cast necessary
• Optimizing KNN Queries with the PostGIS Distance Operator — ordering by distance at scale
• Reading EXPLAIN ANALYZE for Spatial Query Optimization — confirming index assistance

← Back to Coordinate Reference Systems & SRID Handling