deformation: add support for +grids= for GeoTIFF grids

This option is to load grid(s) that have both the horizontal and
vertical velocities in the same file.

Can be tested with the following grid
https://github.com/rouault/sample_proj_gtiff_grids/blob/master/nkgrf03vel_realigned.tif
converted from the original .ct2 and .gtx with
```
gdal_translate nkgrf03vel_realigned.vrt nkgrf03vel_realigned.tif -co COMPRESS=DEFLATE -co PREDICTOR=3 -co BLOCKYSIZE=241 -co INTERLEAVE=BAND
```

where nkgrf03vel_realigned.vrt is
```
<VRTDataset rasterXSize="223" rasterYSize="241">
  <SRS dataAxisToSRSAxisMapping="2,1">GEOGCRS["Unknown based on GRS80",
    DATUM["Unknown based on GRS80",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]]]
</SRS>
  <GeoTransform>  2.9166666666666670e+00,  1.6666666666666666e-01,  0.0000000000000000e+00,  7.3041666666666686e+01,  0.0000000000000000e+00, -8.3333333333333329e-02</GeoTransform>
  <Metadata>
    <MDI key="DESCRIPTION"></MDI>
    <MDI key="area_of_use">Nordic and Baltic countries</MDI>
    <MDI key="AREA_OR_POINT">Point</MDI>
    <MDI key="TIFFTAG_COPYRIGHT">The Nordic Geodetic Commission. Creative Commons Attribution 4.0 https://creativecommons.org/licenses/by/4.0/</MDI>
    <MDI key="TIFFTAG_DATETIME">2019:12:30 00:00:00</MDI>
    <MDI key="TIFFTAG_IMAGEDESCRIPTION">Deformation model covering the Nordic and Baltic countries. Used in transformations between global reference frames and the local realisations of ETRS89 in the Nordic and Baltic countries</MDI>
    <MDI key="TYPE">VELOCITY</MDI>
  </Metadata>
  <VRTRasterBand dataType="Float32" band="1">
    <Description>east_velocity</Description>
    <UnitType>mm/year</UnitType>
    <SimpleSource>
      <SourceFilename relativeToVRT="0">/home/even/proj/proj-datumgrid/europe/nkgrf03vel_realigned_xy.ct2</SourceFilename>
      <SourceBand>2</SourceBand> <!-- GDAL exposes the first physical component of the file (longitude offset normally, here east velocity) as the second band -->
      <SrcRect xOff="0" yOff="0" xSize="223" ySize="241" />
      <DstRect xOff="0" yOff="0" xSize="223" ySize="241" />
    </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="2">
    <Description>north_velocity</Description>
    <UnitType>mm/year</UnitType>
    <SimpleSource>
      <SourceFilename relativeToVRT="0">/home/even/proj/proj-datumgrid/europe/nkgrf03vel_realigned_xy.ct2</SourceFilename>
      <SourceBand>1</SourceBand> <!-- and the second physical component (latitude offset normally, here north velocity) as the first band -->
      <SrcRect xOff="0" yOff="0" xSize="223" ySize="241" />
      <DstRect xOff="0" yOff="0" xSize="223" ySize="241" />
    </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="3">
    <Description>up_velocity</Description>
    <UnitType>mm/year</UnitType>
    <SimpleSource>
      <SourceFilename relativeToVRT="0">/home/even/proj/proj-datumgrid/europe/nkgrf03vel_realigned_z.gtx</SourceFilename>
      <SourceBand>1</SourceBand>
      <SrcRect xOff="0" yOff="0" xSize="223" ySize="241" />
      <DstRect xOff="0" yOff="0" xSize="223" ySize="241" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>
```

7 files changed, 236 insertions, 25 deletions

diff --git a/data/Makefile.am b/data/Makefile.am
index 09fa1989..4bba2649 100644
--- a/data/Makefile.am
+++ b/data/Makefile.am
@@ -82,6 +82,9 @@ EXTRA_DIST = proj.ini GL27 nad.lst nad27 nad83 \
 		tests/subset_of_gr3df97a.tif \
 		tests/egm96_15_uncompressed_truncated.tif \
 		tests/test_vgrid_single_strip_truncated.tif \
