Move 'Coordinate operations' to top level of docs [skip ci]

6 files changed, 707 insertions, 0 deletions

diff --git a/docs/source/operations/transformations/deformation.rst b/docs/source/operations/transformations/deformation.rst
new file mode 100644
index 00000000..2489f18e
--- /dev/null
+++ b/docs/source/operations/transformations/deformation.rst
@@ -0,0 +1,148 @@
+.. _deformation:
+
+================================================================================
+Kinematic datum shifting utilizing a deformation model
+================================================================================
+
+Perform datum shifts means of a deformation/velocity model.
+
++-----------------+--------------------------------------------------------------------+
+| **Input type**  | Cartesian coordinates.                                             |
++-----------------+--------------------------------------------------------------------+
+| **Output type** | Cartesian coordinates.                                             |
++-----------------+--------------------------------------------------------------------+
+| **Options**                                                                          |
++-----------------+--------------------------------------------------------------------+
+| `xy_grids`      | Comma-separated list of grids to load.                             |
++-----------------+--------------------------------------------------------------------+
+| `z_grids`       | Comma-separated list of grids to load.                             |
++-----------------+--------------------------------------------------------------------+
+| `t_epoch`       | Central epoch of transformation. [decimalyear]. Only used in       |
+|                 | spatiotemporal transformations.                                    |
++-----------------+--------------------------------------------------------------------+
+| `t_obs`         | Observation time of coordinate(s). Mostly useful in 2D and 3D      |
+|                 | transformations. [decimalyear]. *Optional*.                        |
++-----------------+--------------------------------------------------------------------+
+
+The deformation operation is used to adjust coordinates for intraplate deformations.
+Usually the transformation parameters for regional plate-fixed reference frames such as
+the ETRS89 does not take intraplate deformation into account. It is assumed that
+tectonic plate of the region is rigid. Often times this is true, but near the plate
+boundary and in areas with post-glacial uplift the assumption breaks. Tntraplate
+deformations can be modelled and then applied to the coordinates so that
+they represent the physical world better. In PROJ this is done with the deformation
+operation.
+
+The horizontal grid is stored in CTable2 format and the vertical grid is stored in the
+GTX format. Both grids are expected to contain grid-values in units of mm/year.
+Details about the formats can be found in the GDAL documentation. GDAL both reads and
+writes both file formats. Using GDAL for construction of new grids is recommended.
+
+Example
+-------------------------------------------------------------------------------
+
+In [Häkli2016]_ coordinate transformation including a deformation model is described.
+The paper describes how coordinates from the global ITRFxx frames are transformed to the
+local Nordic realisations of ERTS89. Scandinavia is an area with significant post-glacial
+rebound. The deformations from the post-glacial uplift is not accounted for in the
+offcial ETRS89 transformations so in order to get accurate transformations in the Nordic
+countries it is necesarry to apply the deformation model. The transformation from ITRF2008
+to the Danish realisation of ETRS89 is in PROJ described as::
+
+    proj =  pipeline ellps = GRS80
+            # ITRF2008@t_obs -> ITRF2000@t_obs
+    step    init = ITRF2008:ITRF2000
+            # ITRF2000@t_obs -> ETRF2000@t_obs
+    step    proj=helmert t_epoch = 2000.0 transpose
+            x =  0.054  rx =  0.000891 drx =  8.1e-05
+            y =  0.051  ry =  0.00539  dry =  0.00049
+            z = -0.048  rz = -0.008712 drz = -0.000792
+            # ETRF2000@t_obs -> NKG_ETRF00@2000.0
+    step    proj = deformation t_epoch = 2000.0
+            xy_grids = ./nkgrf03vel_realigned_xy.ct2
+            z_grids  = ./nkgrf03vel_realigned_z.gtx
+            # NKG_ETRF@2000.0 -> ETRF92@2000.0
+    step    proj=helmert transpose s = -0.009420e
+            x = 0.03863 rx = 0.00617753
+            y = 0.147   ry = 5.064e-05
+            z = 0.02776 rz = 4.729e-05
+            # ETRF92@2000.0 -> ETRF92@1994.704
+    step    proj = deformation t_epoch = 1994.704 t_obs = 2000.0
+            xy_grids = ./nkgrf03vel_realigned_xy.ct2
+            z_grids  = ./nkgrf03vel_realigned_z.gtx
+
+From this we can see that the transformation from ITRF2008 to the Danish realisation of
+ETRS89 is a combination of Helmert transformations and adjustments with a deformation
+model. The first use of the deformation operation is::
+
+    proj = deformation t_epoch = 2000.0
+    xy_grids = ./nkgrf03vel_realigned_xy.ct2
+    z_grids  = ./nkgrf03vel_realigned_z.gtx
+
+Here we set the central epoch of the transformation, 2000.0. The observation epoch
+is expected as part of the input coordinate tuple. The deformation model is
+described by two grids, specified with *xy_grids* and *z_grids*. The first is
+the horizontal part of the model and the second is the vertical component.
+
+Mathematical description
+-------------------------------------------------------------------------------
+
+Mathematically speaking, application of a deformation model is simple. The deformation model is
+represented as a grid of velocities in three dimensions. Coordinate corrections are
+applied in cartesian space. For a given coordinate, :math:`(X, Y, Z)`, velocities
+:math:`(V_X, V_Y, V_Z)` can be interpolated from the gridded model. The time span
+between :math:`t_c` and :math:`t_{obs}` determine the magnitude of the coordinate
+correcton as seen in eq. :eq:`apply_velocity` below.
+
+.. math::
+    :label: apply_velocity
+
+    \begin{align}
+        \begin{pmatrix}
+            X \\
+            Y \\
+            Z \\
+        \end{pmatrix}_B =
+        \begin{pmatrix}
+            X \\
+            Y \\
+            Z \\
+        \end{pmatrix}_A +
+        (t_c - t_{obs})
+        \begin{pmatrix}
+            V_X \\
+            V_Y \\
+            V_Z \\
+        \end{pmatrix}
+    \end{align}
+
+Corrections are done in cartesian space.
+
+Coordinates of the gridded model are in ENU (east, north, up) space because it would
+otherwise require an enormous 3 dimensional grid to handle the corrections in cartesian
+space. Keeping the correction in lat/long space reduces the complexity of the grid
+significantly. Consequentyly though, the input coordinates needs to be converted to
+lat/long space when searching for corrections in the grid. This is done with *cart*
+operation. The converted grid corrections can then be applied to the input coordinates
+in cartesian space. The conversion from ENU space to cartesian space is done in the
+following way:
+
+.. math::
+    :label: enu2xyz
+
+    \begin{align}
+        \begin{pmatrix}
+            X \\
+            Y \\
+            Z \\
+        \end{pmatrix} =
+        \begin{pmatrix}
+            -\sin\phi \cos\lambda N - \sin\lambda E + \cos\phi \cos\lambda U \\
+            -\sin\phi \sin\lambda N + \sin\lambda E + \cos\phi \sin\lambda U \\
+                         \cos\phi N +                 \sin\phi U \\
+        \end{pmatrix}
+    \end{align}
+
+where :math:`\phi` and :math:`\lambda` are the latitude and longitude of the coordinate
+that is searched for in the grid. :math:`(E, N, U)` are the grid values in ENU-space and
+:math:`(X, Y, Z)` are the corrections converted to cartesian space.
diff --git a/docs/source/operations/transformations/helmert.rst b/docs/source/operations/transformations/helmert.rst
new file mode 100644
index 00000000..03e5db15
--- /dev/null
+++ b/docs/source/operations/transformations/helmert.rst
@@ -0,0 +1,388 @@
+.. _helmert:
+
+================================================================================
+Helmert transform
+================================================================================
+
+The Helmert transformation changes coordinates from one reference frame to
+anoether by means of 3-, 4-and 7-parameter shifts, or one of their 6-, 8- and
+14-parameter kinematic counterparts.
+
++-----------------+-------------------------------------------------------------------+
+| **Input type**  | Cartesian coordinates (spatial), decimalyears (temporal).         |
++-----------------+-------------------------------------------------------------------+
+| **Output type** | Cartesian coordinates (spatial), decimalyears (temporal).         |
++-----------------+-------------------------------------------------------------------+
+| **Options**                                                                         |
++----------------+--------------------------------------------------------------------+
+| `x`            | Translation of the x-axis [m]. *Optional*.                         |
++----------------+--------------------------------------------------------------------+
+| `y`            | Translation of the y-axis [m]. *Optional*.                         |
++----------------+--------------------------------------------------------------------+
+| `z`            | Translation of the z-axis. [m]. *Optional*.                        |
++----------------+--------------------------------------------------------------------+
+| `s`            | Scale factor [ppm]. *Optional*.                                    |
++----------------+--------------------------------------------------------------------+
+| `rx`           | X-axis rotation in the 3D Helmert [arc seconds]. *Optional*.       |
++----------------+--------------------------------------------------------------------+
+| `ry`           | Y-axis rotatoin in the 3D Helmert [arc seconds]. *Optional*.       |
++----------------+--------------------------------------------------------------------+
+| `rz`           | Z-axis rotation in the 3D Helmert [arc seconds]. *Optional*.       |
++----------------+--------------------------------------------------------------------+
+| `theta`        | Rotation angle in the 2D Helmert. [arc seconds]. *Optional*.       |
++----------------+--------------------------------------------------------------------+
+| `dx`           | Translation rate of the x-axis. [m/year]. *Optional*.              |
++----------------+--------------------------------------------------------------------+
+| `dy`           | Translation rate of the y-axis. [m/year]. *Optional*.              |
++----------------+--------------------------------------------------------------------+
+| `dz`           | Translation rate of the z-axis. [m/year]. *Optional*.              |
++----------------+--------------------------------------------------------------------+
+| `ds`           | Scale rate factor [ppm/year]. *Optional*.                          |
++----------------+--------------------------------------------------------------------+
+| `drx`          | Rotation rate of the x-axis [arc seconds/year]. *Optional*.        |
++----------------+--------------------------------------------------------------------+
+| `dry`          | Rotation rate of the y-axis [arc seconds/year]. *Optional*.        |
++----------------+--------------------------------------------------------------------+
+| `drz`          | Rotation rate of the y-axis [arc seconds/year]. *Optional*.        |
++----------------+--------------------------------------------------------------------+
+| `t_epoch`      | Central epoch of transformation. [decimalyear]. Only used in       |
+|                | spatiotemporal transformations. *Optional*.                        |
++----------------+--------------------------------------------------------------------+
+| `t_obs`        | Observation time of coordinate(s). Mostly useful in 2D and 3D      |
+|                | transformations. [decimalyear]. *Optional*.                        |
++----------------+--------------------------------------------------------------------+
+| `exact`        | Use exact transformation equations. *Optional*.                    |
++----------------+--------------------------------------------------------------------+
+| `transpose`    | Transpose rotation matrix. *Optional*.                             |
++----------------+--------------------------------------------------------------------+
+
+The Helmert transform, in all it's various incarnations, is used to perform reference
+frame shifts. The transformation operates in cartesian space. It can be used to transform
+planar coordinates from one datum to another, transform 3D cartesian
+coordinates from one static reference frame to another or it can be used to do fully
+kinematic transformations from global reference frames to local static frames.
+
+All of the parameters described in the table above are marked as optional. This is true
+as long as at least one parameter is defined in the setup of the transformation.
+The behaviour of the transformation depends on which parameters are used in the setup.
+For instance, if a rate of change paramater is specified a kinematic version of the
+transformation is used.
+
+The kinematic transformations require an observation time of the coordinate, as well
+as a central epoch for the transformation. The latter is usually documented
+alongside the rest of the transformation parameters for a given transformation.
+The central eopch is controlled with the parameter `t_epoch`. The observation
+time can either by stated as part of the coordinate when using PROJ's
+4D-functionality or it can be controlled in the transformation setup by the
+parameter `t_obs`. When `t_obs` is specified, all transformed coordinates are
+treated as if they have the same observation time.
+
+Examples
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+Transforming coordinates from NAD72 to NAD83 using the 4 parameter 2D Helmert:
+
+::
+
+    proj=helmert ellps=GRS80 x=-9597.3572 y=.6112
+                 s=0.304794780637 theta=-1.244048
+
+Simplified transformations from ITRF2008/IGS08 to ETRS89 using 7 parameters:
+
+::
+
+    proj=helmert ellps=GRS80 x=0.67678    y=0.65495   z=-0.52827
+                rx=-0.022742 ry=0.012667 rz=0.022704  s=-0.01070
+
+Transformation from `ITRF2000@2017.0`  to `ITRF93@2017.0` using 15 parameters:
+
+::
+
+    proj=helmert ellps=GRS80
+         x=0.0127     y=0.0065     z=-0.0209  s=0.00195
+         dx=-0.0029   dy=-0.0002   dz=-0.0006 ds=0.00001
+         rx=-0.00039  ry=0.00080   rz=-0.00114
+         drx=-0.00011 dry=-0.00019 drz=0.00007
+         epoch=1988.0 tobs=2017.0    transpose
+
+Mathematical description
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+In the notation used below, :math:`\hat{P}` is the rate of change of a given transformation
+parameter :math:`P`. :math:`\dot{P}` is the kinematically adjusted version of :math:`P`,
+described by
+
+.. math::
+    :label: propagation
+
+    \dot{P}= P + \hat{P}\left(t - t_{central}\right)
+
+where :math:`t` is the observation time of the coordinate and :math:`t_{central}` is
+the central epoch of the transformation. Equation :eq:`propagation` can be used to
+propagate all transformation parameters in time.
+
+Superscripts of vectors denote the reference frame the coordinates in the vector belong to.
+
+
+2D Helmert
+-------------------------------------------------------------------------------
+
+The simplest version of the Helmert transform is the 2D case. In the 2-dimensional
+case only the horizontal coordinates are changed. The coordinates can be
+translated, rotated and scale. Translation is controlled with the `x` and `y`
+parameters. The rotation is determined by `theta` and the scale is controlled with
+the `s` paramaters.
+
+.. note::
+
+    The scaling parameter `s` is unitless for the 2D Helmert, as opposed to the
+    3D version where the scaling parameter is given in units of ppm.
+
+Mathematically the 2D Helmert is described as:
+
+.. math::
+    :label: 4param
+
+    \begin{align}
+        \begin{bmatrix}
+            X \\
+            Y \\
+        \end{bmatrix}^B =
+        \begin{bmatrix}
+            T_x \\
+            T_y \\
+        \end{bmatrix} +
+        s
+        \begin{bmatrix}
+            \hphantom{-}\cos \theta & \sin \theta \\
+            -\sin \theta & \cos \theta \\
+        \end{bmatrix}
+        \begin{bmatrix}
+            X \\
+            Y \\
+        \end{bmatrix}^A
+    \end{align}
+
+
+:eq:`4param` can be extended to a time-varying kinematic version by
+adjusting the parameters with :eq:`propagation` to :eq:`4param`, which yields
+the kinematic 2D Helmert transform:
+
+.. math::
+    :label: 8param
+
+    \begin{align}
+        \begin{bmatrix}
+            X \\
+            Y \\
+        \end{bmatrix}^B =
+        \begin{bmatrix}
+            \dot{T_x} \\
+            \dot{T_y} \\
+        \end{bmatrix} +
+        s(t)
+        \begin{bmatrix}
+             \hphantom{-}\cos \dot{\theta} & \sin \dot{\theta}  \\
+                        -\sin\ \dot{\theta} & \cos \dot{\theta} \\
+        \end{bmatrix}
+        \begin{bmatrix}
+            X \\
+            Y \\
+        \end{bmatrix}^A
+    \end{align}
+
+All parameters in :eq:`8param` are determined by the use of :eq:`propagation`,
+which applies the rate of change to each individual parameter for a given
+timespan between :math:`t` and :math:`t_{central}`.
+
+
+3D Helmert
+-------------------------------------------------------------------------------
+
+The general form of the 3D Helmert is
+
+.. math::
+    :label: general-helmert
+
+
+    \begin{align}
+        V^B = T + \left(1 + s \times 10^{-6}\right) \mathbf{R} V^A
+    \end{align}
+
+Where :math:`T` is a vector consisting of the three translation parameters, :math:`s`
+is the scaling factor and :math:`\mathbf{R}` is a rotation matrix. :math:`V^A` and
+:math:`V^B` are coordinate vectors, with :math:`V^A` being the input coordinate and
+:math:`V^B` is the output coordinate.
+
+The rotation matrix is composed of three rotation matrices, one for each axis:
+
+.. math::
+
+    \begin{align}
+        \mathbf{R}_X &= \begin{bmatrix} 1 & 0 & 0\\ 0 & \cos\phi & -\sin\phi\\ 0 & \sin\phi & \cos\phi \end{bmatrix}
+    \end{align}
+
+.. math::
+
+    \begin{align}
+        \mathbf{R}_Y &= \begin{bmatrix} \cos\theta & 0 & \sin\theta\\ 0 & 1 & 0\\ -\sin\theta & 0 & \cos\theta \end{bmatrix}
+    \end{align}
+
+.. math::
+
+    \begin{align}
+        \mathbf{R}_Z &= \begin{bmatrix} \cos\psi & -\sin\psi & 0\\ \sin\psi & \cos\psi & 0\\ 0 & 0 & 1 \end{bmatrix}
+    \end{align}
+
+The three rotation matrices can be combined in one:
+
+.. math::
+
+    \begin{align}
+        \mathbf{R} = \mathbf{R_X} \mathbf{R_Y} \mathbf{R_Y}
+    \end{align}
+
+For :math:`\mathbf{R}`, this yields:
+
+.. math::
+    :label: rot_exact
+
+    \begin{bmatrix}
+      \cos\theta \cos\psi &  -\cos\phi \sin\psi + \sin\phi \sin\theta \cos\psi &   \sin\phi \sin\psi + \cos\phi \sin\theta \cos\psi \\
+      \cos\theta\sin\psi &  \cos\phi \cos\psi + \sin\phi \sin\theta \sin\psi &  - \sin\phi \cos\psi + \cos\phi \sin\theta \sin\psi \\
+      -\sin\theta             &  \sin\phi \cos\theta                                          &   \cos\phi \cos\theta \\
+     \end{bmatrix}
+
+
+Using the small angle approxition the rotation matrix can be simplified to
+
+.. math::
+    :label: rot_approx
+
+    \begin{align} \mathbf{R} =
+        \begin{bmatrix}
+             1  & -R_z  &  R_y \\
+             Rz &  1    & -R_x \\
+            -Ry &  R_x  &  1   \\
+        \end{bmatrix}
+    \end{align}
+
+Which allow us to express the most common version of the Helmert transform,
+using the approximated rotation matrix:
+
+
+.. math::
+    :label: 7param
+
+    \begin{align}
+        \begin{bmatrix}
+            X \\
+            Y \\
+            Z \\
+        \end{bmatrix}^B =
+        \begin{bmatrix}
+            T_x \\
+            T_y \\
+            T_z \\
+        \end{bmatrix} +
+        \left(1 + s \times 10^{-6}\right)
+        \begin{bmatrix}
+             1  & -R_z  &  R_y \\
+             Rz &  1    & -R_x \\
+            -Ry &  R_x  &  1   \\
+        \end{bmatrix}
+        \begin{bmatrix}
+            X \\
+            Y \\
+            Z \\
+        \end{bmatrix}^A
+    \end{align}
+
+If the rotation matrix is transposed the transformation is effectively reversed.
+This is cause for some confusion since there is no correct way of defining the
+rotation matrix. Two conventions exists and they seem to be equally popular.
+In PROJ the rotation matrix can be transposed by adding the `transpose` flag in
+the transformation setup.
+
+Applying :eq:`propagation` we get the kinematic version of the approximated
+3D Helmert:
+
+.. math::
+    :label: 14param
+
+    \begin{align}
+        \begin{bmatrix}
+            X \\
+            Y \\
+            Z \\
+        \end{bmatrix}^B =
+        \begin{bmatrix}
+            \dot{T_x} \\
+            \dot{T_y} \\
+            \dot{T_z} \\
+        \end{bmatrix} +
+        \left(1 + \dot{s} \times 10^{-6}\right)
+        \begin{bmatrix}
+             1         & -\dot{R_z}  &  \dot{R_y} \\
+             \dot{R_z} &  1          & -\dot{R_x} \\
+            -\dot{R_y} &  \dot{R_x}  &  1      \\
+        \end{bmatrix}
+        \begin{bmatrix}
+            X \\
+            Y \\
+            Z \\
+        \end{bmatrix}^A
+    \end{align}
+
+
+
+
+The Helmert transformation can be applied without using the rotation parameters,
+in which case it becomes a simlpe translation of the origin of the coordinate
+system. When using the Helmert in this version equation :eq:`general-helmert`
+simplies to:
+
+.. math::
+    :label: 3param
+
+    \begin{align}
+        \begin{bmatrix}
+            X \\
+            Y \\
+            Z \\
+        \end{bmatrix}^B =
+        \begin{bmatrix}
+            T_x \\
+            T_y \\
+            T_z \\
+        \end{bmatrix} +
+        \begin{bmatrix}
+            X \\
+            Y \\
+            Z \\
+        \end{bmatrix}^A
+    \end{align}
+
+That after application of :eq:`propagation` has the following kinematic
+countpart:
+
+.. math::
+    :label: 6param
+
+    \begin{align}
+        \begin{bmatrix}
+            X \\
+            Y \\
+            Z \\
+        \end{bmatrix}^B =
+        \begin{bmatrix}
+            \dot{T_x} \\
+            \dot{T_y} \\
+            \dot{T_z} \\
+        \end{bmatrix} +
+        \begin{bmatrix}
+            X \\
+            Y \\
+            Z \\
+        \end{bmatrix}^A
+    \end{align}
diff --git a/docs/source/operations/transformations/hgridshift.rst b/docs/source/operations/transformations/hgridshift.rst
new file mode 100644
index 00000000..7f3c9895
--- /dev/null
+++ b/docs/source/operations/transformations/hgridshift.rst
@@ -0,0 +1,39 @@
+.. _hgridshift:
+
+================================================================================
+Horizontal grid shift
+================================================================================
+
+Change of horizontal datum by grid shift.
+
++-----------------+--------------------------------------------------------------------+
+| **Input type**  | Geodetic coordinates.                                              |
++-----------------+--------------------------------------------------------------------+
+| **Output type** | Geodetic coordinates.                                              |
++-----------------+--------------------------------------------------------------------+
+| **Options**                                                                          |
++-----------------+--------------------------------------------------------------------+
+| `+grids`        | Comma-separated list of grids to load.                             |
++-----------------+--------------------------------------------------------------------+
+
+The horizontal grid shift is done by offsetting the planar input coordinates by
+a specific amount determined by the loaded grids. The simplest use case of the
+horizontal grid shift is applying a single grid::
+
+    +hgridshift +grids=nzgr2kgrid0005.gsb
+
+
+More than one grid can be loaded at the same time, for instance in case the
+dataset needs to be transformed spans several countries. In this example grids
+of the continental US, Alaska and Canada is loaded at the same time::
+
+    +hgridshift +grids=@conus,@alaska,@ntv2_0.gsb,@ntv_can.dat
+
+The ``@`` in the above example states that the grid is optional, in case the grid
+is not found in the PROJ search path. The list of grids is prioritized so that
+grids in the start of the list takes presedence over the grids in the back of the
+list.
+
+PROJ supports CTable2, NTv1 and NTv2 files for horizontal grid corrections. Details
+about all three formats can be found in the GDAL documentation. GDAL reads and
+writes all three formats. Using GDAL for construction of new grids is recommended.
diff --git a/docs/source/operations/transformations/index.rst b/docs/source/operations/transformations/index.rst
new file mode 100644
index 00000000..5112ffa2
--- /dev/null
+++ b/docs/source/operations/transformations/index.rst
@@ -0,0 +1,17 @@
+.. _transformation_list:
+
+================================================================================
+Transformations
+================================================================================
+
+Transformations coordinate operation in which the two coordinate reference
+systems are based on different datums.
+
+.. toctree::
+   :maxdepth: 1
+
+   deformation
+   helmert
+   molodensky
+   hgridshift
+   vgridshift
diff --git a/docs/source/operations/transformations/molodensky.rst b/docs/source/operations/transformations/molodensky.rst
new file mode 100644
index 00000000..d4ff3e79
--- /dev/null
+++ b/docs/source/operations/transformations/molodensky.rst
@@ -0,0 +1,74 @@
+.. _molodensky:
+
+================================================================================
+Molodensky transform
+================================================================================
+
+The Molodensky transformation resembles a :ref:`Helmert` with zero
+rotations and a scale of unity, but converts directly from geodetic coordinates to
+geodetic coordinates, without the intermediate shifts to and from cartesian
+geocentric coordinates, associated with the Helmert transformation.
+The Molodensky transformation is simple to implement and to parameterize, requiring
+only the 3 shifts between the input and output frame, and the corresponding
+differences between the semimajor axes and flattening parameters of the reference
+ellipsoids. Due to its algorithmic simplicity, it was popular prior to the
+ubiquity of digital computers. Today, it is mostly interesting for historical
+reasons, but nevertheless indispensable due to the large amount of data that has
+already been transformed that way [EversKnudsen2017]_.
+
++---------------------+----------------------------------------------------------+
+| **Input type**      | Geodetic coordinates.                                    |
++---------------------+----------------------------------------------------------+
+| **Output type**     | Geodetic coordinates.                                    |
++---------------------+----------------------------------------------------------+
+| **Options**                                                                    |
++---------------------+----------------------------------------------------------+
+| `+da`               | Difference in semimajor axis of the defining ellipsoids. |
++---------------------+----------------------------------------------------------+
+| `+df`               | Difference in flattening of the defining ellipsoids.     |
++---------------------+----------------------------------------------------------+
+| `+dx`               | Offset of the X-axes of the defining ellipsoids.         |
++---------------------+----------------------------------------------------------+
+| `+dy`               | Offset of the Y-axes of the defining ellipsoids.         |
++---------------------+----------------------------------------------------------+
+| `+dz`               | Offset of the Z-axes of the defining ellipsoids.         |
++---------------------+----------------------------------------------------------+
+| `+ellps`            | Ellipsoid definition of source coordinates.              |
+|                     | Any specification can be used (e.g. `+a`, `+rf`, etc).   |
+|                     | If not specified, default ellipsoid is used.             |
++---------------------+----------------------------------------------------------+
+| `+abridged`         | Use the abridged version of the Molodensky transform.    |
+|                     | Optional.                                                |
++---------------------+----------------------------------------------------------+
+
+The Molodensky transform can be used to perform a datum shift from coordinate
+:math:`(\phi_1, \lambda_1, h_1)` to :math:`(\phi_2, \lambda_2, h_2)` where the two
+coordinates are referenced to different ellipsoids. This is based on three
+assumptions:
+
+    1. The cartesian axes, :math:`X, Y, Z`,  of the two ellipsoids are parallel.
+    2. The offset, :math:`\delta X, \delta Y, \delta Z`,  between the two ellipsoid
+       are known.
+    3. The characteristics of the two ellipsoids, expressed as the difference in
+       semimajor axis (:math:`\delta a`) and flattening (:math:`\delta f`),  are known.
+
+The Molodensky transform is mostly used for transforming between old systems
+dating back to the time before computers. The advantage of the Molodensky transform
+is that it is fairly simple to compute by hand. The ease of computation come at the
+cost of limited accuracy.
+
+A derivation of the mathematical formulas for the Molodensky transform can be found
+in [Deakin2004]_.
+
+
+Examples
+###############################################################################
+
+The abridged Molodensky::
+
+    proj=molodensky a=6378160 rf=298.25 da=-23 df=-8.120449e-8  dx=-134 dy=-48 dz=149 abridged
+
+The same transformation using the standard Molodensky::
+
+    proj=molodensky a=6378160 rf=298.25 da=-23 df=-8.120449e-8  dx=-134 dy=-48 dz=149
+
diff --git a/docs/source/operations/transformations/vgridshift.rst b/docs/source/operations/transformations/vgridshift.rst
new file mode 100644
index 00000000..e9f78310
--- /dev/null
+++ b/docs/source/operations/transformations/vgridshift.rst
@@ -0,0 +1,41 @@
+.. _vgridshift:
+
+================================================================================
+Vertical grid shift
+================================================================================
+
+Change Vertical datum change by grid shift
+
++-----------------+--------------------------------------------------------------------+
+| **Input type**  | Geodetic coordinates (horizontal), meters (vertical).              |
++-----------------+--------------------------------------------------------------------+
+| **Output type** | Geodetic coordinates (horizontal), meters (vertical).              |
++-----------------+--------------------------------------------------------------------+
+| **Options**                                                                          |
++-----------------+--------------------------------------------------------------------+
+| `+grids`        | Comma-separated list of grids to load.                             |
++-----------------+--------------------------------------------------------------------+
+
+The vertical grid shift is done by offsetting the vertical input coordinates by
+a specific amount determined by the loaded grids. The simplest use case of the
+horizontal grid shift is applying a single grid. Here we change the vertical
+reference from the ellipsoid to the global geoid model, EGM96::
+
+    +vgridshift +grids=egm96_16.gtx
+
+
+More than one grid can be loaded at the same time, for instance in the case where
+a better geoid model than the global is available for a certain area. Here the
+gridshift is set up so that the local DVR90 geoid model takes presedence over
+the global model::
+
+    +vgridshift +grids=@dvr90.gtx,egm96_16.gtx
+
+The ``@`` in the above example states that the grid is optional, in case the grid
+is not found in the PROJ search path. The list of grids is prioritized so that
+grids in the start of the list takes presedence over the grids in the back of the
+list.
+
+PROJ supports the GTX file format for vertical grid corrections. Details
+about all the format can be found in the GDAL documentation. GDAL both reads and
+writes the format. Using GDAL for construction of new grids is recommended.