horizontal grid shift: rework to no longer load whole grid into memory

16 files changed, 982 insertions, 1850 deletions

diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index ed6a62d6..c994df63 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -1553,12 +1553,9 @@ 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) {
-
+    {
         const auto gridSet = NS_PROJ::VerticalShiftGridSet::open(ctx, gridname);
         if( gridSet )
         {
@@ -1591,44 +1588,49 @@ PJ_GRID_INFO proj_grid_info(const char *gridname) {
                 grinfo.upperright.lam = extent.eastLon;
                 grinfo.upperright.phi = extent.northLat;
 
-                pj_gridinfo_free(ctx, gridinfo);
-
                 return grinfo;
             }
         }
-
-        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 */
-
-    /* name of grid */
-    strncpy (grinfo.gridname, gridname, sizeof(grinfo.gridname) - 1);
+    {
+        const auto gridSet = NS_PROJ::HorizontalShiftGridSet::open(ctx, gridname);
+        if( gridSet )
+        {
+            const auto& grids = gridSet->grids();
+            if( !grids.empty() )
+            {
+                const auto& grid = grids.front();
+                const auto& extent = grid->extentAndRes();
 
-    /* full path of grid */
-    pj_find_file(ctx, gridname, grinfo.filename, sizeof(grinfo.filename) - 1);
+                /* name of grid */
+                strncpy (grinfo.gridname, gridname, sizeof(grinfo.gridname) - 1);
 
-    /* grid format */
-    strncpy (grinfo.format, gridinfo->format, sizeof(grinfo.format) - 1);
+                /* full path of grid */
+                pj_find_file(ctx, gridname, grinfo.filename, sizeof(grinfo.filename) - 1);
 
-    /* grid size */
-    grinfo.n_lon = gridinfo->ct->lim.lam;
-    grinfo.n_lat = gridinfo->ct->lim.phi;
+                /* grid format */
+                strncpy (grinfo.format, gridSet->format().c_str(), sizeof(grinfo.format) - 1);
 
-    /* cell size */
-    grinfo.cs_lon = gridinfo->ct->del.lam;
-    grinfo.cs_lat = gridinfo->ct->del.phi;
+                /* grid size */
+                grinfo.n_lon = grid->width();
+                grinfo.n_lat = grid->height();
 
-    /* bounds of grid */
-    grinfo.lowerleft  = gridinfo->ct->ll;
-    grinfo.upperright.lam = grinfo.lowerleft.lam + (grinfo.n_lon-1)*grinfo.cs_lon;
-    grinfo.upperright.phi = grinfo.lowerleft.phi + (grinfo.n_lat-1)*grinfo.cs_lat;
+                /* cell size */
+                grinfo.cs_lon = extent.resLon;
+                grinfo.cs_lat = extent.resLat;
 
-    pj_gridinfo_free(ctx, gridinfo);
+                /* bounds of grid */
+                grinfo.lowerleft.lam  = extent.westLon;
+                grinfo.lowerleft.phi  = extent.southLat;
+                grinfo.upperright.lam = extent.eastLon;
+                grinfo.upperright.phi = extent.northLat;
 
+                return grinfo;
+            }
+        }
+    }
+    strcpy(grinfo.format, "missing");
     return grinfo;
 }
 
diff --git a/src/Makefile.am b/src/Makefile.am
index 3d0868a7..a12de4e1 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -195,9 +195,8 @@ libproj_la_SOURCES = \
 	release.cpp gauss.cpp \
 	fileapi.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 \
 	\
diff --git a/src/apply_gridshift.cpp b/src/apply_gridshift.cpp
index 2f46ca11..5a5404da 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,231 +38,62 @@
 
 #include "proj.h"
 #include "proj_internal.h"
+#include "proj/internal/internal.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.                                    */
-/************************************************************************/
+using namespace NS_PROJ;
 
-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 std::vector<std::unique_ptr<HorizontalShiftGridSet>>& 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->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 std::vector<std::unique_ptr<NS_PROJ::HorizontalShiftGridSet>>
+                                        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 )
+    std::vector<std::unique_ptr<NS_PROJ::HorizontalShiftGridSet>> 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) {
+int proj_hgrid_init(PJ* P, const char *gridkey) {
 /**********************************************
 
   Initizalize and populate list of horizontal
@@ -275,29 +108,158 @@ 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 0;
 
-    if (P->gridlist == nullptr) {
-        P->gridlist = pj_gridlist_from_nadgrids(
-            P->ctx,
-            pj_param(P->ctx, P->params, sgrids).s,
-            &(P->gridlist_count)
-        );
+    P->hgrids = getListOfGridSets(P->ctx, grids);
+    return static_cast<int>(P->hgrids.size());
+}
 
-        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;
 }
 
 /********************************************/
@@ -306,22 +268,22 @@ 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 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(P->hgrids, 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);
@@ -330,33 +292,89 @@ 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 std::vector<std::unique_ptr<HorizontalShiftGridSet>>& 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, PJ_LP lp, PJ_DIRECTION direction) {
+    return proj_hgrid_apply_internal(P->ctx, lp, direction, P->hgrids);
+}
+
+/************************************************************************/
+/*                         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 = 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/gridinfo.cpp b/src/gridinfo.cpp
deleted file mode 100644
index 01d392b3..00000000
--- a/src/gridinfo.cpp
+++ /dev/null
@@ -1,834 +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;
-    }
-
-    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()                          */
-/*                                                                      */
-/*      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( 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
index 1a18f975..54dae555 100644
--- a/src/grids.cpp
+++ b/src/grids.cpp
@@ -302,4 +302,582 @@ const VerticalShiftGrid *VerticalShiftGridSet::gridAt(double lon,
     return nullptr;
 }
 
+// ---------------------------------------------------------------------------
+
+HorizontalShiftGrid::HorizontalShiftGrid(int widthIn, int heightIn,
+                                         const ExtentAndRes &extentIn)
+    : m_width(widthIn), m_height(heightIn), m_extent(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
index 0e3ccffb..0c4af09f 100644
--- a/src/grids.hpp
+++ b/src/grids.hpp
@@ -98,6 +98,57 @@ class VerticalShiftGridSet {
     const VerticalShiftGrid *gridAt(double lon, double lat) const;
 };
 
+// ---------------------------------------------------------------------------
+
+class HorizontalShiftGrid {
+  protected:
+    int m_width;
+    int m_height;
+    ExtentAndRes m_extent;
+    std::vector<std::unique_ptr<HorizontalShiftGrid>> m_children{};
+
+  public:
+    HorizontalShiftGrid(int widthIn, int heightIn,
+                        const ExtentAndRes &extentIn);
+    virtual ~HorizontalShiftGrid();
+
+    int width() const { return m_width; }
+    int height() const { return m_height; }
+    const ExtentAndRes &extentAndRes() const { return m_extent; }
+
+    const HorizontalShiftGrid *gridAt(double lon, double lat) const;
+
+    virtual bool isNullGrid() const { return false; }
+
+    // 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;
+};
+
 NS_PROJ_END
 
 #endif // GRIDS_HPP_INCLUDED
diff --git a/src/init.cpp b/src/init.cpp
index 69b4ec5f..19fcf47b 100644
--- a/src/init.cpp
+++ b/src/init.cpp
@@ -646,9 +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;
-
     /* 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 1cc3b90c..f6112aef 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -237,8 +237,6 @@ set(SRC_LIBPROJ_CORE
   geocent.cpp
   geocent.h
   geodesic.c
-  gridinfo.cpp
-  gridlist.cpp
   init.cpp
   initcache.cpp
   internal.cpp
@@ -250,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
diff --git a/src/malloc.cpp b/src/malloc.cpp
index 050721e6..b5fe76fc 100644
--- a/src/malloc.cpp
+++ b/src/malloc.cpp
@@ -225,9 +225,6 @@ 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 );
-
     /* We used to call pj_dalloc( P->catalog ), but this will leak */
     /* memory. The safe way to clear catalog and grid is to call */
     /* pj_gc_unloadall(pj_get_default_ctx()); and pj_deallocate_grids(); */
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 1e43e066..b7a84f87 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -479,8 +479,7 @@ 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;
+    std::vector<std::unique_ptr<NS_PROJ::HorizontalShiftGridSet>> hgrids{};
 
     int     has_geoid_vgrids = 0;      /* TODO: Description needed */
     std::vector<std::unique_ptr<NS_PROJ::VerticalShiftGridSet>> vgrids{};
@@ -752,35 +751,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;
-
 /* procedure prototypes */
 double PROJ_DLL dmstor(const char *, char **);
 double dmstor_ctx(projCtx_t *ctx, const char *, char **);
@@ -827,32 +797,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_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 * );
-
 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 863da605..f6002b70 100644
--- a/src/transform.cpp
+++ b/src/transform.cpp
@@ -867,6 +867,58 @@ 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.empty() )
+    {
+        proj_hgrid_init(defn, "nadgrids");
+    }
+
+    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, 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()                         */
 /*                                                                      */
@@ -1107,3 +1159,8 @@ static int adjust_axis( projCtx ctx,
 
     return 0;
 }
+// ---------------------------------------------------------------------------
+
+void pj_deallocate_grids()
+{
+}
diff --git a/src/transformations/hgridshift.cpp b/src/transformations/hgridshift.cpp
index 90633939..465bfe98 100644
--- a/src/transformations/hgridshift.cpp
+++ b/src/transformations/hgridshift.cpp
@@ -20,7 +20,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
     PJ_COORD point = {{0,0,0,0}};
     point.lpz = lpz;
 
-    if (P->gridlist != nullptr) {
+    if (!P->hgrids.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);
@@ -34,7 +34,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
     PJ_COORD point = {{0,0,0,0}};
     point.xyz = xyz;
 
-    if (P->gridlist != nullptr) {
+    if (!P->hgrids.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);