Getting Started

This section demonstrates recommended usage of the three main classes in pyproj:

New to coordinate reference systems, EPSG codes, or PROJ? Read the Concepts page first. Unfamiliar terms are defined in the Glossary. There are also more examples of usage within the API documentation and tests.

Also see: Gotchas/FAQ

Using the CRS class

A pyproj.crs.CRS object holds the definition of a Coordinate Reference System (CRS). For more usage examples and documentation see pyproj.crs.CRS.

Creating a CRS

The pyproj.crs.CRS class can be initialized from any of the common ways of describing a CRS. Here are some examples; all four produce the same CRS, WGS 84 latitude/longitude.

>>> from pyproj import CRS
>>> crs = CRS.from_epsg(4326)  # EPSG code as an integer
>>> crs = CRS.from_string("EPSG:4326")  # authority string "AUTHORITY:CODE"
>>> crs = CRS.from_proj4("+proj=latlon")  # PROJ string
>>> crs = CRS.from_user_input(4326)  # any of the above (and more)

4326 is an authority code in the EPSG registry and +proj=latlon is a PROJ string; see Concepts for what these mean. pyproj.crs.CRS.from_user_input() lists every kind of input that is accepted, including Well-Known Text (WKT) and PROJ JSON, and is what the CRS(...) constructor and most other pyproj functions use to interpret CRS arguments.

Converting CRS to a different format

A CRS can be exported as an authority code, Well-Known Text (WKT) in several versions, a PROJ string, or a CF grid mapping dictionary.

Warning

You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems

>>> from pyproj import CRS
>>> crs = CRS.from_epsg(4326)
>>> crs.to_epsg()
4326
>>> crs.to_authority()
('EPSG', '4326')
>>> crs = CRS.from_proj4("+proj=omerc +lat_0=-36 +lonc=147 +alpha=-54 +k=1 +x_0=0 +y_0=0 +gamma=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0")
>>> crs
<Bound CRS: +proj=omerc +lat_0=-36 +lonc=147 +alpha=-54 +k=1 + ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: Transformation from unknown to WGS84
- method: Position Vector transformation (geog2D domain)
Datum: Unknown based on WGS84 ellipsoid
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
Source CRS: unknown

>>> print(crs.to_wkt(pretty=True))
BOUNDCRS[
    SOURCECRS[
        PROJCRS["unknown",
            BASEGEOGCRS["unknown",
                DATUM["Unknown based on WGS84 ellipsoid",
                    ELLIPSOID["WGS 84",6378137,298.257223563,
                        LENGTHUNIT["metre",1],
                        ID["EPSG",7030]]],
...
        PARAMETER["Z-axis rotation",0,
            ID["EPSG",8610]],
        PARAMETER["Scale difference",1,
            ID["EPSG",8611]]]]

>>> from pyproj.enums import WktVersion
>>> print(crs.to_wkt(WktVersion.WKT1_GDAL, pretty=True))
PROJCS["unknown",
    GEOGCS["unknown",
        DATUM["Unknown_based_on_WGS84_ellipsoid",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]]],
    PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
    PARAMETER["latitude_of_center",-36],
    PARAMETER["longitude_of_center",147],
    PARAMETER["azimuth",-54],
    PARAMETER["rectified_grid_angle",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]

>>> from pprint import pprint
>>> pprint(crs.to_cf())
{'azimuth_of_central_line': -54,
'crs_wkt': 'BOUNDCRS[SOURCECRS[PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unknown '
...
            'difference",1,ID["EPSG",8611]]]]',
'false_easting': 0.0,
'false_northing': 0.0,
'grid_mapping_name': 'oblique_mercator',
'horizontal_datum_name': 'Unknown based on WGS84 ellipsoid',
'inverse_flattening': 298.257223563,
'latitude_of_projection_origin': -36.0,
'longitude_of_prime_meridian': 0.0,
'longitude_of_projection_origin': 147.0,
'prime_meridian_name': 'Greenwich',
'reference_ellipsoid_name': 'WGS 84',
'scale_factor_at_projection_origin': 1.0,
'semi_major_axis': 6378137.0,
'semi_minor_axis': 6356752.314245179,
'towgs84': [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]}

Extracting attributes from CRS

There are many attributes you can pull from the pyproj.crs.CRS. This is just a small subset of what is available. The example below uses an OGC URN that combines two EPSG codes (a projected CRS and a vertical CRS) into a single compound CRS.

>>> crs = CRS("urn:ogc:def:crs,crs:EPSG::2393,crs:EPSG::5717")
>>> crs
<Compound CRS: EPSG:3901>
Name: KKJ / Finland Uniform Coordinate System + N60 height
Axis Info [cartesian|vertical]:
- X[north]: Northing (metre)
- Y[east]: Easting (metre)
- H[up]: Gravity-related height (metre)
Area of Use:
- undefined
Datum: Kartastokoordinaattijarjestelma (1966)
- Ellipsoid: International 1924
- Prime Meridian: Greenwich
Sub CRS:
- KKJ / Finland Uniform Coordinate System
- N60 height
>>> crs.sub_crs_list
[<Projected CRS: EPSG:2393>
Name: KKJ / Finland Uniform Coordinate System
Axis Info [cartesian]:
- X[north]: Northing (metre)
- Y[east]: Easting (metre)
Area of Use:
- name: Finland - 25.5°E to 28.5°E onshore. Also all country.
- bounds: (19.24, 59.75, 31.59, 70.09)
Coordinate Operation:
- name: Finland Uniform Coordinate System
- method: Transverse Mercator
Datum: Kartastokoordinaattijarjestelma (1966)
- Ellipsoid: International 1924
- Prime Meridian: Greenwich
, <Vertical CRS: EPSG:5717>
Name: N60 height
Axis Info [vertical]:
- H[up]: Gravity-related height (metre)
Area of Use:
- name: Finland - onshore.
- bounds: (19.24, 59.75, 31.59, 70.09)
Datum: Helsinki 1960
- Ellipsoid: undefined
- Prime Meridian: undefined
]
>>> cop = crs.sub_crs_list[0].coordinate_operation
>>> print(cop.to_wkt(pretty=True))
CONVERSION["Finland Uniform Coordinate System",
    METHOD["Transverse Mercator",
        ID["EPSG",9807]],
    PARAMETER["Latitude of natural origin",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8801]],
    PARAMETER["Longitude of natural origin",27,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8802]],
    PARAMETER["Scale factor at natural origin",1,
        SCALEUNIT["unity",1],
        ID["EPSG",8805]],
    PARAMETER["False easting",3500000,
        LENGTHUNIT["metre",1],
        ID["EPSG",8806]],
    PARAMETER["False northing",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8807]]]
