Merge pull request #1138 from rouault/geogoffset

Add geographic offset transformation method.

10 files changed, 530 insertions, 2 deletions

diff --git a/docs/source/operations/transformations/affine.rst b/docs/source/operations/transformations/affine.rst
new file mode 100644
index 00000000..8b1681b3
--- /dev/null
+++ b/docs/source/operations/transformations/affine.rst
@@ -0,0 +1,123 @@
+.. _affine:
+
+================================================================================
+Affine transformation
+================================================================================
+
+.. versionadded:: 6.0.0
+
+The affine transformation applies translation and scaling/rotation terms on the
+x,y,z coordinates, and translation and scaling on the temporal cordinate.
+
++---------------------+----------------------------------------------------------+
+| **Alias**           | affine                                                   |
++---------------------+----------------------------------------------------------+
+| **Domain**          | 4D                                                       |
++---------------------+----------------------------------------------------------+
+| **Input type**      | XYZT                                                     |
++---------------------+----------------------------------------------------------+
+| **output type**     | XYZT                                                     |
++---------------------+----------------------------------------------------------+
+
+By default, the parameters are set for an identity transforms. The transformation
+is reversible unless the determinant of the sji matrix is 0, or `tscale` is 0
+
+
+Parameters
+################################################################################
+
+Optional
+-------------------------------------------------------------------------------
+
+.. option:: +xoff=<value>
+
+    Offset in X. Default value: 0
+
+.. option:: +yoff=<value>
+
+    Offset in Y. Default value: 0
+
+.. option:: +zoff=<value>
+
+    Offset in Z. Default value: 0
+
+.. option:: +toff=<value>
+
+    Offset in T. Default value: 0
+
+.. option:: +s11=<value>
+
+    Rotation/scaling term. Default value: 1
+
+.. option:: +s12=<value>
+
+    Rotation/scaling term. Default value: 0
+
+.. option:: +s13=<value>
+
+    Rotation/scaling term. Default value: 0
+
+.. option:: +s21=<value>
+
+    Rotation/scaling term. Default value: 0
+
+.. option:: +s22=<value>
+
+    Rotation/scaling term. Default value: 1
+
+.. option:: +s23=<value>
+
+    Rotation/scaling term. Default value: 0
+
+.. option:: +s31=<value>
+
+    Rotation/scaling term. Default value: 0
+
+.. option:: +s32=<value>
+
+    Rotation/scaling term. Default value: 0
+
+.. option:: +s33=<value>
+
+    Rotation/scaling term. Default value: 1
+
+.. option:: +tscale=<value>
+
+    Time scaling term. Default value: 1
+
+
+
+Mathematical description
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+.. math::
+    :label: formula
+
+    \begin{align}
+        \begin{bmatrix}
+            X \\
+            Y \\
+            Z \\
+            T \\
+        \end{bmatrix}^{dest} =
+        \begin{bmatrix}
+            xoff \\
+            yoff \\
+            zoff \\
+            toff \\
+        \end{bmatrix} +
+        \begin{bmatrix}
+             s11 & s12 & s13 & 0 \\
+             s21 & s22 & s23 & 0 \\
+             s31 & s32 & s33 & 0 \\
+             0   & 0   & 0   & tscale \\
+        \end{bmatrix}
+        \begin{bmatrix}
+            X \\
+            Y \\
+            Z \\
+            T \\
+        \end{bmatrix}^{source}
+    \end{align}
+
+
diff --git a/docs/source/operations/transformations/geogoffset.rst b/docs/source/operations/transformations/geogoffset.rst
new file mode 100644
index 00000000..f643485e
--- /dev/null
+++ b/docs/source/operations/transformations/geogoffset.rst
@@ -0,0 +1,70 @@
+.. _geogoffset:
+
+================================================================================
+Geographic offsets
+================================================================================
+
+.. versionadded:: 6.0.0
+
+The Geographic offsets transformation adds an offset to the geographic longitude,
+latitude coordinates, and an offset to the ellipsoidal height.
+
++---------------------+----------------------------------------------------------+
+| **Alias**           | geogoffset                                               |
++---------------------+----------------------------------------------------------+
+| **Domain**          | 3D                                                       |
++---------------------+----------------------------------------------------------+
+| **Input type**      | Geodetic coordinates (horizontal), meters (vertical)     |
++---------------------+----------------------------------------------------------+
+| **output type**     | Geodetic coordinates (horizontal), meters (vertical)     |
++---------------------+----------------------------------------------------------+
+
+This method is normally only used when low accuracy is tolerated. It is documented
+as coordinate operation method code 9619 (for geographic 2D) and 9660 (for
+geographic 3D) in the EPSG dataset (:cite:`EPSGGuidanceNumber7Part2`)
+
+It can also be used to implement the method Geographic2D with Height Offsets
+(code 9618) by noting that the input vertical component is a gravity-related
+height and the output vertical component is the ellispoid height (dh being
+the geoid undulation).
+
+It can also be used to implement the method Vertical offset (code 9616)
+
+The reverse transformation simply consists in subtracting the offsets.
+
+This method is a conveniency wrapper for the more general :ref:`affine`.
+
+Examples
+###############################################################################
+
+Geographic offset from the old Greek geographic 2D CRS to the newer GGRS87 CRS::
+
+    proj=geogoffset dlon=0.28 dlat=-5.86
+
+Conversion from Tokyo + JSLD69 height to WGS 84::
+
+    proj=geogoffset dlon=-13.97 dlat=7.94 dh=26.9
+
+Conversion from Baltic 1977 height to Black Sea height::
+
+    proj=geogoffset dh=0.4
+
+
+Parameters
+################################################################################
+
+Optional
+-------------------------------------------------------------------------------
+
+.. option:: +dlon=<value>
+
+    Offset in longitude, expressed in arc-second, to add.
+
+.. option:: +dlat=<value>
+
+    Offset in latitude, expressed in arc-second, to add.
+
+.. option:: +dh=<value>
+
+    Offset in height, expressed in meter, to add.
+
diff --git a/docs/source/operations/transformations/index.rst b/docs/source/operations/transformations/index.rst
index 5112ffa2..121215f7 100644
--- a/docs/source/operations/transformations/index.rst
+++ b/docs/source/operations/transformations/index.rst
@@ -10,7 +10,9 @@ systems are based on different datums.
 .. toctree::
    :maxdepth: 1
 
