Add geographic offset transformation method.

The Geographic offsets transformation adds an offset to the geographic longitude,
latitude coordinates, and an offset to the ellipsoidal height.
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.

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)

It is used for example to transform:
- from the old Greek geographic 2D CRS to the newer GGRS87 CRS
- from Tokyo + JSLD69 height to WGS 84
- from Baltic 1977 height to Black Sea height

It is also useful to document the implicit zero-offset transformation
we do in pipelines such as

+proj=pipeline +step +inv +proj=longlat +ellps=A
               +step +proj=longlat +ellps=B

that can be explicited as

+proj=pipeline +step +inv +proj=longlat +ellps=A
               +step +proj=geogoffset [+dlon=0 +dlat=0 +dh=0]
               +step +proj=longlat +ellps=B

9 files changed, 224 insertions, 2 deletions

diff --git a/docs/source/operations/transformations/geogoffset.rst b/docs/source/operations/transformations/geogoffset.rst
new file mode 100644
index 00000000..91c91c83
--- /dev/null
+++ b/docs/source/operations/transformations/geogoffset.rst
@@ -0,0 +1,68 @@
+.. _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.
+
+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..c2fe6baa 100644
--- a/docs/source/operations/transformations/index.rst
+++ b/docs/source/operations/transformations/index.rst
@@ -11,6 +11,7 @@ systems are based on different datums.
    :maxdepth: 1
 
    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..2320e812 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_geogoffset.c
 
 install-exec-local:
 	rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT)
diff --git a/src/PJ_geogoffset.c b/src/PJ_geogoffset.c
new file mode 100644
index 00000000..7f917fd7
--- /dev/null
+++ b/src/PJ_geogoffset.c
@@ -0,0 +1,111 @@
+/************************************************************************
+* 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.h"
+#include "projects.h"
+
+PROJ_HEAD(geogoffset, "Geographic Offset");
+
+struct pj_opaque_geogoffset {
+    double dlon;
+    double dlat;
+    double dh;
+};
+
+
+static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
+    struct pj_opaque_geogoffset *Q = (struct pj_opaque_geogoffset *) P->opaque;
+    /* offset coordinate */
+    obs.lpz.lam += Q->dlon;
+    obs.lpz.phi += Q->dlat;
+    obs.lpz.z   += Q->dh;
+    return obs;
+}
+
+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) {
+    struct pj_opaque_geogoffset *Q = (struct pj_opaque_geogoffset *) P->opaque;
+    /* offset coordinate */
+    obs.lpz.lam -= Q->dlon;
+    obs.lpz.phi -= Q->dlat;
+    obs.lpz.z   -= Q->dh;
+    return obs;
+}
+
+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;
+}
+
+
+/* Arcsecond to radians */
+#define ARCSEC_TO_RAD (DEG_TO_RAD / 3600.0)
+
+
+PJ *TRANSFORMATION(geogoffset,0 /* no need for ellipsoid */) {
+    struct pj_opaque_geogoffset *Q = pj_calloc(1, sizeof(struct pj_opaque_geogoffset));
+    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->dlon = pj_param(P->ctx, P->params, "ddlon").f * ARCSEC_TO_RAD;
+    Q->dlat = pj_param(P->ctx, P->params, "ddlat").f * ARCSEC_TO_RAD;
+    Q->dh = pj_param(P->ctx, P->params, "ddh").f;
+
+    return P;
+}
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 0779568a..f1e76ff3 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -72,6 +72,7 @@ SET(SRC_LIBPROJ_PJ
         PJ_fouc_s.c
         PJ_gall.c
         PJ_geoc.c
+        PJ_geogoffset.c
         PJ_geos.c
         PJ_gins8.c
         PJ_gnom.c
diff --git a/src/makefile.vc b/src/makefile.vc
index cf9d878c..21db1362 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_geogoffset.obj
 
 geodesic = geodesic.obj
 
diff --git a/src/pj_list.h b/src/pj_list.h
index dfcfe30a..6d94f4fd 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -46,6 +46,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..d3fe64f6 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -557,6 +557,38 @@ expect    0 -90
 
 -------------------------------------------------------------------------------
 
+
+-------------------------------------------------------------------------------
+Test for PJ_geogoffset
+-------------------------------------------------------------------------------
+-------------------------------------------------------------------------------
+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
+
 -------------------------------------------------------------------------------
 run the few gie-builtin tests, which are currently either awkward or impossible
 to express in the gie command set