aboutsummaryrefslogtreecommitdiff
path: root/src/apply_vgridshift.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/apply_vgridshift.cpp')
-rw-r--r--src/apply_vgridshift.cpp217
1 files changed, 0 insertions, 217 deletions
diff --git a/src/apply_vgridshift.cpp b/src/apply_vgridshift.cpp
deleted file mode 100644
index b0136e5c..00000000
--- a/src/apply_vgridshift.cpp
+++ /dev/null
@@ -1,217 +0,0 @@
-/******************************************************************************
- * Project: PROJ.4
- * Purpose: Apply vertical datum shifts based on grid shift files, normally
- * geoid grids mapping WGS84 to NAVD88 or something similar.
- * Author: Frank Warmerdam, warmerdam@pobox.com
- *
- ******************************************************************************
- * Copyright (c) 2010, Frank Warmerdam <warmerdam@pobox.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
-
-#include <assert.h>
-#include <stdio.h>
-#include <string.h>
-
-#include <math.h>
-#include "proj_internal.h"
-#include "proj/internal/internal.hpp"
-#include "grids.hpp"
-
-NS_PROJ_START
-
-static double read_vgrid_value(const ListOfVGrids& grids, PJ_LP input, double vmultiplier) {
-
- /* do not deal with NaN coordinates */
- /* cppcheck-suppress duplicateExpression */
- if( isnan(input.phi) || isnan(input.lam) )
- {
- return HUGE_VAL;
- }
-
- const VerticalShiftGrid* grid = nullptr;
- for( const auto& gridset: grids )
- {
- grid = gridset->gridAt(input.lam, input.phi);
- if( grid )
- break;
- }
- if( !grid )
- {
- return HUGE_VAL;
- }
-
- const auto& extent = grid->extentAndRes();
-
- /* Interpolation a location within the grid */
- double grid_x = (input.lam - extent.westLon) / extent.resLon;
- if( extent.fullWorldLongitude() ) {
- // The first fmod goes to ]-lim, lim[ range
- // So we add lim again to be in ]0, 2*lim[ and fmod again
- grid_x = fmod(
- fmod(grid_x + grid->width(), grid->width()) + grid->width(),
- grid->width());
- }
- double grid_y = (input.phi - extent.southLat) / extent.resLat;
- int grid_ix = static_cast<int>(lround(floor(grid_x)));
- assert(grid_ix >= 0 && grid_ix < grid->width());
- int grid_iy = static_cast<int>(lround(floor(grid_y)));
- assert(grid_iy >= 0 && grid_iy < grid->height());
- grid_x -= grid_ix;
- grid_y -= grid_iy;
-
- int grid_ix2 = grid_ix + 1;
- if( grid_ix2 >= grid->width() ) {
- if( extent.fullWorldLongitude() ) {
- grid_ix2 = 0;
- } else {
- grid_ix2 = grid->width() - 1;
- }
- }
- int grid_iy2 = grid_iy + 1;
- if( grid_iy2 >= grid->height() )
- grid_iy2 = grid->height() - 1;
-
- float value_a = 0;
- float value_b = 0;
- float value_c = 0;
- float value_d = 0;
- if( !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) )
- {
- 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 proj_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 )
- {
- pj_ctx_set_errno( P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
- return {};
- }
- }
- else
- {
- grids.emplace_back(std::move(gridSet));
- }
- }
-
- return grids;
-}
-
-/***********************************************/
-double proj_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(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;
-}
-
-NS_PROJ_END
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 17:57:39 +0000