diff options
Diffstat (limited to 'src/apply_gridshift.cpp')
| -rw-r--r-- | src/apply_gridshift.cpp | 366 |
1 files changed, 0 insertions, 366 deletions
diff --git a/src/apply_gridshift.cpp b/src/apply_gridshift.cpp deleted file mode 100644 index 39e7c3b5..00000000 --- a/src/apply_gridshift.cpp +++ /dev/null @@ -1,366 +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. - *****************************************************************************/ - -#define PJ_LIB__ - -#include <math.h> -#include <stdio.h> -#include <string.h> - -#include "proj.h" -#include "proj_internal.h" - -/************************************************************************/ -/* 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 ) - -{ - PJ_GRIDINFO **gridlist; - int grid_count; - int ret; - - gridlist = pj_gridlist_from_nadgrids( ctx, nadgrids, &grid_count ); - - if( gridlist == nullptr || grid_count == 0 ) - { - pj_dalloc( gridlist ); - return ctx->last_errno; - } - - ret = pj_apply_gridshift_3( ctx, gridlist, grid_count, inverse, - point_count, point_offset, x, y, z ); - - /* - ** Note this frees the array of grid list pointers, but not the grids - ** which is as intended. The grids themselves live on. - */ - pj_dalloc( gridlist ); - - return ret; -} - -/************************************************************************/ -/* pj_apply_gridshift_2() */ -/* */ -/* This implementation uses the gridlist from a coordinate */ -/* system definition. If the gridlist has not yet been */ -/* populated in the coordinate system definition we set it up */ -/* now. */ -/************************************************************************/ - -int pj_apply_gridshift_2( PJ *defn, int inverse, - long point_count, int point_offset, - double *x, double *y, double *z ) - -{ - if( defn->catalog_name != nullptr ) - return pj_gc_apply_gridshift( defn, inverse, point_count, point_offset, - x, y, z ); - - if( defn->gridlist == nullptr ) - { - defn->gridlist = - pj_gridlist_from_nadgrids( pj_get_ctx( defn ), - pj_param(defn->ctx, defn->params,"snadgrids").s, - &(defn->gridlist_count) ); - - if( defn->gridlist == nullptr || defn->gridlist_count == 0 ) - return defn->ctx->last_errno; - } - - return pj_apply_gridshift_3( pj_get_ctx( defn ), - defn->gridlist, defn->gridlist_count, inverse, - point_count, point_offset, x, y, z ); -} - -/************************************************************************/ -/* find_ctable() */ -/* */ -/* Determine which grid is the correct given an input coordinate. */ -/************************************************************************/ - -struct CTABLE* find_ctable(projCtx ctx, PJ_LP input, int grid_count, PJ_GRIDINFO **tables) { - int itable; - - /* keep trying till we find a table that works */ - for( itable = 0; itable < grid_count; itable++ ) - { - - PJ_GRIDINFO *gi = tables[itable]; - struct CTABLE *ct = gi->ct; - double epsilon = (fabs(ct->del.phi)+fabs(ct->del.lam))/10000.0; - /* skip tables that don't match our point at all. */ - if ( ct->ll.phi - epsilon > input.phi - || ct->ll.lam - epsilon > input.lam - || (ct->ll.phi + (ct->lim.phi-1) * ct->del.phi + epsilon < input.phi) - || (ct->ll.lam + (ct->lim.lam-1) * ct->del.lam + epsilon < input.lam) ) { - continue; - } - - /* If we have child nodes, check to see if any of them apply. */ - while( gi->child ) - { - PJ_GRIDINFO *child; - - for( child = gi->child; child != nullptr; child = child->next ) - { - struct CTABLE *ct1 = child->ct; - epsilon = (fabs(ct1->del.phi)+fabs(ct1->del.lam))/10000.0; - - if( ct1->ll.phi - epsilon > input.phi - || ct1->ll.lam - epsilon > input.lam - || (ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi + epsilon < input.phi) - || (ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam + epsilon < input.lam) ) { - continue; - } - break; - } - - /* If we didn't find a child then nothing more to do */ - if( child == nullptr ) break; - - /* Otherwise use the child, first checking it's children */ - gi = child; - ct = child->ct; - } - /* load the grid shift info if we don't have it. */ - if( ct->cvs == nullptr) { - if (!pj_gridinfo_load( ctx, gi ) ) { - pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - return nullptr; - } - } - /* if we get this far we have found a suitable grid */ - return ct; - } - - return nullptr; -} - -/************************************************************************/ -/* pj_apply_gridshift_3() */ -/* */ -/* This is the real workhorse, given a gridlist. */ -/************************************************************************/ - -int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **gridlist, int gridlist_count, - int inverse, long point_count, int point_offset, - double *x, double *y, double *z ) -{ - int i; - struct CTABLE *ct; - static int debug_count = 0; - (void) z; - - if( gridlist== nullptr || gridlist_count == 0 ) - { - pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); - return PJD_ERR_FAILED_TO_LOAD_GRID; - } - - ctx->last_errno = 0; - - for( i = 0; i < point_count; i++ ) - { - long io = i * point_offset; - PJ_LP input, output; - int itable; - - input.phi = y[io]; - input.lam = x[io]; - output.phi = HUGE_VAL; - output.lam = HUGE_VAL; - - ct = find_ctable(ctx, input, gridlist_count, gridlist); - if( ct != nullptr ) - { - output = nad_cvt( ctx, input, inverse, ct, gridlist_count, gridlist); - - if ( output.lam != HUGE_VAL && debug_count++ < 20 ) - pj_log( ctx, PJ_LOG_DEBUG_MINOR, "pj_apply_gridshift(): used %s", ct->id ); - } - - if ( output.lam == HUGE_VAL ) - { - 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 ); - for( itable = 0; itable < gridlist_count; itable++ ) - { - PJ_GRIDINFO *gi = gridlist[itable]; - if( itable == 0 ) - pj_log( ctx, PJ_LOG_DEBUG_MAJOR, " tried: %s", gi->gridname ); - else - pj_log( ctx, PJ_LOG_DEBUG_MAJOR, ",%s", gi->gridname ); - } - } - - /* - * We don't actually have any machinery currently to set the - * following macro, so this is mostly kept here to make it clear - * how we ought to operate if we wanted to make it super clear - * that an error has occurred when points are outside our available - * datum shift areas. But if this is on, we will find that "low - * value" points on the fringes of some datasets will completely - * fail causing lots of problems when it is more or less ok to - * just not apply a datum shift. So rather than deal with - * that we just fallback to no shift. (see also bug #45). - */ -#ifdef ERR_GRID_AREA_TRANSIENT_SEVERE - y[io] = HUGE_VAL; - x[io] = HUGE_VAL; -#else - /* leave x/y unshifted. */ -#endif - } - else - { - y[io] = output.phi; - x[io] = output.lam; - } - } - - return 0; -} - -/**********************************************/ -int proj_hgrid_init(PJ* P, const char *grids) { -/********************************************** - - 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. - -***********************************************/ - - /* prepend "s" to the "grids" string to allow usage with pj_param */ - char *sgrids = (char *) pj_malloc( (strlen(grids)+1+1) *sizeof(char) ); - sprintf(sgrids, "%s%s", "s", grids); - - if (P->gridlist == nullptr) { - P->gridlist = pj_gridlist_from_nadgrids( - P->ctx, - pj_param(P->ctx, P->params, sgrids).s, - &(P->gridlist_count) - ); - - if( P->gridlist == nullptr || P->gridlist_count == 0 ) { - pj_dealloc(sgrids); - return 0; - } - } - - if (P->gridlist_count == 0) { - proj_errno_set(P, PJD_ERR_FAILED_TO_LOAD_GRID); - } - - pj_dealloc(sgrids); - return P->gridlist_count; -} - -/********************************************/ -/* proj_hgrid_value() */ -/* */ -/* Return coordinate offset in grid */ -/********************************************/ -PJ_LP proj_hgrid_value(PJ *P, PJ_LP lp) { - struct CTABLE *ct; - PJ_LP out = proj_coord_error().lp; - - ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist); - if (ct == nullptr) { - pj_ctx_set_errno( P->ctx, PJD_ERR_GRID_AREA); - return out; - } - - /* normalize input to ll origin */ - lp.lam -= ct->ll.lam; - lp.phi -= ct->ll.phi; - - lp.lam = adjlon(lp.lam - M_PI) + M_PI; - - out = nad_intr(lp, ct); - - if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) { - pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA); - } - - return out; -} - -PJ_LP proj_hgrid_apply(PJ *P, PJ_LP lp, PJ_DIRECTION direction) { - struct CTABLE *ct; - int inverse; - PJ_LP out; - - out.lam = HUGE_VAL; out.phi = HUGE_VAL; - - ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist); - - if (ct == nullptr || ct->cvs == nullptr) { - if( P->gridlist_count == 1 && - strcmp(P->gridlist[0]->gridname, "null") == 0) { - // TODO: remove this particular case that is put there just to be - // able to handle longitudes outside of -180,180 - out = lp; - } else { - pj_ctx_set_errno( P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); - } - return out; - } - - inverse = direction == PJ_FWD ? 0 : 1; - out = nad_cvt(P->ctx, lp, inverse, ct, P->gridlist_count, P->gridlist); - - if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) - pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA); - - return out; - -} |