>>> cop.method_code
'9807'
>>> cop.method_name
'Transverse Mercator'
>>> cop.params
[Param(name=Latitude of natural origin, auth_name=EPSG, code=8801, value=0.0, unit_name=degree, unit_auth_name=, unit_code=, unit_category=angular),
 ...
 Param(name=False northing, auth_name=EPSG, code=8807, value=0.0, unit_name=metre, unit_auth_name=, unit_code=, unit_category=linear)]

Find UTM CRS by Latitude and Longitude

Note

For more database methods see: Database.

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

utm_crs_list = query_utm_crs_info(
    datum_name="WGS 84",
    area_of_interest=AreaOfInterest(
        west_lon_degree=-93.581543,
        south_lat_degree=42.032974,
        east_lon_degree=-93.581543,
        north_lat_degree=42.032974,
    ),
)
utm_crs = CRS.from_epsg(utm_crs_list[0].code)

Transformations from CRS to CRS

A pyproj.transformer.Transformer converts coordinates from one CRS to another, including any change of datum between them. See Concepts for background on what happens during a transformation.

Step 1: Inspect CRS definition to ensure proper area of use and axis order

For more options available for inspection, usage examples, and documentation see pyproj.crs.CRS.

>>> from pyproj import CRS
>>> crs_4326 = CRS.from_epsg(4326)
>>> crs_4326
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

>>> crs_26917 = CRS.from_epsg(26917)
>>> crs_26917
<Projected CRS: EPSG:26917>
Name: NAD83 / UTM zone 17N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: North America - 84°W to 78°W and NAD83 by country
- bounds: (-84.0, 23.81, -78.0, 84.0)
Coordinate Operation:
- name: UTM zone 17N
- method: Transverse Mercator
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

Note that crs_4326 has the latitude (north) axis first and the crs_26917 has the easting axis first. This means that in the transformation, we will need to input the data with latitude first and longitude second. Also, note that the second projection is a UTM projection with bounds (-84.0, 23.81, -78.0, 84.0) which are in the form (min_x, min_y, max_x, max_y), so the transformation input/output should be within those bounds for best results.

Step 2: Create Transformer to convert from CRS to CRS

The pyproj.transformer.Transformer can be initialized with anything supported by pyproj.crs.CRS.from_user_input(). There are a couple of examples added here for demonstration. For more usage examples and documentation, see pyproj.transformer.Transformer.

>>> from pyproj import Transformer
>>> transformer = Transformer.from_crs(crs_4326, crs_26917)
>>> transformer = Transformer.from_crs(4326, 26917)
>>> transformer = Transformer.from_crs("EPSG:4326", "EPSG:26917")
>>> transformer
<Unknown Transformer: unknown>
Inverse of NAD83 to WGS 84 (1) + UTM zone 17N
>>> transformer.transform(50, -80)
(571666.4475041276, 5539109.815175673)

