[BREAKING] Hermert: add +convention=position_vector/coordinate_frame, forbids +transpose (fixes #1091)

As identified in #1091, Helmert implementation in PROJ 5.0 and 5.1 is confusing.
It happens that by default it used the coordinate_frame convention, contrary to
the position_vector convention used traditionaly for +towgs84. The documentation
of Helmert was also wrongly specifying that the default convention was
position_vector.

This commit:
- bans the confusing +transpose parameter
- removes the concept of a default convention, since in practice both are
  equally found, and requires +convention as soon as a rotational term parameter
  is present.
  For translation only, convention is ignored and optional, as having no effect.
- fixes all the identified uses of proj=helmert in code, doc and tests

This is obviously a breaking change:
- users will have to adapt their pipeline expressions
- in particular, init files that would use helmert must be adapted

However, as designed, the break will be explicit, and not silent.

4 files changed, 61 insertions, 33 deletions

diff --git a/docs/source/operations/transformations/deformation.rst b/docs/source/operations/transformations/deformation.rst
index 79d9b7b8..87b08453 100644
--- a/docs/source/operations/transformations/deformation.rst
+++ b/docs/source/operations/transformations/deformation.rst
@@ -51,7 +51,7 @@ to the Danish realisation of ETRS89 is in PROJ described as::
             # ITRF2008@t_obs -> ITRF2000@t_obs
     step    init = ITRF2008:ITRF2000
             # ITRF2000@t_obs -> ETRF2000@t_obs
-    step    proj=helmert t_epoch = 2000.0 transpose
+    step    proj=helmert t_epoch = 2000.0 convention=position_vector
             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
@@ -60,7 +60,7 @@ to the Danish realisation of ETRS89 is in PROJ described as::
             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
+    step    proj=helmert convention=position_vector s = -0.009420e
             x = 0.03863 rx = 0.00617753
             y = 0.147   ry = 5.064e-05
             z = 0.02776 rz = 4.729e-05
diff --git a/docs/source/operations/transformations/helmert.rst b/docs/source/operations/transformations/helmert.rst
index 69776285..6adeaa17 100644
--- a/docs/source/operations/transformations/helmert.rst
+++ b/docs/source/operations/transformations/helmert.rst
@@ -53,25 +53,25 @@ Transforming coordinates from NAD72 to NAD83 using the 4 parameter 2D Helmert:
 
 ::
 
-    proj=helmert x=-9597.3572 y=.6112 s=0.304794780637 theta=-1.244048
+    proj=helmert convention=coordinate_frame x=-9597.3572 y=.6112 s=0.304794780637 theta=-1.244048
 
 Simplified transformations from ITRF2008/IGS08 to ETRS89 using 7 parameters:
 
 ::
 
-    proj=helmert x=0.67678    y=0.65495   z=-0.52827
+    proj=helmert convention=coordinate_frame 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
+    proj=helmert convention=position_vector 
          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
-         t_epoch=1988.0 t_obs=2017.0    transpose
+         t_epoch=1988.0 t_obs=2017.0
 
 Parameters
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -81,6 +81,36 @@ Parameters
     All parameters are optional but at least one should be used, otherwise the
     operation will return the coordinates unchanged.
 
+.. option:: +convention=coordinate_frame/position_vector
+
+    .. versionadded:: 5.2.0
+
+    Indicates the convention to express the rotational terms when a 3D-Helmert /
+    7-parameter more transform is involved. As soon as a rotational parameter
+    is specified (one of ``rx``, ``ry``, ``rz``, ``drx``, ``dry``, ``drz``),
+    ``convention`` is required.
+
+    The two conventions are equally popular and a frequent source of confusion.
+    The coordinate frame convention is also described as an clockwise
+    rotation of the coordinate frame. It corresponds to EPSG method code
+    1032 (in the geocentric domain) or 9607 (in the geographic domain)
+    The position vector convention is also described as an anticlockwise
+    (counter-clockwise) rotation of the coordinate frame.
+    It corresponds to as EPSG method code 1033 (in the geocentric domain) or
+    9606 (in the geographic domain).
+
+    This parameter is ignored when only a 3-parameter
+    (translation terms only: ``x``, ``y``, ``z``) , 4-parameter (3-parameter
+    and ``theta``) or 6-parameter (3-parameter and their derivative terms)
+    is used.
+
+    The result obtained with parameters specified in a given convention
+    can be obtained in the other convention by negating the rotational parameters
+    (``rx``, ``ry``, ``rz``, ``drx``, ``dry``, ``drz``)
+
+    .. note:: This parameter obsoletes ``transpose`` which was present in
+              PROJ 5.0 and 5.1, and is forbidden starting with PROJ 5.2
+
 .. option:: +x=<value>
 
     Translation of the x-axis given in meters.
@@ -161,13 +191,6 @@ Parameters
 
     See :eq:`rot_exact`
 
-.. option:: +transpose
-
-    Transpose rotation matrix and follow the **Cordinate Frame** rotation
-    convention. If :option:`+transpose` is not added the **Position Vector**
-    rotation convention is used.
-
-
 
 Mathematical description
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -278,24 +301,30 @@ is the scaling factor and :math:`\mathbf{R}` is a rotation matrix. :math:`V^A` a
 :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:
+In the *Position Vector* convention, we define :math:`R_x = radians \left( rx \right)`,
+:math:`R_z = radians \left( ry \right)` and :math:`R_z = radians \left( rz \right)`
+
+In the *Coordinate Frame* convention, :math:`R_x = - radians \left( rx \right)`,
+:math:`R_z = - radians \left( ry \right)` and :math:`R_z = - radians \left( rz \right)`
+
+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}
+        \mathbf{R}_X &= \begin{bmatrix} 1 & 0 & 0\\ 0 & \cos R_x  & -\sin R_x \\ 0 & \sin R_x  & \cos R_x  \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}
+        \mathbf{R}_Y &= \begin{bmatrix} \cos R_y & 0 & \sin R_y\\ 0 & 1 & 0\\ -\sin R_y & 0 & \cos R_y \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}
+        \mathbf{R}_Z &= \begin{bmatrix} \cos R_z  & -\sin R_z  & 0\\ \sin R_z  & \cos R_z  & 0\\ 0 & 0 & 1 \end{bmatrix}
     \end{align}
 
 The three rotation matrices can be combined in one:
