geodetic-engine – Quickstart Examples¶
This notebook shows how to use the geodetic_engine.geodesy package: a pyproj wrapper that resolves coordinate reference systems, applies a specific EPSG coordinate operation, and returns coordinates together with their provenance.
Two conventions apply throughout:
Declared axis order is what the EPSG dataset says a CRS’s axes are.
EPSG:4326declares(Lat, Lon).Coordinate value order is always
xy: longitude/easting first, then latitude/northing, then height. The two are reported separately; see Section 3.
A transformation never returns a ballpark approximation, silently substitutes a different requested EPSG operation, or drops a required coordinate epoch. These conditions raise specific exceptions.
Table of Contents¶
1. Setup¶
The package resolves CRSs and builds PROJ transformers lazily, so importing it does no PROJ work by itself.
from geodetic_engine.geodesy import (
CoordinateReferenceSystem,
Transformation,
transform,
)
print("geodetic_engine.geodesy is ready")
geodetic_engine.geodesy is ready
2. Four Quick Examples¶
Four different kinds of transformation, each exercising a different part of the package: a same-datum conversion, a named datum shift, a vertical transformation through a grid, and an epoch-dependent transformation between dynamic reference frames.
2.1 One-Shot Conversion (Same Datum)¶
EPSG:4326 (WGS 84, geographic) to EPSG:3395 (WGS 84 / World Mercator, projected). Both CRSs share the WGS 84 datum, so this is a projection and unit change only – no coordinate operation needs to be named, and PROJ’s choice is recorded rather than assumed.
transform() is the one-shot entry point: it resolves the transformation, applies it, and returns immediately.
oslo = (10.7522, 59.9139) # (longitude, latitude) -- values are always xy
result = transform("EPSG:4326", "EPSG:3395", [oslo])
print("coordinates (xy):", result.coordinates)
print("operation route :", result.operation.route)
print("operation name :", result.operation.name)
print("target axes :", result.target_axes, result.target_units)
coordinates (xy): ((1196929.428907436, 8343586.513829198),)
operation route : proj_default
operation name : World Mercator
target axes : ('E', 'N') ('metre', 'metre')
2.2 Named Datum Shift¶
EPSG:4230 (ED50) to EPSG:4326 (WGS 84): this is a datum change, so the coordinate operation to apply must be named. Several valid ED50-to-WGS84 shifts exist with different accuracy and area of use; geodetic_engine will not pick one on your behalf, and once you name one it verifies that PROJ actually applied it rather than a different candidate.
Omitting operation here raises AmbiguousOperationError. The one exception is a CRS that names its own transformation – see Section 7.
datum_shift = Transformation("EPSG:4230", "EPSG:4326", operation="EPSG:1133")
result = datum_shift.transform([oslo]) # oslo, expressed as an ED50 lon/lat
print("coordinates (xy) :", result.coordinates)
print("applied operation:", datum_shift.operation.authority_code)
print("method :", datum_shift.operation.method_name)
coordinates (xy) : ((10.750769179490943, 59.91344872175593),)
applied operation: EPSG:1133
method : Geocentric translations (geog2D domain)
2.3 Vertical (Geoid) Height Transformation¶
EPSG:4979 (WGS 84, 3D geographic) to EPSG:3855 (EGM2008 height): converts an ellipsoidal height to an orthometric one through a geoid grid. The target CRS declares a single height axis, so the result carries exactly one value per point even though PROJ computes three coordinates internally. The grid the operation depends on is reported on the Transformation, not just used silently.
geoid = Transformation("EPSG:4979", "EPSG:3855", operation="EPSG:3858")
result = geoid.transform([(-144.0, 72.0, 548.4082)]) # lon, lat, ellipsoidal height (m)
print("orthometric height (m):", result.coordinates)
print("target axes :", result.target_axes, result.target_units)
print("grids used :", [(grid.name, grid.available) for grid in geoid.grids])
orthometric height (m): ((556.3834421379089,),)
target axes : ('H',) ('metre',)
grids used : [('us_nga_egm08_25.tif', True)]
2.4 Dynamic Reference Frame With a Coordinate Epoch¶
EPSG:4896 (ITRF2008, geocentric, dynamic) to EPSG:4938 (GDA94, geocentric): the operation reads rates of change, so a coordinate epoch is required. Without it the result would silently be displaced by the motion between the true and assumed epochs – Transformation.requires_epoch tells you this before you transform anything.
dynamic = Transformation("EPSG:4896", "EPSG:4938", operation="EPSG:6277")
print("requires a coordinate epoch:", dynamic.requires_epoch)
result = dynamic.transform(
[(-2593197.524, 5656917.6189, -1394397.8828)],
coordinate_epoch=1993.0, # decimal year the coordinates were observed at
)
print("coordinates (xy):", result.coordinates)
requires a coordinate epoch: True
coordinates (xy): ((-2593197.589192494, 5656917.670923312, -1394397.8240364394),)
3. Inspecting a CRS’s Declared Axes and Units¶
CoordinateReferenceSystem reports the axis roles, directions and units the EPSG dataset declares – axis_abbreviations and axis_units are always in declared order, (Lat, Lon) for EPSG:4326, regardless of the xy value convention.
value_axis_order is the bridge between the two: it gives, for each position in a coordinate value, which declared axis it corresponds to. value_axis_abbreviations renders that as labels.
crs = CoordinateReferenceSystem.from_user_input("EPSG:4326")
print("declared axes :", crs.axis_abbreviations) # EPSG order: latitude first
print("declared units :", crs.axis_units)
print("dimension :", crs.dimension)
print("value axis order :", crs.value_axis_order)
print("value abbrevs :", crs.value_axis_abbreviations) # ('Lon', 'Lat')
print()
for axis in crs.axes:
print(
f" {axis.abbrev:>4} {axis.name:<20} "
f"direction={axis.direction:<6} unit={axis.unit_name}"
)
declared axes : ('Lat', 'Lon')
declared units : ('degree', 'degree')
dimension : 2
value axis order : (1, 0)
value abbrevs : ('Lon', 'Lat')
Lat Geodetic latitude direction=north unit=degree
Lon Geodetic longitude direction=east unit=degree
4. Batch Transforms¶
Coordinates are passed as a plain list of points – a list of tuples, a list of lists, or a 2D numpy array – and reshaped one list per axis internally, so a whole batch crosses into PROJ in a single call. There is no wrapper type to construct first.
cities = {
"Oslo": (10.7522, 59.9139),
"Bergen": (5.3221, 60.3913),
"Tromso": (18.9553, 69.6492),
}
result = Transformation("EPSG:4326", "EPSG:3395").transform(list(cities.values()))
for name, xy in zip(cities, result.coordinates, strict=True):
print(f"{name:8}-> {xy}")
Oslo -> (1196929.428907436, 8343586.513829198)
Bergen -> (592453.4619508812, 8450191.559971118)
Tromso -> (2110094.3438337385, 10915376.196796743)
4.1 NumPy Arrays Work Too¶
A 2D array of shape (n_points, n_axes) – one row per point, values in xy order – works everywhere a sequence of points is accepted: transform() and Transformation.transform(). Nothing needs to be converted to a list of tuples first.
A row may also carry one value more than the CRS declares: a height alongside a 2D horizontal CRS. PROJ carries it through unchanged rather than consuming it, matching pyproj.Transformer.transform(xx, yy, zz). Two extra values, or too few, are rejected.
import numpy as np
rows = np.array([[10, 60], [11, 61]]) # integer dtype works too
print(rows.dtype, rows.shape)
result = transform("EPSG:4326", "EPSG:3395", rows)
print(result.coordinates)
# a height alongside the 2D pair passes through unchanged
with_height = np.array([[10.0, 60.0, 100.0], [11.0, 61.0, 200.0]])
result_with_height = transform("EPSG:4326", "EPSG:3395", with_height)
print(result_with_height.coordinates)
int64 (2, 2)
((1113194.9079327357, 8362698.548500747), (1224514.3987260093, 8588415.031895077))
((1113194.9079327357, 8362698.548500747, 100.0), (1224514.3987260093, 8588415.031895077, 200.0))
5. Reading the Provenance of a Result¶
A TransformationResult is not a bare tuple of floats. to_json_dict() renders every provenance field alongside the coordinates: which operation was applied, its accuracy, the grids it used, and the axis/unit contract for both CRSs – everything needed to answer, after the fact, what produced this number. to_json() renders the same fields as a JSON string, indented by default.
result = transform(
"EPSG:4979", "EPSG:3855", [(-144.0, 72.0, 548.4082)], operation="EPSG:3858"
)
print(result.to_json())
result
{
"coordinates": [
[
556.3834421379089
]
],
"coordinate_order": "xy",
"coordinate_epoch": null,
"source_crs": "EPSG:4979",
"target_crs": "EPSG:3855",
"source_crs_wkt": "GEOGCRS[\"WGS 84\",ENSEMBLE[\"World Geodetic System 1984 ensemble\",MEMBER[\"World Geodetic System 1984 (Transit)\"],MEMBER[\"World Geodetic System 1984 (G730)\"],MEMBER[\"World Geodetic System 1984 (G873)\"],MEMBER[\"World Geodetic System 1984 (G1150)\"],MEMBER[\"World Geodetic System 1984 (G1674)\"],MEMBER[\"World Geodetic System 1984 (G1762)\"],MEMBER[\"World Geodetic System 1984 (G2139)\"],MEMBER[\"World Geodetic System 1984 (G2296)\"],ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ENSEMBLEACCURACY[2.0]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],CS[ellipsoidal,3],AXIS[\"geodetic latitude (Lat)\",north,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433]],AXIS[\"geodetic longitude (Lon)\",east,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433]],AXIS[\"ellipsoidal height (h)\",up,ORDER[3],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Geodesy. Navigation and positioning using GPS satellite system.\"],AREA[\"World.\"],BBOX[-90,-180,90,180]],ID[\"EPSG\",4979]]",
"target_crs_wkt": "VERTCRS[\"EGM2008 height\",VDATUM[\"EGM2008 geoid\"],CS[vertical,1],AXIS[\"gravity-related height (H)\",up,LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Geodesy.\"],AREA[\"World.\"],BBOX[-90,-180,90,180]],ID[\"EPSG\",3855]]",
"source_axes": [
"Lat",
"Lon",
"h"
],
"source_units": [
"degree",
"degree",
"metre"
],
"target_axes": [
"H"
],
"target_units": [
"metre"
],
"operation": {
"requested": "EPSG:3858",
"applied": "EPSG:3858",
"name": "WGS 84 to EGM2008 height (1)",
"method": "Geographic3D to GravityRelatedHeight (EGM2008)",
"accuracy_m": 0.113,
"route": "transformer_group",
"steps": [
"WGS 84 to EGM2008 height (1)",
"axis order change (2D)"
],
"ballpark": false,
"requires_epoch": false,
"execution_direction": "FORWARD",
"bound_operations": [],
"axis_order_corrected": false,
"definition": {
"type": "Transformation",
"name": "WGS 84 to EGM2008 height (1)",
"source_crs": {
"type": "GeographicCRS",
"name": "WGS 84",
"datum_ensemble": {
"name": "World Geodetic System 1984 ensemble",
"members": [
{
"name": "World Geodetic System 1984 (Transit)",
"id": {
"authority": "EPSG",
"code": 1166
}
},
{
"name": "World Geodetic System 1984 (G730)",
"id": {
"authority": "EPSG",
"code": 1152
}
},
{
"name": "World Geodetic System 1984 (G873)",
"id": {
"authority": "EPSG",
"code": 1153
}
},
{
"name": "World Geodetic System 1984 (G1150)",
"id": {
"authority": "EPSG",
"code": 1154
}
},
{
"name": "World Geodetic System 1984 (G1674)",
"id": {
"authority": "EPSG",
"code": 1155
}
},
{
"name": "World Geodetic System 1984 (G1762)",
"id": {
"authority": "EPSG",
"code": 1156
}
},
{
"name": "World Geodetic System 1984 (G2139)",
"id": {
"authority": "EPSG",
"code": 1309
}
},
{
"name": "World Geodetic System 1984 (G2296)",
"id": {
"authority": "EPSG",
"code": 1383
}
}
],
"ellipsoid": {
"name": "WGS 84",
"semi_major_axis": 6378137,
"inverse_flattening": 298.257223563
},
"accuracy": "2.0",
"id": {
"authority": "EPSG",
"code": 6326
}
},
"coordinate_system": {
"subtype": "ellipsoidal",
"axis": [
{
"name": "Geodetic latitude",
"abbreviation": "Lat",
"direction": "north",
"unit": "degree"
},
{
"name": "Geodetic longitude",
"abbreviation": "Lon",
"direction": "east",
"unit": "degree"
},
{
"name": "Ellipsoidal height",
"abbreviation": "h",
"direction": "up",
"unit": "metre"
}
]
},
"id": {
"authority": "EPSG",
"code": 4979
}
},
"target_crs": {
"type": "VerticalCRS",
"name": "EGM2008 height",
"datum": {
"type": "VerticalReferenceFrame",
"name": "EGM2008 geoid"
},
"coordinate_system": {
"subtype": "vertical",
"axis": [
{
"name": "Gravity-related height",
"abbreviation": "H",
"direction": "up",
"unit": "metre"
}
]
},
"id": {
"authority": "EPSG",
"code": 3855
}
},
"method": {
"name": "Geographic3D to GravityRelatedHeight (EGM2008)",
"id": {
"authority": "EPSG",
"code": 1025
}
},
"parameters": [
{
"name": "Geoid (height correction) model file",
"value": "us_nga_egm08_25.tif",
"id": {
"authority": "EPSG",
"code": 8666
}
}
],
"accuracy": "0.113",
"id": {
"authority": "DERIVED_FROM(EPSG)",
"code": 3858
},
"remarks": "Replaces WGS 84 to EGM96 height (1) (CT code 10084). Grid spacing is 2.5 arc-minutes. For a file with smaller (1 arc-minute) grid spacing (in principle more exactly representing the spherical harmonics) see CT code 3859. For reversible alternative see CT code 9704. An executable using spherical harmonics is available from the NGA site."
}
},
"grids": [
{
"name": "us_nga_egm08_25.tif",
"available": true,
"package": "",
"url": "https://cdn.proj.org/us_nga_egm08_25.tif",
"full_name": "/usr/local/share/proj/us_nga_egm08_25.tif",
"open_license": true,
"direct_download": true
}
],
"pipeline": "proj=pipeline step proj=unitconvert xy_in=deg xy_out=rad step inv proj=vgridshift grids=us_nga_egm08_25.tif multiplier=1 step proj=unitconvert xy_in=rad xy_out=deg step proj=axisswap order=2,1",
"database_fingerprints": {
"/usr/local/share/proj/proj.db": "3d3d8b279d100544b5737cd08a732335170fff09ee91166b4137919d6d019198"
}
}
TransformationResult(coordinates=((556.3834421379089,),), source_crs=CoordinateReferenceSystem('EPSG:4979'), target_crs=CoordinateReferenceSystem('EPSG:3855'), operation=AppliedOperation(requested='EPSG:3858', auth_name='EPSG', code='3858', name='WGS 84 to EGM2008 height (1)', method_name='Geographic3D to GravityRelatedHeight (EGM2008)', accuracy=0.113, route=<OperationRoute.TRANSFORMER_GROUP: 'transformer_group'>, ballpark=False, requires_epoch=False, steps=('WGS 84 to EGM2008 height (1)', 'axis order change (2D)'), execution_direction=<TransformDirection.FORWARD: 'FORWARD'>, bound_operations=(), axis_order_corrected=False), grids=(GridUsage(name='us_nga_egm08_25.tif', full_name='/usr/local/share/proj/us_nga_egm08_25.tif', package_name='', url='https://cdn.proj.org/us_nga_egm08_25.tif', available=True, open_license=True, direct_download=True),), coordinate_epoch=None, coordinate_order='xy', pipeline='proj=pipeline step proj=unitconvert xy_in=deg xy_out=rad step inv proj=vgridshift grids=us_nga_egm08_25.tif multiplier=1 step proj=unitconvert xy_in=rad xy_out=deg step proj=axisswap order=2,1', database_fingerprints=(('/usr/local/share/proj/proj.db', '3d3d8b279d100544b5737cd08a732335170fff09ee91166b4137919d6d019198'),))
6. Reusing a Transformation for Many Batches¶
Resolving a transformation – finding the operation, verifying it, building the PROJ transformer – happens once, in Transformation.__init__. Prefer building one Transformation and calling .transform() on it repeatedly over calling the one-shot transform() function in a loop, which re-resolves every time it sees a CRS pair for the first time.
mercator = Transformation("EPSG:4326", "EPSG:3395")
for name, lonlat in cities.items():
print(name, mercator.transform([lonlat]).coordinates)
Oslo ((1196929.428907436, 8343586.513829198),)
Bergen ((592453.4619508812, 8450191.559971118),)
Tromso ((2110094.3438337385, 10915376.196796743),)
7. Bound CRSs: a CRS That Names Its Own Transformation¶
Section 2.2 showed that a datum change without a named operation is refused. A bound CRS is the exception that proves the rule rather than an escape from it: it is a CRS packaged together with the single transformation that ties it to a hub, usually WGS 84. The operation is part of the CRS definition, so it was named – by whoever defined the CRS – and PROJ is left with exactly one candidate.
This is early binding. Transformation applies it without complaint and reports route = bound, alongside the EPSG code of the operation that was embedded, so the result is exactly as auditable as one where you named the operation yourself.
Bound CRSs are what geodetic-projdb writes into text_definition on the geodetic_crs and projected_crs tables when your register defines them. See the README for how they are built and stored.
from pyproj import CRS as PyprojCRS
from pyproj.crs import BoundCRS, CoordinateOperation
ed50_to_wgs84 = CoordinateOperation.from_authority("EPSG", 1133)
# ED50 bound to WGS 84 through EPSG:1133, as a register would define it. The hub
# is the transformation's own target, so WGS 84 geographic rather than a UTM CRS.
ed50_bound = BoundCRS(
PyprojCRS.from_epsg(4230), PyprojCRS.from_epsg(4326), ed50_to_wgs84
)
bound = Transformation(ed50_bound, "EPSG:4326")
print("route :", bound.operation.route)
print("operation :", bound.operation.authority_code, "|", bound.operation.name)
print("requested :", bound.operation.requested, "(the CRS named it, not the caller)")
result = bound.transform([oslo])
print("coordinates:", result.coordinates)
# Naming the same operation explicitly must give the same answer.
named = Transformation("EPSG:4230", "EPSG:4326", operation="EPSG:1133")
agrees = result.coordinates == named.transform([oslo]).coordinates
print("agrees with the named operation:", agrees)
route : bound
operation : EPSG:1133 | ED50 to WGS 84 (1)
requested : None (the CRS named it, not the caller)
coordinates: ((10.750769179490943, 59.91344872175593),)
agrees with the named operation: True
A bound CRS is resolved through the same transformer group as a named operation, rather than being handed to PROJ to build a transformer by itself. That matters because most bound CRSs in a real register have a projected base: ED50 / UTM zone 32N bound to WGS 84, not just ED50. Such a transformation has to unproject, apply the datum shift, and reproject, and going through the group is what supplies those steps while keeping the embedded operation identifiable.
# ED50 / UTM zone 32N bound to WGS 84 -- a projected base, the common real case.
utm_bound = BoundCRS(
PyprojCRS.from_epsg(23032), PyprojCRS.from_epsg(4326), ed50_to_wgs84
)
statfjord = (600000.0, 6643000.0) # ED50 / UTM 32N, offshore Norway
for target in ("EPSG:4326", "EPSG:32632"):
step = Transformation(utm_bound, target)
print(
f"-> {target}: route={step.operation.route} "
f"op={step.operation.authority_code} "
f"coords={step.transform([statfjord]).coordinates[0]}"
)
# The datum shift is still EPSG:1133; the projection steps are wrapped around it.
projected = Transformation(utm_bound, "EPSG:32632")
print("\npipeline steps:", projected.operation.steps)
-> EPSG:4326: route=bound op=EPSG:1133 coords=(10.786654370599862, 59.91050486273897)
-> EPSG:32632: route=bound op=EPSG:1133 coords=(599916.4249825485, 6642792.470067851)
pipeline steps: ('ED50 to WGS 84 (1)', 'Inverse of UTM zone 32N', 'UTM zone 32N')
7.1 When the Transformation Is a Chain¶
PROJ cannot embed a chain of operations in a bound CRS. EPSG:8047 (ED50 to WGS 84 (15)) is two Helmert steps through ED87, so a register defining a bound CRS over it would produce something PROJ refuses to build.
collapse_concatenated rewrites such a chain as one equivalent Helmert. Two Helmerts compose exactly, because each is an affine map on geocentric coordinates and the intermediate conversions cancel. But the algebra is only the proposal: EPSG linearises the rotation matrix, and states rotations in different units for different operations (EPSG:1147 uses microradians, not arc-seconds). So every collapse is verified against PROJ’s own rendering of the original chain and refused if it moves a coordinate by more than a millimetre.
A chain that is genuinely not one Helmert – a grid step, Molodensky-Badekas, time-dependent – is not approximated. It raises NotCollapsibleError.
from geodetic_engine.geodesy import NotCollapsibleError
from geodetic_engine.geodesy.utils import collapse_concatenated
chain = CoordinateOperation.from_authority("EPSG", 8047)
steps = chain.to_json_dict()["steps"]
print("chain :", chain.name)
for number, step in enumerate(steps, start=1):
print(f" step {number}: {step['name']} ({step['method']['name']})")
# PROJ will not embed a chain in a bound CRS.
try:
BoundCRS(PyprojCRS.from_epsg(4230), PyprojCRS.from_epsg(4326), chain)
except Exception as error:
print("\nembedding the chain directly:", type(error).__name__)
collapsed = collapse_concatenated(chain)
print("collapsed:", collapsed.name)
print("method :", collapsed.method_name)
print("params :", [round(value, 6) for value in collapsed.towgs84])
bound_chain = BoundCRS(PyprojCRS.from_epsg(4230), PyprojCRS.from_epsg(4326), collapsed)
over_chain = Transformation(bound_chain, "EPSG:4326")
print("coordinates:", over_chain.transform([oslo]).coordinates)
# A chain that is not equivalent to a single Helmert is refused, not approximated.
try:
collapse_concatenated(CoordinateOperation.from_authority("EPSG", 3896))
except NotCollapsibleError as error:
print("\nNotCollapsibleError:", error)
chain : ED50 to WGS 84 (15)
step 1: ED50 to ED87 (2) (Position Vector transformation (geog2D domain))
step 2: ED87 to WGS 84 (1) (Position Vector transformation (geog2D domain))
embedding the chain directly: CRSError
collapsed: ED50 to WGS 84 (15) (collapsed to a single step)
method : Position Vector transformation (geog2D domain)
params : [-84.491, -100.559002, -114.208998, -0.495159, -0.110703, -0.489714, 0.2947]
coordinates: ((10.750821356614876, 59.913468329177206),)
NotCollapsibleError: EPSG:3896 is a concatenated operation, and a bound CRS can carry only a single transformation, so its steps have to compose into one Helmert; step 1 applies 'Longitude rotation', which is not a Helmert and so does not compose
7.2 OSDU Bound CRS References¶
An OSDU persistableReference carries its CRS and transformation definitions together. It can be read without a custom catalogue database.
This example creates an illustrative early-bound payload from the stock EPSG definitions, then reads and applies the payload. ED50 is bound to WGS 84 through EPSG:1133. No producer authority code is invented during export.
from pyproj import CRS as PyprojCRS
from pyproj.crs import BoundCRS, CoordinateOperation
from geodetic_engine.geodesy import CoordinateReferenceSystem, Transformation
from geodetic_engine.persistablereference import (
parse_persistable_reference,
to_persistable_reference,
)
bound = BoundCRS(
PyprojCRS.from_epsg(4230),
PyprojCRS.from_epsg(4326),
CoordinateOperation.from_authority("EPSG", "1133"),
)
source_payload = to_persistable_reference(bound)
source_reference = parse_persistable_reference(source_payload)
print(source_reference.kind, source_reference.name)
print(source_reference.operation.method_names)
EBC ED50
('Geocentric_Translation',)
points = (10.0, 60.0)
trg_crs = "EPSG:4326"
crs = CoordinateReferenceSystem.from_persistable_reference(source_payload)
ct = Transformation(crs, trg_crs)
print(f"OSDU reference: {source_reference.name}")
print(f" base={crs.crs.name} bound={crs.crs.is_bound}")
print(f" axes={crs.axis_abbreviations} route={ct.operation.route}")
print(f" applied={ct.operation.name}")
print(f" {points} -> {ct.transform([points]).coordinates[0]}")
ct.operation
OSDU reference: ED50
base=ED50 bound=True
axes=('lon', 'lat') route=bound
applied=ED50 to WGS 84 (1)
(10.0, 60.0) -> (9.998541177251479, 59.999543827718284)
AppliedOperation(requested=None, auth_name=None, code=None, name='ED50 to WGS 84 (1)', method_name='Geocentric translations (geog2D domain)', accuracy=None, route=<OperationRoute.BOUND: 'bound'>, ballpark=False, requires_epoch=False, steps=('ED50 to WGS 84 (1)', 'axis order change (2D)'), execution_direction=<TransformDirection.FORWARD: 'FORWARD'>, bound_operations=('ED50 to WGS 84 (1)',), axis_order_corrected=False)
Exporting the WKT¶
Both halves export as WKT, and they are not the same object.
crs.crs.to_wkt() gives the BOUNDCRS exactly as it sits in text_definition: the base CRS under SOURCECRS, the hub under TARGETCRS, and the datum shift as an ABRIDGEDTRANSFORMATION.
operation.to_wkt() gives just the shift, as a full COORDINATEOPERATION. It renders from what PROJ actually built, not from a lookup by code, so it also works for an operation EPSG does not define – the collapsed chain in 7.1 has no authority code and could not be fetched back any other way. The trade is that the registry’s VERSION, USAGE and REMARK are not carried: take parameters from here, take scope and area of validity from the EPSG dataset via authority_code.
Abridged versus full is a real distinction rather than a formatting one. Compare Scale difference in the two outputs: 1.0000012 in the BOUNDCRS, against 1.2 with SCALEUNIT["parts per million"] in the operation. The abridged form carries no units at all – metres and arc-seconds are implied, and scale is stored as \(1 + s \cdot 10^{-6}\). Same quantity, two conventions, and reading one as the other is a 1.2 ppm error waiting to happen.
print(crs.crs.to_wkt(pretty=True))
# The shift on its own, rendered from what PROJ built rather than looked up.
print("\n" + (ct.operation.to_wkt(pretty=True) or "no WKT available"))
BOUNDCRS[
SOURCECRS[
GEOGCRS["ED50",
DATUM["European Datum 1950",
ELLIPSOID["International 1924",6378388,297,
LENGTHUNIT["metre",1]],
ID["EPSG",6230]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]],
CS[ellipsoidal,2],
AXIS["longitude",east,
ORDER[1],
ANGLEUNIT["Degree",0.0174532925199433]],
AXIS["latitude",north,
ORDER[2],
ANGLEUNIT["Degree",0.0174532925199433]]]],
TARGETCRS[
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
MEMBER["World Geodetic System 1984 (G2296)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]],
ABRIDGEDTRANSFORMATION["ED50 to WGS 84 (1)",
METHOD["Geocentric translations (geog2D domain)",
ID["EPSG",9603]],
PARAMETER["X-axis translation",-87,
ID["EPSG",8605]],
PARAMETER["Y-axis translation",-98,
ID["EPSG",8606]],
PARAMETER["Z-axis translation",-121,
ID["EPSG",8607]]]]
COORDINATEOPERATION["ED50 to WGS 84 (1)",
SOURCECRS[
GEOGCRS["ED50",
DATUM["European Datum 1950",
ELLIPSOID["International 1924",6378388,297,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4230]]],
TARGETCRS[
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
MEMBER["World Geodetic System 1984 (G2296)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]],
METHOD["Geocentric translations (geog2D domain)",
ID["EPSG",9603]],
PARAMETER["X-axis translation",-87,
LENGTHUNIT["metre",1],
ID["EPSG",8605]],
PARAMETER["Y-axis translation",-98,
LENGTHUNIT["metre",1],
ID["EPSG",8606]],
PARAMETER["Z-axis translation",-121,
LENGTHUNIT["metre",1],
ID["EPSG",8607]]]