+   affine
    deformation
+   geogoffset
    helmert
    molodensky
    hgridshift
diff --git a/docs/source/references.bib b/docs/source/references.bib
index bd3089ee..0fb077ec 100644
--- a/docs/source/references.bib
+++ b/docs/source/references.bib
@@ -82,6 +82,14 @@
   Url                      = {http://www.calcofi.org/publications/calcofireports/v20/Vol_20_Eber___Hewitt.pdf}
 }
 
+@Manual{EPSGGuidanceNumber7Part2,
+  Title                    = {Geomatics Guidance Note 7, part 2 - Coordinate Conversions & Transformations including Formulas},
+  Institution              = {International Association For Oil And Gas Producers},
+  Year                     = {2018},
+
+  Url                      = {https://www.iogp.org/bookstore/product/coordinate-conversions-and-transformation-including-formulas/}
+}
+
 @Manual{Evenden2005,
   Title                    = {libproj4: A Comprehensive Library of Cartographic Projection Functions (Preliminary Draft)},
   Author                   = {G. I. Evenden},
diff --git a/src/Makefile.am b/src/Makefile.am
index 84846c1c..c1fbc12e 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -89,7 +89,7 @@ libproj_la_SOURCES = \
 	\
 	proj_4D_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c \
 	PJ_vgridshift.c PJ_hgridshift.c PJ_unitconvert.c PJ_molodensky.c \
-	PJ_deformation.c pj_internal.c PJ_axisswap.c
+	PJ_deformation.c pj_internal.c PJ_axisswap.c PJ_affine.c
 
 install-exec-local:
 	rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT)
diff --git a/src/PJ_affine.c b/src/PJ_affine.c
new file mode 100644
index 00000000..0d8b6374
--- /dev/null
+++ b/src/PJ_affine.c
@@ -0,0 +1,246 @@
+/************************************************************************
+* Copyright (c) 2018, Even Rouault <even.rouault at spatialys.com>
+*
+* Permission is hereby granted, free of charge, to any person obtaining a
+* copy of this software and associated documentation files (the "Software"),
+* to deal in the Software without restriction, including without limitation
+* the rights to use, copy, modify, merge, publish, distribute, sublicense,
+* and/or sell copies of the Software, and to permit persons to whom the
+* Software is furnished to do so, subject to the following conditions:
+*
+* The above copyright notice and this permission notice shall be included
+* in all copies or substantial portions of the Software.
+*
+* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+* DEALINGS IN THE SOFTWARE.
+*
+***********************************************************************/
+#define PJ_LIB__
+
+#include <errno.h>
+#include <math.h>
+
+#include "proj_internal.h"
+#include "proj.h"
+#include "projects.h"
+
+PROJ_HEAD(affine, "Affine transformation");
+PROJ_HEAD(geogoffset, "Geographic Offset");
+
+struct pj_affine_coeffs {
+    double s11;
+    double s12;
+    double s13;
+    double s21;
+    double s22;
+    double s23;
+    double s31;
+    double s32;
+    double s33;
+    double tscale;
+};
+
+struct pj_opaque_affine {
+    double xoff;
+    double yoff;
+    double zoff;
+    double toff;
+    struct pj_affine_coeffs forward;
+    struct pj_affine_coeffs reverse;
+};
+
+
+static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
+    PJ_COORD newObs;
+    const struct pj_opaque_affine *Q = (const struct pj_opaque_affine *) P->opaque;
+    const struct pj_affine_coeffs *C = &(Q->forward);
+    newObs.xyzt.x = Q->xoff + C->s11 * obs.xyzt.x + C->s12 * obs.xyzt.y + C->s13 * obs.xyzt.z;
+    newObs.xyzt.y = Q->yoff + C->s21 * obs.xyzt.x + C->s22 * obs.xyzt.y + C->s23 * obs.xyzt.z;
+    newObs.xyzt.z = Q->zoff + C->s31 * obs.xyzt.x + C->s32 * obs.xyzt.y + C->s33 * obs.xyzt.z;
+    newObs.xyzt.t = Q->toff + C->tscale * obs.xyzt.t;
+    return newObs;
+}
+
+static XYZ forward_3d(LPZ lpz, PJ *P) {
+    PJ_COORD point = {{0,0,0,0}};
+    point.lpz = lpz;
+    return forward_4d(point, P).xyz;
+}
+
+
+static XY forward_2d(LP lp, PJ *P) {
+    PJ_COORD point = {{0,0,0,0}};
+    point.lp = lp;
+    return forward_4d(point, P).xy;
+}
+
+
+static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
+    PJ_COORD newObs;
+    const struct pj_opaque_affine *Q = (const struct pj_opaque_affine *) P->opaque;
+    const struct pj_affine_coeffs *C = &(Q->reverse);
+    obs.xyzt.x -= Q->xoff;
+    obs.xyzt.y -= Q->yoff;
+    obs.xyzt.z -= Q->zoff;
+    newObs.xyzt.x = C->s11 * obs.xyzt.x + C->s12 * obs.xyzt.y + C->s13 * obs.xyzt.z;
+    newObs.xyzt.y = C->s21 * obs.xyzt.x + C->s22 * obs.xyzt.y + C->s23 * obs.xyzt.z;
+    newObs.xyzt.z = C->s31 * obs.xyzt.x + C->s32 * obs.xyzt.y + C->s33 * obs.xyzt.z;
+    newObs.xyzt.t = C->tscale * (obs.xyzt.t - Q->toff);
+    return newObs;
+}
+
+static LPZ reverse_3d(XYZ xyz, PJ *P) {
+    PJ_COORD point = {{0,0,0,0}};
+    point.xyz = xyz;
+    return reverse_4d(point, P).lpz;
+}
+
+static LP reverse_2d(XY xy, PJ *P) {
+    PJ_COORD point = {{0,0,0,0}};
+    point.xy = xy;
+    return reverse_4d(point, P).lp;
+}
+
+static struct pj_opaque_affine * initQ() {
+    struct pj_opaque_affine *Q = pj_calloc(1, sizeof(struct pj_opaque_affine));
+    if (0==Q)
+        return 0;
+
+    /* default values */
+    Q->forward.s11 = 1.0;
+    Q->forward.s22 = 1.0;
+    Q->forward.s33 = 1.0;
+    Q->forward.tscale = 1.0;
+
+    Q->reverse.s11 = 1.0;
+    Q->reverse.s22 = 1.0;
+    Q->reverse.s33 = 1.0;
+    Q->reverse.tscale = 1.0;
+
+    return Q;
+}
+
+static void computeReverseParameters(PJ* P)
+{
+    struct pj_opaque_affine *Q = (struct pj_opaque_affine *) P->opaque;
+
+    /* cf https://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_3_%C3%97_3_matrices */
+    const double a = Q->forward.s11;
+    const double b = Q->forward.s12;
+    const double c = Q->forward.s13;
+    const double d = Q->forward.s21;
+    const double e = Q->forward.s22;
+    const double f = Q->forward.s23;
+    const double g = Q->forward.s31;
+    const double h = Q->forward.s32;
+    const double i = Q->forward.s33;
+    const double A = e * i - f * h;
+    const double B = -(d * i - f * g);
+    const double C = (d * h - e * g);
+    const double D = -(b * i - c * h);
+    const double E = (a * i - c * g);
+    const double F = -(a * h - b * g);
+    const double G = b * f - c * e;
+    const double H = -(a * f - c * d);
+    const double I = a * e - b * d;
+    const double det = a * A + b * B + c * C;
+    if( det == 0.0 || Q->forward.tscale == 0.0 ) {
+        if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) {
+            proj_log_debug(P, "Affine: matrix non invertible");
+        }
+        P->inv4d = NULL;
+        P->inv3d = NULL;
+        P->inv = NULL;
+    } else {
+        Q->reverse.s11 = A / det;
+        Q->reverse.s12 = D / det;
+        Q->reverse.s13 = G / det;
+        Q->reverse.s21 = B / det;
+        Q->reverse.s22 = E / det;
+        Q->reverse.s23 = H / det;
+        Q->reverse.s31 = C / det;
+        Q->reverse.s32 = F / det;
+        Q->reverse.s33 = I / det;
+        Q->reverse.tscale = 1.0 / Q->forward.tscale;
+    }
+}
+
+PJ *TRANSFORMATION(affine,0 /* no need for ellipsoid */) {
+    struct pj_opaque_affine *Q = initQ();
+    if (0==Q)
+        return pj_default_destructor(P, ENOMEM);
+    P->opaque = (void *) Q;
+
+    P->fwd4d = forward_4d;
+    P->inv4d = reverse_4d;
+    P->fwd3d  = forward_3d;
+    P->inv3d  = reverse_3d;
+    P->fwd    = forward_2d;
+    P->inv    = reverse_2d;
+
+    P->left   = PJ_IO_UNITS_WHATEVER;
+    P->right  = PJ_IO_UNITS_WHATEVER;
+
+    /* read args */
+    Q->xoff = pj_param(P->ctx, P->params, "dxoff").f;
+    Q->yoff = pj_param(P->ctx, P->params, "dyoff").f;
+    Q->zoff = pj_param(P->ctx, P->params, "dzoff").f;
+    Q->toff = pj_param(P->ctx, P->params, "dtoff").f;
+
+    if(pj_param (P->ctx, P->params, "ts11").i) {
+        Q->forward.s11 = pj_param(P->ctx, P->params, "ds11").f;
+    }
+    Q->forward.s12 = pj_param(P->ctx, P->params, "ds12").f;
+    Q->forward.s13 = pj_param(P->ctx, P->params, "ds13").f;
+    Q->forward.s21 = pj_param(P->ctx, P->params, "ds21").f;
+    if(pj_param (P->ctx, P->params, "ts22").i) {
+        Q->forward.s22 = pj_param(P->ctx, P->params, "ds22").f;
+    }
+    Q->forward.s23 = pj_param(P->ctx, P->params, "ds23").f;
+    Q->forward.s31 = pj_param(P->ctx, P->params, "ds31").f;
+    Q->forward.s32 = pj_param(P->ctx, P->params, "ds32").f;
+    if(pj_param (P->ctx, P->params, "ts33").i) {
+        Q->forward.s33 = pj_param(P->ctx, P->params, "ds33").f;
+    }
+    if(pj_param (P->ctx, P->params, "ttscale").i) {
+        Q->forward.tscale = pj_param(P->ctx, P->params, "dtscale").f;
+    }
+
+    computeReverseParameters(P);
+
+    return P;
+}
+
+
+/* Arcsecond to radians */
+#define ARCSEC_TO_RAD (DEG_TO_RAD / 3600.0)
+
+
+PJ *TRANSFORMATION(geogoffset,0 /* no need for ellipsoid */) {
+    struct pj_opaque_affine *Q = initQ();
+    if (0==Q)
+        return pj_default_destructor(P, ENOMEM);
+    P->opaque = (void *) Q;
+
+    P->fwd4d = forward_4d;
+    P->inv4d = reverse_4d;
+    P->fwd3d  = forward_3d;
+    P->inv3d  = reverse_3d;
+    P->fwd    = forward_2d;
+    P->inv    = reverse_2d;
+
+    P->left   = PJ_IO_UNITS_ANGULAR;
+    P->right  = PJ_IO_UNITS_ANGULAR;
+
+    /* read args */
+    Q->xoff = pj_param(P->ctx, P->params, "ddlon").f * ARCSEC_TO_RAD;
+    Q->yoff = pj_param(P->ctx, P->params, "ddlat").f * ARCSEC_TO_RAD;
+    Q->zoff = pj_param(P->ctx, P->params, "ddh").f;
+
+    return P;
+}
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 0779568a..a359369e 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -40,6 +40,7 @@ SET(SRC_LIBPROJ_PJ
         nad_init.c
         PJ_aea.c
         PJ_aeqd.c
+        PJ_affine.c
         PJ_airy.c
         PJ_aitoff.c
         PJ_august.c
diff --git a/src/makefile.vc b/src/makefile.vc
index cf9d878c..73c68c0f 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -61,7 +61,7 @@ support = \
 pipeline = \
 	proj_4D_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj \
 	PJ_vgridshift.obj PJ_hgridshift.obj PJ_unitconvert.obj PJ_molodensky.obj \
-	PJ_deformation.obj PJ_axisswap.obj
+	PJ_deformation.obj PJ_axisswap.obj PJ_affine.obj
 
 geodesic = geodesic.obj
 
diff --git a/src/pj_list.h b/src/pj_list.h
index dfcfe30a..72a9af7f 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -7,6 +7,7 @@ static const char PJ_LIST_H_ID[] = "@(#)pj_list.h	4.5	95/08/09	GIE	REL";
 */
 PROJ_HEAD(aea, "Albers Equal Area")
 PROJ_HEAD(aeqd, "Azimuthal Equidistant")
+PROJ_HEAD(affine, "Affine transformation")
 PROJ_HEAD(airy, "Airy")
 PROJ_HEAD(aitoff, "Aitoff")
 PROJ_HEAD(alsk, "Mod. Stererographics of Alaska")
@@ -46,6 +47,7 @@ PROJ_HEAD(fouc_s, "Foucaut Sinusoidal")
 PROJ_HEAD(gall, "Gall (Gall Stereographic)")
 PROJ_HEAD(geoc, "Geocentric Latitude")
 PROJ_HEAD(geocent, "Geocentric")
+PROJ_HEAD(geogoffset, "Geographic Offset")
 PROJ_HEAD(geos, "Geostationary Satellite View")
 PROJ_HEAD(gins8, "Ginsburg VIII (TsNIIGAiK)")
 PROJ_HEAD(gn_sinu, "General Sinusoidal Series")
diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
index 96212cba..10feea70 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -557,6 +557,82 @@ expect    0 -90
 
 -------------------------------------------------------------------------------
 
+
+-------------------------------------------------------------------------------
+Test for PJ_affine
+-------------------------------------------------------------------------------
+-------------------------------------------------------------------------------
+operation +proj=geogoffset
+-------------------------------------------------------------------------------
+direction forward
+tolerance 1mm
+
+accept    10 20
+expect    10 20
+roundtrip  1
+
+-------------------------------------------------------------------------------
+operation +proj=geogoffset +dlon=3600 +dlat=-3600 +dh=3
+-------------------------------------------------------------------------------
+direction forward
+tolerance 1mm
+
+accept    10 20
+expect    11 19
+roundtrip  1
+
+accept    10 20 30
+expect    11 19 33
+roundtrip  1
+
+accept    10 20 30 40
+expect    11 19 33 40
+roundtrip  1
+
+-------------------------------------------------------------------------------
+operation +proj=affine
+-------------------------------------------------------------------------------
+direction forward
+tolerance 1mm
+
+accept    10 20 30 40
+expect    10 20 30 40
+roundtrip  1
+
+-------------------------------------------------------------------------------
+operation +proj=affine +xoff=1 +yoff=2 +zoff=3 +toff=4 +s11=11 +s12=12 +s13=13 +s21=21 +s22=22 +s23=23 +s31=-31 +s32=32 +s33=33 +tscale=34
+-------------------------------------------------------------------------------
+direction forward
+tolerance 1mm
+
+accept    2 49 10 100
+expect    741.0000      1352.0000     1839.0000     3404.0000
+roundtrip  1
+
+accept    2 49 10
+expect    741.0000      1352.0000     1839.0000
+roundtrip  1
+
+accept    2 49
+expect    611.0000      1122.0000
+roundtrip  1
+
+-------------------------------------------------------------------------------
+# Non invertible
+operation +proj=affine +s11=0 +s22=0 +s23=0
+-------------------------------------------------------------------------------
+direction reverse
+accept    0 0 0 0
+expect    failure
+
+-------------------------------------------------------------------------------
+# Non invertible
+operation +proj=affine +tscale=0
+-------------------------------------------------------------------------------
+direction reverse
+accept    0 0 0 0
+expect    failure
+
 -------------------------------------------------------------------------------
 run the few gie-builtin tests, which are currently either awkward or impossible
 to express in the gie command set