Merge pull request #1777 from rouault/rfc4_impl

[RFC4_dev] Grid loading refactoring

22 files changed, 1742 insertions, 2977 deletions

diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index f37594b5..efb4a86a 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -47,6 +47,7 @@
 #include "proj_internal.h"
 #include <math.h>
 #include "geodesic.h"
+#include "grids.hpp"
 
 #include "proj/common.hpp"
 #include "proj/coordinateoperation.hpp"
@@ -1553,43 +1554,65 @@ PJ_GRID_INFO proj_grid_info(const char *gridname) {
 
     /*PJ_CONTEXT *ctx = proj_context_create(); */
     PJ_CONTEXT *ctx = pj_get_default_ctx();
-    PJ_GRIDINFO *gridinfo = pj_gridinfo_init(ctx, gridname);
     memset(&grinfo, 0, sizeof(PJ_GRID_INFO));
 
-    /* in case the grid wasn't found */
-    if (gridinfo->filename == nullptr || gridinfo->ct == nullptr) {
-        pj_gridinfo_free(ctx, gridinfo);
-        strcpy(grinfo.format, "missing");
-        return grinfo;
-    }
-
-    /* The string copies below are automatically null-terminated due to */
-    /* the memset above, so strncpy is safe */
+    const auto fillGridInfo = [&grinfo, ctx, gridname]
+                        (const NS_PROJ::Grid& grid, const std::string& format)
+    {
+        const auto& extent = grid.extentAndRes();
 
-    /* name of grid */
-    strncpy (grinfo.gridname, gridname, sizeof(grinfo.gridname) - 1);
+        /* name of grid */
+        strncpy (grinfo.gridname, gridname, sizeof(grinfo.gridname) - 1);
 
-    /* full path of grid */
-    pj_find_file(ctx, gridname, grinfo.filename, sizeof(grinfo.filename) - 1);
+        /* full path of grid */
+        pj_find_file(ctx, gridname, grinfo.filename, sizeof(grinfo.filename) - 1);
 
-    /* grid format */
-    strncpy (grinfo.format, gridinfo->format, sizeof(grinfo.format) - 1);
+        /* grid format */
+        strncpy (grinfo.format, format.c_str(), sizeof(grinfo.format) - 1);
 
-    /* grid size */
-    grinfo.n_lon = gridinfo->ct->lim.lam;
-    grinfo.n_lat = gridinfo->ct->lim.phi;
+        /* grid size */
+        grinfo.n_lon = grid.width();
+        grinfo.n_lat = grid.height();
 
-    /* cell size */
-    grinfo.cs_lon = gridinfo->ct->del.lam;
-    grinfo.cs_lat = gridinfo->ct->del.phi;
+        /* cell size */
+        grinfo.cs_lon = extent.resLon;
+        grinfo.cs_lat = extent.resLat;
 
-    /* bounds of grid */
-    grinfo.lowerleft  = gridinfo->ct->ll;
-    grinfo.upperright.lam = grinfo.lowerleft.lam + grinfo.n_lon*grinfo.cs_lon;
-    grinfo.upperright.phi = grinfo.lowerleft.phi + grinfo.n_lat*grinfo.cs_lat;
+        /* bounds of grid */
+        grinfo.lowerleft.lam  = extent.westLon;
+        grinfo.lowerleft.phi  = extent.southLat;
+        grinfo.upperright.lam = extent.eastLon;
+        grinfo.upperright.phi = extent.northLat;
+    };
 
-    pj_gridinfo_free(ctx, gridinfo);
+    {
+        const auto gridSet = NS_PROJ::VerticalShiftGridSet::open(ctx, gridname);
+        if( gridSet )
+        {
+            const auto& grids = gridSet->grids();
+            if( !grids.empty() )
+            {
+                const auto& grid = grids.front();
+                fillGridInfo(*grid, gridSet->format());
+                return grinfo;
+            }
+        }
+    }
 
+    {
+        const auto gridSet = NS_PROJ::HorizontalShiftGridSet::open(ctx, gridname);
+        if( gridSet )
+        {
+            const auto& grids = gridSet->grids();
+            if( !grids.empty() )
+            {
+                const auto& grid = grids.front();
+                fillGridInfo(*grid, gridSet->format());
+                return grinfo;
+            }
+        }
+    }
+    strcpy(grinfo.format, "missing");
     return grinfo;
 }
 
diff --git a/src/Makefile.am b/src/Makefile.am
index 2e14d1dd..a12de4e1 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -195,10 +195,8 @@ libproj_la_SOURCES = \
 	release.cpp gauss.cpp \
 	fileapi.cpp \
 	\
-	gc_reader.cpp gridcatalog.cpp \
-	nad_cvt.cpp nad_init.cpp nad_intr.cpp \
 	apply_gridshift.cpp datums.cpp datum_set.cpp transform.cpp \
-	geocent.cpp geocent.h utils.cpp gridinfo.cpp gridlist.cpp \
+	geocent.cpp geocent.h utils.cpp \
 	jniproj.cpp mutex.cpp initcache.cpp apply_vgridshift.cpp geodesic.c \
 	strtod.cpp \
 	\
@@ -213,7 +211,10 @@ libproj_la_SOURCES = \
 	proj_json_streaming_writer.hpp \
 	proj_json_streaming_writer.cpp \
 	\
-	tracing.cpp
+	tracing.cpp \
+	\
+	grids.hpp \
+	grids.cpp
 
 
 # The sed hack is to please MSVC
diff --git a/src/apply_gridshift.cpp b/src/apply_gridshift.cpp
index e6c7b2b1..a0ffa394 100644
--- a/src/apply_gridshift.cpp
+++ b/src/apply_gridshift.cpp
@@ -28,7 +28,9 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
-#define PJ_LIB__
+#ifndef FROM_PROJ_CPP
+#define FROM_PROJ_CPP
+#endif
 
 #include <math.h>
 #include <stdio.h>
@@ -36,235 +38,60 @@
 
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj/internal/internal.hpp"
+#include "grids.hpp"
 
-/************************************************************************/
-/*                         pj_apply_gridshift()                         */
-/*                                                                      */
-/*      This is the externally callable interface - part of the         */
-/*      public API - though it is not used internally any more and I    */
-/*      doubt it is used by any other applications.  But we preserve    */
-/*      it to honour our public api.                                    */
-/************************************************************************/
+NS_PROJ_START
 
-int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse,
-                        long point_count, int point_offset,
-                        double *x, double *y, double *z )
+// ---------------------------------------------------------------------------
 
+static const HorizontalShiftGrid* findGrid(const ListOfHGrids& grids,
+                                           PJ_LP input)
 {
-    PJ_GRIDINFO **gridlist;
-    int           grid_count;
-    int           ret;
-
-    gridlist = pj_gridlist_from_nadgrids( ctx, nadgrids, &grid_count );
-
-    if( gridlist == nullptr || grid_count == 0 )
-    {
-        pj_dalloc( gridlist );
-        return ctx->last_errno;
-    }
-
-    ret = pj_apply_gridshift_3( ctx, gridlist, grid_count, inverse,
-                                point_count, point_offset, x, y, z );
-
-    /*
-    ** Note this frees the array of grid list pointers, but not the grids
-    ** which is as intended.  The grids themselves live on.
-    */
-    pj_dalloc( gridlist );
-
-    return ret;
-}
-
-/************************************************************************/
-/*                        pj_apply_gridshift_2()                        */
-/*                                                                      */
-/*      This implementation uses the gridlist from a coordinate         */
-/*      system definition.  If the gridlist has not yet been            */
-/*      populated in the coordinate system definition we set it up      */
-/*      now.                                                            */
-/************************************************************************/
-
-int pj_apply_gridshift_2( PJ *defn, int inverse,
-                          long point_count, int point_offset,
-                          double *x, double *y, double *z )
-
-{
-    if( defn->catalog_name != nullptr )
-        return pj_gc_apply_gridshift( defn, inverse, point_count, point_offset,
-                                      x, y, z );
-
-    if( defn->gridlist == nullptr )
-    {
-        defn->gridlist =
-            pj_gridlist_from_nadgrids( pj_get_ctx( defn ),
-                                       pj_param(defn->ctx, defn->params,"snadgrids").s,
-                                       &(defn->gridlist_count) );
-
-        if( defn->gridlist == nullptr || defn->gridlist_count == 0 )
-            return defn->ctx->last_errno;
-    }
-
-    return pj_apply_gridshift_3( pj_get_ctx( defn ),
-                                 defn->gridlist, defn->gridlist_count, inverse,
-                                 point_count, point_offset, x, y, z );
-}
-
-/************************************************************************/
-/*                             find_ctable()                            */
-/*                                                                      */
-/*    Determine which grid is the correct given an input coordinate.    */
-/************************************************************************/
-
-static struct CTABLE* find_ctable(projCtx ctx, PJ_LP input, int grid_count, PJ_GRIDINFO **tables) {
-    int itable;
-
-    /* keep trying till we find a table that works */
-    for( itable = 0; itable < grid_count; itable++ )
+    for( const auto& gridset: grids )
     {
-
-        PJ_GRIDINFO *gi = tables[itable];
-        struct CTABLE *ct = gi->ct;
-        double epsilon = (fabs(ct->del.phi)+fabs(ct->del.lam))/10000.0;
-        /* skip tables that don't match our point at all.  */
-        if ( ct->ll.phi - epsilon > input.phi
-                || ct->ll.lam - epsilon > input.lam
-                || (ct->ll.phi + (ct->lim.phi-1) * ct->del.phi + epsilon < input.phi)
-                || (ct->ll.lam + (ct->lim.lam-1) * ct->del.lam + epsilon < input.lam) ) {
-            continue;
-        }
-
-        /* If we have child nodes, check to see if any of them apply. */
-        while( gi->child )
-        {
-            PJ_GRIDINFO *child;
-
-            for( child = gi->child; child != nullptr; child = child->next )
-            {
-                struct CTABLE *ct1 = child->ct;
-                epsilon = (fabs(ct1->del.phi)+fabs(ct1->del.lam))/10000.0;
-
-                if( ct1->ll.phi - epsilon > input.phi
-                    || ct1->ll.lam - epsilon > input.lam
-                    || (ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi + epsilon < input.phi)
-                    || (ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam + epsilon < input.lam) ) {
-                    continue;
-                }
-                break;
-            }
-
-            /* If we didn't find a child then nothing more to do */
-            if( child == nullptr ) break;
-
-            /* Otherwise use the child, first checking it's children */
-            gi = child;
-            ct = child->ct;
-        }
-        /* load the grid shift info if we don't have it. */
-        if( ct->cvs == nullptr) {
-            if (!pj_gridinfo_load( ctx, gi ) ) {
-                pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-                return nullptr;
-            }
-        }
-        /* if we get this far we have found a suitable grid */
-        return ct;
+        auto grid = gridset->gridAt(input.lam, input.phi);
+        if( grid )
+            return grid;
     }
-
     return nullptr;
 }
 
-/************************************************************************/
-/*                        pj_apply_gridshift_3()                        */
-/*                                                                      */
-/*      This is the real workhorse, given a gridlist.                   */
-/************************************************************************/
+// ---------------------------------------------------------------------------
 
-int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **gridlist, int gridlist_count,
-                          int inverse, long point_count, int point_offset,
-                          double *x, double *y, double *z )
+static ListOfHGrids getListOfGridSets(PJ_CONTEXT* ctx, const char* grids)
 {
-    int  i;
-    struct CTABLE *ct;
-    static int debug_count = 0;
-    (void) z;
-
-    if( gridlist== nullptr || gridlist_count == 0 )
+    ListOfHGrids list;
+    auto listOfGrids = internal::split(std::string(grids), ',');
+    for( const auto& grid: listOfGrids )
     {
-        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
-        return PJD_ERR_FAILED_TO_LOAD_GRID;
-    }
-
-    ctx->last_errno = 0;
-
-    for( i = 0; i < point_count; i++ )
-    {
-        long io = i * point_offset;
-        PJ_LP   input, output;
-        int  itable;
-
-        input.phi = y[io];
-        input.lam = x[io];
-        output.phi = HUGE_VAL;
-        output.lam = HUGE_VAL;
-
-        ct = find_ctable(ctx, input, gridlist_count, gridlist);
-        if( ct != nullptr )
+        const char* gridname = grid.c_str();
+        bool canFail = false;
+        if( gridname[0] == '@' )
         {
-            output = nad_cvt( input, inverse, ct );
-
-            if ( output.lam != HUGE_VAL && debug_count++ < 20 )
-                pj_log( ctx, PJ_LOG_DEBUG_MINOR, "pj_apply_gridshift(): used %s", ct->id );
+            canFail = true;
+            gridname ++;
         }
-
-        if ( output.lam == HUGE_VAL )
+        auto gridSet = HorizontalShiftGridSet::open(ctx, gridname);
+        if( !gridSet )
         {
-            if( ctx->debug_level >= PJ_LOG_DEBUG_MAJOR )
+            if( !canFail )
             {
-                pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
-                    "pj_apply_gridshift(): failed to find a grid shift table for\n"
-                    "                      location (%.7fdW,%.7fdN)",
-                    x[io] * RAD_TO_DEG,
-                    y[io] * RAD_TO_DEG );
-                for( itable = 0; itable < gridlist_count; itable++ )
-                {
-                    PJ_GRIDINFO *gi = gridlist[itable];
-                    if( itable == 0 )
-                        pj_log( ctx, PJ_LOG_DEBUG_MAJOR, "   tried: %s", gi->gridname );
-                    else
-                        pj_log( ctx, PJ_LOG_DEBUG_MAJOR, ",%s", gi->gridname );
-                }
+                pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
+                return {};
             }
-
-            /*
-             * We don't actually have any machinery currently to set the
-             * following macro, so this is mostly kept here to make it clear
-             * how we ought to operate if we wanted to make it super clear
-             * that an error has occurred when points are outside our available
-             * datum shift areas.  But if this is on, we will find that "low
-             * value" points on the fringes of some datasets will completely
-             * fail causing lots of problems when it is more or less ok to
-             * just not apply a datum shift.  So rather than deal with
-             * that we just fallback to no shift. (see also bug #45).
-             */
-#ifdef ERR_GRID_AREA_TRANSIENT_SEVERE
-            y[io] = HUGE_VAL;
-            x[io] = HUGE_VAL;
-#else
-            /* leave x/y unshifted. */
-#endif
         }
         else
         {
-            y[io] = output.phi;
-            x[io] = output.lam;
+            list.emplace_back(std::move(gridSet));
         }
     }
-
-    return 0;
+    return list;
 }
 
+
 /**********************************************/
