pyproj
Cython bindings around PROJ, the C library that nearly all open source geospatial software uses for coordinate reference systems and datum transformations. The pieces you touch are CRS for describing a reference system, Transformer for converting coordinates between two of them, and Geod for geodesic distance, azimuth and area on an ellipsoid. It operates on plain numbers and NumPy arrays rather than geometry objects, which is why GeoPandas, rasterio, cartopy and Shapely-based stacks all sit on top of it instead of reimplementing the maths.
The right and effectively only choice when Python code has to speak the same coordinate reference language as the rest of the geospatial world. Budget an afternoon for the two things that bite everyone: axis order, and the fact that accurate datum transformations need grid data the wheels do not ship.
Use it if
- You need to convert coordinates between reference systems and want the same answers the rest of the geospatial world gets, because it is the same PROJ underneath
- You are transforming large arrays, since Transformer.transform accepts NumPy arrays and does the work in one call rather than per point
- You need geodesic calculations on an ellipsoid: distances, forward and inverse azimuths, intermediate points, and polygon area and perimeter through Geod
- You want to inspect a CRS programmatically, including its axis order, datum, area of use, units and whether the EPSG code has been deprecated
- Your coordinates are all WGS84 and you only need approximate distances, where a haversine function in a few lines avoids a compiled dependency and a PROJ data directory entirely
- You want geometry operations. There is no clipping, buffering, intersection or reprojection of shapes here; that is shapely and geopandas territory, and pyproj is the layer they call
- You need accurate datum shifts out of the box: wheels since pyproj 3 ship no transformation grids at all, so a NAD27 or vertical datum conversion silently falls back to a lower accuracy path unless you enable PROJ network access or download grids yourself
- Your deployment cannot reach cdn.proj.org and cannot carry an extra 500 MB to 1 GB of grid data, because those are the two ways to get the accurate transformations
- You are stuck on an older Python or an older PROJ: 3.7.2 requires Python 3.11 or newer and PROJ 9.4 or newer, and building from source against a mismatched system PROJ is a common way to get an import that fails at runtime
- You expect axis order to match your intuition. EPSG:4326 is officially latitude then longitude, so Transformer.from_crs('EPSG:4326', 'EPSG:3857') without always_xy=True takes (lat, lon), and forgetting that produces coordinates in the wrong hemisphere rather than an error
Setup reality
pip install pyproj gives you binary wheels for the common platforms with PROJ built in, and its only Python dependency is certifi. That covers most people. What the wheels do not carry is transformation grids, and the documentation states this outright: for datum shifts that need a grid, PROJ either downloads it at runtime (set PROJ_NETWORK=ON or call pyproj.network.set_network_enabled(True)) or falls back to a less accurate ballpark transformation without telling you. The pyproj sync command line tool exists to pre-download grids into the user data directory for machines that cannot fetch at runtime, and the full set is 500 MB to 1 GB. Building from source is a separate project: you need a system PROJ that satisfies the version matrix (9.4 or newer for 3.7.x), and PROJ_DIR or PROJ_VERSION pointing at it, otherwise the import fails or the CRS database is missing. Mixing pip and conda in the same environment causes a distinct failure mode where the Python extension finds a PROJ data directory belonging to a different PROJ build. The axis order trap catches almost everyone once: authority-defined order for EPSG:4326 is latitude, longitude, so pass always_xy=True to Transformer.from_crs if your data is in longitude, latitude order like most GeoJSON and shapefiles. Also build the Transformer once and reuse it: from_crs does a database lookup and operation search each time, so calling it inside a loop is the usual reason a transformation feels slow.
Patterns
Convert between two reference systemstransform-coordinates
from pyproj import Transformer
transformer = Transformer.from_crs('EPSG:4326', 'EPSG:3857', always_xy=True)
x, y = transformer.transform(72.8777, 19.0760) # lon, lat
print(x, y)always_xy=True is what makes the call take longitude then latitude. Without it EPSG:4326 uses its authority order of latitude then longitude and you get a plausible but wrong answer.
Build the Transformer oncereuse-transformer
transformer = Transformer.from_crs(4326, 32643, always_xy=True)
def project_all(points):
return [transformer.transform(lon, lat) for lon, lat in points]from_crs searches the PROJ database for a transformation pipeline. Constructing it inside a loop repeats that search per point and dominates the runtime.
Transform whole arrays at oncetransform-numpy-arrays
import numpy as np
lons = np.array([72.87, 77.20, 88.36])
lats = np.array([19.07, 28.61, 22.57])
xs, ys = transformer.transform(lons, lats)Array input goes through one call into PROJ rather than one per point. Values that cannot be transformed come back as inf rather than raising, so check with np.isfinite.
Read what a CRS actually saysinspect-crs
from pyproj import CRS
crs = CRS.from_epsg(32643)
print(crs.name)
print(crs.axis_info)
print(crs.area_of_use)
print(crs.to_authority())
print(crs.is_projected, crs.is_geographic)axis_info is how you confirm the coordinate order before transforming. area_of_use tells you when a projected CRS is being applied outside the region it was defined for.
Accept WKT, PROJ strings or codesparse-any-crs-definition
crs = CRS.from_user_input('EPSG:4326')
crs = CRS.from_user_input('+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs')
crs = CRS.from_wkt(wkt_string)
print(crs.to_wkt(version='WKT2_2019'))from_user_input takes codes, WKT, PROJ strings and dicts. Prefer emitting WKT2 rather than a PROJ string when storing a CRS, since PROJ strings lose datum detail.
Distance and azimuth on an ellipsoidgeodesic-distance
from pyproj import Geod
geod = Geod(ellps='WGS84')
fwd_az, back_az, dist_m = geod.inv(72.8777, 19.0760, 77.2090, 28.6139)
print(round(dist_m / 1000, 1), 'km')inv returns metres on the ellipsoid, which differs from a spherical haversine result by a few tenths of a percent. Argument order is lon, lat for both points.
Area and perimeter of a polygongeodesic-area
lons = [72.80, 72.95, 72.95, 72.80]
lats = [18.90, 18.90, 19.10, 19.10]
area_m2, perimeter_m = geod.polygon_area_perimeter(lons, lats)
print(abs(area_m2))The sign of the area encodes winding direction, so take abs() unless you need orientation. Do not close the ring by repeating the first point.
Turn on PROJ network grid accessenable-grid-download
import pyproj
pyproj.network.set_network_enabled(True)
print(pyproj.network.is_network_enabled())
# or before the process starts
# export PROJ_NETWORK=ONWithout this, a transformation needing a grid quietly uses a lower accuracy fallback instead of failing. Enabling it means each transform may make an HTTPS call to cdn.proj.org.
Fetch grids ahead of time for offline usepredownload-grids
# list what would be downloaded for a bounding box
pyproj sync --bbox 68,6,98,38 --list-files
# download them into the user data directory
pyproj sync --bbox 68,6,98,38This is the option for machines with no outbound network at runtime. Downloading everything is 500 MB to 1 GB, so scope it with --bbox or --source-id.
See every candidate transformationcompare-transformation-options
from pyproj.transformer import TransformerGroup
group = TransformerGroup('EPSG:4326', 'EPSG:2154')
for t in group.transformers:
print(t.description, t.accuracy)
print(group.unavailable_operations)unavailable_operations lists the pipelines that need grids you do not have locally. This is the direct way to find out whether your accurate option is missing.
Reproject a bounding box safelytransform-bounds
left, bottom, right, top = transformer.transform_bounds(
68.0, 6.0, 98.0, 38.0, densify_pts=21
)Transforming only the four corners understates the extent when the projection curves edges. transform_bounds samples along the sides, and densify_pts controls how finely.
Confirm which PROJ you are actually usingcheck-versions
import pyproj
print(pyproj.__version__, pyproj.proj_version_str)
print(pyproj.datadir.get_data_dir())
pyproj.show_versions()Run this first when transformations behave oddly. A data directory pointing at a different PROJ install is the usual cause in mixed pip and conda environments.
Alternatives
| Package | Registry | Pick it when |
|---|---|---|
| geopandas | PyPI | You are working with tables of geometries and want to_crs on a DataFrame rather than raw coordinate arrays |
| geopy | PyPI | You mainly need geocoding and simple great-circle or geodesic distances, without a compiled PROJ dependency |
| cartopy | PyPI | The goal is drawing maps in matplotlib, where projection handling should be part of the plotting layer |