+		tests/nkgrf03vel_realigned_extract.tif \
+		tests/nkgrf03vel_realigned_xy_extract.ct2 \
+		tests/nkgrf03vel_realigned_z_extract.gtx \
 		null \
 		generate_all_sql_in.cmake sql_filelist.cmake \
 		$(SQL_ORDERED_LIST)
diff --git a/data/tests/nkgrf03vel_realigned_extract.tif b/data/tests/nkgrf03vel_realigned_extract.tif
new file mode 100644
index 00000000..70e79e93
--- /dev/null
+++ b/data/tests/nkgrf03vel_realigned_extract.tif
diff --git a/data/tests/nkgrf03vel_realigned_xy_extract.ct2 b/data/tests/nkgrf03vel_realigned_xy_extract.ct2
new file mode 100644
index 00000000..89232b9f
--- /dev/null
+++ b/data/tests/nkgrf03vel_realigned_xy_extract.ct2
diff --git a/data/tests/nkgrf03vel_realigned_z_extract.gtx b/data/tests/nkgrf03vel_realigned_z_extract.gtx
new file mode 100644
index 00000000..5ea8aac7
--- /dev/null
+++ b/data/tests/nkgrf03vel_realigned_z_extract.gtx
diff --git a/docs/source/operations/transformations/deformation.rst b/docs/source/operations/transformations/deformation.rst
index 3a9d025c..b45c5e2b 100644
--- a/docs/source/operations/transformations/deformation.rst
+++ b/docs/source/operations/transformations/deformation.rst
@@ -35,6 +35,9 @@ GTX format. Both grids are expected to contain grid-values in units of
 mm/year. GDAL both reads and writes both file formats. Using GDAL for
 construction of new grids is recommended.
 
+Starting with PROJ 7.0, use of a GeoTIFF format is recommended to store both
+the horizontal and vertical velocities.
+
 Example
 -------------------------------------------------------------------------------
 
@@ -93,7 +96,7 @@ Parameters
     grid is considered optional and PROJ will the not complain if the grid is
     not available.
 
-    Grids for the horizontla component of a deformation model is expected to be
+    Grids for the horizontal component of a deformation model is expected to be
     in CTable2 format.
 
 .. option:: +z_grids=<list>
@@ -105,6 +108,34 @@ Parameters
     Grids for the vertical component of a deformation model is expected to be
     in either GTX format.
 
+.. option:: +grids=<list>
+
+    .. versionadded:: 7.0.0
+
+    Comma-separated list of grids to load. If a grid is prefixed by an `@` the
+    grid is considered optional and PROJ will the not complain if the grid is
+    not available.
+
+    Grids should be in GeoTIFF format with the first 3 components being
+    respectively the easting, northing and up velocities in mm/year.
+    Setting the Description and Unit Type GDAL band metadata items is strongly
+    recommended, so that gdalinfo reports:
+
+    ::
+
+        Band 1 Block=... Type=Float32, ColorInterp=Gray
+            Description = east_velocity
+            Unit Type: mm/year
+        Band 2 Block=... Type=Float32, ColorInterp=Undefined
+            Description = north_velocity
+            Unit Type: mm/year
+        Band 3 Block=... Type=Float32, ColorInterp=Undefined
+            Description = up_velocity
+            Unit Type: mm/year
+
+    .. note:: :option:`+grids` is mutually exclusive with :option:`+xy_grids`
+              and :option:`+z_grids`
+
 .. option:: +t_epoch=<value>
 
     Central epoch of transformation given in decimalyears. Will be used in
diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp
index f5468775..031e9eba 100644
--- a/src/transformations/deformation.cpp
+++ b/src/transformations/deformation.cpp
@@ -58,6 +58,8 @@ grid-values in units of mm/year in ENU-space.
 #include <math.h>
 #include "grids.hpp"
 
