Reprojecting CRS for Accurate Distance Calculations

TL;DR: Never compute distances, buffers, or nearest-neighbour joins on EPSG:4326 (longitude/latitude degrees). Detect the right metric CRS with utm = gdf.estimate_utm_crs(), reproject with gdf_m = gdf.to_crs(utm), then confirm gdf_m.crs.axis_info[0].unit_name == "metre". All subsequent .distance(), .buffer(), and sjoin_nearest(max_distance=...) calls now return and accept metres.

Why This Matters #

A degree of longitude is roughly 111 km at the equator and zero at the poles, so treating latitude/longitude as if they were planar coordinates makes every Euclidean distance wrong — and wrong by a factor that changes with latitude. Buffers become ellipses of unintended size, sjoin_nearest matches the wrong features, and any distance-based spatial weight matrix or variogram is built on nonsense lags. This page is a core step in the GeoPandas Data Preparation guide, and correct units underpin everything downstream in the Python Workflows for Spatial Modeling & Regression stack — kriging ranges, cross-validation buffers, and neighbourhood definitions all assume metric coordinates.

The fix is one reprojection at the start of the pipeline, done once, verified once. Getting it right removes a whole class of silent errors that otherwise surface only as inexplicably poor model performance much later.

The Problem, Formally #

Euclidean distance between two points assumes orthogonal axes in the same unit:

$$d = \sqrt{(\Delta x)^2 + (\Delta y)^2}.$$

In EPSG:4326, $\Delta x$ is a longitude difference in degrees and $\Delta y$ a latitude difference in degrees, but a degree of longitude spans

$$111{,}320 \cdot \cos(\phi)\ \text{metres}$$

at latitude $\phi$, while a degree of latitude is a near-constant $\approx 110{,}540$ m. The two axes are in incompatible physical units, so the formula above is invalid. A projected CRS maps the ellipsoid onto a plane where both axes are metres, restoring the identity.

Environment #

pip install \
  geopandas==0.14.4 \
  pyproj==3.6.1 \
  shapely==2.0.4 \
  pyogrio==0.7.2

import geopandas as gpd
from shapely.geometry import Point

gpd.options.io_engine = "pyogrio"

Step-by-Step Implementation #

Step 1 — Inspect the Current CRS and Its Units #

Before anything, find out whether the data are in degrees. crs.axis_info exposes the axis unit; crs.is_geographic is True for lat/lon systems.

gdf = gpd.read_file("stations.geojson")   # commonly EPSG:4326
print(gdf.crs)
print("geographic?", gdf.crs.is_geographic)
print("unit:", gdf.crs.axis_info[0].unit_name)
# EPSG:4326 -> geographic? True   unit: degree

If is_geographic is True, any distance you compute now is invalid until you reproject.

Step 2 — Estimate the Correct UTM Zone #

estimate_utm_crs reads the data extent and returns the UTM CRS covering it, with the correct hemisphere. This avoids manually looking up zone numbers.

utm_crs = gdf.estimate_utm_crs()
print(utm_crs)                    # e.g. EPSG:32633 (UTM zone 33N)
print("unit:", utm_crs.axis_info[0].unit_name)   # metre

UTM is the right default when your study area spans up to a few hundred kilometres and you care about local distance fidelity.

Step 3 — Reproject with to_crs #

to_crs transforms every geometry into the target CRS. It returns a new GeoDataFrame; the original is untouched.

gdf_m = gdf.to_crs(utm_crs)
print(gdf_m.crs.axis_info[0].unit_name)   # metre

For continental or global area calculations, target an equal-area CRS instead, for example an Albers or Lambert azimuthal equal-area definition appropriate to your region:

# Example: Europe Lambert azimuthal equal-area
gdf_ea = gdf.to_crs("EPSG:3035")

Step 4 — Verify Units and a Known Distance #

Trust nothing until a known distance checks out. Compute the separation between two points whose real-world distance you can sanity-check.

