# Geodetic transformation¶

PROJ can do everything from the most simple projection to very complex transformations across many reference frames. While originally developed as a tool for cartographic projections, PROJ has over time evolved into a powerful generic coordinate transformation engine that makes it possible to do both large scale cartographic projections as well as coordinate transformation at a geodetic high precision level. This chapter delves into the details of how geodetic transformations of varying complexity can be performed.

In PROJ, two frameworks for geodetic transformations exists, the
*PROJ 4.x/5.x* / cs2cs / `pj_transform()`

framework and the *transformation pipelines* framework. The first is the original,
and limited, framework for doing geodetic transforms in PROJ. The latter is a
newer addition that aims to be a more complete transformation framework. Both are
described in the sections below. Large portions of the text are based on
[EversKnudsen2017].

Before describing the details of the two frameworks, let us first remark that most cases of geodetic transformations can be expressed as a series of elementary operations, the output of one operation being the input of the next. E.g. when going from UTM zone 32, datum ED50, to UTM zone 32, datum ETRS89, one must, in the simplest case, go through 5 steps:

Back-project the UTM coordinates to geographic coordinates

Convert the geographic coordinates to 3D cartesian geocentric coordinates

Apply a Helmert transformation from ED50 to ETRS89

Convert back from cartesian to geographic coordinates

Finally project the geographic coordinates to UTM zone 32 planar coordinates.

## Transformation pipelines¶

The homology between the above steps and a Unix shell style pipeline is evident. It is there the main architectural inspiration behind the transformation pipeline framework. The pipeline framework is realized by utilizing a special "projection", that takes as its user supplied arguments, a series of elementary operations, which it strings together in order to implement the full transformation needed. Additionally, a number of elementary geodetic operations, including Helmert transformations, general high order polynomial shifts and the Molodensky transformation are available as part of the pipeline framework. In anticipation of upcoming support for full time-varying transformations, we also introduce a 4D spatiotemporal data type, and a programming interface (API) for handling this.

The Molodensky transformation converts directly from geodetic coordinates
in one datum, to geodetic coordinates in another datum, while the (typically more
accurate) Helmert transformation converts from 3D cartesian to 3D cartesian
coordinates. So when using the Helmert transformation one typically needs to do an
initial conversion from geodetic to cartesian coordinates, and a final conversion
the other way round, to arrive at the desired result. Fortunately, this three-step
compound transformation has the attractive characteristic that each step depends
only on the output of the immediately preceding step. Hence, we can build a
geodetic-to-geodetic Helmert transformation by tying together the outputs and inputs
of 3 steps (geodetic-to-cartesian → Helmert → cartesian-to-geodetic), pipeline style.
The pipeline driver, makes this kind of chained transformations possible.
The implementation is compact, consisting of just one pseudo-projection, called
`pipeline`

, which takes as its arguments strings of elementary projections
(note: "projection" is the, slightly misleading, PROJ term used for any kind of
transformation).
The pipeline pseudo projection is supplemented by a number of elementary
transformations, all in all providing a framework for building high accuracy
solutions for a wide spectrum of geodetic tasks.

As a first example, let us take a look at the iconic
*geodetic → Cartesian → Helmert → geodetic* case (steps 2 to 4 in the example in
the introduction).
(Note, for the rest of this book the "+"s are removed from the examples.)
In PROJ it can be implemented as

```
proj=pipeline
step proj=cart ellps=intl
step proj=helmert convention=coordinate_frame
x=-81.0703 y=-89.3603 z=-115.7526
rx=-0.48488 ry=-0.02436 rz=-0.41321 s=-0.540645
step proj=cart inv ellps=GRS80
```

