Merge RFC4 (#1865)

This commit is the result of the squashing of rfc4_dev branch in a single
commit. It implements mostly RFC 4 related work.

* Grid handling:
  - remove obsolete and presumably unfinished implementation of grid catalog functionality
  - all grid functionality is in grids.cpp/.hpp
  - vertical and horizontal grid shift: rework to no longer load whole grid into memory
  - remove hgrids and vgrids member from PJ structure, and store them in hgridshift/vgridshift/deformation structures
  - build systems: add optional libtiff dependency. Must be explicitly disabled if not desired
  - add support for horizontal and vertical grids in GeoTIFF, if libtiff is available
  - add GenericShiftGridSet and GenericShiftGrid classes, relying on TIFF grids, that can be used for generic purpose grid-based adjustment
  - add a +proj=xyzgridshift method to perform geocentric translation by grid. Used for French NTF to RGF93 transformation using gr3df97a.tif grid
  - deformation: add support for +grids= for GeoTIFF grids
  - horizontal grid shift: fix failures on points slightly outside a subgrid (fixes #209)

* File management:
  - add a filemanager.cpp/.hpp to deal with file related work
  - test for legacy proj_api.h fileapi
  - proj.h: add proj_context_set_fileapi() and proj_context_set_sqlite3_vfs_name() (fixes #866)
  - add capability to read resource files from the user writable directory

* Network access:
  - build systems: add optional curl dependency
  - add a curl-based default implementation for network related functionality
  - proj.h: add C API to control network functionality, and optionaly provide network callbacks
  - add data/proj.ini with default settings
  - add a SQLite3 local cache of downloaded chunks
  - add proj_is_download_needed() and proj_download_file()

* Use Win32 Unicode APIs and expect all strings to be UTF-8 (fixes #1765)
  For backward compatibility, if PROJ_LIB content is found to be not UTF-8 or
  pointing to a non existing directory, then an attempt at interpretating it
  in the ANSI page encoding is done.
  proj_context_set_search_paths() now assumes strings to be in UTF-8, and
  functions returning paths will also return values in UTF-8.

4 files changed, 615 insertions, 89 deletions

diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp
index f1311a54..8aee50c9 100644
--- a/src/transformations/deformation.cpp
+++ b/src/transformations/deformation.cpp
@@ -56,20 +56,90 @@ grid-values in units of mm/year in ENU-space.
 #include "proj.h"
 #include "proj_internal.h"
 #include <math.h>
+#include "grids.hpp"
+
+#include <algorithm>
 
 PROJ_HEAD(deformation, "Kinematic grid shift");
 
 #define TOL 1e-8
 #define MAX_ITERATIONS 10
 
+using namespace NS_PROJ;
+
 namespace { // anonymous namespace
-struct pj_opaque {
-    double dt;
-    double t_epoch;
-    PJ *cart;
+struct deformationData {
+    double dt = 0;
+    double t_epoch = 0;
+    PJ *cart = nullptr;
+    ListOfGenericGrids grids{};
+    ListOfHGrids hgrids{};
+    ListOfVGrids vgrids{};
 };
 } // anonymous namespace
 
+// ---------------------------------------------------------------------------
+
+static bool get_grid_values(PJ* P,
+                            deformationData* Q,
+                            const PJ_LP& lp,
+                            double& vx,
+                            double& vy,
+                            double& vz)
+{
+    GenericShiftGridSet* gridset = nullptr;
+    auto grid = pj_find_generic_grid(Q->grids, lp, gridset);
+    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 != "millimetres per year" ) {
+        proj_log_error(P, "deformation: Only unit=millimetres per year currently handled");
+        return false;
+    }
+
+    bool must_retry = false;
+    if( !pj_bilinear_interpolation_three_samples(grid, lp,
+                                                 sampleE, sampleN, sampleU,
+                                                 vx, vy, vz,
+                                                 must_retry) )
+    {
+        if( must_retry )
+            return get_grid_values( P, Q, lp, vx, vy, vz);
+        return false;
+    }
+    // divide by 1000 to get m/year
+    vx /= 1000;
+    vy /= 1000;
+    vz /= 1000;
+    return true;
+}
+
 /********************************************************************************/
 static PJ_XYZ get_grid_shift(PJ* P, const PJ_XYZ& cartesian) {
 /********************************************************************************
@@ -85,22 +155,39 @@ static PJ_XYZ get_grid_shift(PJ* P, const PJ_XYZ& cartesian) {
     PJ_COORD geodetic, shift, temp;
     double sp, cp, sl, cl;
     int previous_errno = proj_errno_reset(P);
+    auto Q = static_cast<deformationData*>(P->opaque);
 
     /* cartesian to geodetic */
-    geodetic.lpz = pj_inv3d(cartesian, static_cast<struct pj_opaque*>(P->opaque)->cart);
+    geodetic.lpz = pj_inv3d(cartesian, Q->cart);
 
     /* look up correction values in grids */
-    shift.lp    = proj_hgrid_value(P, geodetic.lp);
-    shift.enu.u = proj_vgrid_value(P, 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));
-
-    /* grid values are stored as mm/yr, we need m/yr */
-    shift.xyz.x /= 1000;
-    shift.xyz.y /= 1000;
-    shift.xyz.z /= 1000;
+    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    = pj_hgrid_value(P, Q->hgrids, geodetic.lp);
+        shift.enu.u = pj_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));
+
+        /* 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);
@@ -130,6 +217,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;
@@ -163,7 +253,7 @@ static PJ_XYZ reverse_shift(PJ *P, PJ_XYZ input, double dt) {
 }
 
 static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
-    struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+    struct deformationData *Q = (struct deformationData *) P->opaque;
     PJ_COORD out, in;
     PJ_XYZ shift;
     in.lpz = lpz;
@@ -176,6 +266,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;
@@ -186,7 +279,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
 
 
 static PJ_COORD forward_4d(PJ_COORD in, PJ *P) {
-    struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+    struct deformationData *Q = (struct deformationData *) P->opaque;
     double dt;
     PJ_XYZ shift;
     PJ_COORD out = in;
@@ -209,7 +302,7 @@ static PJ_COORD forward_4d(PJ_COORD in, PJ *P) {
 
 
 static PJ_LPZ reverse_3d(PJ_XYZ in, PJ *P) {
-    struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+    struct deformationData *Q = (struct deformationData *) P->opaque;
     PJ_COORD out;
     out.xyz = in;
 
@@ -225,7 +318,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ in, PJ *P) {
 }
 
 static PJ_COORD reverse_4d(PJ_COORD in, PJ *P) {
-    struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+    struct deformationData *Q = (struct deformationData *) P->opaque;
     PJ_COORD out = in;
     double dt;
 
@@ -244,23 +337,23 @@ static PJ *destructor(PJ *P, int errlev) {
     if (nullptr==P)
         return nullptr;
 
-    if (nullptr==P->opaque)
-        return pj_default_destructor (P, errlev);
-
-    if (static_cast<struct pj_opaque*>(P->opaque)->cart)
-        static_cast<struct pj_opaque*>(P->opaque)->cart->destructor (static_cast<struct pj_opaque*>(P->opaque)->cart, errlev);
+    auto Q = static_cast<struct deformationData*>(P->opaque);
+    if( Q )
+    {
+        if (Q->cart)
+            Q->cart->destructor (Q->cart, errlev);
+        delete Q;
+    }
+    P->opaque = nullptr;
 
     return pj_default_destructor(P, errlev);
 }
 
 
 PJ *TRANSFORMATION(deformation,1) {
-    int has_xy_grids = 0;
-    int has_z_grids  = 0;
-    struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
-    if (nullptr==Q)
-        return destructor(P, ENOMEM);
+    auto Q = new deformationData;
     P->opaque = (void *) Q;
+    P->destructor = destructor;
 
     // Pass a dummy ellipsoid definition that will be overridden just afterwards
     Q->cart = proj_create(P->ctx, "+proj=cart +a=1");
@@ -270,25 +363,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 );
     }
 
-    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 = pj_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);
+        }
     }
-
-    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 = pj_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 = pj_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;
@@ -325,7 +431,6 @@ PJ *TRANSFORMATION(deformation,1) {
 
     P->left  = PJ_IO_UNITS_CARTESIAN;
     P->right = PJ_IO_UNITS_CARTESIAN;
-    P->destructor = destructor;
 
     return P;
 }
diff --git a/src/transformations/hgridshift.cpp b/src/transformations/hgridshift.cpp
index 90633939..122a7ab2 100644
--- a/src/transformations/hgridshift.cpp
+++ b/src/transformations/hgridshift.cpp
@@ -6,24 +6,38 @@
 #include <time.h>
 
 #include "proj_internal.h"
+#include "grids.hpp"
 
 PROJ_HEAD(hgridshift, "Horizontal grid shift");
 
+using namespace NS_PROJ;
+
 namespace { // anonymous namespace
-struct pj_opaque_hgridshift {
-    double t_final;
-    double t_epoch;
+struct hgridshiftData {
+    double t_final = 0;
+    double t_epoch = 0;
+    ListOfHGrids grids{};
+    bool defer_grid_opening = false;
 };
 } // anonymous namespace
 
 static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
+    auto Q = static_cast<hgridshiftData*>(P->opaque);
     PJ_COORD point = {{0,0,0,0}};
     point.lpz = lpz;
 
-    if (P->gridlist != nullptr) {
+    if ( Q->defer_grid_opening ) {
+        Q->defer_grid_opening = false;
+        Q->grids = pj_hgrid_init(P, "grids");
+        if ( proj_errno(P) ) {
+            return proj_coord_error().xyz;
+        }
+    }
+
+    if (!Q->grids.empty()) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        point.lp = proj_hgrid_apply(P, point.lp, PJ_FWD);
+        point.lp = pj_hgrid_apply(P->ctx, Q->grids, point.lp, PJ_FWD);
     }
 
     return point.xyz;
@@ -31,20 +45,29 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
 
 
 static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
+    auto Q = static_cast<hgridshiftData*>(P->opaque);
     PJ_COORD point = {{0,0,0,0}};
     point.xyz = xyz;
 
-    if (P->gridlist != nullptr) {
+    if ( Q->defer_grid_opening ) {
+        Q->defer_grid_opening = false;
+        Q->grids = pj_hgrid_init(P, "grids");
+        if ( proj_errno(P) ) {
+            return proj_coord_error().lpz;
+        }
+    }
+
+    if (!Q->grids.empty()) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        point.lp = proj_hgrid_apply(P, point.lp, PJ_INV);
+        point.lp = pj_hgrid_apply(P->ctx, Q->grids, point.lp, PJ_INV);
     }
 
     return point.lpz;
 }
 
 static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
-    struct pj_opaque_hgridshift *Q = (struct pj_opaque_hgridshift *) P->opaque;
+    struct hgridshiftData *Q = (struct hgridshiftData *) P->opaque;
     PJ_COORD point = obs;
 
     /* If transformation is not time restricted, we always call it */
@@ -62,7 +85,7 @@ static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
 }
 
 static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
-    struct pj_opaque_hgridshift *Q = (struct pj_opaque_hgridshift *) P->opaque;
+    struct hgridshiftData *Q = (struct hgridshiftData *) P->opaque;
     PJ_COORD point = obs;
 
     /* If transformation is not time restricted, we always call it */
@@ -78,12 +101,29 @@ static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
     return point;
 }
 