assert gdf_m.crs.axis_info[0].unit_name == "metre", "Not a metric CRS"

# Distance between the first two stations, now in metres
d = gdf_m.geometry.iloc[0].distance(gdf_m.geometry.iloc[1])
print(f"Separation: {d:,.1f} m")

# Compare against the geodesic (true ellipsoidal) distance as a cross-check
geod = gdf.crs.get_geod() if gdf.crs.get_geod() else None

A projected distance within a fraction of a percent of the geodesic distance confirms the projection is suitable for the area. Large discrepancies mean the extent is too big for a single UTM zone.

Step 5 — Run Distance-Dependent Operations Safely #

With metric coordinates, buffers and nearest joins now take and return metres.

# 500 m buffer around each station
buffers = gdf_m.copy()
buffers["geometry"] = gdf_m.geometry.buffer(500)

# Nearest-neighbour join with a real 1 km distance cap
from_gdf = gdf_m
matched = gpd.sjoin_nearest(from_gdf, gdf_m, max_distance=1000, how="left")

Interpreting the Output #

After Step 4, gdf_m.crs.axis_info[0].unit_name must read metre and a spot-checked separation must match reality. If a coastal station 30 km from its neighbour reports 0.27 in the projected CRS, you are still in degrees; if it reports 29,980 m, you are correct. The geodesic cross-check matters most near a UTM zone edge, where projected distance distortion grows: UTM error stays under about 0.1% within its 6-degree-wide zone but climbs beyond it, which is the practical signal that your area needs a custom projection or a wider equal-area CRS.

Once units are confirmed metric, every downstream metre-valued parameter — a kriging range, a buffered-CV exclusion radius, a DistanceBand threshold — is finally interpretable and portable.

Critical Best Practices #

Reproject Once, Early, and Assert #

Do the reprojection at ingestion and immediately assert not gdf.crs.is_geographic. Threading degree-CRS data deep into a pipeline before someone computes a distance is how silent errors survive to production. Systematic harmonisation across all inputs is the job of the parent GeoPandas Data Preparation workflow.

Align Every Layer to the Same CRS Before Joining #

sjoin, sjoin_nearest, and overlay silently mis-locate features if layers differ in CRS, and GeoPandas will warn but still run. Reproject all participating layers to the one chosen metric CRS first; never mix a degree layer with a metre layer in a spatial operation.

Choose the Projection by Task, Not Habit #

UTM preserves local distance and angle but distorts area far from its central meridian; equal-area projections preserve area but distort distance. Pick by what you will measure: distances and buffers favour UTM; density, area-normalised rates, and zonal statistics favour an equal-area CRS. Do not reuse a web-mapping EPSG:3857 (Web Mercator) for measurement — it distorts distance badly away from the equator.

Reproject Vectors, Not Rasters, When Sampling #

When aligning points to a raster, reproject the points to the raster CRS rather than resampling the raster to the points. Reprojecting a raster resamples every pixel and changes values; moving the vector layer is exact and cheap. This mirrors the guidance in Chunked Raster Processing with Dask-GeoPandas.

Watch Antimeridian and Polar Edge Cases #

estimate_utm_crs fails gracefully for most areas but a study region straddling the antimeridian or a pole needs a deliberate polar or custom CRS. Check the returned zone against your known geography rather than trusting it blindly for global datasets.

Troubleshooting #

Next Steps #

With coordinates in metres, distance-based joins scale better too; see Optimizing GeoPandas Spatial Joins for Large Datasets for index strategies, and return to the parent GeoPandas Data Preparation guide for the full cleaning and harmonisation sequence.

Related:

• Optimizing GeoPandas Spatial Joins for Large Datasets — distance-capped joins that require a metric CRS
• Spatial Weight Matrices — distance-band and KNN weights that assume metric coordinates
• Chunked Raster Processing with Dask-GeoPandas — why to reproject vectors, not rasters, when sampling

← Back to GeoPandas Data Preparation