geodetic-engine examples¶
Examples for CRS inspection, explicitly selected coordinate transformations, operation provenance, and coordinate export.
The examples use the standard PROJ database and do not require this repository’s built proj.db. Sections 2.1 and 2.2 also require the Norwegian height grid no_kv_href2008a.tif, which is installed in this repository’s dev container but is not bundled with the standard database. Outside the container, install it explicitly with python -m pyproj sync --file no_kv_href2008a.tif before running those sections. Without it, those transformations raise MissingGridError; they do not substitute a less accurate operation.
Contents¶
1. Horizontal transformations¶
These examples transform geographic and projected coordinates while keeping the applied operation explicit and traceable.
1.1 ED50 to WGS 84 with an explicit operation¶
Transform an ED50 geographic point to WGS 84 using a named EPSG operation. The result includes the transformed coordinates, operation provenance, JSON representations, and WKT.
import numpy as np
from geodetic_engine.geodesy import Transformation
src_crs = "EPSG:4230"
trg_crs = "EPSG:4326"
operation = "EPSG:1612"
longitude, latitude, height = 2.5, 63.5, 100
ct = Transformation(source_crs=src_crs, target_crs=trg_crs, operation=operation)
results = ct.transform(np.array([longitude, latitude, height]))
trans_lon, trans_lat, trans_h = results.coordinates[0]
print(trans_lon, trans_lat, trans_h)
print(results.to_json())
print(results.to_json(pretty=False))
print(results.operation.to_wkt())
2.4981894919143337 63.49961375156072 100.0
{
"coordinates": [
[
2.4981894919143337,
63.49961375156072,
100.0
]
],
"coordinate_order": "xy",
"coordinate_epoch": null,
"source_crs": "EPSG:4230",
"target_crs": "EPSG:4326",
"source_crs_wkt": "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]],USAGE[SCOPE[\"Geodesy.\"],AREA[\"Europe - west: Andorra; Cyprus; Denmark - onshore and offshore; Faroe Islands - onshore; France - offshore; Germany - offshore North Sea; Gibraltar; Greece - offshore; Israel - offshore; Italy including San Marino and Vatican City State; Ireland offshore; Malta; Netherlands - offshore; North Sea; Norway including Svalbard - onshore and offshore; Portugal - mainland - offshore; Spain - onshore; T\u00fcrkiye (Turkey) - onshore and offshore; United Kingdom - UKCS offshore east of 6\u00b0W including Channel Islands (Guernsey and Jersey). Egypt - Western Desert; Iraq - onshore; Jordan.\"],BBOX[25.71,-16.1,84.73,48.61]],ID[\"EPSG\",4230]]",
"target_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,2],AXIS[\"geodetic latitude (Lat)\",north,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433]],AXIS[\"geodetic longitude (Lon)\",east,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433]],USAGE[SCOPE[\"Horizontal component of 3D system.\"],AREA[\"World.\"],BBOX[-90,-180,90,180]],ID[\"EPSG\",4326]]",
"source_axes": [
"Lat",
"Lon"
],
"source_units": [
"degree",
"degree"
],
"target_axes": [
"Lat",
"Lon"
],
"target_units": [
"degree",
"degree"
],
"operation": {
"requested": "EPSG:1612",
"applied": "EPSG:1612",
"name": "ED50 to WGS 84 (23)",
"method": "Position Vector transformation (geog2D domain)",
"accuracy_m": 1.0,
"route": "transformer_group",
"steps": [
"ED50 to WGS 84 (23)",
"axis order change (2D)"
],
"ballpark": false,
"requires_epoch": false,
"execution_direction": "FORWARD",
"bound_operations": [],
"axis_order_corrected": false,
"definition": {
"type": "Transformation",
"name": "ED50 to WGS 84 (23)",
"source_crs": {
"type": "GeographicCRS",
"name": "ED50",
"datum": {
"type": "GeodeticReferenceFrame",
"name": "European Datum 1950",
"ellipsoid": {
"name": "International 1924",
"semi_major_axis": 6378388,
"inverse_flattening": 297
}
},
"coordinate_system": {
"subtype": "ellipsoidal",
"axis": [
{
"name": "Geodetic latitude",
"abbreviation": "Lat",
"direction": "north",
"unit": "degree"
},
{
"name": "Geodetic longitude",
"abbreviation": "Lon",
"direction": "east",
"unit": "degree"
}
]
},
"id": {
"authority": "EPSG",
"code": 4230
}
},
"target_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"
}
]
},
"id": {
"authority": "EPSG",
"code": 4326
}
},
"method": {
"name": "Position Vector transformation (geog2D domain)",
"id": {
"authority": "EPSG",
"code": 9606
}
},
"parameters": [
{
"name": "X-axis translation",
"value": -116.641,
"unit": "metre",
"id": {
"authority": "EPSG",
"code": 8605
}
},
{
"name": "Y-axis translation",
"value": -56.931,
"unit": "metre",
"id": {
"authority": "EPSG",
"code": 8606
}
},
{
"name": "Z-axis translation",
"value": -110.559,
"unit": "metre",
"id": {
"authority": "EPSG",
"code": 8607
}
},
{
"name": "X-axis rotation",
"value": 0.893,
"unit": {
"type": "AngularUnit",
"name": "arc-second",
"conversion_factor": 4.84813681109536e-06
},
"id": {
"authority": "EPSG",
"code": 8608
}
},
{
"name": "Y-axis rotation",
"value": 0.921,
"unit": {
"type": "AngularUnit",
"name": "arc-second",
"conversion_factor": 4.84813681109536e-06
},
"id": {
"authority": "EPSG",
"code": 8609
}
},
{
"name": "Z-axis rotation",
"value": -0.917,
"unit": {
"type": "AngularUnit",
"name": "arc-second",
"conversion_factor": 4.84813681109536e-06
},
"id": {
"authority": "EPSG",
"code": 8610
}
},
{
"name": "Scale difference",
"value": -3.52,
"unit": {
"type": "ScaleUnit",
"name": "parts per million",
"conversion_factor": 1e-06
},
"id": {
"authority": "EPSG",
"code": 8611
}
}
],
"accuracy": "1.0",
"id": {
"authority": "EPSG",
"code": 1612
},
"remarks": "Parameter values are taken from ED50 to ETRS89 (1), code 1588, (rotations converted from \u03bcrad). Adopted for ED50 to WGS 84 CTs offshore Norway north of 62\u00b0N from April 2001, replacing code 1590. Included in Statens Kartverk programme wsktrans from v4.0."
}
},
"grids": [],
"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",
"database_fingerprints": {
"/usr/local/share/proj/proj.db": "3d3d8b279d100544b5737cd08a732335170fff09ee91166b4137919d6d019198"
}
}
{"coordinates": [[2.4981894919143337, 63.49961375156072, 100.0]], "coordinate_order": "xy", "coordinate_epoch": null, "source_crs": "EPSG:4230", "target_crs": "EPSG:4326", "source_crs_wkt": "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]],USAGE[SCOPE[\"Geodesy.\"],AREA[\"Europe - west: Andorra; Cyprus; Denmark - onshore and offshore; Faroe Islands - onshore; France - offshore; Germany - offshore North Sea; Gibraltar; Greece - offshore; Israel - offshore; Italy including San Marino and Vatican City State; Ireland offshore; Malta; Netherlands - offshore; North Sea; Norway including Svalbard - onshore and offshore; Portugal - mainland - offshore; Spain - onshore; T\u00fcrkiye (Turkey) - onshore and offshore; United Kingdom - UKCS offshore east of 6\u00b0W including Channel Islands (Guernsey and Jersey). Egypt - Western Desert; Iraq - onshore; Jordan.\"],BBOX[25.71,-16.1,84.73,48.61]],ID[\"EPSG\",4230]]", "target_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,2],AXIS[\"geodetic latitude (Lat)\",north,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433]],AXIS[\"geodetic longitude (Lon)\",east,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433]],USAGE[SCOPE[\"Horizontal component of 3D system.\"],AREA[\"World.\"],BBOX[-90,-180,90,180]],ID[\"EPSG\",4326]]", "source_axes": ["Lat", "Lon"], "source_units": ["degree", "degree"], "target_axes": ["Lat", "Lon"], "target_units": ["degree", "degree"], "operation": {"requested": "EPSG:1612", "applied": "EPSG:1612", "name": "ED50 to WGS 84 (23)", "method": "Position Vector transformation (geog2D domain)", "accuracy_m": 1.0, "route": "transformer_group", "steps": ["ED50 to WGS 84 (23)", "axis order change (2D)"], "ballpark": false, "requires_epoch": false, "execution_direction": "FORWARD", "bound_operations": [], "axis_order_corrected": false, "definition": {"type": "Transformation", "name": "ED50 to WGS 84 (23)", "source_crs": {"type": "GeographicCRS", "name": "ED50", "datum": {"type": "GeodeticReferenceFrame", "name": "European Datum 1950", "ellipsoid": {"name": "International 1924", "semi_major_axis": 6378388, "inverse_flattening": 297}}, "coordinate_system": {"subtype": "ellipsoidal", "axis": [{"name": "Geodetic latitude", "abbreviation": "Lat", "direction": "north", "unit": "degree"}, {"name": "Geodetic longitude", "abbreviation": "Lon", "direction": "east", "unit": "degree"}]}, "id": {"authority": "EPSG", "code": 4230}}, "target_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"}]}, "id": {"authority": "EPSG", "code": 4326}}, "method": {"name": "Position Vector transformation (geog2D domain)", "id": {"authority": "EPSG", "code": 9606}}, "parameters": [{"name": "X-axis translation", "value": -116.641, "unit": "metre", "id": {"authority": "EPSG", "code": 8605}}, {"name": "Y-axis translation", "value": -56.931, "unit": "metre", "id": {"authority": "EPSG", "code": 8606}}, {"name": "Z-axis translation", "value": -110.559, "unit": "metre", "id": {"authority": "EPSG", "code": 8607}}, {"name": "X-axis rotation", "value": 0.893, "unit": {"type": "AngularUnit", "name": "arc-second", "conversion_factor": 4.84813681109536e-06}, "id": {"authority": "EPSG", "code": 8608}}, {"name": "Y-axis rotation", "value": 0.921, "unit": {"type": "AngularUnit", "name": "arc-second", "conversion_factor": 4.84813681109536e-06}, "id": {"authority": "EPSG", "code": 8609}}, {"name": "Z-axis rotation", "value": -0.917, "unit": {"type": "AngularUnit", "name": "arc-second", "conversion_factor": 4.84813681109536e-06}, "id": {"authority": "EPSG", "code": 8610}}, {"name": "Scale difference", "value": -3.52, "unit": {"type": "ScaleUnit", "name": "parts per million", "conversion_factor": 1e-06}, "id": {"authority": "EPSG", "code": 8611}}], "accuracy": "1.0", "id": {"authority": "EPSG", "code": 1612}, "remarks": "Parameter values are taken from ED50 to ETRS89 (1), code 1588, (rotations converted from \u03bcrad). Adopted for ED50 to WGS 84 CTs offshore Norway north of 62\u00b0N from April 2001, replacing code 1590. Included in Statens Kartverk programme wsktrans from v4.0."}}, "grids": [], "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", "database_fingerprints": {"/usr/local/share/proj/proj.db": "3d3d8b279d100544b5737cd08a732335170fff09ee91166b4137919d6d019198"}}
COORDINATEOPERATION["ED50 to WGS 84 (23)",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["Position Vector transformation (geog2D domain)",ID["EPSG",9606]],PARAMETER["X-axis translation",-116.641,LENGTHUNIT["metre",1],ID["EPSG",8605]],PARAMETER["Y-axis translation",-56.931,LENGTHUNIT["metre",1],ID["EPSG",8606]],PARAMETER["Z-axis translation",-110.559,LENGTHUNIT["metre",1],ID["EPSG",8607]],PARAMETER["X-axis rotation",0.893,ANGLEUNIT["arc-second",4.84813681109536E-06],ID["EPSG",8608]],PARAMETER["Y-axis rotation",0.921,ANGLEUNIT["arc-second",4.84813681109536E-06],ID["EPSG",8609]],PARAMETER["Z-axis rotation",-0.917,ANGLEUNIT["arc-second",4.84813681109536E-06],ID["EPSG",8610]],PARAMETER["Scale difference",-3.52,SCALEUNIT["parts per million",1E-06],ID["EPSG",8611]],OPERATIONACCURACY[1.0],ID["EPSG",1612],REMARK["Parameter values are taken from ED50 to ETRS89 (1), code 1588, (rotations converted from μrad). Adopted for ED50 to WGS 84 CTs offshore Norway north of 62°N from April 2001, replacing code 1590. Included in Statens Kartverk programme wsktrans from v4.0."]]
1.2 OSDU BoundCRS to UTM 32N¶
Transform an OSDU BoundCRS from longitude, latitude, and ellipsoidal height into UTM zone 32N, then apply the inverse transformation using the generated easting and northing values. The bound CRS is stated as a persistableReference payload (the same kind used throughout section 6), so this needs no OSDU catalogue database, custom or otherwise.
from geodetic_engine.geodesy import Transformation
# OSDU record BoundGeographic2D:EPSG::4230_EPSG::1612, stated entirely by its
# persistableReference payload rather than looked up by an "OSDU:" code -- see
# section 6 for more on this. This way the example needs no catalogue
# database, custom or otherwise.
BOUND_ED50_VIA_1612 = r"""
{"authCode":{"auth":"OSDU",
"code":"4230023"},
"lateBoundCRS":{"authCode":{"auth":"EPSG","code":"4230"},
"name":"GCS_European_1950",
"type":"LBC",
"ver":"PE_10_9_1",
"wkt":"GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433],AUTHORITY[\"EPSG\",4230]]"},
"name":"ED50 * EPSG-Nor N62 2001 [4230,1612]",
"singleCT":{"authCode":{"auth":"EPSG","code":"1612"},
"name":"ED_1950_To_WGS_1984_23",
"type":"ST",
"ver":"PE_10_9_1",
"wkt":"GEOGTRAN[\"ED_1950_To_WGS_1984_23\",GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],METHOD[\"Position_Vector\"],PARAMETER[\"X_Axis_Translation\",-116.641],PARAMETER[\"Y_Axis_Translation\",-56.931],PARAMETER[\"Z_Axis_Translation\",-110.559],PARAMETER[\"X_Axis_Rotation\",0.893],PARAMETER[\"Y_Axis_Rotation\",0.921],PARAMETER[\"Z_Axis_Rotation\",-0.917],PARAMETER[\"Scale_Difference\",-3.52],OPERATIONACCURACY[1.0],AUTHORITY[\"EPSG\",1612]]"},
"type":"EBC",
"ver":"PE_10_9_1"}
"""
bound_crs = BOUND_ED50_VIA_1612 # name: ED50 * EPSG-Nor N62 2001 [4230,1612]
trg_crs = "EPSG:32632"
lon, lat, h = 10, 60, 100
ct = Transformation(source_crs=bound_crs, target_crs=trg_crs)
results = ct.transform(lon, lat, h)
trans_lon, trans_lat, trans_h = results.coordinates[0]
print(trans_lon, trans_lat, trans_h)
# See the results in JSON format and WKT
print(results.to_json())
print(results.to_json(pretty=False))
print(results.operation.to_wkt())
# Inverse direction: feed the UTM point the forward step just produced, not the
# original lon/lat -- the source CRS here is EPSG:32632, in metres.
east, north = trans_lon, trans_lat
ct = Transformation(source_crs=trg_crs, target_crs=bound_crs)
results = ct.transform(east, north, h)
trans_lon, trans_lat, trans_h = results.coordinates[0]
print(trans_lon, trans_lat, trans_h)
555699.3111382048 6651781.958444249 100.0
{
"coordinates": [
[
555699.3111382048,
6651781.958444249,
100.0
]
],
"coordinate_order": "xy",
"coordinate_epoch": null,
"source_crs": "ED50",
"target_crs": "EPSG:32632",
"source_crs_wkt": "BOUNDCRS[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]]],ABRIDGEDTRANSFORMATION[\"ED_1950_To_WGS_1984_23\",METHOD[\"Position Vector transformation (geog2D domain)\",ID[\"EPSG\",9606]],PARAMETER[\"X-axis translation\",-116.641,ID[\"EPSG\",8605]],PARAMETER[\"Y-axis translation\",-56.931,ID[\"EPSG\",8606]],PARAMETER[\"Z-axis translation\",-110.559,ID[\"EPSG\",8607]],PARAMETER[\"X-axis rotation\",0.893,ID[\"EPSG\",8608]],PARAMETER[\"Y-axis rotation\",0.921,ID[\"EPSG\",8609]],PARAMETER[\"Z-axis rotation\",-0.917,ID[\"EPSG\",8610]],PARAMETER[\"Scale difference\",0.99999648,ID[\"EPSG\",8611]]]]",
"target_crs_wkt": "PROJCRS[\"WGS 84 / UTM zone 32N\",BASEGEOGCRS[\"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]],ID[\"EPSG\",4326]],CONVERSION[\"UTM zone 32N\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",9,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Navigation and medium accuracy spatial referencing.\"],AREA[\"Between 6\u00b0E and 12\u00b0E, northern hemisphere between equator and 84\u00b0N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.\"],BBOX[0,6,84,12]],ID[\"EPSG\",32632]]",
"source_axes": [
"Lat",
"Lon"
],
"source_units": [
"degree",
"degree"
],
"target_axes": [
"E",
"N"
],
"target_units": [
"metre",
"metre"
],
"operation": {
"requested": null,
"applied": null,
"name": "ED_1950_To_WGS_1984_23",
"method": "Position Vector transformation (geog2D domain)",
"accuracy_m": null,
"route": "bound",
"steps": [
"ED_1950_To_WGS_1984_23",
"UTM zone 32N",
"axis order change (2D)"
],
"ballpark": false,
"requires_epoch": false,
"execution_direction": "FORWARD",
"bound_operations": [
"ED_1950_To_WGS_1984_23"
],
"axis_order_corrected": false,
"definition": {
"type": "Transformation",
"name": "ED_1950_To_WGS_1984_23",
"source_crs": {
"type": "GeographicCRS",
"name": "ED50",
"datum": {
"type": "GeodeticReferenceFrame",
"name": "European Datum 1950",
"ellipsoid": {
"name": "International 1924",
"semi_major_axis": 6378388,
"inverse_flattening": 297
}
},
"coordinate_system": {
"subtype": "ellipsoidal",
"axis": [
{
"name": "Geodetic latitude",
"abbreviation": "Lat",
"direction": "north",
"unit": "degree"
},
{
"name": "Geodetic longitude",
"abbreviation": "Lon",
"direction": "east",
"unit": "degree"
}
]
},
"id": {
"authority": "EPSG",
"code": 4230
}
},
"target_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"
}
]
},
"id": {
"authority": "EPSG",
"code": 4326
}
},
"method": {
"name": "Position Vector transformation (geog2D domain)",
"id": {
"authority": "EPSG",
"code": 9606
}
},
"parameters": [
{
"name": "X-axis translation",
"value": -116.641,
"unit": "metre",
"id": {
"authority": "EPSG",
"code": 8605
}
},
{
"name": "Y-axis translation",
"value": -56.931,
"unit": "metre",
"id": {
"authority": "EPSG",
"code": 8606
}
},
{
"name": "Z-axis translation",
"value": -110.559,
"unit": "metre",
"id": {
"authority": "EPSG",
"code": 8607
}
},
{
"name": "X-axis rotation",
"value": 0.893,
"unit": {
"type": "AngularUnit",
"name": "arc-second",
"conversion_factor": 4.84813681109536e-06
},
"id": {
"authority": "EPSG",
"code": 8608
}
},
{
"name": "Y-axis rotation",
"value": 0.921,
"unit": {
"type": "AngularUnit",
"name": "arc-second",
"conversion_factor": 4.84813681109536e-06
},
"id": {
"authority": "EPSG",
"code": 8609
}
},
{
"name": "Z-axis rotation",
"value": -0.917,
"unit": {
"type": "AngularUnit",
"name": "arc-second",
"conversion_factor": 4.84813681109536e-06
},
"id": {
"authority": "EPSG",
"code": 8610
}
},
{
"name": "Scale difference",
"value": -3.52,
"unit": {
"type": "ScaleUnit",
"name": "parts per million",
"conversion_factor": 1e-06
},
"id": {
"authority": "EPSG",
"code": 8611
}
}
]
}
},
"grids": [],
"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=utm zone=32 ellps=WGS84",
"database_fingerprints": {
"/usr/local/share/proj/proj.db": "3d3d8b279d100544b5737cd08a732335170fff09ee91166b4137919d6d019198"
}
}
{"coordinates": [[555699.3111382048, 6651781.958444249, 100.0]], "coordinate_order": "xy", "coordinate_epoch": null, "source_crs": "ED50", "target_crs": "EPSG:32632", "source_crs_wkt": "BOUNDCRS[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]]],ABRIDGEDTRANSFORMATION[\"ED_1950_To_WGS_1984_23\",METHOD[\"Position Vector transformation (geog2D domain)\",ID[\"EPSG\",9606]],PARAMETER[\"X-axis translation\",-116.641,ID[\"EPSG\",8605]],PARAMETER[\"Y-axis translation\",-56.931,ID[\"EPSG\",8606]],PARAMETER[\"Z-axis translation\",-110.559,ID[\"EPSG\",8607]],PARAMETER[\"X-axis rotation\",0.893,ID[\"EPSG\",8608]],PARAMETER[\"Y-axis rotation\",0.921,ID[\"EPSG\",8609]],PARAMETER[\"Z-axis rotation\",-0.917,ID[\"EPSG\",8610]],PARAMETER[\"Scale difference\",0.99999648,ID[\"EPSG\",8611]]]]", "target_crs_wkt": "PROJCRS[\"WGS 84 / UTM zone 32N\",BASEGEOGCRS[\"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]],ID[\"EPSG\",4326]],CONVERSION[\"UTM zone 32N\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",9,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Navigation and medium accuracy spatial referencing.\"],AREA[\"Between 6\u00b0E and 12\u00b0E, northern hemisphere between equator and 84\u00b0N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.\"],BBOX[0,6,84,12]],ID[\"EPSG\",32632]]", "source_axes": ["Lat", "Lon"], "source_units": ["degree", "degree"], "target_axes": ["E", "N"], "target_units": ["metre", "metre"], "operation": {"requested": null, "applied": null, "name": "ED_1950_To_WGS_1984_23", "method": "Position Vector transformation (geog2D domain)", "accuracy_m": null, "route": "bound", "steps": ["ED_1950_To_WGS_1984_23", "UTM zone 32N", "axis order change (2D)"], "ballpark": false, "requires_epoch": false, "execution_direction": "FORWARD", "bound_operations": ["ED_1950_To_WGS_1984_23"], "axis_order_corrected": false, "definition": {"type": "Transformation", "name": "ED_1950_To_WGS_1984_23", "source_crs": {"type": "GeographicCRS", "name": "ED50", "datum": {"type": "GeodeticReferenceFrame", "name": "European Datum 1950", "ellipsoid": {"name": "International 1924", "semi_major_axis": 6378388, "inverse_flattening": 297}}, "coordinate_system": {"subtype": "ellipsoidal", "axis": [{"name": "Geodetic latitude", "abbreviation": "Lat", "direction": "north", "unit": "degree"}, {"name": "Geodetic longitude", "abbreviation": "Lon", "direction": "east", "unit": "degree"}]}, "id": {"authority": "EPSG", "code": 4230}}, "target_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"}]}, "id": {"authority": "EPSG", "code": 4326}}, "method": {"name": "Position Vector transformation (geog2D domain)", "id": {"authority": "EPSG", "code": 9606}}, "parameters": [{"name": "X-axis translation", "value": -116.641, "unit": "metre", "id": {"authority": "EPSG", "code": 8605}}, {"name": "Y-axis translation", "value": -56.931, "unit": "metre", "id": {"authority": "EPSG", "code": 8606}}, {"name": "Z-axis translation", "value": -110.559, "unit": "metre", "id": {"authority": "EPSG", "code": 8607}}, {"name": "X-axis rotation", "value": 0.893, "unit": {"type": "AngularUnit", "name": "arc-second", "conversion_factor": 4.84813681109536e-06}, "id": {"authority": "EPSG", "code": 8608}}, {"name": "Y-axis rotation", "value": 0.921, "unit": {"type": "AngularUnit", "name": "arc-second", "conversion_factor": 4.84813681109536e-06}, "id": {"authority": "EPSG", "code": 8609}}, {"name": "Z-axis rotation", "value": -0.917, "unit": {"type": "AngularUnit", "name": "arc-second", "conversion_factor": 4.84813681109536e-06}, "id": {"authority": "EPSG", "code": 8610}}, {"name": "Scale difference", "value": -3.52, "unit": {"type": "ScaleUnit", "name": "parts per million", "conversion_factor": 1e-06}, "id": {"authority": "EPSG", "code": 8611}}]}}, "grids": [], "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=utm zone=32 ellps=WGS84", "database_fingerprints": {"/usr/local/share/proj/proj.db": "3d3d8b279d100544b5737cd08a732335170fff09ee91166b4137919d6d019198"}}
COORDINATEOPERATION["ED_1950_To_WGS_1984_23",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["Position Vector transformation (geog2D domain)",ID["EPSG",9606]],PARAMETER["X-axis translation",-116.641,LENGTHUNIT["metre",1],ID["EPSG",8605]],PARAMETER["Y-axis translation",-56.931,LENGTHUNIT["metre",1],ID["EPSG",8606]],PARAMETER["Z-axis translation",-110.559,LENGTHUNIT["metre",1],ID["EPSG",8607]],PARAMETER["X-axis rotation",0.893,ANGLEUNIT["arc-second",4.84813681109536E-06],ID["EPSG",8608]],PARAMETER["Y-axis rotation",0.921,ANGLEUNIT["arc-second",4.84813681109536E-06],ID["EPSG",8609]],PARAMETER["Z-axis rotation",-0.917,ANGLEUNIT["arc-second",4.84813681109536E-06],ID["EPSG",8610]],PARAMETER["Scale difference",-3.52,SCALEUNIT["parts per million",1E-06],ID["EPSG",8611]]]
10.00000000073917 59.99999999956904 100.0
1.3 Naming a list of operations¶
A transformation route can contain several named EPSG operations. The example compares a two-step route with its published concatenated operation and checks the reverse direction.
from geodetic_engine.geodesy import Transformation, available_operations
src_crs = "EPSG:4230" # ED50 Geographic 2D CRS
trg_crs = "EPSG:4326" # WGS 84 Geographic 2D CRS
points = [[4.12789451, 63.58496782, 100]]
ct = Transformation(
source_crs=src_crs, target_crs=trg_crs, operation=["EPSG:1147", "EPSG:1146"]
)
results = ct.transform(*points)
trans_lon, trans_lat, trans_h = results.coordinates[0]
print(f"Forward: {trans_lon}, {trans_lat}, {trans_h}")
# Shall produce the same result as the previous transformation
ct = Transformation(source_crs=src_crs, target_crs=trg_crs, operation=["EPSG:8047"])
results = ct.transform(*points)
trans_lon, trans_lat, trans_h = results.coordinates[0]
print(f"Forward: {trans_lon}, {trans_lat}, {trans_h}")
# # Reversed direction
ct_rev = Transformation(
source_crs=trg_crs, target_crs=src_crs, operation=["EPSG:1147", "EPSG:1146"]
)
results_rev = ct_rev.transform(*results.coordinates[0])
rev_lon, rev_lat, rev_h = results_rev.coordinates[0]
print(f"Reversed: {rev_lon:.8f}, {rev_lat:.8f}, {rev_h:.8f}")
Forward: 4.126139926951677, 63.584613419996394, 100.0
Forward: 4.126139926951677, 63.584613419996394, 100.0
Reversed: 4.12789452, 63.58496782, 100.00000000
2. Vertical transformations¶
The following examples combine horizontal and vertical components, including a height transformation to the NN54 vertical reference system.
2.1 Vertical transformation within one horizontal datum¶
Apply a named vertical operation from ellipsoidal height to the NN54 height reference while transforming between the associated CRS definitions.
from geodetic_engine.geodesy import Transformation
bound_crs = "EPSG:4937"
trg_crs = "EPSG:6172" # ETRS89 / UTM zone 32N + NN54 height
operation = "EPSG:11559"
points = [[11.12789451, 63.58496782, 100]]
cts = available_operations(source_crs=bound_crs, target_crs=trg_crs)
print(cts)
ct = Transformation(source_crs=bound_crs, target_crs=trg_crs, operation=operation)
results = ct.transform(*points)
trans_lon, trans_lat, trans_h = results.coordinates[0]
print(trans_lon, trans_lat, trans_h) # shall be [[605606.253, 7052523.904, 61.742]]
# Reversed direction
ct_rev = Transformation(source_crs=trg_crs, target_crs=bound_crs, operation=operation)
results_rev = ct_rev.transform(*results.coordinates[0])
rev_lon, rev_lat, rev_h = results_rev.coordinates[0]
print(rev_lon, rev_lat, rev_h) # shall be [[11.12789451, 63.58496782, 100]]
(OperationCandidate(auth_name=None, code=None, name='ETRS89-NOR [EUREF89] to NN54 height (1) + UTM zone 32N', method_name=None, accuracy=0.02, area_of_use=AreaOfUse(west=4.39, south=57.9, east=31.32, north=71.24, name='Norway - onshore.'), ballpark=False, requires_epoch=False, grids=(GridUsage(name='no_kv_href2008a.tif', full_name='/usr/local/share/proj/no_kv_href2008a.tif', package_name='', url='https://cdn.proj.org/no_kv_href2008a.tif', available=True, open_license=True, direct_download=True),), usable=True, steps=(OperationStep(auth_name=None, code=None, name='ETRS89-NOR [EUREF89] to NN54 height (1)', method_name='PROJ-based operation method (approximate): +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=no_kv_href2008a.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1'), OperationStep(auth_name='EPSG', code='16032', name='UTM zone 32N', method_name='Transverse Mercator'))), OperationCandidate(auth_name=None, code=None, name='Inverse of Transformation from NN54 height to ETRS89 (ballpark vertical transformation, without ellipsoid height to vertical height correction) + UTM zone 32N', method_name=None, accuracy=None, area_of_use=AreaOfUse(west=-180.0, south=-90.0, east=180.0, north=90.0, name='World'), ballpark=False, requires_epoch=False, grids=(), usable=True, steps=(OperationStep(auth_name=None, code=None, name='Inverse of Transformation from NN54 height to ETRS89 (ballpark vertical transformation, without ellipsoid height to vertical height correction)', method_name='PROJ-based operation method (approximate): +proj=noop'), OperationStep(auth_name='EPSG', code='16032', name='UTM zone 32N', method_name='Transverse Mercator'))))
605606.253274944 7052523.904230434 61.741534204810534
11.12789451 63.58496782000001 100.0
2.2 Vertical transformation with a horizontal datum change¶
This route combines a horizontal datum transformation with a vertical height transformation. Both component operations are named explicitly so that no part of the route is selected silently.
from geodetic_engine.geodesy import Transformation, available_operations
bound_crs = "EPSG:4979"
trg_crs = "EPSG:6172"
points = [(11.12789451, 63.58496782, 100)]
operations = ["EPSG:11028", "EPSG:11559"]
cts = available_operations(source_crs=bound_crs, target_crs=trg_crs)
print(cts)
ct = Transformation(source_crs=bound_crs, target_crs=trg_crs, operation=operations)
results = ct.transform(points)
trans_lon, trans_lat, trans_h = results.coordinates[0]
print(f"Forward: {trans_lon}, {trans_lat}, {trans_h}")
ct_rev = Transformation(
source_crs=trg_crs, target_crs=bound_crs, operation=list(reversed(operations))
)
results_rev = ct_rev.transform(results.coordinates)
rev_lon, rev_lat, rev_h = results_rev.coordinates[0]
print(f"Reversed: {rev_lon}, {rev_lat}, {rev_h}")
(OperationCandidate(auth_name=None, code=None, name='Inverse of ETRS89-NOR [EUREF89] to WGS 84 (1) + ETRS89-NOR [EUREF89] to NN54 height (1) + UTM zone 32N', method_name=None, accuracy=1.02, area_of_use=AreaOfUse(west=4.39, south=57.9, east=31.32, north=71.24, name='Norway - onshore.'), ballpark=False, requires_epoch=False, grids=(GridUsage(name='no_kv_href2008a.tif', full_name='/usr/local/share/proj/no_kv_href2008a.tif', package_name='', url='https://cdn.proj.org/no_kv_href2008a.tif', available=True, open_license=True, direct_download=True),), usable=True, steps=(OperationStep(auth_name=None, code=None, name='Inverse of ETRS89-NOR [EUREF89] to WGS 84 (1) + ETRS89-NOR [EUREF89] to NN54 height (1)', method_name='PROJ-based operation method (approximate): +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=no_kv_href2008a.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1'), OperationStep(auth_name='EPSG', code='16032', name='UTM zone 32N', method_name='Transverse Mercator'))), OperationCandidate(auth_name=None, code=None, name='Inverse of Transformation from NN54 height to WGS 84 (ballpark vertical transformation, without ellipsoid height to vertical height correction) + Inverse of Ballpark geographic offset from ETRS89-NOR [EUREF89] to WGS 84 + UTM zone 32N', method_name=None, accuracy=None, area_of_use=AreaOfUse(west=-180.0, south=-90.0, east=180.0, north=90.0, name='World'), ballpark=False, requires_epoch=False, grids=(), usable=True, steps=(OperationStep(auth_name=None, code=None, name='Inverse of Transformation from NN54 height to WGS 84 (ballpark vertical transformation, without ellipsoid height to vertical height correction) + Inverse of Ballpark geographic offset from ETRS89-NOR [EUREF89] to WGS 84', method_name='PROJ-based operation method (approximate): +proj=noop'), OperationStep(auth_name='EPSG', code='16032', name='UTM zone 32N', method_name='Transverse Mercator'))), OperationCandidate(auth_name=None, code=None, name='Inverse of Transformation from NN54 height to WGS 84 (ballpark vertical transformation, without ellipsoid height to vertical height correction) + Inverse of ETRS89-NOR [EUREF89] to WGS 84 (1) + UTM zone 32N', method_name=None, accuracy=None, area_of_use=AreaOfUse(west=-3.35, south=56.08, east=38.01, north=84.73, name='Norway including Svalbard - onshore and offshore.'), ballpark=False, requires_epoch=False, grids=(), usable=True, steps=(OperationStep(auth_name=None, code=None, name='Inverse of Transformation from NN54 height to WGS 84 (ballpark vertical transformation, without ellipsoid height to vertical height correction) + Inverse of ETRS89-NOR [EUREF89] to WGS 84 (1)', method_name='PROJ-based operation method (approximate): +proj=noop'), OperationStep(auth_name='EPSG', code='16032', name='UTM zone 32N', method_name='Transverse Mercator'))))
Forward: 605606.253274944, 7052523.904230434, 61.741534204810534
Reversed: 11.12789451, 63.58496782000001, 100.0
3. Engineering CRS transformations¶
Engineering CRSs use local coordinate systems and may rely on similarity or other specialized transformations. This example compares a named operation with a usable candidate discovered through available_operations() and selected explicitly.
3.1 Explicit and discovered operation selection¶
Use a registered engineering CRS pair with an explicit EPSG operation, compare it with the operation discovered through available_operations(), and inspect the transformed coordinates in both directions.
from geodetic_engine.geodesy import (
CoordinateReferenceSystem,
Transformation,
available_operations,
)
src_crs = CoordinateReferenceSystem.from_user_input("EPSG:5817") # Tombak LNG plant
trg_crs = CoordinateReferenceSystem.from_user_input(
"EPSG:3307"
) # Nakhl-e Ghanem / UTM zone 39N
transformation = "EPSG:15747"
x, y = 20000.0, 10000.0
# Explicit EPSG operation.
tombak_result = Transformation(src_crs, trg_crs, operation=transformation).transform(
x, y
)
coords = tombak_result.coordinates
print(f"Result: {coords[0][0]:.3f}, {coords[0][1]:.3f}") # 618336.748, 3067774.210
# Reversed transformation
tombak_reversed_result = Transformation(
trg_crs, src_crs, operation=transformation
).transform(coords[0][0], coords[0][1])
reversed_coords = tombak_reversed_result.coordinates
print(
f"Reversed Result: {reversed_coords[0][0]:.3f}, {reversed_coords[0][1]:.3f}"
) # 20000.000, 10000.000
# A datum change is never resolved silently -- allow_any_operation no longer
# does that. available_operations() is how a candidate is discovered
# instead; here there is only one usable candidate, so there is nothing left
# to choose between.
candidates = [c for c in available_operations(src_crs, trg_crs) if c.usable]
assert [c.authority_code for c in candidates] == [transformation]
auto_result = Transformation(src_crs, trg_crs, operation=candidates[0]).transform(x, y)
auto_coords = auto_result.coordinates
print(f"Discovered Result: {auto_coords[0][0]:.3f}, {auto_coords[0][1]:.3f}")
print(
f"Using the: {tombak_result.operation.name} "
f"({tombak_result.operation.auth_name}:{tombak_result.operation.code})"
)
Result: 618336.748, 3067774.210
Reversed Result: 20000.000, 10000.000
Discovered Result: 618336.748, 3067774.210
Using the: Tombak LNG Plant Grid to Nakhl-e Ghanem / UTM zone 39N (1) (EPSG:15747)
4. Operation selection and coordinate export¶
These examples show how to discover and explicitly select an applicable operation, then export the resulting coordinates.
4.1 Discover and explicitly select an operation¶
A datum change is never resolved silently. available_operations() discovers candidates; the example filters them by usability and area of use, then explicitly requests the most accurate applicable one. allow_any_operation does not enable automatic selection. The result records the operation actually used.
from geodetic_engine.geodesy import Transformation, available_operations
src_crs = "EPSG:4230" # ED50 Geographic 2D CRS
trg_crs = "EPSG:32632" # WGS 84 / UTM zone 32N
points = [[11.12789451, 63.58496782, 100]]
lon, lat, _ = points[0]
def _covers(area, lon, lat):
"""Whether a candidate's declared area of use contains the point."""
if area is None:
return False
west, south, east, north = area.bounds
if east < west: # crosses the antimeridian
return (lon >= west or lon <= east) and south <= lat <= north
return west <= lon <= east and south <= lat <= north
candidates = sorted(
(
c
for c in available_operations(src_crs, trg_crs)
if c.usable and _covers(c.area_of_use, lon, lat)
),
key=lambda c: c.accuracy,
)
best = candidates[0]
print(f"Selected: {best.name} (accuracy {best.accuracy} m, area: {best.area_of_use})")
ct = Transformation(source_crs=src_crs, target_crs=trg_crs, operation=best)
results = ct.transform(*points)
trans_lon, trans_lat, trans_h = results.coordinates[0]
print(f"Forward: {trans_lon}, {trans_lat}, {trans_h}")
# See the results in JSON format and WKT
print(results.to_json(pretty=True))
Selected: ED50 to WGS 84 (23) + UTM zone 32N (accuracy 1.0 m, area: Norway - offshore north of 62°N. Also Svalbard - onshore and offshore.)
Forward: 605532.3193044778, 7052489.142253298, 100.0
{
"coordinates": [
[
605532.3193044778,
7052489.142253298,
100.0
]
],
"coordinate_order": "xy",
"coordinate_epoch": null,
"source_crs": "EPSG:4230",
"target_crs": "EPSG:32632",
"source_crs_wkt": "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]],USAGE[SCOPE[\"Geodesy.\"],AREA[\"Europe - west: Andorra; Cyprus; Denmark - onshore and offshore; Faroe Islands - onshore; France - offshore; Germany - offshore North Sea; Gibraltar; Greece - offshore; Israel - offshore; Italy including San Marino and Vatican City State; Ireland offshore; Malta; Netherlands - offshore; North Sea; Norway including Svalbard - onshore and offshore; Portugal - mainland - offshore; Spain - onshore; T\u00fcrkiye (Turkey) - onshore and offshore; United Kingdom - UKCS offshore east of 6\u00b0W including Channel Islands (Guernsey and Jersey). Egypt - Western Desert; Iraq - onshore; Jordan.\"],BBOX[25.71,-16.1,84.73,48.61]],ID[\"EPSG\",4230]]",
"target_crs_wkt": "PROJCRS[\"WGS 84 / UTM zone 32N\",BASEGEOGCRS[\"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]],ID[\"EPSG\",4326]],CONVERSION[\"UTM zone 32N\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",9,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Navigation and medium accuracy spatial referencing.\"],AREA[\"Between 6\u00b0E and 12\u00b0E, northern hemisphere between equator and 84\u00b0N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.\"],BBOX[0,6,84,12]],ID[\"EPSG\",32632]]",
"source_axes": [
"Lat",
"Lon"
],
"source_units": [
"degree",
"degree"
],
"target_axes": [
"E",
"N"
],
"target_units": [
"metre",
"metre"
],
"operation": {
"requested": "EPSG:1612 + EPSG:16032",
"applied": null,
"name": "ED50 to WGS 84 (23) + UTM zone 32N",
"method": null,
"accuracy_m": 1.0,
"route": "transformer_group",
"steps": [
"ED50 to WGS 84 (23)",
"UTM zone 32N",
"axis order change (2D)"
],
"ballpark": false,
"requires_epoch": false,
"execution_direction": "FORWARD",
"bound_operations": [],
"axis_order_corrected": false,
"definition": {
"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json",
"type": "ConcatenatedOperation",
"name": "ED50 to WGS 84 (23) + UTM zone 32N (with axis order normalized for visualization)",
"source_crs": {
"type": "GeographicCRS",
"name": "ED50 (with axis order normalized for visualization)",
"datum": {
"type": "GeodeticReferenceFrame",
"name": "European Datum 1950",
"ellipsoid": {
"name": "International 1924",
"semi_major_axis": 6378388,
"inverse_flattening": 297
},
"id": {
"authority": "EPSG",
"code": 6230
}
},
"coordinate_system": {
"subtype": "ellipsoidal",
"axis": [
{
"name": "Geodetic longitude",
"abbreviation": "Lon",
"direction": "east",
"unit": "degree"
},
{
"name": "Geodetic latitude",
"abbreviation": "Lat",
"direction": "north",
"unit": "degree"
}
]
},
"remarks": "Axis order reversed compared to EPSG:4230"
},
"target_crs": {
"type": "ProjectedCRS",
"name": "WGS 84 / UTM zone 32N",
"base_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"
}
]
},
"id": {
"authority": "EPSG",
"code": 4326
}
},
"conversion": {
"name": "UTM zone 32N",
"method": {
"name": "Transverse Mercator",
"id": {
"authority": "EPSG",
"code": 9807
}
},
"parameters": [
{
"name": "Latitude of natural origin",
"value": 0,
"unit": "degree",
"id": {
"authority": "EPSG",
"code": 8801
}
},
{
"name": "Longitude of natural origin",
"value": 9,
"unit": "degree",
"id": {
"authority": "EPSG",
"code": 8802
}
},
{
"name": "Scale factor at natural origin",
"value": 0.9996,
"unit": "unity",
"id": {
"authority": "EPSG",
"code": 8805
}
},
{
"name": "False easting",
"value": 500000,
"unit": "metre",
"id": {
"authority": "EPSG",
"code": 8806
}
},
{
"name": "False northing",
"value": 0,
"unit": "metre",
"id": {
"authority": "EPSG",
"code": 8807
}
}
]
},
"coordinate_system": {
"subtype": "Cartesian",
"axis": [
{
"name": "Easting",
"abbreviation": "E",
"direction": "east",
"unit": "metre"
},
{
"name": "Northing",
"abbreviation": "N",
"direction": "north",
"unit": "metre"
}
]
},
"id": {
"authority": "EPSG",
"code": 32632
}
},
"steps": [
{
"type": "Conversion",
"name": "axis order change (2D)",
"method": {
"name": "Axis Order Reversal (2D)",
"id": {
"authority": "EPSG",
"code": 9843
}
},
"id": {
"authority": "EPSG",
"code": 15498
}
},
{
"type": "Transformation",
"name": "ED50 to WGS 84 (23)",
"source_crs": {
"type": "GeographicCRS",
"name": "ED50",
"datum": {
"type": "GeodeticReferenceFrame",
"name": "European Datum 1950",
"ellipsoid": {
"name": "International 1924",
"semi_major_axis": 6378388,
"inverse_flattening": 297
}
},
"coordinate_system": {
"subtype": "ellipsoidal",
"axis": [
{
"name": "Geodetic latitude",
"abbreviation": "Lat",
"direction": "north",
"unit": "degree"
},
{
"name": "Geodetic longitude",
"abbreviation": "Lon",
"direction": "east",
"unit": "degree"
}
]
},
"id": {
"authority": "EPSG",
"code": 4230
}
},
"target_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"
}
]
},
"id": {
"authority": "EPSG",
"code": 4326
}
},
"method": {
"name": "Position Vector transformation (geog2D domain)",
"id": {
"authority": "EPSG",
"code": 9606
}
},
"parameters": [
{
"name": "X-axis translation",
"value": -116.641,
"unit": "metre",
"id": {
"authority": "EPSG",
"code": 8605
}
},
{
"name": "Y-axis translation",
"value": -56.931,
"unit": "metre",
"id": {
"authority": "EPSG",
"code": 8606
}
},
{
"name": "Z-axis translation",
"value": -110.559,
"unit": "metre",
"id": {
"authority": "EPSG",
"code": 8607
}
},
{
"name": "X-axis rotation",
"value": 0.893,
"unit": {
"type": "AngularUnit",
"name": "arc-second",
"conversion_factor": 4.84813681109536e-06
},
"id": {
"authority": "EPSG",
"code": 8608
}
},
{
"name": "Y-axis rotation",
"value": 0.921,
"unit": {
"type": "AngularUnit",
"name": "arc-second",
"conversion_factor": 4.84813681109536e-06
},
"id": {
"authority": "EPSG",
"code": 8609
}
},
{
"name": "Z-axis rotation",
"value": -0.917,
"unit": {
"type": "AngularUnit",
"name": "arc-second",
"conversion_factor": 4.84813681109536e-06
},
"id": {
"authority": "EPSG",
"code": 8610
}
},
{
"name": "Scale difference",
"value": -3.52,
"unit": {
"type": "ScaleUnit",
"name": "parts per million",
"conversion_factor": 1e-06
},
"id": {
"authority": "EPSG",
"code": 8611
}
}
],
"accuracy": "1.0",
"id": {
"authority": "EPSG",
"code": 1612
},
"remarks": "Parameter values are taken from ED50 to ETRS89 (1), code 1588, (rotations converted from \u03bcrad). Adopted for ED50 to WGS 84 CTs offshore Norway north of 62\u00b0N from April 2001, replacing code 1590. Included in Statens Kartverk programme wsktrans from v4.0."
},
{
"type": "Conversion",
"name": "UTM zone 32N",
"method": {
"name": "Transverse Mercator",
"id": {
"authority": "EPSG",
"code": 9807
}
},
"parameters": [
{
"name": "Latitude of natural origin",
"value": 0,
"unit": "degree",
"id": {
"authority": "EPSG",
"code": 8801
}
},
{
"name": "Longitude of natural origin",
"value": 9,
"unit": "degree",
"id": {
"authority": "EPSG",
"code": 8802
}
},
{
"name": "Scale factor at natural origin",
"value": 0.9996,
"unit": "unity",
"id": {
"authority": "EPSG",
"code": 8805
}
},
{
"name": "False easting",
"value": 500000,
"unit": "metre",
"id": {
"authority": "EPSG",
"code": 8806
}
},
{
"name": "False northing",
"value": 0,
"unit": "metre",
"id": {
"authority": "EPSG",
"code": 8807
}
}
],
"id": {
"authority": "EPSG",
"code": 16032
}
}
],
"accuracy": "1",
"area": "Norway - offshore north of 62\u00b0N. Also Svalbard - onshore and offshore.",
"bbox": {
"south_latitude": 62,
"west_longitude": -3.35,
"north_latitude": 84.73,
"east_longitude": 38.01
}
}
},
"grids": [],
"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=utm zone=32 ellps=WGS84",
"database_fingerprints": {
"/usr/local/share/proj/proj.db": "3d3d8b279d100544b5737cd08a732335170fff09ee91166b4137919d6d019198"
}
}
4.2 Export coordinates¶
TransformationResult.coordinates can be exported as plain Python lists, NumPy arrays, or a pandas DataFrame. DataFrame columns follow the target CRS axis abbreviations.
from geodetic_engine.geodesy import Transformation
ct = Transformation(
source_crs="EPSG:4258", target_crs="EPSG:25832", operation="EPSG:16032"
)
result = ct.transform([(10.75, 59.91), (5.32, 60.39)])
# Plain Python
print(result.coordinates.to_list())
# NumPy array, shape (n_points, n_axes)
print(result.coordinates.to_numpy())
# pandas DataFrame, columns named after the target CRS's declared axes
result.coordinates.to_dataframe()
[[597868.3810645777, 6642681.510038425], [297230.2202071025, 6700510.175130839]]
[[ 597868.38106458 6642681.51003843]
[ 297230.2202071 6700510.17513084]]
| E | N | |
|---|---|---|
| 0 | 597868.381065 | 6.642682e+06 |
| 1 | 297230.220207 | 6.700510e+06 |
5. Available operations¶
Use available_operations() to inspect the candidate operations for a CRS pair, including authority code, stated accuracy, usability, and area of use – the same tool used in section 4.1 to discover an operation for a datum change.
from geodetic_engine.geodesy import Transformation, available_operations
src_crs = "EPSG:4230"
trg_crs = "EPSG:4326"
accuracy = None # 10
authority = "any" # "EPSG"
# Without an explicit operation, this datum change raises AmbiguousOperationError.
ct_list = available_operations(
src_crs, trg_crs, authority=authority, accuracy=accuracy, allow_ballpark=False
)
print(f"Number of available operations: {len(ct_list)}\n")
for c in ct_list:
print(
f"{c.authority_code:10} {c.name:30} accuracy={c.accuracy!s:>5} "
f"usable={c.usable} area={c.area_of_use}"
)
# candidates = available_operations(src_crs, trg_crs)
# ct = Transformation(src_crs, trg_crs, operation=candidates[0].authority_code)
# ct
Number of available operations: 42
EPSG:1133 ED50 to WGS 84 (1) accuracy= 10.0 usable=True area=Austria; Belgium; Denmark; Finland; Faroe islands; France; Germany (west); Gibraltar; Greece; Italy; Luxembourg; Netherlands; Norway; Portugal; Spain; Sweden; Switzerland.
ESRI:108335 ED_1950_To_WGS_1984_NGA_7PAR accuracy= 10.0 usable=True area=Austria; Belgium; Denmark; Finland; Faroe islands; France; Germany (west); Gibraltar; Greece; Italy; Luxembourg; Netherlands; Norway; Portugal; Spain; Sweden; Switzerland.
EPSG:1612 ED50 to WGS 84 (23) accuracy= 1.0 usable=True area=Norway - offshore north of 62°N. Also Svalbard - onshore and offshore.
EPSG:1311 ED50 to WGS 84 (18) accuracy= 1.0 usable=True area=Denmark - offshore North Sea; Ireland - offshore; Netherlands - offshore; United Kingdom - UKCS offshore.
EPSG:1134 ED50 to WGS 84 (2) accuracy= 6.0 usable=True area=Austria; Denmark; France; Germany (west); Netherlands; Switzerland.
EPSG:8569 ED50 to WGS 84 (21) accuracy= 1.5 usable=True area=Norway - offshore north of 65°N. Also Svalbard.
EPSG:8047 ED50 to WGS 84 (15) accuracy= 1.5 usable=True area=Norway - offshore north of 65°N. Also Svalbard.
EPSG:1139 ED50 to WGS 84 (7) accuracy= 7.0 usable=True area=Finland and Norway - onshore.
EPSG:1275 ED50 to WGS 84 (17) accuracy= 2.0 usable=True area=France - onshore and offshore, mainland and Corsica (France métropolitaine including Corsica).
EPSG:1784 ED50 to WGS 84 (30) accuracy= 2.0 usable=True area=Türkiye (Turkey) - onshore and offshore.
EPSG:1075 ED50 to WGS 84 (38) accuracy= 10.0 usable=True area=Israel and Palestine - offshore Mediterranean Sea; Türkiye (Turkey) - offshore Black Sea.
EPSG:1145 ED50 to WGS 84 (13) accuracy= 9.0 usable=True area=Portugal; Spain - mainland.
EPSG:15933 ED50 to WGS 84 (41) accuracy= 1.0 usable=True area=Spain - mainland, Balearic Islands and Ceuta - onshore.
EPSG:1138 ED50 to WGS 84 (6) accuracy= 6.0 usable=True area=Channel islands - onshore. Ireland - onshore. Isle of Man - onshore. United Kingdom (UK) - onshore - England; Scotland; Wales; Northern Ireland.
ESRI:108068 ED_1950_To_WGS_1984_NTv2_PENR2009 accuracy= 0.9 usable=True area=Spain - mainland and Ceuta - onshore.
EPSG:1633 ED50 to WGS 84 (28) accuracy= 1.5 usable=True area=Spain - onshore mainland except northwest (north of 41°30'N and west of 4°30'W).
EPSG:1140 ED50 to WGS 84 (8) accuracy= 44.0 usable=True area=Greece - onshore.
EPSG:1245 ED50 to WGS 84 (16) accuracy= 44.0 usable=True area=Tunisia - onshore and offshore.
EPSG:15964 ED50 to WGS 84 (42) accuracy= 5.0 usable=True area=Portugal - mainland - offshore.
EPSG:1613 ED50 to WGS 84 (24) accuracy= 1.0 usable=True area=Norway - offshore south of 62°N - North Sea.
EPSG:1137 ED50 to WGS 84 (5) accuracy= 13.0 usable=True area=Egypt - Western Desert.
EPSG:1810 ED50 to WGS 84 (31) accuracy= 15.0 usable=True area=Egypt - Western Desert.
EPSG:1087 ED50 to WGS 84 (37) accuracy= 2.5 usable=True area=Jordan.
EPSG:1136 ED50 to WGS 84 (4) accuracy= 26.0 usable=True area=Cyprus - onshore and offshore.
EPSG:1989 ED50 to WGS 84 (34) accuracy= 1.0 usable=True area=Portugal - mainland - onshore.
EPSG:1985 ED50 to WGS 84 (33) accuracy= 5.0 usable=True area=Portugal - mainland - onshore.
EPSG:1627 ED50 to WGS 84 (25) accuracy= 1.0 usable=True area=Denmark - onshore.
EPSG:3904 ED50 to WGS 84 (32) accuracy= 5.0 usable=True area=Netherlands - offshore North Sea.
EPSG:1635 ED50 to WGS 84 (29) accuracy= 1.5 usable=True area=Spain - onshore mainland north of 41°30'N and west of 4°30' W.
EPSG:1998 ED50 to WGS 84 (36) accuracy= 1.0 usable=True area=Germany - offshore North Sea.
ESRI:108359 ED_1950_To_WGS_1984_NTv2_Catalonia accuracy= 0.05 usable=True area=Spain - Catalonia onshore.
EPSG:1143 ED50 to WGS 84 (11) accuracy= 35.0 usable=True area=Italy - Sicily onshore.
ESRI:108069 ED_1950_To_WGS_1984_NTv2_BALR2009 accuracy= 0.9 usable=True area=Spain - Balearic Islands.
EPSG:1631 ED50 to WGS 84 (27) accuracy= 1.5 usable=True area=Spain - Balearic Islands.
EPSG:1853 ED50 to WGS 84 (39) accuracy= 5.0 usable=True area=Ireland - offshore - Corrib and Errigal fields.
EPSG:1142 ED50 to WGS 84 (10) accuracy= 44.0 usable=True area=Italy - Sardinia onshore.
EPSG:1144 ED50 to WGS 84 (12) accuracy= 44.0 usable=True area=Malta - onshore.
EPSG:1629 ED50 to WGS 84 (26) accuracy= 1.0 usable=True area=Gibraltar - onshore and offshore.
EPSG:1135 ED50 to WGS 84 (3) accuracy= None usable=True area=Iraq, Israel, Jordan, Lebanon, Kuwait, Saudi Arabia and Syria - onshore.
EPSG:1440 ED50 to WGS 84 (19) accuracy= None usable=True area=Greece - onshore.
ESRI:108121 ED_1950_To_WGS_1984_PT7 accuracy= None usable=True area=Portugal - mainland - onshore.
ESRI:108111 ED_1950_To_WGS_1984_PT3 accuracy= None usable=True area=Portugal - mainland - onshore.
6. OSDU persistableReference¶
OSDU identifies the CRS, transformation or unit of a record with a persistableReference: a JSON envelope wrapping ESRI WKT, often URL-encoded and often embedded as a string inside another document.
It is self-contained. It carries the full definition rather than a code to look up, and geodetic_engine.persistablereference reads it that way — the parameters are what is built from, and the authority code beside them is provenance. A payload whose code is wrong, or names a register this machine has never seen, still resolves to exactly the CRS its sender meant.
Every cell below is standalone: each one carries the payloads it uses, taken verbatim from an OSDU reference-data catalogue, and can be run on its own in any order.
6.1 ED50 / UTM zone 32N to WGS 84, stated entirely by reference¶
Both ends of this transformation are persistableReference payloads. Nothing is named by code.
The source is an early bound (EBC) payload: ED50 / UTM zone 32N packaged with EPSG:1133, the three-parameter shift to WGS 84, so the datum change is part of the CRS definition rather than a choice made later. The target is a late bound (LBC) payload for WGS 84 geographic.
The same CRS named by its plain code raises AmbiguousOperationError, because EPSG publishes several shifts from ED50 and this package will not pick one on your behalf. Binding is what settles it.
from geodetic_engine.geodesy import transform
from geodetic_engine.persistablereference import parse_persistable_reference
# OSDU record BoundProjected:EPSG::23032_EPSG::1133
ED50_UTM32N_VIA_1133 = r"""
{"authCode":{"auth":"OSDU",
"code":"23032001"},
"lateBoundCRS":{"authCode":{"auth":"EPSG","code":"23032"},
"name":"ED_1950_UTM_Zone_32N",
"type":"LBC",
"ver":"PE_10_9_1",
"wkt":"PROJCS[\"ED_1950_UTM_Zone_32N\",GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"False_Easting\",500000.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",9.0],PARAMETER[\"Scale_Factor\",0.9996],PARAMETER[\"Latitude_Of_Origin\",0.0],UNIT[\"Meter\",1.0],AUTHORITY[\"EPSG\",23032]]"},
"name":"ED50 * DMA-mean / UTM zone 32N [23032,1133]",
"singleCT":{"authCode":{"auth":"EPSG","code":"1133"},
"name":"ED_1950_To_WGS_1984_1",
"type":"ST",
"ver":"PE_10_9_1",
"wkt":"GEOGTRAN[\"ED_1950_To_WGS_1984_1\",GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],METHOD[\"Geocentric_Translation\"],PARAMETER[\"X_Axis_Translation\",-87.0],PARAMETER[\"Y_Axis_Translation\",-98.0],PARAMETER[\"Z_Axis_Translation\",-121.0],OPERATIONACCURACY[10.0],AUTHORITY[\"EPSG\",1133]]"},
"type":"EBC",
"ver":"PE_10_9_1"}
"""
# OSDU record Geographic2D:EPSG::4326
WGS84_GEOGRAPHIC = r"""
{"authCode":{"auth":"EPSG",
"code":"4326"},
"name":"GCS_WGS_1984",
"type":"LBC",
"ver":"PE_10_9_1",
"wkt":"GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433],AUTHORITY[\"EPSG\",4326]]"}
"""
easting, northing = 500000.0, 6600000.0
source = parse_persistable_reference(ED50_UTM32N_VIA_1133)
target = parse_persistable_reference(WGS84_GEOGRAPHIC)
print(f"source {source.kind.name} {source.name}")
print(f" base {source.late_bound.authority_code} {source.late_bound.name}")
print(
f" bound by {source.operation.authority_code} {source.operation.method_names[0]}"
)
print(f"target {target.kind.name} {target.name}")
result = transform(ED50_UTM32N_VIA_1133, WGS84_GEOGRAPHIC, [(easting, northing)])
longitude, latitude = result.coordinates[0]
print(f"\n{easting:.1f} E, {northing:.1f} N -> {longitude:.9f}, {latitude:.9f}")
print(f"operation applied {result.operation.name}")
print(
f"target axes {result.target_axes}, "
f"values in {result.coordinate_order} order"
)
source EARLY_BOUND_CRS ED50 * DMA-mean / UTM zone 32N [23032,1133]
base EPSG:23032 ED_1950_UTM_Zone_32N
bound by EPSG:1133 Geocentric_Translation
target LATE_BOUND_CRS GCS_WGS_1984
500000.0 E, 6600000.0 N -> 8.998529777, 59.536487235
operation applied ED_1950_To_WGS_1984_1
target axes ('Lat', 'Lon'), values in xy order
6.2 The same transformation from a geographic CRS, into a projected CRS¶
The second example applies the same EPSG:1133 transformation, but the bound CRS wraps ED50 as a geographic CRS and the target is WGS 84 / UTM zone 32N. Again both ends are payloads.
Reading it back the other way round returns the point it started from, which is the check worth running on any early bound definition: the binding has to be applied in both directions, not just outwards to the hub.
from geodetic_engine.geodesy import transform
from geodetic_engine.persistablereference import parse_persistable_reference
# OSDU record BoundGeographic2D:EPSG::4230_EPSG::1133
ED50_GEOGRAPHIC_VIA_1133 = r"""
{"authCode":{"auth":"OSDU",
"code":"4230001"},
"lateBoundCRS":{"authCode":{"auth":"EPSG","code":"4230"},
"name":"GCS_European_1950",
"type":"LBC",
"ver":"PE_10_9_1",
"wkt":"GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433],AUTHORITY[\"EPSG\",4230]]"},
"name":"ED50 * DMA-mean [4230,1133]",
"singleCT":{"authCode":{"auth":"EPSG","code":"1133"},
"name":"ED_1950_To_WGS_1984_1",
"type":"ST",
"ver":"PE_10_9_1",
"wkt":"GEOGTRAN[\"ED_1950_To_WGS_1984_1\",GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],METHOD[\"Geocentric_Translation\"],PARAMETER[\"X_Axis_Translation\",-87.0],PARAMETER[\"Y_Axis_Translation\",-98.0],PARAMETER[\"Z_Axis_Translation\",-121.0],OPERATIONACCURACY[10.0],AUTHORITY[\"EPSG\",1133]]"},
"type":"EBC",
"ver":"PE_10_9_1"}
"""
# OSDU record Projected:EPSG::32632
WGS84_UTM32N = r"""
{"authCode":{"auth":"EPSG",
"code":"32632"},
"name":"WGS_1984_UTM_Zone_32N",
"type":"LBC",
"ver":"PE_10_9_1",
"wkt":"PROJCS[\"WGS_1984_UTM_Zone_32N\",GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"False_Easting\",500000.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",9.0],PARAMETER[\"Scale_Factor\",0.9996],PARAMETER[\"Latitude_Of_Origin\",0.0],UNIT[\"Meter\",1.0],AUTHORITY[\"EPSG\",32632]]"}
"""
source = parse_persistable_reference(ED50_GEOGRAPHIC_VIA_1133)
print(f"source {source.kind.name} {source.name} (base is geographic)")
print(
f"bound by {source.operation.authority_code} {source.operation.method_names[0]}"
)
longitude, latitude = 9.0, 59.5
forward = transform(ED50_GEOGRAPHIC_VIA_1133, WGS84_UTM32N, [(longitude, latitude)])
easting, northing = forward.coordinates[0]
print(f"\n{longitude}, {latitude} (ED50) -> {easting:.4f} E, {northing:.4f} N")
print(f"operation applied {forward.operation.name}")
back = transform(WGS84_UTM32N, ED50_GEOGRAPHIC_VIA_1133, [(easting, northing)])
print(f"and back {back.coordinates[0][0]:.9f}, {back.coordinates[0][1]:.9f}")
source EARLY_BOUND_CRS ED50 * DMA-mean [4230,1133] (base is geographic)
bound by EPSG:1133 Geocentric_Translation
9.0, 59.5 (ED50) -> 499916.8500 E, 6595675.2754 N
operation applied ED_1950_To_WGS_1984_1
and back 9.000000006, 59.500000002
6.4 Projected to projected, across a datum change¶
Both ends are projected, both arrive as payloads, and the datum change between them is the seven-parameter shift EPSG:1612 that the source payload carries.
Two projections and a datum shift are three steps, and the result reports the one that mattered. The points move about 220 m - that is the difference between ED50 and WGS 84 in the North Sea, not a rounding artefact, and it is why the shift has to be stated rather than assumed.
from geodetic_engine.geodesy import transform
# OSDU record BoundProjected:EPSG::23032_EPSG::1612 -- a seven-parameter shift
ED50_UTM32N_VIA_1612 = r"""
{"authCode":{"auth":"OSDU",
"code":"23032023"},
"lateBoundCRS":{"authCode":{"auth":"EPSG","code":"23032"},
"name":"ED_1950_UTM_Zone_32N",
"type":"LBC",
"ver":"PE_10_9_1",
"wkt":"PROJCS[\"ED_1950_UTM_Zone_32N\",GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"False_Easting\",500000.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",9.0],PARAMETER[\"Scale_Factor\",0.9996],PARAMETER[\"Latitude_Of_Origin\",0.0],UNIT[\"Meter\",1.0],AUTHORITY[\"EPSG\",23032]]"},
"name":"ED50 * EPSG-Nor N62 2001 / UTM zone 32N [23032,1612]",
"singleCT":{"authCode":{"auth":"EPSG","code":"1612"},
"name":"ED_1950_To_WGS_1984_23",
"type":"ST",
"ver":"PE_10_9_1",
"wkt":"GEOGTRAN[\"ED_1950_To_WGS_1984_23\",GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],METHOD[\"Position_Vector\"],PARAMETER[\"X_Axis_Translation\",-116.641],PARAMETER[\"Y_Axis_Translation\",-56.931],PARAMETER[\"Z_Axis_Translation\",-110.559],PARAMETER[\"X_Axis_Rotation\",0.893],PARAMETER[\"Y_Axis_Rotation\",0.921],PARAMETER[\"Z_Axis_Rotation\",-0.917],PARAMETER[\"Scale_Difference\",-3.52],OPERATIONACCURACY[1.0],AUTHORITY[\"EPSG\",1612]]"},
"type":"EBC",
"ver":"PE_10_9_1"}
"""
# OSDU record Projected:EPSG::32632
WGS84_UTM32N = r"""
{"authCode":{"auth":"EPSG",
"code":"32632"},
"name":"WGS_1984_UTM_Zone_32N",
"type":"LBC",
"ver":"PE_10_9_1",
"wkt":"PROJCS[\"WGS_1984_UTM_Zone_32N\",GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"False_Easting\",500000.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",9.0],PARAMETER[\"Scale_Factor\",0.9996],PARAMETER[\"Latitude_Of_Origin\",0.0],UNIT[\"Meter\",1.0],AUTHORITY[\"EPSG\",32632]]"}
"""
survey = [(500000.0, 6600000.0), (612345.6, 6750000.0)]
result = transform(ED50_UTM32N_VIA_1612, WGS84_UTM32N, survey)
print(f"operation applied {result.operation.name}\n")
print(f"{'ED50 / UTM 32N':<28}{'WGS 84 / UTM 32N':<28}moved")
for (east, north), (new_east, new_north) in zip(
survey, result.coordinates, strict=True
):
moved = ((new_east - east) ** 2 + (new_north - north) ** 2) ** 0.5
print(
f"{east:>11.2f} {north:>13.2f} "
f"{new_east:>11.2f} {new_north:>13.2f} {moved:6.2f} m"
)
operation applied ED_1950_To_WGS_1984_23
ED50 / UTM 32N WGS 84 / UTM 32N moved
500000.00 6600000.00 499920.40 6599793.97 220.87 m
612345.60 6750000.00 612265.78 6749793.22 221.65 m
6.6 A concatenated transformation¶
Some datum shifts are published as a chain through an intermediate frame. EPSG:8047 takes ED50 to WGS 84 by way of ED87, and OSDU states it as a compoundCT holding both steps.
A bound CRS carries one transformation, so the chain has to become a single equivalent step before it can be bound. The two Helmerts compose exactly, and the composition is checked against PROJ’s own rendering of the original chain over the area of use before it is accepted — a collapse that does not reproduce the chain is refused rather than returned. The result says collapsed to a single step so that this is visible rather than implied.
from geodetic_engine.geodesy import transform
from geodetic_engine.persistablereference import parse_persistable_reference
# OSDU record BoundGeographic2D:EPSG::4230_EPSG::8047 -- ED50 to WGS 84 via ED87
ED50_VIA_8047 = r"""
{"authCode":{"auth":"OSDU",
"code":"4230015"},
"compoundCT":{"authCode":{"auth":"EPSG","code":"8047"},
"cts":[{"authCode":{"auth":"EPSG","code":"1147"},"name":"ED_1950_To_ED_1987_2","type":"ST","ver":"PE_10_9_1","wkt":"GEOGTRAN[\"ED_1950_To_ED_1987_2\",GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],GEOGCS[\"GCS_European_1987\",DATUM[\"D_European_1987\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],METHOD[\"Position_Vector\"],PARAMETER[\"X_Axis_Translation\",-1.51],PARAMETER[\"Y_Axis_Translation\",-0.84],PARAMETER[\"Z_Axis_Translation\",-3.5],PARAMETER[\"X_Axis_Rotation\",-0.3904592782257534],PARAMETER[\"Y_Axis_Rotation\",-0.1417039218917552],PARAMETER[\"Z_Axis_Rotation\",-0.5701159244669742],PARAMETER[\"Scale_Difference\",0.609],OPERATIONACCURACY[1.0],AUTHORITY[\"EPSG\",1147]]"},{"authCode":{"auth":"EPSG","code":"1146"},"name":"ED_1987_To_WGS_1984_1","type":"ST","ver":"PE_10_9_1","wkt":"GEOGTRAN[\"ED_1987_To_WGS_1984_1\",GEOGCS[\"GCS_European_1987\",DATUM[\"D_European_1987\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],METHOD[\"Position_Vector\"],PARAMETER[\"X_Axis_Translation\",-82.981],PARAMETER[\"Y_Axis_Translation\",-99.719],PARAMETER[\"Z_Axis_Translation\",-110.709],PARAMETER[\"X_Axis_Rotation\",-0.1047000156510261],PARAMETER[\"Y_Axis_Rotation\",0.03100160037893858],PARAMETER[\"Z_Axis_Rotation\",0.08040202147511816],PARAMETER[\"Scale_Difference\",-0.3143],OPERATIONACCURACY[0.8],AUTHORITY[\"EPSG\",1146]]"}],
"name":"ED50 to WGS 84 (15)",
"policy":"Concatenated",
"type":"CT",
"ver":"PE_10_9_1"},
"lateBoundCRS":{"authCode":{"auth":"EPSG","code":"4230"},
"name":"GCS_European_1950",
"type":"LBC",
"ver":"PE_10_9_1",
"wkt":"GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433],AUTHORITY[\"EPSG\",4230]]"},
"name":"ED50 * NMA-Nor N65 1991 [4230,8047]",
"type":"EBC",
"ver":"PE_10_9_1"}
"""
# OSDU record Geographic2D:EPSG::4326
WGS84_GEOGRAPHIC = r"""
{"authCode":{"auth":"EPSG",
"code":"4326"},
"name":"GCS_WGS_1984",
"type":"LBC",
"ver":"PE_10_9_1",
"wkt":"GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433],AUTHORITY[\"EPSG\",4326]]"}
"""
chain = parse_persistable_reference(ED50_VIA_8047).operation
print(f"stated as {' then '.join(chain.method_names)}")
for step in chain.to_operation().operations:
print(f" {step.name:<24} accuracy {step.accuracy} m")
point = (4.12789451, 63.58496782)
result = transform(ED50_VIA_8047, WGS84_GEOGRAPHIC, [point])
print(f"\napplied as {result.operation.name}")
print(f"ED50 {point}")
print(f" -> WGS 84 {result.coordinates[0]}")
# The same chain named by its published code, which PROJ applies step by step.
by_code = transform("EPSG:4230", "EPSG:4326", [point], operation="EPSG:8047")
print(f"\nthe chain by code {by_code.coordinates[0]}")
print(
f"collapsed matches to {
max(
abs(first - second)
for first, second in zip(
result.coordinates[0], by_code.coordinates[0], strict=True
)
):.2e} degrees"
)
stated as Position_Vector then Position_Vector
ED_1950_To_ED_1987_2 accuracy 1.0 m
ED_1987_To_WGS_1984_1 accuracy 0.8 m
applied as ED50 to WGS 84 (15) (collapsed to a single step)
ED50 (4.12789451, 63.58496782)
-> WGS 84 (4.126139897300084, 63.584613414179096)
the chain by code (4.126139897255748, 63.584613414122096)
collapsed matches to 5.70e-11 degrees
6.7 Binding a CRS to a transformation manually¶
An OSDU catalogue publishes a CRS and a coordinate transformation as records of their own. A data record that uses them names both, and CrsReference.to_bound_crs() puts the two together, giving the same thing an early bound payload states in one piece.
This is where the choice actually lives. EPSG publishes several shifts from ED50 to WGS 84, and they disagree by metres — which is exactly why this package refuses to pick one for you. Binding is how you record the one your data was surveyed against.
from pyproj import Geod
from geodetic_engine.geodesy import transform
from geodetic_engine.persistablereference import parse_persistable_reference
# OSDU record Geographic2D:EPSG::4230 -- the CRS, on its own
ED50_GEOGRAPHIC = r"""
{"authCode":{"auth":"EPSG",
"code":"4230"},
"name":"GCS_European_1950",
"type":"LBC",
"ver":"PE_10_9_1",
"wkt":"GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433],AUTHORITY[\"EPSG\",4230]]"}
"""
# OSDU record EPSG::1612 -- ED50 to WGS 84 (23), seven parameters
SHIFT_1612 = r"""
{"authCode":{"auth":"EPSG",
"code":"1612"},
"name":"ED_1950_To_WGS_1984_23",
"type":"ST",
"ver":"PE_10_9_1",
"wkt":"GEOGTRAN[\"ED_1950_To_WGS_1984_23\",GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],METHOD[\"Position_Vector\"],PARAMETER[\"X_Axis_Translation\",-116.641],PARAMETER[\"Y_Axis_Translation\",-56.931],PARAMETER[\"Z_Axis_Translation\",-110.559],PARAMETER[\"X_Axis_Rotation\",0.893],PARAMETER[\"Y_Axis_Rotation\",0.921],PARAMETER[\"Z_Axis_Rotation\",-0.917],PARAMETER[\"Scale_Difference\",-3.52],OPERATIONACCURACY[1.0],AUTHORITY[\"EPSG\",1612]]"}
"""
# OSDU record EPSG::1133 -- ED50 to WGS 84 (1), three parameters
SHIFT_1133 = r"""
{"authCode":{"auth":"EPSG",
"code":"1133"},
"name":"ED_1950_To_WGS_1984_1",
"type":"ST",
"ver":"PE_10_9_1",
"wkt":"GEOGTRAN[\"ED_1950_To_WGS_1984_1\",GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],METHOD[\"Geocentric_Translation\"],PARAMETER[\"X_Axis_Translation\",-87.0],PARAMETER[\"Y_Axis_Translation\",-98.0],PARAMETER[\"Z_Axis_Translation\",-121.0],OPERATIONACCURACY[10.0],AUTHORITY[\"EPSG\",1133]]"}
"""
crs = parse_persistable_reference(ED50_GEOGRAPHIC)
point = (4.12789451, 63.58496782)
landed = {}
for payload in (SHIFT_1612, SHIFT_1133):
shift = parse_persistable_reference(payload)
bound = crs.to_bound_crs(shift)
result = transform(bound, "EPSG:4326", [point])
landed[str(shift.authority_code)] = result.coordinates[0]
print(
f"{shift.authority_code} {shift.name:<24} {shift.method_names[0]:<24} "
f"accuracy {shift.to_operation().accuracy} m"
)
print(f" -> {result.coordinates[0]}")
(first, second) = landed.values()
_, _, apart = Geod(ellps="WGS84").inv(first[0], first[1], second[0], second[1])
print(f"\nthe two published shifts disagree by {apart:.2f} m at this point")
EPSG:1612 ED_1950_To_WGS_1984_23 Position_Vector accuracy 1.0 m
-> (4.126133774855332, 63.58460392067997)
EPSG:1133 ED_1950_To_WGS_1984_1 Geocentric_Translation accuracy 10.0 m
-> (4.126052210412462, 63.58458429493627)
the two published shifts disagree by 4.60 m at this point
6.8 The package’s own types, built from a payload¶
Everything above works through parse_persistable_reference(), which returns a CrsReference, an OperationReference or a UnitReference — the payload’s own metadata, with to_crs() and to_operation() handing back pyproj objects.
Most of the time you want this package’s types instead, and a payload is accepted wherever a CRS definition is:
CoordinateReferenceSystem.from_persistable_reference()builds the package’s CRS, which adds the axis-order reporting a payload cannot state for itself.Transformation(source_crs=..., target_crs=...)takes payloads directly for either end, so an OSDU record transforms without being converted to codes first.
The three routes are layers over the same definition, not alternatives: the reference is the payload read, the CRS is that wrapped, and the transformation is two of those put to work.
Two things in the output below are worth reading twice. The reference is stamped OSDU:23032023 while the CRS reports authority_code None — that is the point of the whole module, not a gap: the code is provenance carried on the payload, and the CRS was built from the parameters rather than looked up, so it claims no register identity. And the target’s axes are declared Lat, Lon but its values are given as Lon, Lat, which is the sort of thing ESRI WKT cannot tell you and a caller should never have to guess.
from geodetic_engine.geodesy import CoordinateReferenceSystem, Transformation
# OSDU record BoundProjected:EPSG::23032_EPSG::1612
ED50_UTM32N_VIA_1612 = r"""
{"authCode":{"auth":"OSDU",
"code":"23032023"},
"lateBoundCRS":{"authCode":{"auth":"EPSG","code":"23032"},
"name":"ED_1950_UTM_Zone_32N",
"type":"LBC",
"ver":"PE_10_9_1",
"wkt":"PROJCS[\"ED_1950_UTM_Zone_32N\",GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"False_Easting\",500000.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",9.0],PARAMETER[\"Scale_Factor\",0.9996],PARAMETER[\"Latitude_Of_Origin\",0.0],UNIT[\"Meter\",1.0],AUTHORITY[\"EPSG\",23032]]"},
"name":"ED50 * EPSG-Nor N62 2001 / UTM zone 32N [23032,1612]",
"singleCT":{"authCode":{"auth":"EPSG","code":"1612"},
"name":"ED_1950_To_WGS_1984_23",
"type":"ST",
"ver":"PE_10_9_1",
"wkt":"GEOGTRAN[\"ED_1950_To_WGS_1984_23\",GEOGCS[\"GCS_European_1950\",DATUM[\"D_European_1950\",SPHEROID[\"International_1924\",6378388.0,297.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],METHOD[\"Position_Vector\"],PARAMETER[\"X_Axis_Translation\",-116.641],PARAMETER[\"Y_Axis_Translation\",-56.931],PARAMETER[\"Z_Axis_Translation\",-110.559],PARAMETER[\"X_Axis_Rotation\",0.893],PARAMETER[\"Y_Axis_Rotation\",0.921],PARAMETER[\"Z_Axis_Rotation\",-0.917],PARAMETER[\"Scale_Difference\",-3.52],OPERATIONACCURACY[1.0],AUTHORITY[\"EPSG\",1612]]"},
"type":"EBC",
"ver":"PE_10_9_1"}
"""
# OSDU record Geographic2D:EPSG::4326
WGS84_GEOGRAPHIC = r"""
{"authCode":{"auth":"EPSG",
"code":"4326"},
"name":"GCS_WGS_1984",
"type":"LBC",
"ver":"PE_10_9_1",
"wkt":"GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433],AUTHORITY[\"EPSG\",4326]]"}
"""
# 1. The same payload as this package's CRS.
source = CoordinateReferenceSystem.from_persistable_reference(ED50_UTM32N_VIA_1612)
target = CoordinateReferenceSystem.from_persistable_reference(WGS84_GEOGRAPHIC)
print("\nCoordinateReferenceSystem:")
print(f" name {source.name}")
print(f" authority_code {source.authority_code}")
# ESRI WKT declares no axes, so the order a payload's values arrive in is the
# thing worth asking about rather than assuming.
print(f"\n{'':<22}{'declared':<18}{'values given as'}")
for label, crs in (("source", source), ("target", target)):
print(f" {label:<19}{crs.axis_abbreviations!s:<18}{crs.value_axis_abbreviations}")
# 2. Two payloads put to work. Either end may be a payload, a code, or a CRS object.
survey = [(500000.0, 6600000.0), (612345.6, 6750000.0)]
ct = Transformation(source_crs=ED50_UTM32N_VIA_1612, target_crs=WGS84_GEOGRAPHIC)
result = ct.transform(survey)
print(f"\nTransformation -> {type(ct).__name__}")
print(f" operation applied {result.operation.name}")
for (east, north), point in zip(survey, result.coordinates, strict=True):
print(f" {east:>10.1f} E {north:>11.1f} N -> {point[0]:.9f}, {point[1]:.9f}")
# The same two payloads handed to CRS objects instead give the same answer.
same = Transformation(source_crs=source, target_crs=target).transform(survey)
print(
"\n via CRS objects rather than payloads: "
f"{same.coordinates == result.coordinates}"
)
CoordinateReferenceSystem:
name ED50 / UTM zone 32N
authority_code None
declared values given as
source ('E', 'N') ('E', 'N')
target ('Lat', 'Lon') ('Lon', 'Lat')
Transformation -> Transformation
operation applied ED_1950_To_WGS_1984_23
500000.0 E 6600000.0 N -> 8.998592587, 59.536498959
612345.6 E 6750000.0 N -> 11.067324637, 60.867452799
via CRS objects rather than payloads: True
6.9 Plain ESRI WKT, without the OSDU envelope¶
A persistableReference is ESRI WKT inside a JSON envelope. Where a register supplies WKT on its own, it can be given directly: a GEOGCS for either end and a GEOGTRAN as the operation.
The shift is built from the seven parameters below, not looked up by its authority code. ESRI’s GEOGTRAN is not standard WKT and PROJ will not read one directly, so this package translates it.
from geodetic_engine.geodesy import transform
# ED50 and WGS 84 as ESRI writes them: no authority code, nothing to look up.
ED50_WKT = (
'GEOGCS["GCS_European_1950",DATUM["D_European_1950",'
'SPHEROID["International_1924",6378388.0,297.0]],'
'PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]'
)
WGS84_WKT = (
'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",'
'SPHEROID["WGS_1984",6378137.0,298.257223563]],'
'PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]'
)
# ED50 to WGS 84 (23), EPSG:1612 -- seven parameters, as an ESRI GEOGTRAN.
EPSG_1612_WKT = (
'GEOGTRAN["ED_1950_To_WGS_1984_23",' + ED50_WKT + "," + WGS84_WKT + ","
'METHOD["Position_Vector"],'
'PARAMETER["X_Axis_Translation",-116.641],'
'PARAMETER["Y_Axis_Translation",-56.931],'
'PARAMETER["Z_Axis_Translation",-110.559],'
'PARAMETER["X_Axis_Rotation",0.893],'
'PARAMETER["Y_Axis_Rotation",0.921],'
'PARAMETER["Z_Axis_Rotation",-0.917],'
'PARAMETER["Scale_Difference",-3.52],'
'OPERATIONACCURACY[1.0],AUTHORITY["EPSG",1612]]'
)
# Inside EPSG:1612's area of use: Norway offshore, north of 62N.
point = (2.5, 63.5)
result = transform(ED50_WKT, WGS84_WKT, point, operation=EPSG_1612_WKT)
longitude, latitude = result.coordinates[0]
print(f"{point[0]}, {point[1]} -> {longitude:.8f}, {latitude:.8f}")
print(f"operation applied {result.operation.name}")
print(f"method {result.operation.method_name}")
2.5, 63.5 -> 2.49818948, 63.49961375
operation applied ED_1950_To_WGS_1984_23
method Position Vector transformation (geog2D domain)
6.10 The same transformation between two projected CRSs¶
from geodetic_engine.geodesy import transform
# The geographic frames, needed inside each PROJCS and inside the GEOGTRAN.
ED50_WKT = (
'GEOGCS["GCS_European_1950",DATUM["D_European_1950",'
'SPHEROID["International_1924",6378388.0,297.0]],'
'PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]'
)
WGS84_WKT = (
'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",'
'SPHEROID["WGS_1984",6378137.0,298.257223563]],'
'PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]'
)
# UTM zone 31N over each: same projection, different datum.
ED50_UTM31N_WKT = (
'PROJCS["ED_1950_UTM_Zone_31N",' + ED50_WKT + ","
'PROJECTION["Transverse_Mercator"],'
'PARAMETER["False_Easting",500000.0],'
'PARAMETER["False_Northing",0.0],'
'PARAMETER["Central_Meridian",3.0],'
'PARAMETER["Scale_Factor",0.9996],'
'PARAMETER["Latitude_Of_Origin",0.0],'
'UNIT["Meter",1.0]]'
)
WGS84_UTM31N_WKT = (
'PROJCS["WGS_1984_UTM_Zone_31N",' + WGS84_WKT + ","
'PROJECTION["Transverse_Mercator"],'
'PARAMETER["False_Easting",500000.0],'
'PARAMETER["False_Northing",0.0],'
'PARAMETER["Central_Meridian",3.0],'
'PARAMETER["Scale_Factor",0.9996],'
'PARAMETER["Latitude_Of_Origin",0.0],'
'UNIT["Meter",1.0]]'
)
# ED50 to WGS 84 (23), EPSG:1612 -- between the geographic frames, as before.
EPSG_1612_WKT = (
'GEOGTRAN["ED_1950_To_WGS_1984_23",' + ED50_WKT + "," + WGS84_WKT + ","
'METHOD["Position_Vector"],'
'PARAMETER["X_Axis_Translation",-116.641],'
'PARAMETER["Y_Axis_Translation",-56.931],'
'PARAMETER["Z_Axis_Translation",-110.559],'
'PARAMETER["X_Axis_Rotation",0.893],'
'PARAMETER["Y_Axis_Rotation",0.921],'
'PARAMETER["Z_Axis_Rotation",-0.917],'
'PARAMETER["Scale_Difference",-3.52],'
'OPERATIONACCURACY[1.0],AUTHORITY["EPSG",1612]]'
)
# Eastings and northings in zone 31N, north of 62N: inside EPSG:1612's area.
survey = [(475000.0, 7040000.0), (512500.0, 7125000.0)]
result = transform(ED50_UTM31N_WKT, WGS84_UTM31N_WKT, survey, operation=EPSG_1612_WKT)
print(f"operation applied {result.operation.name}")
print(f"method {result.operation.method_name}\n")
print(f"{'ED50 / UTM 31N':<28}{'WGS 84 / UTM 31N':<28}moved")
for (east, north), (new_east, new_north) in zip(
survey, result.coordinates, strict=True
):
moved = ((new_east - east) ** 2 + (new_north - north) ** 2) ** 0.5
print(
f"{east:>11.4f} {north:>13.4f} "
f"{new_east:>11.4f} {new_north:>13.4f} {moved:6.4f} m"
)
operation applied ED_1950_To_WGS_1984_23
method Position Vector transformation (geog2D domain)
ED50 / UTM 31N WGS 84 / UTM 31N moved
475000.0000 7040000.0000 474910.7896 7039785.0135 232.7610 m
512500.0000 7125000.0000 512410.8083 7124784.5996 233.1362 m