If you prefer to always have the axis order in the x,y or lon,lat order, you can use the always_xy option when creating the transformer.

>>> from pyproj import Transformer
>>> transformer = Transformer.from_crs("EPSG:4326", "EPSG:26917", always_xy=True)
>>> transformer.transform(-80, 50)
(571666.4475041276, 5539109.815175673)

Converting between geographic and projection coordinates within one datum

Step 1: Retrieve the geodetic CRS based on original CRS

>>> from pyproj import CRS
>>> crs = CRS.from_epsg(3857)
>>> crs
<Projected CRS: EPSG:3857>
Name: WGS 84 / Pseudo-Mercator
Axis Info [cartesian]:
- X[east]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: World - 85°S to 85°N
- bounds: (-180.0, -85.06, 180.0, 85.06)
Coordinate Operation:
- name: Popular Visualisation Pseudo-Mercator
- method: Popular Visualisation Pseudo Mercator
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

>>> crs.geodetic_crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Step 2: Create Transformer to convert from geodetic CRS to CRS

>>> proj = Transformer.from_crs(crs.geodetic_crs, crs)
>>> proj
<Conversion Transformer: pipeline>
Popular Visualisation Pseudo-Mercator
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
>>> proj.transform(12, 15)
(1669792.3618991035, 1345708.4084091093)

4D Transformations with Time

Note

If you are doing a transformation with a CRS that is time based, it is recommended to include the time in the transformation operation.

>>> transformer = Transformer.from_crs(7789, 8401)
>>> transformer
<Transformation Transformer: helmert>
ITRF2014 to ETRF2014 (1)
>>> transformer.transform(xx=3496737.2679, yy=743254.4507, zz=5264462.9620, tt=2019.0)
(3496737.757717311, 743253.9940103051, 5264462.701132784, 2019.0)

Geodesic calculations

Geodesic calculations are useful if you need to calculate the distance between two points or the area of a geometry on Earth’s surface without projecting it onto a plane.

For more examples of usage and documentation, see pyproj.Geod.

Creating Geod class

This example demonstrates creating a pyproj.Geod using an ellipsoid name as well as deriving one using a pyproj.crs.CRS.

>>> from pyproj import CRS, Geod
>>> geod_clrk = Geod(ellps='clrk66') # Use Clarke 1866 ellipsoid.
>>> geod_clrk
Geod(ellps='clrk66')
>>> geod_wgs84 = CRS("EPSG:4326").get_geod()
>>> geod_wgs84
Geod('+a=6378137 +f=0.0033528106647475126')

Geodesic line length

Calculate the geodesic length of a line (See: pyproj.Geod.line_length()):

>>> from pyproj import Geod
>>> lats = [-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7,
...         -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7]
>>> lons = [-74, -102, -102, -131, -163, 163, 172, 140, 113,
...         88, 59, 25, -4, -14, -33, -46, -61]
>>> geod = Geod(ellps="WGS84")
>>> total_length = geod.line_length(lons, lats)
>>> f"{total_length:.3f}"
'14259605.611'

Calculate the geodesic length of a shapely geometry (See: pyproj.Geod.geometry_length()):

>>> from pyproj import Geod
>>> from shapely.geometry import Point, LineString
>>> line_string = LineString([Point(1, 2), Point(3, 4)]))
>>> geod = Geod(ellps="WGS84")
>>> total_length = geod.geometry_length(line_string)
>>> f"{total_length:.3f}"
'313588.397'

Geodesic area

Calculate the geodesic area and perimeter of a polygon (See: pyproj.Geod.polygon_area_perimeter()):

>>> from pyproj import Geod
>>> geod = Geod('+a=6378137 +f=0.0033528106647475126')
>>> lats = [-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7,
...         -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7]
>>> lons = [-74, -102, -102, -131, -163, 163, 172, 140, 113,
...         88, 59, 25, -4, -14, -33, -46, -61]
>>> poly_area, poly_perimeter = geod.polygon_area_perimeter(lons, lats)
>>> f"{poly_area:.3f} {poly_perimeter:.3f}"
'13376856682207.406 14710425.407'

Calculate the geodesic area and perimeter of a shapely polygon (See: pyproj.Geod.geometry_area_perimeter()):

>>> from pyproj import Geod
>>> from shapely.geometry import LineString, Point, Polygon
>>> geod = Geod('+a=6378137 +f=0.0033528106647475126')
>>> poly_area, poly_perimeter = geod.geometry_area_perimeter(
        Polygon(
            LineString([Point(1, 1), Point(1, 10), Point(10, 10), Point(10, 1)]),
            holes=[LineString([Point(1, 2), Point(3, 4), Point(5, 2)])],
        )
    )
>>> f"{poly_area:.3f} {poly_perimeter:.3f}"
'-944373881400.339 3979008.036'