The pipeline can be expanded at both ends to accommodate whatever coordinate type is needed for input and output: In the example below, we transform from the deprecated Danish System 45, a 2D system with some tension in the original defining network, to UTM zone 33, ETRS89. The tension is reduced using a polynomial transformation (the init=./s45b... step, s45b.pol is a file containing the polynomial coefficients), taking the S45 coordinates to a technical coordinate system (TC32), defined to represent "UTM zone 32 coordinates, as they would look if the Helmert transformation between ED50 and ETRS89 was perfect". The TC32 coordinates are then converted back to geodetic(ED50) coordinates, using an inverse UTM projection, further to cartesian(ED50), then to cartesian(ETRS89), using the relevant Helmert transformation, and back to geodetic(ETRS89), before finally being projected onto the UTM zone 33, ETRS89 system. All in all a 6 step pipeline, implementing a transformation with centimeter level accuracy from a deprecated system with decimeter level tensions.

```
proj=pipeline
step init=./s45b.pol:s45b_tc32
step proj=utm inv ellps=intl zone=32
step proj=cart ellps=intl
step proj=helmert convention=coordinate_frame
x=-81.0703 y=-89.3603 z=-115.7526
rx=-0.48488 ry=-0.02436 rz=-0.41321 s=-0.540645
step proj=cart inv ellps=GRS80
step proj=utm ellps=GRS80 zone=33
```

With the pipeline framework spatiotemporal transformation is possible. This is possible by leveraging the time dimension in PROJ that enables 4D coordinates (three spatial components and one temporal component) to be passed through a transformation pipeline. In the example below a transformation from ITRF93 to ITRF2000 is defined. The temporal component is given as GPS weeks in the input data, but the 14-parameter Helmert transform expects temporal units in decimalyears. Hence the first step in the pipeline is the unitconvert pseudo-projection that makes sure the correct units are passed along to the Helmert transform. Most parameters of the Helmert transform are taken from [Altamimi2002], except the epoch which is the epoch of the transformation. The last step in the pipeline is converting the coordinate timestamps back to GPS weeks.

```
proj=pipeline
step proj=unitconvert t_in=gps_week t_out=decimalyear
step proj=helmert convention=coordinate_frame
x=0.0127 y=0.0065 z=-0.0209 s=0.00195
rx=0.00039 ry=-0.00080 rz=0.00114
dx=-0.0029 dy=-0.0002 dz=-0.0006 ds=0.00001
drx=0.00011 dry=0.00019 drz=-0.00007
t_epoch=1988.0
step proj=unitconvert t_in=decimalyear t_out=gps_week
```

## PROJ 4.x/5.x paradigm¶

Parameter

Description

+datum

Datum name (see

`proj -ld`

)+geoidgrids

Filename of GTX grid file to use for vertical datum transforms

+geoid_crs

(Added in PROJ 9.1) Value can be "WGS84" or "horizontal_crs". Only used if geoidgrids is specified.

+nadgrids

Filename of NTv2 grid file to use for datum transforms

+towgs84

3 or 7 term datum transform parameters

+to_meter

Multiplier to convert map units to 1.0m

+vto_meter

Vertical conversion to meters

Warning

This section documents the behavior of PROJ 4.x and 5.x. In PROJ 6.x,
cs2cs has been reworked to use `proj_create_crs_to_crs()`

internally,
with *late binding* capabilities, and thus is no longer constrained to using
WGS84 as a pivot (also called as *early binding* method).
When cs2cs of PROJ 6 is used with PROJ.4 expanded strings to describe the CRS,
including `+towgs84`

, `+nadgrids`

and `+geoidgrids`

, it will generally give
the same results as earlier PROJ versions. When used with AUTHORITY:CODE
CRS descriptions, it may return different results.

The *cs2cs* framework in PROJ 4 and 5 delivers a subset of the geodetic transformations available
with the *pipeline* framework. Coordinate transformations done in this framework
were transformed in a two-step process with WGS84 as a pivot datum. That is, the
input coordinates are transformed to WGS84 geodetic coordinates and then transformed
from WGS84 coordinates to the specified output coordinate reference system, by
utilizing either the Helmert transform, datum shift grids or a combination of both.
Datum shifts can be described in a proj-string with the parameters `+towgs84`

,
`+nadgrids`

