A CRS is part of the data model, not display metadata. An incorrect or missing CRS can make numerically plausible results geographically wrong.
from pyproj import CRS
if gdf.crs is None:
raise ValueError("CRS is missing; recover it from authoritative source metadata")
crs = CRS.from_user_input(gdf.crs)
Do not infer EPSG:4326 because values resemble longitude/latitude. Coordinate ranges are not evidence of datum, axis interpretation, units, or epoch.
set_crs versus to_crs# Assign metadata only; coordinate numbers do not change.
gdf = gdf.set_crs("EPSG:4326")
# Transform coordinates; the source CRS must already be set.
projected = gdf.to_crs("EPSG:32631")
set_crs() only when the coordinate values are already expressed in that
CRS and metadata is absent or demonstrably wrong.allow_override=True; record why.to_crs() to transform the active geometry column. Transform each
additional geometry column explicitly.gdf.crs = ...; manual override is deprecated.EPSG definitions can be latitude/longitude while GIS coordinate arrays are normally x/y (longitude/latitude). Inspect the authority axes:
for axis in crs.axis_info:
print(axis.name, axis.abbrev, axis.direction, axis.unit_name)
For explicit pyproj transformations, request traditional GIS x/y order:
from pyproj import Transformer
transformer = Transformer.from_crs(
source_crs,
target_crs,
always_xy=True,
allow_ballpark=False,
only_best=True,
)
x_out, y_out = transformer.transform(x_in, y_in, errcheck=True)
Record always_xy=True; it changes the API coordinate order, not the CRS
definition. GeoParquet 1.1 explicitly stores WKB/native coordinates as x/y even
when the CRS authority uses another axis order.
GeoPandas and Shapely compute planar Cartesian geometry and ignore Z:
.area returns squared coordinate units;.length, .distance, .buffer, sjoin_nearest(max_distance=...),
sjoin(predicate="dwithin", distance=...), precision grids, simplify
tolerances, and gap widths use coordinate units.axis_units = [(axis.unit_name, axis.unit_conversion_factor) for axis in crs.axis_info]
if crs.is_geographic:
raise ValueError("Planar measurement on angular coordinates is not accepted")
Do not label an output metres merely because a CRS is projected. Convert units
using the CRS axis metadata and document the conversion. Web Mercator
(EPSG:3857) is for web display, not general area/distance analysis.
Choose based on the operation, study extent, datum, and required accuracy:
estimate_utm_crs() is a convenience based on dataset bounds, not a proof of
suitability. It can be poor for multi-zone, polar, antimeridian-crossing, or
very large datasets.
When projection distortion is unacceptable, use the ellipsoid associated with
the CRS through pyproj.Geod, not Shapely's planar distance:
geod = crs.get_geod()
azimuth_fwd, azimuth_back, metres = geod.inv(lon1, lat1, lon2, lat2)
area_m2, perimeter_m = geod.geometry_area_perimeter(polygon)
Ensure inputs are longitude/latitude on the intended geodetic datum. Geodesic area is signed according to ring orientation and has documented limitations for very large polygons; normalize orientation and test known controls.
The same source/target CRS pair can have several operations. Selection depends on area of interest, installed grids, authority, and accuracy:
from pyproj.aoi import AreaOfInterest
from pyproj.transformer import TransformerGroup
from pyproj import network
network.set_network_enabled(False)
group = TransformerGroup(
source_crs,
target_crs,
always_xy=True,
area_of_interest=AreaOfInterest(west, south, east, north),
allow_ballpark=False,
)
if not group.best_available:
raise RuntimeError("Best transformation unavailable; inspect missing grids")
candidate = group.transformers[0]
print(candidate.description, candidate.accuracy, candidate.area_of_use)
Privacy note: an exact area of interest can reveal a sensitive study location. Do not log it; retain only an approved coarse region or protected audit record.
Operational rules:
allow_ballpark=False for accuracy-sensitive work.only_best=True with Transformer.from_crs() when failure is preferable
to silently selecting a lower-quality operation.accuracy (-1 means unknown) and area_of_use.TransformerGroup.unavailable_operations for missing grids.The bundled scripts/crs_reprojection_plan.py performs this inspection without
transforming coordinates or enabling network access.
GeoDataFrame.to_crs() transforms every existing vertex. It does not interpret
a segment as a geodesic arc:
out = gdf.to_crs(target_crs)
If source linework is sparse, densify according to a documented geodesic or source-space tolerance before projection when shape fidelity matters. Validate the resulting topology and bounds.
to_crs() warns that objects crossing the dateline or another projection
boundary have undesirable behavior. A naive line from 179°E to 179°W can be
treated as spanning almost the whole map.
Safe workflow:
[-180, 180] or [0, 360)).For transformed bounds, use Transformer.transform_bounds(..., densify_pts=...).
When geographic output returns right < left, pyproj documents that the bounds
cross the antimeridian and should be represented as two polygons. Do not sort
the numbers and erase that meaning.
Compare semantic CRS objects, not raw WKT strings:
left_crs = CRS.from_user_input(left.crs)
right_crs = CRS.from_user_input(right.crs)
if not left_crs.equals(right_crs):
right = right.to_crs(left_crs)
Equivalent does not mean equally appropriate for the analysis. Before concat, join, overlay, or clip, require matching CRS and confirm both datasets use the same coordinate epoch/realization where relevant.
Record:
always_xy);