How to Measure Distance Accurately in Python: Geodesic vs Projected

On this page

Problem statement

Four ways to measure the flight from Leeds Bradford Airport to Paris Charles de Gaulle, four different answers (airport positions from OurAirports; the setup is under The data used on this page):

import geopandas as gpd
from pyproj import Geod

leeds, paris = AIRPORTS["EGNM"], AIRPORTS["LFPG"]

print(f"1. in degrees      {leeds.distance(paris):.3f}")

pts = gpd.GeoSeries([leeds, paris], crs=4326).to_crs(32630)
print(f"2. UTM zone 30N    {pts.iloc[0].distance(pts.iloc[1]) / 1000:,.1f} km")

pts = gpd.GeoSeries([leeds, paris], crs=4326).to_crs(3857)
print(f"3. Web Mercator    {pts.iloc[0].distance(pts.iloc[1]) / 1000:,.1f} km")

_, _, m = Geod(ellps="WGS84").inv(leeds.x, leeds.y, paris.x, paris.y)
print(f"4. geodesic        {m / 1000:,.1f} km")

1. in degrees      6.431
2. UTM zone 30N    614.7 km
3. Web Mercator    987.1 km
4. geodesic        614.5 km

The geodesic figure is the correct one. The first is not a distance at all. The UTM figure is close but a few hundred metres out, because Paris is outside zone 30. Web Mercator is wrong by well over half.

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

Scale and accuracy requirement decide the method. There is no single right 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	a few tenths of a percent (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.

What you get — from 1. Never measure in degrees

    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.30 km   1° lat = 111.32 km
   60°  1° lon =  55.66 km   1° lat = 111.32 km
   70°  1° lon =  38.07 km   1° lat = 111.32 km

No degrees, Geod for exact, projected for local

1. Never measure in degrees

print(leeds.distance(paris))               # a number of degrees — step 1 above

That is Euclidean distance over longitude and latitude treated as if they were the same unit. They are not: at Leeds' 53.8°N one degree of longitude is about 66 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.30 km   1° lat = 111.32 km
   60°  1° lon =  55.66 km   1° lat = 111.32 km
   70°  1° lon =  38.07 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(leeds.x, leeds.y, paris.x, paris.y)
print(f"{metres:,.1f} m   bearing {azimuth_fwd:.1f}°")

614,465.0 m   bearing 149.9°

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(leeds.x, leeds.y, 45.0, 100_000)
print(f"100 km northeast of the airport: {lon:.4f}, {lat:.4f}")

100 km northeast of the airport: -0.5692, 54.4963

3. 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

import time

t0 = time.perf_counter()
sites["dist_to_hq_m"] = geodesic_distance(sites, hq)
print(f"{len(sites):,} stops in {time.perf_counter() - t0:.3f} s; "
      f"median {sites['dist_to_hq_m'].median() / 1000:.1f} km from Leeds station")

16,605 stops in 0.010 s; median 16.3 km from Leeds station

Fast 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.036%   worst   0.037%   (18 m on the worst pair)
32630      median  0.031%   worst   0.035%   (14 m on the worst pair)
3857       median 68.915%   worst  69.708%   (34,297 m on the worst pair)

A national grid or UTM zone is accurate to a few hundredths of a percent over West Yorkshire — 18 m on the worst of these pairs, which run to tens of kilometres, and 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

Mercator's scale factor is 1/cos(latitude). At 60°N every distance is doubled.

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:

import numpy as np

for lat, place in [(0, "equator"), (51.5, "London"), (53.8, "Leeds"), (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.6%
Leeds      53.8°   distances overstated by   69.3%
Oslo         60°   distances overstated by  100.0%
Tromsø       71°   distances overstated by  207.2%

A 500 m buffer computed in EPSG:3857 at London's latitude is actually about 311 m on the ground, and at Leeds' about 295 m. 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, a few tenths of a percent off 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, (a, b) in {
    "Leeds–Paris":     ("EGNM", "LFPG"),
    "London–New York": ("EGLL", "KJFK"),
    "Sydney–Santiago": ("YSSY", "SCEL"),
}.items():
    lon1, lat1, lon2, lat2 = AIRPORTS[a].x, AIRPORTS[a].y, AIRPORTS[b].x, AIRPORTS[b].y
    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}%")

Leeds–Paris        haversine     613.8 km   exact     614.5 km    0.11%
London–New York    haversine   5,539.8 km   exact   5,554.7 km    0.27%
Sydney–Santiago    haversine  11,340.4 km   exact  11,362.6 km    0.20%

Haversine models the earth as a sphere, so it is off by up to a few tenths of a percent — the earth's flattening, which is about 0.34%. 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.

NaPTAN stops in West Yorkshire measured from Leeds station

Everything else in the function follows from this one threshold.

The data used on this page

sites.gpkg is every active bus stop, station entrance and ferry point in West Yorkshire from the Department for Transport's NaPTAN, and hq is Leeds station's main entrance from the same file. AIRPORTS holds six airports' positions from OurAirports' open data, and flight_routes.gpkg joins some of them with straight lines — plus Leeds to Bradford Interchange by rail station — to have lines of three very different lengths. Everything is saved in WGS 84, the way it arrives.

import io

import geopandas as gpd
import pandas as pd
import requests
from shapely.geometry import LineString, Point

# NaPTAN area 450 as saved on 9 October 2026 (the register changes weekly); live:
# naptan.api.dft.gov.uk/v1/access-nodes?dataFormat=csv&atcoAreaCodes=450
naptan = pd.read_csv(io.StringIO(requests.get("https://data.spatialworkflow.io/naptan-area-450-2026-10-09.csv", timeout=300).text),
                     low_memory=False)
naptan = naptan[naptan["Status"] == "active"]
sites = gpd.GeoDataFrame(naptan[["ATCOCode", "CommonName", "StopType"]].reset_index(drop=True), crs=4326,
                         geometry=gpd.points_from_xy(naptan["Longitude"], naptan["Latitude"]))
sites.to_file("sites.gpkg")
hq = sites.loc[sites["ATCOCode"] == "4500LEEDS0", "geometry"].iloc[0]          # Leeds station
bradford = sites.loc[sites["ATCOCode"] == "4500BRADIN0", "geometry"].iloc[0]   # Bradford Interchange

airports = pd.read_csv("https://davidmegginson.github.io/ourairports-data/airports.csv",
                       usecols=["ident", "name", "latitude_deg", "longitude_deg"])
airports = airports[airports["ident"].isin(["EGNM", "LFPG", "EGLL", "KJFK", "YSSY", "SCEL"])]
AIRPORTS = {r.ident: Point(r.longitude_deg, r.latitude_deg) for r in airports.itertuples()}

routes = gpd.GeoDataFrame({"name": ["Leeds–Bradford (rail)", "Leeds–Paris", "London–New York"]}, crs=4326,
                          geometry=[LineString([hq, bradford]),
                                    LineString([AIRPORTS["EGNM"], AIRPORTS["LFPG"]]),
                                    LineString([AIRPORTS["EGLL"], AIRPORTS["KJFK"]])])
routes.to_file("flight_routes.gpkg")
print(f"{len(sites):,} stops, {len(AIRPORTS)} airports, {len(routes)} routes")

16,605 stops, 6 airports, 3 routes

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.030% apart, worst 0.032%

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 within 0.04% 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")          # EPSG:4326
print(f"one UTM zone for the whole layer: {routes.estimate_utm_crs().name}\n")
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))