and `+geoidgrids`

.
An inverse transform exists for all three and is applied if
specified in the input proj-string. The most common is `+towgs84`

, which is used to
define a 3- or 7-parameter Helmert shift from the input reference frame to WGS84.
Exactly which realization of WGS84 is not specified, hence a fair amount of
uncertainty is introduced in this step of the transformation. With the +nadgrids
parameter a non-linear planar correction derived from interpolation in a
correction grid can be applied. Originally this was implemented as a means to
transform coordinates between the North American datums NAD27 and NAD83, but
corrections can be applied for any datum for which a correction grid exists. The
inverse transform for the horizontal grid shift is "dumb", in the sense that the
correction grid is applied verbatim without taking into account that the inverse
operation is non-linear. Similar to the horizontal grid correction, `+geoidgrids`

can be used to perform grid corrections in the vertical component.
Both grid correction methods allow inclusion of more than one grid in the same
transformation

The `+geoid_crs=WGS84/horizontal_crs`

parameter was added in PROJ 9.1 to solve
an ambiguity that has been rampant for a long time. Before PROJ 6, the interpolation
CRS of the geoid transformation was assumed to be the one of the horizontal CRS.
Starting with PROJ 6, this was inadvertently changed to be WGS84.

When transforming from a CRS that specifies a horizontal datum shift
with `+towgs84`

or `+nadgrids`

, and a vertical one with `+geoidgrids`

with
`+geoid_crs=horizontal_crs`

, to WGS84, the vertical transformation is done before
the horizontal one.

```
$ projinfo -s "+proj=longlat +ellps=GRS80 +nadgrids=@foo +geoidgrids=@bar +geoid_crs=horizontal_crs +type=crs" -t EPSG:4979 -o PROJ -q
+proj=pipeline
+step +proj=unitconvert +xy_in=deg +xy_out=rad
+step +proj=vgridshift +grids=@bar +multiplier=1
+step +proj=hgridshift +grids=@foo
+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m
+step +proj=axisswap +order=2,1
```

Conversely, when `+geoid_crs=WGS84`

(or omitting it with PROJ >= 6),
the horizontal transformation is done before the vertical one.

```
$ projinfo -s "+proj=longlat +ellps=GRS80 +nadgrids=@foo +geoidgrids=@bar +geoid_crs=WGS84 +type=crs" -t EPSG:4979 -o PROJ -q
+proj=pipeline
+step +proj=unitconvert +xy_in=deg +xy_out=rad
+step +proj=hgridshift +grids=@foo
+step +proj=vgridshift +grids=@bar +multiplier=1
+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1
```

In contrast to the *transformation pipeline* framework, transformations with the
*cs2cs* framework in PROJ 4 and 5 were expressed as two separate proj-strings. One proj-string *to*
WGS84 and one *from* WGS84. Together they form the mapping from the source
coordinate reference system to the destination coordinate reference system.
When used with the cs2cs the source and destination CRS's are separated by the
special `+to`

parameter.

The following example demonstrates converting from the Greek GGRS87 datum
to WGS84 with the `+towgs84`

parameter.

```
cs2cs +proj=latlong +ellps=GRS80 +towgs84=-199.87,74.79,246.62
+to +proj=latlong +datum=WGS84
20 35
20d0'5.467"E 35d0'9.575"N 0.000
```

With PROJ 6, you can simply use the following:

Note

With PROJ 6, the order of coordinates for EPSG geographic coordinate reference systems is latitude first, longitude second.

```
cs2cs "GGRS87" "WGS 84"
35 20
35d0'9.575"N 20d0'5.467"E 0.000
cs2cs EPSG:4121 EPSG:4326
35 20
35d0'9.575"N 20d0'5.467"E 0.000
```

The EPSG database provides this example for transforming from WGS72 to WGS84 using an approximated 7 parameter transformation.

```
cs2cs +proj=latlong +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.219 \
+to +proj=latlong +datum=WGS84
4 55
4d0'0.554"E 55d0'0.09"N 0.000
```

