# 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 *cs2cs*
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). 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
```

## cs2cs paradigm¶

Parameter Description +datum Datum name (see `proj -ld`

)+geoidgrids Filename of GTX grid file to use for vertical datum transforms +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

The *cs2cs* framework delivers a subset of the geodetic transformations available
with the *pipeline* framework. Coordinate transformations done in this framework
are 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

In contrast to the *transformation pipeline* framework, transformations with the
*cs2cs* framework are 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 8.570
```

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 3.223
```

## 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 shift with releases after 4.4.6 (not
inclusive). 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 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 shifted - There are additional grids for shifting between NAD83 and various HPGN versions of the NAD83 datum. Use of these haven’t been tried recently so you may encounter problems. The FL.lla, WO.lla, MD.lla, TN.lla and WI.lla are examples of high precision grid shifts. Take care!
- Additional 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
- The
*cs2cs*framework always assumes that grids contain a shift**to**NAD83 (essentially WGS84). Other types of grids can be used with the*pipeline*framework.