🔧 Fix APIs & tile services CRS & projections Serving Spatial Data · step 19 of 20 Updated

Fixing Coordinates in the Wrong Order in an API Response

Problem statement

The API returns valid GeoJSON, the numbers are in plausible ranges, and every point is in the wrong place:

{"type": "Point", "coordinates": [51.5072, -0.1276]}

That is London written latitude-first. Read as GeoJSON — longitude first — it is a point in the Indian Ocean off Somalia, about 7,500 km away.

Nothing raises. The document is well-formed, the values are within range, and every consumer that follows the specification puts the point in the sea. The bug is a convention mismatch, and it enters through one of four doors:

  • EPSG:4326's official axis order is latitude first, and GeoJSON's is longitude first
  • ST_Transform in DuckDB follows the authority order unless told otherwise
  • PostGIS ST_X / ST_Y are unambiguous, but a hand-built coordinate array can transpose them
  • A client library that takes (lat, lon) — Leaflet does — feeding a serialiser that expects (lon, lat)

Quick answer

Assert a known coordinate. One test catches the entire class:

KNOWN = {
    "London":   (-0.1276, 51.5072),
    "New York": (-74.0060, 40.7128),
    "Sydney":   (151.2093, -33.8688),
}


def test_coordinate_order(client):
    for name, (lon, lat) in KNOWN.items():
        feature = client.get(f"/features?name={name}").json()["features"][0]
        got_lon, got_lat = feature["geometry"]["coordinates"][:2]
        assert abs(got_lon - lon) < 0.01, (
            f"{name}: got longitude {got_lon}, expected {lon} — "
            f"the pair is transposed")
        assert abs(got_lat - lat) < 0.01

Sydney is in the set on purpose: a southern-hemisphere, eastern-longitude point catches sign errors that London does not.

Grid of six libraries and specifications with the coordinate order each uses.
Every one is right within its own frame, which is what makes the mismatch so easy.

Step-by-step solution

1. Establish which convention each layer uses 

GeoJSON specification        [longitude, latitude]     always
GeoJSON "CRS84"              [longitude, latitude]     the OGC name for it
EPSG:4326 (authority order)  [latitude, longitude]     what the registry says
Shapely / GEOS               (x, y) = (lon, lat)
PostGIS ST_X / ST_Y          X = longitude
Leaflet                      L.latLng(lat, lon)
MapLibre / Mapbox GL         [lon, lat]

The two rows that cause the damage are the second and third: OGC named the longitude-first variant CRS84 precisely because "EPSG:4326" is ambiguous in practice.

2. Test the transform, not the data

In DuckDB, ST_Transform follows the authority axis order by default:

select st_astext(st_transform(st_point(-0.1276, 51.5072),
                              'EPSG:4326', 'EPSG:3857'));
POINT (5733755.276187301 -14204.378766821068)          -- wrong

select st_astext(st_transform(st_point(-0.1276, 51.5072),
                              'EPSG:4326', 'EPSG:3857', always_xy := true));
POINT (-14204.367025221705 6711506.705400525)          -- correct

pyproj and GeoPandas default to always_xy=True, so the same transform in two tools on one machine gives two answers. If a pipeline crosses between them, this is where to look first.

3. Check the range, which catches most cases free

A latitude outside ±90 is impossible, so a transposed pair is detectable whenever the longitude exceeds 90:

def looks_transposed(lon, lat):
    if abs(lat) > 90:
        return True, "latitude out of range — certainly transposed"
    if abs(lon) <= 90 and abs(lat) <= 90:
        return None, "ambiguous — both values are valid latitudes"
    return False, "plausible"

The ambiguous case is real and common: anywhere within 90° of the meridian, both orderings produce valid-looking coordinates. That is why a known-location assertion beats a range check.

4. Check where the points land

The strongest test uses data you already have. Join the API's output to a country layer and count what falls outside the expected country:

import geopandas as gpd


