Fixing Coordinates in the Wrong Order in an API Response
On this page
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 some 600 km off the Somali coast, 7,484 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_Transformin DuckDB follows the authority order unless told otherwise- PostGIS
ST_X/ST_Yare 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.1, (
f"{name}: got longitude {got_lon}, expected {lon} — "
f"the pair is transposed")
assert abs(got_lat - lat) < 0.1for label, client in (("good", good), ("transposed", transposed)):
try:
test_coordinate_order(client)
print(f"{label:10} passed")
except AssertionError as exc:
print(f"{label:10} failed: {exc}")good passed
transposed failed: London: got longitude 51.5019406, expected -0.1276 — the pair is transposed
The tolerance is 0.1°, not 0.01°, for a reason found by running it: with 0.01° this test failed the correct API, because Natural Earth places New York 0.0103° east of the value above. Two gazetteers disagree about where a city is by that much; a swap moves it by tens of degrees, so the tolerance only has to sit between the two.
The run also shows why a name is a weak key. Natural Earth has two Sydneys, and before the app sorted its matches by population the test picked up Sydney, Nova Scotia, at longitude −60.18 and reported a transposition that was not there. Query your reference points by something unambiguous — an id, or a name plus a country.
Sydney is in the set on purpose: a southern-hemisphere, eastern-longitude point catches sign errors that London does not.
Know each layer's convention, test the transform
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) -- correctpyproj 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"for name, (lon, lat) in {"London": (51.5072, -0.1276), "Sydney": (-33.8688, 151.2093),
"New York": (40.7128, -74.0060)}.items():
print(f"{name:9} given ({lon}, {lat}):", looks_transposed(lon, lat))London given (51.5072, -0.1276): (None, 'ambiguous — both values are valid latitudes')
Sydney given (-33.8688, 151.2093): (True, 'latitude out of range — certainly transposed')
New York given (40.7128, -74.006): (None, 'ambiguous — both values are valid latitudes')
Only Sydney, whose longitude is over 90, is caught by the range test. London and New York, transposed, are still valid coordinates — which is why the next check asks where the points land.
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 wrongfrom shapely.geometry import shape
def as_points(client):
body = client.get("/features").json()
return gpd.GeoDataFrame({"place": [f["properties"]["name"] for f in body["features"]]},
geometry=[shape(f["geometry"]) for f in body["features"]], crs=4326)
world = countries.rename(columns={"name": "name"})
uk = places[places.ADM0NAME == "United Kingdom"].NAME.tolist()
for label, client in (("good", good), ("transposed", transposed)):
pts = as_points(client)
print(label, end=": ")
landing_check(pts[pts.place.isin(uk)], "United Kingdom", world)good: 0 of 7 points are not in United Kingdom
transposed: 7 of 7 points are not in United Kingdom
! more than half — this is a systematic transposition, not scattered errors
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.
Validating and normalising world places
The data used on this page
Natural Earth's 1:50m populated places and countries (public domain), served by a small FastAPI app so the tests have a real API to call. The app is built twice: good writes each place as [lon, lat], and transposed writes [lat, lon], the way a service does when it serialises EPSG:4326 in its authority order. Both are exercised with FastAPI's TestClient.
import io
import zipfile
from pathlib import Path
import geopandas as gpd
import requests
from fastapi import FastAPI
from fastapi.testclient import TestClient
def natural_earth(name):
if not Path(f"{name}.shp").exists():
zipfile.ZipFile(io.BytesIO(requests.get(
f"https://naciscdn.org/naturalearth/50m/cultural/{name}.zip", timeout=120).content)).extractall(".")
return gpd.read_file(f"{name}.shp")
places = (natural_earth("ne_50m_populated_places")[["NAME", "ADM0NAME", "POP_MAX", "geometry"]]
.sort_values("POP_MAX", ascending=False)) # a name search returns the biggest match first
countries = natural_earth("ne_50m_admin_0_countries")[["NAME", "geometry"]].rename(columns={"NAME": "name"})
countries.to_file("countries.gpkg")
def make_app(lat_first: bool) -> FastAPI:
app = FastAPI()
@app.get("/features")
def features(name: str | None = None, limit: int = 1000):
rows = places if name is None else places[places.NAME == name]
feats = [{"type": "Feature", "properties": {"name": r.NAME, "country": r.ADM0NAME},
"geometry": {"type": "Point",
"coordinates": [r.geometry.y, r.geometry.x] if lat_first
else [r.geometry.x, r.geometry.y]}}
for r in rows.head(limit).itertuples()]
return {"type": "FeatureCollection", "features": feats}
return app
good, transposed = TestClient(make_app(False)), TestClient(make_app(True))
print(len(places), "places,", len(countries), "countries;",
good.get("/features?name=London").json()["features"][0]["geometry"],
transposed.get("/features?name=London").json()["features"][0]["geometry"])1251 places, 242 countries; {'type': 'Point', 'coordinates': [-0.1186677, 51.5019406]} {'type': 'Point', 'coordinates': [51.5019406, -0.1186677]}
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]for label, client in (("good", good), ("transposed", transposed)):
body = client.get("/features").json()
body["features"] = [f for f in body["features"]
if f["properties"]["country"] == "United Kingdom"]
print(label, "->", validate_coordinate_order(body, expected_region="gb") or "no problems")good -> no problems
! 7 features fall inside gb when swapped, against 0 as given — the pair is transposed
transposed -> ['7 features fall inside gb when swapped, against 0 as given — the pair is transposed']
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]wfs_style = Boundary(inbound=EPSG4326_AUTHORITY)
print(wfs_style.parse_bbox("50.8,-0.2,50.9,-0.1"), wfs_style.emit_point(-0.1276, 51.5072))(-0.2, 50.8, -0.1, 50.9) [-0.1276, 51.5072]
Example 3 — a regression test against real geometry
import geopandas as gpd
from shapely.geometry import Point, 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()
geoms = [shape(f["geometry"]).representative_point() for f in body["features"]]
def in_no_country(points):
gdf = gpd.GeoDataFrame(geometry=points, crs="EPSG:4326")
joined = gpd.sjoin(gdf, countries, how="left", predicate="within")
return int(joined["name"].isna().sum())
as_given = in_no_country(geoms)
swapped = in_no_country([Point(p.y, p.x) for p in geoms])
assert as_given < swapped, (
f"{as_given} of {len(geoms)} features are in no country as given, "
f"{swapped} when swapped — coordinates are probably transposed")for label, client in (("good", good), ("transposed", transposed)):
try:
test_features_land_in_the_right_country(client, "countries.gpkg")
print(f"{label:10} passed")
except AssertionError as exc:
print(f"{label:10} failed: {exc}")good passed
transposed failed: 144 of 200 features are in no country as given, 13 when swapped — 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.
Why the error is continents, not metres
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.
DuckDB transforms, Leaflet and bounding boxes
always_xy := trueon every DuckDBST_Transform. Its default follows the authority order.pyprojand GeoPandas default toalways_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.
Internal links
- My points plot in the ocean: latitude and longitude swapped — the same bug in a data pipeline
- CRS in DuckDB: why ST_Transform moves your data to the wrong place — the axis-order default
- OGC API Features explained — why CRS84 is named explicitly
- How to test a spatial API with pytest and httpx — the known-coordinate test
- Projected versus geographic CRS explained — the wider CRS picture
- EPSG codes explained — what the registry says
- How to clean CSV coordinates into a GeoDataFrame — the ingestion side
- Fixing a CRS mismatch in GeoPandas — the related failure
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.