-int proj_hgrid_init(PJ* P, const char *grids) {
+ListOfHGrids proj_hgrid_init(PJ* P, const char *gridkey) {
 /**********************************************
 
   Initizalize and populate list of horizontal
@@ -279,29 +106,157 @@ int proj_hgrid_init(PJ* P, const char *grids) {
 
 ***********************************************/
 
-    /* prepend "s" to the "grids" string to allow usage with pj_param */
-    char *sgrids = (char *) pj_malloc( (strlen(grids)+1+1) *sizeof(char) );
-    sprintf(sgrids, "%s%s", "s", grids);
+    std::string key("s");
+    key += gridkey;
+    const char* grids = pj_param(P->ctx, P->params, key.c_str()).s;
+    if( grids == nullptr )
+        return {};
 
-    if (P->gridlist == nullptr) {
-        P->gridlist = pj_gridlist_from_nadgrids(
-            P->ctx,
-            pj_param(P->ctx, P->params, sgrids).s,
-            &(P->gridlist_count)
-        );
+    return getListOfGridSets(P->ctx, grids);
+}
 
-        if( P->gridlist == nullptr || P->gridlist_count == 0 ) {
-            pj_dealloc(sgrids);
-            return 0;
+typedef struct { pj_int32 lam, phi; } ILP;
+
+static PJ_LP nad_intr(PJ_LP t, const HorizontalShiftGrid* grid) {
+	PJ_LP val, frct;
+	ILP indx;
+	int in;
+
+        const auto& extent = grid->extentAndRes();
+	t.lam /= extent.resLon;
+	indx.lam = isnan(t.lam) ? 0 : (pj_int32)lround(floor(t.lam));
+	t.phi /= extent.resLat;
+	indx.phi = isnan(t.phi) ? 0 : (pj_int32)lround(floor(t.phi));
+
+	frct.lam = t.lam - indx.lam;
+	frct.phi = t.phi - indx.phi;
+	val.lam = val.phi = HUGE_VAL;
+	if (indx.lam < 0) {
+		if (indx.lam == -1 && frct.lam > 0.99999999999) {
+			++indx.lam;
+			frct.lam = 0.;
+		} else
+			return val;
+	} else if ((in = indx.lam + 1) >= grid->width()) {
+		if (in == grid->width() && frct.lam < 1e-11) {
+			--indx.lam;
+			frct.lam = 1.;
+		} else
+			return val;
+	}
+	if (indx.phi < 0) {
+		if (indx.phi == -1 && frct.phi > 0.99999999999) {
+			++indx.phi;
+			frct.phi = 0.;
+		} else
+			return val;
+	} else if ((in = indx.phi + 1) >= grid->height()) {
+		if (in == grid->height() && frct.phi < 1e-11) {
+			--indx.phi;
+			frct.phi = 1.;
+		} else
+			return val;
+	}
+	
+	float f00Lon = 0, f00Lat = 0;
+        float f10Lon = 0, f10Lat = 0;
+        float f01Lon = 0, f01Lat = 0;
+        float f11Lon = 0, f11Lat = 0;
+	if( !grid->valueAt(indx.lam, indx.phi, f00Lon, f00Lat) ||
+            !grid->valueAt(indx.lam + 1, indx.phi, f10Lon, f10Lat) ||
+            !grid->valueAt(indx.lam, indx.phi + 1, f01Lon, f01Lat) ||
+            !grid->valueAt(indx.lam + 1, indx.phi + 1, f11Lon, f11Lat) )
+        {
+            return val;
         }
+
+        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;
+	val.lam = m00 * f00Lon + m10 * f10Lon +
+			  m01 * f01Lon + m11 * f11Lon;
+	val.phi = m00 * f00Lat + m10 * f10Lat +
+			  m01 * f01Lat + m11 * f11Lat;
+	return val;
+}
+
+
+#define MAX_ITERATIONS 10
+#define TOL 1e-12
+
+static
+PJ_LP nad_cvt(PJ_LP in, int inverse, const HorizontalShiftGrid* grid) {
+    PJ_LP t, tb,del, dif;
+    int i = MAX_ITERATIONS;
+    const double toltol = TOL*TOL;
+
+    if (in.lam == HUGE_VAL)
+        return in;
+
+    /* normalize input to ll origin */
+    tb = in;
+    const auto& extent = grid->extentAndRes();
+    tb.lam -= extent.westLon;
+    tb.phi -= extent.southLat;
+
+    tb.lam = adjlon (tb.lam - M_PI) + M_PI;
+
+    t = nad_intr (tb, grid);
+    if (t.lam == HUGE_VAL)
+        return t;
+
+    if (!inverse) {
+        in.lam -= t.lam;
+        in.phi += t.phi;
+        return in;
     }
 
-    if (P->gridlist_count == 0) {
-        proj_errno_set(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+    t.lam = tb.lam + t.lam;
+    t.phi = tb.phi - t.phi;
+
+    do {
+        del = nad_intr(t, grid);
+
+        /* This case used to return failure, but I have
+           changed it to return the first order approximation
+           of the inverse shift.  This avoids cases where the
+           grid shift *into* this grid came from another grid.
+           While we aren't returning optimally correct results
+           I feel a close result in this case is better than
+           no result.  NFW
+           To demonstrate use -112.5839956 49.4914451 against
+           the NTv2 grid shift file from Canada. */
+        if (del.lam == HUGE_VAL)
+            break;
+
+        dif.lam = t.lam - del.lam - tb.lam;
+           dif.phi = t.phi + del.phi - tb.phi;
+        t.lam -= dif.lam;
+        t.phi -= dif.phi;
+
+    } while (--i && (dif.lam*dif.lam + dif.phi*dif.phi > toltol)); /* prob. slightly faster than hypot() */
+
+    if (i==0) {
+        /* If we had access to a context, this should go through pj_log, and we should set ctx->errno */
+        if (getenv ("PROJ_DEBUG"))
+            fprintf( stderr, "Inverse grid shift iterator failed to converge.\n" );
+        t.lam = t.phi = HUGE_VAL;
+        return t;
     }
 
-    pj_dealloc(sgrids);
-    return P->gridlist_count;
+    /* and again: pj_log and ctx->errno */
+    if (del.lam==HUGE_VAL && getenv ("PROJ_DEBUG"))
+        fprintf (stderr, "Inverse grid shift iteration failed, presumably at grid edge.\nUsing first approximation.\n");
+
+    in.lam = adjlon (t.lam + extent.westLon);
+    in.phi = t.phi + extent.southLat;
+    return in;
 }
 
 /********************************************/
@@ -309,23 +264,23 @@ int proj_hgrid_init(PJ* P, const char *grids) {
 /*                                          */
 /*    Return coordinate offset in grid      */
 /********************************************/
-PJ_LP proj_hgrid_value(PJ *P, PJ_LP lp) {
-    struct CTABLE *ct;
+PJ_LP proj_hgrid_value(PJ *P, const ListOfHGrids& grids, PJ_LP lp) {
     PJ_LP out = proj_coord_error().lp;
 
-    ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist);
-    if (ct == nullptr) {
-        pj_ctx_set_errno( P->ctx, PJD_ERR_GRID_AREA);
+    const auto grid = findGrid(grids, lp);
+    if( !grid ) {
+        pj_ctx_set_errno( P->ctx, PJD_ERR_GRID_AREA );
         return out;
     }
 
     /* normalize input to ll origin */
-    lp.lam -= ct->ll.lam;
-    lp.phi -= ct->ll.phi;
+    const auto& extent = grid->extentAndRes();
+    lp.lam -= extent.westLon;
+    lp.phi -= extent.southLat;
 
     lp.lam = adjlon(lp.lam - M_PI) + M_PI;
 
-    out = nad_intr(lp, ct);
+    out = nad_intr(lp, grid);
 
     if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) {
         pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA);
@@ -334,33 +289,92 @@ PJ_LP proj_hgrid_value(PJ *P, PJ_LP lp) {
     return out;
 }
 
-PJ_LP proj_hgrid_apply(PJ *P, PJ_LP lp, PJ_DIRECTION direction) {
-    struct CTABLE *ct;
-    int inverse;
+static
+PJ_LP proj_hgrid_apply_internal(PJ_CONTEXT* ctx,
+                       PJ_LP lp,
+                       PJ_DIRECTION direction,
+                       const ListOfHGrids& grids)
+{
     PJ_LP out;
 
-    out.lam = HUGE_VAL; out.phi = HUGE_VAL;
+    out.lam = HUGE_VAL;
+    out.phi = HUGE_VAL;
 
-    ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist);
-
-    if (ct == nullptr || ct->cvs == nullptr) {
-        if( P->gridlist_count == 1 &&
-            strcmp(P->gridlist[0]->gridname, "null") == 0) {
-          // TODO: remove this particular case that is put there just to be
-          // able to handle longitudes outside of -180,180
-          out = lp;
-        } else {
-          pj_ctx_set_errno( P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        }
+    const auto grid = findGrid(grids, lp);
+    if( !grid ) {
+        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
         return out;
     }
+    if( grid->isNullGrid() )
+    {
+        return lp;
+    }
 
-    inverse = direction == PJ_FWD ? 0 : 1;
-    out = nad_cvt(lp, inverse, ct);
+    int inverse = direction == PJ_FWD ? 0 : 1;
+    out = nad_cvt(lp, inverse, grid);
 
     if (out.lam == HUGE_VAL || out.phi == HUGE_VAL)
-        pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA);
+        pj_ctx_set_errno(ctx, PJD_ERR_GRID_AREA);
 
     return out;
+}
+
+PJ_LP proj_hgrid_apply(PJ *P, const ListOfHGrids& grids, PJ_LP lp, PJ_DIRECTION direction) {
+    return proj_hgrid_apply_internal(P->ctx, lp, direction, grids);
+}
+
+
+NS_PROJ_END
+
+/************************************************************************/
+/*                         pj_apply_gridshift()                         */
+/*                                                                      */
+/*      This is the externally callable interface - part of the         */
+/*      public API - though it is not used internally any more and I    */
+/*      doubt it is used by any other applications.  But we preserve    */
+/*      it to honour our public api.                                    */
+/************************************************************************/
+
+int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse,
+                        long point_count, int point_offset,
+                        double *x, double *y, double * /*z */ )
+
+{
+    auto hgrids = NS_PROJ::getListOfGridSets(ctx, nadgrids);
+    if( hgrids.empty() )
+    {
+        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
+        return 1;
+    }
+
+
+    for( long i = 0; i < point_count; i++ )
+    {
+        PJ_LP   input;
+
+        long io = i * point_offset;
+        input.phi = y[io];
+        input.lam = x[io];
+
+        auto output = proj_hgrid_apply_internal(ctx, input, inverse ? PJ_INV : PJ_FWD, hgrids);
 
+        if ( output.lam != HUGE_VAL )
+        {
+            y[io] = output.phi;
+            x[io] = output.lam;
+        }
+        else
+        {
+            if( ctx->debug_level >= PJ_LOG_DEBUG_MAJOR )
+            {
+                pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
+                    "pj_apply_gridshift(): failed to find a grid shift table for\n"
+                    "                      location (%.7fdW,%.7fdN)",
+                    x[io] * RAD_TO_DEG,
+                    y[io] * RAD_TO_DEG );
+            }
+        }
+    }
+
+    return 0;
 }
diff --git a/src/apply_vgridshift.cpp b/src/apply_vgridshift.cpp
index daa44858..b0136e5c 100644
--- a/src/apply_vgridshift.cpp
+++ b/src/apply_vgridshift.cpp
@@ -26,7 +26,9 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
-#define PJ_LIB__
+#ifndef FROM_PROJ_CPP
+#define FROM_PROJ_CPP
+#endif
 
 #include <assert.h>
 #include <stdio.h>
@@ -34,273 +36,117 @@
 
 #include <math.h>
 #include "proj_internal.h"
+#include "proj/internal/internal.hpp"
+#include "grids.hpp"
 
-static int is_nodata(float value, double vmultiplier)
-{
-    /* nodata?  */
-    /* GTX official nodata value if  -88.88880f, but some grids also */
-    /* use other  big values for nodata (e.g naptrans2008.gtx has */
-    /* nodata values like -2147479936), so test them too */
-    return value * vmultiplier > 1000 || value * vmultiplier < -1000 || value == -88.88880f;
-}
+NS_PROJ_START
+
+static double read_vgrid_value(const ListOfVGrids& grids, PJ_LP input, double vmultiplier) {
 
-static double read_vgrid_value( PJ *defn, PJ_LP input, double vmultiplier, int *gridlist_count_p, PJ_GRIDINFO **tables, struct CTABLE *ct) {
-    int  itable = 0;
-    double value = HUGE_VAL;
-    double grid_x, grid_y;
-    long   grid_ix, grid_iy;
-    long   grid_ix2, grid_iy2;
-    float  *cvs;
     /* do not deal with NaN coordinates */
     /* cppcheck-suppress duplicateExpression */
     if( isnan(input.phi) || isnan(input.lam) )
-        itable = *gridlist_count_p;
-
-    /* keep trying till we find a table that works */
-    for ( ; itable < *gridlist_count_p; itable++ )
     {
-        PJ_GRIDINFO *gi = tables[itable];
-
-        ct = gi->ct;
-
-        /* skip tables that don't match our point at all (latitude check).  */
-        if( ct->ll.phi > input.phi
-            || ct->ll.phi + (ct->lim.phi-1) * ct->del.phi < input.phi )
-            continue;
-
-        bool fullWorldLongExtent = false;
-        if( fabs(ct->lim.lam * ct->del.lam - 2 * M_PI) < 1e-10 )
-        {
-            fullWorldLongExtent = true;
-        }
-
-        /* skip tables that don't match our point at all (longitude check).  */
-        else if( ct->ll.lam > input.lam
-            || ct->ll.lam + (ct->lim.lam-1) * ct->del.lam < input.lam )
-            continue;
-
-        /* If we have child nodes, check to see if any of them apply. */
-        while( gi->child != nullptr )
-        {
-            PJ_GRIDINFO *child;
-
-            for( child = gi->child; child != nullptr; child = child->next )
-            {
-                struct CTABLE *ct1 = child->ct;
-
-                fullWorldLongExtent = false;
-
-                if( ct1->ll.phi > input.phi
-                  || ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi < input.phi)
-                    continue;
-
-                if( fabs(ct1->lim.lam * ct1->del.lam - 2 * M_PI) < 1e-10 )
-                {
-                    fullWorldLongExtent = true;
-                }
-                else if( ct1->ll.lam > input.lam
-                  || ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam < input.lam)
-                    continue;
-
-                break;
-            }
-
-            /* we didn't find a more refined child node to use, so go with current grid */
-            if( child == nullptr )
-            {
-                break;
-            }
-
-            /* Otherwise let's try for childrens children .. */
-            gi = child;
-            ct = child->ct;
-        }
-
-        /* load the grid shift info if we don't have it. */
-        if( ct->cvs == nullptr )
-        {
-            pj_gridinfo_load( pj_get_ctx(defn), gi );
-        }
-        if( ct->cvs == nullptr )
-        {
-            pj_ctx_set_errno( defn->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-            return PJD_ERR_FAILED_TO_LOAD_GRID;
-        }
+        return HUGE_VAL;
+    }
 
-        /* Interpolation a location within the grid */
-        grid_x = (input.lam - ct->ll.lam) / ct->del.lam;
-        if( fullWorldLongExtent ) {
-            // The first fmod goes to ]-lim, lim[ range
-            // So we add lim again to be in ]0, 2*lim[ and fmod again
-            grid_x = fmod(
-                fmod(grid_x + ct->lim.lam, ct->lim.lam) + ct->lim.lam,
-                ct->lim.lam);
-        }
-        grid_y = (input.phi - ct->ll.phi) / ct->del.phi;
-        grid_ix = lround(floor(grid_x));
-        assert(grid_ix >= 0 && grid_ix < ct->lim.lam);
-        grid_iy = lround(floor(grid_y));
-        assert(grid_iy >= 0 && grid_iy < ct->lim.phi);
-        grid_x -= grid_ix;
-        grid_y -= grid_iy;
+    const VerticalShiftGrid* grid = nullptr;
+    for( const auto& gridset: grids )
+    {
+        grid = gridset->gridAt(input.lam, input.phi);
+        if( grid )
+            break;
+    }
+    if( !grid )
+    {
+        return HUGE_VAL;
+    }
 
-        grid_ix2 = grid_ix + 1;
-        if( grid_ix2 >= ct->lim.lam ) {
-            if( fullWorldLongExtent ) {
-              grid_ix2 = 0;
-            } else {
-              grid_ix2 = ct->lim.lam - 1;
-            }
-        }
-        grid_iy2 = grid_iy + 1;
-        if( grid_iy2 >= ct->lim.phi )
-            grid_iy2 = ct->lim.phi - 1;
+    const auto& extent = grid->extentAndRes();
 
-        cvs = (float *) ct->cvs;
-        {
-            float value_a = cvs[grid_ix + grid_iy * ct->lim.lam];
-            float value_b = cvs[grid_ix2 + grid_iy * ct->lim.lam];
-            float value_c = cvs[grid_ix + grid_iy2 * ct->lim.lam];
-            float value_d = cvs[grid_ix2 + grid_iy2 * ct->lim.lam];
-            double total_weight = 0.0;
-            int n_weights = 0;
-            value = 0.0f;
-            if( !is_nodata(value_a, vmultiplier) )
-            {
-                double weight = (1.0-grid_x) * (1.0-grid_y);
-                value += value_a * weight;
-                total_weight += weight;
-                n_weights ++;
-            }
-            if( !is_nodata(value_b, vmultiplier) )
-            {
-                double weight = (grid_x) * (1.0-grid_y);
-                value += value_b * weight;
-                total_weight += weight;
-                n_weights ++;
-            }
-            if( !is_nodata(value_c, vmultiplier) )
-            {
-                double weight = (1.0-grid_x) * (grid_y);
-                value += value_c * weight;
-                total_weight += weight;
-                n_weights ++;
-            }
-            if( !is_nodata(value_d, vmultiplier) )
-            {
-                double weight = (grid_x) * (grid_y);
-                value += value_d * weight;
-                total_weight += weight;
-                n_weights ++;
-            }
-            if( n_weights == 0 )
-                value = HUGE_VAL;
-            else if( n_weights != 4 )
-                value /= total_weight;
+    /* Interpolation a location within the grid */
+    double grid_x = (input.lam - extent.westLon) / extent.resLon;
+    if( extent.fullWorldLongitude() ) {
+        // The first fmod goes to ]-lim, lim[ range
+        // So we add lim again to be in ]0, 2*lim[ and fmod again
+        grid_x = fmod(
+            fmod(grid_x + grid->width(), grid->width()) + grid->width(),
+            grid->width());
+    }
+    double grid_y = (input.phi - extent.southLat) / extent.resLat;
+    int grid_ix = static_cast<int>(lround(floor(grid_x)));
+    assert(grid_ix >= 0 && grid_ix < grid->width());
+    int grid_iy = static_cast<int>(lround(floor(grid_y)));
+    assert(grid_iy >= 0 && grid_iy < grid->height());
+    grid_x -= grid_ix;
+    grid_y -= grid_iy;
+
+    int grid_ix2 = grid_ix + 1;
+    if( grid_ix2 >= grid->width() ) {
+        if( extent.fullWorldLongitude() ) {
+            grid_ix2 = 0;
+        } else {
+            grid_ix2 = grid->width() - 1;
         }
-
+    }
+    int grid_iy2 = grid_iy + 1;
+    if( grid_iy2 >= grid->height() )
+        grid_iy2 = grid->height() - 1;
+
+    float value_a = 0;
+    float value_b = 0;
+    float value_c = 0;
+    float value_d = 0;
+    if( !grid->valueAt(grid_ix, grid_iy, value_a) ||
+        !grid->valueAt(grid_ix2, grid_iy, value_b) ||
+        !grid->valueAt(grid_ix, grid_iy2, value_c) ||
+        !grid->valueAt(grid_ix2, grid_iy2, value_d) )
+    {
+        return HUGE_VAL;
     }
 
-    return value * vmultiplier;
-}
-
-/************************************************************************/
-/*                        pj_apply_vgridshift()                         */
-/*                                                                      */
-/*      This implementation takes uses the gridlist from a coordinate   */
-/*      system definition.  If the gridlist has not yet been            */
-/*      populated in the coordinate system definition we set it up      */
-/*      now.                                                            */
-/************************************************************************/
-int pj_apply_vgridshift( PJ *defn, const char *listname,
-                         PJ_GRIDINFO ***gridlist_p,
-                         int *gridlist_count_p,
-                         int inverse,
-                         long point_count, int point_offset,
-                         double *x, double *y, double *z )
-
-{
-    int  i;
-    static int debug_count = 0;
-    PJ_GRIDINFO **tables;
-    struct CTABLE ct;
+    double total_weight = 0.0;
+    int n_weights = 0;
+    double value = 0.0f;
 
-    if( *gridlist_p == nullptr )
+    if( !grid->isNodata(value_a, vmultiplier) )
     {
-        *gridlist_p =
-            pj_gridlist_from_nadgrids( pj_get_ctx(defn),
-                                       pj_param(defn->ctx,defn->params,listname).s,
-                                       gridlist_count_p );
-
-        if( *gridlist_p == nullptr || *gridlist_count_p == 0 )
-            return defn->ctx->last_errno;
+        double weight = (1.0-grid_x) * (1.0-grid_y);
+        value += value_a * weight;
+        total_weight += weight;
+        n_weights ++;
     }
-
-    if( *gridlist_count_p == 0 )
+    if( !grid->isNodata(value_b, vmultiplier) )
     {
-        pj_ctx_set_errno( defn->ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
-        return PJD_ERR_FAILED_TO_LOAD_GRID;
+        double weight = (grid_x) * (1.0-grid_y);
+        value += value_b * weight;
+        total_weight += weight;
+        n_weights ++;
     }
-
-    tables = *gridlist_p;
-    defn->ctx->last_errno = 0;
-
-    for( i = 0; i < point_count; i++ )
+    if( !grid->isNodata(value_c, vmultiplier) )
     {
-        double value;
-        long io = i * point_offset;
-        PJ_LP   input;
-
-        input.phi = y[io];
-        input.lam = x[io];
-
-        value = read_vgrid_value(defn, input, 1.0, gridlist_count_p, tables, &ct);
-
-        if( inverse )
-            z[io] -= value;
-        else
-            z[io] += value;
-        if( value != HUGE_VAL )
-        {
-            if( debug_count++ < 20 ) {
-                proj_log_trace(defn, "pj_apply_gridshift(): used %s", ct.id);
-                break;
-            }
-        }
-
-        if( value == HUGE_VAL )
-        {
-            int itable;
-            std::string gridlist;
-
-            proj_log_debug(defn,
-                "pj_apply_vgridshift(): failed to find a grid shift table for\n"
-                "                       location (%.7fdW,%.7fdN)",
-                x[io] * RAD_TO_DEG,
-                y[io] * RAD_TO_DEG );
-
-            for( itable = 0; itable < *gridlist_count_p; itable++ )
-            {
-                PJ_GRIDINFO *gi = tables[itable];
-                if( itable == 0 )
-                    gridlist += "   tried: ";
-                else
-                    gridlist += ',';
-                gridlist += gi->gridname;
-            }
-
-            proj_log_debug(defn, "%s", gridlist.c_str());
-            pj_ctx_set_errno( defn->ctx, PJD_ERR_GRID_AREA );
-
-            return PJD_ERR_GRID_AREA;
-        }
+        double weight = (1.0-grid_x) * (grid_y);
+        value += value_c * weight;
+        total_weight += weight;
+        n_weights ++;
     }
+    if( !grid->isNodata(value_d, vmultiplier) )
+    {
+        double weight = (grid_x) * (grid_y);
+        value += value_d * weight;
+        total_weight += weight;
+        n_weights ++;
+    }
+    if( n_weights == 0 )
+        value = HUGE_VAL;
+    else if( n_weights != 4 )
+        value /= total_weight;
 
-    return 0;
+    return value * vmultiplier;
 }
 
 /**********************************************/
-int proj_vgrid_init(PJ* P, const char *grids) {
+ListOfVGrids proj_vgrid_init(PJ* P, const char *gridkey) {
 /**********************************************
 
   Initizalize and populate gridlist.
@@ -314,33 +160,43 @@ int proj_vgrid_init(PJ* P, const char *grids) {
 
 ***********************************************/
 
-    /* prepend "s" to the "grids" string to allow usage with pj_param */
-    char *sgrids = (char *) pj_malloc( (strlen(grids)+1+1) *sizeof(char) );
-    sprintf(sgrids, "%s%s", "s", grids);
+    std::string key("s");
+    key += gridkey;
+    const char* gridnames = pj_param(P->ctx, P->params, key.c_str()).s;
+    if( gridnames == nullptr )
+        return {};
 
-    if (P->vgridlist_geoid == nullptr) {
-        P->vgridlist_geoid = pj_gridlist_from_nadgrids(
-            P->ctx,
-            pj_param(P->ctx, P->params, sgrids).s,
-            &(P->vgridlist_geoid_count)
-        );
-
-        if( P->vgridlist_geoid == nullptr || P->vgridlist_geoid_count == 0 ) {
-            pj_dealloc(sgrids);
-            return 0;
+    auto listOfGridNames = internal::split(std::string(gridnames), ',');
+    ListOfVGrids grids;
+    for( const auto& gridnameStr: listOfGridNames )
+    {
+        const char* gridname = gridnameStr.c_str();
+        bool canFail = false;
+        if( gridname[0] == '@' )
+        {
+            canFail = true;
+            gridname ++;
+        }
+        auto gridSet = VerticalShiftGridSet::open(P->ctx, gridname);
+        if( !gridSet )
+        {
+            if( !canFail )
+            {
+                pj_ctx_set_errno( P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
+                return {};
+            }
+        }
+        else
+        {
+            grids.emplace_back(std::move(gridSet));
         }
     }
 
-    if (P->vgridlist_geoid_count == 0) {
-        proj_errno_set(P, PJD_ERR_FAILED_TO_LOAD_GRID);
-    }
-
-    pj_dealloc(sgrids);
-    return P->vgridlist_geoid_count;
+    return grids;
 }
 
 /***********************************************/
-double proj_vgrid_value(PJ *P, PJ_LP lp, double vmultiplier){
+double proj_vgrid_value(PJ *P, const ListOfVGrids& grids, PJ_LP lp, double vmultiplier){
 /***********************************************
 
   Read grid value at position lp in grids loaded
@@ -350,12 +206,12 @@ double proj_vgrid_value(PJ *P, PJ_LP lp, double vmultiplier){
 
 ************************************************/
 
-    struct CTABLE used_grid;
     double value;
-    memset(&used_grid, 0, sizeof(struct CTABLE));
 
-    value = read_vgrid_value(P, lp, vmultiplier, &(P->vgridlist_geoid_count), P->vgridlist_geoid, &used_grid);
+    value = read_vgrid_value(grids, lp, vmultiplier);
     proj_log_trace(P, "proj_vgrid_value: (%f, %f) = %f", lp.lam*RAD_TO_DEG, lp.phi*RAD_TO_DEG, value);
 
     return value;
 }
+
+NS_PROJ_END
diff --git a/src/datum_set.cpp b/src/datum_set.cpp
index 873d7be5..15d51613 100644
--- a/src/datum_set.cpp
+++ b/src/datum_set.cpp
@@ -41,7 +41,7 @@
 int pj_datum_set(projCtx ctx, paralist *pl, PJ *projdef)
 
 {
-    const char *name, *towgs84, *nadgrids, *catalog;
+    const char *name, *towgs84, *nadgrids;
 
     projdef->datum_type = PJD_UNKNOWN;
 
@@ -118,25 +118,6 @@ int pj_datum_set(projCtx ctx, paralist *pl, PJ *projdef)
     }
 
 /* -------------------------------------------------------------------- */
-/*      Check for grid catalog parameter, and optional date.            */
-/* -------------------------------------------------------------------- */
-    else if( (catalog = pj_param(ctx, pl,"scatalog").s) != nullptr )
-    {
-        const char *date;
-
-        projdef->datum_type = PJD_GRIDSHIFT;
-        projdef->catalog_name = pj_strdup(catalog);
-        if (!projdef->catalog_name) {
-            pj_ctx_set_errno(ctx, ENOMEM);
-            return 1;
-        }
-
-        date = pj_param(ctx, pl, "sdate").s;
-        if( date != nullptr) 
-            projdef->datum_date = pj_gc_parsedate( ctx, date);
-    }
-
-/* -------------------------------------------------------------------- */
 /*      Check for towgs84 parameter.                                    */
 /* -------------------------------------------------------------------- */
     else if( (towgs84 = pj_param(ctx, pl,"stowgs84").s) != nullptr )
diff --git a/src/gc_reader.cpp b/src/gc_reader.cpp
deleted file mode 100644
index def52a11..00000000
--- a/src/gc_reader.cpp
+++ /dev/null
@@ -1,248 +0,0 @@
-/******************************************************************************
- * Project:  PROJ.4
- * Purpose:  Code to read a grid catalog from a .cvs file.
- * Author:   Frank Warmerdam, warmerdam@pobox.com
- *
- ******************************************************************************
- * Copyright (c) 2012, Frank Warmerdam <warmerdam@pobox.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 <ctype.h>
-#include <errno.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include "proj.h"
-#include "proj_internal.h"
-
-static int gc_readentry(projCtx ctx, PAFile fid, PJ_GridCatalogEntry *entry);
-
-/************************************************************************/
-/*                         pj_gc_readcatalog()                          */
-/*                                                                      */
-/*      Read a grid catalog from a .csv file.                           */
-/************************************************************************/
-
-PJ_GridCatalog *pj_gc_readcatalog( projCtx ctx, const char *catalog_name )
-{
-    PAFile fid;
-    PJ_GridCatalog *catalog;
-    int entry_max;
-    char line[302];
-
-    fid = pj_open_lib( ctx, catalog_name, "r" );
-    if (fid == nullptr) 
-        return nullptr;
-
-    /* discard title line */
-    pj_ctx_fgets(ctx, line, sizeof(line)-1, fid);
-
-    catalog = (PJ_GridCatalog *) calloc(1,sizeof(PJ_GridCatalog));
-    if( !catalog )
-    {
-        pj_ctx_set_errno(ctx, ENOMEM);
-        pj_ctx_fclose(ctx, fid);
-        return nullptr;
-    }
-    
-    catalog->catalog_name = pj_strdup(catalog_name);
-    if (!catalog->catalog_name) {
-        pj_ctx_set_errno(ctx, ENOMEM);
-        free(catalog);
-        pj_ctx_fclose(ctx, fid);
-        return nullptr;
-    }
-    
-    entry_max = 10;
-    catalog->entries = (PJ_GridCatalogEntry *) 
-        malloc(entry_max * sizeof(PJ_GridCatalogEntry));
-    if (!catalog->entries) {
-        pj_ctx_set_errno(ctx, ENOMEM);
-        free(catalog->catalog_name);
-        free(catalog);
-        pj_ctx_fclose(ctx, fid);
-        return nullptr;
-    }
-    
-    while( gc_readentry( ctx, fid, 
-                            catalog->entries+catalog->entry_count) == 0)
-    {
-        catalog->entry_count++;
-        
-        if( catalog->entry_count == entry_max ) 
-        {
-            PJ_GridCatalogEntry* new_entries;
-            entry_max = entry_max * 2;
-            new_entries = (PJ_GridCatalogEntry *) 
-                realloc(catalog->entries, 
-                        entry_max * sizeof(PJ_GridCatalogEntry));
-            if (new_entries == nullptr )
-            {
-                int i;
-                for( i = 0; i < catalog->entry_count; i++ )
-                    free( catalog->entries[i].definition );
-                free( catalog->entries );
-                free( catalog->catalog_name );
-                free( catalog );
-                pj_ctx_fclose(ctx, fid);
-                return nullptr;
-            }
-            catalog->entries = new_entries;
-        }
-    }
-
-    pj_ctx_fclose(ctx, fid);
-
-    return catalog;
-}
-
-/************************************************************************/
-/*                        gc_read_csv_line()                         */
-/*                                                                      */
-/*      Simple csv line splitter with fixed maximum line size and       */
-/*      token count.                                                    */
-/************************************************************************/
-
-static int gc_read_csv_line( projCtx ctx, PAFile fid, 
-                                char **tokens, int max_tokens ) 
-{
-    char line[302];
-   
-    while( pj_ctx_fgets(ctx, line, sizeof(line)-1, fid) != nullptr )
-    {
-        char *next = line;
-        int token_count = 0;
-        
-        while( isspace(*next) ) 
-            next++;
-        
-        /* skip blank and comment lines */
-        if( next[0] == '#' || next[0] == '\0' )
-            continue;
-        
-        while( token_count < max_tokens && *next != '\0' ) 
-        {
-            const char *start = next;
-            char* token;
-            
-            while( *next != '\0' && *next != ',' ) 
-                next++;
-            
-            if( *next == ',' )
-            {
-                *next = '\0';
-                next++;
-            }
-            
-            token = pj_strdup(start);
-            if (!token) {
-                while (token_count > 0)
-                    free(tokens[--token_count]);
-                pj_ctx_set_errno(ctx, ENOMEM);
-                return 0;
-            }
-            tokens[token_count++] = token;
-        }
-
-        return token_count;
-    }
-    
-    return 0; 
-}
-
-/************************************************************************/
-/*                          pj_gc_parsedate()                           */
-/*                                                                      */
-/*      Parse a date into a floating point year value.  Acceptable      */
-/*      values are "yyyy.fraction" and "yyyy-mm-dd".  Anything else     */
-/*      returns 0.0.                                                    */
-/************************************************************************/
-
-double pj_gc_parsedate( projCtx ctx, const char *date_string )
-{
-    (void) ctx;
-
-    if( strlen(date_string) == 10 
-        && date_string[4] == '-' && date_string[7] == '-' ) 
-    {
-        int year = atoi(date_string);
-        int month = atoi(date_string+5);
-        int day = atoi(date_string+8);
-
-        /* simplified calculation so we don't need to know all about months */
-        return year + ((month-1) * 31 + (day-1)) / 372.0;
-    }
-    else 
-    {
-        return pj_atof(date_string);
-    }
-}
-
-
-/************************************************************************/
-/*                          gc_readentry()                           */
-/*                                                                      */
-/*      Read one catalog entry from the file                            */
-/*                                                                      */
-/*      Format:                                                         */
-/*        gridname,ll_long,ll_lat,ur_long,ur_lat,priority,date          */
-/************************************************************************/
-
-static int gc_readentry(projCtx ctx, PAFile fid, PJ_GridCatalogEntry *entry) 
-{
-#define MAX_TOKENS 30
-    char *tokens[MAX_TOKENS];
-    int token_count, i;
-    int error = 0;
-
-    memset( entry, 0, sizeof(PJ_GridCatalogEntry) );
-    
-    token_count = gc_read_csv_line( ctx, fid, tokens, MAX_TOKENS );
-    if( token_count < 5 )
-    {
-        error = 1; /* TODO: need real error codes */
-        if( token_count != 0 )
-            pj_log( ctx, PJ_LOG_ERROR, "Short line in grid catalog." );
-    }
-    else
-    {
-        entry->definition = tokens[0];
-        tokens[0] = nullptr;   /* We take ownership of tokens[0] */
-        entry->region.ll_long = dmstor_ctx( ctx, tokens[1], nullptr );
-        entry->region.ll_lat = dmstor_ctx( ctx, tokens[2], nullptr );
-        entry->region.ur_long = dmstor_ctx( ctx, tokens[3], nullptr );
-        entry->region.ur_lat = dmstor_ctx( ctx, tokens[4], nullptr );
-        if( token_count > 5 )
-            entry->priority = atoi( tokens[5] ); /* defaults to zero */
-        if( token_count > 6 )
-            entry->date = pj_gc_parsedate( ctx, tokens[6] );
-    }
-
-    for( i = 0; i < token_count; i++ )
-        free( tokens[i] );
-
-    return error;
-}
-
-
-
diff --git a/src/gridcatalog.cpp b/src/gridcatalog.cpp
deleted file mode 100644
index 15d81dd7..00000000
--- a/src/gridcatalog.cpp
+++ /dev/null
@@ -1,306 +0,0 @@
-/******************************************************************************
- * Project:  PROJ.4
- * Purpose:  Code in support of grid catalogs
- * Author:   Frank Warmerdam, warmerdam@pobox.com
- *
- ******************************************************************************
- * Copyright (c) 2012, Frank Warmerdam <warmerdam@pobox.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 <assert.h>
-#include <math.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include "proj.h"
-#include "proj_internal.h"
-
-static PJ_GridCatalog *grid_catalog_list = nullptr;
-
-static
-PJ_GRIDINFO *pj_gc_findgrid( projCtx_t *ctx,
-                             PJ_GridCatalog *catalog, int after,
-                             PJ_LP location, double date,
-                             PJ_Region *optional_region,
-                             double *grid_date );
-
-/************************************************************************/
-/*                          pj_gc_unloadall()                           */
-/*                                                                      */
-/*      Deallocate all the grid catalogs (but not the referenced        */
-/*      grids).                                                         */
-/************************************************************************/
-
-void pj_gc_unloadall( projCtx ctx )
-{
-    (void) ctx;
-
-    while( grid_catalog_list != nullptr )
-    {
-        int i;
-        PJ_GridCatalog *catalog = grid_catalog_list;
-        grid_catalog_list = grid_catalog_list->next;
-
-        for( i = 0; i < catalog->entry_count; i++ )
-        {
-            /* we don't own gridinfo - do not free here */
-            free( catalog->entries[i].definition );
-        }
-        free( catalog->entries );
-        free( catalog->catalog_name );
-        free( catalog );
-    }
-}
-
-/************************************************************************/
-/*                         pj_gc_findcatalog()                          */
-/************************************************************************/
-
-PJ_GridCatalog *pj_gc_findcatalog( projCtx ctx, const char *name )
-
-{
-    PJ_GridCatalog *catalog;
-
-    pj_acquire_lock();
-
-    for( catalog=grid_catalog_list; catalog != nullptr; catalog = catalog->next ) 
-    {
-        if( strcmp(catalog->catalog_name, name) == 0 )
-        {
-            pj_release_lock();
-            return catalog;
-        }
-    }
-
-    pj_release_lock();
-
-    catalog = pj_gc_readcatalog( ctx, name );
-    if( catalog == nullptr )
-        return nullptr;
-
-    pj_acquire_lock();
-    catalog->next = grid_catalog_list;
-    grid_catalog_list = catalog;
-    pj_release_lock();
-
-    return catalog;
-}
-
-/************************************************************************/
-/*                       pj_gc_apply_gridshift()                        */
-/************************************************************************/
-
-int pj_gc_apply_gridshift( PJ *defn, int inverse, 
-                           long point_count, int point_offset, 
-                           double *x, double *y, double *z )
-
-{
-    int i;
-    (void) z;
-
-    if( defn->catalog == nullptr ) 
-    {
-        defn->catalog = pj_gc_findcatalog( defn->ctx, defn->catalog_name );
-        if( defn->catalog == nullptr )
-            return defn->ctx->last_errno;
-    }
-
-    defn->ctx->last_errno = 0;
-
-    for( i = 0; i < point_count; i++ )
-    {
-        long io = i * point_offset;
-        PJ_LP   input, output_after, output_before;
-        double mix_ratio;
-        PJ_GRIDINFO *gi;
-
-        input.phi = y[io];
-        input.lam = x[io];
-
-        /* make sure we have appropriate "after" shift file available */
-        if( defn->last_after_grid == nullptr
-            || input.lam < defn->last_after_region.ll_long
-            || input.lam > defn->last_after_region.ur_long
-            || input.phi < defn->last_after_region.ll_lat
-            || input.phi > defn->last_after_region.ll_lat ) {
-            defn->last_after_grid = 
-                pj_gc_findgrid( defn->ctx, defn->catalog, 
-                                1, input, defn->datum_date, 
-                                &(defn->last_after_region), 
-                                &(defn->last_after_date));
-            if( defn->last_after_grid == nullptr )
-            {
-                pj_ctx_set_errno( defn->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-                return PJD_ERR_FAILED_TO_LOAD_GRID;
-            }
-        }
-        gi = defn->last_after_grid;
-        assert( gi->child == nullptr );
-
-        /* load the grid shift info if we don't have it. */
-        if( gi->ct->cvs == nullptr && !pj_gridinfo_load( defn->ctx, gi ) )
-        {
-            pj_ctx_set_errno( defn->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-            return PJD_ERR_FAILED_TO_LOAD_GRID;
-        }
-            
-        output_after = nad_cvt( input, inverse, gi->ct );
-        if( output_after.lam == HUGE_VAL )
-        {
-            if( defn->ctx->debug_level >= PJ_LOG_DEBUG_MAJOR )
-            {
-                pj_log( defn->ctx, PJ_LOG_DEBUG_MAJOR,
-                        "pj_apply_gridshift(): failed to find a grid shift table for\n"
-                        "                      location (%.7fdW,%.7fdN)",
-                        x[io] * RAD_TO_DEG, 
-                        y[io] * RAD_TO_DEG );
-            }
-            continue;
-        }
-
-        if( defn->datum_date == 0.0 ) 
-        {
-            y[io] = output_after.phi;
-            x[io] = output_after.lam;
-            continue;
-        }
-
-        /* make sure we have appropriate "before" shift file available */
-        if( defn->last_before_grid == nullptr
-            || input.lam < defn->last_before_region.ll_long
-            || input.lam > defn->last_before_region.ur_long
-            || input.phi < defn->last_before_region.ll_lat
-            || input.phi > defn->last_before_region.ll_lat ) {
-            defn->last_before_grid = 
-                pj_gc_findgrid( defn->ctx, defn->catalog, 
-                                0, input, defn->datum_date, 
-                                &(defn->last_before_region), 
-                                &(defn->last_before_date));
-            if( defn->last_before_grid == nullptr )
-            {
-                pj_ctx_set_errno( defn->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-                return PJD_ERR_FAILED_TO_LOAD_GRID;
-            }
-        }
-
-        gi = defn->last_before_grid;
-        assert( gi->child == nullptr );
-
-        /* load the grid shift info if we don't have it. */
-        if( gi->ct->cvs == nullptr && !pj_gridinfo_load( defn->ctx, gi ) )
-        {
-            pj_ctx_set_errno( defn->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-            return PJD_ERR_FAILED_TO_LOAD_GRID;
-        }
-            
-        output_before = nad_cvt( input, inverse, gi->ct );
-        if( output_before.lam == HUGE_VAL )
-        {
-            if( defn->ctx->debug_level >= PJ_LOG_DEBUG_MAJOR )
-            {
-                pj_log( defn->ctx, PJ_LOG_DEBUG_MAJOR,
-                        "pj_apply_gridshift(): failed to find a grid shift table for\n"
-                        "                      location (%.7fdW,%.7fdN)",
-                        x[io] * RAD_TO_DEG, 
-                        y[io] * RAD_TO_DEG );
-            }
-            continue;
-        }
-
-        mix_ratio = (defn->datum_date - defn->last_before_date) 
-            / (defn->last_after_date - defn->last_before_date);
-
-        y[io] = mix_ratio * output_after.phi 
-            + (1.0-mix_ratio) * output_before.phi;
-        x[io] = mix_ratio * output_after.lam 
-            + (1.0-mix_ratio) * output_before.lam;
-    }
-
-    return 0;
-}
-
-/************************************************************************/
-/*                           pj_c_findgrid()                            */
-/************************************************************************/
-
-static
-PJ_GRIDINFO *pj_gc_findgrid( projCtx ctx, PJ_GridCatalog *catalog, int after,
-                             PJ_LP location, double date,
-                             PJ_Region *optional_region,
-                             double *grid_date ) 
-{
-    int iEntry;
-    PJ_GridCatalogEntry *entry = nullptr;
-
-    for( iEntry = 0; iEntry < catalog->entry_count; iEntry++ ) 
-    {
-        entry = catalog->entries + iEntry;
-
-        if( (after && entry->date < date) 
-            || (!after && entry->date > date) )
-            continue;
-
-        if( location.lam < entry->region.ll_long
-            || location.lam > entry->region.ur_long
-            || location.phi < entry->region.ll_lat
-            || location.phi > entry->region.ur_lat )
-            continue;
-
-        if( entry->available == -1 )
-            continue;
-
-        break;
-    }
-
-    if( entry == nullptr )
-    {
-        if( grid_date )
-            *grid_date = 0.0;
-        if( optional_region != nullptr )
-            memset( optional_region, 0, sizeof(PJ_Region));
-        return nullptr;
-    }
-
-    if( grid_date )
-        *grid_date = entry->date;
-
-    if( optional_region )
-    {
-        
-    }
-
-    if( entry->gridinfo == nullptr )
-    {
-        PJ_GRIDINFO **gridlist = nullptr;
-        int grid_count = 0;
-        gridlist = pj_gridlist_from_nadgrids( ctx, entry->definition, 
-                                              &grid_count);
-        // FIXME: this leaks gridlist itself, and memory ownership of
-        // entry->gridinfo is also confusing. Coverity CID 193539
-        if( grid_count == 1 )
-            entry->gridinfo = gridlist[0];
-    }
-    
-    return entry->gridinfo;
-}
-                             
diff --git a/src/gridinfo.cpp b/src/gridinfo.cpp
deleted file mode 100644
index 7f7d930f..00000000
--- a/src/gridinfo.cpp
+++ /dev/null
@@ -1,988 +0,0 @@
-/******************************************************************************
- * Project:  PROJ.4
- * Purpose:  Functions for handling individual PJ_GRIDINFO's.  Includes
- *           loaders for all formats but CTABLE (in nad_init.c).
- * Author:   Frank Warmerdam, warmerdam@pobox.com
- *
- ******************************************************************************
- * Copyright (c) 2000, Frank Warmerdam <warmerdam@pobox.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 <stddef.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include "proj_internal.h"
-
-/************************************************************************/
-/*                             swap_words()                             */
-/*                                                                      */
-/*      Convert the byte order of the given word(s) in place.           */
-/************************************************************************/
-
-static const int  byte_order_test = 1;
-#define IS_LSB	(1 == ((const unsigned char *) (&byte_order_test))[0])
-
-static void swap_words( unsigned char *data, int word_size, int word_count )
-
-{
-    int	word;
-
-    for( word = 0; word < word_count; word++ )
-    {
-        int	i;
-
-        for( i = 0; i < word_size/2; i++ )
-        {
-            unsigned char	t;
-
-            t = data[i];
-            data[i] = data[word_size-i-1];
-            data[word_size-i-1] = t;
-        }
-
-        data += word_size;
-    }
-}
-
-/************************************************************************/
-/*                             to_double()                              */
-/*                                                                      */
-/*      Returns a double from the pointed data.                         */
-/************************************************************************/
-
-static double to_double( unsigned char* data )
-{
-    double d;
-    memcpy(&d, data, sizeof(d));
-    return d;
-}
-
-/************************************************************************/
-/*                          pj_gridinfo_free()                          */
-/************************************************************************/
-
-void pj_gridinfo_free( projCtx ctx, PJ_GRIDINFO *gi )
-
-{
-    if( gi == nullptr )
-        return;
-
-    if( gi->child != nullptr )
-    {
-        PJ_GRIDINFO *child, *next;
-
-        for( child = gi->child; child != nullptr; child=next)
-        {
-            next=child->next;
-            pj_gridinfo_free( ctx, child );
-        }
-    }
-
-    if( gi->ct != nullptr )
-        nad_free( gi->ct );
-
-    free( gi->gridname );
-    if( gi->filename != nullptr )
-        free( gi->filename );
-
-    pj_dalloc( gi );
-}
-
-/************************************************************************/
-/*                          pj_gridinfo_load()                          */
-/*                                                                      */
-/*      This function is intended to implement delayed loading of       */
-/*      the data contents of a grid file.  The header and related       */
-/*      stuff are loaded by pj_gridinfo_init().                         */
-/************************************************************************/
-
-int pj_gridinfo_load( projCtx_t* ctx, PJ_GRIDINFO *gi )
-
-{
-    struct CTABLE ct_tmp;
-
-    if( gi == nullptr || gi->ct == nullptr )
-        return 0;
-
-    pj_acquire_lock();
-    if( gi->ct->cvs != nullptr )
-    {
-        pj_release_lock();
-        return 1;
-    }
-
-    memcpy(&ct_tmp, gi->ct, sizeof(struct CTABLE));
-
-/* -------------------------------------------------------------------- */
-/*      Original platform specific CTable format.                       */
-/* -------------------------------------------------------------------- */
-    if( strcmp(gi->format,"ctable") == 0 )
-    {
-        PAFile fid;
-        int result;
-
-        fid = pj_open_lib( ctx, gi->filename, "rb" );
-
-        if( fid == nullptr )
-        {
-            pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-            pj_release_lock();
-            return 0;
-        }
-
-        result = nad_ctable_load( ctx, &ct_tmp, (struct projFileAPI_t*)fid );
-
-        pj_ctx_fclose( ctx, fid );
-
-        gi->ct->cvs = ct_tmp.cvs;
-        pj_release_lock();
-
-        return result;
-    }
-
-/* -------------------------------------------------------------------- */
-/*      CTable2 format.                                                 */
-/* -------------------------------------------------------------------- */
-    else if( strcmp(gi->format,"ctable2") == 0 )
-    {
-        PAFile fid;
-        int result;
-
-        fid = pj_open_lib( ctx, gi->filename, "rb" );
-
-        if( fid == nullptr )
-        {
-            pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-            pj_release_lock();
-            return 0;
-        }
-
-        result = nad_ctable2_load( ctx, &ct_tmp, (struct projFileAPI_t*)fid );
-
-        pj_ctx_fclose( ctx, fid );
-
-        gi->ct->cvs = ct_tmp.cvs;
-
-        pj_release_lock();
-        return result;
-    }
-
-/* -------------------------------------------------------------------- */
-/*      NTv1 format.                                                    */
-/*      We process one line at a time.  Note that the array storage     */
-/*      direction (e-w) is different in the NTv1 file and what          */
-/*      the CTABLE is supposed to have.  The phi/lam are also           */
-/*      reversed, and we have to be aware of byte swapping.             */
-/* -------------------------------------------------------------------- */
-    else if( strcmp(gi->format,"ntv1") == 0 )
-    {
-        double	*row_buf;
-        int	row;
-        PAFile fid;
-
-        fid = pj_open_lib( ctx, gi->filename, "rb" );
-
-        if( fid == nullptr )
-        {
-            pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-            pj_release_lock();
-            return 0;
-        }
-
-        pj_ctx_fseek( ctx, fid, gi->grid_offset, SEEK_SET );
-
-        row_buf = (double *) pj_malloc(gi->ct->lim.lam * sizeof(double) * 2);
-        ct_tmp.cvs = (FLP *) pj_malloc(gi->ct->lim.lam*gi->ct->lim.phi*sizeof(FLP));
-        if( row_buf == nullptr || ct_tmp.cvs == nullptr )
-        {
-            pj_dalloc( row_buf );
-            pj_dalloc( ct_tmp.cvs );
-            pj_ctx_set_errno( ctx, ENOMEM );
-            pj_release_lock();
-            return 0;
-        }
-
-        for( row = 0; row < gi->ct->lim.phi; row++ )
-        {
-            int	    i;
-            FLP     *cvs;
-            double  *diff_seconds;
-
-            if( pj_ctx_fread( ctx, row_buf,
-                              sizeof(double), gi->ct->lim.lam * 2, fid )
-                != (size_t)( 2 * gi->ct->lim.lam ) )
-            {
-                pj_dalloc( row_buf );
-                pj_dalloc( ct_tmp.cvs );
-                pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-                pj_release_lock();
-                return 0;
-            }
-
-            if( IS_LSB )
-                swap_words( (unsigned char *) row_buf, 8, gi->ct->lim.lam*2 );
-
-            /* convert seconds to radians */
-            diff_seconds = row_buf;
-
-            for( i = 0; i < gi->ct->lim.lam; i++ )
-            {
-                cvs = ct_tmp.cvs + (row) * gi->ct->lim.lam
-                    + (gi->ct->lim.lam - i - 1);
-
-                cvs->phi = (float) (*(diff_seconds++) * ((M_PI/180.0) / 3600.0));
-                cvs->lam = (float) (*(diff_seconds++) * ((M_PI/180.0) / 3600.0));
-            }
-        }
-
-        pj_dalloc( row_buf );
-
-        pj_ctx_fclose( ctx, fid );
-
-        gi->ct->cvs = ct_tmp.cvs;
-        pj_release_lock();
-
-        return 1;
-    }
-
-/* -------------------------------------------------------------------- */
-/*      NTv2 format.                                                    */
-/*      We process one line at a time.  Note that the array storage     */
-/*      direction (e-w) is different in the NTv2 file and what          */
-/*      the CTABLE is supposed to have.  The phi/lam are also           */
-/*      reversed, and we have to be aware of byte swapping.             */
-/* -------------------------------------------------------------------- */
-    else if( strcmp(gi->format,"ntv2") == 0 )
-    {
-        float	*row_buf;
-        int	row;
-        PAFile fid;
-
-        pj_log( ctx, PJ_LOG_DEBUG_MINOR,
-                "NTv2 - loading grid %s", gi->ct->id );
-
-        fid = pj_open_lib( ctx, gi->filename, "rb" );
-
-        if( fid == nullptr )
-        {
-            pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-            pj_release_lock();
-            return 0;
-        }
-
-        pj_ctx_fseek( ctx, fid, gi->grid_offset, SEEK_SET );
-
-        row_buf = (float *) pj_malloc(gi->ct->lim.lam * sizeof(float) * 4);
-        ct_tmp.cvs = (FLP *) pj_malloc(gi->ct->lim.lam*gi->ct->lim.phi*sizeof(FLP));
-        if( row_buf == nullptr || ct_tmp.cvs == nullptr )
-        {
-            pj_dalloc( row_buf );
-            pj_dalloc( ct_tmp.cvs );
-            pj_ctx_set_errno( ctx, ENOMEM );
-            pj_release_lock();
-            return 0;
-        }
-
-        for( row = 0; row < gi->ct->lim.phi; row++ )
-        {
-            int	    i;
-            FLP     *cvs;
-            float   *diff_seconds;
-
-            if( pj_ctx_fread( ctx, row_buf, sizeof(float),
-                              gi->ct->lim.lam*4, fid )
-                != (size_t)( 4 * gi->ct->lim.lam ) )
-            {
-                pj_dalloc( row_buf );
-                pj_dalloc( ct_tmp.cvs );
-                pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-                pj_release_lock();
-                return 0;
-            }
-
-            if( gi->must_swap )
-                swap_words( (unsigned char *) row_buf, 4,
-                            gi->ct->lim.lam*4 );
-
-            /* convert seconds to radians */
-            diff_seconds = row_buf;
-
-            for( i = 0; i < gi->ct->lim.lam; i++ )
-            {
-                cvs = ct_tmp.cvs + (row) * gi->ct->lim.lam
-                    + (gi->ct->lim.lam - i - 1);
-
-                cvs->phi = (float) (*(diff_seconds++) * ((M_PI/180.0) / 3600.0));
-                cvs->lam = (float) (*(diff_seconds++) * ((M_PI/180.0) / 3600.0));
-                diff_seconds += 2; /* skip accuracy values */
-            }
-        }
-
-        pj_dalloc( row_buf );
-
-        pj_ctx_fclose( ctx, fid );
-
-        gi->ct->cvs = ct_tmp.cvs;
-
-        pj_release_lock();
-        return 1;
-    }
-
-/* -------------------------------------------------------------------- */
-/*      GTX format.                                                     */
-/* -------------------------------------------------------------------- */
-    else if( strcmp(gi->format,"gtx") == 0 )
-    {
-        int   words = gi->ct->lim.lam * gi->ct->lim.phi;
-        PAFile fid;
-
-        fid = pj_open_lib( ctx, gi->filename, "rb" );
-
-        if( fid == nullptr )
-        {
-            pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-            pj_release_lock();
-            return 0;
-        }
-
-        pj_ctx_fseek( ctx, fid, gi->grid_offset, SEEK_SET );
-
-        ct_tmp.cvs = (FLP *) pj_malloc(words*sizeof(float));
-        if( ct_tmp.cvs == nullptr )
-        {
-            pj_ctx_set_errno( ctx, ENOMEM );
-            pj_release_lock();
-            return 0;
-        }
-
-        if( pj_ctx_fread( ctx, ct_tmp.cvs, sizeof(float), words, fid )
-            != (size_t)words )
-        {
-            pj_dalloc( ct_tmp.cvs );
-            pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-            pj_release_lock();
-            return 0;
-        }
-
-        if( IS_LSB )
-            swap_words( (unsigned char *) ct_tmp.cvs, 4, words );
-
-        pj_ctx_fclose( ctx, fid );
-        gi->ct->cvs = ct_tmp.cvs;
-        pj_release_lock();
-        return 1;
-    }
-
-    else
-    {
-        pj_release_lock();
-        return 0;
-    }
-}
-
-/************************************************************************/
-/*                        gridinfo_parent()                          */
-/*                                                                      */
-/*      Seek a parent grid file by name from a grid list                */
-/************************************************************************/
-
-static PJ_GRIDINFO* gridinfo_parent( PJ_GRIDINFO *gilist,
-        const char *name, int length )
-{
-    while( gilist )
-    {
-        if( strncmp(gilist->ct->id,name,length) == 0 ) return gilist;
-        if( gilist->child )
-        {
-            PJ_GRIDINFO *parent=gridinfo_parent( gilist->child, name, length );
-            if( parent ) return parent;
-        }
-        gilist=gilist->next;
-    }
-    return gilist;
-}
-
-/************************************************************************/
-/*                       pj_gridinfo_init_ntv2()                        */
-/*                                                                      */
-/*      Load a ntv2 (.gsb) file.                                        */
-/************************************************************************/
-
-static int pj_gridinfo_init_ntv2( projCtx ctx, PAFile fid, PJ_GRIDINFO *gilist )
-{
-    unsigned char header[11*16];
-    int num_subfiles, subfile;
-    int must_swap;
-
-    /* cppcheck-suppress sizeofCalculation */
-    STATIC_ASSERT( sizeof(pj_int32) == 4 );
-    /* cppcheck-suppress sizeofCalculation */
-    STATIC_ASSERT( sizeof(double) == 8 );
-
-/* -------------------------------------------------------------------- */
-/*      Read the overview header.                                       */
-/* -------------------------------------------------------------------- */
-    if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 )
-    {
-        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        return 0;
-    }
-
-    if( header[8] == 11 )
-        must_swap = !IS_LSB;
-    else
-        must_swap = IS_LSB;
-
-/* -------------------------------------------------------------------- */
-/*      Byte swap interesting fields if needed.                         */
-/* -------------------------------------------------------------------- */
-    if( must_swap )
-    {
-        swap_words( header+8, 4, 1 );
-        swap_words( header+8+16, 4, 1 );
-        swap_words( header+8+32, 4, 1 );
-        swap_words( header+8+7*16, 8, 1 );
-        swap_words( header+8+8*16, 8, 1 );
-        swap_words( header+8+9*16, 8, 1 );
-        swap_words( header+8+10*16, 8, 1 );
-    }
-
-/* -------------------------------------------------------------------- */
-/*      Get the subfile count out ... all we really use for now.        */
-/* -------------------------------------------------------------------- */
-    memcpy( &num_subfiles, header+8+32, 4 );
-
-/* ==================================================================== */
-/*      Step through the subfiles, creating a PJ_GRIDINFO for each.     */
-/* ==================================================================== */
-    for( subfile = 0; subfile < num_subfiles; subfile++ )
-    {
-        struct CTABLE *ct;
-        PJ_LP ur;
-        int gs_count;
-        PJ_GRIDINFO *gi;
-
-/* -------------------------------------------------------------------- */
-/*      Read header.                                                    */
-/* -------------------------------------------------------------------- */
-        if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 )
-        {
-            pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-            return 0;
-        }
-
-        if( strncmp((const char *) header,"SUB_NAME",8) != 0 )
-        {
-            pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-            return 0;
-        }
-
-/* -------------------------------------------------------------------- */
-/*      Byte swap interesting fields if needed.                         */
-/* -------------------------------------------------------------------- */
-        if( must_swap )
-        {
-            swap_words( header+8+16*4, 8, 1 );
-            swap_words( header+8+16*5, 8, 1 );
-            swap_words( header+8+16*6, 8, 1 );
-            swap_words( header+8+16*7, 8, 1 );
-            swap_words( header+8+16*8, 8, 1 );
-            swap_words( header+8+16*9, 8, 1 );
-            swap_words( header+8+16*10, 4, 1 );
-        }
-
-/* -------------------------------------------------------------------- */
-/*      Initialize a corresponding "ct" structure.                      */
-/* -------------------------------------------------------------------- */
-        ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE));
-        if (!ct) {
-            pj_ctx_set_errno(ctx, ENOMEM);
-            return 0;
-        }
-        strncpy( ct->id, (const char *) header + 8, 8 );
-        ct->id[8] = '\0';
-
-        ct->ll.lam = - to_double(header+7*16+8); /* W_LONG */
-        ct->ll.phi = to_double(header+4*16+8);   /* S_LAT */
-
-        ur.lam = - to_double(header+6*16+8);     /* E_LONG */
-        ur.phi = to_double(header+5*16+8);       /* N_LAT */
-
-        ct->del.lam = to_double(header+9*16+8);
-        ct->del.phi = to_double(header+8*16+8);
-
-        ct->lim.lam = (pj_int32) (fabs(ur.lam-ct->ll.lam)/ct->del.lam + 0.5) + 1;
-        ct->lim.phi = (pj_int32) (fabs(ur.phi-ct->ll.phi)/ct->del.phi + 0.5) + 1;
-
-        pj_log( ctx, PJ_LOG_DEBUG_MINOR,
-                "NTv2 %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)",
-                ct->id,
-                ct->lim.lam, ct->lim.phi,
-                ct->ll.lam/3600.0, ct->ll.phi/3600.0,
-                ur.lam/3600.0, ur.phi/3600.0 );
-
-        ct->ll.lam *= DEG_TO_RAD/3600.0;
-        ct->ll.phi *= DEG_TO_RAD/3600.0;
-        ct->del.lam *= DEG_TO_RAD/3600.0;
-        ct->del.phi *= DEG_TO_RAD/3600.0;
-
-        memcpy( &gs_count, header + 8 + 16*10, 4 );
-        if( gs_count != ct->lim.lam * ct->lim.phi )
-        {
-            pj_log( ctx, PJ_LOG_ERROR,
-                    "GS_COUNT(%d) does not match expected cells (%dx%d=%d)",
-                    gs_count, ct->lim.lam, ct->lim.phi,
-                    ct->lim.lam * ct->lim.phi );
-            pj_dalloc(ct);
-            pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-            return 0;
-        }
-
-        ct->cvs = nullptr;
-
-/* -------------------------------------------------------------------- */
-/*      Create a new gridinfo for this if we aren't processing the      */
-/*      1st subfile, and initialize our grid info.                      */
-/* -------------------------------------------------------------------- */
-        if( subfile == 0 )
-            gi = gilist;
-        else
-        {
-            gi = (PJ_GRIDINFO *) pj_calloc(1, sizeof(PJ_GRIDINFO));
-            if (!gi) {
-                pj_dalloc(ct);
-                pj_gridinfo_free(ctx, gilist);
-                pj_ctx_set_errno(ctx, ENOMEM);
-                return 0;
-            }
-
-            gi->gridname = pj_strdup( gilist->gridname );
-            gi->filename = pj_strdup( gilist->filename );
-            if (!gi->gridname || !gi->filename) {
-                pj_gridinfo_free(ctx, gi);
-                pj_dalloc(ct);
-                pj_gridinfo_free(ctx, gilist);
-                pj_ctx_set_errno(ctx, ENOMEM);
-                return 0;
-            }
-            gi->next = nullptr;
-        }
-
-        gi->must_swap = must_swap;
-        gi->ct = ct;
-        gi->format = "ntv2";
-        gi->grid_offset = pj_ctx_ftell( ctx, fid );
-
-/* -------------------------------------------------------------------- */
-/*      Attach to the correct list or sublist.                          */
-/* -------------------------------------------------------------------- */
-        if( strncmp((const char *)header+24,"NONE",4) == 0 )
-        {
-            if( gi != gilist )
-            {
-                PJ_GRIDINFO *lnk;
-
-                for( lnk = gilist; lnk->next != nullptr; lnk = lnk->next ) {}
-                lnk->next = gi;
-            }
-        }
-
-        else
-        {
-            PJ_GRIDINFO *lnk;
-            PJ_GRIDINFO *gp = gridinfo_parent(gilist,
-                                                 (const char*)header+24,8);
-
-            if( gp == nullptr )
-            {
-                pj_log( ctx, PJ_LOG_ERROR,
-                        "pj_gridinfo_init_ntv2(): "
-                        "failed to find parent %8.8s for %s.",
-                        (const char *) header+24, gi->ct->id );
-
-                for( lnk = gilist; lnk->next != nullptr; lnk = lnk->next ) {}
-                lnk->next = gi;
-            }
-            else
-            {
-                if( gp->child == nullptr )
-                {
-                    gp->child = gi;
-                }
-                else
-                {
-                    for( lnk = gp->child; lnk->next != nullptr; lnk = lnk->next ) {}
-                    lnk->next = gi;
-                }
-            }
-        }
-
-/* -------------------------------------------------------------------- */
-/*      Seek past the data.                                             */
-/* -------------------------------------------------------------------- */
-        pj_ctx_fseek( ctx, fid, gs_count * 16, SEEK_CUR );
-    }
-
-    return 1;
-}
-
-/************************************************************************/
-/*                       pj_gridinfo_init_ntv1()                        */
-/*                                                                      */
-/*      Load an NTv1 style Canadian grid shift file.                    */
-/************************************************************************/
-
-static int pj_gridinfo_init_ntv1( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi )
-
-{
-    unsigned char header[192]; /* 12 records of 16 bytes */
-    struct CTABLE *ct;
-    PJ_LP		ur;
-
-    /* cppcheck-suppress sizeofCalculation */
-    STATIC_ASSERT( sizeof(pj_int32) == 4 );
-    /* cppcheck-suppress sizeofCalculation */
-    STATIC_ASSERT( sizeof(double) == 8 );
-
-/* -------------------------------------------------------------------- */
-/*      Read the header.                                                */
-/* -------------------------------------------------------------------- */
-    if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 )
-    {
-        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        return 0;
-    }
-
-/* -------------------------------------------------------------------- */
-/*      Regularize fields of interest.                                  */
-/* -------------------------------------------------------------------- */
-    if( IS_LSB )
-    {
-        swap_words( header+8, 4, 1 );
-        swap_words( header+24, 8, 1 );
-        swap_words( header+40, 8, 1 );
-        swap_words( header+56, 8, 1 );
-        swap_words( header+72, 8, 1 );
-        swap_words( header+88, 8, 1 );
-        swap_words( header+104, 8, 1 );
-    }
-
-    if( *((int *) (header+8)) != 12 )
-    {
-        pj_log( ctx, PJ_LOG_ERROR,
-                "NTv1 grid shift file has wrong record count, corrupt?" );
-        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        return 0;
-    }
-
-/* -------------------------------------------------------------------- */
-/*      Fill in CTABLE structure.                                       */
-/* -------------------------------------------------------------------- */
-    ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE));
-    if (!ct) {
-        pj_ctx_set_errno(ctx, ENOMEM);
-        return 0;
-    }
-    strcpy( ct->id, "NTv1 Grid Shift File" );
-
-    ct->ll.lam = - to_double(header+72);
-    ct->ll.phi = to_double(header+24);
-    ur.lam = - to_double(header+56);
-    ur.phi = to_double(header+40);
-    ct->del.lam = to_double(header+104);
-    ct->del.phi = to_double(header+88);
-    ct->lim.lam = (pj_int32) (fabs(ur.lam-ct->ll.lam)/ct->del.lam + 0.5) + 1;
-    ct->lim.phi = (pj_int32) (fabs(ur.phi-ct->ll.phi)/ct->del.phi + 0.5) + 1;
-
-    pj_log( ctx, PJ_LOG_DEBUG_MINOR,
-            "NTv1 %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)",
-            ct->lim.lam, ct->lim.phi,
-            ct->ll.lam, ct->ll.phi, ur.lam, ur.phi );
-
-    ct->ll.lam *= DEG_TO_RAD;
-    ct->ll.phi *= DEG_TO_RAD;
-    ct->del.lam *= DEG_TO_RAD;
-    ct->del.phi *= DEG_TO_RAD;
-    ct->cvs = nullptr;
-
-    gi->ct = ct;
-    gi->grid_offset = (long) sizeof(header);
-    gi->format = "ntv1";
-
-    return 1;
-}
-
-/************************************************************************/
-/*                       pj_gridinfo_init_gtx()                         */
-/*                                                                      */
-/*      Load a NOAA .gtx vertical datum shift file.                     */
-/************************************************************************/
-
-static int pj_gridinfo_init_gtx( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi )
-
-{
-    unsigned char header[40];
-    struct CTABLE *ct;
-    double      xorigin,yorigin,xstep,ystep;
-    int         rows, columns;
-
-    /* cppcheck-suppress sizeofCalculation */
-    STATIC_ASSERT( sizeof(pj_int32) == 4 );
-    /* cppcheck-suppress sizeofCalculation */
-    STATIC_ASSERT( sizeof(double) == 8 );
-
-/* -------------------------------------------------------------------- */
-/*      Read the header.                                                */
-/* -------------------------------------------------------------------- */
-    if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 )
-    {
-        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        return 0;
-    }
-
-/* -------------------------------------------------------------------- */
-/*      Regularize fields of interest and extract.                      */
-/* -------------------------------------------------------------------- */
-    if( IS_LSB )
-    {
-        swap_words( header+0, 8, 4 );
-        swap_words( header+32, 4, 2 );
-    }
-
-    memcpy( &yorigin, header+0, 8 );
-    memcpy( &xorigin, header+8, 8 );
-    memcpy( &ystep, header+16, 8 );
-    memcpy( &xstep, header+24, 8 );
-
-    memcpy( &rows, header+32, 4 );
-    memcpy( &columns, header+36, 4 );
-
-    if( xorigin < -360 || xorigin > 360
-        || yorigin < -90 || yorigin > 90 )
-    {
-        pj_log( ctx, PJ_LOG_ERROR,
-                "gtx file header has invalid extents, corrupt?");
-        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        return 0;
-    }
-
-/* -------------------------------------------------------------------- */
-/*      Fill in CTABLE structure.                                       */
-/* -------------------------------------------------------------------- */
-    ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE));
-    if (!ct) {
-        pj_ctx_set_errno(ctx, ENOMEM);
-        return 0;
-    }
-    strcpy( ct->id, "GTX Vertical Grid Shift File" );
-
-    ct->ll.lam = xorigin;
-    ct->ll.phi = yorigin;
-    ct->del.lam = xstep;
-    ct->del.phi = ystep;
-    ct->lim.lam = columns;
-    ct->lim.phi = rows;
-
-    /* some GTX files come in 0-360 and we shift them back into the
-       expected -180 to 180 range if possible.  This does not solve
-       problems with grids spanning the dateline. */
-    if( ct->ll.lam >= 180.0 )
-        ct->ll.lam -= 360.0;
-
-    if( ct->ll.lam >= 0.0 && ct->ll.lam + ct->del.lam * ct->lim.lam > 180.0 )
-    {
-        pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
-                "This GTX spans the dateline!  This will cause problems." );
-    }
-
-    pj_log( ctx, PJ_LOG_DEBUG_MINOR,
-            "GTX %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)",
-            ct->lim.lam, ct->lim.phi,
-            ct->ll.lam, ct->ll.phi,
-            ct->ll.lam + (columns-1)*xstep, ct->ll.phi + (rows-1)*ystep);
-
-    ct->ll.lam *= DEG_TO_RAD;
-    ct->ll.phi *= DEG_TO_RAD;
-    ct->del.lam *= DEG_TO_RAD;
-    ct->del.phi *= DEG_TO_RAD;
-    ct->cvs = nullptr;
-
-    gi->ct = ct;
-    gi->grid_offset = 40;
-    gi->format = "gtx";
-
-    return 1;
-}
-
-/************************************************************************/
-/*                          pj_gridinfo_init()                          */
-/*                                                                      */
-/*      Open and parse header details from a datum gridshift file       */
-/*      returning a list of PJ_GRIDINFOs for the grids in that          */
-/*      file.  This superceeds use of nad_init() for modern             */
-/*      applications.                                                   */
-/************************************************************************/
-
-PJ_GRIDINFO *pj_gridinfo_init( projCtx ctx, const char *gridname )
-
-{
-    PJ_GRIDINFO *gilist;
-    PAFile 	fp;
-    char	header[160];
-    size_t      header_size = 0;
-
-    errno = pj_errno = 0;
-    ctx->last_errno = 0;
-
-/* -------------------------------------------------------------------- */
-/*      Initialize a GRIDINFO with stub info we would use if it         */
-/*      cannot be loaded.                                               */
-/* -------------------------------------------------------------------- */
-    gilist = (PJ_GRIDINFO *) pj_calloc(1, sizeof(PJ_GRIDINFO));
-    if (!gilist) {
-        pj_ctx_set_errno(ctx, ENOMEM);
-        return nullptr;
-    }
-
-    gilist->gridname = pj_strdup( gridname );
-    if (!gilist->gridname) {
-        pj_dalloc(gilist);
-        pj_ctx_set_errno(ctx, ENOMEM);
-        return nullptr;
-    }
-    gilist->filename = nullptr;
-    gilist->format = "missing";
-    gilist->grid_offset = 0;
-    gilist->ct = nullptr;
-    gilist->next = nullptr;
-
-/* -------------------------------------------------------------------- */
-/*      Open the file using the usual search rules.                     */
-/* -------------------------------------------------------------------- */
-    if (!(fp = pj_open_lib(ctx, gridname, "rb"))) {
-        ctx->last_errno = 0; /* don't treat as a persistent error */
-        return gilist;
-    }
-
-    gilist->filename = pj_strdup(gridname);
-    if (!gilist->filename) {
-        pj_dalloc(gilist->gridname);
-        pj_dalloc(gilist);
-        pj_ctx_set_errno(ctx, ENOMEM);
-        return nullptr;
-    }
-
-/* -------------------------------------------------------------------- */
-/*      Load a header, to determine the file type.                      */
-/* -------------------------------------------------------------------- */
-    if( (header_size = pj_ctx_fread( ctx, header, 1,
-                                     sizeof(header), fp ) ) != sizeof(header) )
-    {
-        /* some files may be smaller that sizeof(header), eg 160, so */
-        ctx->last_errno = 0; /* don't treat as a persistent error */
-        pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
-                "pj_gridinfo_init: short header read of %d bytes",
-                (int)header_size );
-    }
-
-    pj_ctx_fseek( ctx, fp, SEEK_SET, 0 );
-
-/* -------------------------------------------------------------------- */
-/*      Determine file type.                                            */
-/* -------------------------------------------------------------------- */
-    if( header_size >= 144 + 16
-        && strncmp(header + 0, "HEADER", 6) == 0
-        && strncmp(header + 96, "W GRID", 6) == 0
-        && strncmp(header + 144, "TO      NAD83   ", 16) == 0 )
-    {
-        pj_gridinfo_init_ntv1( ctx, fp, gilist );
-    }
-
-    else if( header_size >= 48 + 7
-             && strncmp(header + 0, "NUM_OREC", 8) == 0
-             && strncmp(header + 48, "GS_TYPE", 7) == 0 )
-    {
-        pj_gridinfo_init_ntv2( ctx, fp, gilist );
-    }
-
-    else if( strlen(gridname) > 4
-             && (strcmp(gridname+strlen(gridname)-3,"gtx") == 0
-                 || strcmp(gridname+strlen(gridname)-3,"GTX") == 0) )
-    {
-        pj_gridinfo_init_gtx( ctx, fp, gilist );
-    }
-
-    else if( header_size >= 9 && strncmp(header + 0,"CTABLE V2",9) == 0 )
-    {
-        struct CTABLE *ct = nad_ctable2_init( ctx, (struct projFileAPI_t*)fp );
-
-        gilist->format = "ctable2";
-        gilist->ct = ct;
-
-        if (ct == nullptr)
-        {
-            pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
-                    "CTABLE V2 ct is NULL.");
-        }
-        else
-        {
-            pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
-                    "Ctable2 %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)",
-                    ct->id,
-                    ct->lim.lam, ct->lim.phi,
-                    ct->ll.lam * RAD_TO_DEG, ct->ll.phi * RAD_TO_DEG,
-                    (ct->ll.lam + (ct->lim.lam-1)*ct->del.lam) * RAD_TO_DEG,
-                    (ct->ll.phi + (ct->lim.phi-1)*ct->del.phi) * RAD_TO_DEG );
-        }
-    }
-
-    else
-    {
-        struct CTABLE *ct = nad_ctable_init( ctx, (struct projFileAPI_t*)fp );
-        if (ct == nullptr)
-        {
-            pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
-                    "CTABLE ct is NULL.");
-        } else
-        {
-            gilist->format = "ctable";
-            gilist->ct = ct;
-
-            pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
-                    "Ctable %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)",
-                    ct->id,
-                    ct->lim.lam, ct->lim.phi,
-                    ct->ll.lam * RAD_TO_DEG, ct->ll.phi * RAD_TO_DEG,
-                    (ct->ll.lam + (ct->lim.lam-1)*ct->del.lam) * RAD_TO_DEG,
-                    (ct->ll.phi + (ct->lim.phi-1)*ct->del.phi) * RAD_TO_DEG );
-        }
-    }
-
-    pj_ctx_fclose(ctx, fp);
-
-    return gilist;
-}
diff --git a/src/gridlist.cpp b/src/gridlist.cpp
deleted file mode 100644
index c540b134..00000000
--- a/src/gridlist.cpp
+++ /dev/null
@@ -1,220 +0,0 @@
-/******************************************************************************
- * Project:  PROJ.4
- * Purpose:  Code to manage the list of currently loaded (cached) PJ_GRIDINFOs
- *           See pj_gridinfo.c for details of loading individual grids.
- * Author:   Frank Warmerdam, warmerdam@pobox.com
- *
- ******************************************************************************
- * Copyright (c) 2000, Frank Warmerdam <warmerdam@pobox.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 <stddef.h>
-#include <string.h>
-
-#include "proj.h"
-#include "proj_internal.h"
-
-static PJ_GRIDINFO *grid_list = nullptr;
-#define PJ_MAX_PATH_LENGTH 1024
-
-/************************************************************************/
-/*                        pj_deallocate_grids()                         */
-/*                                                                      */
-/*      Deallocate all loaded grids.                                    */
-/************************************************************************/
-
-void pj_deallocate_grids()
-
-{
-    while( grid_list != nullptr )
-    {
-        PJ_GRIDINFO *item = grid_list;
-        grid_list = grid_list->next;
-        item->next = nullptr;
-
-        pj_gridinfo_free( pj_get_default_ctx(), item );
-    }
-}
-
-/************************************************************************/
-/*                       pj_gridlist_merge_grid()                       */
-/*                                                                      */
-/*      Find/load the named gridfile and merge it into the              */
-/*      last_nadgrids_list.                                             */
-/************************************************************************/
-
-static int pj_gridlist_merge_gridfile( projCtx ctx, 
-                                       const char *gridname,
-                                       PJ_GRIDINFO ***p_gridlist,
-                                       int *p_gridcount, 
-                                       int *p_gridmax )
-
-{
-    int got_match=0;
-    PJ_GRIDINFO *this_grid, *tail = nullptr;
-
-/* -------------------------------------------------------------------- */
-/*      Try to find in the existing list of loaded grids.  Add all      */
-/*      matching grids as with NTv2 we can get many grids from one      */
-/*      file (one shared gridname).                                     */
-/* -------------------------------------------------------------------- */
-    for( this_grid = grid_list; this_grid != nullptr; this_grid = this_grid->next)
-    {
-        if( strcmp(this_grid->gridname,gridname) == 0 )
-        {
-            got_match = 1;
-
-            /* don't add to the list if it is invalid. */
-            if( this_grid->ct == nullptr )
-                return 0;
-
-            /* do we need to grow the list? */
-            if( *p_gridcount >= *p_gridmax - 2 )
-            {
-                PJ_GRIDINFO **new_list;
-                int new_max = *p_gridmax + 20;
-
-                new_list = (PJ_GRIDINFO **) pj_calloc(new_max, sizeof(void *));
-                if (!new_list) {
-                    pj_ctx_set_errno( ctx, ENOMEM );
-                    return 0;
-                }
-                if( *p_gridlist != nullptr )
-                {
-                    memcpy( new_list, *p_gridlist,
-                            sizeof(void *) * (*p_gridmax) );
-                    pj_dalloc( *p_gridlist );
-                }
-
-                *p_gridlist = new_list;
-                *p_gridmax = new_max;
-            }
-
-            /* add to the list */
-            (*p_gridlist)[(*p_gridcount)++] = this_grid;
-            (*p_gridlist)[*p_gridcount] = nullptr;
-        }
-
-        tail = this_grid;
-    }
-
-    if( got_match )
-        return 1;
-
-/* -------------------------------------------------------------------- */
-/*      Try to load the named grid.                                     */
-/* -------------------------------------------------------------------- */
-    this_grid = pj_gridinfo_init( ctx, gridname );
-
-    if( this_grid == nullptr )
-    {
-        return 0;
-    }
-    
-    if( tail != nullptr )
-        tail->next = this_grid;
-    else
-        grid_list = this_grid;
-
-/* -------------------------------------------------------------------- */
-/*      Recurse to add the grid now that it is loaded.                  */
-/* -------------------------------------------------------------------- */
-    return pj_gridlist_merge_gridfile( ctx, gridname, p_gridlist, 
-                                       p_gridcount, p_gridmax );
-}
-
-/************************************************************************/
-/*                     pj_gridlist_from_nadgrids()                      */
-/*                                                                      */
-/*      This functions loads the list of grids corresponding to a       */
-/*      particular nadgrids string into a list, and returns it.  The    */
-/*      list is kept around till a request is made with a different     */
-/*      string in order to cut down on the string parsing cost, and     */
-/*      the cost of building the list of tables each time.              */
-/************************************************************************/
-
-PJ_GRIDINFO **pj_gridlist_from_nadgrids( projCtx ctx, const char *nadgrids, 
-                                         int *grid_count)
-
-{
-    const char *s;
-    PJ_GRIDINFO **gridlist = nullptr;
-    int grid_max = 0;
-
-    pj_errno = 0;
-    *grid_count = 0;
-
-    pj_acquire_lock();
-
-/* -------------------------------------------------------------------- */
-/*      Loop processing names out of nadgrids one at a time.            */
-/* -------------------------------------------------------------------- */
-    for( s = nadgrids; *s != '\0'; )
-    {
-        size_t end_char;
-        int    required = 1;
-        char   name[PJ_MAX_PATH_LENGTH];
-
-        if( *s == '@' )
-        {
-            required = 0;
-            s++;
-        }
-
-        for( end_char = 0; 
-             s[end_char] != '\0' && s[end_char] != ','; 
-             end_char++ ) {}
-
-        if( end_char >= sizeof(name) )
-        {
-            pj_dalloc( gridlist );
-            pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-            pj_release_lock();
-            return nullptr;
-        }
-        
-        strncpy( name, s, end_char );
-        name[end_char] = '\0';
-
-        s += end_char;
-        if( *s == ',' )
-            s++;
-
-        if( !pj_gridlist_merge_gridfile( ctx, name, &gridlist, grid_count, 
-                                         &grid_max) 
-            && required )
-        {
-            pj_dalloc( gridlist );
-            pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-            pj_release_lock();
-            return nullptr;
-        }
-        else
-            pj_errno = 0;
-    }
-
-    pj_release_lock();
-
-    return gridlist;
-}
diff --git a/src/grids.cpp b/src/grids.cpp
new file mode 100644
index 00000000..7d19b1f7
--- /dev/null
+++ b/src/grids.cpp
@@ -0,0 +1,895 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Grid management
+ * Author:   Even Rouault, <even.rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2000, Frank Warmerdam <warmerdam@pobox.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.
+ *****************************************************************************/
+
+#ifndef FROM_PROJ_CPP
+#define FROM_PROJ_CPP
+#endif
+
+#include "grids.hpp"
+#include "proj/internal/internal.hpp"
+#include "proj_internal.h"
+
+NS_PROJ_START
+
+using namespace internal;
+
+/************************************************************************/
+/*                             swap_words()                             */
+/*                                                                      */
+/*      Convert the byte order of the given word(s) in place.           */
+/************************************************************************/
+
+static const int byte_order_test = 1;
+#define IS_LSB (1 == ((const unsigned char *)(&byte_order_test))[0])
+
+static void swap_words(void *dataIn, size_t word_size, size_t word_count)
+
+{
+    unsigned char *data = static_cast<unsigned char *>(dataIn);
+    for (size_t word = 0; word < word_count; word++) {
+        for (size_t i = 0; i < word_size / 2; i++) {
+            unsigned char t;
+
+            t = data[i];
+            data[i] = data[word_size - i - 1];
+            data[word_size - i - 1] = t;
+        }
+
+        data += word_size;
+    }
+}
+
+bool ExtentAndRes::fullWorldLongitude() const {
+    return eastLon - westLon + resLon >= 2 * M_PI - 1e-10;
+}
+
+// ---------------------------------------------------------------------------
+
+Grid::Grid(int widthIn, int heightIn, const ExtentAndRes &extentIn)
+    : m_width(widthIn), m_height(heightIn), m_extent(extentIn) {}
+
+// ---------------------------------------------------------------------------
+
+Grid::~Grid() = default;
+
+// ---------------------------------------------------------------------------
+
+VerticalShiftGrid::VerticalShiftGrid(int widthIn, int heightIn,
+                                     const ExtentAndRes &extentIn)
+    : Grid(widthIn, heightIn, extentIn) {}
+
+// ---------------------------------------------------------------------------
+
+bool VerticalShiftGrid::isNodata(float /*val*/, double /* multiplier */) const {
+    return false;
+}
+
+// ---------------------------------------------------------------------------
+
+static ExtentAndRes globalExtent() {
+    ExtentAndRes extent;
+    extent.westLon = -M_PI;
+    extent.southLat = -M_PI / 2;
+    extent.eastLon = M_PI;
+    extent.northLat = M_PI / 2;
+    extent.resLon = M_PI;
+    extent.resLat = M_PI / 2;
+    return extent;
+}
+
+// ---------------------------------------------------------------------------
+
+class NullVerticalShiftGrid : public VerticalShiftGrid {
+
+  public:
+    NullVerticalShiftGrid() : VerticalShiftGrid(3, 3, globalExtent()) {}
+
+    bool isNullGrid() const override { return true; }
+    bool valueAt(int, int, float &out) const override;
+    bool isNodata(float, double) const override { return false; }
+};
+
+// ---------------------------------------------------------------------------
+
+bool NullVerticalShiftGrid::valueAt(int, int, float &out) const {
+    out = 0.0f;
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+class GTXVerticalShiftGrid : public VerticalShiftGrid {
+    PJ_CONTEXT *m_ctx;
+    PAFile m_fp;
+
+    GTXVerticalShiftGrid(const GTXVerticalShiftGrid &) = delete;
+    GTXVerticalShiftGrid &operator=(const GTXVerticalShiftGrid &) = delete;
+
+  public:
+    GTXVerticalShiftGrid(PJ_CONTEXT *ctx, PAFile fp, int widthIn, int heightIn,
+                         const ExtentAndRes &extentIn)
+        : VerticalShiftGrid(widthIn, heightIn, extentIn), m_ctx(ctx), m_fp(fp) {
+    }
+
+    ~GTXVerticalShiftGrid() override;
+
+    bool valueAt(int x, int y, float &out) const override;
+    bool isNodata(float val, double multiplier) const override;
+
+    static GTXVerticalShiftGrid *open(PJ_CONTEXT *ctx, PAFile fp);
+};
+
+// ---------------------------------------------------------------------------
+
+GTXVerticalShiftGrid::~GTXVerticalShiftGrid() { pj_ctx_fclose(m_ctx, m_fp); }
+
+// ---------------------------------------------------------------------------
+
+GTXVerticalShiftGrid *GTXVerticalShiftGrid::open(PJ_CONTEXT *ctx, PAFile fp) {
+    unsigned char header[40];
+
+    /* -------------------------------------------------------------------- */
+    /*      Read the header.                                                */
+    /* -------------------------------------------------------------------- */
+    if (pj_ctx_fread(ctx, header, sizeof(header), 1, fp) != 1) {
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return nullptr;
+    }
+
+    /* -------------------------------------------------------------------- */
+    /*      Regularize fields of interest and extract.                      */
+    /* -------------------------------------------------------------------- */
+    if (IS_LSB) {
+        swap_words(header + 0, 8, 4);
+        swap_words(header + 32, 4, 2);
+    }
+
+    double xorigin, yorigin, xstep, ystep;
+    int rows, columns;
+
+    memcpy(&yorigin, header + 0, 8);
+    memcpy(&xorigin, header + 8, 8);
+    memcpy(&ystep, header + 16, 8);
+    memcpy(&xstep, header + 24, 8);
+
+    memcpy(&rows, header + 32, 4);
+    memcpy(&columns, header + 36, 4);
+
+    if (xorigin < -360 || xorigin > 360 || yorigin < -90 || yorigin > 90) {
+        pj_log(ctx, PJ_LOG_ERROR,
+               "gtx file header has invalid extents, corrupt?");
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return nullptr;
+    }
+
+    /* some GTX files come in 0-360 and we shift them back into the
+       expected -180 to 180 range if possible.  This does not solve
+       problems with grids spanning the dateline. */
+    if (xorigin >= 180.0)
+        xorigin -= 360.0;
+
+    if (xorigin >= 0.0 && xorigin + xstep * columns > 180.0) {
+        pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+               "This GTX spans the dateline!  This will cause problems.");
+    }
+
+    ExtentAndRes extent;
+    extent.westLon = xorigin * DEG_TO_RAD;
+    extent.southLat = yorigin * DEG_TO_RAD;
+    extent.resLon = xstep * DEG_TO_RAD;
+    extent.resLat = ystep * DEG_TO_RAD;
+    extent.eastLon = (xorigin + xstep * (columns - 1)) * DEG_TO_RAD;
+    extent.northLat = (yorigin + ystep * (rows - 1)) * DEG_TO_RAD;
+
+    return new GTXVerticalShiftGrid(ctx, fp, columns, rows, extent);
+}
+
+// ---------------------------------------------------------------------------
+
+bool GTXVerticalShiftGrid::valueAt(int x, int y, float &out) const {
+    assert(x >= 0 && y >= 0 && x < m_width && y < m_height);
+
+    pj_ctx_fseek(m_ctx, m_fp, 40 + sizeof(float) * (y * m_width + x), SEEK_SET);
+    if (pj_ctx_fread(m_ctx, &out, sizeof(out), 1, m_fp) != 1) {
+        pj_ctx_set_errno(m_ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return false;
+    }
+    if (IS_LSB) {
+        swap_words(&out, sizeof(float), 1);
+    }
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+bool GTXVerticalShiftGrid::isNodata(float val, double multiplier) const {
+    /* nodata?  */
+    /* GTX official nodata value if  -88.88880f, but some grids also */
+    /* use other  big values for nodata (e.g naptrans2008.gtx has */
+    /* nodata values like -2147479936), so test them too */
+    return val * multiplier > 1000 || val * multiplier < -1000 ||
+           val == -88.88880f;
+}
+
+// ---------------------------------------------------------------------------
+
+VerticalShiftGridSet::VerticalShiftGridSet() = default;
+
+// ---------------------------------------------------------------------------
+
+VerticalShiftGridSet::~VerticalShiftGridSet() = default;
+
+// ---------------------------------------------------------------------------
+
+std::unique_ptr<VerticalShiftGridSet>
+VerticalShiftGridSet::open(PJ_CONTEXT *ctx, const std::string &filename) {
+    if (filename == "null") {
+        auto set =
+            std::unique_ptr<VerticalShiftGridSet>(new VerticalShiftGridSet());
+        set->m_name = filename;
+        set->m_format = "null";
+        set->m_grids.push_back(std::unique_ptr<NullVerticalShiftGrid>(
+            new NullVerticalShiftGrid()));
+        return set;
+    }
+
+    PAFile fp;
+    if (!(fp = pj_open_lib(ctx, filename.c_str(), "rb"))) {
+        ctx->last_errno = 0; /* don't treat as a persistent error */
+        return nullptr;
+    }
+    if (ends_with(filename, "gtx") || ends_with(filename, "GTX")) {
+        auto grid = GTXVerticalShiftGrid::open(ctx, fp);
+        if (!grid) {
+            pj_ctx_fclose(ctx, fp);
+            return nullptr;
+        }
+        auto set =
+            std::unique_ptr<VerticalShiftGridSet>(new VerticalShiftGridSet());
+        set->m_name = filename;
+        set->m_format = "gtx";
+        set->m_grids.push_back(std::unique_ptr<VerticalShiftGrid>(grid));
+        return set;
+    }
+
+    pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Unrecognized vertical grid format");
+    pj_ctx_fclose(ctx, fp);
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+const VerticalShiftGrid *VerticalShiftGrid::gridAt(double lon,
+                                                   double lat) const {
+    for (const auto &child : m_children) {
+        const auto &extentChild = child->extentAndRes();
+        if ((extentChild.fullWorldLongitude() ||
+             (lon >= extentChild.westLon && lon <= extentChild.eastLon)) &&
+            lat >= extentChild.southLat && lat <= extentChild.northLat) {
+            return child->gridAt(lon, lat);
+        }
+    }
+    return this;
+}
+// ---------------------------------------------------------------------------
+
+const VerticalShiftGrid *VerticalShiftGridSet::gridAt(double lon,
+                                                      double lat) const {
+    for (const auto &grid : m_grids) {
+        if (dynamic_cast<NullVerticalShiftGrid *>(grid.get())) {
+            return grid.get();
+        }
+        const auto &extent = grid->extentAndRes();
+        if ((extent.fullWorldLongitude() ||
+             (lon >= extent.westLon && lon <= extent.eastLon)) &&
+            lat >= extent.southLat && lat <= extent.northLat) {
+            return grid->gridAt(lon, lat);
+        }
+    }
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+HorizontalShiftGrid::HorizontalShiftGrid(int widthIn, int heightIn,
+                                         const ExtentAndRes &extentIn)
+    : Grid(widthIn, heightIn, extentIn) {}
+
+// ---------------------------------------------------------------------------
+
+HorizontalShiftGrid::~HorizontalShiftGrid() = default;
+
+// ---------------------------------------------------------------------------
+
+HorizontalShiftGridSet::HorizontalShiftGridSet() = default;
+
+// ---------------------------------------------------------------------------
+
+HorizontalShiftGridSet::~HorizontalShiftGridSet() = default;
+
+// ---------------------------------------------------------------------------
+
+class NullHorizontalShiftGrid : public HorizontalShiftGrid {
+
+  public:
+    NullHorizontalShiftGrid() : HorizontalShiftGrid(3, 3, globalExtent()) {}
+
+    bool isNullGrid() const override { return true; }
+
+    bool valueAt(int, int, float &lonShift, float &latShift) const override;
+};
+
+// ---------------------------------------------------------------------------
+
+bool NullHorizontalShiftGrid::valueAt(int, int, float &lonShift,
+                                      float &latShift) const {
+    lonShift = 0.0f;
+    latShift = 0.0f;
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+static double to_double(const void *data) {
+    double d;
+    memcpy(&d, data, sizeof(d));
+    return d;
+}
+
+// ---------------------------------------------------------------------------
+
+class NTv1Grid : public HorizontalShiftGrid {
+    PJ_CONTEXT *m_ctx;
+    PAFile m_fp;
+
+    NTv1Grid(const NTv1Grid &) = delete;
+    NTv1Grid &operator=(const NTv1Grid &) = delete;
+
+  public:
+    NTv1Grid(PJ_CONTEXT *ctx, PAFile fp, int widthIn, int heightIn,
+             const ExtentAndRes &extentIn)
+        : HorizontalShiftGrid(widthIn, heightIn, extentIn), m_ctx(ctx),
+          m_fp(fp) {}
+
+    ~NTv1Grid() override;
+
+    bool valueAt(int, int, float &lonShift, float &latShift) const override;
+
+    static NTv1Grid *open(PJ_CONTEXT *ctx, PAFile fp,
+                          const std::string &filename);
+};
+
+// ---------------------------------------------------------------------------
+
+NTv1Grid::~NTv1Grid() { pj_ctx_fclose(m_ctx, m_fp); }
+
+// ---------------------------------------------------------------------------
+
+NTv1Grid *NTv1Grid::open(PJ_CONTEXT *ctx, PAFile fp,
+                         const std::string &filename) {
+    unsigned char header[192];
+
+    /* -------------------------------------------------------------------- */
+    /*      Read the header.                                                */
+    /* -------------------------------------------------------------------- */
+    if (pj_ctx_fread(ctx, header, sizeof(header), 1, fp) != 1) {
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return nullptr;
+    }
+
+    /* -------------------------------------------------------------------- */
+    /*      Regularize fields of interest.                                  */
+    /* -------------------------------------------------------------------- */
+    if (IS_LSB) {
+        swap_words(header + 8, sizeof(int), 1);
+        swap_words(header + 24, sizeof(double), 1);
+        swap_words(header + 40, sizeof(double), 1);
+        swap_words(header + 56, sizeof(double), 1);
+        swap_words(header + 72, sizeof(double), 1);
+        swap_words(header + 88, sizeof(double), 1);
+        swap_words(header + 104, sizeof(double), 1);
+    }
+
+    if (*((int *)(header + 8)) != 12) {
+        pj_log(ctx, PJ_LOG_ERROR,
+               "NTv1 grid shift file has wrong record count, corrupt?");
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return nullptr;
+    }
+
+    ExtentAndRes extent;
+    extent.westLon = -to_double(header + 72) * DEG_TO_RAD;
+    extent.southLat = to_double(header + 24) * DEG_TO_RAD;
+    extent.eastLon = -to_double(header + 56) * DEG_TO_RAD;
+    extent.northLat = to_double(header + 40) * DEG_TO_RAD;
+    extent.resLon = to_double(header + 104) * DEG_TO_RAD;
+    extent.resLat = to_double(header + 88) * DEG_TO_RAD;
+    if (!(fabs(extent.westLon) <= 4 * M_PI &&
+          fabs(extent.eastLon) <= 4 * M_PI &&
+          fabs(extent.northLat) <= M_PI + 1e-5 &&
+          fabs(extent.southLat) <= M_PI + 1e-5 &&
+          extent.westLon < extent.eastLon &&
+          extent.southLat < extent.northLat && extent.resLon > 1e-10 &&
+          extent.resLat > 1e-10)) {
+        pj_log(ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s",
+               filename.c_str());
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return nullptr;
+    }
+    const int columns = static_cast<int>(
+        fabs((extent.eastLon - extent.westLon) / extent.resLon + 0.5) + 1);
+    const int rows = static_cast<int>(
+        fabs((extent.northLat - extent.southLat) / extent.resLat + 0.5) + 1);
+
+    return new NTv1Grid(ctx, fp, columns, rows, extent);
+}
+
+// ---------------------------------------------------------------------------
+
+bool NTv1Grid::valueAt(int x, int y, float &lonShift, float &latShift) const {
+    assert(x >= 0 && y >= 0 && x < m_width && y < m_height);
+
+    double two_doubles[2];
+    // NTv1 is organized from east to west !
+    pj_ctx_fseek(m_ctx, m_fp,
+                 192 + 2 * sizeof(double) * (y * m_width + m_width - 1 - x),
+                 SEEK_SET);
+    if (pj_ctx_fread(m_ctx, &two_doubles[0], sizeof(two_doubles), 1, m_fp) !=
+        1) {
+        pj_ctx_set_errno(m_ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return false;
+    }
+    if (IS_LSB) {
+        swap_words(&two_doubles[0], sizeof(double), 2);
+    }
+    /* convert seconds to radians */
+    latShift = static_cast<float>(two_doubles[0] * ((M_PI / 180.0) / 3600.0));
+    lonShift = static_cast<float>(two_doubles[1] * ((M_PI / 180.0) / 3600.0));
+
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+class CTable2Grid : public HorizontalShiftGrid {
+    PJ_CONTEXT *m_ctx;
+    PAFile m_fp;
+
+    CTable2Grid(const CTable2Grid &) = delete;
+    CTable2Grid &operator=(const CTable2Grid &) = delete;
+
+  public:
+    CTable2Grid(PJ_CONTEXT *ctx, PAFile fp, int widthIn, int heightIn,
+                const ExtentAndRes &extentIn)
+        : HorizontalShiftGrid(widthIn, heightIn, extentIn), m_ctx(ctx),
+          m_fp(fp) {}
+
+    ~CTable2Grid() override;
+
+    bool valueAt(int, int, float &lonShift, float &latShift) const override;
+
+    static CTable2Grid *open(PJ_CONTEXT *ctx, PAFile fp,
+                             const std::string &filename);
+};
+
+// ---------------------------------------------------------------------------
+
+CTable2Grid::~CTable2Grid() { pj_ctx_fclose(m_ctx, m_fp); }
+
+// ---------------------------------------------------------------------------
+
+CTable2Grid *CTable2Grid::open(PJ_CONTEXT *ctx, PAFile fp,
+                               const std::string &filename) {
+    unsigned char header[160];
+
+    /* -------------------------------------------------------------------- */
+    /*      Read the header.                                                */
+    /* -------------------------------------------------------------------- */
+    if (pj_ctx_fread(ctx, header, sizeof(header), 1, fp) != 1) {
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return nullptr;
+    }
+
+    /* -------------------------------------------------------------------- */
+    /*      Regularize fields of interest.                                  */
+    /* -------------------------------------------------------------------- */
+    if (!IS_LSB) {
+        swap_words(header + 96, sizeof(double), 4);
+        swap_words(header + 128, sizeof(int), 2);
+    }
+
+    ExtentAndRes extent;
+    static_assert(sizeof(extent.westLon) == 8, "wrong sizeof");
+    static_assert(sizeof(extent.southLat) == 8, "wrong sizeof");
+    static_assert(sizeof(extent.resLon) == 8, "wrong sizeof");
+    static_assert(sizeof(extent.resLat) == 8, "wrong sizeof");
+    memcpy(&extent.westLon, header + 96, 8);
+    memcpy(&extent.southLat, header + 104, 8);
+    memcpy(&extent.resLon, header + 112, 8);
+    memcpy(&extent.resLat, header + 120, 8);
+    if (!(fabs(extent.westLon) <= 4 * M_PI &&
+          fabs(extent.southLat) <= M_PI + 1e-5 && extent.resLon > 1e-10 &&
+          extent.resLat > 1e-10)) {
+        pj_log(ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s",
+               filename.c_str());
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return nullptr;
+    }
+    int width;
+    int height;
+    memcpy(&width, header + 128, 4);
+    memcpy(&height, header + 132, 4);
+    if (width <= 0 || height <= 0) {
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return nullptr;
+    }
+    extent.eastLon = extent.westLon + (width - 1) * extent.resLon;
+    extent.northLat = extent.southLat + (height - 1) * extent.resLon;
+
+    return new CTable2Grid(ctx, fp, width, height, extent);
+}
+
+// ---------------------------------------------------------------------------
+
+bool CTable2Grid::valueAt(int x, int y, float &lonShift,
+                          float &latShift) const {
+    assert(x >= 0 && y >= 0 && x < m_width && y < m_height);
+
+    float two_floats[2];
+    pj_ctx_fseek(m_ctx, m_fp, 160 + 2 * sizeof(float) * (y * m_width + x),
+                 SEEK_SET);
+    if (pj_ctx_fread(m_ctx, &two_floats[0], sizeof(two_floats), 1, m_fp) != 1) {
+        pj_ctx_set_errno(m_ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return false;
+    }
+    if (!IS_LSB) {
+        swap_words(&two_floats[0], sizeof(float), 2);
+    }
+
+    latShift = two_floats[1];
+    lonShift = two_floats[0];
+
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+class NTv2GridSet : public HorizontalShiftGridSet {
+    PJ_CONTEXT *m_ctx;
+    PAFile m_fp;
+
+    NTv2GridSet(const NTv2GridSet &) = delete;
+    NTv2GridSet &operator=(const NTv2GridSet &) = delete;
+
+    NTv2GridSet(PJ_CONTEXT *ctx, PAFile fp) : m_ctx(ctx), m_fp(fp) {}
+
+  public:
+    ~NTv2GridSet() override;
+
+    static std::unique_ptr<NTv2GridSet> open(PJ_CONTEXT *ctx, PAFile fp,
+                                             const std::string &filename);
+};
+
+// ---------------------------------------------------------------------------
+
+class NTv2Grid : public HorizontalShiftGrid {
+    friend class NTv2GridSet;
+
+    std::string m_name;
+    PJ_CONTEXT *m_ctx; // owned by the parent NTv2GridSet
+    PAFile m_fp;       // owned by the parent NTv2GridSet
+    unsigned long long m_offset;
+    bool m_mustSwap;
+
+    NTv2Grid(const NTv2Grid &) = delete;
+    NTv2Grid &operator=(const NTv2Grid &) = delete;
+
+  public:
+    NTv2Grid(const std::string &nameIn, PJ_CONTEXT *ctx, PAFile fp,
+             unsigned long long offsetIn, bool mustSwapIn, int widthIn,
+             int heightIn, const ExtentAndRes &extentIn)
+        : HorizontalShiftGrid(widthIn, heightIn, extentIn), m_name(nameIn),
+          m_ctx(ctx), m_fp(fp), m_offset(offsetIn), m_mustSwap(mustSwapIn) {}
+
+    bool valueAt(int, int, float &lonShift, float &latShift) const override;
+};
+
+// ---------------------------------------------------------------------------
+
+bool NTv2Grid::valueAt(int x, int y, float &lonShift, float &latShift) const {
+    assert(x >= 0 && y >= 0 && x < m_width && y < m_height);
+
+    float two_float[2];
+    // NTv2 is organized from east to west !
+    // there are 4 components: lat shift, lon shift, lat error, lon error
+    pj_ctx_fseek(
+        m_ctx, m_fp,
+        // FIXME when fseek support unsigned long long
+        static_cast<long>(m_offset +
+                          4 * sizeof(float) * (y * m_width + m_width - 1 - x)),
+        SEEK_SET);
+    if (pj_ctx_fread(m_ctx, &two_float[0], sizeof(two_float), 1, m_fp) != 1) {
+        pj_ctx_set_errno(m_ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return false;
+    }
+    if (m_mustSwap) {
+        swap_words(&two_float[0], sizeof(float), 2);
+    }
+    /* convert seconds to radians */
+    latShift = static_cast<float>(two_float[0] * ((M_PI / 180.0) / 3600.0));
+    lonShift = static_cast<float>(two_float[1] * ((M_PI / 180.0) / 3600.0));
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+NTv2GridSet::~NTv2GridSet() { pj_ctx_fclose(m_ctx, m_fp); }
+
+// ---------------------------------------------------------------------------
+
+std::unique_ptr<NTv2GridSet> NTv2GridSet::open(PJ_CONTEXT *ctx, PAFile fp,
+                                               const std::string &filename) {
+    auto set = std::unique_ptr<NTv2GridSet>(new NTv2GridSet(ctx, fp));
+    set->m_name = filename;
+    set->m_format = "ntv2";
+
+    char header[11 * 16];
+
+    /* -------------------------------------------------------------------- */
+    /*      Read the header.                                                */
+    /* -------------------------------------------------------------------- */
+    if (pj_ctx_fread(ctx, header, sizeof(header), 1, fp) != 1) {
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return nullptr;
+    }
+
+    const bool must_swap = (header[8] == 11) ? !IS_LSB : IS_LSB;
+    if (must_swap) {
+        // swap_words( header+8, 4, 1 );
+        // swap_words( header+8+16, 4, 1 );
+        swap_words(header + 8 + 32, 4, 1);
+        // swap_words( header+8+7*16, 8, 1 );
+        // swap_words( header+8+8*16, 8, 1 );
+        // swap_words( header+8+9*16, 8, 1 );
+        // swap_words( header+8+10*16, 8, 1 );
+    }
+
+    /* -------------------------------------------------------------------- */
+    /*      Get the subfile count out ... all we really use for now.        */
+    /* -------------------------------------------------------------------- */
+    unsigned int num_subfiles;
+    memcpy(&num_subfiles, header + 8 + 32, 4);
+
+    std::map<std::string, NTv2Grid *> mapGrids;
+
+    /* ==================================================================== */
+    /*      Step through the subfiles, creating a grid for each.            */
+    /* ==================================================================== */
+    for (unsigned subfile = 0; subfile < num_subfiles; subfile++) {
+        // Read header
+        if (pj_ctx_fread(ctx, header, sizeof(header), 1, fp) != 1) {
+            pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+            return nullptr;
+        }
+
+        if (strncmp(header, "SUB_NAME", 8) != 0) {
+            pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+            return nullptr;
+        }
+
+        // Byte swap interesting fields if needed.
+        if (must_swap) {
+            swap_words(header + 8 + 16 * 4, sizeof(double), 6);
+            swap_words(header + 8 + 16 * 10, sizeof(int), 1);
+        }
+
+        std::string gridName;
+        gridName.append(header + 8, 8);
+
+        ExtentAndRes extent;
+        extent.westLon =
+            -to_double(header + 7 * 16 + 8) * DEG_TO_RAD / 3600.0; /* W_LONG */
+        extent.southLat =
+            to_double(header + 4 * 16 + 8) * DEG_TO_RAD / 3600.0; /* S_LAT */
+        extent.eastLon =
+            -to_double(header + 6 * 16 + 8) * DEG_TO_RAD / 3600.0; /* E_LONG */
+        extent.northLat =
+            to_double(header + 5 * 16 + 8) * DEG_TO_RAD / 3600.0; /* N_LAT */
+        extent.resLon = to_double(header + 9 * 16 + 8) * DEG_TO_RAD / 3600.0;
+        extent.resLat = to_double(header + 8 * 16 + 8) * DEG_TO_RAD / 3600.0;
+
+        if (!(fabs(extent.westLon) <= 4 * M_PI &&
+              fabs(extent.eastLon) <= 4 * M_PI &&
+              fabs(extent.northLat) <= M_PI + 1e-5 &&
+              fabs(extent.southLat) <= M_PI + 1e-5 &&
+              extent.westLon < extent.eastLon &&
+              extent.southLat < extent.northLat && extent.resLon > 1e-10 &&
+              extent.resLat > 1e-10)) {
+            pj_log(ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s",
+                   filename.c_str());
+            pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+            return nullptr;
+        }
+        const int columns = static_cast<int>(
+            fabs((extent.eastLon - extent.westLon) / extent.resLon + 0.5) + 1);
+        const int rows = static_cast<int>(
+            fabs((extent.northLat - extent.southLat) / extent.resLat + 0.5) +
+            1);
+
+        pj_log(ctx, PJ_LOG_DEBUG_MINOR,
+               "NTv2 %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)", gridName.c_str(),
+               columns, rows, extent.westLon * RAD_TO_DEG,
+               extent.southLat * RAD_TO_DEG, extent.eastLon * RAD_TO_DEG,
+               extent.northLat * RAD_TO_DEG);
+
+        unsigned int gs_count;
+        memcpy(&gs_count, header + 8 + 16 * 10, 4);
+        if (gs_count / columns != static_cast<unsigned>(rows)) {
+            pj_log(ctx, PJ_LOG_ERROR,
+                   "GS_COUNT(%u) does not match expected cells (%dx%d)",
+                   gs_count, columns, rows);
+            pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+            return nullptr;
+        }
+
+        auto offset = pj_ctx_ftell(ctx, fp);
+        auto grid = std::unique_ptr<NTv2Grid>(new NTv2Grid(
+            gridName, ctx, fp, offset, must_swap, columns, rows, extent));
+        std::string parentName;
+        parentName.assign(header + 24, 8);
+        auto iter = mapGrids.find(parentName);
+        auto gridPtr = grid.get();
+        if (iter == mapGrids.end()) {
+            set->m_grids.emplace_back(std::move(grid));
+        } else {
+            iter->second->m_children.emplace_back(std::move(grid));
+        }
+        mapGrids[gridName] = gridPtr;
+
+        // Skip grid data. 4 components of size float
+        pj_ctx_fseek(ctx, fp, gs_count * 4 * 4, SEEK_CUR);
+    }
+    return set;
+}
+
+// ---------------------------------------------------------------------------
+
+std::unique_ptr<HorizontalShiftGridSet>
+HorizontalShiftGridSet::open(PJ_CONTEXT *ctx, const std::string &filename) {
+    if (filename == "null") {
+        auto set = std::unique_ptr<HorizontalShiftGridSet>(
+            new HorizontalShiftGridSet());
+        set->m_name = filename;
+        set->m_format = "null";
+        set->m_grids.push_back(std::unique_ptr<NullHorizontalShiftGrid>(
+            new NullHorizontalShiftGrid()));
+        return set;
+    }
+
+    PAFile fp;
+    if (!(fp = pj_open_lib(ctx, filename.c_str(), "rb"))) {
+        ctx->last_errno = 0; /* don't treat as a persistent error */
+        return nullptr;
+    }
+
+    char header[160];
+    /* -------------------------------------------------------------------- */
+    /*      Load a header, to determine the file type.                      */
+    /* -------------------------------------------------------------------- */
+    size_t header_size;
+    if ((header_size = pj_ctx_fread(ctx, header, 1, sizeof(header), fp)) !=
+        sizeof(header)) {
+        /* some files may be smaller that sizeof(header), eg 160, so */
+        ctx->last_errno = 0; /* don't treat as a persistent error */
+        pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+               "pj_gridinfo_init: short header read of %d bytes",
+               (int)header_size);
+    }
+
+    pj_ctx_fseek(ctx, fp, SEEK_SET, 0);
+
+    /* -------------------------------------------------------------------- */
+    /*      Determine file type.                                            */
+    /* -------------------------------------------------------------------- */
+    if (header_size >= 144 + 16 && strncmp(header + 0, "HEADER", 6) == 0 &&
+        strncmp(header + 96, "W GRID", 6) == 0 &&
+        strncmp(header + 144, "TO      NAD83   ", 16) == 0) {
+        auto grid = NTv1Grid::open(ctx, fp, filename);
+        if (!grid) {
+            pj_ctx_fclose(ctx, fp);
+            return nullptr;
+        }
+        auto set = std::unique_ptr<HorizontalShiftGridSet>(
+            new HorizontalShiftGridSet());
+        set->m_name = filename;
+        set->m_format = "ntv1";
+        set->m_grids.push_back(std::unique_ptr<HorizontalShiftGrid>(grid));
+        return set;
+    } else if (header_size >= 9 && strncmp(header + 0, "CTABLE V2", 9) == 0) {
+        auto grid = CTable2Grid::open(ctx, fp, filename);
+        if (!grid) {
+            pj_ctx_fclose(ctx, fp);
+            return nullptr;
+        }
+        auto set = std::unique_ptr<HorizontalShiftGridSet>(
+            new HorizontalShiftGridSet());
+        set->m_name = filename;
+        set->m_format = "ctable2";
+        set->m_grids.push_back(std::unique_ptr<HorizontalShiftGrid>(grid));
+        return set;
+    } else if (header_size >= 48 + 7 &&
+               strncmp(header + 0, "NUM_OREC", 8) == 0 &&
+               strncmp(header + 48, "GS_TYPE", 7) == 0) {
+        return NTv2GridSet::open(ctx, fp, filename);
+    }
+
+    pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Unrecognized horizontal grid format");
+    pj_ctx_fclose(ctx, fp);
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+const HorizontalShiftGrid *HorizontalShiftGrid::gridAt(double lon,
+                                                       double lat) const {
+    for (const auto &child : m_children) {
+        const auto &extentChild = child->extentAndRes();
+        const double epsilon =
+            (extentChild.resLon + extentChild.resLat) / 10000.0;
+        if ((extentChild.fullWorldLongitude() ||
+             (lon + epsilon >= extentChild.westLon &&
+              lon - epsilon <= extentChild.eastLon)) &&
+            lat + epsilon >= extentChild.southLat &&
+            lat - epsilon <= extentChild.northLat) {
+            return child->gridAt(lon, lat);
+        }
+    }
+    return this;
+}
+// ---------------------------------------------------------------------------
+
+const HorizontalShiftGrid *HorizontalShiftGridSet::gridAt(double lon,
+                                                          double lat) const {
+    for (const auto &grid : m_grids) {
+        if (dynamic_cast<NullHorizontalShiftGrid *>(grid.get())) {
+            return grid.get();
+        }
+        const auto &extent = grid->extentAndRes();
+        const double epsilon = (extent.resLon + extent.resLat) / 10000.0;
+        if ((extent.fullWorldLongitude() ||
+             (lon + epsilon >= extent.westLon &&
+              lon - epsilon <= extent.eastLon)) &&
+            lat + epsilon >= extent.southLat &&
+            lat - epsilon <= extent.northLat) {
+            return grid->gridAt(lon, lat);
+        }
+    }
+    return nullptr;
+}
+
+NS_PROJ_END
diff --git a/src/grids.hpp b/src/grids.hpp
new file mode 100644
index 00000000..65bf502b
--- /dev/null
+++ b/src/grids.hpp
@@ -0,0 +1,167 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Grid management
+ * 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.
+ *****************************************************************************/
+
+#ifndef GRIDS_HPP_INCLUDED
+#define GRIDS_HPP_INCLUDED
+
+#include <memory>
+#include <vector>
+
+#include "proj.h"
+#include "proj/util.hpp"
+
+NS_PROJ_START
+
+struct ExtentAndRes {
+    double westLon;  // in radian
+    double southLat; // in radian
+    double eastLon;  // in radian
+    double northLat; // in radian
+    double resLon;   // in radian
+    double resLat;   // in radian
+
+    bool fullWorldLongitude() const;
+};
+
+// ---------------------------------------------------------------------------
+
+class Grid {
+  protected:
+    int m_width;
+    int m_height;
+    ExtentAndRes m_extent;
+
+    Grid(int widthIn, int heightIn, const ExtentAndRes &extentIn);
+
+  public:
+    virtual ~Grid();
+
+    int width() const { return m_width; }
+    int height() const { return m_height; }
+    const ExtentAndRes &extentAndRes() const { return m_extent; }
+
+    virtual bool isNullGrid() const { return false; }
+};
+
+// ---------------------------------------------------------------------------
+
+class VerticalShiftGrid : public Grid {
+  protected:
+    std::vector<std::unique_ptr<VerticalShiftGrid>> m_children{};
+
+  public:
+    VerticalShiftGrid(int widthIn, int heightIn, const ExtentAndRes &extentIn);
+
+    const VerticalShiftGrid *gridAt(double lon, double lat) const;
+
+    virtual bool isNodata(float /*val*/, double /* multiplier */) const;
+
+    // x = 0 is western-most column, y = 0 is southern-most line
+    virtual bool valueAt(int x, int y, float &out) const = 0;
+};
+
+// ---------------------------------------------------------------------------
+
+class VerticalShiftGridSet {
+    std::string m_name{};
+    std::string m_format{};
+    std::vector<std::unique_ptr<VerticalShiftGrid>> m_grids{};
+
+    VerticalShiftGridSet();
+
+  public:
+    virtual ~VerticalShiftGridSet();
+
+    static std::unique_ptr<VerticalShiftGridSet>
+    open(PJ_CONTEXT *ctx, const std::string &filename);
+
+    const std::string &name() const { return m_name; }
+    const std::string &format() const { return m_format; }
+    const std::vector<std::unique_ptr<VerticalShiftGrid>> &grids() const {
+        return m_grids;
+    }
+    const VerticalShiftGrid *gridAt(double lon, double lat) const;
+};
+
+// ---------------------------------------------------------------------------
+
+class HorizontalShiftGrid : public Grid {
+  protected:
+    std::vector<std::unique_ptr<HorizontalShiftGrid>> m_children{};
+
+  public:
+    HorizontalShiftGrid(int widthIn, int heightIn,
+                        const ExtentAndRes &extentIn);
+    ~HorizontalShiftGrid() override;
+
+    const HorizontalShiftGrid *gridAt(double lon, double lat) const;
+
+    // x = 0 is western-most column, y = 0 is southern-most line
+    virtual bool valueAt(int x, int y, float &lonShift,
+                         float &latShift) const = 0;
+};
+
+// ---------------------------------------------------------------------------
+
+class HorizontalShiftGridSet {
+  protected:
+    std::string m_name{};
+    std::string m_format{};
+    std::vector<std::unique_ptr<HorizontalShiftGrid>> m_grids{};
+
+    HorizontalShiftGridSet();
+
+  public:
+    virtual ~HorizontalShiftGridSet();
+
+    static std::unique_ptr<HorizontalShiftGridSet>
+    open(PJ_CONTEXT *ctx, const std::string &filename);
+
+    const std::string &name() const { return m_name; }
+    const std::string &format() const { return m_format; }
+    const std::vector<std::unique_ptr<HorizontalShiftGrid>> &grids() const {
+        return m_grids;
+    }
+    const HorizontalShiftGrid *gridAt(double lon, double lat) const;
+};
+
+// ---------------------------------------------------------------------------
+
+typedef std::vector<std::unique_ptr<HorizontalShiftGridSet>> ListOfHGrids;
+typedef std::vector<std::unique_ptr<VerticalShiftGridSet>> ListOfVGrids;
+
+ListOfVGrids proj_vgrid_init(PJ *P, const char *grids);
+ListOfHGrids proj_hgrid_init(PJ *P, const char *grids);
+double proj_vgrid_value(PJ *P, const ListOfVGrids &, PJ_LP lp,
+                        double vmultiplier);
+PJ_LP proj_hgrid_value(PJ *P, const ListOfHGrids &, PJ_LP lp);
+PJ_LP proj_hgrid_apply(PJ *P, const ListOfHGrids &, PJ_LP lp,
+                       PJ_DIRECTION direction);
+
+NS_PROJ_END
+
+#endif // GRIDS_HPP_INCLUDED
diff --git a/src/init.cpp b/src/init.cpp
index 1ed46e5a..19fcf47b 100644
--- a/src/init.cpp
+++ b/src/init.cpp
@@ -646,12 +646,6 @@ pj_init_ctx_with_allow_init_epsg(projCtx ctx, int argc, char **argv, int allow_i
     PIN->long_wrap_center = 0.0;
     strcpy( PIN->axis, "enu" );
 
-    PIN->gridlist = nullptr;
-    PIN->gridlist_count = 0;
-
-    PIN->vgridlist_geoid = nullptr;
-    PIN->vgridlist_geoid_count = 0;
-
     /* Set datum parameters. Similarly to +init parameters we want to expand    */
     /* +datum parameters as late as possible when dealing with pipelines.       */
     /* otherwise only the first occurrence of +datum will be expanded and that */
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index d1bc8836..f6112aef 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -234,13 +234,9 @@ set(SRC_LIBPROJ_CORE
   fileapi.cpp
   fwd.cpp
   gauss.cpp
-  gc_reader.cpp
   geocent.cpp
   geocent.h
   geodesic.c
-  gridcatalog.cpp
-  gridinfo.cpp
-  gridlist.cpp
   init.cpp
   initcache.cpp
   internal.cpp
@@ -252,9 +248,6 @@ set(SRC_LIBPROJ_CORE
   mlfn.cpp
   msfn.cpp
   mutex.cpp
-  nad_cvt.cpp
-  nad_init.cpp
-  nad_intr.cpp
   open_lib.cpp
   param.cpp
   phi2.cpp
@@ -286,6 +279,8 @@ set(SRC_LIBPROJ_CORE
   proj_json_streaming_writer.hpp
   proj_json_streaming_writer.cpp
   tracing.cpp
+  grids.hpp
+  grids.cpp
   ${CMAKE_CURRENT_BINARY_DIR}/proj_config.h
 )
 
diff --git a/src/malloc.cpp b/src/malloc.cpp
index 393437e3..3fd3699f 100644
--- a/src/malloc.cpp
+++ b/src/malloc.cpp
@@ -48,6 +48,9 @@
 
 #include "proj.h"
 #include "proj_internal.h"
+#include "grids.hpp"
+
+using namespace NS_PROJ;
 
 /**********************************************************************/
 void *pj_malloc(size_t size) {
@@ -225,10 +228,8 @@ PJ *pj_default_destructor (PJ *P, int errlev) {   /* Destructor */
     pj_dealloc(P->def_spherification);
     pj_dealloc(P->def_ellps);
 
-    /* free grid lists */
-    pj_dealloc( P->gridlist );
-    pj_dealloc( P->vgridlist_geoid );
-    pj_dealloc( P->catalog_name );
+    delete static_cast<ListOfHGrids*>(P->hgrids_legacy);
+    delete static_cast<ListOfVGrids*>(P->vgrids_legacy);
 
     /* We used to call pj_dalloc( P->catalog ), but this will leak */
     /* memory. The safe way to clear catalog and grid is to call */
diff --git a/src/nad_cvt.cpp b/src/nad_cvt.cpp
deleted file mode 100644
index 79441d0a..00000000
--- a/src/nad_cvt.cpp
+++ /dev/null
@@ -1,77 +0,0 @@
-#define PJ_LIB__
-
-#include <stdio.h>
-#include <stdlib.h>
-
-#include "proj.h"
-#include "proj_internal.h"
-#include <math.h>
-
-#define MAX_ITERATIONS 10
-#define TOL 1e-12
-
-PJ_LP nad_cvt(PJ_LP in, int inverse, struct CTABLE *ct) {
-    PJ_LP t, tb,del, dif;
-    int i = MAX_ITERATIONS;
-    const double toltol = TOL*TOL;
-
-    if (in.lam == HUGE_VAL)
-        return in;
-
-    /* normalize input to ll origin */
-    tb = in;
-    tb.lam -= ct->ll.lam;
-    tb.phi -= ct->ll.phi;
-    tb.lam = adjlon (tb.lam - M_PI) + M_PI;
-
-    t = nad_intr (tb, ct);
-    if (t.lam == HUGE_VAL)
-        return t;
-
-    if (!inverse) {
-        in.lam -= t.lam;
-        in.phi += t.phi;
-        return in;
-    }
-
-    t.lam = tb.lam + t.lam;
-    t.phi = tb.phi - t.phi;
-
-    do {
-        del = nad_intr(t, ct);
-
-        /* This case used to return failure, but I have
-           changed it to return the first order approximation
-           of the inverse shift.  This avoids cases where the
-           grid shift *into* this grid came from another grid.
-           While we aren't returning optimally correct results
-           I feel a close result in this case is better than
-           no result.  NFW
-           To demonstrate use -112.5839956 49.4914451 against
-           the NTv2 grid shift file from Canada. */
-        if (del.lam == HUGE_VAL)
-            break;
-
-        dif.lam = t.lam - del.lam - tb.lam;
-           dif.phi = t.phi + del.phi - tb.phi;
-        t.lam -= dif.lam;
-        t.phi -= dif.phi;
-
-    } while (--i && (dif.lam*dif.lam + dif.phi*dif.phi > toltol)); /* prob. slightly faster than hypot() */
-
-    if (i==0) {
-        /* If we had access to a context, this should go through pj_log, and we should set ctx->errno */
-        if (getenv ("PROJ_DEBUG"))
-            fprintf( stderr, "Inverse grid shift iterator failed to converge.\n" );
-        t.lam = t.phi = HUGE_VAL;
-        return t;
-    }
-
-    /* and again: pj_log and ctx->errno */
-    if (del.lam==HUGE_VAL && getenv ("PROJ_DEBUG"))
-        fprintf (stderr, "Inverse grid shift iteration failed, presumably at grid edge.\nUsing first approximation.\n");
-
-    in.lam = adjlon (t.lam + ct->ll.lam);
-    in.phi = t.phi + ct->ll.phi;
-    return in;
-}
diff --git a/src/nad_init.cpp b/src/nad_init.cpp
deleted file mode 100644
index 315318be..00000000
--- a/src/nad_init.cpp
+++ /dev/null
@@ -1,308 +0,0 @@
-/******************************************************************************
- * Project:  PROJ.4
- * Purpose:  Load datum shift files into memory.
- * Author:   Frank Warmerdam, warmerdam@pobox.com
- *
- ******************************************************************************
- * Copyright (c) 2000, Frank Warmerdam
- *
- * 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 <stddef.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include "proj.h"
-#include "proj_internal.h"
-
-/************************************************************************/
-/*                             swap_words()                             */
-/*                                                                      */
-/*      Convert the byte order of the given word(s) in place.           */
-/************************************************************************/
-
-static const int  byte_order_test = 1;
-#define IS_LSB	(((const unsigned char *) (&byte_order_test))[0] == 1)
-
-static void swap_words( void *data_in, int word_size, int word_count )
-
-{
-    int	word;
-    unsigned char *data = (unsigned char *) data_in;
-
-    for( word = 0; word < word_count; word++ )
-    {
-        int	i;
-        
-        for( i = 0; i < word_size/2; i++ )
-        {
-            unsigned char	t;
-            
-            t = data[i];
-            data[i] = data[word_size-i-1];
-            data[word_size-i-1] = t;
-        }
-        
-        data += word_size;
-    }
-}
-
-/************************************************************************/
-/*                          nad_ctable_load()                           */
-/*                                                                      */
-/*      Load the data portion of a ctable formatted grid.               */
-/************************************************************************/
-
-int nad_ctable_load( projCtx ctx, struct CTABLE *ct, struct projFileAPI_t* fileapi )
-
-{
-    PAFile fid = (PAFile)fileapi;
-    size_t a_size;
-
-    pj_ctx_fseek( ctx, fid, sizeof(struct CTABLE), SEEK_SET );
-
-    /* read all the actual shift values */
-    a_size = ct->lim.lam * ct->lim.phi;
-    ct->cvs = (FLP *) pj_malloc(sizeof(FLP) * a_size);
-    if( ct->cvs == nullptr 
-        || pj_ctx_fread(ctx, ct->cvs, sizeof(FLP), a_size, fid) != a_size )
-    {
-        pj_dalloc( ct->cvs );
-        ct->cvs = nullptr;
-
-        pj_log( ctx, PJ_LOG_ERROR, 
-                "ctable loading failed on fread() - binary incompatible?" );
-        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        return 0;
-    }
-
-    return 1;
-} 
-
-/************************************************************************/
-/*                          nad_ctable_init()                           */
-/*                                                                      */
-/*      Read the header portion of a "ctable" format grid.              */
-/************************************************************************/
-
-struct CTABLE *nad_ctable_init( projCtx ctx, struct projFileAPI_t* fileapi )
-{
-    PAFile fid = (PAFile)fileapi;
-    struct CTABLE *ct;
-    int		id_end;
-
-    /* read the table header */
-    ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE));
-    if( ct == nullptr 
-        || pj_ctx_fread( ctx, ct, sizeof(struct CTABLE), 1, fid ) != 1 )
-    {
-        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        pj_dalloc( ct );
-        return nullptr;
-    }
-
-    /* do some minimal validation to ensure the structure isn't corrupt */
-    if( ct->lim.lam < 1 || ct->lim.lam > 100000 
-        || ct->lim.phi < 1 || ct->lim.phi > 100000 )
-    {
-        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        pj_dalloc( ct );
-        return nullptr;
-    }
-    
-    /* trim white space and newlines off id */
-    for( id_end = (int)strlen(ct->id)-1; id_end > 0; id_end-- )
-    {
-        if( ct->id[id_end] == '\n' || ct->id[id_end] == ' ' )
-            ct->id[id_end] = '\0';
-        else
-            break;
-    }
-
-    ct->cvs = nullptr;
-
-    return ct;
-}
-
-/************************************************************************/
-/*                          nad_ctable2_load()                          */
-/*                                                                      */
-/*      Load the data portion of a ctable2 formatted grid.              */
-/************************************************************************/
-
-int nad_ctable2_load( projCtx ctx, struct CTABLE *ct, struct projFileAPI_t* fileapi )
-
-{
-    PAFile fid = (PAFile)fileapi;
-    size_t a_size;
-
-    pj_ctx_fseek( ctx, fid, 160, SEEK_SET );
-
-    /* read all the actual shift values */
-    a_size = ct->lim.lam * ct->lim.phi;
-    ct->cvs = (FLP *) pj_malloc(sizeof(FLP) * a_size);
-    if( ct->cvs == nullptr 
-        || pj_ctx_fread(ctx, ct->cvs, sizeof(FLP), a_size, fid) != a_size )
-    {
-        pj_dalloc( ct->cvs );
-        ct->cvs = nullptr;
-
-        if( getenv("PROJ_DEBUG") != nullptr )
-        {
-            fprintf( stderr,
-            "ctable2 loading failed on fread() - binary incompatible?\n" );
-        }
-
-        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        return 0;
-    }
-
-    if( !IS_LSB )
-    {
-        swap_words( ct->cvs, 4, (int)a_size * 2 );
-    }
-
-    return 1;
-} 
-
-/************************************************************************/
-/*                          nad_ctable2_init()                          */
-/*                                                                      */
-/*      Read the header portion of a "ctable2" format grid.             */
-/************************************************************************/
-
-struct CTABLE *nad_ctable2_init( projCtx ctx, struct projFileAPI_t* fileapi )
-{
-    PAFile fid = (PAFile)fileapi;
-    struct CTABLE *ct;
-    int		id_end;
-    char        header[160];
-
-    if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 )
-    {
-        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        return nullptr;
-    }
-
-    if( !IS_LSB )
-    {
-        swap_words( header +  96, 8, 4 );
-        swap_words( header + 128, 4, 2 );
-    }
-
-    if( strncmp(header,"CTABLE V2",9) != 0 )
-    {
-        pj_log( ctx, PJ_LOG_ERROR, "ctable2 - wrong header!" );
-        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        return nullptr;
-    }
-
-    /* read the table header */
-    ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE));
-    if( ct == nullptr )
-    {
-        pj_ctx_set_errno( ctx, ENOMEM );
-        return nullptr;
-    }
-
-    memcpy( ct->id,       header +  16, 80 );
-    memcpy( &ct->ll.lam,  header +  96, 8 );
-    memcpy( &ct->ll.phi,  header + 104, 8 );
-    memcpy( &ct->del.lam, header + 112, 8 );
-    memcpy( &ct->del.phi, header + 120, 8 );
-    memcpy( &ct->lim.lam, header + 128, 4 );
-    memcpy( &ct->lim.phi, header + 132, 4 );
-
-    /* do some minimal validation to ensure the structure isn't corrupt */
-    if( ct->lim.lam < 1 || ct->lim.lam > 100000 
-        || ct->lim.phi < 1 || ct->lim.phi > 100000 )
-    {
-        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        pj_dalloc( ct );
-        return nullptr;
-    }
-    
-    /* trim white space and newlines off id */
-    for( id_end = (int)strlen(ct->id)-1; id_end > 0; id_end-- )
-    {
-        if( ct->id[id_end] == '\n' || ct->id[id_end] == ' ' )
-            ct->id[id_end] = '\0';
-        else
-            break;
-    }
-
-    ct->cvs = nullptr;
-
-    return ct;
-}
-
-/************************************************************************/
-/*                              nad_init()                              */
-/*                                                                      */
-/*      Read a datum shift file in any of the supported binary formats. */
-/************************************************************************/
-
-struct CTABLE *nad_init(projCtx ctx, char *name) 
-{
-    struct CTABLE *ct;
-    PAFile      fid;
-
-    ctx->last_errno = 0;
-
-/* -------------------------------------------------------------------- */
-/*      Open the file using the usual search rules.                     */
-/* -------------------------------------------------------------------- */
-    if (!(fid = pj_open_lib(ctx, name, "rb"))) {
-        return nullptr;
-    }
-
-    ct = nad_ctable_init( ctx, (struct projFileAPI_t*)fid );
-    if( ct != nullptr )
-    {
-        if( !nad_ctable_load( ctx, ct, (struct projFileAPI_t*)fid ) )
-        {
-            nad_free( ct );
-            ct = nullptr;
-        }
-    }
-
-    pj_ctx_fclose(ctx, fid);
-    return ct;
-}
-
-/************************************************************************/
-/*                              nad_free()                              */
-/*                                                                      */
-/*      Free a CTABLE grid shift structure produced by nad_init().      */
-/************************************************************************/
-
-void nad_free(struct CTABLE *ct) 
-{
-    if (ct) {
-        if( ct->cvs != nullptr )
-            pj_dalloc(ct->cvs);
-
-        pj_dalloc(ct);
-    }
-}
diff --git a/src/nad_intr.cpp b/src/nad_intr.cpp
deleted file mode 100644
index 36ab0a99..00000000
--- a/src/nad_intr.cpp
+++ /dev/null
@@ -1,67 +0,0 @@
-/* Determine nad table correction value */
-#define PJ_LIB__
-#include "proj.h"
-#include "proj_internal.h"
-#include <math.h>
-
-PJ_LP nad_intr(PJ_LP t, struct CTABLE *ct) {
-	PJ_LP val, frct;
-	ILP indx;
-	double m00, m10, m01, m11;
-	FLP *f00, *f10, *f01, *f11;
-	long index;
-	int in;
-
-	t.lam /= ct->del.lam;
-	indx.lam = isnan(t.lam) ? 0 : (pj_int32)lround(floor(t.lam));
-	t.phi /= ct->del.phi;
-	indx.phi = isnan(t.phi) ? 0 : (pj_int32)lround(floor(t.phi));
-
-	frct.lam = t.lam - indx.lam;
-	frct.phi = t.phi - indx.phi;
-	val.lam = val.phi = HUGE_VAL;
-	if (indx.lam < 0) {
-		if (indx.lam == -1 && frct.lam > 0.99999999999) {
-			++indx.lam;
-			frct.lam = 0.;
-		} else
-			return val;
-	} else if ((in = indx.lam + 1) >= ct->lim.lam) {
-		if (in == ct->lim.lam && frct.lam < 1e-11) {
-			--indx.lam;
-			frct.lam = 1.;
-		} else
-			return val;
-	}
-	if (indx.phi < 0) {
-		if (indx.phi == -1 && frct.phi > 0.99999999999) {
-			++indx.phi;
-			frct.phi = 0.;
-		} else
-			return val;
-	} else if ((in = indx.phi + 1) >= ct->lim.phi) {
-		if (in == ct->lim.phi && frct.phi < 1e-11) {
-			--indx.phi;
-			frct.phi = 1.;
-		} else
-			return val;
-	}
-	index = indx.phi * ct->lim.lam + indx.lam;
-	f00 = ct->cvs + index++;
-	f10 = ct->cvs + index;
-	index += ct->lim.lam;
-	f11 = ct->cvs + index--;
-	f01 = ct->cvs + index;
-	m11 = m10 = frct.lam;
-	m00 = m01 = 1. - frct.lam;
-	m11 *= frct.phi;
-	m01 *= frct.phi;
-	frct.phi = 1. - frct.phi;
-	m00 *= frct.phi;
-	m10 *= frct.phi;
-	val.lam = m00 * f00->lam + m10 * f10->lam +
-			  m01 * f01->lam + m11 * f11->lam;
-	val.phi = m00 * f00->phi + m10 * f10->phi +
-			  m01 * f01->phi + m11 * f11->phi;
-	return val;
-}
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 3ca927a3..7d826414 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -195,14 +195,6 @@ PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P);
 PJ_COORD PROJ_DLL pj_approx_2D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo);
 PJ_COORD PROJ_DLL pj_approx_3D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo);
 
-
-/* Grid functionality */
-int             proj_vgrid_init(PJ *P, const char *grids);
-int             proj_hgrid_init(PJ *P, const char *grids);
-double          proj_vgrid_value(PJ *P, PJ_LP lp, double vmultiplier);
-PJ_LP           proj_hgrid_value(PJ *P, PJ_LP lp);
-PJ_LP           proj_hgrid_apply(PJ *P, PJ_LP lp, PJ_DIRECTION direction);
-
 void PROJ_DLL proj_log_error (PJ *P, const char *fmt, ...);
 void proj_log_debug (PJ *P, const char *fmt, ...);
 void proj_log_trace (PJ *P, const char *fmt, ...);
@@ -478,32 +470,16 @@ struct PJconsts {
 
     int     datum_type = PJD_UNKNOWN;  /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */
     double  datum_params[7] = {0,0,0,0,0,0,0}; /* Parameters for 3PARAM and 7PARAM */
-    struct _pj_gi **gridlist = nullptr;     /* TODO: Description needed */
-    int     gridlist_count = 0;
 
-    int     has_geoid_vgrids = 0;      /* TODO: Description needed */
-    struct _pj_gi **vgridlist_geoid = nullptr;   /* TODO: Description needed */
-    int     vgridlist_geoid_count = 0;
+    int     has_geoid_vgrids = 0;      /* used by legacy transform.cpp */
+    void*   hgrids_legacy = nullptr;   /* used by legacy transform.cpp. Is a pointer to a ListOfHGrids* */ 
+    void*   vgrids_legacy = nullptr;   /* used by legacy transform.cpp. Is a pointer to a ListOfVGrids* */ 
 
     double  from_greenwich = 0.0;       /* prime meridian offset (in radians) */
     double  long_wrap_center = 0.0;     /* 0.0 for -180 to 180, actually in radians*/
     int     is_long_wrap_set = 0;
     char    axis[4] = {0,0,0,0};        /* Axis order, pj_transform/pj_adjust_axis */
 
-    /* New Datum Shift Grid Catalogs */
-    char   *catalog_name = nullptr;
-    struct _PJ_GridCatalog *catalog = nullptr;
-
-    double  datum_date = 0.0;           /* TODO: Description needed */
-
-    struct _pj_gi *last_before_grid = nullptr;    /* TODO: Description needed */
-    PJ_Region     last_before_region = {0,0,0,0}; /* TODO: Description needed */
-    double        last_before_date = 0.0;         /* TODO: Description needed */
-
-    struct _pj_gi *last_after_grid = nullptr;     /* TODO: Description needed */
-    PJ_Region     last_after_region = {0,0,0,0};  /* TODO: Description needed */
-    double        last_after_date = 0.0;          /* TODO: Description needed */
-
     /*************************************************************************************
      ISO-19111 interface
     **************************************************************************************/
@@ -766,56 +742,6 @@ PJ *pj_projection_specific_setup_##name (PJ *P)
 
 #endif /* def PJ_LIB__ */
 
-
-#define MAX_TAB_ID 80
-typedef struct { float lam, phi; } FLP;
-typedef struct { pj_int32 lam, phi; } ILP;
-
-struct CTABLE {
-    char id[MAX_TAB_ID];    /* ascii info */
-    PJ_LP ll;               /* lower left corner coordinates */
-    PJ_LP del;              /* size of cells */
-    ILP lim;                /* limits of conversion matrix */
-    FLP *cvs;               /* conversion matrix */
-};
-
-typedef struct _pj_gi {
-    char *gridname;     /* identifying name of grid, eg "conus" or ntv2_0.gsb */
-    char *filename;     /* full path to filename */
-
-    const char *format; /* format of this grid, ie "ctable", "ntv1",
-                           "ntv2" or "missing". */
-
-    long   grid_offset;  /* offset in file, for delayed loading */
-    int   must_swap;    /* only for NTv2 */
-
-    struct CTABLE *ct;
-
-    struct _pj_gi *next;
-    struct _pj_gi *child;
-} PJ_GRIDINFO;
-
-typedef struct {
-    PJ_Region region;
-    int  priority;      /* higher used before lower */
-    double date;        /* year.fraction */
-    char *definition;   /* usually the gridname */
-
-    PJ_GRIDINFO  *gridinfo;
-    int available;      /* 0=unknown, 1=true, -1=false */
-} PJ_GridCatalogEntry;
-
-typedef struct _PJ_GridCatalog {
-    char *catalog_name;
-
-    PJ_Region region;   /* maximum extent of catalog data */
-
-    int entry_count;
-    PJ_GridCatalogEntry *entries;
-
-    struct _PJ_GridCatalog *next;
-} PJ_GridCatalog;
-
 /* procedure prototypes */
 double PROJ_DLL dmstor(const char *, char **);
 double dmstor_ctx(projCtx_t *ctx, const char *, char **);
@@ -862,50 +788,6 @@ COMPLEX pj_zpolyd1(COMPLEX, const COMPLEX *, int, COMPLEX *);
 int pj_deriv(PJ_LP, double, const PJ *, struct DERIVS *);
 int pj_factors(PJ_LP, const PJ *, double, struct FACTORS *);
 
-/* nadcon related protos */
-PJ_LP             nad_intr(PJ_LP, struct CTABLE *);
-PJ_LP             nad_cvt(PJ_LP, int, struct CTABLE *);
-struct CTABLE *nad_init(projCtx_t *ctx, char *);
-struct CTABLE *nad_ctable_init( projCtx_t *ctx, struct projFileAPI_t* fid );
-int            nad_ctable_load( projCtx_t *ctx, struct CTABLE *, struct projFileAPI_t* fid );
-struct CTABLE *nad_ctable2_init( projCtx_t *ctx, struct projFileAPI_t* fid );
-int            nad_ctable2_load( projCtx_t *ctx, struct CTABLE *, struct projFileAPI_t* fid );
-void           nad_free(struct CTABLE *);
-
-/* higher level handling of datum grid shift files */
-
-int pj_apply_vgridshift( PJ *defn, const char *listname,
-                         PJ_GRIDINFO ***gridlist_p,
-                         int *gridlist_count_p,
-                         int inverse,
-                         long point_count, int point_offset,
-                         double *x, double *y, double *z );
-int pj_apply_gridshift_2( PJ *defn, int inverse,
-                          long point_count, int point_offset,
-                          double *x, double *y, double *z );
-int pj_apply_gridshift_3( projCtx_t *ctx,
-                          PJ_GRIDINFO **gridlist, int gridlist_count,
-                          int inverse, long point_count, int point_offset,
-                          double *x, double *y, double *z );
-
-PJ_GRIDINFO **pj_gridlist_from_nadgrids( projCtx_t *, const char *, int * );
-
-PJ_GRIDINFO *pj_gridinfo_init( projCtx_t *, const char * );
-int          pj_gridinfo_load( projCtx_t *, PJ_GRIDINFO * );
-void         pj_gridinfo_free( projCtx_t *, PJ_GRIDINFO * );
-
-PJ_GridCatalog *pj_gc_findcatalog( projCtx_t *, const char * );
-PJ_GridCatalog *pj_gc_readcatalog( projCtx_t *, const char * );
-void pj_gc_unloadall( projCtx_t *);
-int pj_gc_apply_gridshift( PJ *defn, int inverse,
-                           long point_count, int point_offset,
-                           double *x, double *y, double *z );
-int pj_gc_apply_gridshift( PJ *defn, int inverse,
-                           long point_count, int point_offset,
-                           double *x, double *y, double *z );
-
-double pj_gc_parsedate( projCtx_t *, const char * );
-
 void  *proj_mdist_ini(double);
 double proj_mdist(double, double, double, const void *);
 double proj_inv_mdist(projCtx_t *ctx, double, const void *);
diff --git a/src/transform.cpp b/src/transform.cpp
index d111d835..020d62ea 100644
--- a/src/transform.cpp
+++ b/src/transform.cpp
@@ -34,6 +34,9 @@
 #include "proj.h"
 #include "proj_internal.h"
 #include "geocent.h"
+#include "grids.hpp"
+
+using namespace NS_PROJ;
 
 static int adjust_axis( projCtx ctx, const char *axis, int denormalize_flag,
                            long point_count, int point_offset,
@@ -431,6 +434,78 @@ static int height_unit (PJ *P, PJ_DIRECTION dir, long n, int dist, double *z) {
 }
 
 
+/************************************************************************/
+/*                        pj_apply_vgridshift()                         */
+/*                                                                      */
+/*      This implementation takes uses the gridlist from a coordinate   */
+/*      system definition.  If the gridlist has not yet been            */
+/*      populated in the coordinate system definition we set it up      */
+/*      now.                                                            */
+/************************************************************************/
+static int pj_apply_vgridshift( PJ *defn,
+                         int inverse,
+                         long point_count, int point_offset,
+                         double *x, double *y, double *z )
+
+{
+    if( defn->vgrids_legacy == nullptr )
+    {
+        defn->vgrids_legacy = new ListOfVGrids;
+        auto vgrids = proj_vgrid_init(defn, "geoidgrids");
+        if( vgrids.empty() )
+            return 0;
+        *static_cast<ListOfVGrids*>(defn->vgrids_legacy) = std::move(vgrids);
+    }
+    if( static_cast<ListOfVGrids*>(defn->vgrids_legacy)->empty() )
+    {
+        return 0;
+    }
+
+    for( int i = 0; i < point_count; i++ )
+    {
+        double value;
+        long io = i * point_offset;
+        PJ_LP   input;
+
+        input.phi = y[io];
+        input.lam = x[io];
+
+        value = proj_vgrid_value(defn, *static_cast<ListOfVGrids*>(defn->vgrids_legacy), input, 1.0);
+
+        if( inverse )
+            z[io] -= value;
+        else
+            z[io] += value;
+
+        if( value == HUGE_VAL )
+        {
+            std::string gridlist;
+
+            proj_log_debug(defn,
+                "pj_apply_vgridshift(): failed to find a grid shift table for\n"
+                "                       location (%.7fdW,%.7fdN)",
+                x[io] * RAD_TO_DEG,
+                y[io] * RAD_TO_DEG );
+
+            for( const auto& gridset: *static_cast<ListOfVGrids*>(defn->vgrids_legacy) )
+            {
+                if( gridlist.empty() )
+                    gridlist += "   tried: ";
+                else
+                    gridlist += ',';
+                gridlist += gridset->name();
+            }
+
+            proj_log_debug(defn, "%s", gridlist.c_str());
+            pj_ctx_set_errno( defn->ctx, PJD_ERR_GRID_AREA );
+
+            return PJD_ERR_GRID_AREA;
+        }
+    }
+
+    return 0;
+}
+
 
 /* -------------------------------------------------------------------- */
 /*           Transform to ellipsoidal heights if needed                 */
@@ -441,10 +516,7 @@ static int geometric_to_orthometric (PJ *P, PJ_DIRECTION dir, long n, int dist,
         return 0;
     if (z==nullptr)
         return PJD_ERR_GEOCENTRIC;
-    err = pj_apply_vgridshift (P, "sgeoidgrids",
-              &(P->vgridlist_geoid),
-              &(P->vgridlist_geoid_count),
-              dir==PJ_FWD ? 1 : 0, n, dist, x, y, z );
+    err = pj_apply_vgridshift (P, dir==PJ_FWD ? 1 : 0, n, dist, x, y, z );
     if (err)
         return pj_ctx_get_errno(P->ctx);
     return 0;
@@ -806,6 +878,66 @@ int pj_geocentric_from_wgs84( PJ *defn,
     return 0;
 }
 
+
+/************************************************************************/
+/*                        pj_apply_gridshift_2()                        */
+/*                                                                      */
+/*      This implementation uses the gridlist from a coordinate         */
+/*      system definition.  If the gridlist has not yet been            */
+/*      populated in the coordinate system definition we set it up      */
+/*      now.                                                            */
+/************************************************************************/
+static
+int pj_apply_gridshift_2( PJ *defn, int inverse,
+                          long point_count, int point_offset,
+                          double *x, double *y, double * /*z*/ )
+
+{
+    if( defn->hgrids_legacy == nullptr )
+    {
+        defn->hgrids_legacy = new ListOfHGrids;
+        auto hgrids = proj_hgrid_init(defn, "nadgrids");
+        if( hgrids.empty() )
+            return 0;
+        *static_cast<ListOfHGrids*>(defn->hgrids_legacy) = std::move(hgrids);
+    }
+    if( static_cast<ListOfHGrids*>(defn->hgrids_legacy)->empty() )
+    {
+        return 0;
+    }
+
+    for( long i = 0; i < point_count; i++ )
+    {
+        PJ_LP   input;
+
+        long io = i * point_offset;
+        input.phi = y[io];
+        input.lam = x[io];
+
+        auto output = proj_hgrid_apply(defn, *static_cast<ListOfHGrids*>(defn->hgrids_legacy), input, inverse ? PJ_INV : PJ_FWD);
+
+        if ( output.lam != HUGE_VAL )
+        {
+            y[io] = output.phi;
+            x[io] = output.lam;
+        }
+        else
+        {
+            if( defn->ctx->debug_level >= PJ_LOG_DEBUG_MAJOR )
+            {
+                pj_log( defn->ctx, PJ_LOG_DEBUG_MAJOR,
+                    "pj_apply_gridshift(): failed to find a grid shift table for\n"
+                    "                      location (%.7fdW,%.7fdN)",
+                    x[io] * RAD_TO_DEG,
+                    y[io] * RAD_TO_DEG );
+            }
+        }
+    }
+
+    return 0;
+}
+
+
 /************************************************************************/
 /*                         pj_datum_transform()                         */
 /*                                                                      */
@@ -1046,3 +1178,8 @@ static int adjust_axis( projCtx ctx,
 
     return 0;
 }
+// ---------------------------------------------------------------------------
+
+void pj_deallocate_grids()
+{
+}
diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp
index f1311a54..f5468775 100644
--- a/src/transformations/deformation.cpp
+++ b/src/transformations/deformation.cpp
@@ -56,17 +56,22 @@ grid-values in units of mm/year in ENU-space.
 #include "proj.h"
 #include "proj_internal.h"
 #include <math.h>
+#include "grids.hpp"
 
 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;
+    ListOfHGrids hgrids{};
+    ListOfVGrids vgrids{};
 };
 } // anonymous namespace
 
@@ -85,13 +90,14 @@ 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);
+    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",
@@ -163,7 +169,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;
@@ -186,7 +192,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 +215,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 +231,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,11 +250,14 @@ 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);
 }
@@ -257,10 +266,9 @@ static PJ *destructor(PJ *P, int 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");
@@ -279,13 +287,13 @@ PJ *TRANSFORMATION(deformation,1) {
         return destructor(P, PJD_ERR_NO_ARGS );
     }
 
-    proj_hgrid_init(P, "xy_grids");
+    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);
     }
 
-    proj_vgrid_init(P, "z_grids");
+    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);
@@ -325,7 +333,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..e9983df6 100644
--- a/src/transformations/hgridshift.cpp
+++ b/src/transformations/hgridshift.cpp
@@ -6,24 +6,29 @@
 #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{};
 };
 } // 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->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 = proj_hgrid_apply(P, Q->grids, point.lp, PJ_FWD);
     }
 
     return point.xyz;
@@ -31,20 +36,21 @@ 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->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 = proj_hgrid_apply(P, 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 +68,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 +84,20 @@ 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);
+}
 
 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->fwd4d  = forward_4d;
     P->inv4d  = reverse_4d;
@@ -97,7 +111,7 @@ 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  */
@@ -121,11 +135,11 @@ PJ *TRANSFORMATION(hgridshift,0) {
         Q->t_epoch = pj_param (P->ctx, P->params, "dt_epoch").f;
 
 
-    proj_hgrid_init(P, "grids");
+    Q->grids = 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);
+        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..b964f45b 100644
--- a/src/transformations/vgridshift.cpp
+++ b/src/transformations/vgridshift.cpp
@@ -6,26 +6,30 @@
 #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{};
 };
 } // anonymous namespace
 
 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->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 += proj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier);
     }
 
     return point.xyz;
@@ -33,14 +37,14 @@ 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->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 -= proj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier);
     }
 
     return point.lpz;
@@ -48,7 +52,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 +70,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 +86,24 @@ 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);
+}
 
 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;
 
    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  */
@@ -121,12 +133,12 @@ PJ *TRANSFORMATION(vgridshift,0) {
     }
 
     /* Build gridlist. P->vgridlist_geoid can be empty if +grids only ask for optional grids. */
-    proj_vgrid_init(P, "grids");
+    Q->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);
+        return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
     }
 
     P->fwd4d = forward_4d;