one UTM zone for the whole layer: WGS 84 / UTM zone 25N

                 name  geodesic_km  utm_km  utm_err_%  merc_err_%
Leeds–Bradford (rail)         13.4    14.1        5.1        68.9
          Leeds–Paris        614.5   654.5        6.5        60.6
      London–New York       5554.7  5758.1        3.7        50.3

The UTM column is the lesson, and it is not the one you might expect. Every route is wrong by several percent — even the 13 km line between two Leeds and Bradford stations — because estimate_utm_crs() picks one zone for the whole layer, from its centre, and the centre of a layer reaching from Leeds to New York is out over the Atlantic in zone 25. Leeds is 26 degrees of longitude from that zone's central meridian. Projected the right way, one zone per local feature, the short line would be as accurate as the stops were in step 4; the transatlantic one has no zone that fits. Web Mercator is wrong on all three, by between half and two-thirds.

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          nan m
mercator   effective radius        295.7 m
correct    effective radius        499.7 m

Three outcomes. Buffering by 500 in degrees produces a polygon so far outside valid coordinates that even its geodesic area comes back as NaN — obvious the moment it is plotted. Buffering by 500 in Web Mercator produces about 296 m at Leeds' latitude, which is not obvious at all: it looks like a 500 m buffer, plots correctly against a basemap, and is 41% 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.

Distance depends on the surface

The shortest path on an ellipsoid is not a straight line in any flat representation.

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 a few tenths of a percent — 0.11 to 0.27% on the three routes in step 6 — 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.

Example 2’s routes: even the 13 km line is 5% out in the one zone chosen for the whole layer.

Argument order, geometry_length and Mercator stretch

Geod.inv takes lon, lat, lon, lat — x before y, the opposite of Folium's convention.
Geod.inv is vectorised. Pass NumPy arrays; 184,000 distances take about 0.4 s.
Geod.geometry_length measures a LineString geodesically, segment by segment.
geod.fwd solves 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() not centroid for polygon distances — a centroid can fall outside.
Geod(ellps="GRS80") versus "WGS84" differ by under a millimetre; the choice rarely matters.
sjoin_nearest measures in the frame's CRS, so reproject before using max_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
How to compute distance to the coast for many points — the same question at scale

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 a few tenths of a percent — up to 0.27% on the routes measured here. 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.