+#include <algorithm>
+
 PROJ_HEAD(deformation, "Kinematic grid shift");
 
 #define TOL 1e-8
@@ -70,11 +72,124 @@ struct deformationData {
     double dt = 0;
     double t_epoch = 0;
     PJ *cart = nullptr;
+    ListOfGenericGrids grids{};
     ListOfHGrids hgrids{};
     ListOfVGrids vgrids{};
 };
 } // anonymous namespace
 
+// ---------------------------------------------------------------------------
+
+static const GenericShiftGrid* findGrid(const ListOfGenericGrids& grids,
+                                        const PJ_LP& input)
+{
+    for( const auto& gridset: grids )
+    {
+        auto grid = gridset->gridAt(input.lam, input.phi);
+        if( grid )
+            return grid;
+    }
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+static bool get_grid_values(PJ* P,
+                            deformationData* Q,
+                            const PJ_LP& lp,
+                            double& vx,
+                            double& vy,
+                            double& vz)
+{
+    auto grid = findGrid(Q->grids, lp);
+    if( !grid ) {
+        return false;
+    }
+    if( grid->isNullGrid() ) {
+        vx = 0;
+        vy = 0;
+        vz = 0;
+        return true;
+    }
+    const auto samplesPerPixel = grid->samplesPerPixel();
+    if( samplesPerPixel < 3 ) {
+        proj_log_error(P, "deformation: grid has not enough samples");
+        return false;
+    }
+    int sampleE = 0;
+    int sampleN = 1;
+    int sampleU = 2;
+    for( int i = 0; i < samplesPerPixel; i++ )
+    {
+        const auto desc = grid->description(i);
+        if( desc == "east_velocity") {
+            sampleE = i;
+        } else if( desc == "north_velocity") {
+            sampleN = i;
+        } else if( desc == "up_velocity") {
+            sampleU = i;
+        }
+    }
+    const auto unit = grid->unit(sampleE);
+    if( !unit.empty() && unit != "mm/year" ) {
+        proj_log_error(P, "deformation: Only unit=mm/year currently handled");
+        return false;
+    }
+
+    const auto& extent = grid->extentAndRes();
+    double grid_x = (lp.lam - extent.westLon) / extent.resLon;
+    double grid_y = (lp.phi - extent.southLat) / extent.resLat;
+    int ix = static_cast<int>(grid_x);
+    int iy = static_cast<int>(grid_y);
+    int ix2 = std::min(ix + 1, grid->width() - 1);
+    int iy2 = std::min(iy + 1, grid->height() - 1);
+
+    float dx1, dy1, dz1;
+    if( !grid->valueAt(ix, iy, sampleE, dx1) ||
+        !grid->valueAt(ix, iy, sampleN, dy1) ||
+        !grid->valueAt(ix, iy, sampleU, dz1) ) {
+        return false;
+    }
+
+    float dx2, dy2, dz2;
+    if( !grid->valueAt(ix2, iy, sampleE, dx2) ||
+        !grid->valueAt(ix2, iy, sampleN, dy2) ||
+        !grid->valueAt(ix2, iy, sampleU, dz2) ) {
+        return false;
+    }
+
+    float dx3, dy3, dz3;
+    if( !grid->valueAt(ix, iy2, sampleE, dx3) ||
+        !grid->valueAt(ix, iy2, sampleN, dy3) ||
+        !grid->valueAt(ix, iy2, sampleU, dz3) ) {
+        return false;
+    }
+
+    float dx4, dy4, dz4;
+    if( !grid->valueAt(ix2, iy2, sampleE, dx4) ||
+        !grid->valueAt(ix2, iy2, sampleN, dy4) ||
+        !grid->valueAt(ix2, iy2, sampleU, dz4) ) {
+        return false;
+    }
+
+    double frct_lam = grid_x - ix;
+    double frct_phi = grid_y - iy;
+    double m10 = frct_lam;
+    double m11 = m10;
+    double m01 = 1. - frct_lam;
+    double m00 = m01;
+    m11 *= frct_phi;
+    m01 *= frct_phi;
+    frct_phi = 1. - frct_phi;
+    m00 *= frct_phi;
+    m10 *= frct_phi;
+    // divide by 1000 to get m/year
+    vx = (m00 * dx1 + m10 * dx2 + m01 * dx3 + m11 * dx4) / 1000;
+    vy = (m00 * dy1 + m10 * dy2 + m01 * dy3 + m11 * dy4) / 1000;
+    vz = (m00 * dz1 + m10 * dz2 + m01 * dz3 + m11 * dz4) / 1000;
+    return true;
+}
+
 /********************************************************************************/
 static PJ_XYZ get_grid_shift(PJ* P, const PJ_XYZ& cartesian) {
 /********************************************************************************
@@ -96,17 +211,33 @@ static PJ_XYZ get_grid_shift(PJ* P, const PJ_XYZ& cartesian) {
     geodetic.lpz = pj_inv3d(cartesian, Q->cart);
 
     /* look up correction values in grids */
-    shift.lp    = proj_hgrid_value(P, Q->hgrids, geodetic.lp);
-    shift.enu.u = proj_vgrid_value(P, Q->vgrids, geodetic.lp, 1.0);
+    if( !Q->grids.empty() )
+    {
+        double vx = 0;
+        double vy = 0;
+        double vz = 0;
+        if( !get_grid_values(P, Q, geodetic.lp, vx, vy, vz) )
+        {
+            return proj_coord_error().xyz;
+        }
+        shift.xyz.x = vx;
+        shift.xyz.y = vy;
+        shift.xyz.z = vz;
+    }
+    else
+    {
+        shift.lp    = proj_hgrid_value(P, Q->hgrids, geodetic.lp);
+        shift.enu.u = proj_vgrid_value(P, Q->vgrids, geodetic.lp, 1.0);
 
-    if (proj_errno(P) == PJD_ERR_GRID_AREA)
-        proj_log_debug(P, "deformation: coordinate (%.3f, %.3f) outside deformation model",
-                       proj_todeg(geodetic.lpz.lam), proj_todeg(geodetic.lpz.phi));
+        if (proj_errno(P) == PJD_ERR_GRID_AREA)
+            proj_log_debug(P, "deformation: coordinate (%.3f, %.3f) outside deformation model",
+                        proj_todeg(geodetic.lpz.lam), proj_todeg(geodetic.lpz.phi));
 
-    /* grid values are stored as mm/yr, we need m/yr */
-    shift.xyz.x /= 1000;
-    shift.xyz.y /= 1000;
-    shift.xyz.z /= 1000;
+        /* grid values are stored as mm/yr, we need m/yr */
+        shift.xyz.x /= 1000;
+        shift.xyz.y /= 1000;
+        shift.xyz.z /= 1000;
+    }
 
     /* pre-calc cosines and sines */
     sp = sin(geodetic.lpz.phi);
@@ -136,6 +267,9 @@ static PJ_XYZ reverse_shift(PJ *P, PJ_XYZ input, double dt) {
     int i = MAX_ITERATIONS;
 
     delta = get_grid_shift(P, input);
+    if (delta.x == HUGE_VAL) {
+        return delta;
+    }
 
     /* Store the origial z shift for later application */
     z0 = delta.z;
@@ -182,6 +316,9 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
     }
 
     shift = get_grid_shift(P, in.xyz);
+    if (shift.x == HUGE_VAL) {
+        return shift;
+    }
 
     out.xyz.x += Q->dt * shift.x;
     out.xyz.y += Q->dt * shift.y;
@@ -264,8 +401,6 @@ static PJ *destructor(PJ *P, int errlev) {
 
 
 PJ *TRANSFORMATION(deformation,1) {
-    int has_xy_grids = 0;
-    int has_z_grids  = 0;
     auto Q = new deformationData;
     P->opaque = (void *) Q;
     P->destructor = destructor;
@@ -278,25 +413,38 @@ PJ *TRANSFORMATION(deformation,1) {
     /* inherit ellipsoid definition from P to Q->cart */
     pj_inherit_ellipsoid_def (P, Q->cart);
 
-    has_xy_grids = pj_param(P->ctx, P->params, "txy_grids").i;
-    has_z_grids  = pj_param(P->ctx, P->params, "tz_grids").i;
+    int has_xy_grids = pj_param(P->ctx, P->params, "txy_grids").i;
+    int has_z_grids  = pj_param(P->ctx, P->params, "tz_grids").i;
+    int has_grids  = pj_param(P->ctx, P->params, "tgrids").i;
 
     /* Build gridlists. Both horizontal and vertical grids are mandatory. */
-    if (!has_xy_grids || !has_z_grids) {
-        proj_log_error(P, "deformation: Both +xy_grids and +z_grids should be specified.");
+    if ( !has_grids && (!has_xy_grids || !has_z_grids)) {
+        proj_log_error(P, "deformation: Either +grids or (+xy_grids and +z_grids) should be specified.");
         return destructor(P, PJD_ERR_NO_ARGS );
     }
 
-    Q->hgrids = proj_hgrid_init(P, "xy_grids");
-    if (proj_errno(P)) {
-        proj_log_error(P, "deformation: could not find requested xy_grid(s).");
-        return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+    if( has_grids )
+    {
+        Q->grids = proj_generic_grid_init(P, "grids");
+        /* Was gridlist compiled properly? */
+        if ( proj_errno(P) ) {
+            proj_log_error(P, "deformation: could not find required grid(s).");
+            return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+        }
     }
-
-    Q->vgrids = proj_vgrid_init(P, "z_grids");
-    if (proj_errno(P)) {
-        proj_log_error(P, "deformation: could not find requested z_grid(s).");
-        return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+    else
+    {
+        Q->hgrids = proj_hgrid_init(P, "xy_grids");
+        if (proj_errno(P)) {
+            proj_log_error(P, "deformation: could not find requested xy_grid(s).");
+            return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+        }
+
+        Q->vgrids = proj_vgrid_init(P, "z_grids");
+        if (proj_errno(P)) {
+            proj_log_error(P, "deformation: could not find requested z_grid(s).");
+            return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+        }
     }
 
     Q->dt = HUGE_VAL;
diff --git a/test/gie/deformation.gie b/test/gie/deformation.gie
index 37f62e6e..848b9e89 100644
--- a/test/gie/deformation.gie
+++ b/test/gie/deformation.gie
@@ -12,6 +12,35 @@ The input coordinate is located at lon=60, lam=-160 - somewhere in Alaska.
 <gie>
 
 -------------------------------------------------------------------------------
+Test with an extract of nkgrf03vel_realigned with ctable2+gtx
+-------------------------------------------------------------------------------
+operation    +proj=pipeline
+             +step +proj=cart +ellps=GRS80
+             +step +proj=deformation
+             +xy_grids=tests/nkgrf03vel_realigned_xy_extract.ct2
+             +z_grids=tests/nkgrf03vel_realigned_z_extract.gtx +ellps=GRS80 +dt=1
+             +step +proj=cart +ellps=GRS80 +inv
+-------------------------------------------------------------------------------
+tolerance   0.05 mm
+accept      21.5            63                   0
+expect      21.5000000049   62.9999999937        0.0083
+roundtrip   5
+
+-------------------------------------------------------------------------------
+Test with an extract of nkgrf03vel_realigned with GeoTIFF
+-------------------------------------------------------------------------------
+operation    +proj=pipeline
+             +step +proj=cart +ellps=GRS80
+             +step +proj=deformation
+             +grids=tests/nkgrf03vel_realigned_extract.tif +ellps=GRS80 +dt=1
+             +step +proj=cart +ellps=GRS80 +inv
+-------------------------------------------------------------------------------
+tolerance   0.05 mm
+accept      21.5            63                   0
+expect      21.5000000049   62.9999999937        0.0083
+roundtrip   5
+
+-------------------------------------------------------------------------------
 Test the +dt parameter
 -------------------------------------------------------------------------------
 operation   +proj=deformation +xy_grids=alaska +z_grids=egm96_15.gtx