Merge RFC4 (#1865)

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

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

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

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

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

1 files changed, 3402 insertions, 0 deletions

diff --git a/src/grids.cpp b/src/grids.cpp
new file mode 100644
index 00000000..68bf8600
--- /dev/null
+++ b/src/grids.cpp
@@ -0,0 +1,3402 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Grid management
+ * Author:   Even Rouault, <even.rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2000, Frank Warmerdam <warmerdam@pobox.com>
+ * Copyright (c) 2019, Even Rouault, <even.rouault at spatialys.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#ifndef FROM_PROJ_CPP
+#define FROM_PROJ_CPP
+#endif
+#define LRU11_DO_NOT_DEFINE_OUT_OF_CLASS_METHODS
+
+#include "grids.hpp"
+#include "filemanager.hpp"
+#include "proj/internal/internal.hpp"
+#include "proj/internal/lru_cache.hpp"
+#include "proj_internal.h"
+
+#ifdef TIFF_ENABLED
+#include "tiffio.h"
+#endif
+
+#include <algorithm>
+#include <cmath>
+
+NS_PROJ_START
+
+using namespace internal;
+
+/************************************************************************/
+/*                             swap_words()                             */
+/*                                                                      */
+/*      Convert the byte order of the given word(s) in place.           */
+/************************************************************************/
+
+static const int byte_order_test = 1;
+#define IS_LSB (1 == ((const unsigned char *)(&byte_order_test))[0])
+
+static void swap_words(void *dataIn, size_t word_size, size_t word_count)
+
+{
+    unsigned char *data = static_cast<unsigned char *>(dataIn);
+    for (size_t word = 0; word < word_count; word++) {
+        for (size_t i = 0; i < word_size / 2; i++) {
+            unsigned char t;
+
+            t = data[i];
+            data[i] = data[word_size - i - 1];
+            data[word_size - i - 1] = t;
+        }
+
+        data += word_size;
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+bool ExtentAndRes::fullWorldLongitude() const {
+    return eastLon - westLon + resLon >= 2 * M_PI - 1e-10;
+}
+
+// ---------------------------------------------------------------------------
+
+bool ExtentAndRes::contains(const ExtentAndRes &other) const {
+    return other.westLon >= westLon && other.eastLon <= eastLon &&
+           other.southLat >= southLat && other.northLat <= northLat;
+}
+
+// ---------------------------------------------------------------------------
+
+bool ExtentAndRes::intersects(const ExtentAndRes &other) const {
+    return other.westLon < eastLon && westLon <= other.westLon &&
+           other.southLat < northLat && southLat <= other.northLat;
+}
+
+// ---------------------------------------------------------------------------
+
+Grid::Grid(const std::string &name, int widthIn, int heightIn,
+           const ExtentAndRes &extentIn)
+    : m_name(name), m_width(widthIn), m_height(heightIn), m_extent(extentIn) {}
+
+// ---------------------------------------------------------------------------
+
+Grid::~Grid() = default;
+
+// ---------------------------------------------------------------------------
+
+VerticalShiftGrid::VerticalShiftGrid(const std::string &nameIn, int widthIn,
+                                     int heightIn, const ExtentAndRes &extentIn)
+    : Grid(nameIn, widthIn, heightIn, extentIn) {}
+
+// ---------------------------------------------------------------------------
+
+VerticalShiftGrid::~VerticalShiftGrid() = default;
+
+// ---------------------------------------------------------------------------
+
+static ExtentAndRes globalExtent() {
+    ExtentAndRes extent;
+    extent.westLon = -M_PI;
+    extent.southLat = -M_PI / 2;
+    extent.eastLon = M_PI;
+    extent.northLat = M_PI / 2;
+    extent.resLon = M_PI;
+    extent.resLat = M_PI / 2;
+    return extent;
+}
+
+// ---------------------------------------------------------------------------
+
+class NullVerticalShiftGrid : public VerticalShiftGrid {
+
+  public:
+    NullVerticalShiftGrid() : VerticalShiftGrid("null", 3, 3, globalExtent()) {}
+
+    bool isNullGrid() const override { return true; }
+    bool valueAt(int, int, float &out) const override;
+    bool isNodata(float, double) const override { return false; }
+    void reassign_context(PJ_CONTEXT *) override {}
+    bool hasChanged() const override { return false; }
+};
+
+// ---------------------------------------------------------------------------
+
+bool NullVerticalShiftGrid::valueAt(int, int, float &out) const {
+    out = 0.0f;
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+class GTXVerticalShiftGrid : public VerticalShiftGrid {
+    PJ_CONTEXT *m_ctx;
+    std::unique_ptr<File> m_fp;
+
+    GTXVerticalShiftGrid(const GTXVerticalShiftGrid &) = delete;
+    GTXVerticalShiftGrid &operator=(const GTXVerticalShiftGrid &) = delete;
+
+  public:
+    explicit GTXVerticalShiftGrid(PJ_CONTEXT *ctx, std::unique_ptr<File> &&fp,
+                                  const std::string &nameIn, int widthIn,
+                                  int heightIn, const ExtentAndRes &extentIn)
+        : VerticalShiftGrid(nameIn, widthIn, heightIn, extentIn), m_ctx(ctx),
+          m_fp(std::move(fp)) {}
+
+    ~GTXVerticalShiftGrid() override;
+
+    bool valueAt(int x, int y, float &out) const override;
+    bool isNodata(float val, double multiplier) const override;
+
+    static GTXVerticalShiftGrid *open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
+                                      const std::string &name);
+
+    void reassign_context(PJ_CONTEXT *ctx) override {
+        m_ctx = ctx;
+        m_fp->reassign_context(ctx);
+    }
+
+    bool hasChanged() const override { return m_fp->hasChanged(); }
+};
+
+// ---------------------------------------------------------------------------
+
+GTXVerticalShiftGrid::~GTXVerticalShiftGrid() = default;
+
+// ---------------------------------------------------------------------------
+
+GTXVerticalShiftGrid *GTXVerticalShiftGrid::open(PJ_CONTEXT *ctx,
+                                                 std::unique_ptr<File> fp,
+                                                 const std::string &name) {
+    unsigned char header[40];
+
+    /* -------------------------------------------------------------------- */
+    /*      Read the header.                                                */
+    /* -------------------------------------------------------------------- */
+    if (fp->read(header, sizeof(header)) != sizeof(header)) {
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return nullptr;
+    }
+
+    /* -------------------------------------------------------------------- */
+    /*      Regularize fields of interest and extract.                      */
+    /* -------------------------------------------------------------------- */
+    if (IS_LSB) {
+        swap_words(header + 0, 8, 4);
+        swap_words(header + 32, 4, 2);
+    }
+
+    double xorigin, yorigin, xstep, ystep;
+    int rows, columns;
+
+    memcpy(&yorigin, header + 0, 8);
+    memcpy(&xorigin, header + 8, 8);
+    memcpy(&ystep, header + 16, 8);
+    memcpy(&xstep, header + 24, 8);
+
+    memcpy(&rows, header + 32, 4);
+    memcpy(&columns, header + 36, 4);
+
+    if (xorigin < -360 || xorigin > 360 || yorigin < -90 || yorigin > 90) {
+        pj_log(ctx, PJ_LOG_ERROR,
+               "gtx file header has invalid extents, corrupt?");
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return nullptr;
+    }
+
+    /* some GTX files come in 0-360 and we shift them back into the
+       expected -180 to 180 range if possible.  This does not solve
+       problems with grids spanning the dateline. */
+    if (xorigin >= 180.0)
+        xorigin -= 360.0;
+
+    if (xorigin >= 0.0 && xorigin + xstep * columns > 180.0) {
+        pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+               "This GTX spans the dateline!  This will cause problems.");
+    }
+
+    ExtentAndRes extent;
+    extent.westLon = xorigin * DEG_TO_RAD;
+    extent.southLat = yorigin * DEG_TO_RAD;
+    extent.resLon = xstep * DEG_TO_RAD;
+    extent.resLat = ystep * DEG_TO_RAD;
+    extent.eastLon = (xorigin + xstep * (columns - 1)) * DEG_TO_RAD;
+    extent.northLat = (yorigin + ystep * (rows - 1)) * DEG_TO_RAD;
+
+    return new GTXVerticalShiftGrid(ctx, std::move(fp), name, columns, rows,
+                                    extent);
+}
+
+// ---------------------------------------------------------------------------
+
+bool GTXVerticalShiftGrid::valueAt(int x, int y, float &out) const {
+    assert(x >= 0 && y >= 0 && x < m_width && y < m_height);
+
+    m_fp->seek(40 + sizeof(float) * (y * m_width + x));
+    if (m_fp->read(&out, sizeof(out)) != sizeof(out)) {
+        pj_ctx_set_errno(m_ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return false;
+    }
+    if (IS_LSB) {
+        swap_words(&out, sizeof(float), 1);
+    }
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+bool GTXVerticalShiftGrid::isNodata(float val, double multiplier) const {
+    /* nodata?  */
+    /* GTX official nodata value if  -88.88880f, but some grids also */
+    /* use other  big values for nodata (e.g naptrans2008.gtx has */
+    /* nodata values like -2147479936), so test them too */
+    return val * multiplier > 1000 || val * multiplier < -1000 ||
+           val == -88.88880f;
+}
+
+// ---------------------------------------------------------------------------
+
+VerticalShiftGridSet::VerticalShiftGridSet() = default;
+
+// ---------------------------------------------------------------------------
+
+VerticalShiftGridSet::~VerticalShiftGridSet() = default;
+
+// ---------------------------------------------------------------------------
+
+static bool IsTIFF(size_t header_size, const unsigned char *header) {
+    // Test combinations of signature for ClassicTIFF/BigTIFF little/big endian
+    return header_size >= 4 && (((header[0] == 'I' && header[1] == 'I') ||
+                                 (header[0] == 'M' && header[1] == 'M')) &&
+                                ((header[2] == 0x2A && header[3] == 0) ||
+                                 (header[3] == 0x2A && header[2] == 0) ||
+                                 (header[2] == 0x2B && header[3] == 0) ||
+                                 (header[3] == 0x2B && header[2] == 0)));
+}
+
+#ifdef TIFF_ENABLED
+
+// ---------------------------------------------------------------------------
+
+enum class TIFFDataType { Int16, UInt16, Int32, UInt32, Float32, Float64 };
+
+// ---------------------------------------------------------------------------
+
+constexpr uint16 TIFFTAG_GEOPIXELSCALE = 33550;
+constexpr uint16 TIFFTAG_GEOTIEPOINTS = 33922;
+constexpr uint16 TIFFTAG_GEOTRANSMATRIX = 34264;
+constexpr uint16 TIFFTAG_GEOKEYDIRECTORY = 34735;
+constexpr uint16 TIFFTAG_GEODOUBLEPARAMS = 34736;
+constexpr uint16 TIFFTAG_GEOASCIIPARAMS = 34737;
+constexpr uint16 TIFFTAG_GDAL_METADATA = 42112;
+constexpr uint16 TIFFTAG_GDAL_NODATA = 42113;
+
+// ---------------------------------------------------------------------------
+
+class BlockCache {
+  public:
+    void insert(uint32 ifdIdx, uint32 blockNumber,
+                const std::vector<unsigned char> &data);
+    std::shared_ptr<std::vector<unsigned char>> get(uint32 ifdIdx,
+                                                    uint32 blockNumber);
+
+  private:
+    struct Key {
+        uint32 ifdIdx;
+        uint32 blockNumber;
+
+        Key(uint32 ifdIdxIn, uint32 blockNumberIn)
+            : ifdIdx(ifdIdxIn), blockNumber(blockNumberIn) {}
+        bool operator==(const Key &other) const {
+            return ifdIdx == other.ifdIdx && blockNumber == other.blockNumber;
+        }
+    };
+
+    struct KeyHasher {
+        std::size_t operator()(const Key &k) const {
+            return k.ifdIdx ^ (k.blockNumber << 16) ^ (k.blockNumber >> 16);
+        }
+    };
+
+    static constexpr int NUM_BLOCKS_AT_CROSSING_TILES = 4;
+    static constexpr int MAX_SAMPLE_COUNT = 3;
+    lru11::Cache<
+        Key, std::shared_ptr<std::vector<unsigned char>>, lru11::NullLock,
+        std::unordered_map<
+            Key,
+            typename std::list<lru11::KeyValuePair<
+                Key, std::shared_ptr<std::vector<unsigned char>>>>::iterator,
+            KeyHasher>>
+        cache_{NUM_BLOCKS_AT_CROSSING_TILES * MAX_SAMPLE_COUNT};
+};
+
+// ---------------------------------------------------------------------------
+
+void BlockCache::insert(uint32 ifdIdx, uint32 blockNumber,
+                        const std::vector<unsigned char> &data) {
+    cache_.insert(Key(ifdIdx, blockNumber),
+                  std::make_shared<std::vector<unsigned char>>(data));
+}
+
+// ---------------------------------------------------------------------------
+
+std::shared_ptr<std::vector<unsigned char>>
+BlockCache::get(uint32 ifdIdx, uint32 blockNumber) {
+    std::shared_ptr<std::vector<unsigned char>> ret;
+    cache_.tryGet(Key(ifdIdx, blockNumber), ret);
+    return ret;
+}
+
+// ---------------------------------------------------------------------------
+
+class GTiffGrid : public Grid {
+    PJ_CONTEXT *m_ctx;   // owned by the belonging GTiffDataset
+    TIFF *m_hTIFF;       // owned by the belonging GTiffDataset
+    BlockCache &m_cache; // owned by the belonging GTiffDataset
+    File *m_fp;          // owned by the belonging GTiffDataset
+    uint32 m_ifdIdx;
+    TIFFDataType m_dt;
+    uint16 m_samplesPerPixel;
+    uint16 m_planarConfig;
+    bool m_bottomUp;
+    toff_t m_dirOffset;
+    bool m_tiled;
+    uint32 m_blockWidth = 0;
+    uint32 m_blockHeight = 0;
+    mutable std::vector<unsigned char> m_buffer{};
+    unsigned m_blocksPerRow = 0;
+    unsigned m_blocksPerCol = 0;
+    std::map<int, double> m_mapOffset{};
+    std::map<int, double> m_mapScale{};
+    std::map<std::pair<int, std::string>, std::string> m_metadata{};
+    bool m_hasNodata = false;
+    float m_noData = 0.0f;
+    uint32 m_subfileType = 0;
+
+    GTiffGrid(const GTiffGrid &) = delete;
+    GTiffGrid &operator=(const GTiffGrid &) = delete;
+
+    void getScaleOffset(double &scale, double &offset, uint16 sample) const;
+
+    template <class T>
+    float readValue(const std::vector<unsigned char> &buffer,
+                    uint32 offsetInBlock, uint16 sample) const;
+
+  public:
+    GTiffGrid(PJ_CONTEXT *ctx, TIFF *hTIFF, BlockCache &cache, File *fp,
+              uint32 ifdIdx, const std::string &nameIn, int widthIn,
+              int heightIn, const ExtentAndRes &extentIn, TIFFDataType dtIn,
+              uint16 samplesPerPixelIn, uint16 planarConfig, bool bottomUpIn);
+
+    ~GTiffGrid() override;
+
+    uint16 samplesPerPixel() const { return m_samplesPerPixel; }
+
+    bool valueAt(uint16 sample, int x, int y, float &out) const;
+
+    bool isNodata(float val) const;
+
+    std::string metadataItem(const std::string &key, int sample = -1) const;
+
+    uint32 subfileType() const { return m_subfileType; }
+
+    void reassign_context(PJ_CONTEXT *ctx) { m_ctx = ctx; }
+
+    bool hasChanged() const override { return m_fp->hasChanged(); }
+};
+
+// ---------------------------------------------------------------------------
+
+GTiffGrid::GTiffGrid(PJ_CONTEXT *ctx, TIFF *hTIFF, BlockCache &cache, File *fp,
+                     uint32 ifdIdx, const std::string &nameIn, int widthIn,
+                     int heightIn, const ExtentAndRes &extentIn,
+                     TIFFDataType dtIn, uint16 samplesPerPixelIn,
+                     uint16 planarConfig, bool bottomUpIn)
+    : Grid(nameIn, widthIn, heightIn, extentIn), m_ctx(ctx), m_hTIFF(hTIFF),
+      m_cache(cache), m_fp(fp), m_ifdIdx(ifdIdx), m_dt(dtIn),
+      m_samplesPerPixel(samplesPerPixelIn), m_planarConfig(planarConfig),
+      m_bottomUp(bottomUpIn), m_dirOffset(TIFFCurrentDirOffset(hTIFF)),
+      m_tiled(TIFFIsTiled(hTIFF) != 0) {
+
+    if (m_tiled) {
+        TIFFGetField(m_hTIFF, TIFFTAG_TILEWIDTH, &m_blockWidth);
+        TIFFGetField(m_hTIFF, TIFFTAG_TILELENGTH, &m_blockHeight);
+    } else {
+        m_blockWidth = widthIn;
+        TIFFGetField(m_hTIFF, TIFFTAG_ROWSPERSTRIP, &m_blockHeight);
+        if (m_blockHeight > static_cast<unsigned>(m_height))
+            m_blockHeight = m_height;
+    }
+
+    TIFFGetField(m_hTIFF, TIFFTAG_SUBFILETYPE, &m_subfileType);
+
+    m_blocksPerRow = (m_width + m_blockWidth - 1) / m_blockWidth;
+    m_blocksPerCol = (m_height + m_blockHeight - 1) / m_blockHeight;
+
+    const char *text = nullptr;
+    // Poor-man XML parsing of TIFFTAG_GDAL_METADATA tag. Hopefully good
+    // enough for our purposes.
+    if (TIFFGetField(m_hTIFF, TIFFTAG_GDAL_METADATA, &text)) {
+        const char *ptr = text;
+        while (true) {
+            ptr = strstr(ptr, "<Item ");
+            if (ptr == nullptr)
+                break;
+            const char *endTag = strchr(ptr, '>');
+            if (endTag == nullptr)
+                break;
+            const char *endValue = strchr(endTag, '<');
+            if (endValue == nullptr)
+                break;
+
+            std::string tag;
+            tag.append(ptr, endTag - ptr);
+
+            std::string value;
+            value.append(endTag + 1, endValue - (endTag + 1));
+
+            std::string name;
+            auto namePos = tag.find("name=\"");
+            if (namePos == std::string::npos)
+                break;
+            {
+                namePos += strlen("name=\"");
+                const auto endQuote = tag.find('"', namePos);
+                if (endQuote == std::string::npos)
+                    break;
+                name = tag.substr(namePos, endQuote - namePos);
+            }
+
+            const auto samplePos = tag.find("sample=\"");
+            int sample = -1;
+            if (samplePos != std::string::npos) {
+                sample = atoi(tag.c_str() + samplePos + strlen("sample=\""));
+            }
+
+            m_metadata[std::pair<int, std::string>(sample, name)] = value;
+
+            auto rolePos = tag.find("role=\"");
+            if (rolePos != std::string::npos) {
+                rolePos += strlen("role=\"");
+                const auto endQuote = tag.find('"', rolePos);
+                if (endQuote == std::string::npos)
+                    break;
+                const auto role = tag.substr(rolePos, endQuote - rolePos);
+                if (role == "offset") {
+                    if (sample >= 0) {
+                        try {
+                            m_mapOffset[sample] = c_locale_stod(value);
+                        } catch (const std::exception &) {
+                        }
+                    }
+                } else if (role == "scale") {
+                    if (sample >= 0) {
+                        try {
+                            m_mapScale[sample] = c_locale_stod(value);
+                        } catch (const std::exception &) {
+                        }
+                    }
+                }
+            }
+
+            ptr = endValue + 1;
+        }
+    }
+
+    if (TIFFGetField(m_hTIFF, TIFFTAG_GDAL_NODATA, &text)) {
+        try {
+            m_noData = static_cast<float>(c_locale_stod(text));
+            m_hasNodata = true;
+        } catch (const std::exception &) {
+        }
+    }
+
+    auto oIter = m_metadata.find(std::pair<int, std::string>(-1, "grid_name"));
+    if (oIter != m_metadata.end()) {
+        m_name += ", " + oIter->second;
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+GTiffGrid::~GTiffGrid() = default;
+
+// ---------------------------------------------------------------------------
+
+void GTiffGrid::getScaleOffset(double &scale, double &offset,
+                               uint16 sample) const {
+    {
+        auto iter = m_mapScale.find(sample);
+        if (iter != m_mapScale.end())
+            scale = iter->second;
+    }
+
+    {
+        auto iter = m_mapOffset.find(sample);
+        if (iter != m_mapOffset.end())
+            offset = iter->second;
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+template <class T>
+float GTiffGrid::readValue(const std::vector<unsigned char> &buffer,
+                           uint32 offsetInBlock, uint16 sample) const {
+    const auto ptr = reinterpret_cast<const T *>(buffer.data());
+    assert(offsetInBlock < buffer.size() / sizeof(T));
+    const auto val = ptr[offsetInBlock];
+    if (!m_hasNodata || static_cast<float>(val) != m_noData) {
+        double scale = 1;
+        double offset = 0;
+        getScaleOffset(scale, offset, sample);
+        return static_cast<float>(val * scale + offset);
+    } else {
+        return static_cast<float>(val);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+bool GTiffGrid::valueAt(uint16 sample, int x, int yFromBottom,
+                        float &out) const {
+    assert(x >= 0 && yFromBottom >= 0 && x < m_width && yFromBottom < m_height);
+    assert(sample < m_samplesPerPixel);
+
+    const int blockX = x / m_blockWidth;
+
+    // All non-TIFF grids have the first rows in the file being the one
+    // corresponding to the southern-most row. In GeoTIFF, the convention is
+    // *generally* different (when m_bottomUp == false), TIFF being an
+    // image-oriented image. If m_bottomUp == true, then we had GeoTIFF hints
+    // that the first row of the image is the southern-most.
+    const int yTIFF = m_bottomUp ? yFromBottom : m_height - 1 - yFromBottom;
+    const int blockY = yTIFF / m_blockHeight;
+
+    uint32 blockId = blockY * m_blocksPerRow + blockX;
+    if (m_planarConfig == PLANARCONFIG_SEPARATE) {
+        blockId += sample * m_blocksPerCol * m_blocksPerRow;
+    }
+
+    auto cachedBuffer = m_cache.get(m_ifdIdx, blockId);
+    std::vector<unsigned char> *pBuffer = &m_buffer;
+    if (cachedBuffer != nullptr) {
+        // Safe as we don't access the cache before pBuffer is used
+        pBuffer = cachedBuffer.get();
+    } else {
+        if (TIFFCurrentDirOffset(m_hTIFF) != m_dirOffset &&
+            !TIFFSetSubDirectory(m_hTIFF, m_dirOffset)) {
+            return false;
+        }
+        if (m_buffer.empty()) {
+            const auto blockSize = static_cast<size_t>(
+                m_tiled ? TIFFTileSize64(m_hTIFF) : TIFFStripSize64(m_hTIFF));
+            try {
+                m_buffer.resize(blockSize);
+            } catch (const std::exception &e) {
+                pj_log(m_ctx, PJ_LOG_ERROR, "Exception %s", e.what());
+                return false;
+            }
+        }
+
+        if (m_tiled) {
+            if (TIFFReadEncodedTile(m_hTIFF, blockId, m_buffer.data(),
+                                    m_buffer.size()) == -1) {
+                return false;
+            }
+        } else {
+            if (TIFFReadEncodedStrip(m_hTIFF, blockId, m_buffer.data(),
+                                     m_buffer.size()) == -1) {
+                return false;
+            }
+        }
+
+        try {
+            m_cache.insert(m_ifdIdx, blockId, m_buffer);
+        } catch (const std::exception &e) {
+            // Should normally not happen
+            pj_log(m_ctx, PJ_LOG_ERROR, "Exception %s", e.what());
+        }
+    }
+
+    uint32 offsetInBlock =
+        (x % m_blockWidth) + (yTIFF % m_blockHeight) * m_blockWidth;
+    if (m_planarConfig == PLANARCONFIG_CONTIG)
+        offsetInBlock = offsetInBlock * m_samplesPerPixel + sample;
+
+    switch (m_dt) {
+    case TIFFDataType::Int16:
+        out = readValue<short>(*pBuffer, offsetInBlock, sample);
+        break;
+
+    case TIFFDataType::UInt16:
+        out = readValue<unsigned short>(*pBuffer, offsetInBlock, sample);
+        break;
+
+    case TIFFDataType::Int32:
+        out = readValue<int>(*pBuffer, offsetInBlock, sample);
+        break;
+
+    case TIFFDataType::UInt32:
+        out = readValue<unsigned int>(*pBuffer, offsetInBlock, sample);
+        break;
+
+    case TIFFDataType::Float32:
+        out = readValue<float>(*pBuffer, offsetInBlock, sample);
+        break;
+
+    case TIFFDataType::Float64:
+        out = readValue<double>(*pBuffer, offsetInBlock, sample);
+        break;
+    }
+
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+bool GTiffGrid::isNodata(float val) const {
+    return (m_hasNodata && val == m_noData) || std::isnan(val);
+}
+
+// ---------------------------------------------------------------------------
+
+std::string GTiffGrid::metadataItem(const std::string &key, int sample) const {
+    auto iter = m_metadata.find(std::pair<int, std::string>(sample, key));
+    if (iter == m_metadata.end()) {
+        return std::string();
+    }
+    return iter->second;
+}
+
+// ---------------------------------------------------------------------------
+
+class GTiffDataset {
+    PJ_CONTEXT *m_ctx;
+    std::unique_ptr<File> m_fp;
+    TIFF *m_hTIFF = nullptr;
+    bool m_hasNextGrid = false;
+    uint32 m_ifdIdx = 0;
+    toff_t m_nextDirOffset = 0;
+    std::string m_filename{};
+    BlockCache m_cache{};
+
+    GTiffDataset(const GTiffDataset &) = delete;
+    GTiffDataset &operator=(const GTiffDataset &) = delete;
+
+    // libtiff I/O routines
+    static tsize_t tiffReadProc(thandle_t fd, tdata_t buf, tsize_t size) {
+        GTiffDataset *self = static_cast<GTiffDataset *>(fd);
+        return self->m_fp->read(buf, size);
+    }
+
+    static tsize_t tiffWriteProc(thandle_t, tdata_t, tsize_t) {
+        assert(false);
+        return 0;
+    }
+
+    static toff_t tiffSeekProc(thandle_t fd, toff_t off, int whence) {
+        GTiffDataset *self = static_cast<GTiffDataset *>(fd);
+        if (self->m_fp->seek(off, whence))
+            return static_cast<toff_t>(self->m_fp->tell());
+        else
+            return static_cast<toff_t>(-1);
+    }
+
+    static int tiffCloseProc(thandle_t) {
+        // done in destructor
+        return 0;
+    }
+
+    static toff_t tiffSizeProc(thandle_t fd) {
+        GTiffDataset *self = static_cast<GTiffDataset *>(fd);
+        const auto old_off = self->m_fp->tell();
+        self->m_fp->seek(0, SEEK_END);
+        const auto file_size = static_cast<toff_t>(self->m_fp->tell());
+        self->m_fp->seek(old_off);
+        return file_size;
+    }
+
+    static int tiffMapProc(thandle_t, tdata_t *, toff_t *) { return (0); }
+
+    static void tiffUnmapProc(thandle_t, tdata_t, toff_t) {}
+
+  public:
+    GTiffDataset(PJ_CONTEXT *ctx, std::unique_ptr<File> &&fp)
+        : m_ctx(ctx), m_fp(std::move(fp)) {}
+    virtual ~GTiffDataset();
+
+    bool openTIFF(const std::string &filename);
+
+    std::unique_ptr<GTiffGrid> nextGrid();
+
+    void reassign_context(PJ_CONTEXT *ctx) {
+        m_ctx = ctx;
+        m_fp->reassign_context(ctx);
+    }
+};
+
+// ---------------------------------------------------------------------------
+
+GTiffDataset::~GTiffDataset() {
+    if (m_hTIFF)
+        TIFFClose(m_hTIFF);
+}
+
+// ---------------------------------------------------------------------------
+class OneTimeTIFFTagInit {
+
+    static TIFFExtendProc ParentExtender;
+
+    // Function called by libtiff when initializing a TIFF directory
+    static void GTiffTagExtender(TIFF *tif) {
+        static const TIFFFieldInfo xtiffFieldInfo[] = {
+            // GeoTIFF tags
+            {TIFFTAG_GEOPIXELSCALE, -1, -1, TIFF_DOUBLE, FIELD_CUSTOM, TRUE,
+             TRUE, const_cast<char *>("GeoPixelScale")},
+            {TIFFTAG_GEOTIEPOINTS, -1, -1, TIFF_DOUBLE, FIELD_CUSTOM, TRUE,
+             TRUE, const_cast<char *>("GeoTiePoints")},
+            {TIFFTAG_GEOTRANSMATRIX, -1, -1, TIFF_DOUBLE, FIELD_CUSTOM, TRUE,
+             TRUE, const_cast<char *>("GeoTransformationMatrix")},
+
+            {TIFFTAG_GEOKEYDIRECTORY, -1, -1, TIFF_SHORT, FIELD_CUSTOM, TRUE,
+             TRUE, const_cast<char *>("GeoKeyDirectory")},
+            {TIFFTAG_GEODOUBLEPARAMS, -1, -1, TIFF_DOUBLE, FIELD_CUSTOM, TRUE,
+             TRUE, const_cast<char *>("GeoDoubleParams")},
+            {TIFFTAG_GEOASCIIPARAMS, -1, -1, TIFF_ASCII, FIELD_CUSTOM, TRUE,
+             FALSE, const_cast<char *>("GeoASCIIParams")},
+
+            // GDAL tags
+            {TIFFTAG_GDAL_METADATA, -1, -1, TIFF_ASCII, FIELD_CUSTOM, TRUE,
+             FALSE, const_cast<char *>("GDALMetadata")},
+            {TIFFTAG_GDAL_NODATA, -1, -1, TIFF_ASCII, FIELD_CUSTOM, TRUE, FALSE,
+             const_cast<char *>("GDALNoDataValue")},
+
+        };
+
+        if (ParentExtender)
+            (*ParentExtender)(tif);
+
+        TIFFMergeFieldInfo(tif, xtiffFieldInfo,
+                           sizeof(xtiffFieldInfo) / sizeof(xtiffFieldInfo[0]));
+    }
+
+  public:
+    OneTimeTIFFTagInit() {
+        assert(ParentExtender == nullptr);
+        // Install our TIFF tag extender
+        ParentExtender = TIFFSetTagExtender(GTiffTagExtender);
+    }
+};
+
+TIFFExtendProc OneTimeTIFFTagInit::ParentExtender = nullptr;
+
+// ---------------------------------------------------------------------------
+
+bool GTiffDataset::openTIFF(const std::string &filename) {
+    static OneTimeTIFFTagInit oneTimeTIFFTagInit;
+    m_hTIFF =
+        TIFFClientOpen(filename.c_str(), "r", static_cast<thandle_t>(this),
+                       GTiffDataset::tiffReadProc, GTiffDataset::tiffWriteProc,
+                       GTiffDataset::tiffSeekProc, GTiffDataset::tiffCloseProc,
+                       GTiffDataset::tiffSizeProc, GTiffDataset::tiffMapProc,
+                       GTiffDataset::tiffUnmapProc);
+
+    m_filename = filename;
+    m_hasNextGrid = true;
+    return m_hTIFF != nullptr;
+}
+// ---------------------------------------------------------------------------
+
+std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() {
+    if (!m_hasNextGrid)
+        return nullptr;
+    if (m_nextDirOffset) {
+        TIFFSetSubDirectory(m_hTIFF, m_nextDirOffset);
+    }
+
+    uint32 width = 0;
+    uint32 height = 0;
+    TIFFGetField(m_hTIFF, TIFFTAG_IMAGEWIDTH, &width);
+    TIFFGetField(m_hTIFF, TIFFTAG_IMAGELENGTH, &height);
+    if (width == 0 || height == 0 || width > INT_MAX || height > INT_MAX) {
+        pj_log(m_ctx, PJ_LOG_ERROR, "Invalid image size");
+        return nullptr;
+    }
+
+    uint16 samplesPerPixel = 0;
+    if (!TIFFGetField(m_hTIFF, TIFFTAG_SAMPLESPERPIXEL, &samplesPerPixel)) {
+        pj_log(m_ctx, PJ_LOG_ERROR, "Missing SamplesPerPixel tag");
+        return nullptr;
+    }
+    if (samplesPerPixel == 0) {
+        pj_log(m_ctx, PJ_LOG_ERROR, "Invalid SamplesPerPixel value");
+        return nullptr;
+    }
+
+    uint16 bitsPerSample = 0;
+    if (!TIFFGetField(m_hTIFF, TIFFTAG_BITSPERSAMPLE, &bitsPerSample)) {
+        pj_log(m_ctx, PJ_LOG_ERROR, "Missing BitsPerSample tag");
+        return nullptr;
+    }
+
+    uint16 planarConfig = 0;
+    if (!TIFFGetField(m_hTIFF, TIFFTAG_PLANARCONFIG, &planarConfig)) {
+        pj_log(m_ctx, PJ_LOG_ERROR, "Missing PlanarConfig tag");
+        return nullptr;
+    }
+
+    uint16 sampleFormat = 0;
+    if (!TIFFGetField(m_hTIFF, TIFFTAG_SAMPLEFORMAT, &sampleFormat)) {
+        pj_log(m_ctx, PJ_LOG_ERROR, "Missing SampleFormat tag");
+        return nullptr;
+    }
+
+    TIFFDataType dt;
+    if (sampleFormat == SAMPLEFORMAT_INT && bitsPerSample == 16)
+        dt = TIFFDataType::Int16;
+    else if (sampleFormat == SAMPLEFORMAT_UINT && bitsPerSample == 16)
+        dt = TIFFDataType::UInt16;
+    else if (sampleFormat == SAMPLEFORMAT_INT && bitsPerSample == 32)
+        dt = TIFFDataType::Int32;
+    else if (sampleFormat == SAMPLEFORMAT_UINT && bitsPerSample == 32)
+        dt = TIFFDataType::UInt32;
+    else if (sampleFormat == SAMPLEFORMAT_IEEEFP && bitsPerSample == 32)
+        dt = TIFFDataType::Float32;
+    else if (sampleFormat == SAMPLEFORMAT_IEEEFP && bitsPerSample == 64)
+        dt = TIFFDataType::Float64;
+    else {
+        pj_log(
+            m_ctx, PJ_LOG_ERROR,
+            "Unsupported combination of SampleFormat and BitsPerSample values");
+        return nullptr;
+    }
+
+    uint16 photometric = PHOTOMETRIC_MINISBLACK;
+    if (!TIFFGetField(m_hTIFF, TIFFTAG_PHOTOMETRIC, &photometric))
+        photometric = PHOTOMETRIC_MINISBLACK;
+    if (photometric != PHOTOMETRIC_MINISBLACK) {
+        pj_log(m_ctx, PJ_LOG_ERROR, "Unsupported Photometric value");
+        return nullptr;
+    }
+
+    uint16 compression = COMPRESSION_NONE;
+    if (!TIFFGetField(m_hTIFF, TIFFTAG_COMPRESSION, &compression))
+        compression = COMPRESSION_NONE;
+
+    if (compression != COMPRESSION_NONE &&
+        !TIFFIsCODECConfigured(compression)) {
+        pj_log(m_ctx, PJ_LOG_ERROR,
+               "Cannot open TIFF file due to missing codec.");
+        return nullptr;
+    }
+    // We really don't want to try dealing with old-JPEG images
+    if (compression == COMPRESSION_OJPEG) {
+        pj_log(m_ctx, PJ_LOG_ERROR, "Unsupported compression method.");
+        return nullptr;
+    }
+
+    const auto blockSize = TIFFIsTiled(m_hTIFF) ? TIFFTileSize64(m_hTIFF)
+                                                : TIFFStripSize64(m_hTIFF);
+    if (blockSize == 0 || blockSize > 64 * 1024 * 2014) {
+        pj_log(m_ctx, PJ_LOG_ERROR, "Unsupported block size.");
+        return nullptr;
+    }
+
+    unsigned short count = 0;
+    unsigned short *geokeys = nullptr;
+    bool pixelIsArea = false;
+    if (!TIFFGetField(m_hTIFF, TIFFTAG_GEOKEYDIRECTORY, &count, &geokeys)) {
+        pj_log(m_ctx, PJ_LOG_DEBUG_MINOR, "No GeoKeys tag");
+    } else {
+        if (count < 4 || (count % 4) != 0) {
+            pj_log(m_ctx, PJ_LOG_ERROR,
+                   "Wrong number of values in GeoKeys tag");
+            return nullptr;
+        }
+
+        if (geokeys[0] != 1) {
+            pj_log(m_ctx, PJ_LOG_ERROR, "Unsupported GeoTIFF major version");
+            return nullptr;
+        }
+        // We only know that we support GeoTIFF 1.0 and 1.1 at that time
+        if (geokeys[1] != 1 || geokeys[2] > 1) {
+            pj_log(m_ctx, PJ_LOG_DEBUG_MINOR,
+                   "GeoTIFF %d.%d possibly not handled", geokeys[1],
+                   geokeys[2]);
+        }
+
+        for (unsigned int i = 4; i + 3 < count; i += 4) {
+            constexpr unsigned short GTModelTypeGeoKey = 1024;
+            constexpr unsigned short ModelTypeGeographic = 2;
+
+            constexpr unsigned short GTRasterTypeGeoKey = 1025;
+            constexpr unsigned short RasterPixelIsArea = 1;
+            // constexpr unsigned short RasterPixelIsPoint = 2;
+
+            if (geokeys[i] == GTModelTypeGeoKey) {
+                if (geokeys[i + 3] != ModelTypeGeographic) {
+                    pj_log(m_ctx, PJ_LOG_ERROR, "Only GTModelTypeGeoKey = "
+                                                "ModelTypeGeographic is "
+                                                "supported");
+                    return nullptr;
+                }
+            } else if (geokeys[i] == GTRasterTypeGeoKey) {
+                if (geokeys[i + 3] == RasterPixelIsArea) {
+                    pixelIsArea = true;
+                }
+            }
+        }
+    }
+
+    double hRes = 0;
+    double vRes = 0;
+    double westLon = 0;
+    double northLat = 0;
+
+    double *matrix = nullptr;
+    if (TIFFGetField(m_hTIFF, TIFFTAG_GEOTRANSMATRIX, &count, &matrix) &&
+        count == 16) {
+        // If using GDAL to produce a bottom-up georeferencing, it will produce
+        // a GeoTransformationMatrix, since negative values in GeoPixelScale
+        // have historically been implementation bugs.
+        if (matrix[1] != 0 || matrix[4] != 0) {
+            pj_log(m_ctx, PJ_LOG_ERROR, "Rotational terms not supported in "
+                                        "GeoTransformationMatrix tag");
+            return nullptr;
+        }
+
+        westLon = matrix[3];
+        hRes = matrix[0];
+        northLat = matrix[7];
+        vRes = -matrix[5]; // negation to simulate GeoPixelScale convention
+    } else {
+        double *geopixelscale = nullptr;
+        if (TIFFGetField(m_hTIFF, TIFFTAG_GEOPIXELSCALE, &count,
+                         &geopixelscale) != 1) {
+            pj_log(m_ctx, PJ_LOG_ERROR, "No GeoPixelScale tag");
+            return nullptr;
+        }
+        if (count != 3) {
+            pj_log(m_ctx, PJ_LOG_ERROR,
+                   "Wrong number of values in GeoPixelScale tag");
+            return nullptr;
+        }
+        hRes = geopixelscale[0];
+        vRes = geopixelscale[1];
+
+        double *geotiepoints = nullptr;
+        if (TIFFGetField(m_hTIFF, TIFFTAG_GEOTIEPOINTS, &count,
+                         &geotiepoints) != 1) {
+            pj_log(m_ctx, PJ_LOG_ERROR, "No GeoTiePoints tag");
+            return nullptr;
+        }
+        if (count != 6) {
+            pj_log(m_ctx, PJ_LOG_ERROR,
+                   "Wrong number of values in GeoTiePoints tag");
+            return nullptr;
+        }
+
+        westLon = geotiepoints[3] - geotiepoints[0] * hRes;
+        northLat = geotiepoints[4] + geotiepoints[1] * vRes;
+    }
+
+    if (pixelIsArea) {
+        westLon += 0.5 * hRes;
+        northLat -= 0.5 * vRes;
+    }
+
+    ExtentAndRes extent;
+    extent.westLon = westLon * DEG_TO_RAD;
+    extent.northLat = northLat * DEG_TO_RAD;
+    extent.resLon = hRes * DEG_TO_RAD;
+    extent.resLat = fabs(vRes) * DEG_TO_RAD;
+    extent.eastLon = (westLon + hRes * (width - 1)) * DEG_TO_RAD;
+    extent.southLat = (northLat - vRes * (height - 1)) * DEG_TO_RAD;
+
+    if (vRes < 0) {
+        std::swap(extent.northLat, extent.southLat);
+    }
+
+    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(m_ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s",
+               m_filename.c_str());
+        return nullptr;
+    }
+
+    auto ret = std::unique_ptr<GTiffGrid>(new GTiffGrid(
+        m_ctx, m_hTIFF, m_cache, m_fp.get(), m_ifdIdx, m_filename, width,
+        height, extent, dt, samplesPerPixel, planarConfig, vRes < 0));
+    m_ifdIdx++;
+    m_hasNextGrid = TIFFReadDirectory(m_hTIFF) != 0;
+    m_nextDirOffset = TIFFCurrentDirOffset(m_hTIFF);
+    return ret;
+}
+
+// ---------------------------------------------------------------------------
+
+class GTiffVGridShiftSet : public VerticalShiftGridSet {
+
+    std::unique_ptr<GTiffDataset> m_GTiffDataset;
+
+    GTiffVGridShiftSet(PJ_CONTEXT *ctx, std::unique_ptr<File> &&fp)
+        : m_GTiffDataset(new GTiffDataset(ctx, std::move(fp))) {}
+
+  public:
+    ~GTiffVGridShiftSet() override;
+
+    static std::unique_ptr<GTiffVGridShiftSet>
+    open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
+         const std::string &filename);
+
+    void reassign_context(PJ_CONTEXT *ctx) override {
+        VerticalShiftGridSet::reassign_context(ctx);
+        if (m_GTiffDataset) {
+            m_GTiffDataset->reassign_context(ctx);
+        }
+    }
+
+    bool reopen(PJ_CONTEXT *ctx) override {
+        pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Grid %s has changed. Re-loading it",
+               m_name.c_str());
+        m_grids.clear();
+        m_GTiffDataset.reset();
+        auto fp = FileManager::open_resource_file(ctx, m_name.c_str());
+        if (!fp) {
+            return false;
+        }
+        auto newGS = open(ctx, std::move(fp), m_name);
+        if (newGS) {
+            m_grids = std::move(newGS->m_grids);
+            m_GTiffDataset = std::move(newGS->m_GTiffDataset);
+        }
+        return !m_grids.empty();
+    }
+};
+
+#endif // TIFF_ENABLED
+
+// ---------------------------------------------------------------------------
+
+template <class GridType, class GenericGridType>
+static void
+insertIntoHierarchy(PJ_CONTEXT *ctx, std::unique_ptr<GridType> &&grid,
+                    const std::string &gridName, const std::string &parentName,
+                    std::vector<std::unique_ptr<GenericGridType>> &topGrids,
+                    std::map<std::string, GridType *> &mapGrids) {
+    const auto &extent = grid->extentAndRes();
+
+    // If we have one or both of grid_name and parent_grid_name, try to use
+    // the names to recreate the hiearchy
+    if (!gridName.empty()) {
+        if (mapGrids.find(gridName) != mapGrids.end()) {
+            pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Several grids called %s found!",
+                   gridName.c_str());
+        }
+        mapGrids[gridName] = grid.get();
+    }
+    bool gridInserted = false;
+    if (!parentName.empty()) {
+        auto iter = mapGrids.find(parentName);
+        if (iter == mapGrids.end()) {
+            pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+                   "Grid %s refers to non-existing parent %s. "
+                   "Using bounding-box method.",
+                   gridName.c_str(), parentName.c_str());
+        } else {
+            if (iter->second->extentAndRes().contains(extent)) {
+                iter->second->m_children.emplace_back(std::move(grid));
+                gridInserted = true;
+            } else {
+                pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+                       "Grid %s refers to parent %s, but its extent is "
+                       "not included in it. Using bounding-box method.",
+                       gridName.c_str(), parentName.c_str());
+            }
+        }
+    } else if (!gridName.empty()) {
+        topGrids.emplace_back(std::move(grid));
+        gridInserted = true;
+    }
+
+    // Fallback to analyzing spatial extents
+    if (!gridInserted) {
+        for (const auto &candidateParent : topGrids) {
+            const auto &candidateParentExtent = candidateParent->extentAndRes();
+            if (candidateParentExtent.contains(extent)) {
+                static_cast<GridType *>(candidateParent.get())
+                    ->insertGrid(ctx, std::move(grid));
+                gridInserted = true;
+                break;
+            } else if (candidateParentExtent.intersects(extent)) {
+                pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+                       "Partially intersecting grids found!");
+            }
+        }
+        if (!gridInserted) {
+            topGrids.emplace_back(std::move(grid));
+        }
+    }
+}
+
+#ifdef TIFF_ENABLED
+// ---------------------------------------------------------------------------
+
+class GTiffVGrid : public VerticalShiftGrid {
+    friend void insertIntoHierarchy<GTiffVGrid, VerticalShiftGrid>(
+        PJ_CONTEXT *ctx, std::unique_ptr<GTiffVGrid> &&grid,
+        const std::string &gridName, const std::string &parentName,
+        std::vector<std::unique_ptr<VerticalShiftGrid>> &topGrids,
+        std::map<std::string, GTiffVGrid *> &mapGrids);
+
+    std::unique_ptr<GTiffGrid> m_grid;
+    uint16 m_idxSample;
+
+  public:
+    GTiffVGrid(std::unique_ptr<GTiffGrid> &&grid, uint16 idxSample);
+
+    ~GTiffVGrid() override;
+
+    bool valueAt(int x, int y, float &out) const override {
+        return m_grid->valueAt(m_idxSample, x, y, out);
+    }
+
+    bool isNodata(float val, double /* multiplier */) const override {
+        return m_grid->isNodata(val);
+    }
+
+    void insertGrid(PJ_CONTEXT *ctx, std::unique_ptr<GTiffVGrid> &&subgrid);
+
+    void reassign_context(PJ_CONTEXT *ctx) override {
+        m_grid->reassign_context(ctx);
+    }
+
+    bool hasChanged() const override { return m_grid->hasChanged(); }
+};
+
+// ---------------------------------------------------------------------------
+
+GTiffVGridShiftSet::~GTiffVGridShiftSet() = default;
+
+// ---------------------------------------------------------------------------
+
+GTiffVGrid::GTiffVGrid(std::unique_ptr<GTiffGrid> &&grid, uint16 idxSample)
+    : VerticalShiftGrid(grid->name(), grid->width(), grid->height(),
+                        grid->extentAndRes()),
+      m_grid(std::move(grid)), m_idxSample(idxSample) {}
+
+// ---------------------------------------------------------------------------
+
+GTiffVGrid::~GTiffVGrid() = default;
+
+// ---------------------------------------------------------------------------
+
+void GTiffVGrid::insertGrid(PJ_CONTEXT *ctx,
+                            std::unique_ptr<GTiffVGrid> &&subgrid) {
+    bool gridInserted = false;
+    const auto &extent = subgrid->extentAndRes();
+    for (const auto &candidateParent : m_children) {
+        const auto &candidateParentExtent = candidateParent->extentAndRes();
+        if (candidateParentExtent.contains(extent)) {
+            static_cast<GTiffVGrid *>(candidateParent.get())
+                ->insertGrid(ctx, std::move(subgrid));
+            gridInserted = true;
+            break;
+        } else if (candidateParentExtent.intersects(extent)) {
+            pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+                   "Partially intersecting grids found!");
+        }
+    }
+    if (!gridInserted) {
+        m_children.emplace_back(std::move(subgrid));
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+std::unique_ptr<GTiffVGridShiftSet>
+GTiffVGridShiftSet::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
+                         const std::string &filename) {
+    auto set = std::unique_ptr<GTiffVGridShiftSet>(
+        new GTiffVGridShiftSet(ctx, std::move(fp)));
+    set->m_name = filename;
+    set->m_format = "gtiff";
+    if (!set->m_GTiffDataset->openTIFF(filename)) {
+        return nullptr;
+    }
+    uint16 idxSample = 0;
+
+    std::map<std::string, GTiffVGrid *> mapGrids;
+    for (int ifd = 0;; ++ifd) {
+        auto grid = set->m_GTiffDataset->nextGrid();
+        if (!grid) {
+            if (ifd == 0) {
+                return nullptr;
+            }
+            break;
+        }
+
+        const auto subfileType = grid->subfileType();
+        if (subfileType != 0 && subfileType != FILETYPE_PAGE) {
+            if (ifd == 0) {
+                pj_log(ctx, PJ_LOG_ERROR, "Invalid subfileType");
+                return nullptr;
+            } else {
+                pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+                       "Ignoring IFD %d as it has a unsupported subfileType",
+                       ifd);
+                continue;
+            }
+        }
+
+        // Identify the index of the geoid_undulation/vertical_offset
+        bool foundDescriptionForAtLeastOneSample = false;
+        bool foundDescriptionForShift = false;
+        for (int i = 0; i < static_cast<int>(grid->samplesPerPixel()); ++i) {
+            const auto desc = grid->metadataItem("DESCRIPTION", i);
+            if (!desc.empty()) {
+                foundDescriptionForAtLeastOneSample = true;
+            }
+            if (desc == "geoid_undulation" || desc == "vertical_offset") {
+                idxSample = static_cast<uint16>(i);
+                foundDescriptionForShift = true;
+            }
+        }
+
+        if (foundDescriptionForAtLeastOneSample) {
+            if (!foundDescriptionForShift) {
+                if (ifd > 0) {
+                    // Assuming that extra IFD without our channel of interest
+                    // can be ignored
+                    // One could imagine to put the accuracy values in separate
+                    // IFD for example
+                    pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+                           "Ignoring IFD %d as it has no "
+                           "geoid_undulation/vertical_offset channel",
+                           ifd);
+                    continue;
+                } else {
+                    pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+                           "IFD 0 has channel descriptions, but no "
+                           "geoid_undulation/vertical_offset channel");
+                    return nullptr;
+                }
+            }
+        }
+
+        if (idxSample >= grid->samplesPerPixel()) {
+            pj_log(ctx, PJ_LOG_ERROR, "Invalid sample index");
+            return nullptr;
+        }
+
+        const std::string gridName = grid->metadataItem("grid_name");
+        const std::string parentName = grid->metadataItem("parent_grid_name");
+
+        auto vgrid =
+            internal::make_unique<GTiffVGrid>(std::move(grid), idxSample);
+
+        insertIntoHierarchy(ctx, std::move(vgrid), gridName, parentName,
+                            set->m_grids, mapGrids);
+    }
+    return set;
+}
+#endif // TIFF_ENABLED
+
+// ---------------------------------------------------------------------------
+
+std::unique_ptr<VerticalShiftGridSet>
+VerticalShiftGridSet::open(PJ_CONTEXT *ctx, const std::string &filename) {
+    if (filename == "null") {
+        auto set =
+            std::unique_ptr<VerticalShiftGridSet>(new VerticalShiftGridSet());
+        set->m_name = filename;
+        set->m_format = "null";
+        set->m_grids.push_back(std::unique_ptr<NullVerticalShiftGrid>(
+            new NullVerticalShiftGrid()));
+        return set;
+    }
+
+    auto fp = FileManager::open_resource_file(ctx, filename.c_str());
+    if (!fp) {
+        return nullptr;
+    }
+    const auto actualName(fp->name());
+    if (ends_with(actualName, "gtx") || ends_with(actualName, "GTX")) {
+        auto grid = GTXVerticalShiftGrid::open(ctx, std::move(fp), actualName);
+        if (!grid) {
+            return nullptr;
+        }
+        auto set =
+            std::unique_ptr<VerticalShiftGridSet>(new VerticalShiftGridSet());
+        set->m_name = actualName;
+        set->m_format = "gtx";
+        set->m_grids.push_back(std::unique_ptr<VerticalShiftGrid>(grid));
+        return set;
+    }
+
+    /* -------------------------------------------------------------------- */
+    /*      Load a header, to determine the file type.                      */
+    /* -------------------------------------------------------------------- */
+    unsigned char header[4];
+    size_t header_size = fp->read(header, sizeof(header));
+    if (header_size != sizeof(header)) {
+        return nullptr;
+    }
+    fp->seek(0);
+
+    if (IsTIFF(header_size, header)) {
+#ifdef TIFF_ENABLED
+        auto set = GTiffVGridShiftSet::open(ctx, std::move(fp), actualName);
+        if (!set)
+            pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return set;
+#else
+        pj_log(ctx, PJ_LOG_ERROR,
+               "TIFF grid, but TIFF support disabled in this build");
+        return nullptr;
+#endif
+    }
+
+    pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Unrecognized vertical grid format");
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+bool VerticalShiftGridSet::reopen(PJ_CONTEXT *ctx) {
+    pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Grid %s has changed. Re-loading it",
+           m_name.c_str());
+    auto newGS = open(ctx, m_name);
+    m_grids.clear();
+    if (newGS) {
+        m_grids = std::move(newGS->m_grids);
+    }
+    return !m_grids.empty();
+}
+
+// ---------------------------------------------------------------------------
+
+const VerticalShiftGrid *VerticalShiftGrid::gridAt(double lon,
+                                                   double lat) const {
+    for (const auto &child : m_children) {
+        const auto &extentChild = child->extentAndRes();
+        if ((extentChild.fullWorldLongitude() ||
+             (lon >= extentChild.westLon && lon <= extentChild.eastLon)) &&
+            lat >= extentChild.southLat && lat <= extentChild.northLat) {
+            return child->gridAt(lon, lat);
+        }
+    }
+    return this;
+}
+// ---------------------------------------------------------------------------
+
+const VerticalShiftGrid *VerticalShiftGridSet::gridAt(double lon,
+                                                      double lat) const {
+    for (const auto &grid : m_grids) {
+        if (dynamic_cast<NullVerticalShiftGrid *>(grid.get())) {
+            return grid.get();
+        }
+        const auto &extent = grid->extentAndRes();
+        if ((extent.fullWorldLongitude() ||
+             (lon >= extent.westLon && lon <= extent.eastLon)) &&
+            lat >= extent.southLat && lat <= extent.northLat) {
+            return grid->gridAt(lon, lat);
+        }
+    }
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+void VerticalShiftGridSet::reassign_context(PJ_CONTEXT *ctx) {
+    for (const auto &grid : m_grids) {
+        grid->reassign_context(ctx);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+HorizontalShiftGrid::HorizontalShiftGrid(const std::string &nameIn, int widthIn,
+                                         int heightIn,
+                                         const ExtentAndRes &extentIn)
+    : Grid(nameIn, widthIn, heightIn, extentIn) {}
+
+// ---------------------------------------------------------------------------
+
+HorizontalShiftGrid::~HorizontalShiftGrid() = default;
+
+// ---------------------------------------------------------------------------
+
+HorizontalShiftGridSet::HorizontalShiftGridSet() = default;
+
+// ---------------------------------------------------------------------------
+
+HorizontalShiftGridSet::~HorizontalShiftGridSet() = default;
+
+// ---------------------------------------------------------------------------
+
+class NullHorizontalShiftGrid : public HorizontalShiftGrid {
+
+  public:
+    NullHorizontalShiftGrid()
+        : HorizontalShiftGrid("null", 3, 3, globalExtent()) {}
+
+    bool isNullGrid() const override { return true; }
+
+    bool valueAt(int, int, bool, float &lonShift,
+                 float &latShift) const override;
+
+    void reassign_context(PJ_CONTEXT *) override {}
+
+    bool hasChanged() const override { return false; }
+};
+
+// ---------------------------------------------------------------------------
+
+bool NullHorizontalShiftGrid::valueAt(int, int, bool, 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;
+    std::unique_ptr<File> m_fp;
+
+    NTv1Grid(const NTv1Grid &) = delete;
+    NTv1Grid &operator=(const NTv1Grid &) = delete;
+
+  public:
+    explicit NTv1Grid(PJ_CONTEXT *ctx, std::unique_ptr<File> &&fp,
+                      const std::string &nameIn, int widthIn, int heightIn,
+                      const ExtentAndRes &extentIn)
+        : HorizontalShiftGrid(nameIn, widthIn, heightIn, extentIn), m_ctx(ctx),
+          m_fp(std::move(fp)) {}
+
+    ~NTv1Grid() override;
+
+    bool valueAt(int, int, bool, float &lonShift,
+                 float &latShift) const override;
+
+    static NTv1Grid *open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
+                          const std::string &filename);
+
+    void reassign_context(PJ_CONTEXT *ctx) override {
+        m_ctx = ctx;
+        m_fp->reassign_context(ctx);
+    }
+
+    bool hasChanged() const override { return m_fp->hasChanged(); }
+};
+
+// ---------------------------------------------------------------------------
+
+NTv1Grid::~NTv1Grid() = default;
+
+// ---------------------------------------------------------------------------
+
+NTv1Grid *NTv1Grid::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
+                         const std::string &filename) {
+    unsigned char header[192];
+
+    /* -------------------------------------------------------------------- */
+    /*      Read the header.                                                */
+    /* -------------------------------------------------------------------- */
+    if (fp->read(header, sizeof(header)) != sizeof(header)) {
+        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, std::move(fp), filename, columns, rows, extent);
+}
+
+// ---------------------------------------------------------------------------
+
+bool NTv1Grid::valueAt(int x, int y, bool compensateNTConvention,
+                       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 !
+    m_fp->seek(192 + 2 * sizeof(double) * (y * m_width + m_width - 1 - x));
+    if (m_fp->read(&two_doubles[0], sizeof(two_doubles)) !=
+        sizeof(two_doubles)) {
+        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));
+    // west longitude positive convention !
+    lonShift = (compensateNTConvention ? -1 : 1) *
+               static_cast<float>(two_doubles[1] * ((M_PI / 180.0) / 3600.0));
+
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+class CTable2Grid : public HorizontalShiftGrid {
+    PJ_CONTEXT *m_ctx;
+    std::unique_ptr<File> m_fp;
+
+    CTable2Grid(const CTable2Grid &) = delete;
+    CTable2Grid &operator=(const CTable2Grid &) = delete;
+
+  public:
+    CTable2Grid(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
+                const std::string &nameIn, int widthIn, int heightIn,
+                const ExtentAndRes &extentIn)
+        : HorizontalShiftGrid(nameIn, widthIn, heightIn, extentIn), m_ctx(ctx),
+          m_fp(std::move(fp)) {}
+
+    ~CTable2Grid() override;
+
+    bool valueAt(int, int, bool, float &lonShift,
+                 float &latShift) const override;
+
+    static CTable2Grid *open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
+                             const std::string &filename);
+
+    void reassign_context(PJ_CONTEXT *ctx) override {
+        m_ctx = ctx;
+        m_fp->reassign_context(ctx);
+    }
+
+    bool hasChanged() const override { return m_fp->hasChanged(); }
+};
+
+// ---------------------------------------------------------------------------
+
+CTable2Grid::~CTable2Grid() = default;
+
+// ---------------------------------------------------------------------------
+
+CTable2Grid *CTable2Grid::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
+                               const std::string &filename) {
+    unsigned char header[160];
+
+    /* -------------------------------------------------------------------- */
+    /*      Read the header.                                                */
+    /* -------------------------------------------------------------------- */
+    if (fp->read(header, sizeof(header)) != sizeof(header)) {
+        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, std::move(fp), filename, width, height, extent);
+}
+
+// ---------------------------------------------------------------------------
+
+bool CTable2Grid::valueAt(int x, int y, bool compensateNTConvention,
+                          float &lonShift, float &latShift) const {
+    assert(x >= 0 && y >= 0 && x < m_width && y < m_height);
+
+    float two_floats[2];
+    m_fp->seek(160 + 2 * sizeof(float) * (y * m_width + x));
+    if (m_fp->read(&two_floats[0], sizeof(two_floats)) != sizeof(two_floats)) {
+        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];
+    // west longitude positive convention !
+    lonShift = (compensateNTConvention ? -1 : 1) * two_floats[0];
+
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+class NTv2GridSet : public HorizontalShiftGridSet {
+    std::unique_ptr<File> m_fp;
+
+    NTv2GridSet(const NTv2GridSet &) = delete;
+    NTv2GridSet &operator=(const NTv2GridSet &) = delete;
+
+    explicit NTv2GridSet(std::unique_ptr<File> &&fp) : m_fp(std::move(fp)) {}
+
+  public:
+    ~NTv2GridSet() override;
+
+    static std::unique_ptr<NTv2GridSet> open(PJ_CONTEXT *ctx,
+                                             std::unique_ptr<File> fp,
+                                             const std::string &filename);
+
+    void reassign_context(PJ_CONTEXT *ctx) override {
+        HorizontalShiftGridSet::reassign_context(ctx);
+        m_fp->reassign_context(ctx);
+    }
+};
+
+// ---------------------------------------------------------------------------
+
+class NTv2Grid : public HorizontalShiftGrid {
+    friend class NTv2GridSet;
+
+    std::string m_name;
+    PJ_CONTEXT *m_ctx; // owned by the parent NTv2GridSet
+    File *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, File *fp,
+             unsigned long long offsetIn, bool mustSwapIn, int widthIn,
+             int heightIn, const ExtentAndRes &extentIn)
+        : HorizontalShiftGrid(nameIn, widthIn, heightIn, extentIn),
+          m_name(nameIn), m_ctx(ctx), m_fp(fp), m_offset(offsetIn),
+          m_mustSwap(mustSwapIn) {}
+
+    bool valueAt(int, int, bool, float &lonShift,
+                 float &latShift) const override;
+
+    void reassign_context(PJ_CONTEXT *ctx) override {
+        m_ctx = ctx;
+        m_fp->reassign_context(ctx);
+    }
+
+    bool hasChanged() const override { return m_fp->hasChanged(); }
+};
+
+// ---------------------------------------------------------------------------
+
+bool NTv2Grid::valueAt(int x, int y, bool compensateNTConvention,
+                       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
+    m_fp->seek(
+        m_offset +
+        4 * sizeof(float) *
+            (static_cast<unsigned long long>(y) * m_width + m_width - 1 - x));
+    if (m_fp->read(&two_float[0], sizeof(two_float)) != sizeof(two_float)) {
+        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));
+    // west longitude positive convention !
+    lonShift = (compensateNTConvention ? -1 : 1) *
+               static_cast<float>(two_float[1] * ((M_PI / 180.0) / 3600.0));
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+NTv2GridSet::~NTv2GridSet() = default;
+
+// ---------------------------------------------------------------------------
+
+std::unique_ptr<NTv2GridSet> NTv2GridSet::open(PJ_CONTEXT *ctx,
+                                               std::unique_ptr<File> fp,
+                                               const std::string &filename) {
+    File *fpRaw = fp.get();
+    auto set = std::unique_ptr<NTv2GridSet>(new NTv2GridSet(std::move(fp)));
+    set->m_name = filename;
+    set->m_format = "ntv2";
+
+    char header[11 * 16];
+
+    /* -------------------------------------------------------------------- */
+    /*      Read the header.                                                */
+    /* -------------------------------------------------------------------- */
+    if (fpRaw->read(header, sizeof(header)) != sizeof(header)) {
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return nullptr;
+    }
+
+    constexpr int OFFSET_GS_TYPE = 56;
+    if (memcmp(header + OFFSET_GS_TYPE, "SECONDS", 7) != 0) {
+        pj_log(ctx, PJ_LOG_ERROR, "Only GS_TYPE=SECONDS is supported");
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return nullptr;
+    }
+
+    const bool must_swap = (header[8] == 11) ? !IS_LSB : IS_LSB;
+    constexpr int OFFSET_NUM_SUBFILES = 8 + 32;
+    if (must_swap) {
+        // swap_words( header+8, 4, 1 );
+        // swap_words( header+8+16, 4, 1 );
+        swap_words(header + OFFSET_NUM_SUBFILES, 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 + OFFSET_NUM_SUBFILES, 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 (fpRaw->read(header, sizeof(header)) != sizeof(header)) {
+            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.
+        constexpr int OFFSET_GS_COUNT = 8 + 16 * 10;
+        constexpr int OFFSET_SOUTH_LAT = 8 + 16 * 4;
+        if (must_swap) {
+            // 6 double values: southLat, northLat, eastLon, westLon, resLat,
+            // resLon
+            swap_words(header + OFFSET_SOUTH_LAT, sizeof(double), 6);
+            swap_words(header + OFFSET_GS_COUNT, sizeof(int), 1);
+        }
+
+        std::string gridName;
+        gridName.append(header + 8, 8);
+
+        ExtentAndRes extent;
+        extent.southLat = to_double(header + OFFSET_SOUTH_LAT) * DEG_TO_RAD /
+                          3600.0; /* S_LAT */
+        extent.northLat = to_double(header + OFFSET_SOUTH_LAT + 16) *
+                          DEG_TO_RAD / 3600.0; /* N_LAT */
+        extent.eastLon = -to_double(header + OFFSET_SOUTH_LAT + 16 * 2) *
+                         DEG_TO_RAD / 3600.0; /* E_LONG */
+        extent.westLon = -to_double(header + OFFSET_SOUTH_LAT + 16 * 3) *
+                         DEG_TO_RAD / 3600.0; /* W_LONG */
+        extent.resLat =
+            to_double(header + OFFSET_SOUTH_LAT + 16 * 4) * DEG_TO_RAD / 3600.0;
+        extent.resLon =
+            to_double(header + OFFSET_SOUTH_LAT + 16 * 5) * 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 + OFFSET_GS_COUNT, 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;
+        }
+
+        const auto offset = fpRaw->tell();
+        auto grid = std::unique_ptr<NTv2Grid>(
+            new NTv2Grid(filename + ", " + gridName, ctx, fpRaw, 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
+        fpRaw->seek(static_cast<unsigned long long>(gs_count) * 4 * 4,
+                    SEEK_CUR);
+    }
+    return set;
+}
+
+#ifdef TIFF_ENABLED
+
+// ---------------------------------------------------------------------------
+
+class GTiffHGridShiftSet : public HorizontalShiftGridSet {
+
+    std::unique_ptr<GTiffDataset> m_GTiffDataset;
+
+    GTiffHGridShiftSet(PJ_CONTEXT *ctx, std::unique_ptr<File> &&fp)
+        : m_GTiffDataset(new GTiffDataset(ctx, std::move(fp))) {}
+
+  public:
+    ~GTiffHGridShiftSet() override;
+
+    static std::unique_ptr<GTiffHGridShiftSet>
+    open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
+         const std::string &filename);
+
+    void reassign_context(PJ_CONTEXT *ctx) override {
+        HorizontalShiftGridSet::reassign_context(ctx);
+        if (m_GTiffDataset) {
+            m_GTiffDataset->reassign_context(ctx);
+        }
+    }
+
+    bool reopen(PJ_CONTEXT *ctx) override {
+        pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Grid %s has changed. Re-loading it",
+               m_name.c_str());
+        m_grids.clear();
+        m_GTiffDataset.reset();
+        auto fp = FileManager::open_resource_file(ctx, m_name.c_str());
+        if (!fp) {
+            return false;
+        }
+        auto newGS = open(ctx, std::move(fp), m_name);
+        if (newGS) {
+            m_grids = std::move(newGS->m_grids);
+            m_GTiffDataset = std::move(newGS->m_GTiffDataset);
+        }
+        return !m_grids.empty();
+    }
+};
+
+// ---------------------------------------------------------------------------
+
+class GTiffHGrid : public HorizontalShiftGrid {
+    friend void insertIntoHierarchy<GTiffHGrid, HorizontalShiftGrid>(
+        PJ_CONTEXT *ctx, std::unique_ptr<GTiffHGrid> &&grid,
+        const std::string &gridName, const std::string &parentName,
+        std::vector<std::unique_ptr<HorizontalShiftGrid>> &topGrids,
+        std::map<std::string, GTiffHGrid *> &mapGrids);
+
+    std::unique_ptr<GTiffGrid> m_grid;
+    uint16 m_idxLatShift;
+    uint16 m_idxLonShift;
+    double m_convFactorToRadian;
+    bool m_positiveEast;
+
+  public:
+    GTiffHGrid(std::unique_ptr<GTiffGrid> &&grid, uint16 idxLatShift,
+               uint16 idxLonShift, double convFactorToRadian,
+               bool positiveEast);
+
+    ~GTiffHGrid() override;
+
+    bool valueAt(int x, int y, bool, float &lonShift,
+                 float &latShift) const override;
+
+    void insertGrid(PJ_CONTEXT *ctx, std::unique_ptr<GTiffHGrid> &&subgrid);
+
+    void reassign_context(PJ_CONTEXT *ctx) override {
+        m_grid->reassign_context(ctx);
+    }
+
+    bool hasChanged() const override { return m_grid->hasChanged(); }
+};
+
+// ---------------------------------------------------------------------------
+
+GTiffHGridShiftSet::~GTiffHGridShiftSet() = default;
+
+// ---------------------------------------------------------------------------
+
+GTiffHGrid::GTiffHGrid(std::unique_ptr<GTiffGrid> &&grid, uint16 idxLatShift,
+                       uint16 idxLonShift, double convFactorToRadian,
+                       bool positiveEast)
+    : HorizontalShiftGrid(grid->name(), grid->width(), grid->height(),
+                          grid->extentAndRes()),
+      m_grid(std::move(grid)), m_idxLatShift(idxLatShift),
+      m_idxLonShift(idxLonShift), m_convFactorToRadian(convFactorToRadian),
+      m_positiveEast(positiveEast) {}
+
+// ---------------------------------------------------------------------------
+
+GTiffHGrid::~GTiffHGrid() = default;
+
+// ---------------------------------------------------------------------------
+
+bool GTiffHGrid::valueAt(int x, int y, bool, float &lonShift,
+                         float &latShift) const {
+    if (!m_grid->valueAt(m_idxLatShift, x, y, latShift) ||
+        !m_grid->valueAt(m_idxLonShift, x, y, lonShift)) {
+        return false;
+    }
+    // From arc-seconds to radians
+    latShift = static_cast<float>(latShift * m_convFactorToRadian);
+    lonShift = static_cast<float>(lonShift * m_convFactorToRadian);
+    if (!m_positiveEast) {
+        lonShift = -lonShift;
+    }
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+void GTiffHGrid::insertGrid(PJ_CONTEXT *ctx,
+                            std::unique_ptr<GTiffHGrid> &&subgrid) {
+    bool gridInserted = false;
+    const auto &extent = subgrid->extentAndRes();
+    for (const auto &candidateParent : m_children) {
+        const auto &candidateParentExtent = candidateParent->extentAndRes();
+        if (candidateParentExtent.contains(extent)) {
+            static_cast<GTiffHGrid *>(candidateParent.get())
+                ->insertGrid(ctx, std::move(subgrid));
+            gridInserted = true;
+            break;
+        } else if (candidateParentExtent.intersects(extent)) {
+            pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+                   "Partially intersecting grids found!");
+        }
+    }
+    if (!gridInserted) {
+        m_children.emplace_back(std::move(subgrid));
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+std::unique_ptr<GTiffHGridShiftSet>
+GTiffHGridShiftSet::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
+                         const std::string &filename) {
+    auto set = std::unique_ptr<GTiffHGridShiftSet>(
+        new GTiffHGridShiftSet(ctx, std::move(fp)));
+    set->m_name = filename;
+    set->m_format = "gtiff";
+    if (!set->m_GTiffDataset->openTIFF(filename)) {
+        return nullptr;
+    }
+
+    // Defaults inspired from NTv2
+    uint16 idxLatShift = 0;
+    uint16 idxLonShift = 1;
+    constexpr double ARC_SECOND_TO_RADIAN = (M_PI / 180.0) / 3600.0;
+    double convFactorToRadian = ARC_SECOND_TO_RADIAN;
+    bool positiveEast = true;
+
+    std::map<std::string, GTiffHGrid *> mapGrids;
+    for (int ifd = 0;; ++ifd) {
+        auto grid = set->m_GTiffDataset->nextGrid();
+        if (!grid) {
+            if (ifd == 0) {
+                return nullptr;
+            }
+            break;
+        }
+
+        const auto subfileType = grid->subfileType();
+        if (subfileType != 0 && subfileType != FILETYPE_PAGE) {
+            if (ifd == 0) {
+                pj_log(ctx, PJ_LOG_ERROR, "Invalid subfileType");
+                return nullptr;
+            } else {
+                pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+                       "Ignoring IFD %d as it has a unsupported subfileType",
+                       ifd);
+                continue;
+            }
+        }
+
+        if (grid->samplesPerPixel() < 2) {
+            if (ifd == 0) {
+                pj_log(ctx, PJ_LOG_ERROR,
+                       "At least 2 samples per pixel needed");
+                return nullptr;
+            } else {
+                pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+                       "Ignoring IFD %d as it has not at least 2 samples", ifd);
+                continue;
+            }
+        }
+
+        // Identify the index of the latitude and longitude offset channels
+        bool foundDescriptionForAtLeastOneSample = false;
+        bool foundDescriptionForLatOffset = false;
+        bool foundDescriptionForLonOffset = false;
+        for (int i = 0; i < static_cast<int>(grid->samplesPerPixel()); ++i) {
+            const auto desc = grid->metadataItem("DESCRIPTION", i);
+            if (!desc.empty()) {
+                foundDescriptionForAtLeastOneSample = true;
+            }
+            if (desc == "latitude_offset") {
+                idxLatShift = static_cast<uint16>(i);
+                foundDescriptionForLatOffset = true;
+            } else if (desc == "longitude_offset") {
+                idxLonShift = static_cast<uint16>(i);
+                foundDescriptionForLonOffset = true;
+            }
+        }
+
+        if (foundDescriptionForAtLeastOneSample) {
+            if (!foundDescriptionForLonOffset &&
+                !foundDescriptionForLatOffset) {
+                if (ifd > 0) {
+                    // Assuming that extra IFD without
+                    // longitude_offset/latitude_offset can be ignored
+                    // One could imagine to put the accuracy values in separate
+                    // IFD for example
+                    pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+                           "Ignoring IFD %d as it has no "
+                           "longitude_offset/latitude_offset channel",
+                           ifd);
+                    continue;
+                } else {
+                    pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+                           "IFD 0 has channel descriptions, but no "
+                           "longitude_offset/latitude_offset channel");
+                    return nullptr;
+                }
+            }
+        }
+        if (foundDescriptionForLatOffset && !foundDescriptionForLonOffset) {
+            pj_log(ctx, PJ_LOG_ERROR,
+                   "Found latitude_offset channel, but not longitude_offset");
+            return nullptr;
+        } else if (foundDescriptionForLonOffset &&
+                   !foundDescriptionForLatOffset) {
+            pj_log(ctx, PJ_LOG_ERROR,
+                   "Found longitude_offset channel, but not latitude_offset");
+            return nullptr;
+        }
+
+        if (idxLatShift >= grid->samplesPerPixel() ||
+            idxLonShift >= grid->samplesPerPixel()) {
+            pj_log(ctx, PJ_LOG_ERROR, "Invalid sample index");
+            return nullptr;
+        }
+
+        if (foundDescriptionForLonOffset) {
+            const std::string positiveValue =
+                grid->metadataItem("positive_value", idxLonShift);
+            if (!positiveValue.empty()) {
+                if (positiveValue == "west") {
+                    positiveEast = false;
+                } else if (positiveValue == "east") {
+                    positiveEast = true;
+                } else {
+                    pj_log(ctx, PJ_LOG_ERROR,
+                           "Unsupported value %s for 'positive_value'",
+                           positiveValue.c_str());
+                    return nullptr;
+                }
+            }
+        }
+
+        // Identify their unit
+        {
+            const auto unitLatShift =
+                grid->metadataItem("UNITTYPE", idxLatShift);
+            const auto unitLonShift =
+                grid->metadataItem("UNITTYPE", idxLonShift);
+            if (unitLatShift != unitLonShift) {
+                pj_log(ctx, PJ_LOG_ERROR,
+                       "Different unit for longitude and latitude offset");
+                return nullptr;
+            }
+            if (!unitLatShift.empty()) {
+                if (unitLatShift == "arc-second") {
+                    convFactorToRadian = ARC_SECOND_TO_RADIAN;
+                } else if (unitLatShift == "radian") {
+                    convFactorToRadian = 1.0;
+                } else if (unitLatShift == "degree") {
+                    convFactorToRadian = M_PI / 180.0;
+                } else {
+                    pj_log(ctx, PJ_LOG_ERROR, "Unsupported unit %s",
+                           unitLatShift.c_str());
+                    return nullptr;
+                }
+            }
+        }
+
+        const std::string gridName = grid->metadataItem("grid_name");
+        const std::string parentName = grid->metadataItem("parent_grid_name");
+
+        auto hgrid = internal::make_unique<GTiffHGrid>(
+            std::move(grid), idxLatShift, idxLonShift, convFactorToRadian,
+            positiveEast);
+
+        insertIntoHierarchy(ctx, std::move(hgrid), gridName, parentName,
+                            set->m_grids, mapGrids);
+    }
+    return set;
+}
+#endif // TIFF_ENABLED
+
+// ---------------------------------------------------------------------------
+
+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;
+    }
+
+    auto fp = FileManager::open_resource_file(ctx, filename.c_str());
+    if (!fp) {
+        return nullptr;
+    }
+    const auto actualName(fp->name());
+
+    char header[160];
+    /* -------------------------------------------------------------------- */
+    /*      Load a header, to determine the file type.                      */
+    /* -------------------------------------------------------------------- */
+    size_t header_size = fp->read(header, sizeof(header));
+    if (header_size != 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);
+    }
+    fp->seek(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, std::move(fp), actualName);
+        if (!grid) {
+            return nullptr;
+        }
+        auto set = std::unique_ptr<HorizontalShiftGridSet>(
+            new HorizontalShiftGridSet());
+        set->m_name = actualName;
+        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, std::move(fp), actualName);
+        if (!grid) {
+            return nullptr;
+        }
+        auto set = std::unique_ptr<HorizontalShiftGridSet>(
+            new HorizontalShiftGridSet());
+        set->m_name = actualName;
+        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, std::move(fp), actualName);
+    } else if (IsTIFF(header_size,
+                      reinterpret_cast<const unsigned char *>(header))) {
+#ifdef TIFF_ENABLED
+        auto set = GTiffHGridShiftSet::open(ctx, std::move(fp), actualName);
+        if (!set)
+            pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return set;
+#else
+        pj_log(ctx, PJ_LOG_ERROR,
+               "TIFF grid, but TIFF support disabled in this build");
+        return nullptr;
+#endif
+    }
+
+    pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Unrecognized horizontal grid format");
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+bool HorizontalShiftGridSet::reopen(PJ_CONTEXT *ctx) {
+    pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Grid %s has changed. Re-loading it",
+           m_name.c_str());
+    auto newGS = open(ctx, m_name);
+    m_grids.clear();
+    if (newGS) {
+        m_grids = std::move(newGS->m_grids);
+    }
+    return !m_grids.empty();
+}
+
+// ---------------------------------------------------------------------------
+
+#define REL_TOLERANCE_HGRIDSHIFT 1e-5
+
+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) *
+                               REL_TOLERANCE_HGRIDSHIFT;
+        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) * REL_TOLERANCE_HGRIDSHIFT;
+        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;
+}
+
+// ---------------------------------------------------------------------------
+
+void HorizontalShiftGridSet::reassign_context(PJ_CONTEXT *ctx) {
+    for (const auto &grid : m_grids) {
+        grid->reassign_context(ctx);
+    }
+}
+
+#ifdef TIFF_ENABLED
+// ---------------------------------------------------------------------------
+
+class GTiffGenericGridShiftSet : public GenericShiftGridSet {
+
+    std::unique_ptr<GTiffDataset> m_GTiffDataset;
+
+    GTiffGenericGridShiftSet(PJ_CONTEXT *ctx, std::unique_ptr<File> &&fp)
+        : m_GTiffDataset(new GTiffDataset(ctx, std::move(fp))) {}
+
+  public:
+    ~GTiffGenericGridShiftSet() override;
+
+    static std::unique_ptr<GTiffGenericGridShiftSet>
+    open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
+         const std::string &filename);
+
+    void reassign_context(PJ_CONTEXT *ctx) override {
+        GenericShiftGridSet::reassign_context(ctx);
+        if (m_GTiffDataset) {
+            m_GTiffDataset->reassign_context(ctx);
+        }
+    }
+
+    bool reopen(PJ_CONTEXT *ctx) override {
+        pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Grid %s has changed. Re-loading it",
+               m_name.c_str());
+        m_grids.clear();
+        m_GTiffDataset.reset();
+        auto fp = FileManager::open_resource_file(ctx, m_name.c_str());
+        if (!fp) {
+            return false;
+        }
+        auto newGS = open(ctx, std::move(fp), m_name);
+        if (newGS) {
+            m_grids = std::move(newGS->m_grids);
+            m_GTiffDataset = std::move(newGS->m_GTiffDataset);
+        }
+        return !m_grids.empty();
+    }
+};
+
+// ---------------------------------------------------------------------------
+
+class GTiffGenericGrid : public GenericShiftGrid {
+    friend void insertIntoHierarchy<GTiffGenericGrid, GenericShiftGrid>(
+        PJ_CONTEXT *ctx, std::unique_ptr<GTiffGenericGrid> &&grid,
+        const std::string &gridName, const std::string &parentName,
+        std::vector<std::unique_ptr<GenericShiftGrid>> &topGrids,
+        std::map<std::string, GTiffGenericGrid *> &mapGrids);
+
+    std::unique_ptr<GTiffGrid> m_grid;
+
+  public:
+    GTiffGenericGrid(std::unique_ptr<GTiffGrid> &&grid);
+
+    ~GTiffGenericGrid() override;
+
+    bool valueAt(int x, int y, int sample, float &out) const override;
+
+    int samplesPerPixel() const override { return m_grid->samplesPerPixel(); }
+
+    std::string unit(int sample) const override {
+        return m_grid->metadataItem("UNITTYPE", sample);
+    }
+
+    std::string description(int sample) const override {
+        return m_grid->metadataItem("DESCRIPTION", sample);
+    }
+
+    std::string metadataItem(const std::string &key,
+                             int sample = -1) const override {
+        return m_grid->metadataItem(key, sample);
+    }
+
+    void insertGrid(PJ_CONTEXT *ctx,
+                    std::unique_ptr<GTiffGenericGrid> &&subgrid);
+
+    void reassign_context(PJ_CONTEXT *ctx) override {
+        m_grid->reassign_context(ctx);
+    }
+
+    bool hasChanged() const override { return m_grid->hasChanged(); }
+};
+
+// ---------------------------------------------------------------------------
+
+GTiffGenericGridShiftSet::~GTiffGenericGridShiftSet() = default;
+
+// ---------------------------------------------------------------------------
+
+GTiffGenericGrid::GTiffGenericGrid(std::unique_ptr<GTiffGrid> &&grid)
+    : GenericShiftGrid(grid->name(), grid->width(), grid->height(),
+                       grid->extentAndRes()),
+      m_grid(std::move(grid)) {}
+
+// ---------------------------------------------------------------------------
+
+GTiffGenericGrid::~GTiffGenericGrid() = default;
+
+// ---------------------------------------------------------------------------
+
+bool GTiffGenericGrid::valueAt(int x, int y, int sample, float &out) const {
+    if (sample < 0 ||
+        static_cast<unsigned>(sample) >= m_grid->samplesPerPixel())
+        return false;
+    return m_grid->valueAt(static_cast<uint16>(sample), x, y, out);
+}
+
+// ---------------------------------------------------------------------------
+
+void GTiffGenericGrid::insertGrid(PJ_CONTEXT *ctx,
+                                  std::unique_ptr<GTiffGenericGrid> &&subgrid) {
+    bool gridInserted = false;
+    const auto &extent = subgrid->extentAndRes();
+    for (const auto &candidateParent : m_children) {
+        const auto &candidateParentExtent = candidateParent->extentAndRes();
+        if (candidateParentExtent.contains(extent)) {
+            static_cast<GTiffGenericGrid *>(candidateParent.get())
+                ->insertGrid(ctx, std::move(subgrid));
+            gridInserted = true;
+            break;
+        } else if (candidateParentExtent.intersects(extent)) {
+            pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+                   "Partially intersecting grids found!");
+        }
+    }
+    if (!gridInserted) {
+        m_children.emplace_back(std::move(subgrid));
+    }
+}
+#endif // TIFF_ENABLED
+
+// ---------------------------------------------------------------------------
+
+class NullGenericShiftGrid : public GenericShiftGrid {
+
+  public:
+    NullGenericShiftGrid() : GenericShiftGrid("null", 3, 3, globalExtent()) {}
+
+    bool isNullGrid() const override { return true; }
+    bool valueAt(int, int, int, float &out) const override;
+
+    int samplesPerPixel() const override { return 0; }
+
+    std::string unit(int) const override { return std::string(); }
+
+    std::string description(int) const override { return std::string(); }
+
+    std::string metadataItem(const std::string &, int) const override {
+        return std::string();
+    }
+
+    void reassign_context(PJ_CONTEXT *) override {}
+
+    bool hasChanged() const override { return false; }
+};
+
+// ---------------------------------------------------------------------------
+
+bool NullGenericShiftGrid::valueAt(int, int, int, float &out) const {
+    out = 0.0f;
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+#ifdef TIFF_ENABLED
+
+std::unique_ptr<GTiffGenericGridShiftSet>
+GTiffGenericGridShiftSet::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
+                               const std::string &filename) {
+    auto set = std::unique_ptr<GTiffGenericGridShiftSet>(
+        new GTiffGenericGridShiftSet(ctx, std::move(fp)));
+    set->m_name = filename;
+    set->m_format = "gtiff";
+    if (!set->m_GTiffDataset->openTIFF(filename)) {
+        return nullptr;
+    }
+
+    std::map<std::string, GTiffGenericGrid *> mapGrids;
+    for (int ifd = 0;; ++ifd) {
+        auto grid = set->m_GTiffDataset->nextGrid();
+        if (!grid) {
+            if (ifd == 0) {
+                return nullptr;
+            }
+            break;
+        }
+
+        const auto subfileType = grid->subfileType();
+        if (subfileType != 0 && subfileType != FILETYPE_PAGE) {
+            if (ifd == 0) {
+                pj_log(ctx, PJ_LOG_ERROR, "Invalid subfileType");
+                return nullptr;
+            } else {
+                pj_log(ctx, PJ_LOG_DEBUG_MAJOR,
+                       "Ignoring IFD %d as it has a unsupported subfileType",
+                       ifd);
+                continue;
+            }
+        }
+
+        const std::string gridName = grid->metadataItem("grid_name");
+        const std::string parentName = grid->metadataItem("parent_grid_name");
+
+        auto hgrid = internal::make_unique<GTiffGenericGrid>(std::move(grid));
+
+        insertIntoHierarchy(ctx, std::move(hgrid), gridName, parentName,
+                            set->m_grids, mapGrids);
+    }
+    return set;
+}
+#endif // TIFF_ENABLED
+
+// ---------------------------------------------------------------------------
+
+GenericShiftGrid::GenericShiftGrid(const std::string &nameIn, int widthIn,
+                                   int heightIn, const ExtentAndRes &extentIn)
+    : Grid(nameIn, widthIn, heightIn, extentIn) {}
+
+// ---------------------------------------------------------------------------
+
+GenericShiftGrid::~GenericShiftGrid() = default;
+
+// ---------------------------------------------------------------------------
+
+GenericShiftGridSet::GenericShiftGridSet() = default;
+
+// ---------------------------------------------------------------------------
+
+GenericShiftGridSet::~GenericShiftGridSet() = default;
+
+// ---------------------------------------------------------------------------
+
+std::unique_ptr<GenericShiftGridSet>
+GenericShiftGridSet::open(PJ_CONTEXT *ctx, const std::string &filename) {
+    if (filename == "null") {
+        auto set =
+            std::unique_ptr<GenericShiftGridSet>(new GenericShiftGridSet());
+        set->m_name = filename;
+        set->m_format = "null";
+        set->m_grids.push_back(
+            std::unique_ptr<NullGenericShiftGrid>(new NullGenericShiftGrid()));
+        return set;
+    }
+
+    auto fp = FileManager::open_resource_file(ctx, filename.c_str());
+    if (!fp) {
+        return nullptr;
+    }
+    const auto actualName(fp->name());
+
+    /* -------------------------------------------------------------------- */
+    /*      Load a header, to determine the file type.                      */
+    /* -------------------------------------------------------------------- */
+    unsigned char header[4];
+    size_t header_size = fp->read(header, sizeof(header));
+    if (header_size != sizeof(header)) {
+        return nullptr;
+    }
+    fp->seek(0);
+
+    if (IsTIFF(header_size, header)) {
+#ifdef TIFF_ENABLED
+        auto set =
+            GTiffGenericGridShiftSet::open(ctx, std::move(fp), actualName);
+        if (!set)
+            pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return set;
+#else
+        pj_log(ctx, PJ_LOG_ERROR,
+               "TIFF grid, but TIFF support disabled in this build");
+        return nullptr;
+#endif
+    }
+
+    pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Unrecognized generic grid format");
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+bool GenericShiftGridSet::reopen(PJ_CONTEXT *ctx) {
+    pj_log(ctx, PJ_LOG_DEBUG_MAJOR, "Grid %s has changed. Re-loading it",
+           m_name.c_str());
+    auto newGS = open(ctx, m_name);
+    m_grids.clear();
+    if (newGS) {
+        m_grids = std::move(newGS->m_grids);
+    }
+    return !m_grids.empty();
+}
+
+// ---------------------------------------------------------------------------
+
+const GenericShiftGrid *GenericShiftGrid::gridAt(double lon, double lat) const {
+    for (const auto &child : m_children) {
+        const auto &extentChild = child->extentAndRes();
+        if ((extentChild.fullWorldLongitude() ||
+             (lon >= extentChild.westLon && lon <= extentChild.eastLon)) &&
+            lat >= extentChild.southLat && lat <= extentChild.northLat) {
+            return child->gridAt(lon, lat);
+        }
+    }
+    return this;
+}
+
+// ---------------------------------------------------------------------------
+
+const GenericShiftGrid *GenericShiftGridSet::gridAt(double lon,
+                                                    double lat) const {
+    for (const auto &grid : m_grids) {
+        if (dynamic_cast<NullGenericShiftGrid *>(grid.get())) {
+            return grid.get();
+        }
+        const auto &extent = grid->extentAndRes();
+        if ((extent.fullWorldLongitude() ||
+             (lon >= extent.westLon && lon <= extent.eastLon)) &&
+            lat >= extent.southLat && lat <= extent.northLat) {
+            return grid->gridAt(lon, lat);
+        }
+    }
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+void GenericShiftGridSet::reassign_context(PJ_CONTEXT *ctx) {
+    for (const auto &grid : m_grids) {
+        grid->reassign_context(ctx);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+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;
+    if (gridnames == nullptr)
+        return {};
+
+    auto listOfGridNames = internal::split(std::string(gridnames), ',');
+    ListOfGenericGrids grids;
+    for (const auto &gridnameStr : listOfGridNames) {
+        const char *gridname = gridnameStr.c_str();
+        bool canFail = false;
+        if (gridname[0] == '@') {
+            canFail = true;
+            gridname++;
+        }
+        auto gridSet = GenericShiftGridSet::open(P->ctx, gridname);
+        if (!gridSet) {
+            if (!canFail) {
+                if (proj_context_errno(P->ctx) != PJD_ERR_NETWORK_ERROR) {
+                    pj_ctx_set_errno(P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+                }
+                return {};
+            }
+            pj_ctx_set_errno(P->ctx, 0); // don't treat as a persistent error
+        } else {
+            grids.emplace_back(std::move(gridSet));
+        }
+    }
+
+    return grids;
+}
+
+// ---------------------------------------------------------------------------
+
+static const HorizontalShiftGrid *
+findGrid(const ListOfHGrids &grids, const PJ_LP &input,
+         HorizontalShiftGridSet *&gridSetOut) {
+    for (const auto &gridset : grids) {
+        auto grid = gridset->gridAt(input.lam, input.phi);
+        if (grid) {
+            gridSetOut = gridset.get();
+            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) {
+                if (proj_context_errno(ctx) != PJD_ERR_NETWORK_ERROR) {
+                    pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+                }
+                return {};
+            }
+            pj_ctx_set_errno(ctx, 0); // don't treat as a persistent error
+        } 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;
+
+// Apply bilinear interpolation for horizontal shift grids
+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 > 1 - 10 * REL_TOLERANCE_HGRIDSHIFT) {
+            ++indx.lam;
+            frct.lam = 0.;
+        } else
+            return val;
+    } else if ((in = indx.lam + 1) >= grid->width()) {
+        if (in == grid->width() && frct.lam < 10 * REL_TOLERANCE_HGRIDSHIFT) {
+            --indx.lam;
+            frct.lam = 1.;
+        } else
+            return val;
+    }
+    if (indx.phi < 0) {
+        if (indx.phi == -1 && frct.phi > 1 - 10 * REL_TOLERANCE_HGRIDSHIFT) {
+            ++indx.phi;
+            frct.phi = 0.;
+        } else
+            return val;
+    } else if ((in = indx.phi + 1) >= grid->height()) {
+        if (in == grid->height() && frct.phi < 10 * REL_TOLERANCE_HGRIDSHIFT) {
+            --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,
+                                     HorizontalShiftGridSet *gridset,
+                                     const ListOfHGrids &grids,
+                                     bool &shouldRetry) {
+    PJ_LP t, tb, del, dif;
+    int i = MAX_ITERATIONS;
+    const double toltol = TOL * TOL;
+
+    shouldRetry = false;
+    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;
+
+    t = pj_hgrid_interpolate(tb, grid, true);
+    if (grid->hasChanged()) {
+        shouldRetry = gridset->reopen(ctx);
+        return t;
+    }
+    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);
+        if (grid->hasChanged()) {
+            shouldRetry = gridset->reopen(ctx);
+            return t;
+        }
+
+        /* 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, gridset);
+            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;
+            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;
+
+    while (true) {
+        HorizontalShiftGridSet *gridset = nullptr;
+        const auto grid = findGrid(grids, lp, gridset);
+        if (!grid) {
+            pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+            return out;
+        }
+        if (grid->isNullGrid()) {
+            return lp;
+        }
+
+        bool shouldRetry = false;
+        out = pj_hgrid_apply_internal(ctx, lp, direction, grid, gridset, grids,
+                                      shouldRetry);
+        if (!shouldRetry) {
+            break;
+        }
+    }
+
+    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;
+
+    HorizontalShiftGridSet *gridset = nullptr;
+    const auto grid = findGrid(grids, lp, gridset);
+    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 (grid->hasChanged()) {
+        if (gridset->reopen(P->ctx)) {
+            return pj_hgrid_value(P, grids, lp);
+        }
+        out.lam = HUGE_VAL;
+        out.phi = HUGE_VAL;
+    }
+
+    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(PJ_CONTEXT *ctx, const ListOfVGrids &grids,
+                               const PJ_LP &input, const double vmultiplier) {
+
+    /* do not deal with NaN coordinates */
+    /* cppcheck-suppress duplicateExpression */
+    if (std::isnan(input.phi) || std::isnan(input.lam)) {
+        return HUGE_VAL;
+    }
+
+    VerticalShiftGridSet *curGridset = nullptr;
+    const VerticalShiftGrid *grid = nullptr;
+    for (const auto &gridset : grids) {
+        grid = gridset->gridAt(input.lam, input.phi);
+        if (grid) {
+            curGridset = gridset.get();
+            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;
+    bool error = (!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));
+    if (grid->hasChanged()) {
+        if (curGridset->reopen(ctx)) {
+            return read_vgrid_value(ctx, grids, input, vmultiplier);
+        }
+        error = true;
+    }
+
+    if (error) {
+        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) {
+                if (proj_context_errno(P->ctx) != PJD_ERR_NETWORK_ERROR) {
+                    pj_ctx_set_errno(P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+                }
+                return {};
+            }
+            pj_ctx_set_errno(P->ctx, 0); // don't treat as a persistent error
+        } 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(P->ctx, 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;
+}
+
+// ---------------------------------------------------------------------------
+
+const GenericShiftGrid *pj_find_generic_grid(const ListOfGenericGrids &grids,
+                                             const PJ_LP &input,
+                                             GenericShiftGridSet *&gridSetOut) {
+    for (const auto &gridset : grids) {
+        auto grid = gridset->gridAt(input.lam, input.phi);
+        if (grid) {
+            gridSetOut = gridset.get();
+            return grid;
+        }
+    }
+    return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+// Used by +proj=deformation and +proj=xyzgridshift to do bilinear interpolation
+// on 3 sample values per node.
+bool pj_bilinear_interpolation_three_samples(const GenericShiftGrid *grid,
+                                             const PJ_LP &lp, int idx1,
+                                             int idx2, int idx3, double &v1,
+                                             double &v2, double &v3,
+                                             bool &must_retry) {
+    must_retry = false;
+
+    const auto &extent = grid->extentAndRes();
+    double grid_x = (lp.lam - extent.westLon) / extent.resLon;
+    double grid_y = (lp.phi - extent.southLat) / extent.resLat;
+    int ix = static_cast<int>(grid_x);
+    int iy = static_cast<int>(grid_y);
+    int ix2 = std::min(ix + 1, grid->width() - 1);
+    int iy2 = std::min(iy + 1, grid->height() - 1);
+
+    float dx1 = 0.0f, dy1 = 0.0f, dz1 = 0.0f;
+    float dx2 = 0.0f, dy2 = 0.0f, dz2 = 0.0f;
+    float dx3 = 0.0f, dy3 = 0.0f, dz3 = 0.0f;
+    float dx4 = 0.0f, dy4 = 0.0f, dz4 = 0.0f;
+    bool error = (!grid->valueAt(ix, iy, idx1, dx1) ||
+                  !grid->valueAt(ix, iy, idx2, dy1) ||
+                  !grid->valueAt(ix, iy, idx3, dz1) ||
+                  !grid->valueAt(ix2, iy, idx1, dx2) ||
+                  !grid->valueAt(ix2, iy, idx2, dy2) ||
+                  !grid->valueAt(ix2, iy, idx3, dz2) ||
+                  !grid->valueAt(ix, iy2, idx1, dx3) ||
+                  !grid->valueAt(ix, iy2, idx2, dy3) ||
+                  !grid->valueAt(ix, iy2, idx3, dz3) ||
+                  !grid->valueAt(ix2, iy2, idx1, dx4) ||
+                  !grid->valueAt(ix2, iy2, idx2, dy4) ||
+                  !grid->valueAt(ix2, iy2, idx3, dz4));
+    if (grid->hasChanged()) {
+        must_retry = true;
+        return false;
+    }
+    if (error) {
+        return false;
+    }
+
+    double frct_lam = grid_x - ix;
+    double frct_phi = grid_y - iy;
+    double m10 = frct_lam;
+    double m11 = m10;
+    double m01 = 1. - frct_lam;
+    double m00 = m01;
+    m11 *= frct_phi;
+    m01 *= frct_phi;
+    frct_phi = 1. - frct_phi;
+    m00 *= frct_phi;
+    m10 *= frct_phi;
+
+    v1 = m00 * dx1 + m10 * dx2 + m01 * dx3 + m11 * dx4;
+    v2 = m00 * dy1 + m10 * dy2 + m01 * dy3 + m11 * dy4;
+    v3 = m00 * dz1 + m10 * dz2 + m01 * dz3 + m11 * dz4;
+    return true;
+}
+
+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;
+}