diff options
Diffstat (limited to 'src/apply_gridshift.cpp')
| -rw-r--r-- | src/apply_gridshift.cpp | 397 |
1 files changed, 0 insertions, 397 deletions
diff --git a/src/apply_gridshift.cpp b/src/apply_gridshift.cpp deleted file mode 100644 index 4ef86fc0..00000000 --- a/src/apply_gridshift.cpp +++ /dev/null @@ -1,397 +0,0 @@ -/****************************************************************************** - * Project: PROJ.4 - * Purpose: Apply datum shifts based on grid shift files (normally NAD27 to - * NAD83 or the reverse). This module is responsible for keeping - * a list of loaded grids, and calling with each one that is - * allowed for a given datum (expressed as the nadgrids= parameter). - * Author: Frank Warmerdam, warmerdam@pobox.com - * - ****************************************************************************** - * Copyright (c) 2000, 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 <math.h> -#include <stdio.h> -#include <string.h> - -#include <limits> - -#include "proj.h" -#include "proj_internal.h" -#include "proj/internal/internal.hpp" -#include "grids.hpp" - -NS_PROJ_START - -// --------------------------------------------------------------------------- - -static const HorizontalShiftGrid* findGrid(const ListOfHGrids& grids, - PJ_LP input) -{ - for( const auto& gridset: grids ) - { - auto grid = gridset->gridAt(input.lam, input.phi); - if( grid ) - 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 ) - { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return {}; - } - } - else - { - list.emplace_back(std::move(gridSet)); - } - } - return list; -} - - -/**********************************************/ -ListOfHGrids proj_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; - -static PJ_LP nad_intr(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 = isnan(t.lam) ? 0 : (pj_int32)lround(floor(t.lam)); - t.phi /= extent.resLat; - indx.phi = 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 > 0.99999999999) { - ++indx.lam; - frct.lam = 0.; - } else - return val; - } else if ((in = indx.lam + 1) >= grid->width()) { - if (in == grid->width() && frct.lam < 1e-11) { - --indx.lam; - frct.lam = 1.; - } else - return val; - } - if (indx.phi < 0) { - if (indx.phi == -1 && frct.phi > 0.99999999999) { - ++indx.phi; - frct.phi = 0.; - } else - return val; - } else if ((in = indx.phi + 1) >= grid->height()) { - if (in == grid->height() && frct.phi < 1e-11) { - --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 nad_cvt(projCtx ctx, PJ_LP in, int inverse, const HorizontalShiftGrid* grid, - const ListOfHGrids& grids) { - PJ_LP t, tb,del, dif; - int i = MAX_ITERATIONS; - const double toltol = TOL*TOL; - - if (in.lam == HUGE_VAL) - return in; - - /* normalize input to ll origin */ - tb = in; - const auto* extent = &(grid->extentAndRes()); - tb.lam -= extent->westLon; - tb.phi -= extent->southLat; - - tb.lam = adjlon (tb.lam - M_PI) + M_PI; - - t = nad_intr (tb, grid, true); - if (t.lam == HUGE_VAL) - return t; - - if (!inverse) { - in.lam += t.lam; - in.phi += t.phi; - return in; - } - - t.lam = tb.lam - t.lam; - t.phi = tb.phi - t.phi; - - do { - del = nad_intr(t, grid, true); - - /* 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); - 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; - tb.phi -= extent->southLat; - tb.lam = adjlon (tb.lam - M_PI) + M_PI; - 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; -} - -/********************************************/ -/* proj_hgrid_value() */ -/* */ -/* Return coordinate offset in grid */ -/********************************************/ -PJ_LP proj_hgrid_value(PJ *P, const ListOfHGrids& grids, PJ_LP lp) { - PJ_LP out = proj_coord_error().lp; - - const auto grid = findGrid(grids, lp); - if( !grid ) { - pj_ctx_set_errno( P->ctx, PJD_ERR_GRID_AREA ); - return out; - } - - /* normalize input to ll origin */ - const auto& extent = grid->extentAndRes(); - lp.lam -= extent.westLon; - lp.phi -= extent.southLat; - - lp.lam = adjlon(lp.lam - M_PI) + M_PI; - - out = nad_intr(lp, grid, false); - - if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) { - pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA); - } - - return out; -} - -static -PJ_LP proj_hgrid_apply_internal(PJ_CONTEXT* ctx, - PJ_LP lp, - PJ_DIRECTION direction, - const ListOfHGrids& grids) -{ - PJ_LP out; - - out.lam = HUGE_VAL; - out.phi = HUGE_VAL; - - const auto grid = findGrid(grids, lp); - if( !grid ) { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return out; - } - if( grid->isNullGrid() ) - { - return lp; - } - - int inverse = direction == PJ_FWD ? 0 : 1; - out = nad_cvt(ctx, lp, inverse, grid, grids); - - if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) - pj_ctx_set_errno(ctx, PJD_ERR_GRID_AREA); - - return out; -} - -PJ_LP proj_hgrid_apply(PJ *P, const ListOfHGrids& grids, PJ_LP lp, PJ_DIRECTION direction) { - return proj_hgrid_apply_internal(P->ctx, lp, direction, grids); -} - - -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 = proj_hgrid_apply_internal(ctx, input, inverse ? PJ_INV : PJ_FWD, hgrids); - - 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; -} |