| 3161 | // --------------------------------------------------------------------------- |
| 3162 | |
| 3163 | static double read_vgrid_value(PJ_CONTEXT *ctx, const ListOfVGrids &grids, |
| 3164 | const PJ_LP &input, const double vmultiplier) { |
| 3165 | |
| 3166 | /* do not deal with NaN coordinates */ |
| 3167 | /* cppcheck-suppress duplicateExpression */ |
| 3168 | if (std::isnan(input.phi) || std::isnan(input.lam)) { |
| 3169 | return HUGE_VAL; |
| 3170 | } |
| 3171 | |
| 3172 | VerticalShiftGridSet *curGridset = nullptr; |
| 3173 | const VerticalShiftGrid *grid = nullptr; |
| 3174 | for (const auto &gridset : grids) { |
| 3175 | grid = gridset->gridAt(input.lam, input.phi); |
| 3176 | if (grid) { |
| 3177 | curGridset = gridset.get(); |
| 3178 | break; |
| 3179 | } |
| 3180 | } |
| 3181 | if (!grid) { |
| 3182 | proj_context_errno_set(ctx, PROJ_ERR_COORD_TRANSFM_OUTSIDE_GRID); |
| 3183 | return HUGE_VAL; |
| 3184 | } |
| 3185 | if (grid->isNullGrid()) { |
| 3186 | return 0; |
| 3187 | } |
| 3188 | |
| 3189 | const auto &extent = grid->extentAndRes(); |
| 3190 | if (!extent.isGeographic) { |
| 3191 | pj_log(ctx, PJ_LOG_ERROR, |
| 3192 | _("Can only handle grids referenced in a geographic CRS")); |
| 3193 | proj_context_errno_set(ctx, |
| 3194 | PROJ_ERR_INVALID_OP_FILE_NOT_FOUND_OR_INVALID); |
| 3195 | return HUGE_VAL; |
| 3196 | } |
| 3197 | |
| 3198 | /* Interpolation of a location within the grid */ |
| 3199 | double grid_x = (input.lam - extent.west) / extent.resX; |
| 3200 | if (input.lam < extent.west) { |
| 3201 | if (extent.fullWorldLongitude()) { |
| 3202 | // The first fmod goes to ]-lim, lim[ range |
| 3203 | // So we add lim again to be in ]0, 2*lim[ and fmod again |
| 3204 | grid_x = fmod(fmod(grid_x + grid->width(), grid->width()) + |
| 3205 | grid->width(), |
| 3206 | grid->width()); |
| 3207 | } else { |
| 3208 | grid_x = (input.lam + 2 * M_PI - extent.west) / extent.resX; |
| 3209 | } |
| 3210 | } else if (input.lam > extent.east) { |
| 3211 | if (extent.fullWorldLongitude()) { |
| 3212 | // The first fmod goes to ]-lim, lim[ range |
| 3213 | // So we add lim again to be in ]0, 2*lim[ and fmod again |
| 3214 | grid_x = fmod(fmod(grid_x + grid->width(), grid->width()) + |
| 3215 | grid->width(), |
| 3216 | grid->width()); |
| 3217 | } else { |
| 3218 | grid_x = (input.lam - 2 * M_PI - extent.west) / extent.resX; |
| 3219 | } |
| 3220 | } |
no test coverage detected