@@ -313,9 +342,11 @@ For :math:`\mathbf{R}`, this yields:
     :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 \\
+      \cos R_y \cos R_z  &  -\cos R_x  \sin R_z  +       &   \sin R_x  \sin R_z  +   \\
+                         &  \sin R_x  \sin R_y \cos R_z  &   \cos R_x  \sin R_y \cos R_z  \\
+      \cos R_y\sin R_z   &  \cos R_x  \cos R_z  +        &  - \sin R_x  \cos R_z +  \\
+                         &  \sin R_x  \sin R_y \sin R_z  &   \cos R_x  \sin R_y \sin R_z  \\
+      -\sin R_y          &  \sin R_x  \cos R_y            &   \cos R_x  \cos R_y \\
      \end{bmatrix}
 
 
@@ -363,12 +394,10 @@ using the approximated rotation matrix:
         \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. PROJ
-uses the **Position Vector** rotation convention. The rotation matrix can be transposed by
-adding the :option:`+transpose` flag in the transformation setup which makes PROJ
-follow the **Coordinate Frame** rotation convention.
+If the rotation matrix is transposed, or the sign of the rotation terms negated,
+the rotational part of the transformation is effectively reversed.
+This is what happens when switching between the 2 conventions ``position_vector``
+and ``coordinate_frame``
 
 Applying :eq:`propagation` we get the kinematic version of the approximated
 3D Helmert:
diff --git a/docs/source/resource_files.rst b/docs/source/resource_files.rst
index 10cbcb7b..c0321b18 100644
--- a/docs/source/resource_files.rst
+++ b/docs/source/resource_files.rst
@@ -412,7 +412,7 @@ which then expands to
 
     +proj=helmert +x=-0.0001 +y=0.0008 +z=0.0058 +s=-0.0004
     +dx=0.0002 +dy=-0.0001 +dz=0.0018 +ds=-0.000008
-    +t_epoch=2000.0 +transpose
+    +t_epoch=2000.0 +convention=position_vector
     +t_obs=2010.5
 
 Below is a list of the init files that are packaged with PROJ.
diff --git a/docs/source/usage/transformation.rst b/docs/source/usage/transformation.rst
index acec0bb5..b4ed7031 100644
--- a/docs/source/usage/transformation.rst
+++ b/docs/source/usage/transformation.rst
@@ -75,7 +75,7 @@ the introduction). In PROJ it can be implemented as
 
     proj=pipeline
     step proj=cart ellps=intl
-    step proj=helmert
+    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
@@ -101,7 +101,7 @@ deprecated system with decimeter level tensions.
     step init=./s45b.pol:s45b_tc32
     step proj=utm inv ellps=intl zone=32
     step proj=cart ellps=intl
-    step proj=helmert
+    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
@@ -116,16 +116,15 @@ data, but the 14-parameter Helmert transform expects temporal units in decimalye
 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 :cite:`Altamimi2002`,
-except the epoch which is the epoch of the transformation. The default setting is to
-use “coordinate frame” convention of the Helmert transform, but “position vector”
-convention can also be used. The last step in the pipeline is converting the
+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
+    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