With PROJ 6, you can simply use the following (note the reversed order for latitude and longitude)

```
cs2cs "WGS 72" "WGS 84"
55 4
55d0'0.09"N 4d0'0.554"E 0.000
cs2cs EPSG:4322 EPSG:4326
55 4
55d0'0.09"N 4d0'0.554"E 0.000
```

## Grid Based Datum Adjustments¶

In many places (notably North America and Australia) national geodetic organizations provide grid shift files for converting between different datums, such as NAD27 to NAD83. These grid shift files include a shift to be applied at each grid location. Actually grid shifts are normally computed based on an interpolation between the containing four grid points.

PROJ supports use of grid files for shifting between various reference frames.
The grid shift table formats are CTable, NTv1 (the old Canadian format), and NTv2
(`.gsb`

- the new Canadian and Australian format).

The text in this section is based on the *cs2cs* framework. Gridshifting is off
course also possible with the *pipeline* framework. The major difference between the
two is that the *cs2cs* framework is limited to grid mappings to WGS84, whereas with
*transformation pipelines* it is possible to perform grid shifts between any two
reference frames, as long as a grid exists.

Use of grid shifts with cs2cs is specified using the `+nadgrids`

keyword in a coordinate system definition. For example:

```
% cs2cs +proj=latlong +ellps=clrk66 +nadgrids=ntv1_can.dat \
+to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF
-111 50
EOF
111d0'2.952"W 50d0'0.111"N 0.000
```

In this case the `/usr/local/share/proj/ntv1_can.dat`

grid shift file was
loaded, and used to get a grid shift value for the selected point.

It is possible to list multiple grid shift files, in which case each will be tried in turn till one is found that contains the point being transformed.

```
cs2cs +proj=latlong +ellps=clrk66 \
+nadgrids=conus,alaska,hawaii,stgeorge,stlrnc,stpaul \
+to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF
-111 44
EOF
111d0'2.788"W 43d59'59.725"N 0.000
```

### Skipping Missing Grids¶

The special prefix `@`

may be prefixed to a grid to make it optional. If it
not found, the search will continue to the next grid. Normally any grid not
found will cause an error. For instance, the following would use the
`ntv2_0.gsb`

file if available, otherwise it would
fallback to using the `ntv1_can.dat`

file.

```
cs2cs +proj=latlong +ellps=clrk66 +nadgrids=@ntv2_0.gsb,ntv1_can.dat \
+to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF
-111 50
EOF
111d0'3.006"W 50d0'0.103"N 0.000
```

### The null Grid¶

A special `null`

grid shift file is distributed with PROJ.
This file provides a zero shift for the whole world. It may be
listed at the end of a nadgrids file list if you want a zero shift to be
applied to points outside the valid region of all the other grids. Normally if
no grid is found that contains the point to be transformed an error will occur.

```
cs2cs +proj=latlong +ellps=clrk66 +nadgrids=conus,null \
+to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF
-111 45
EOF
111d0'3.006"W 50d0'0.103"N 0.000
cs2cs +proj=latlong +ellps=clrk66 +nadgrids=conus,null \
+to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF
-111 44
-111 55
EOF
111d0'2.788"W 43d59'59.725"N 0.000
111dW 55dN 0.000
```

For more information see the chapter on Other transformation grids.

### Caveats¶

Where grids overlap (such as conus and

`ntv1_can.dat`

for instance) the first found for a point will be used regardless of whether it is appropriate or not. So, for instance,`+nadgrids=ntv1_can.dat`

,conus would result in the Canadian data being used for some areas in the northern United States even though the conus data is the approved data to use for the area. Careful selection of files and file order is necessary. In some cases border spanning datasets may need to be pre-segmented into Canadian and American points so they can be properly grid shiftedAdditional detail on the grid shift being applied can be found by setting the PROJ_DEBUG environment variable to a value. This will result in output to stderr on what grid is used to shift points, the bounds of the various grids loaded and so forth