def landing_check(points, expected_country, countries):
    joined = gpd.sjoin(points.to_crs(countries.crs),
                       countries[["name", "geometry"]],
                       how="left", predicate="within")
    wrong = joined[joined["name"] != expected_country]
    print(f"{len(wrong):,} of {len(points):,} points are not in {expected_country}")
    if len(wrong) > len(points) * 0.5:
        print("  ! more than half — this is a systematic transposition, "
              "not scattered errors")
    return wrong

Systematic is the tell. A handful of points in the wrong country is a geocoding problem; nearly all of them is an axis-order problem.

5. Fix it at the boundary, not in the middle

The right place to normalise is where data enters and leaves the service, in one function each:

def to_geojson_coords(x, y, source_order="xy"):
    """Everything leaving this service is [lon, lat]. One place, one rule."""
    return [x, y] if source_order == "xy" else [y, x]

Fixing it per endpoint guarantees that a new endpoint written next month will get it wrong again.

6. Say which convention you use, in the response

For an OGC API Features service the CRS is declared explicitly:

{"crs": ["http://www.opengis.net/def/crs/OGC/1.3/CRS84"]}

CRS84 means longitude first, unambiguously. Declaring it costs one field and removes the question from every client integration.

Four detection steps for transposed coordinates, from a range check to a spatial join.
A range check passes on transposed data for a large fraction of the world.

Code examples

Example 1 — a validator for a feature collection 

KNOWN_BOUNDS = {"gb": (-8.6, 49.9, 1.8, 60.9), "fr": (-5.2, 41.3, 9.6, 51.1)}


def validate_coordinate_order(collection, expected_region=None, sample=200):
    """Three checks, cheapest first."""
    features = collection["features"][:sample]
    problems = []

    impossible = [f for f in features
                  if abs(coords_of(f)[1]) > 90]
    if impossible:
        problems.append(f"{len(impossible)} features have |latitude| > 90 — "
                        f"certainly transposed")

    if expected_region:
        x0, y0, x1, y1 = KNOWN_BOUNDS[expected_region]
        inside = sum(1 for f in features
                     if x0 <= coords_of(f)[0] <= x1 and y0 <= coords_of(f)[1] <= y1)
        swapped_inside = sum(1 for f in features
                             if x0 <= coords_of(f)[1] <= x1
                             and y0 <= coords_of(f)[0] <= y1)
        if swapped_inside > inside:
            problems.append(
                f"{swapped_inside} features fall inside {expected_region} when "
                f"swapped, against {inside} as given — the pair is transposed")

    for problem in problems:
        print("  !", problem)
    return problems


def coords_of(feature):
    geometry = feature["geometry"]
    coordinates = geometry["coordinates"]
    while isinstance(coordinates[0], (list, tuple)):
        coordinates = coordinates[0]
    return coordinates[0], coordinates[1]

The swapped-inside comparison is the decisive test: if more points land in the expected region after swapping than before, the answer is not in doubt.

Example 2 — normalising at the service boundary 

from dataclasses import dataclass


@dataclass(frozen=True)
class CoordinateConvention:
    name: str
    lon_first: bool


CRS84 = CoordinateConvention("CRS84", lon_first=True)
EPSG4326_AUTHORITY = CoordinateConvention("EPSG:4326 (authority)", lon_first=False)


class Boundary:
    """One place where coordinates enter and one where they leave."""

    def __init__(self, inbound=CRS84, outbound=CRS84):
        self.inbound, self.outbound = inbound, outbound

    def parse_bbox(self, raw: str) -> tuple[float, float, float, float]:
        a, b, c, d = (float(v) for v in raw.split(","))
        if self.inbound.lon_first:
            return a, b, c, d
        return b, a, d, c

    def emit_point(self, x: float, y: float) -> list[float]:
        return [x, y] if self.outbound.lon_first else [y, x]

Example 3 — a regression test against real geometry 

import geopandas as gpd
from shapely.geometry import shape


