diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Makefile.am | 2 | ||||
| -rw-r--r-- | src/pj_apply_vgridshift.c | 200 | ||||
| -rw-r--r-- | src/pj_gridinfo.c | 140 | ||||
| -rw-r--r-- | src/pj_init.c | 7 | ||||
| -rw-r--r-- | src/pj_param.c | 4 | ||||
| -rw-r--r-- | src/pj_transform.c | 28 | ||||
| -rw-r--r-- | src/projects.h | 7 |
7 files changed, 382 insertions, 6 deletions
diff --git a/src/Makefile.am b/src/Makefile.am index 676ae916..0dcf782f 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -60,7 +60,7 @@ libproj_la_SOURCES = \ nad_cvt.c nad_init.c nad_intr.c emess.c emess.h \ pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \ geocent.c geocent.h pj_utils.c pj_gridinfo.c pj_gridlist.c \ - jniproj.c pj_mutex.c pj_initcache.c + jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c install-exec-local: 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 <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 <projects.h> +#include <string.h> +#include <math.h> + +/************************************************************************/ +/* 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; +} + diff --git a/src/pj_gridinfo.c b/src/pj_gridinfo.c index fee51bf4..8dcde5ee 100644 --- a/src/pj_gridinfo.c +++ b/src/pj_gridinfo.c @@ -291,6 +291,46 @@ int pj_gridinfo_load( PJ_GRIDINFO *gi ) return 1; } +/* -------------------------------------------------------------------- */ +/* GTX format. */ +/* -------------------------------------------------------------------- */ + else if( strcmp(gi->format,"gtx") == 0 ) + { + int words = gi->ct->lim.lam * gi->ct->lim.phi; + FILE *fid; + + fid = pj_open_lib( gi->filename, "rb" ); + + if( fid == NULL ) + { + pj_errno = -38; + return 0; + } + + fseek( fid, gi->grid_offset, SEEK_SET ); + + gi->ct->cvs = (FLP *) pj_malloc(words*sizeof(float)); + if( gi->ct->cvs == NULL ) + { + pj_errno = -38; + return 0; + } + + if( fread( gi->ct->cvs, sizeof(float), words, fid ) != words ) + { + pj_dalloc( gi->ct->cvs ); + gi->ct->cvs = NULL; + pj_errno = -38; + return 0; + } + + if( IS_LSB ) + swap_words( (char *) gi->ct->cvs, 4, words ); + + fclose( fid ); + return 1; + } + else { return 0; @@ -592,6 +632,97 @@ static int pj_gridinfo_init_ntv1( FILE * fid, PJ_GRIDINFO *gi ) } /************************************************************************/ +/* pj_gridinfo_init_gtx() */ +/* */ +/* Load a NOAA .gtx vertical datum shift file. */ +/************************************************************************/ + +static int pj_gridinfo_init_gtx( FILE * fid, PJ_GRIDINFO *gi ) + +{ + unsigned char header[40]; + struct CTABLE *ct; + double xorigin,yorigin,xstep,ystep; + int rows, columns; + + assert( sizeof(int) == 4 ); + assert( sizeof(double) == 8 ); + if( sizeof(int) != 4 || sizeof(double) != 8 ) + { + fprintf( stderr, + "basic types of inappropraiate size in nad_load_gtx()\n" ); + pj_errno = -38; + return 0; + } + +/* -------------------------------------------------------------------- */ +/* Read the header. */ +/* -------------------------------------------------------------------- */ + if( fread( header, sizeof(header), 1, fid ) != 1 ) + { + pj_errno = -38; + return 0; + } + +/* -------------------------------------------------------------------- */ +/* Regularize fields of interest and extract. */ +/* -------------------------------------------------------------------- */ + if( IS_LSB ) + { + swap_words( header+0, 8, 4 ); + swap_words( header+32, 4, 2 ); + } + + memcpy( &yorigin, header+0, 8 ); + memcpy( &xorigin, header+8, 8 ); + memcpy( &ystep, header+16, 8 ); + memcpy( &xstep, header+24, 8 ); + + memcpy( &rows, header+32, 4 ); + memcpy( &columns, header+36, 4 ); + + if( xorigin < -360 || xorigin > 360 + || yorigin < -90 || yorigin > 90 ) + { + pj_errno = -38; + printf("gtx file header has invalid extents, corrupt?\n"); + return 0; + } + +/* -------------------------------------------------------------------- */ +/* Fill in CTABLE structure. */ +/* -------------------------------------------------------------------- */ + ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE)); + strcpy( ct->id, "GTX Vertical Grid Shift File" ); + + ct->ll.lam = xorigin; + ct->ll.phi = yorigin; + ct->del.lam = xstep; + ct->del.phi = ystep; + ct->lim.lam = columns; + ct->lim.phi = rows; + + if( getenv("PROJ_DEBUG") != NULL ) + fprintf( stderr, + "GTX %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)\n", + ct->lim.lam, ct->lim.phi, + ct->ll.lam, ct->ll.phi, + ct->ll.lam + columns*xstep, ct->ll.phi + rows*ystep); + + ct->ll.lam *= DEG_TO_RAD; + ct->ll.phi *= DEG_TO_RAD; + ct->del.lam *= DEG_TO_RAD; + ct->del.phi *= DEG_TO_RAD; + ct->cvs = NULL; + + gi->ct = ct; + gi->grid_offset = 40; + gi->format = "gtx"; + + return 1; +} + +/************************************************************************/ /* pj_gridinfo_init() */ /* */ /* Open and parse header details from a datum gridshift file */ @@ -662,7 +793,14 @@ PJ_GRIDINFO *pj_gridinfo_init( const char *gridname ) { pj_gridinfo_init_ntv2( fp, gilist ); } - + + else if( strlen(gridname) > 4 + && (strcmp(gridname+strlen(gridname)-3,"gtx") == 0 + || strcmp(gridname+strlen(gridname)-3,"GTX") == 0) ) + { + pj_gridinfo_init_gtx( fp, gilist ); + } + else { struct CTABLE *ct = nad_ctable_init( fp ); diff --git a/src/pj_init.c b/src/pj_init.c index 87d5ca8d..91230a70 100644 --- a/src/pj_init.c +++ b/src/pj_init.c @@ -271,9 +271,13 @@ pj_init(int argc, char **argv) { PIN->is_long_wrap_set = 0; PIN->long_wrap_center = 0.0; strcpy( PIN->axis, "enu" ); + PIN->gridlist = NULL; PIN->gridlist_count = 0; + PIN->vgridlist_geoid = NULL; + PIN->vgridlist_geoid_count = 0; + /* set datum parameters */ if (pj_datum_set(start, PIN)) goto bum_call; @@ -307,6 +311,9 @@ pj_init(int argc, char **argv) { PIN->over = pj_param(start, "bover").i; /* longitude center for wrapping */ + PIN->has_geoid_vgrids = pj_param(start, "tgeoidgrids").i; + + /* longitude center for wrapping */ PIN->is_long_wrap_set = pj_param(start, "tlon_wrap").i; if (PIN->is_long_wrap_set) PIN->long_wrap_center = pj_param(start, "rlon_wrap").f; diff --git a/src/pj_param.c b/src/pj_param.c index a95010d3..c18554a0 100644 --- a/src/pj_param.c +++ b/src/pj_param.c @@ -34,7 +34,7 @@ pj_mkparam(char *str) { /************************************************************************/ PVALUE /* test for presence or get parameter value */ -pj_param(paralist *pl, char *opt) { +pj_param(paralist *pl, const char *opt) { int type; unsigned l; PVALUE value; @@ -63,7 +63,7 @@ pj_param(paralist *pl, char *opt) { value.f = dmstor(opt, 0); break; case 's': /* char string */ - value.s = opt; + value.s = (char *) opt; break; case 'b': /* boolean */ switch (*opt) { diff --git a/src/pj_transform.c b/src/pj_transform.c index 94f7d6ed..33c1dcd1 100644 --- a/src/pj_transform.c +++ b/src/pj_transform.c @@ -200,6 +200,19 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset, } /* -------------------------------------------------------------------- */ +/* Do we need to translate from geoid to ellipsoidal vertical */ +/* datum? */ +/* -------------------------------------------------------------------- */ + if( srcdefn->has_geoid_vgrids ) + { + if( !pj_apply_vgridshift( srcdefn, "sgeoidgrids", + &(srcdefn->vgridlist_geoid), + &(srcdefn->vgridlist_geoid_count), + 0, point_count, point_offset, x, y, z ) ) + return pj_errno; + } + +/* -------------------------------------------------------------------- */ /* Convert datums if needed, and possible. */ /* -------------------------------------------------------------------- */ if( pj_datum_transform( srcdefn, dstdefn, point_count, point_offset, @@ -207,6 +220,19 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset, return pj_errno; /* -------------------------------------------------------------------- */ +/* Do we need to translate from geoid to ellipsoidal vertical */ +/* datum? */ +/* -------------------------------------------------------------------- */ + if( dstdefn->has_geoid_vgrids ) + { + if( !pj_apply_vgridshift( dstdefn, "sgeoidgrids", + &(dstdefn->vgridlist_geoid), + &(dstdefn->vgridlist_geoid_count), + 0, point_count, point_offset, x, y, z ) ) + return pj_errno; + } + +/* -------------------------------------------------------------------- */ /* But if they are staying lat long, adjust for the prime */ /* meridian if there is one in effect. */ /* -------------------------------------------------------------------- */ @@ -618,7 +644,7 @@ int pj_datum_transform( PJ *srcdefn, PJ *dstdefn, dst_a = SRS_WGS84_SEMIMAJOR; dst_es = SRS_WGS84_ESQUARED; } - + /* ==================================================================== */ /* Do we need to go through geocentric coordinates? */ /* ==================================================================== */ diff --git a/src/projects.h b/src/projects.h index fcc9e48d..6a91ebda 100644 --- a/src/projects.h +++ b/src/projects.h @@ -236,6 +236,11 @@ typedef struct PJconsts { double datum_params[7]; struct _pj_gi **gridlist; int gridlist_count; + + int has_geoid_vgrids; + struct _pj_gi **vgridlist_geoid; + int vgridlist_geoid_count; + double from_greenwich; /* prime meridian offset (in radians) */ double long_wrap_center; /* 0.0 for -180 to 180, actually in radians*/ int is_long_wrap_set; @@ -342,7 +347,7 @@ void set_rtodms(int, int); char *rtodms(char *, double, int, int); double adjlon(double); double aacos(double), aasin(double), asqrt(double), aatan2(double, double); -PVALUE pj_param(paralist *, char *); +PVALUE pj_param(paralist *, const char *); paralist *pj_mkparam(char *); int pj_ell_set(paralist *, double *, double *); int pj_datum_set(paralist *, PJ *); |