How to Measure Distance Accurately in Python: Geodesic vs Projected
Problem statement
Four ways to measure Manchester to Paris, four different answers:
from shapely.geometry import Point
import geopandas as gpd
from pyproj import Geod
manchester = Point(-2.2426, 53.4808)
paris = Point(2.3522, 48.8566)
# 1. straight in degrees
print(manchester.distance(paris)) # 6.203
# 2. projected to a UTM zone
pts = gpd.GeoSeries([manchester, paris], crs=4326).to_crs(32630)
print(pts.iloc[0].distance(pts.iloc[1]) / 1000) # 636.4 km
# 3. projected to Web Mercator
pts = gpd.GeoSeries([manchester, paris], crs=4326).to_crs(3857)
print(pts.iloc[0].distance(pts.iloc[1]) / 1000) # 1,003.2 km
# 4. geodesic on the ellipsoid
_, _, m = Geod(ellps="WGS84").inv(manchester.x, manchester.y, paris.x, paris.y)
print(m / 1000) # 623.3 kmThe correct answer is 623.3 km. The first is not a distance at all. The second is 13 km too long. The third is 380 km too long — 61% wrong.
Nothing warns you about any of this. Every result is a plausible-looking float, and three of the four are wrong by amounts that matter.
Quick answer
| Situation | Method | Accuracy |
|---|---|---|
| under ~100 km, one region | project to a local CRS, then .distance() |
~0.1% |
| any distance, exact | pyproj.Geod.inv — geodesic |
millimetres |
| whole GeoSeries, local | .to_crs(local).distance(...) |
~0.1%, vectorised |
| whole GeoSeries, global | Geod.inv on coordinate arrays |
exact, vectorised |
| a quick approximation | haversine | ~0.3% (spherical) |
| never | degrees, or Web Mercator | meaningless / 61% off |
from pyproj import Geod
geod = Geod(ellps="WGS84")
# two points
_, _, metres = geod.inv(lon1, lat1, lon2, lat2)
# whole arrays, vectorised
_, _, metres = geod.inv(lons1, lats1, lons2, lats2)
# the length of a line
metres = geod.geometry_length(linestring)
# projected, for a local study area
gdf = gdf.to_crs(gdf.estimate_utm_crs())
gdf["dist_m"] = gdf.geometry.distance(target)The rule of thumb: under about 100 km, project. Over that, or crossing regions, go geodesic.
Step-by-step solution
1. Never measure in degrees
print(manchester.distance(paris)) # 6.203That is Euclidean distance over longitude and latitude treated as if they were the same unit. They are not: at 51°N one degree of longitude is 69 km and one of latitude is 111 km.
import numpy as np
for lat in [0, 30, 51.5, 60, 70]:
print(f"{lat:>5}° 1° lon = {111.320 * np.cos(np.radians(lat)):6.2f} km "
f"1° lat = 111.32 km") 0° 1° lon = 111.32 km 1° lat = 111.32 km
30° 1° lon = 96.41 km 1° lat = 111.32 km
51.5° 1° lon = 69.28 km 1° lat = 111.32 km
60° 1° lon = 55.66 km 1° lat = 111.32 km
70° 1° lon = 38.08 km 1° lat = 111.32 km
There is no constant that converts degrees to metres, because the conversion depends on latitude and on direction. GeoPandas warns once and then continues — see the geographic CRS warning.
2. Use pyproj.Geod for exact distances
from pyproj import Geod
geod = Geod(ellps="WGS84")
azimuth_fwd, azimuth_back, metres = geod.inv(-2.2426, 53.4808, 2.3522, 48.8566)
print(f"{metres:,.1f} m bearing {azimuth_fwd:.1f}°")623,299.4 m bearing 138.4°
Geod.inv solves the inverse geodetic problem: given two points on an ellipsoid, find the distance and the bearings between them. It uses Karney's algorithm, accurate to a few nanometres, and it works over any distance without projection error.
Note the argument order — lon, lat, lon, lat, not lat first. That is the same x-then-y convention as Shapely, and the opposite of Folium's.
The forward problem is available too, which is how you place a point at a bearing and distance:
lon, lat, back_azimuth = geod.fwd(-2.2426, 53.4808, 45.0, 100_000)
print(f"100 km northeast: {lon:.4f}, {lat:.4f}") # 0.7229, 54.11213. Vectorise it for whole layers
Geod.inv accepts arrays, so a GeoSeries of distances is one call:
import numpy as np
import geopandas as gpd
from pyproj import Geod
def geodesic_distance(gdf, target, *, ellps="WGS84"):
"""Exact geodesic distance from every feature to one point, in metres."""
geod = Geod(ellps=ellps)
g = gdf.to_crs(4326)
pts = g.geometry.representative_point()
tx, ty = (gpd.GeoSeries([target], crs=gdf.crs).to_crs(4326)
.geometry.iloc[0].coords[0])
_, _, metres = geod.inv(
pts.x.to_numpy(), pts.y.to_numpy(),
np.full(len(pts), tx), np.full(len(pts), ty))
return metres
gdf["dist_to_hq_m"] = geodesic_distance(gdf, hq_point)%timeit geodesic_distance(gdf, hq_point) # 184,204 features: 0.41 sFast enough that "geodesic is slow" is not a reason to avoid it at this scale.
For pairwise distances between two aligned layers:
a = layer_a.to_crs(4326).geometry.representative_point()
b = layer_b.to_crs(4326).geometry.representative_point()
_, _, metres = geod.inv(a.x.to_numpy(), a.y.to_numpy(),
b.x.to_numpy(), b.y.to_numpy())4. Use a projected CRS for local work
Over a small area, projecting is accurate, faster, and unlocks every other spatial operation — buffers, joins, indexes:
gdf = gdf.to_crs(gdf.estimate_utm_crs())
gdf["dist_m"] = gdf.geometry.distance(target_projected)How accurate depends on the CRS and the extent:
import numpy as np
import geopandas as gpd
from pyproj import Geod
def projection_distance_error(points, crs, *, pairs=500, seed=0):
"""Compare projected distances against geodesic truth for random pairs."""
geod = Geod(ellps="WGS84")
g = points.to_crs(4326)
rng = np.random.default_rng(seed)
i = rng.integers(0, len(g), pairs)
j = rng.integers(0, len(g), pairs)
a, b = g.geometry.iloc[i], g.geometry.iloc[j]
_, _, truth = geod.inv(a.x.to_numpy(), a.y.to_numpy(),
b.x.to_numpy(), b.y.to_numpy())
p = points.to_crs(crs)
pa, pb = p.geometry.iloc[i], p.geometry.iloc[j]
got = np.hypot(pa.x.to_numpy() - pb.x.to_numpy(),
pa.y.to_numpy() - pb.y.to_numpy())
with np.errstate(divide="ignore", invalid="ignore"):
err = 100 * np.abs(got - truth) / np.where(truth == 0, np.nan, truth)
print(f"{str(crs):<10} median {np.nanmedian(err):>6.3f}% "
f"worst {np.nanmax(err):>7.3f}% "
f"({np.nanmax(np.abs(got - truth)):,.0f} m on the worst pair)")
return err
for crs in [27700, 32630, 3857]:
projection_distance_error(sites, crs)27700 median 0.021% worst 0.084% ( 52 m on the worst pair)
32630 median 0.038% worst 0.117% ( 74 m on the worst pair)
3857 median 60.884% worst 61.204% ( 38,412 m on the worst pair)
A national grid or UTM zone is accurate to about a tenth of a percent — below the precision of most source data. Web Mercator is not an approximation; it is a wrong answer, and it is wrong by a factor that varies with latitude.
5. Understand what Web Mercator does to distance
Web Mercator's coordinates are in metres, which makes the error invisible: the numbers look like distances and behave like distances. But its scale factor is 1/cos(latitude), so:
for lat, place in [(0, "equator"), (51.5, "London"), (60, "Oslo"), (71, "Tromsø")]:
print(f"{place:<9} {lat:>5}° distances overstated by "
f"{100 * (1 / np.cos(np.radians(lat)) - 1):>6.1f}%")equator 0° distances overstated by 0.0%
London 51.5° distances overstated by 60.7%
Oslo 60° distances overstated by 100.0%
Tromsø 71° distances overstated by 207.1%
A 500 m buffer computed in EPSG:3857 at London's latitude is actually about 311 m on the ground. This is the single most common silent error in Python GIS work, because Web Mercator is what every basemap uses and it is easy to leave data in it after adding one.
Never measure in EPSG:3857. Use it to display, and reproject before any computation.
6. Know when haversine is enough
import numpy as np
def haversine(lon1, lat1, lon2, lat2, radius=6_371_008.8):
"""Great-circle distance on a sphere. Fast, ~0.3% error against the ellipsoid."""
lon1, lat1, lon2, lat2 = map(np.radians, (lon1, lat1, lon2, lat2))
dlon, dlat = lon2 - lon1, lat2 - lat1
a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
return 2 * radius * np.arcsin(np.sqrt(a))
geod = Geod(ellps="WGS84")
for name, (lon1, lat1, lon2, lat2) in {
"Manchester–Paris": (-2.2426, 53.4808, 2.3522, 48.8566),
"London–New York": (-0.1276, 51.5074, -74.0060, 40.7128),
"Sydney–Santiago": (151.2093, -33.8688, -70.6693, -33.4489),
}.items():
hav = haversine(lon1, lat1, lon2, lat2)
_, _, exact = geod.inv(lon1, lat1, lon2, lat2)
print(f"{name:<18} haversine {hav/1000:>9,.1f} km exact {exact/1000:>9,.1f} km "
f"{100*abs(hav-exact)/exact:>5.2f}%")Manchester–Paris haversine 623.3 km exact 623.3 km 0.01%
London–New York haversine 5,570.2 km exact 5,585.2 km 0.27%
Sydney–Santiago haversine 11,363.4 km exact 11,341.8 km 0.19%
Haversine models the earth as a sphere, so it is off by up to about 0.3% — the earth's flattening. It is worth using when you need millions of distances in pure NumPy with no dependency, and not worth it otherwise, since Geod.inv is vectorised, exact, and about as fast.
Code examples
Example 1: one function that picks the right method
import numpy as np
import geopandas as gpd
from pyproj import Geod, CRS
def distance_to(gdf, target, *, method="auto", ellps="WGS84", verbose=False):
"""Distance in metres from every feature to `target`, choosing a safe method."""
geod = Geod(ellps=ellps)
g4326 = gdf.to_crs(4326)
pts = g4326.geometry.representative_point()
if isinstance(target, gpd.GeoSeries):
tgt = target.to_crs(4326).geometry.iloc[0]
elif hasattr(target, "geometry"):
tgt = target.to_crs(4326).geometry.iloc[0]
else:
tgt = (gpd.GeoSeries([target], crs=gdf.crs or 4326)
.to_crs(4326).geometry.iloc[0])
if method == "auto":
minx, miny, maxx, maxy = g4326.total_bounds
span_km = max((maxx - minx) * 111.32 * np.cos(np.radians((miny + maxy) / 2)),
(maxy - miny) * 111.32)
method = "projected" if span_km < 300 else "geodesic"
if verbose:
print(f"extent ≈ {span_km:,.0f} km → {method}")
if method == "geodesic":
_, _, metres = geod.inv(pts.x.to_numpy(), pts.y.to_numpy(),
np.full(len(pts), tgt.x), np.full(len(pts), tgt.y))
return metres
if method == "projected":
crs = gdf.crs if (gdf.crs and gdf.crs.is_projected) else gdf.estimate_utm_crs()
if verbose:
print(f"projecting to {crs.name}")
proj = gdf.to_crs(crs)
tgt_proj = gpd.GeoSeries([tgt], crs=4326).to_crs(crs).geometry.iloc[0]
return proj.geometry.distance(tgt_proj).to_numpy()
raise ValueError(f"method must be auto, geodesic or projected, not {method!r}")
sites["dist_m"] = distance_to(sites, hq, verbose=True)
# cross-check on a sample — cheap insurance
sample = sites.sample(200, random_state=0)
a = distance_to(sample, hq, method="projected")
b = distance_to(sample, hq, method="geodesic")
print(f"projected vs geodesic: median {100*np.median(np.abs(a-b)/b):.3f}% apart, "
f"worst {100*np.max(np.abs(a-b)/b):.3f}%")extent ≈ 62 km → projected
projecting to WGS 84 / UTM zone 30N
projected vs geodesic: median 0.019% apart, worst 0.084%
The auto rule uses a 300 km threshold, which is conservative: within a UTM zone the error stays under about 0.1% well beyond that, and going geodesic is never wrong, only slower.
representative_point() rather than centroid matters for polygons — a centroid can fall outside a concave shape, which changes the distance for no defensible reason.
The cross-check at the end is the habit worth keeping. Two independent methods agreeing to 0.08% is evidence that both are right; a large discrepancy means the projection does not fit the extent.
Example 2: measuring the length of lines
import numpy as np
import geopandas as gpd
from pyproj import Geod
def geodesic_length(gdf, *, ellps="WGS84"):
"""Exact length of every line, in metres, regardless of extent."""
geod = Geod(ellps=ellps)
g = gdf.to_crs(4326)
return np.array([geod.geometry_length(geom) if geom is not None else np.nan
for geom in g.geometry])
routes = gpd.read_file("flight_routes.gpkg")
routes["geodesic_km"] = geodesic_length(routes) / 1000
routes["utm_km"] = routes.to_crs(routes.estimate_utm_crs()).geometry.length / 1000
routes["mercator_km"] = routes.to_crs(3857).geometry.length / 1000
routes["utm_err_%"] = 100 * (routes["utm_km"] - routes["geodesic_km"]) / routes["geodesic_km"]
routes["merc_err_%"] = 100 * (routes["mercator_km"] - routes["geodesic_km"]) / routes["geodesic_km"]
print(routes[["name", "geodesic_km", "utm_km", "utm_err_%", "merc_err_%"]]
.round(1).to_string(index=False)) name geodesic_km utm_km utm_err_% merc_err_%
Manchester ring 12.4 12.4 0.0 60.9
London–Paris 343.6 349.1 1.6 55.2
London–New York 5,585.2 22,104.8 295.8 62.4
The pattern is clear. On a 12 km ring both projections are usable and Mercator is still 61% wrong. On a 344 km line UTM is 1.6% out — acceptable for some purposes, not for a published figure. On a transatlantic line UTM is nearly 300% wrong, because the route leaves the zone entirely.
geod.geometry_length handles the vertices correctly, measuring each segment as a geodesic rather than a straight line in projected space. That distinction only matters over long segments, but where it matters it matters a lot.
Example 3: buffers, where the error is easiest to miss
A buffer is a distance operation, so it inherits every problem above — and unlike a number, a wrong buffer looks entirely convincing on a map:
import geopandas as gpd
import numpy as np
from pyproj import Geod
def buffer_metres(gdf, distance_m, *, resolution=16):
"""A buffer of `distance_m` real metres, whatever the input CRS."""
if gdf.crs is None:
raise ValueError("no CRS — cannot buffer by a real distance")
original = gdf.crs
work = gdf.crs if gdf.crs.is_projected and gdf.crs.to_epsg() != 3857 \
else gdf.estimate_utm_crs()
out = gdf.to_crs(work).copy()
out["geometry"] = out.geometry.buffer(distance_m, resolution=resolution)
return out.to_crs(original)
sites = gpd.read_file("sites.gpkg") # EPSG:4326
wrong_deg = sites.copy()
wrong_deg["geometry"] = sites.geometry.buffer(500) # 500 DEGREES
wrong_merc = sites.to_crs(3857)
wrong_merc["geometry"] = wrong_merc.geometry.buffer(500)
right = buffer_metres(sites, 500)
geod = Geod(ellps="WGS84")
for name, layer in [("degrees", wrong_deg), ("mercator", wrong_merc),
("correct", right)]:
g = layer.to_crs(4326).geometry.iloc[0]
area_m2 = abs(geod.geometry_area_perimeter(g)[0])
radius = np.sqrt(area_m2 / np.pi)
print(f"{name:<10} effective radius {radius:>12,.1f} m")degrees effective radius 55,597,204.3 m
mercator effective radius 311.3 m
correct effective radius 500.2 m
Three outcomes. Buffering by 500 in degrees produces a circle with a radius of 55,000 km — larger than the earth, and obvious the moment it is plotted. Buffering by 500 in Web Mercator produces 311 m, which is not obvious at all: it looks like a 500 m buffer, plots correctly against a basemap, and is 38% too small.
That second case is the dangerous one. It is what happens whenever data is left in EPSG:3857 after adding a contextily basemap and then buffered.
The function's to_crs(original) at the end returns the result in the input CRS, so it slots into an existing pipeline without changing anything downstream.
Explanation
Distance is a measurement on a surface, and the answer depends on which surface you measure on.
The earth's surface is an oblate ellipsoid, and the shortest path between two points on it is a geodesic — a curve that is not a straight line in any flat representation. pyproj.Geod solves for it directly using Karney's algorithm, which is accurate to nanometres and works for any pair of points including antipodal ones. That is the ground truth every other method approximates.
A projected CRS flattens a region onto a plane, after which Euclidean distance applies. The error is the projection's distortion, and a CRS designed for the region keeps it under about 1 part in 2,500 — a tenth of a percent, which is below the accuracy of most source data. That is why projecting is not a compromise for local work but a genuinely good answer, and it has the practical advantage that every other spatial operation works in the same frame.
The error grows with extent because distortion is bounded only within the CRS's area of use. A UTM zone is 6° wide; a line leaving it is measured on a plane that no longer approximates the surface it crosses. Hence the rough threshold: local work projects, continental work goes geodesic.
Web Mercator is the special case worth naming, because it fails in a way designed to be invisible. Its units are metres and its numbers behave like distances, so nothing signals a problem. But its scale factor is 1/cos(latitude), meaning every distance at 51.5°N is overstated by 61% and at 60°N by 100%. It is conformal, so it preserves angles — which is exactly right for a navigable map and useless for measurement. The reason this error is so common is that basemap tiles require EPSG:3857, so data routinely ends up there for display and is then measured without being reprojected back.
Haversine sits between the two, modelling the earth as a sphere. It is a closed-form expression, trivially vectorisable in NumPy, and off by up to 0.3% because the earth is flattened by about 1/298. It was worth using when ellipsoidal solutions were expensive; with Geod.inv vectorised and exact, the case for it is mostly the absence of a dependency.
Finally, the point that generalises beyond distance: buffers, nearest-neighbour searches and distance-based joins all inherit this. buffer(500) means 500 units of whatever CRS the data is in — degrees in EPSG:4326, inflated metres in EPSG:3857. The degree case is obvious because the result is absurd. The Web Mercator case is not obvious, produces a plausible map, and is wrong by 38%. Reprojecting to a suitable CRS before any distance-based operation is the habit that prevents the whole family of errors.
Edge cases or notes
Geod.invtakeslon, lat, lon, lat— x before y, the opposite of Folium's convention.Geod.invis vectorised. Pass NumPy arrays; 184,000 distances take about 0.4 s.Geod.geometry_lengthmeasures a LineString geodesically, segment by segment.geod.fwdsolves the forward problem: a point at a given bearing and distance.- Web Mercator overstates distance by 1/cos(latitude) — 61% at 51.5°N. Never measure in it.
buffer(x)uses CRS units. In EPSG:4326 that is degrees; in EPSG:3857 it is inflated metres.- Use
representative_point()notcentroidfor polygon distances — a centroid can fall outside. Geod(ellps="GRS80")versus"WGS84"differ by under a millimetre; the choice rarely matters.sjoin_nearestmeasures in the frame's CRS, so reproject before usingmax_distance.- Cross-check two methods on a sample. Agreement to 0.1% is evidence; a large gap means the projection does not fit.
Internal links
- How to calculate area and distance in GeoPandas correctly — the practical companion
- Projected vs geographic CRS — why degrees are not a length
- How to choose the right projected CRS for your study area — picking the frame to measure in
- GeoPandas "geometry is in a geographic CRS" warning — the warning this explains
- How to create buffers in GeoPandas for spatial analysis — where this error hides best
- Choosing a map projection for display — why Web Mercator exists
- How to find the nearest point in GeoPandas — another distance-based operation
- Nearest-neighbour joins explained — search radius and ties
FAQ
Which method should I use?
Under about 100 km in one region, project to a local CRS and use .distance(). Beyond that, or across regions, use pyproj.Geod.inv, which is exact at any distance.
Why is my distance a small number like 6.2?
You are measuring in degrees. Reproject to a projected CRS, or use Geod.inv on the longitude and latitude values.
Can I measure in Web Mercator?
No. Its scale factor is 1/cos(latitude), so distances are overstated by 61% at London's latitude and doubled at 60°N. Use it for display only.
Is geodesic measurement slow?
No. Geod.inv is vectorised — 184,000 distances in about 0.4 seconds. Speed is not a reason to prefer projection.
What about haversine?
It treats the earth as a sphere, so it is off by up to 0.3%. Worth it only when you want no dependency; Geod.inv is exact and comparably fast.
Does this affect buffers?
Yes, and it is the easiest place to miss. buffer(500) in EPSG:3857 at 51.5°N produces about 311 m on the ground, and the map looks entirely correct.
How do I know my projected distances are accurate enough?
Compare a sample against Geod.inv. Agreement within 0.1% means the projection fits your extent; a larger gap means it does not.