+static PJ *destructor (PJ *P, int errlev) {
+    if (nullptr==P)
+        return nullptr;
+
+    delete static_cast<struct hgridshiftData*>(P->opaque);
+    P->opaque = nullptr;
+
+    return pj_default_destructor(P, errlev);
+}
+
+static void reassign_context( PJ* P, PJ_CONTEXT* ctx )
+{
+    auto Q = (struct hgridshiftData *) P->opaque;
+    for( auto& grid: Q->grids ) {
+        grid->reassign_context(ctx);
+    }
+}
 
 PJ *TRANSFORMATION(hgridshift,0) {
-    struct pj_opaque_hgridshift *Q = static_cast<struct pj_opaque_hgridshift*>(pj_calloc (1, sizeof (struct pj_opaque_hgridshift)));
-    if (nullptr==Q)
-        return pj_default_destructor (P, ENOMEM);
+    auto Q = new hgridshiftData;
     P->opaque = (void *) Q;
+    P->destructor = destructor;
+    P->reassign_context = reassign_context;
 
     P->fwd4d  = forward_4d;
     P->inv4d  = reverse_4d;
@@ -97,12 +137,12 @@ PJ *TRANSFORMATION(hgridshift,0) {
 
     if (0==pj_param(P->ctx, P->params, "tgrids").i) {
         proj_log_error(P, "hgridshift: +grids parameter missing.");
-        return pj_default_destructor (P, PJD_ERR_NO_ARGS);
+        return destructor (P, PJD_ERR_NO_ARGS);
     }
 
-   /* TODO: Refactor into shared function that can be used  */
-   /*       by both vgridshift and hgridshift               */
-   if (pj_param(P->ctx, P->params, "tt_final").i) {
+    /* TODO: Refactor into shared function that can be used  */
+    /*       by both vgridshift and hgridshift               */
+    if (pj_param(P->ctx, P->params, "tt_final").i) {
         Q->t_final = pj_param (P->ctx, P->params, "dt_final").f;
         if (Q->t_final == 0) {
             /* a number wasn't passed to +t_final, let's see if it was "now" */
@@ -117,16 +157,21 @@ PJ *TRANSFORMATION(hgridshift,0) {
         }
     }
 
-   if (pj_param(P->ctx, P->params, "tt_epoch").i)
+    if (pj_param(P->ctx, P->params, "tt_epoch").i)
         Q->t_epoch = pj_param (P->ctx, P->params, "dt_epoch").f;
 
 
-    proj_hgrid_init(P, "grids");
-    /* Was gridlist compiled properly? */
-    if ( proj_errno(P) ) {
-        proj_log_error(P, "hgridshift: could not find required grid(s).");
-        return pj_default_destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+    if( P->ctx->defer_grid_opening ) {
+        Q->defer_grid_opening = true;
     }
+    else {
+        Q->grids = pj_hgrid_init(P, "grids");
+        /* Was gridlist compiled properly? */
+        if ( proj_errno(P) ) {
+            proj_log_error(P, "hgridshift: could not find required grid(s).");
+            return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+        }
+     }
 
     return P;
 }
diff --git a/src/transformations/vgridshift.cpp b/src/transformations/vgridshift.cpp
index de0cdd8c..121b795a 100644
--- a/src/transformations/vgridshift.cpp
+++ b/src/transformations/vgridshift.cpp
@@ -6,26 +6,67 @@
 #include <time.h>
 
 #include "proj_internal.h"
+#include "grids.hpp"
 
 PROJ_HEAD(vgridshift, "Vertical grid shift");
 
+using namespace NS_PROJ;
+
 namespace { // anonymous namespace
-struct pj_opaque_vgridshift {
-    double t_final;
-    double t_epoch;
-    double forward_multiplier;
+struct vgridshiftData {
+    double t_final = 0;
+    double t_epoch = 0;
+    double forward_multiplier = 0;
+    ListOfVGrids grids{};
+    bool defer_grid_opening = false;
 };
 } // anonymous namespace
 
+static void deal_with_vertcon_gtx_hack(PJ *P)
+{
+    struct vgridshiftData *Q = (struct vgridshiftData *) P->opaque;
+    // The .gtx VERTCON files stored millimetres, but the .tif files
+    // are in metres.
+    if( Q->forward_multiplier != 0.001 ) {
+        return;
+    }
+    const char* gridname = pj_param(P->ctx, P->params, "sgrids").s;
+    if( !gridname ) {
+        return;
+    }
+    if( strcmp(gridname, "vertconw.gtx") != 0 &&
+        strcmp(gridname, "vertconc.gtx") != 0 &&
+        strcmp(gridname, "vertcone.gtx") != 0 ) {
+        return;
+    }
+    if( Q->grids.empty() ) {
+        return;
+    }
+    const auto& grids = Q->grids[0]->grids();
+    if( !grids.empty() &&
+        grids[0]->name().find(".tif") != std::string::npos ) {
+        Q->forward_multiplier = 1.0;
+    }
+}
+
 static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
-    struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) P->opaque;
+    struct vgridshiftData *Q = (struct vgridshiftData *) P->opaque;
     PJ_COORD point = {{0,0,0,0}};
     point.lpz = lpz;
 
-    if (P->vgridlist_geoid != nullptr) {
+    if ( Q->defer_grid_opening ) {
+        Q->defer_grid_opening = false;
+        Q->grids = pj_vgrid_init(P, "grids");
+        deal_with_vertcon_gtx_hack(P);
+        if ( proj_errno(P) ) {
+            return proj_coord_error().xyz;
+        }
+    }
+
+    if (!Q->grids.empty()) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        point.xyz.z += proj_vgrid_value(P, point.lp, Q->forward_multiplier);
+        point.xyz.z += pj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier);
     }
 
     return point.xyz;
@@ -33,14 +74,23 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
 
 
 static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
-    struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) P->opaque;
+    struct vgridshiftData *Q = (struct vgridshiftData *) P->opaque;
     PJ_COORD point = {{0,0,0,0}};
     point.xyz = xyz;
 
-    if (P->vgridlist_geoid != nullptr) {
+    if ( Q->defer_grid_opening ) {
+        Q->defer_grid_opening = false;
+        Q->grids = pj_vgrid_init(P, "grids");
+        deal_with_vertcon_gtx_hack(P);
+        if ( proj_errno(P) ) {
+            return proj_coord_error().lpz;
+        }
+    }
+
+    if (!Q->grids.empty()) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        point.xyz.z -= proj_vgrid_value(P, point.lp, Q->forward_multiplier);
+        point.xyz.z -= pj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier);
     }
 
     return point.lpz;
@@ -48,7 +98,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
 
 
 static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
-    struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) P->opaque;
+    struct vgridshiftData *Q = (struct vgridshiftData *) P->opaque;
     PJ_COORD point = obs;
 
     /* If transformation is not time restricted, we always call it */