def test_features_land_in_the_right_country(client, countries_gpkg):
    countries = gpd.read_file(countries_gpkg)[["name", "geometry"]]
    body = client.get("/features?limit=200").json()

    points = gpd.GeoDataFrame(
        {"name": [f["properties"].get("name") for f in body["features"]]},
        geometry=[shape(f["geometry"]).representative_point()
                  for f in body["features"]],
        crs="EPSG:4326")

    joined = gpd.sjoin(points, countries, how="left", predicate="within")
    unmatched = joined["name_right"].isna().sum()
    assert unmatched < len(points) * 0.1, (
        f"{unmatched} of {len(points)} features are in no country — "
        f"coordinates are probably transposed")

Falling in no country is the signature: transposed European coordinates land in the Indian Ocean, and the ocean belongs to nobody.

Explanation

Why the ambiguity exists

The EPSG registry defines coordinate systems including their axis order, and for EPSG:4326 that order is latitude then longitude, because that is how geodesy writes them.

GIS software, GeoJSON, most web APIs and every (x, y) convention put longitude first, because that is how a Cartesian plane works and how a map is drawn.

Both are correct within their own frame. PROJ supports both, which is why always_xy exists, and OGC coined CRS84 as the unambiguous name for the longitude-first variant.

Why the error is thousands of kilometres rather than a small offset

Transposing London’s coordinates gives latitude −0.1276 and longitude 51.5072 — a point in the Indian Ocean, 7,484 km away. The distance is enormous because the two values are unrelated numbers, not a small perturbation of each other.

That is the one mercy in this bug: it is instantly visible on any map. The dangerous variant is a dataset near the diagonal — around 45°N, 45°E — where the transposed point is still on land and merely in the wrong country.

Why range checks are not enough

Anywhere within 90° of the prime meridian, both orderings produce a valid latitude and a valid longitude. Most of Europe, Africa and the Atlantic fall in that band, so a range check passes on transposed data for a large fraction of the world.

The check that always works is a known location, and the check that scales is a spatial join against a country layer: transposed European coordinates land in the ocean, and the ocean is not in any country.

Why fixing it at the boundary is the only durable fix

Transposition bugs recur because the conversion is scattered: one endpoint builds a coordinate array by hand, another passes through a transform with different defaults, a third uses a client library with the opposite convention.

A single function on the way in and a single function on the way out gives one place to be right, one place to test, and a new endpoint that inherits the correct behaviour for free.

Flow showing coordinate normalisation at the inbound and outbound service boundaries.
Declaring CRS84 costs one field and removes the question from every integration.

Edge cases or notes

  • always_xy := true on every DuckDB ST_Transform. Its default follows the authority order.
  • pyproj and GeoPandas default to always_xy=True, so cross-tool pipelines disagree.
  • Leaflet takes (lat, lon); MapLibre takes [lon, lat]. Both are correct in their own API.
  • Declare CRS84, not "EPSG:4326", in an API's metadata.
  • A bbox is minx,miny,maxx,maxy — also longitude first.
  • Near 45°N 45°E a transposed point stays on land, which is the hard case.
  • WKT and WKB are X then Y, so they are unambiguous in practice.
  • Test with a southern, eastern point such as Sydney to catch sign errors.

FAQ

Which order does GeoJSON use?

Longitude first: [lon, lat]. That is the specification, and it is what OGC calls CRS84.

Why does EPSG:4326 mean latitude first?

Because the EPSG registry defines it that way, following geodetic convention. GIS software uses longitude-first almost universally, which is the whole source of the ambiguity.

How far wrong is a transposed coordinate?

For London, about 7,500 km — the point lands in the Indian Ocean. The values are unrelated, so the error is enormous rather than small.

How do I detect it automatically?

Assert a known location, and join a sample of the output to a country layer. If most points fall in no country, or if more fall inside the expected region after swapping, the pair is transposed.

Why does DuckDB transform differently from GeoPandas?

DuckDB's ST_Transform follows the authority axis order by default; pyproj and GeoPandas default to always_xy=True. Pass always_xy := true in DuckDB.

Where should I fix it?

At the service boundary — one function for input, one for output. Fixing it per endpoint guarantees the next endpoint gets it wrong.