Quickstart

This page runs during every documentation build, so the output below is what the code actually produces with PROJ 9.9.0.

Transform one point

Transform a point in Oslo from ETRS89 geographic coordinates (EPSG:4258) to ETRS89 / UTM zone 32N (EPSG:25832). Input is longitude, latitude in degrees. Output is easting, northing in metres.

from geodetic_engine.geodesy import transform

result = transform("EPSG:4258", "EPSG:25832", (10.7522, 59.9139))
result.coordinates

((597979.9028839065, 6643118.991370024),)

Both CRSs use the same datum, so this is a map projection and nothing had to be chosen. Coordinate values are always given in xy order: longitude before latitude, easting before northing, then height. This holds even though EPSG declares EPSG:4258 latitude-first. See Axis order.

Change datum, and say how

A datum change can be done more than one way. PROJ’s EPSG dataset holds dozens of transformations from ED50 to WGS 84, and they disagree by metres. This package will not pick one for you:

transform("EPSG:4230", "EPSG:4326", (2.5, 63.5))

AmbiguousOperationError: EPSG:4230 to EPSG:4326 involves a datum change; every datum operation must be named explicitly. allow_any_operation no longer permits automatic selection or ballpark results

Name the operation instead. EPSG:1612 is ED50 to WGS 84 (23), a seven-parameter Helmert for Norwegian waters north of 62°N:

result = transform("EPSG:4230", "EPSG:4326", (2.5, 63.5), operation="EPSG:1612")
result.coordinates

((2.4981894757094416, 63.49961374934551),)

Read how the answer was produced

Every result records which operation was applied, how it was chosen, its stated accuracy, and the PROJ pipeline that ran:

op = result.operation
print("operation :", op.authority_code, "-", op.name)
print("method    :", op.method_name)
print("accuracy  :", op.accuracy, "m")
print("route     :", op.route)
print("values    :", result.coordinate_order, result.target_crs.value_axis_abbreviations)
print("declared  :", result.target_axes, result.target_units)

operation : EPSG:1612 - ED50 to WGS 84 (23)
method    : Position Vector transformation (geog2D domain)
accuracy  : 1.0 m
route     : transformer_group
values    : xy ('Lon', 'Lat')
declared  : ('Lat', 'Lon') ('degree', 'degree')

The values are xy (Lon, Lat). target_axes gives the order EPSG declares for the target CRS, latitude first, so you can see that the two differ.

print(result.pipeline)

proj=pipeline step proj=unitconvert xy_in=deg xy_out=rad step proj=push v_3 step proj=cart ellps=intl step proj=helmert x=-116.641 y=-56.931 z=-110.559 rx=0.893 ry=0.921 rz=-0.917 s=-3.52 convention=position_vector step inv proj=cart ellps=WGS84 step proj=pop v_3 step proj=unitconvert xy_in=rad xy_out=deg

The pipeline can be given to pyproj.Transformer.from_pipeline to reproduce the result without this package. result.to_json() serialises all of this, including SHA-256 fingerprints of the proj.db that answered. See Results and provenance.

Many points, one resolution

To transform many batches, resolve the transformation once and reuse it:

import numpy as np

from geodetic_engine.geodesy import Transformation

tfm = Transformation("EPSG:4258", "EPSG:25832", operation="EPSG:16032")
points = np.array([[10.75, 59.91], [5.32, 60.39], [7.99, 58.15]])  # lon, lat
tfm.transform(points).coordinates.to_dataframe()
E N
0 597868.381065 6.642682e+06
1 297230.220207 6.700510e+06
2 440550.920462 6.445855e+06

tfm.transform accepts a single point, a list of tuples, a 2D numpy array, or separate x, y[, z] sequences as pyproj does. See Input and output formats.

Where next

  • Core concepts: the vocabulary used in the rest of the documentation.

  • Geodesy: everything the geodesy module does.

  • Examples: longer worked examples.