/******************************************************************************
 * 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();
}

// ---------------------------------------------------------------------------

static bool isPointInExtent(double lon, double lat, const ExtentAndRes &extent,
                            double eps = 0) {
    if (!(lat + eps >= extent.southLat && lat - eps <= extent.northLat))
        return false;
    if (extent.fullWorldLongitude())
        return true;
    if (lon + eps < extent.westLon)
        lon += 2 * M_PI;
    else if (lon - eps > extent.eastLon)
        lon -= 2 * M_PI;
    if (!(lon + eps >= extent.westLon && lon - eps <= extent.eastLon))
        return false;
    return true;
}

// ---------------------------------------------------------------------------

const VerticalShiftGrid *VerticalShiftGrid::gridAt(double lon,
                                                   double lat) const {
    for (const auto &child : m_children) {
        const auto &extentChild = child->extentAndRes();
        if (isPointInExtent(lon, lat, extentChild)) {
            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 (isPointInExtent(lon, lat, extent)) {
            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 (isPointInExtent(lon, lat, extentChild, epsilon)) {
            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 (isPointInExtent(lon, lat, extent, epsilon)) {
            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 (isPointInExtent(lon, lat, extentChild)) {
            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 (isPointInExtent(lon, lat, extent)) {
            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());
    const double epsilon =
        (extent->resLon + extent->resLat) * REL_TOLERANCE_HGRIDSHIFT;
    tb.lam -= extent->westLon;
    if (tb.lam + epsilon < 0)
        tb.lam += 2 * M_PI;
    else if (tb.lam - epsilon > extent->eastLon - extent->westLon)
        tb.lam -= 2 * M_PI;
    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;
            if (tb.lam + epsilon < 0)
                tb.lam += 2 * M_PI;
            else if (tb.lam - epsilon > extent->eastLon - extent->westLon)
                tb.lam -= 2 * M_PI;
            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();
    const double epsilon =
        (extent.resLon + extent.resLat) * REL_TOLERANCE_HGRIDSHIFT;
    lp.lam -= extent.westLon;
    if (lp.lam + epsilon < 0)
        lp.lam += 2 * M_PI;
    else if (lp.lam - epsilon > extent.eastLon - extent.westLon)
        lp.lam -= 2 * M_PI;
    lp.phi -= extent.southLat;

    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) {
        pj_ctx_set_errno(ctx, PJD_ERR_GRID_AREA);
        return HUGE_VAL;
    }
    if (grid->isNullGrid()) {
        return 0;
    }

    const auto &extent = grid->extentAndRes();

    /* Interpolation of a location within the grid */
    double grid_x = (input.lam - extent.westLon) / extent.resLon;
    if (input.lam < extent.westLon) {
        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());
        } else {
            grid_x = (input.lam + 2 * M_PI - extent.westLon) / extent.resLon;
        }
    } else if (input.lam > extent.eastLon) {
        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());
        } else {
            grid_x = (input.lam - 2 * M_PI - extent.westLon) / extent.resLon;
        }
    }
    double grid_y = (input.phi - extent.southLat) / extent.resLat;
    int grid_ix = static_cast<int>(lround(floor(grid_x)));
    if (!(grid_ix >= 0 && grid_ix < grid->width())) {
        // in the unlikely case we end up here...
        pj_log(ctx, PJ_LOG_ERROR, "grid_ix not in grid");
        pj_ctx_set_errno(ctx, PJD_ERR_GRID_AREA);
        return HUGE_VAL;
    }
    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;
    if (grid->isNullGrid()) {
        v1 = 0.0;
        v2 = 0.0;
        v3 = 0.0;
        return true;
    }

    const auto &extent = grid->extentAndRes();
    double grid_x = (lp.lam - extent.westLon) / extent.resLon;
    if (lp.lam < extent.westLon) {
        grid_x = (lp.lam + 2 * M_PI - extent.westLon) / extent.resLon;
    } else if (lp.lam > extent.eastLon) {
        grid_x = (lp.lam - 2 * M_PI - 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;
}
