Choosing the Right UTM Zone Automatically in Python

Metric analysis needs a projected CRS, and for most local-to-regional work that means the correct UTM zone — but hard-coding EPSG:32633 breaks the moment your data moves one zone east. This guide picks the right UTM zone automatically from the data's own extent. It is for anyone who buffers, measures, or computes areas and is tired of guessing EPSG codes. It sits under Coordinate Systems with PyProj in Mastering Core Geospatial Python Libraries.

Why This Approach / What Goes Wrong

UTM divides the world into 60 six-degree zones, each with its own EPSG code (326XX north, 327XX south). Pick the wrong zone and distances stretch the further your data sits from the zone's central meridian; pick a single zone for data spanning several and the distortion is worse still. GeoPandas' estimate_utm_crs() solves the common case by reading the data's centroid and returning the matching zone — no manual lookup. The trap is calling it on data with no CRS (it can't estimate without knowing the input) or on data spanning many zones, where no single UTM zone is appropriate and you should use an equal-area projection instead.

Prerequisites

conda install -c conda-forge "geopandas=0.14.*" "pyproj=3.4.*"

Step-by-Step Implementation

1. Load data and confirm it has a CRS — estimation needs to know the input.

import geopandas as gpd

survey_points = gpd.read_file("survey_points.gpkg")
if survey_points.crs is None:
    survey_points = survey_points.set_crs(epsg=4326)   # assert the known input CRS

2. Let GeoPandas estimate the UTM zone from the extent.

utm_crs = survey_points.estimate_utm_crs()
print(utm_crs.to_epsg(), "-", utm_crs.name)
# 32633 - WGS 84 / UTM zone 33N

3. Reproject and do metric work in the estimated zone.

survey_utm = survey_points.to_crs(utm_crs)
survey_utm["buffer_50m"] = survey_utm.geometry.buffer(50)   # 50 metres, correct

4. For a manual check or a pure-pyproj pipeline, derive the zone yourself from longitude.

from pyproj import CRS
from pyproj.aoi import AreaOfInterest
from pyproj.database import query_utm_crs_info

minx, miny, maxx, maxy = survey_points.to_crs(epsg=4326).total_bounds
utm_info = query_utm_crs_info(
    datum_name="WGS 84",
    area_of_interest=AreaOfInterest(minx, miny, maxx, maxy),
)
manual_crs = CRS.from_epsg(utm_info[0].code)
print(manual_crs.to_epsg())   # 32633

Verification

Confirm the estimate is a metric CRS and that an east–west span doesn't exceed one zone.

# The chosen CRS must use metres, not degrees
assert survey_utm.crs.axis_info[0].unit_name == "metre"
assert str(survey_utm.crs.to_epsg()).startswith(("326", "327"))

# Warn if the data spans more than ~6° of longitude (one UTM zone)
minx, _, maxx, _ = survey_points.to_crs(epsg=4326).total_bounds
span = maxx - minx
print(f"Longitude span: {span:.2f}°")          # Longitude span: 1.84°
assert span < 6, "Data spans multiple UTM zones — use an equal-area CRS instead"

Edge Cases & Debugging

Related