From b08809dfede44f0e32bcee930d408f9f5e6bd7cf Mon Sep 17 00:00:00 2001 From: Frank Warmerdam Date: Tue, 11 May 2010 04:05:04 +0000 Subject: preliminary addition of vertical datum shifting capability git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@1839 4e78687f-474d-0410-85f9-8d5e500ac6b2 --- src/pj_apply_vgridshift.c | 200 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 200 insertions(+) create mode 100644 src/pj_apply_vgridshift.c (limited to 'src/pj_apply_vgridshift.c') diff --git a/src/pj_apply_vgridshift.c b/src/pj_apply_vgridshift.c new file mode 100644 index 00000000..7696ea19 --- /dev/null +++ b/src/pj_apply_vgridshift.c @@ -0,0 +1,200 @@ +/****************************************************************************** + * $Id: pj_apply_gridshift.c 1831 2010-03-16 12:44:36Z warmerdam $ + * + * 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 + * + * 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 +#include +#include + +/************************************************************************/ +/* pj_apply_vgridshift() */ +/* */ +/* This implmentation takes 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_vgridshift( PJ *defn, const char *listname, + PJ_GRIDINFO ***gridlist_p, + int *gridlist_count_p, + int inverse, + long point_count, int point_offset, + double *x, double *y, double *z ) + +{ + int i; + int debug_flag = getenv( "PROJ_DEBUG" ) != NULL; + static int debug_count = 0; + PJ_GRIDINFO **tables; + + if( *gridlist_p == NULL ) + { + *gridlist_p = + pj_gridlist_from_nadgrids( pj_param(defn->params,listname).s, + gridlist_count_p ); + + if( *gridlist_p == NULL || *gridlist_count_p == 0 ) + return pj_errno; + } + + if( tables == NULL || *gridlist_count_p == 0 ) + { + pj_errno = -38; + return -38; + } + + tables = *gridlist_p; + pj_errno = 0; + + for( i = 0; i < point_count; i++ ) + { + long io = i * point_offset; + LP input, output; + int itable; + double value = HUGE_VAL; + + input.phi = y[io]; + input.lam = x[io]; + + /* keep trying till we find a table that works */ + for( itable = 0; itable < *gridlist_count_p; itable++ ) + { + PJ_GRIDINFO *gi = tables[itable]; + struct CTABLE *ct = gi->ct; + double grid_x, grid_y; + int grid_ix, grid_iy; + float *cvs; + + /* skip tables that don't match our point at all. */ + if( ct->ll.phi > input.phi || ct->ll.lam > input.lam + || ct->ll.phi + (ct->lim.phi-1) * ct->del.phi < input.phi + || ct->ll.lam + (ct->lim.lam-1) * ct->del.lam < input.lam ) + continue; + + /* If we have child nodes, check to see if any of them apply. */ + if( gi->child != NULL ) + { + PJ_GRIDINFO *child; + + for( child = gi->child; child != NULL; child = child->next ) + { + struct CTABLE *ct1 = child->ct; + + if( ct1->ll.phi > input.phi || ct1->ll.lam > input.lam + || ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi < input.phi + || ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam < input.lam) + continue; + + break; + } + + /* we found a more refined child node to use */ + if( child != NULL ) + { + gi = child; + ct = child->ct; + } + } + + /* load the grid shift info if we don't have it. */ + if( ct->cvs == NULL && !pj_gridinfo_load( gi ) ) + { + pj_errno = -38; + return -38; + } + + /* Interpolation a location within the grid */ + grid_x = (input.lam - ct->ll.lam) / ct->del.lam; + grid_y = (input.phi - ct->ll.phi) / ct->del.phi; + grid_ix = (int) floor(grid_x); + grid_iy = (int) floor(grid_y); + grid_x -= grid_ix; + grid_y -= grid_iy; + + cvs = (float *) ct->cvs; + value = cvs[grid_ix + grid_iy * ct->lim.lam] + * (1.0-grid_x) * (1.0-grid_y) + + cvs[grid_ix + 1 + grid_iy * ct->lim.lam] + * (grid_x) * (1.0-grid_y) + + cvs[grid_ix + (grid_iy+1) * ct->lim.lam] + * (1.0-grid_x) * (grid_y) + + cvs[grid_ix + 1 + (grid_iy+1) * ct->lim.lam] + * (grid_x) * (grid_y); + + if( value > 1000 || value < -1000 ) /* nodata? */ + value = HUGE_VAL; + else + { + if( inverse ) + z[io] += value; + else + z[io] -= value; + } + + if( value != HUGE_VAL ) + { + if( debug_flag && debug_count++ < 20 ) + fprintf( stderr, + "pj_apply_gridshift(): used %s\n", + ct->id ); + break; + } + } + + if( value == HUGE_VAL ) + { + if( debug_flag ) + { + fprintf( stderr, + "pj_apply_vgridshift(): failed to find a grid shift table for\n" + " location (%.7fdW,%.7fdN)\n", + x[io] * RAD_TO_DEG, + y[io] * RAD_TO_DEG ); + for( itable = 0; itable < *gridlist_count_p; itable++ ) + { + PJ_GRIDINFO *gi = tables[itable]; + if( itable == 0 ) + fprintf( stderr, " tried: %s", gi->gridname ); + else + fprintf( stderr, ",%s", gi->gridname ); + } + fprintf( stderr, "\n" ); + } + + pj_errno = PJD_ERR_GRID_AREA; + return PJD_ERR_GRID_AREA; + } + } + + return 0; +} + -- cgit v1.2.3