Conversions and transformations¶
Every example on this page runs during the documentation build, with PROJ
9.9.0 and the proj-data grids installed in the devcontainer
image.
A conversion stays on one datum, such as a map projection. There is one answer, so no operation needs naming.
A transformation changes datum. Several are usually published, so you must say which one (Core concepts).
Conversions¶
Map projection¶
WGS 84 geographic (lon, lat in degrees) to WGS 84 / UTM zone 32N (E, N in
metres). Same datum, so no operation is named. The result records the
conversion PROJ chose, EPSG:16032 UTM zone 32N:
from geodetic_engine.geodesy import transform
result = transform("EPSG:4326", "EPSG:32632", (10.7522, 59.9139))
print(result.coordinates)
print(result.operation.authority_code, result.operation.name, "/", result.operation.route)
((597979.9028826989, 6643118.991493065),)
EPSG:16032 UTM zone 32N / proj_default
Inverse projection, and projection to projection¶
A projected CRS can be the source, and two projections on the same datum convert directly, through the geographic CRS they share:
# ETRS89 / UTM 32N -> ETRS89 geographic (lon, lat)
print(transform("EPSG:25832", "EPSG:4258", (597979.90, 6643118.99)).coordinates)
# ETRS89 / UTM 32N -> ETRS89 / UTM 33N: same datum, different zone
print(transform("EPSG:25832", "EPSG:25833", (597979.90, 6643118.99)).coordinates)
((10.752199947806304, 59.91389998838889),)
((262560.4791834648, 6649443.582868823),)
Geographic to geocentric¶
A 3D geographic CRS (lon, lat, ellipsoidal height) to Earth-centred X, Y, Z in metres, on the same datum:
transform("EPSG:4979", "EPSG:4978", (10.7522, 59.9139, 100.0)).coordinates
((3149181.0076383473, 598015.9749511877, 5495761.153737682),)
Unit change¶
EPSG:2278 is NAD83 / Texas South Central in US survey feet. EPSG:32139 is
the same projection in metres. Only the units differ:
feet = transform("EPSG:32139", "EPSG:2278", (900000.0, 4200000.0))
print(feet.coordinates, feet.target_units)
((2249815.656301068, 17731323.118244793),) ('US survey foot', 'US survey foot')
Transformations¶
Named operation¶
ED50 to WGS 84 in Norwegian waters north of 62°N, using EPSG:1612, a
seven-parameter Helmert with a stated accuracy of 1 m. Input and output are
(lon, lat, h), in degrees and metres:
result = transform("EPSG:4230", "EPSG:4326", (2.5, 63.5, 100.0), operation="EPSG:1612")
print(result.coordinates)
print(result.operation.name, "|", result.operation.method_name, "|", result.operation.accuracy, "m")
((2.4981894919143337, 63.49961375156072, 100.0),)
ED50 to WGS 84 (23) | Position Vector transformation (geog2D domain) | 1.0 m
An operation can be named as "EPSG:1612", 1612, an OGC URN, or an
OperationCandidate returned by
available_operations()
(Choosing an operation).
A datum change and a projection in one call¶
Named operations are wrapped in whatever conversions the two CRSs need. Here ED50 geographic goes to WGS 84 / UTM zone 32N. The result reports the datum shift, not the projection, because the shift is what you asked for and what determines the accuracy:
result = transform("EPSG:4230", "EPSG:32632", (11.12789451, 63.58496782), operation="EPSG:1612")
print(result.coordinates)
print(result.operation.authority_code, result.operation.route)
print(result.operation.steps)
((605532.3187915211, 7052489.142085133),)
EPSG:1612 transformer_group
('ED50 to WGS 84 (23)', 'UTM zone 32N', 'axis order change (2D)')
Concatenated operations, and naming each step¶
EPSG:8047 (ED50 to WGS 84 (15)) is published as two steps: EPSG:1147 (ED50
to ED87) then EPSG:1146 (ED87 to WGS 84). Naming the concatenated code, or
naming both steps, gives the same result:
point = (4.12789451, 63.58496782, 100.0)
by_code = transform("EPSG:4230", "EPSG:4326", point, operation="EPSG:8047")
by_steps = transform("EPSG:4230", "EPSG:4326", point, operation=["EPSG:1147", "EPSG:1146"])
print(by_code.coordinates)
print(by_steps.coordinates)
print(by_code.operation.authority_code, "vs", by_steps.operation.steps)
((4.126139926951677, 63.584613419996394, 100.0),)
((4.126139926951677, 63.584613419996394, 100.0),)
EPSG:8047 vs ('ED50 to ED87 (2)', 'ED87 to WGS 84 (1)', 'axis order change (2D)')
The steps are a set, not a sequence. Each is checked separately against the pipeline PROJ built, so their order does not matter:
reordered = transform("EPSG:4230", "EPSG:4326", point, operation=["EPSG:1146", "EPSG:1147"])
reordered.coordinates == by_steps.coordinates
True
Grid-based horizontal transformation¶
NAD27 to NAD83 over the conterminous US, using EPSG:1241 (NADCON). PROJ reads
the grid us_noaa_conus.tif. The result names every grid it read, and whether
each was installed:
result = transform("EPSG:4267", "EPSG:4269", (-95.0, 30.0), operation="EPSG:1241")
print(result.coordinates)
for grid in result.grids:
print(grid.name, "available" if grid.available else "MISSING", grid.full_name)
((-95.00020486027408, 30.00021830000477),)
us_noaa_conus.tif available /usr/local/share/proj/us_noaa_conus.tif
The newer NADCON5 transformation EPSG:8555 reads a different grid, and
differs by a few centimetres:
nadcon5 = transform("EPSG:4267", "EPSG:4269", (-95.0, 30.0), operation="EPSG:8555")
print(nadcon5.coordinates, [g.name for g in nadcon5.grids])
((-95.00020430775153, 30.00021849426958),) ['us_noaa_nadcon5_nad27_nad83_1986_conus.tif']
If a grid is not installed, the transformation is refused and the missing file is named (MissingGridError).
Vertical: ellipsoidal height to a geoid height¶
WGS 84 3D (lon, lat, ellipsoidal height in metres) to EGM2008 height, using
EPSG:3858 and the 2.5′ geoid grid. The target is a vertical CRS, so each
point has one output value, the height:
result = transform("EPSG:4979", "EPSG:3855", (-144.0, 72.0, 548.4082), operation="EPSG:3858")
print(result.coordinates, result.target_axes, result.target_units)
print([g.name for g in result.grids])
((556.3834421379089,),) ('H',) ('metre',)
['us_nga_egm08_25.tif']
Compound target: position and height together¶
EPSG:6172 is ETRS89 / UTM zone 32N + NN54 height. From ETRS89 3D
(EPSG:4937), only the vertical part needs a transformation, EPSG:11559,
which reads the Norwegian grid no_kv_href2008a.tif. Output is (E, N, H) in
metres:
result = transform(
"EPSG:4937", "EPSG:6172", (11.12789451, 63.58496782, 100.0), operation="EPSG:11559"
)
print(result.coordinates, result.target_axes)
((605606.253274944, 7052523.904230434, 61.741534204810534),) ('E', 'N', 'H')
From WGS 84 3D (EPSG:4979) there is also a horizontal datum change, so both
operations must be named. PROJ merges the horizontal and vertical steps into
one unidentified operation, so naming only one would leave the other chosen by
PROJ:
result = transform(
"EPSG:4979", "EPSG:6172", (11.12789451, 63.58496782, 100.0),
operation=["EPSG:11028", "EPSG:11559"],
)
print(result.coordinates)
print(result.operation.name)
((605606.253274944, 7052523.904230434, 61.741534204810534),)
Inverse of ETRS89-NOR [EUREF89] to WGS 84 (1) + ETRS89-NOR [EUREF89] to NN54 height (1) + UTM zone 32N
Naming just the vertical one is refused:
transform("EPSG:4979", "EPSG:6172", (11.12789451, 63.58496782, 100.0), operation="EPSG:11559")
OperationNotAvailableError: EPSG:11559 operates between 'ETRS89-NOR [EUREF89] (with axis order normalized for visualization)' and 'NN54 height', neither of which shares a datum with EPSG:4979; it cannot be applied here
Time-dependent: a coordinate epoch¶
EPSG:6277 (ITRF2005 to GDA94) is a 14-parameter Helmert with rates of
change, so its result depends on when the coordinates were observed. Pass the
coordinate epoch as a decimal year. Input and output are geocentric X, Y, Z in
metres:
from geodetic_engine.geodesy import Transformation
tfm = Transformation("EPSG:4896", "EPSG:4938", operation="EPSG:6277")
print("requires epoch:", tfm.requires_epoch)
point = [(-2593197.524, 5656917.6189, -1394397.8828)]
for epoch in (1994.0, 2024.0):
print(epoch, tfm.transform(point, coordinate_epoch=epoch).coordinates)
requires epoch: True
1994.0 ((-2593197.5478785704, 5656917.676734876, -1394397.8797274325),)
2024.0 ((-2593196.308460863, 5656917.8510818165, -1394399.550457241),)
Over 30 years the result moves by about 2 m, which is the Australian plate’s motion relative to ITRF. Leaving the epoch out is refused (MissingCoordinateEpochError).
A dynamic CRS alone does not require an epoch. WGS 72 to WGS 84 (EPSG:1238)
is a static Helmert and gives the same answer at every epoch:
Transformation("EPSG:4322", "EPSG:4326", operation="EPSG:1238").requires_epoch
False
Engineering CRS¶
Engineering CRSs are local plant or site grids. EPSG:5817 is the Tombak LNG
plant grid, with axes x and y in metres. EPSG:15747 is a similarity
transformation to Nakhl-e Ghanem / UTM zone 39N:
result = transform("EPSG:5817", "EPSG:3307", (20000.0, 10000.0), operation="EPSG:15747")
print(result.coordinates, result.source_axes, "->", result.target_axes)
print(result.operation.method_name)
back = transform("EPSG:3307", "EPSG:5817", result.coordinates[0], operation="EPSG:15747")
print(back.coordinates)
((618336.7480000182, 3067774.210000054),) ('x', 'y') -> ('E', 'N')
Similarity transformation
((20000.000000000022, 9999.999999999927),)
A plant grid’s axes need not be east and north. EPSG:5800, the Astra Minas
grid, declares X north and Y west, so its values are given in that declared
order – there is no easting to put first – while the projected result is
xy as always. PROJ on its own gets this pair wrong (see Known issues and workarounds);
the point below lands where it should, at Comodoro Rivadavia:
result = transform("EPSG:5800", "EPSG:22192", (10000.0, 20000.0), operation="EPSG:1035")
print(result.coordinates, result.source_axes, "->", result.target_crs.value_axis_abbreviations)
print(transform("EPSG:5800", "EPSG:4221", (10000.0, 20000.0), operation="EPSG:1035").coordinates)
((2590394.630374936, 4915661.955434943),) ('X', 'Y') -> ('Y', 'X')
((-67.83504442706683, -45.90809129814133),)
Bound CRS: the CRS names its own transformation¶
A bound CRS carries its transformation to a hub, so no operation needs naming. It can be built with pyproj, parsed from an OSDU payload (Persistable references (OSDU)), or looked up by code in a custom database:
from pyproj import CRS
from pyproj.crs import CoordinateOperation
from pyproj.crs.crs import BoundCRS
ed50_utm31_via_1612 = BoundCRS(
source_crs=CRS.from_epsg(23031), # ED50 / UTM zone 31N
target_crs=CRS.from_epsg(4326), # hub: WGS 84
transformation=CoordinateOperation.from_epsg(1612),
)
result = transform(ed50_utm31_via_1612, "EPSG:4326", (500000.0, 7000000.0))
print(result.coordinates)
print(result.operation.authority_code, result.operation.route)
((2.9982275525073256, 63.12741118011052),)
EPSG:1612 bound
The route is bound: the operation came from the CRS definition, not from a
search. Bound CRSs covers bound CRSs over concatenated
operations and why looking them up by code needs a workaround.
A stated operation: OSDU payloads and ESRI GEOGTRAN¶
Instead of a code, the operation can be given in full as an OSDU
persistableReference payload, an ESRI GEOGTRAN WKT string, or a parsed
OperationReference. It is then
applied exactly as written, parameters and all, not looked up in PROJ’s
database. See A stated operation.
Choosing between transform and Transformation¶
transform() caches its resolved transformations,
so repeating a call with the same CRSs and operation does not resolve again.
Transformation makes this explicit. It
resolves when constructed, so a bad operation, missing grid or ambiguous datum
change is raised before any coordinates are passed. The object can then be
reused:
tfm = Transformation("EPSG:4230", "EPSG:4326", operation="EPSG:1612")
print(tfm)
print(tfm.source_crs, "->", tfm.target_crs, "| grids:", tfm.grids, "| epoch:", tfm.requires_epoch)
for batch in ([(2.5, 63.5)], [(3.0, 64.0), (4.0, 65.0)]):
print(tfm.transform(batch).coordinates)
Transformation(EPSG:4230 -> EPSG:4326, operation='EPSG:1612')
CoordinateReferenceSystem('EPSG:4230') -> CoordinateReferenceSystem('EPSG:4326') | grids: () | epoch: False
((2.4981894757094416, 63.49961374934551),)
((2.9981763918698823, 63.99964022048465), (3.9981487744754873, 64.99969354462532))
Inverting a transformation means swapping source and target and naming the same operation. PROJ applies it in reverse.