Grid refactoring: address review comments of https://github.com/OSGeo/PROJ/pull/1777

- Move content of legacy apply_gridshift.cpp and apply_vgridshift.cpp in grids.cpp
- Rename nad_ functions to pj_hgrid_
- Rename internal proj_hgrid_/proj_vgrid_ functions to pj_

11 files changed, 529 insertions, 649 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index d29eb976..39667509 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -196,9 +196,9 @@ libproj_la_SOURCES = \
 	release.cpp gauss.cpp \
 	fileapi.cpp \
 	\
-	apply_gridshift.cpp datums.cpp datum_set.cpp transform.cpp \
+	datums.cpp datum_set.cpp transform.cpp \
 	geocent.cpp geocent.h utils.cpp \
-	jniproj.cpp mutex.cpp initcache.cpp apply_vgridshift.cpp geodesic.c \
+	jniproj.cpp mutex.cpp initcache.cpp geodesic.c \
 	strtod.cpp \
 	\
 	4D_api.cpp pipeline.cpp \
diff --git a/src/apply_gridshift.cpp b/src/apply_gridshift.cpp
deleted file mode 100644
index 4ef86fc0..00000000
--- a/src/apply_gridshift.cpp
+++ /dev/null
@@ -1,397 +0,0 @@
-/******************************************************************************
- * Project:  PROJ.4
- * Purpose:  Apply datum shifts based on grid shift files (normally NAD27 to
- *           NAD83 or the reverse).  This module is responsible for keeping
- *           a list of loaded grids, and calling with each one that is
- *           allowed for a given datum (expressed as the nadgrids= parameter).
- * 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.
- *****************************************************************************/
-
-#ifndef FROM_PROJ_CPP
-#define FROM_PROJ_CPP
-#endif
-
-#include <math.h>
-#include <stdio.h>
-#include <string.h>
-
-#include <limits>
-
-#include "proj.h"
-#include "proj_internal.h"
-#include "proj/internal/internal.hpp"
-#include "grids.hpp"
-
-NS_PROJ_START
-
-// ---------------------------------------------------------------------------
-
-static const HorizontalShiftGrid* findGrid(const ListOfHGrids& grids,
-                                           PJ_LP input)
-{
-    for( const auto& gridset: grids )
-    {
-        auto grid = gridset->gridAt(input.lam, input.phi);
-        if( grid )
-            return grid;
-    }
-    return nullptr;
-}
-
-// ---------------------------------------------------------------------------
-
-static ListOfHGrids getListOfGridSets(PJ_CONTEXT* ctx, const char* grids)
-{
-    ListOfHGrids list;
-    auto listOfGrids = internal::split(std::string(grids), ',');
-    for( const auto& grid: listOfGrids )
-    {
-        const char* gridname = grid.c_str();
-        bool canFail = false;
-        if( gridname[0] == '@' )
-        {
-            canFail = true;
-            gridname ++;
-        }
-        auto gridSet = HorizontalShiftGridSet::open(ctx, gridname);
-        if( !gridSet )
-        {
-            if( !canFail )
-            {
-                pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-                return {};
-            }
-        }
-        else
-        {
-            list.emplace_back(std::move(gridSet));
-        }
-    }
-    return list;
-}
-
-
-/**********************************************/
-ListOfHGrids proj_hgrid_init(PJ* P, const char *gridkey) {
-/**********************************************
-
-  Initizalize and populate list of horizontal
-  grids.
-
-    Takes a PJ-object and the plus-parameter
-    name that is used in the proj-string to
-    specify the grids to load, e.g. "+grids".
-    The + should be left out here.
-
-    Returns the number of loaded grids.
-
-***********************************************/
-
-    std::string key("s");
-    key += gridkey;
-    const char* grids = pj_param(P->ctx, P->params, key.c_str()).s;
-    if( grids == nullptr )
-        return {};
-
-    return getListOfGridSets(P->ctx, grids);
-}
-
-typedef struct { pj_int32 lam, phi; } ILP;
-
-static PJ_LP nad_intr(PJ_LP t, const HorizontalShiftGrid* grid, bool compensateNTConvention) {
-	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, compensateNTConvention, f00Lon, f00Lat) ||
-            !grid->valueAt(indx.lam + 1, indx.phi, compensateNTConvention, f10Lon, f10Lat) ||
-            !grid->valueAt(indx.lam, indx.phi + 1, compensateNTConvention, f01Lon, f01Lat) ||
-            !grid->valueAt(indx.lam + 1, indx.phi + 1, compensateNTConvention, 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(projCtx ctx, PJ_LP in, int inverse, const HorizontalShiftGrid* grid,
-              const ListOfHGrids& grids) {
-    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, true);
-    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, grid, true);
-
-        /* We can possibly go outside of the initial guessed grid, so try */
-        /* to fetch a new grid into which iterate... */
-        if (del.lam == HUGE_VAL)
-        {
-            PJ_LP lp;
-            lp.lam = t.lam + extent->westLon;
-            lp.phi = t.phi + extent->southLat;
-            auto newGrid = findGrid(grids, lp);
-            if( newGrid == nullptr || newGrid == grid || newGrid->isNullGrid() )
-                break;
-            pj_log(ctx, PJ_LOG_DEBUG_MINOR, "Switching from grid %s to grid %s",
-                   grid->name().c_str(),
-                   newGrid->name().c_str());
-            grid = newGrid;
-            extent = &(grid->extentAndRes());
-            t.lam = lp.lam - extent->westLon;
-            t.phi = lp.phi - extent->southLat;
-            tb = in;
-            tb.lam -= extent->westLon;
-            tb.phi -= extent->southLat;
-            tb.lam = adjlon (tb.lam - M_PI) + M_PI;
-            dif.lam = std::numeric_limits<double>::max();
-            dif.phi = std::numeric_limits<double>::max();
-            continue;
-        }
-
-        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 + extent->westLon);
-    in.phi = t.phi + extent->southLat;
-    return in;
-}
-
-/********************************************/
-/*           proj_hgrid_value()             */
-/*                                          */
-/*    Return coordinate offset in grid      */
-/********************************************/
-PJ_LP proj_hgrid_value(PJ *P, const ListOfHGrids& grids, PJ_LP lp) {
-    PJ_LP out = proj_coord_error().lp;
-
-    const auto grid = findGrid(grids, lp);
-    if( !grid ) {
-        pj_ctx_set_errno( P->ctx, PJD_ERR_GRID_AREA );
-        return out;
-    }
-
-    /* normalize input to ll origin */
-    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, grid, false);
-
-    if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) {
-        pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA);
-    }
-
-    return out;
-}
-
-static
-PJ_LP proj_hgrid_apply_internal(PJ_CONTEXT* ctx,
-                       PJ_LP lp,
-                       PJ_DIRECTION direction,
-                       const ListOfHGrids& grids)
-{
-    PJ_LP out;
-
-    out.lam = HUGE_VAL;
-    out.phi = HUGE_VAL;
-
-    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;
-    }
-
-    int inverse = direction == PJ_FWD ? 0 : 1;
-    out = nad_cvt(ctx, lp, inverse, grid, grids);
-
-    if (out.lam == HUGE_VAL || out.phi == HUGE_VAL)
-        pj_ctx_set_errno(ctx, PJD_ERR_GRID_AREA);
-
-    return out;
-}
-
-PJ_LP proj_hgrid_apply(PJ *P, const ListOfHGrids& grids, PJ_LP lp, PJ_DIRECTION direction) {
-    return proj_hgrid_apply_internal(P->ctx, lp, direction, grids);
-}
-
-
-NS_PROJ_END
-
-/************************************************************************/
-/*                         pj_apply_gridshift()                         */
-/*                                                                      */
-/*      This is the externally callable interface - part of the         */
-/*      public API - though it is not used internally any more and I    */
-/*      doubt it is used by any other applications.  But we preserve    */
-/*      it to honour our public api.                                    */
-/************************************************************************/
-
-int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse,
-                        long point_count, int point_offset,
-                        double *x, double *y, double * /*z */ )
-
-{
-    auto hgrids = NS_PROJ::getListOfGridSets(ctx, nadgrids);
-    if( hgrids.empty() )
-    {
-        pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-        return 1;
-    }
-
-
-    for( long i = 0; i < point_count; i++ )
-    {
-        PJ_LP   input;
-
-        long io = i * point_offset;
-        input.phi = y[io];
-        input.lam = x[io];
-
-        auto output = proj_hgrid_apply_internal(ctx, input, inverse ? PJ_INV : PJ_FWD, hgrids);
-
-        if ( output.lam != HUGE_VAL )
-        {
-            y[io] = output.phi;
-            x[io] = output.lam;
-        }
-        else
-        {
-            if( ctx->debug_level >= PJ_LOG_DEBUG_MAJOR )
-            {
-                pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
-                    "pj_apply_gridshift(): failed to find a grid shift table for\n"
-                    "                      location (%.7fdW,%.7fdN)",
-                    x[io] * RAD_TO_DEG,
-                    y[io] * RAD_TO_DEG );
-            }
-        }
-    }
-
-    return 0;
-}
diff --git a/src/apply_vgridshift.cpp b/src/apply_vgridshift.cpp
deleted file mode 100644
index b0136e5c..00000000
--- a/src/apply_vgridshift.cpp
+++ /dev/null
@@ -1,217 +0,0 @@
-/******************************************************************************
- * Project:  PROJ.4
- * Purpose:  Apply vertical datum shifts based on grid shift files, normally
- *           geoid grids mapping WGS84 to NAVD88 or something similar.
- * Author:   Frank Warmerdam, warmerdam@pobox.com
- *
- ******************************************************************************
- * Copyright (c) 2010, 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.
- *****************************************************************************/
-
-#ifndef FROM_PROJ_CPP
-#define FROM_PROJ_CPP
-#endif
-
-#include <assert.h>
-#include <stdio.h>
-#include <string.h>
-
-#include <math.h>
-#include "proj_internal.h"
-#include "proj/internal/internal.hpp"
-#include "grids.hpp"
-
-NS_PROJ_START
-
-static double read_vgrid_value(const ListOfVGrids& grids, PJ_LP input, double vmultiplier) {
-
-    /* do not deal with NaN coordinates */
-    /* cppcheck-suppress duplicateExpression */
-    if( isnan(input.phi) || isnan(input.lam) )
-    {
-        return HUGE_VAL;
-    }
-
-    const VerticalShiftGrid* grid = nullptr;
-    for( const auto& gridset: grids )
-    {
-        grid = gridset->gridAt(input.lam, input.phi);
-        if( grid )
-            break;
-    }
-    if( !grid )
-    {
-        return HUGE_VAL;
-    }
-
-    const auto& extent = grid->extentAndRes();
-
-    /* Interpolation a location within the grid */
-    double grid_x = (input.lam - extent.westLon) / extent.resLon;
-    if( extent.fullWorldLongitude() ) {
-        // The first fmod goes to ]-lim, lim[ range
-        // So we add lim again to be in ]0, 2*lim[ and fmod again
-        grid_x = fmod(
-            fmod(grid_x + grid->width(), grid->width()) + grid->width(),
-            grid->width());
-    }
-    double grid_y = (input.phi - extent.southLat) / extent.resLat;
-    int grid_ix = static_cast<int>(lround(floor(grid_x)));
-    assert(grid_ix >= 0 && grid_ix < grid->width());
-    int grid_iy = static_cast<int>(lround(floor(grid_y)));
-    assert(grid_iy >= 0 && grid_iy < grid->height());
-    grid_x -= grid_ix;
-    grid_y -= grid_iy;
-
-    int grid_ix2 = grid_ix + 1;
-    if( grid_ix2 >= grid->width() ) {
-        if( extent.fullWorldLongitude() ) {
-            grid_ix2 = 0;
-        } else {
-            grid_ix2 = grid->width() - 1;
-        }
-    }
-    int grid_iy2 = grid_iy + 1;
-    if( grid_iy2 >= grid->height() )
-        grid_iy2 = grid->height() - 1;
-
-    float value_a = 0;
-    float value_b = 0;
-    float value_c = 0;
-    float value_d = 0;
-    if( !grid->valueAt(grid_ix, grid_iy, value_a) ||
-        !grid->valueAt(grid_ix2, grid_iy, value_b) ||
-        !grid->valueAt(grid_ix, grid_iy2, value_c) ||
-        !grid->valueAt(grid_ix2, grid_iy2, value_d) )
-    {
-        return HUGE_VAL;
-    }
-
-    double total_weight = 0.0;
-    int n_weights = 0;
-    double value = 0.0f;
-
-    if( !grid->isNodata(value_a, vmultiplier) )
-    {
-        double weight = (1.0-grid_x) * (1.0-grid_y);
-        value += value_a * weight;
-        total_weight += weight;
-        n_weights ++;
-    }
-    if( !grid->isNodata(value_b, vmultiplier) )
-    {
-        double weight = (grid_x) * (1.0-grid_y);
-        value += value_b * weight;
-        total_weight += weight;
-        n_weights ++;
-    }
-    if( !grid->isNodata(value_c, vmultiplier) )
-    {
-        double weight = (1.0-grid_x) * (grid_y);
-        value += value_c * weight;
-        total_weight += weight;
-        n_weights ++;
-    }
-    if( !grid->isNodata(value_d, vmultiplier) )
-    {
-        double weight = (grid_x) * (grid_y);
-        value += value_d * weight;
-        total_weight += weight;
-        n_weights ++;
-    }
-    if( n_weights == 0 )
-        value = HUGE_VAL;
-    else if( n_weights != 4 )
-        value /= total_weight;
-
-    return value * vmultiplier;
-}
-
-/**********************************************/
-ListOfVGrids proj_vgrid_init(PJ* P, const char *gridkey) {
-/**********************************************
-
-  Initizalize and populate gridlist.
-
-    Takes a PJ-object and the plus-parameter
-    name that is used in the proj-string to
-    specify the grids to load, e.g. "+grids".
-    The + should be left out here.
-
-    Returns the number of loaded grids.
-
-***********************************************/
-
-    std::string key("s");
-    key += gridkey;
-    const char* gridnames = pj_param(P->ctx, P->params, key.c_str()).s;
-    if( gridnames == nullptr )
-        return {};
-
-    auto listOfGridNames = internal::split(std::string(gridnames), ',');
-    ListOfVGrids grids;
-    for( const auto& gridnameStr: listOfGridNames )
-    {
-        const char* gridname = gridnameStr.c_str();
-        bool canFail = false;
-        if( gridname[0] == '@' )
-        {
-            canFail = true;
-            gridname ++;
-        }
-        auto gridSet = VerticalShiftGridSet::open(P->ctx, gridname);
-        if( !gridSet )
-        {
-            if( !canFail )
-            {
-                pj_ctx_set_errno( P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
-                return {};
-            }
-        }
-        else
-        {
-            grids.emplace_back(std::move(gridSet));
-        }
-    }
-
-    return grids;
-}
-
-/***********************************************/
-double proj_vgrid_value(PJ *P, const ListOfVGrids& grids, PJ_LP lp, double vmultiplier){
-/***********************************************
-
-  Read grid value at position lp in grids loaded
-  with proj_grid_init.
-
-  Returns the grid value of the given coordinate.
-
-************************************************/
-
-    double value;
-
-    value = read_vgrid_value(grids, lp, vmultiplier);
-    proj_log_trace(P, "proj_vgrid_value: (%f, %f) = %f", lp.lam*RAD_TO_DEG, lp.phi*RAD_TO_DEG, value);
-
-    return value;
-}
-
-NS_PROJ_END
diff --git a/src/grids.cpp b/src/grids.cpp
index 5a99106b..66854188 100644
--- a/src/grids.cpp
+++ b/src/grids.cpp
@@ -2614,7 +2614,7 @@ void GenericShiftGridSet::reassign_context(PJ_CONTEXT *ctx) {
 
 // ---------------------------------------------------------------------------
 
-ListOfGenericGrids proj_generic_grid_init(PJ *P, const char *gridkey) {
+ListOfGenericGrids pj_generic_grid_init(PJ *P, const char *gridkey) {
     std::string key("s");
     key += gridkey;
     const char *gridnames = pj_param(P->ctx, P->params, key.c_str()).s;
@@ -2644,4 +2644,500 @@ ListOfGenericGrids proj_generic_grid_init(PJ *P, const char *gridkey) {
     return grids;
 }
 
+// ---------------------------------------------------------------------------
+
+static const HorizontalShiftGrid *findGrid(const ListOfHGrids &grids,
+                                           PJ_LP input) {
+    for (const auto &gridset : grids) {
+        auto grid = gridset->gridAt(input.lam, input.phi);
+        if (grid)
+            return grid;
+    }
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+static ListOfHGrids getListOfGridSets(PJ_CONTEXT *ctx, const char *grids) {
+    ListOfHGrids list;
+    auto listOfGrids = internal::split(std::string(grids), ',');
+    for (const auto &grid : listOfGrids) {
+        const char *gridname = grid.c_str();
+        bool canFail = false;
+        if (gridname[0] == '@') {
+            canFail = true;
+            gridname++;
+        }
+        auto gridSet = HorizontalShiftGridSet::open(ctx, gridname);
+        if (!gridSet) {
+            if (!canFail) {
+                pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+                return {};
+            }
+        } else {
+            list.emplace_back(std::move(gridSet));
+        }
+    }
+    return list;
+}
+
+/**********************************************/
+ListOfHGrids pj_hgrid_init(PJ *P, const char *gridkey) {
+    /**********************************************
+
+      Initizalize and populate list of horizontal
+      grids.
+
+        Takes a PJ-object and the plus-parameter
+        name that is used in the proj-string to
+        specify the grids to load, e.g. "+grids".
+        The + should be left out here.
+
+        Returns the number of loaded grids.
+
+    ***********************************************/
+
+    std::string key("s");
+    key += gridkey;
+    const char *grids = pj_param(P->ctx, P->params, key.c_str()).s;
+    if (grids == nullptr)
+        return {};
+
+    return getListOfGridSets(P->ctx, grids);
+}
+
+// ---------------------------------------------------------------------------
+
+typedef struct { pj_int32 lam, phi; } ILP;
+
+static PJ_LP pj_hgrid_interpolate(PJ_LP t, const HorizontalShiftGrid *grid,
+                                  bool compensateNTConvention) {
+    PJ_LP val, frct;
+    ILP indx;
+    int in;
+
+    const auto &extent = grid->extentAndRes();
+    t.lam /= extent.resLon;
+    indx.lam = std::isnan(t.lam) ? 0 : (pj_int32)lround(floor(t.lam));
+    t.phi /= extent.resLat;
+    indx.phi = std::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, compensateNTConvention, f00Lon,
+                       f00Lat) ||
+        !grid->valueAt(indx.lam + 1, indx.phi, compensateNTConvention, f10Lon,
+                       f10Lat) ||
+        !grid->valueAt(indx.lam, indx.phi + 1, compensateNTConvention, f01Lon,
+                       f01Lat) ||
+        !grid->valueAt(indx.lam + 1, indx.phi + 1, compensateNTConvention,
+                       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 pj_hgrid_apply_internal(projCtx ctx, PJ_LP in,
+                                     PJ_DIRECTION direction,
+                                     const HorizontalShiftGrid *grid,
+                                     const ListOfHGrids &grids) {
+    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 = pj_hgrid_interpolate(tb, grid, true);
+    if (t.lam == HUGE_VAL)
+        return t;
+
+    if (direction == PJ_FWD) {
+        in.lam += t.lam;
+        in.phi += t.phi;
+        return in;
+    }
+
+    t.lam = tb.lam - t.lam;
+    t.phi = tb.phi - t.phi;
+
+    do {
+        del = pj_hgrid_interpolate(t, grid, true);
+
+        /* We can possibly go outside of the initial guessed grid, so try */
+        /* to fetch a new grid into which iterate... */
+        if (del.lam == HUGE_VAL) {
+            PJ_LP lp;
+            lp.lam = t.lam + extent->westLon;
+            lp.phi = t.phi + extent->southLat;
+            auto newGrid = findGrid(grids, lp);
+            if (newGrid == nullptr || newGrid == grid || newGrid->isNullGrid())
+                break;
+            pj_log(ctx, PJ_LOG_DEBUG_MINOR, "Switching from grid %s to grid %s",
+                   grid->name().c_str(), newGrid->name().c_str());
+            grid = newGrid;
+            extent = &(grid->extentAndRes());
+            t.lam = lp.lam - extent->westLon;
+            t.phi = lp.phi - extent->southLat;
+            tb = in;
+            tb.lam -= extent->westLon;
+            tb.phi -= extent->southLat;
+            tb.lam = adjlon(tb.lam - M_PI) + M_PI;
+            dif.lam = std::numeric_limits<double>::max();
+            dif.phi = std::numeric_limits<double>::max();
+            continue;
+        }
+
+        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 + extent->westLon);
+    in.phi = t.phi + extent->southLat;
+    return in;
+}
+
+// ---------------------------------------------------------------------------
+
+PJ_LP pj_hgrid_apply(PJ_CONTEXT *ctx, const ListOfHGrids &grids, PJ_LP lp,
+                     PJ_DIRECTION direction) {
+    PJ_LP out;
+
+    out.lam = HUGE_VAL;
+    out.phi = HUGE_VAL;
+
+    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;
+    }
+
+    out = pj_hgrid_apply_internal(ctx, lp, direction, grid, grids);
+
+    if (out.lam == HUGE_VAL || out.phi == HUGE_VAL)
+        pj_ctx_set_errno(ctx, PJD_ERR_GRID_AREA);
+
+    return out;
+}
+
+/********************************************/
+/*           proj_hgrid_value()             */
+/*                                          */
+/*    Return coordinate offset in grid      */
+/********************************************/
+PJ_LP pj_hgrid_value(PJ *P, const ListOfHGrids &grids, PJ_LP lp) {
+    PJ_LP out = proj_coord_error().lp;
+
+    const auto grid = findGrid(grids, lp);
+    if (!grid) {
+        pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA);
+        return out;
+    }
+
+    /* normalize input to ll origin */
+    const auto &extent = grid->extentAndRes();
+    lp.lam -= extent.westLon;
+    lp.phi -= extent.southLat;
+
+    lp.lam = adjlon(lp.lam - M_PI) + M_PI;
+
+    out = pj_hgrid_interpolate(lp, grid, false);
+
+    if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) {
+        pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA);
+    }
+
+    return out;
+}
+
+// ---------------------------------------------------------------------------
+
+static double read_vgrid_value(const ListOfVGrids &grids, PJ_LP input,
+                               double vmultiplier) {
+
+    /* do not deal with NaN coordinates */
+    /* cppcheck-suppress duplicateExpression */
+    if (std::isnan(input.phi) || std::isnan(input.lam)) {
+        return HUGE_VAL;
+    }
+
+    const VerticalShiftGrid *grid = nullptr;
+    for (const auto &gridset : grids) {
+        grid = gridset->gridAt(input.lam, input.phi);
+        if (grid)
+            break;
+    }
+    if (!grid) {
+        return HUGE_VAL;
+    }
+
+    const auto &extent = grid->extentAndRes();
+
+    /* Interpolation a location within the grid */
+    double grid_x = (input.lam - extent.westLon) / extent.resLon;
+    if (extent.fullWorldLongitude()) {
+        // The first fmod goes to ]-lim, lim[ range
+        // So we add lim again to be in ]0, 2*lim[ and fmod again
+        grid_x =
+            fmod(fmod(grid_x + grid->width(), grid->width()) + grid->width(),
+                 grid->width());
+    }
+    double grid_y = (input.phi - extent.southLat) / extent.resLat;
+    int grid_ix = static_cast<int>(lround(floor(grid_x)));
+    assert(grid_ix >= 0 && grid_ix < grid->width());
+    int grid_iy = static_cast<int>(lround(floor(grid_y)));
+    assert(grid_iy >= 0 && grid_iy < grid->height());
+    grid_x -= grid_ix;
+    grid_y -= grid_iy;
+
+    int grid_ix2 = grid_ix + 1;
+    if (grid_ix2 >= grid->width()) {
+        if (extent.fullWorldLongitude()) {
+            grid_ix2 = 0;
+        } else {
+            grid_ix2 = grid->width() - 1;
+        }
+    }
+    int grid_iy2 = grid_iy + 1;
+    if (grid_iy2 >= grid->height())
+        grid_iy2 = grid->height() - 1;
+
+    float value_a = 0;
+    float value_b = 0;
+    float value_c = 0;
+    float value_d = 0;
+    if (!grid->valueAt(grid_ix, grid_iy, value_a) ||
+        !grid->valueAt(grid_ix2, grid_iy, value_b) ||
+        !grid->valueAt(grid_ix, grid_iy2, value_c) ||
+        !grid->valueAt(grid_ix2, grid_iy2, value_d)) {
+        return HUGE_VAL;
+    }
+
+    double total_weight = 0.0;
+    int n_weights = 0;
+    double value = 0.0f;
+
+    if (!grid->isNodata(value_a, vmultiplier)) {
+        double weight = (1.0 - grid_x) * (1.0 - grid_y);
+        value += value_a * weight;
+        total_weight += weight;
+        n_weights++;
+    }
+    if (!grid->isNodata(value_b, vmultiplier)) {
+        double weight = (grid_x) * (1.0 - grid_y);
+        value += value_b * weight;
+        total_weight += weight;
+        n_weights++;
+    }
+    if (!grid->isNodata(value_c, vmultiplier)) {
+        double weight = (1.0 - grid_x) * (grid_y);
+        value += value_c * weight;
+        total_weight += weight;
+        n_weights++;
+    }
+    if (!grid->isNodata(value_d, vmultiplier)) {
+        double weight = (grid_x) * (grid_y);
+        value += value_d * weight;
+        total_weight += weight;
+        n_weights++;
+    }
+    if (n_weights == 0)
+        value = HUGE_VAL;
+    else if (n_weights != 4)
+        value /= total_weight;
+
+    return value * vmultiplier;
+}
+
+/**********************************************/
+ListOfVGrids pj_vgrid_init(PJ *P, const char *gridkey) {
+    /**********************************************
+
+      Initizalize and populate gridlist.
+
+        Takes a PJ-object and the plus-parameter
+        name that is used in the proj-string to
+        specify the grids to load, e.g. "+grids".
+        The + should be left out here.
+
+        Returns the number of loaded grids.
+
+    ***********************************************/
+
+    std::string key("s");
+    key += gridkey;
+    const char *gridnames = pj_param(P->ctx, P->params, key.c_str()).s;
+    if (gridnames == nullptr)
+        return {};
+
+    auto listOfGridNames = internal::split(std::string(gridnames), ',');
+    ListOfVGrids grids;
+    for (const auto &gridnameStr : listOfGridNames) {
+        const char *gridname = gridnameStr.c_str();
+        bool canFail = false;
+        if (gridname[0] == '@') {
+            canFail = true;
+            gridname++;
+        }
+        auto gridSet = VerticalShiftGridSet::open(P->ctx, gridname);
+        if (!gridSet) {
+            if (!canFail) {
+                pj_ctx_set_errno(P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+                return {};
+            }
+        } else {
+            grids.emplace_back(std::move(gridSet));
+        }
+    }
+
+    return grids;
+}
+
+/***********************************************/
+double pj_vgrid_value(PJ *P, const ListOfVGrids &grids, PJ_LP lp,
+                      double vmultiplier) {
+    /***********************************************
+
+      Read grid value at position lp in grids loaded
+      with proj_grid_init.
+
+      Returns the grid value of the given coordinate.
+
+    ************************************************/
+
+    double value;
+
+    value = read_vgrid_value(grids, lp, vmultiplier);
+    proj_log_trace(P, "proj_vgrid_value: (%f, %f) = %f", lp.lam * RAD_TO_DEG,
+                   lp.phi * RAD_TO_DEG, value);
+
+    return value;
+}
+
 NS_PROJ_END
+
+/************************************************************************/
+/*                         pj_apply_gridshift()                         */
+/*                                                                      */
+/*      This is the externally callable interface - part of the         */
+/*      public API - though it is not used internally any more and I    */
+/*      doubt it is used by any other applications.  But we preserve    */
+/*      it to honour our public api.                                    */
+/************************************************************************/
+
+int pj_apply_gridshift(projCtx ctx, const char *nadgrids, int inverse,
+                       long point_count, int point_offset, double *x, double *y,
+                       double * /*z */)
+
+{
+    auto hgrids = NS_PROJ::getListOfGridSets(ctx, nadgrids);
+    if (hgrids.empty()) {
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return 1;
+    }
+
+    for (long i = 0; i < point_count; i++) {
+        PJ_LP input;
+
+        long io = i * point_offset;
+        input.phi = y[io];
+        input.lam = x[io];
+
+        auto output =
+            pj_hgrid_apply(ctx, hgrids, input, inverse ? PJ_INV : PJ_FWD);
+
+        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/grids.hpp b/src/grids.hpp
index aa852ef6..fbe4e7f8 100644
--- a/src/grids.hpp
+++ b/src/grids.hpp
@@ -225,15 +225,15 @@ typedef std::vector<std::unique_ptr<HorizontalShiftGridSet>> ListOfHGrids;
 typedef std::vector<std::unique_ptr<VerticalShiftGridSet>> ListOfVGrids;
 typedef std::vector<std::unique_ptr<GenericShiftGridSet>> ListOfGenericGrids;
 
-ListOfVGrids proj_vgrid_init(PJ *P, const char *grids);
-ListOfHGrids proj_hgrid_init(PJ *P, const char *grids);
-ListOfGenericGrids proj_generic_grid_init(PJ *P, const char *grids);
-
-double proj_vgrid_value(PJ *P, const ListOfVGrids &, PJ_LP lp,
-                        double vmultiplier);
-PJ_LP proj_hgrid_value(PJ *P, const ListOfHGrids &, PJ_LP lp);
-PJ_LP proj_hgrid_apply(PJ *P, const ListOfHGrids &, PJ_LP lp,
-                       PJ_DIRECTION direction);
+ListOfVGrids pj_vgrid_init(PJ *P, const char *grids);
+ListOfHGrids pj_hgrid_init(PJ *P, const char *grids);
+ListOfGenericGrids pj_generic_grid_init(PJ *P, const char *grids);
+
+PJ_LP pj_hgrid_value(PJ *P, const ListOfHGrids &grids, PJ_LP lp);
+double pj_vgrid_value(PJ *P, const ListOfVGrids &, PJ_LP lp,
+                      double vmultiplier);
+PJ_LP pj_hgrid_apply(PJ_CONTEXT *ctx, const ListOfHGrids &grids, PJ_LP lp,
+                     PJ_DIRECTION direction);
 
 NS_PROJ_END
 
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 18c91021..fdb59434 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -220,8 +220,6 @@ set(SRC_LIBPROJ_CORE
   4D_api.cpp
   aasincos.cpp
   adjlon.cpp
-  apply_gridshift.cpp
-  apply_vgridshift.cpp
   auth.cpp
   ctx.cpp
   datum_set.cpp
diff --git a/src/transform.cpp b/src/transform.cpp
index 020d62ea..811e2a6a 100644
--- a/src/transform.cpp
+++ b/src/transform.cpp
@@ -451,7 +451,7 @@ static int pj_apply_vgridshift( PJ *defn,
     if( defn->vgrids_legacy == nullptr )
     {
         defn->vgrids_legacy = new ListOfVGrids;
-        auto vgrids = proj_vgrid_init(defn, "geoidgrids");
+        auto vgrids = pj_vgrid_init(defn, "geoidgrids");
         if( vgrids.empty() )
             return 0;
         *static_cast<ListOfVGrids*>(defn->vgrids_legacy) = std::move(vgrids);
@@ -470,7 +470,7 @@ static int pj_apply_vgridshift( PJ *defn,
         input.phi = y[io];
         input.lam = x[io];
 
-        value = proj_vgrid_value(defn, *static_cast<ListOfVGrids*>(defn->vgrids_legacy), input, 1.0);
+        value = pj_vgrid_value(defn, *static_cast<ListOfVGrids*>(defn->vgrids_legacy), input, 1.0);
 
         if( inverse )
             z[io] -= value;
@@ -896,7 +896,7 @@ int pj_apply_gridshift_2( PJ *defn, int inverse,
     if( defn->hgrids_legacy == nullptr )
     {
         defn->hgrids_legacy = new ListOfHGrids;
-        auto hgrids = proj_hgrid_init(defn, "nadgrids");
+        auto hgrids = pj_hgrid_init(defn, "nadgrids");
         if( hgrids.empty() )
             return 0;
         *static_cast<ListOfHGrids*>(defn->hgrids_legacy) = std::move(hgrids);
@@ -914,7 +914,7 @@ int pj_apply_gridshift_2( PJ *defn, int inverse,
         input.phi = y[io];
         input.lam = x[io];
 
-        auto output = proj_hgrid_apply(defn, *static_cast<ListOfHGrids*>(defn->hgrids_legacy), input, inverse ? PJ_INV : PJ_FWD);
+        auto output = pj_hgrid_apply(defn->ctx, *static_cast<ListOfHGrids*>(defn->hgrids_legacy), input, inverse ? PJ_INV : PJ_FWD);
 
         if ( output.lam != HUGE_VAL )
         {
diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp
index eb109826..993647fc 100644
--- a/src/transformations/deformation.cpp
+++ b/src/transformations/deformation.cpp
@@ -226,8 +226,8 @@ static PJ_XYZ get_grid_shift(PJ* P, const PJ_XYZ& cartesian) {
     }
     else
     {
-        shift.lp    = proj_hgrid_value(P, Q->hgrids, geodetic.lp);
-        shift.enu.u = proj_vgrid_value(P, Q->vgrids, geodetic.lp, 1.0);
+        shift.lp    = pj_hgrid_value(P, Q->hgrids, geodetic.lp);
+        shift.enu.u = pj_vgrid_value(P, Q->vgrids, geodetic.lp, 1.0);
 
         if (proj_errno(P) == PJD_ERR_GRID_AREA)
             proj_log_debug(P, "deformation: coordinate (%.3f, %.3f) outside deformation model",
@@ -425,7 +425,7 @@ PJ *TRANSFORMATION(deformation,1) {
 
     if( has_grids )
     {
-        Q->grids = proj_generic_grid_init(P, "grids");
+        Q->grids = pj_generic_grid_init(P, "grids");
         /* Was gridlist compiled properly? */
         if ( proj_errno(P) ) {
             proj_log_error(P, "deformation: could not find required grid(s).");
@@ -434,13 +434,13 @@ PJ *TRANSFORMATION(deformation,1) {
     }
     else
     {
-        Q->hgrids = proj_hgrid_init(P, "xy_grids");
+        Q->hgrids = pj_hgrid_init(P, "xy_grids");
         if (proj_errno(P)) {
             proj_log_error(P, "deformation: could not find requested xy_grid(s).");
             return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
         }
 
-        Q->vgrids = proj_vgrid_init(P, "z_grids");
+        Q->vgrids = pj_vgrid_init(P, "z_grids");
         if (proj_errno(P)) {
             proj_log_error(P, "deformation: could not find requested z_grid(s).");
             return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
diff --git a/src/transformations/hgridshift.cpp b/src/transformations/hgridshift.cpp
index 24da4dde..122a7ab2 100644
--- a/src/transformations/hgridshift.cpp
+++ b/src/transformations/hgridshift.cpp
@@ -28,7 +28,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
 
     if ( Q->defer_grid_opening ) {
         Q->defer_grid_opening = false;
-        Q->grids = proj_hgrid_init(P, "grids");
+        Q->grids = pj_hgrid_init(P, "grids");
         if ( proj_errno(P) ) {
             return proj_coord_error().xyz;
         }
@@ -37,7 +37,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
     if (!Q->grids.empty()) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        point.lp = proj_hgrid_apply(P, Q->grids, point.lp, PJ_FWD);
+        point.lp = pj_hgrid_apply(P->ctx, Q->grids, point.lp, PJ_FWD);
     }
 
     return point.xyz;
@@ -51,7 +51,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
 
     if ( Q->defer_grid_opening ) {
         Q->defer_grid_opening = false;
-        Q->grids = proj_hgrid_init(P, "grids");
+        Q->grids = pj_hgrid_init(P, "grids");
         if ( proj_errno(P) ) {
             return proj_coord_error().lpz;
         }
@@ -60,7 +60,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
     if (!Q->grids.empty()) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        point.lp = proj_hgrid_apply(P, Q->grids, point.lp, PJ_INV);
+        point.lp = pj_hgrid_apply(P->ctx, Q->grids, point.lp, PJ_INV);
     }
 
     return point.lpz;
@@ -165,7 +165,7 @@ PJ *TRANSFORMATION(hgridshift,0) {
         Q->defer_grid_opening = true;
     }
     else {
-        Q->grids = proj_hgrid_init(P, "grids");
+        Q->grids = pj_hgrid_init(P, "grids");
         /* Was gridlist compiled properly? */
         if ( proj_errno(P) ) {
             proj_log_error(P, "hgridshift: could not find required grid(s).");
diff --git a/src/transformations/vgridshift.cpp b/src/transformations/vgridshift.cpp
index 3e7a015e..121b795a 100644
--- a/src/transformations/vgridshift.cpp
+++ b/src/transformations/vgridshift.cpp
@@ -56,7 +56,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
 
     if ( Q->defer_grid_opening ) {
         Q->defer_grid_opening = false;
-        Q->grids = proj_vgrid_init(P, "grids");
+        Q->grids = pj_vgrid_init(P, "grids");
         deal_with_vertcon_gtx_hack(P);
         if ( proj_errno(P) ) {
             return proj_coord_error().xyz;
@@ -66,7 +66,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
     if (!Q->grids.empty()) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        point.xyz.z += proj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier);
+        point.xyz.z += pj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier);
     }
 
     return point.xyz;
@@ -80,7 +80,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
 
     if ( Q->defer_grid_opening ) {
         Q->defer_grid_opening = false;
-        Q->grids = proj_vgrid_init(P, "grids");
+        Q->grids = pj_vgrid_init(P, "grids");
         deal_with_vertcon_gtx_hack(P);
         if ( proj_errno(P) ) {
             return proj_coord_error().lpz;
@@ -90,7 +90,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
     if (!Q->grids.empty()) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        point.xyz.z -= proj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier);
+        point.xyz.z -= pj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier);
     }
 
     return point.lpz;
@@ -193,7 +193,7 @@ PJ *TRANSFORMATION(vgridshift,0) {
     }
     else {
         /* Build gridlist. P->vgridlist_geoid can be empty if +grids only ask for optional grids. */
-        Q->grids = proj_vgrid_init(P, "grids");
+        Q->grids = pj_vgrid_init(P, "grids");
 
         /* Was gridlist compiled properly? */
         if ( proj_errno(P) ) {
diff --git a/src/transformations/xyzgridshift.cpp b/src/transformations/xyzgridshift.cpp
index a76f3255..3ec3863c 100644
--- a/src/transformations/xyzgridshift.cpp
+++ b/src/transformations/xyzgridshift.cpp
@@ -77,7 +77,7 @@ static bool get_grid_values(PJ* P,
 {
     if ( Q->defer_grid_opening ) {
         Q->defer_grid_opening = false;
-        Q->grids = proj_generic_grid_init(P, "grids");
+        Q->grids = pj_generic_grid_init(P, "grids");
         if ( proj_errno(P) ) {
             return false;
         }
@@ -341,7 +341,7 @@ PJ *TRANSFORMATION(xyzgridshift,0) {
         Q->defer_grid_opening = true;
     }
     else {
-        Q->grids = proj_generic_grid_init(P, "grids");
+        Q->grids = pj_generic_grid_init(P, "grids");
         /* Was gridlist compiled properly? */
         if ( proj_errno(P) ) {
             proj_log_error(P, "xyzgridshift: could not find required grid(s).");