Method grid_interpolate

src/transformations/gridshift.cpp:241–580

#define REL_TOLERANCE_HGRIDSHIFT 1e-5

PJ_XYZ gridshiftData::grid_interpolate(PJ_CONTEXT *ctx, const std::string &type,
                                       PJ_XY xy, const GenericShiftGrid *grid,
                                       bool &biquadraticInterpolationOut) {
    PJ_XYZ val;

    val.x = val.y = HUGE_VAL;
    val.z = 0;

    const bool isProjectedCoord = !grid->extentAndRes().isGeographic;
    auto iterCache = m_cacheGridInfo.find(grid);
    if (iterCache == m_cacheGridInfo.end()) {
        bool eastingNorthingOffset = false;
        const auto samplesPerPixel = grid->samplesPerPixel();
        int idxSampleY = -1;
        int idxSampleX = -1;
        int idxSampleZ = -1;
        for (int i = 0; i < samplesPerPixel; i++) {
            const auto desc = grid->description(i);
            if (!isProjectedCoord && desc == "latitude_offset") {
                idxSampleY = i;
                const auto unit = grid->unit(idxSampleY);
                if (!unit.empty() && unit != "arc-second") {
                    pj_log(ctx, PJ_LOG_ERROR,
                           "gridshift: Only unit=arc-second currently handled");
                    return val;
                }
            } else if (!isProjectedCoord && desc == "longitude_offset") {
                idxSampleX = i;
                const auto unit = grid->unit(idxSampleX);
                if (!unit.empty() && unit != "arc-second") {
                    pj_log(ctx, PJ_LOG_ERROR,
                           "gridshift: Only unit=arc-second currently handled");
                    return val;
                }
            } else if (isProjectedCoord && desc == "easting_offset") {
                eastingNorthingOffset = true;
                idxSampleX = i;
                const auto unit = grid->unit(idxSampleX);
                if (!unit.empty() && unit != "metre") {
                    pj_log(ctx, PJ_LOG_ERROR,
                           "gridshift: Only unit=metre currently handled");
                    return val;
                }
            } else if (isProjectedCoord && desc == "northing_offset") {
                eastingNorthingOffset = true;
                idxSampleY = i;
                const auto unit = grid->unit(idxSampleY);
                if (!unit.empty() && unit != "metre") {
                    pj_log(ctx, PJ_LOG_ERROR,
                           "gridshift: Only unit=metre currently handled");
                    return val;
                }
            } else if (desc == "ellipsoidal_height_offset" ||
                       desc == "geoid_undulation" || desc == "hydroid_height" ||
                       desc == "vertical_offset") {
                idxSampleZ = i;
                const auto unit = grid->unit(idxSampleZ);
                if (!unit.empty() && unit != "metre") {