@@ -66,7 +116,7 @@ static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
 }
 
 static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
-    struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) P->opaque;
+    struct vgridshiftData *Q = (struct vgridshiftData *) P->opaque;
     PJ_COORD point = obs;
 
     /* If transformation is not time restricted, we always call it */
@@ -82,16 +132,34 @@ static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
     return point;
 }
 
+static PJ *destructor (PJ *P, int errlev) {
+    if (nullptr==P)
+        return nullptr;
+
+    delete static_cast<struct vgridshiftData*>(P->opaque);
+    P->opaque = nullptr;
+
+    return pj_default_destructor(P, errlev);
+}
+
+static void reassign_context( PJ* P, PJ_CONTEXT* ctx )
+{
+    auto Q = (struct vgridshiftData *) P->opaque;
+    for( auto& grid: Q->grids ) {
+        grid->reassign_context(ctx);
+    }
+}
+
 
 PJ *TRANSFORMATION(vgridshift,0) {
-    struct pj_opaque_vgridshift *Q = static_cast<struct pj_opaque_vgridshift*>(pj_calloc (1, sizeof (struct pj_opaque_vgridshift)));
-    if (nullptr==Q)
-        return pj_default_destructor (P, ENOMEM);
+    auto Q = new vgridshiftData;
     P->opaque = (void *) Q;
+    P->destructor = destructor;
+    P->reassign_context = reassign_context;
 
    if (!pj_param(P->ctx, P->params, "tgrids").i) {
         proj_log_error(P, "vgridshift: +grids parameter missing.");
-        return pj_default_destructor(P, PJD_ERR_NO_ARGS);
+        return destructor(P, PJD_ERR_NO_ARGS);
     }
 
    /* TODO: Refactor into shared function that can be used  */
@@ -120,13 +188,18 @@ PJ *TRANSFORMATION(vgridshift,0) {
         Q->forward_multiplier = pj_param(P->ctx, P->params, "dmultiplier").f;
     }
 
-    /* Build gridlist. P->vgridlist_geoid can be empty if +grids only ask for optional grids. */
-    proj_vgrid_init(P, "grids");
-
-    /* Was gridlist compiled properly? */
-    if ( proj_errno(P) ) {
-        proj_log_error(P, "vgridshift: could not find required grid(s).");
-        return pj_default_destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+    if( P->ctx->defer_grid_opening ) {
+        Q->defer_grid_opening = true;
+    }
+    else {
+        /* Build gridlist. P->vgridlist_geoid can be empty if +grids only ask for optional grids. */
+        Q->grids = pj_vgrid_init(P, "grids");
+
+        /* Was gridlist compiled properly? */
+        if ( proj_errno(P) ) {
+            proj_log_error(P, "vgridshift: could not find required grid(s).");
+            return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+        }
     }
 
     P->fwd4d = forward_4d;
diff --git a/src/transformations/xyzgridshift.cpp b/src/transformations/xyzgridshift.cpp
new file mode 100644
index 00000000..e1c76518
--- /dev/null
+++ b/src/transformations/xyzgridshift.cpp
@@ -0,0 +1,303 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Geocentric translation using a grid
+ * Author:   Even Rouault, <even.rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2019, 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 <string.h>
+#include <stddef.h>
+#include <time.h>
+
+#include "proj_internal.h"
+#include "grids.hpp"
+
+#include <algorithm>
+
+PROJ_HEAD(xyzgridshift, "Geocentric grid shift");
+
+using namespace NS_PROJ;
+
+namespace { // anonymous namespace
+struct xyzgridshiftData {
+    PJ *cart = nullptr;
+    bool grid_ref_is_input = true;
+    ListOfGenericGrids grids{};
+    bool defer_grid_opening = false;
+    double multiplier = 1.0;
+};
+} // anonymous namespace
+
+// ---------------------------------------------------------------------------
+
+static bool get_grid_values(PJ* P,
+                            xyzgridshiftData* Q,
+                            const PJ_LP& lp,
+                            double& dx,
+                            double& dy,
+                            double& dz)
+{
+    if ( Q->defer_grid_opening ) {
+        Q->defer_grid_opening = false;
+        Q->grids = pj_generic_grid_init(P, "grids");
+        if ( proj_errno(P) ) {
+            return false;
+        }
+    }
+
+    GenericShiftGridSet* gridset = nullptr;
+    auto grid = pj_find_generic_grid(Q->grids, lp, gridset);
+    if( !grid ) {
+        return false;
+    }
+    if( grid->isNullGrid() ) {
+        dx = 0;
+        dy = 0;
+        dz = 0;
+        return true;
+    }
+    const auto samplesPerPixel = grid->samplesPerPixel();
+    if( samplesPerPixel < 3 ) {
+        proj_log_error(P, "xyzgridshift: grid has not enough samples");
+        return false;
+    }
+    int sampleX = 0;
+    int sampleY = 1;
+    int sampleZ = 2;
+    for( int i = 0; i < samplesPerPixel; i++ )
+    {
+        const auto desc = grid->description(i);
+        if( desc == "x_translation") {
+            sampleX = i;
+        } else if( desc == "y_translation") {
+            sampleY = i;
+        } else if( desc == "z_translation") {
+            sampleZ = i;
+        }
+    }
+    const auto unit = grid->unit(sampleX);
+    if( !unit.empty() && unit != "metre" ) {
+        proj_log_error(P, "xyzgridshift: Only unit=metre currently handled");
+        return false;
+    }
+
+    bool must_retry = false;
+    if( !pj_bilinear_interpolation_three_samples(grid, lp,
+                                                 sampleX, sampleY, sampleZ,
+                                                 dx, dy, dz,
+                                                 must_retry) )
+    {
+        if( must_retry )
+            return get_grid_values( P, Q, lp, dx, dy, dz);
+        return false;
+    }
+
+    dx *= Q->multiplier;
+    dy *= Q->multiplier;
+    dz *= Q->multiplier;
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+#define SQUARE(x) ((x)*(x))
+
+// ---------------------------------------------------------------------------
+
+static PJ_COORD iterative_adjustment(PJ* P,
+                                     xyzgridshiftData* Q,
+                                     const PJ_COORD& pointInit,
+                                     double factor)
+{
+    PJ_COORD point = pointInit;
+    for(int i = 0; i < 10; i++) {
+        PJ_COORD geodetic;
+        geodetic.lpz = pj_inv3d(point.xyz, Q->cart);
+
+        double dx, dy, dz;
+        if( !get_grid_values(P, Q, geodetic.lp, dx, dy, dz) ) {
+            return proj_coord_error();
+        }
+
+        dx *= factor;
+        dy *= factor;
+        dz *= factor;
+
+        const double err = SQUARE((point.xyz.x - pointInit.xyz.x) - dx) +
+                            SQUARE((point.xyz.y - pointInit.xyz.y) - dy) +
+                            SQUARE((point.xyz.z - pointInit.xyz.z) - dz);
+
+        point.xyz.x = pointInit.xyz.x + dx;
+        point.xyz.y = pointInit.xyz.y + dy;
+        point.xyz.z = pointInit.xyz.z + dz;
+        if( err < 1e-10 ) {
+            break;
+        }
+    }
+    return point;
+}
+
+// ---------------------------------------------------------------------------
+
+static PJ_COORD direct_adjustment(PJ* P,
+                                  xyzgridshiftData* Q,
+                                  PJ_COORD point,
+                                  double factor)
+{
+    PJ_COORD geodetic;
+    geodetic.lpz = pj_inv3d(point.xyz, Q->cart);
+
+    double dx, dy, dz;
+    if( !get_grid_values(P, Q, geodetic.lp, dx, dy, dz) ) {
+        return proj_coord_error();
+    }
+    point.xyz.x += factor * dx;
+    point.xyz.y += factor * dy;
+    point.xyz.z += factor * dz;
+    return point;
+}
+
+// ---------------------------------------------------------------------------
+
+static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
+    auto Q = static_cast<xyzgridshiftData*>(P->opaque);
+    PJ_COORD point = {{0,0,0,0}};
+    point.lpz = lpz;
+
+    if( Q->grid_ref_is_input ) {
+        point = direct_adjustment(P, Q, point, 1.0);
+    }
+    else {
+        point = iterative_adjustment(P, Q, point, 1.0);
+    }
+
+    return point.xyz;
+}
+
+
+static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
+    auto Q = static_cast<xyzgridshiftData*>(P->opaque);
+    PJ_COORD point = {{0,0,0,0}};
+    point.xyz = xyz;
+
+    if( Q->grid_ref_is_input ) {
+        point = iterative_adjustment(P, Q, point, -1.0);
+    }
+    else {
+        point = direct_adjustment(P, Q, point, -1.0);
+    }
+
+    return point.lpz;
+}
+
+static PJ *destructor (PJ *P, int errlev) {
+    if (nullptr==P)
+        return nullptr;
+
+    auto Q = static_cast<struct xyzgridshiftData*>(P->opaque);
+    if( Q )
+    {
+        if (Q->cart)
+            Q->cart->destructor (Q->cart, errlev);
+        delete Q;
+    }
+    P->opaque = nullptr;
+
+    return pj_default_destructor(P, errlev);
+}
+
+static void reassign_context( PJ* P, PJ_CONTEXT* ctx )
+{
+    auto Q = (struct xyzgridshiftData *) P->opaque;
+    for( auto& grid: Q->grids ) {
+        grid->reassign_context(ctx);
+    }
+}
+
+
+PJ *TRANSFORMATION(xyzgridshift,0) {
+    auto Q = new xyzgridshiftData;
+    P->opaque = (void *) Q;
+    P->destructor = destructor;
+    P->reassign_context = reassign_context;
+
+    P->fwd4d  = nullptr;
+    P->inv4d  = nullptr;
+    P->fwd3d  = forward_3d;
+    P->inv3d  = reverse_3d;
+    P->fwd    = nullptr;
+    P->inv    = nullptr;
+
+    P->left  = PJ_IO_UNITS_CARTESIAN;
+    P->right = PJ_IO_UNITS_CARTESIAN;
+
+    // Pass a dummy ellipsoid definition that will be overridden just afterwards
+    Q->cart = proj_create(P->ctx, "+proj=cart +a=1");
+    if (Q->cart == nullptr)
+        return destructor(P, ENOMEM);
+
+    /* inherit ellipsoid definition from P to Q->cart */
+    pj_inherit_ellipsoid_def (P, Q->cart);
+
+    const char* grid_ref = pj_param (P->ctx, P->params, "sgrid_ref").s;
+    if( grid_ref ) {
+        if (strcmp(grid_ref, "input_crs") == 0 ) {
+            // default
+        } else if (strcmp(grid_ref, "output_crs") == 0 ) {
+            // Convention use for example for NTF->RGF93 grid that contains
+            // delta x,y,z from NTF to RGF93, but the grid itself is referenced
+            // in RGF93
+            Q->grid_ref_is_input = false;
+        } else {
+            proj_log_error(P, "xyzgridshift: unusupported value for grid_ref");
+            return destructor (P, PJD_ERR_NO_ARGS);
+        }
+    }
+
+    if (0==pj_param(P->ctx, P->params, "tgrids").i) {
+        proj_log_error(P, "xyzgridshift: +grids parameter missing.");
+        return destructor (P, PJD_ERR_NO_ARGS);
+    }
+
+    /* multiplier for delta x,y,z */
+    if (pj_param(P->ctx, P->params, "tmultiplier").i) {
+        Q->multiplier = pj_param(P->ctx, P->params, "dmultiplier").f;
+    }
+
+    if( P->ctx->defer_grid_opening ) {
+        Q->defer_grid_opening = true;
+    }
+    else {
+        Q->grids = pj_generic_grid_init(P, "grids");
+        /* Was gridlist compiled properly? */
+        if ( proj_errno(P) ) {
+            proj_log_error(P, "xyzgridshift: could not find required grid(s).");
+            return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+        }
+    }
+
+    return P;
+}