| 3590 | // --------------------------------------------------------------------------- |
| 3591 | |
| 3592 | static double read_vgrid_value(PJ_CONTEXT *ctx, const ListOfVGrids &grids, |
| 3593 | const PJ_LP &input, const double vmultiplier) { |
| 3594 | |
| 3595 | /* do not deal with NaN coordinates */ |
| 3596 | /* cppcheck-suppress duplicateExpression */ |
| 3597 | if (std::isnan(input.phi) || std::isnan(input.lam)) { |
| 3598 | return HUGE_VAL; |
| 3599 | } |
| 3600 | |
| 3601 | VerticalShiftGridSet *curGridset = nullptr; |
| 3602 | const VerticalShiftGrid *grid = nullptr; |
| 3603 | for (const auto &gridset : grids) { |
| 3604 | grid = gridset->gridAt(input.lam, input.phi); |
| 3605 | if (grid) { |
| 3606 | curGridset = gridset.get(); |
| 3607 | break; |
| 3608 | } |
| 3609 | } |
| 3610 | if (!grid) { |
| 3611 | proj_context_errno_set(ctx, PROJ_ERR_COORD_TRANSFM_OUTSIDE_GRID); |
| 3612 | return HUGE_VAL; |
| 3613 | } |
| 3614 | if (grid->isNullGrid()) { |
| 3615 | return 0; |
| 3616 | } |
| 3617 | |
| 3618 | const auto &extent = grid->extentAndRes(); |
| 3619 | if (!extent.isGeographic) { |
| 3620 | pj_log(ctx, PJ_LOG_ERROR, |
| 3621 | _("Can only handle grids referenced in a geographic CRS")); |
| 3622 | proj_context_errno_set(ctx, |
| 3623 | PROJ_ERR_INVALID_OP_FILE_NOT_FOUND_OR_INVALID); |
| 3624 | return HUGE_VAL; |
| 3625 | } |
| 3626 | |
| 3627 | /* Interpolation of a location within the grid */ |
| 3628 | double grid_x = (input.lam - extent.west) * extent.invResX; |
| 3629 | if (input.lam < extent.west) { |
| 3630 | if (extent.fullWorldLongitude()) { |
| 3631 | // The first fmod goes to ]-lim, lim[ range |
| 3632 | // So we add lim again to be in ]0, 2*lim[ and fmod again |
| 3633 | grid_x = fmod(fmod(grid_x + grid->width(), grid->width()) + |
| 3634 | grid->width(), |
| 3635 | grid->width()); |
| 3636 | } else { |
| 3637 | grid_x = (input.lam + 2 * M_PI - extent.west) * extent.invResX; |
| 3638 | } |
| 3639 | } else if (input.lam > extent.east) { |
| 3640 | if (extent.fullWorldLongitude()) { |
| 3641 | // The first fmod goes to ]-lim, lim[ range |
| 3642 | // So we add lim again to be in ]0, 2*lim[ and fmod again |
| 3643 | grid_x = fmod(fmod(grid_x + grid->width(), grid->width()) + |
| 3644 | grid->width(), |
| 3645 | grid->width()); |
| 3646 | } else { |
| 3647 | grid_x = (input.lam - 2 * M_PI - extent.west) * extent.invResX; |
| 3648 | } |
| 3649 | } |
no test coverage detected