diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/PJ_hgridshift.c | 52 | ||||
| -rw-r--r-- | src/PJ_vgridshift.c | 25 | ||||
| -rw-r--r-- | src/pj_apply_gridshift.c | 275 | ||||
| -rw-r--r-- | src/pj_apply_vgridshift.c | 305 | ||||
| -rw-r--r-- | src/proj_internal.h | 7 |
5 files changed, 398 insertions, 266 deletions
diff --git a/src/PJ_hgridshift.c b/src/PJ_hgridshift.c index 0adc9e00..659039ab 100644 --- a/src/PJ_hgridshift.c +++ b/src/PJ_hgridshift.c @@ -4,7 +4,6 @@ PROJ_HEAD(hgridshift, "Horizontal grid shift"); - static XYZ forward_3d(LPZ lpz, PJ *P) { PJ_TRIPLET point; point.lpz = lpz; @@ -12,9 +11,7 @@ static XYZ forward_3d(LPZ lpz, PJ *P) { if (P->gridlist != NULL) { /* Only try the gridshift if at least one grid is loaded, * otherwise just pass the coordinate through unchanged. */ - pj_apply_gridshift_3( P->ctx, P->gridlist, - P->gridlist_count, 1, 1, 0, - &point.xyz.x, &point.xyz.y, &point.xyz.z ); + point.lp = proj_hgrid_apply(P, point.lp, PJ_FWD); } return point.xyz; @@ -28,9 +25,7 @@ static LPZ reverse_3d(XYZ xyz, PJ *P) { if (P->gridlist != NULL) { /* Only try the gridshift if at least one grid is loaded, * otherwise just pass the coordinate through unchanged. */ - pj_apply_gridshift_3( P->ctx, P->gridlist, - P->gridlist_count, 0, 1, 0, - &point.xyz.x, &point.xyz.y, &point.xyz.z ); + point.lp = proj_hgrid_apply(P, point.lp, PJ_INV); } return point.lpz; @@ -52,29 +47,8 @@ static PJ_OBS reverse_obs(PJ_OBS obs, PJ *P) { -#if 0 -static XY forward_xy(LP lp, PJ *P) { - PJ_TRIPLET point; - point.lp = lp; - point.lpz.z = 0; - point.xyz = forward_3d (point.lpz, P); - return point.xy; -} - - -static LP reverse_lp(XY xy, PJ *P) { - PJ_TRIPLET point; - point.xy = xy; - point.xyz.z = 0; - point.lpz = reverse_3d (point.xyz, P); - return point.lp; -} -#endif - - - PJ *PROJECTION(hgridshift) { - + P->fwdobs = forward_obs; P->invobs = reverse_obs; P->fwd3d = forward_3d; @@ -90,14 +64,12 @@ PJ *PROJECTION(hgridshift) { return pj_default_destructor (P, PJD_ERR_NO_ARGS); } - /* Build gridlist. P->gridlist can be empty if +grids only ask for optional grids. */ - P->gridlist = pj_gridlist_from_nadgrids( P->ctx, pj_param(P->ctx, P->params, "sgrids").s, - &(P->gridlist_count) ); + proj_hgrid_init(P, "grids"); /* Was gridlist compiled properly? */ - if ( pj_ctx_get_errno(pj_get_ctx(P)) ) { + if ( proj_errno(P) ) { proj_log_error(P, "hgridshift: could not find required grid(s)."); - return pj_default_destructor (P, PJD_ERR_FAILED_TO_LOAD_GRID); + return pj_default_destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); } return P; @@ -119,26 +91,28 @@ int pj_hgridshift_selftest (void) { proj_destroy (P); return 99; } - + /* fail on purpose: open non-existing grid */ P = proj_create(PJ_DEFAULT_CTX, "+proj=hgridshift +grids=@nonexistinggrid.gsb,anothernonexistinggrid.gsb"); if (0!=P) { proj_destroy (P); return 999; } - + /* Failure most likely means the grid is missing */ P = proj_create(PJ_DEFAULT_CTX, "+proj=hgridshift +grids=nzgd2kgrid0005.gsb +ellps=GRS80"); if (0==P) return 10; - + a = proj_obs_null; a.coo.lpz.lam = PJ_TORAD(173); a.coo.lpz.phi = PJ_TORAD(-45); - + dist = proj_roundtrip (P, PJ_FWD, 1, a.coo); - if (dist > 0.00000001) + if (dist > 0.00000001) { + printf("dist: %f\n",dist); return 1; + } expect.coo.lpz.lam = PJ_TORAD(172.999892181021551); expect.coo.lpz.phi = PJ_TORAD(-45.001620431954613); diff --git a/src/PJ_vgridshift.c b/src/PJ_vgridshift.c index ededd544..691f791b 100644 --- a/src/PJ_vgridshift.c +++ b/src/PJ_vgridshift.c @@ -9,14 +9,10 @@ static XYZ forward_3d(LPZ lpz, PJ *P) { PJ_TRIPLET point; point.lpz = lpz; - if (P->gridlist != NULL) { + if (P->vgridlist_geoid != NULL) { /* Only try the gridshift if at least one grid is loaded, * otherwise just pass the coordinate through unchanged. */ - pj_apply_vgridshift( P, "sgrids", - &(P->gridlist), - &(P->gridlist_count), - 1, 1, 0, - &point.xyz.x, &point.xyz.y, &point.xyz.z ); + point.xyz.z -= proj_vgrid_value(P, point.lp); } return point.xyz; @@ -27,14 +23,10 @@ static LPZ reverse_3d(XYZ xyz, PJ *P) { PJ_TRIPLET point; point.xyz = xyz; - if (P->gridlist != NULL) { + if (P->vgridlist_geoid != NULL) { /* Only try the gridshift if at least one grid is loaded, * otherwise just pass the coordinate through unchanged. */ - pj_apply_vgridshift( P, "sgrids", - &(P->gridlist), - &(P->gridlist_count), - 0, 1, 0, - &point.xyz.x, &point.xyz.y, &point.xyz.z ); + point.xyz.z += proj_vgrid_value(P, point.lp); } return point.lpz; @@ -61,14 +53,13 @@ PJ *PROJECTION(vgridshift) { return pj_default_destructor(P, PJD_ERR_NO_ARGS); } - /* Build gridlist. P->gridlist can be empty if +grids only ask for optional grids. */ - P->gridlist = pj_gridlist_from_nadgrids( P->ctx, pj_param(P->ctx, P->params, "sgrids").s, - &(P->gridlist_count) ); + /* Build gridlist. P->vgridlist_geoid can be empty if +grids only ask for optional grids. */ + proj_vgrid_init(P, "grids"); /* Was gridlist compiled properly? */ - if ( pj_ctx_get_errno(P->ctx) ) { + if ( proj_errno(P) ) { proj_log_error(P, "vgridshift: could not find required grid(s)."); - return pj_default_destructor(P, -38); + return pj_default_destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); } P->fwdobs = forward_obs; diff --git a/src/pj_apply_gridshift.c b/src/pj_apply_gridshift.c index 91e2de26..45887abd 100644 --- a/src/pj_apply_gridshift.c +++ b/src/pj_apply_gridshift.c @@ -2,7 +2,7 @@ * 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 + * 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 * @@ -30,9 +30,10 @@ #define PJ_LIB__ -#include <projects.h> #include <string.h> #include <math.h> +#include "proj_internal.h" +#include "projects.h" /************************************************************************/ /* pj_apply_gridshift() */ @@ -43,7 +44,7 @@ /* it to honour our public api. */ /************************************************************************/ -int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse, +int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse, long point_count, int point_offset, double *x, double *y, double *z ) @@ -51,16 +52,16 @@ int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse, PJ_GRIDINFO **gridlist; int grid_count; int ret; - + gridlist = pj_gridlist_from_nadgrids( ctx, nadgrids, &grid_count ); if( gridlist == NULL || grid_count == 0 ) return ctx->last_errno; - ret = pj_apply_gridshift_3( ctx, gridlist, grid_count, inverse, + 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. */ @@ -72,13 +73,13 @@ int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse, /************************************************************************/ /* pj_apply_gridshift_2() */ /* */ -/* This implementation takes uses the gridlist from a coordinate */ +/* 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, +int pj_apply_gridshift_2( PJ *defn, int inverse, long point_count, int point_offset, double *x, double *y, double *z ) @@ -86,10 +87,10 @@ int pj_apply_gridshift_2( PJ *defn, int inverse, if( defn->catalog_name != NULL ) return pj_gc_apply_gridshift( defn, inverse, point_count, point_offset, x, y, z ); - + if( defn->gridlist == NULL ) { - defn->gridlist = + defn->gridlist = pj_gridlist_from_nadgrids( pj_get_ctx( defn ), pj_param(defn->ctx, defn->params,"snadgrids").s, &(defn->gridlist_count) ); @@ -97,12 +98,76 @@ int pj_apply_gridshift_2( PJ *defn, int inverse, if( defn->gridlist == NULL || defn->gridlist_count == 0 ) return defn->ctx->last_errno; } - + return pj_apply_gridshift_3( pj_get_ctx( defn ), - defn->gridlist, defn->gridlist_count, inverse, + 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. */ +/************************************************************************/ + +static struct CTABLE* find_ctable(projCtx ctx, LP input, int grid_count, PJ_GRIDINFO **tables) { + int itable; + double epsilon; + struct CTABLE *ct = NULL; + + /* keep trying till we find a table that works */ + for( itable = 0; itable < grid_count; itable++ ) + { + + PJ_GRIDINFO *gi = tables[itable]; + ct = gi->ct; + 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 != NULL; 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 == NULL ) 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 == NULL) { + if (!pj_gridinfo_load( ctx, gi ) ) { + pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); + return NULL; + } + } + /* if we get this far we have found a suitable grid */ + break; + } + return ct; +} /************************************************************************/ /* pj_apply_gridshift_3() */ @@ -113,16 +178,16 @@ int pj_apply_gridshift_2( PJ *defn, int inverse, int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_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( tables == NULL || grid_count == 0 ) { - pj_ctx_set_errno( ctx, -38); - return -38; + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return PJD_ERR_FAILED_TO_LOAD_GRID; } ctx->last_errno = 0; @@ -138,100 +203,39 @@ int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_count, output.phi = HUGE_VAL; output.lam = HUGE_VAL; - /* 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 != NULL; 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 == NULL ) break; - - /* Otherwise use the child, first checking it's children */ - - gi = child; - ct = child->ct; - } + ct = find_ctable(ctx, input, grid_count, tables); + output = nad_cvt( input, inverse, ct ); - /* load the grid shift info if we don't have it. */ - if( ct->cvs == NULL && !pj_gridinfo_load( ctx, gi ) ) - { - pj_ctx_set_errno( ctx, -38 ); - return -38; - } - - output = nad_cvt( input, inverse, ct ); - if( output.lam != HUGE_VAL ) - { - if( debug_count++ < 20 ) - pj_log( ctx, PJ_LOG_DEBUG_MINOR, - "pj_apply_gridshift(): used %s", ct->id ); - break; - } - } + 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 ( 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, + x[io] * RAD_TO_DEG, y[io] * RAD_TO_DEG ); for( itable = 0; itable < grid_count; itable++ ) { PJ_GRIDINFO *gi = tables[itable]; if( itable == 0 ) - pj_log( ctx, PJ_LOG_DEBUG_MAJOR, - " tried: %s", gi->gridname ); + pj_log( ctx, PJ_LOG_DEBUG_MAJOR, " tried: %s", gi->gridname ); else - pj_log( ctx, PJ_LOG_DEBUG_MAJOR, - ",%s", gi->gridname ); + 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 + /* + * 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 + * 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). */ @@ -252,3 +256,92 @@ int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_count, 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) *sizeof(char) ); + sprintf(sgrids, "%s%s", "s", grids); + + if (P->gridlist == NULL) { + P->gridlist = pj_gridlist_from_nadgrids( + P->ctx, + pj_param(P->ctx, P->params, sgrids).s, + &(P->gridlist_count) + ); + + if( P->gridlist == NULL || 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 */ +/********************************************/ +LP proj_hgrid_value(PJ *P, LP lp) { + struct CTABLE *ct; + LP out; + + ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist); + + /* 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; +} + +LP proj_hgrid_apply(PJ *P, LP lp, PJ_DIRECTION direction) { + struct CTABLE *ct; + int inverse; + LP out; + + out.lam = HUGE_VAL; out.phi = HUGE_VAL; + + ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist); + + if (ct == NULL || ct->cvs == NULL) { + pj_ctx_set_errno( P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); + return out; + } + + inverse = direction == PJ_FWD ? 0 : 1; + out = nad_cvt(lp, inverse, ct); + + if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) + pj_ctx_set_errno(P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + + return out; + +} diff --git a/src/pj_apply_vgridshift.c b/src/pj_apply_vgridshift.c index 35047a19..3c7cc210 100644 --- a/src/pj_apply_vgridshift.c +++ b/src/pj_apply_vgridshift.c @@ -28,23 +28,120 @@ #define PJ_LIB__ -#include <projects.h> #include <string.h> #include <math.h> +#include "proj_internal.h" +#include "projects.h" + +static double pj_read_vgrid_value( PJ *defn, LP input, int *gridlist_count_p, PJ_GRIDINFO **tables, struct CTABLE *ct) { + int itable = 0; + double value = HUGE_VAL; + double grid_x, grid_y; + int grid_ix, grid_iy; + int grid_ix2, grid_iy2; + float *cvs; + /* do not deal with NaN coordinates */ + if( input.phi != input.phi || input.lam != input.lam ) + itable = *gridlist_count_p; + + /* keep trying till we find a table that works */ + for ( ; itable < *gridlist_count_p; itable++ ) + { + PJ_GRIDINFO *gi = tables[itable]; + + ct = gi->ct; + + /* 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. */ + while( 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 didn't find a more refined child node to use, so go with current grid */ + if( child == NULL ) + { + break; + } + + /* Otherwise let's try for childrens children .. */ + gi = child; + ct = child->ct; + } + + /* load the grid shift info if we don't have it. */ + if( ct->cvs == NULL && !pj_gridinfo_load( pj_get_ctx(defn), gi ) ) + { + pj_ctx_set_errno( defn->ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); + return PJD_ERR_FAILED_TO_LOAD_GRID; + } + + } + + /* 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; + + grid_ix2 = grid_ix + 1; + if( grid_ix2 >= ct->lim.lam ) + grid_ix2 = ct->lim.lam - 1; + grid_iy2 = grid_iy + 1; + if( grid_iy2 >= ct->lim.phi ) + grid_iy2 = ct->lim.phi - 1; + + cvs = (float *) ct->cvs; + value = cvs[grid_ix + grid_iy * ct->lim.lam] + * (1.0-grid_x) * (1.0-grid_y) + + cvs[grid_ix2 + grid_iy * ct->lim.lam] + * (grid_x) * (1.0-grid_y) + + cvs[grid_ix + grid_iy2 * ct->lim.lam] + * (1.0-grid_x) * (grid_y) + + cvs[grid_ix2 + grid_iy2 * ct->lim.lam] + * (grid_x) * (grid_y); + + /* nodata? */ + /* GTX official nodata value if -88.88880f, but some grids also */ + /* use other big values for nodata (e.g naptrans2008.gtx has */ + /* nodata values like -2147479936), so test them too */ + if( value > 1000 || value < -1000 || value == -88.88880f ) + value = HUGE_VAL; + + + return value; +} /************************************************************************/ /* pj_apply_vgridshift() */ /* */ -/* This implementation takes uses the gridlist from a coordinate */ +/* This implementation 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, + PJ_GRIDINFO ***gridlist_p, int *gridlist_count_p, - int inverse, + int inverse, long point_count, int point_offset, double *x, double *y, double *z ) @@ -52,22 +149,23 @@ int pj_apply_vgridshift( PJ *defn, const char *listname, int i; static int debug_count = 0; PJ_GRIDINFO **tables; + struct CTABLE ct; if( *gridlist_p == NULL ) { - *gridlist_p = - pj_gridlist_from_nadgrids( pj_get_ctx(defn), + *gridlist_p = + pj_gridlist_from_nadgrids( pj_get_ctx(defn), pj_param(defn->ctx,defn->params,listname).s, gridlist_count_p ); if( *gridlist_p == NULL || *gridlist_count_p == 0 ) return defn->ctx->last_errno; } - + if( *gridlist_count_p == 0 ) { - pj_ctx_set_errno( defn->ctx, -38); - return -38; + pj_ctx_set_errno( defn->ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return PJD_ERR_FAILED_TO_LOAD_GRID; } tables = *gridlist_p; @@ -75,127 +173,37 @@ int pj_apply_vgridshift( PJ *defn, const char *listname, for( i = 0; i < point_count; i++ ) { + double value; long io = i * point_offset; LP input; - int itable = 0; - double value = HUGE_VAL; input.phi = y[io]; input.lam = x[io]; - /* do not deal with NaN coordinates */ - if( input.phi != input.phi || input.lam != input.lam ) - itable = *gridlist_count_p; + value = pj_read_vgrid_value(defn, input, gridlist_count_p, tables, &ct); - /* keep trying till we find a table that works */ - for( ; itable < *gridlist_count_p; itable++ ) + if( inverse ) + z[io] -= value; + else + z[io] += value; + if( value != HUGE_VAL ) { - PJ_GRIDINFO *gi = tables[itable]; - struct CTABLE *ct = gi->ct; - double grid_x, grid_y; - int grid_ix, grid_iy; - int grid_ix2, grid_iy2; - 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. */ - while( 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 didn't find a more refined child node to use, so go with current grid */ - if( child == NULL ) - { - break; - } - - /* Otherwise let's try for childrens children .. */ - gi = child; - ct = child->ct; - } - - /* load the grid shift info if we don't have it. */ - if( ct->cvs == NULL && !pj_gridinfo_load( pj_get_ctx(defn), gi ) ) - { - pj_ctx_set_errno( defn->ctx, -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; - - grid_ix2 = grid_ix + 1; - if( grid_ix2 >= ct->lim.lam ) - grid_ix2 = ct->lim.lam - 1; - grid_iy2 = grid_iy + 1; - if( grid_iy2 >= ct->lim.phi ) - grid_iy2 = ct->lim.phi - 1; - - cvs = (float *) ct->cvs; - value = cvs[grid_ix + grid_iy * ct->lim.lam] - * (1.0-grid_x) * (1.0-grid_y) - + cvs[grid_ix2 + grid_iy * ct->lim.lam] - * (grid_x) * (1.0-grid_y) - + cvs[grid_ix + grid_iy2 * ct->lim.lam] - * (1.0-grid_x) * (grid_y) - + cvs[grid_ix2 + grid_iy2 * ct->lim.lam] - * (grid_x) * (grid_y); - - /* nodata? */ - /* GTX official nodata value if -88.88880f, but some grids also */ - /* use other big values for nodata (e.g naptrans2008.gtx has */ - /* nodata values like -2147479936), so test them too */ - if( value > 1000 || value < -1000 || value == -88.88880f ) - value = HUGE_VAL; - else - { - if( inverse ) - z[io] -= value; - else - z[io] += value; - } - - if( value != HUGE_VAL ) - { - if( debug_count++ < 20 ) - pj_log( defn->ctx, PJ_LOG_DEBUG_MINOR, - "pj_apply_gridshift(): used %s", - ct->id ); + if( debug_count++ < 20 ) { + proj_log_trace(defn, "pj_apply_gridshift(): used %s", ct.id); break; } } if( value == HUGE_VAL ) { + int itable; char gridlist[3000]; - pj_log( defn->ctx, PJ_LOG_DEBUG_MAJOR, - "pj_apply_vgridshift(): failed to find a grid shift table for\n" - " location (%.7fdW,%.7fdN)", - x[io] * RAD_TO_DEG, - y[io] * RAD_TO_DEG ); + proj_log_debug(defn, + "pj_apply_vgridshift(): failed to find a grid shift table for\n" + " location (%.7fdW,%.7fdN)", + x[io] * RAD_TO_DEG, + y[io] * RAD_TO_DEG ); gridlist[0] = '\0'; for( itable = 0; itable < *gridlist_count_p; itable++ ) @@ -212,10 +220,10 @@ int pj_apply_vgridshift( PJ *defn, const char *listname, else sprintf( gridlist+strlen(gridlist), ",%s", gi->gridname ); } - pj_log( defn->ctx, PJ_LOG_DEBUG_MAJOR, - "%s", gridlist ); - + + proj_log_debug(defn, "%s", gridlist); pj_ctx_set_errno( defn->ctx, PJD_ERR_GRID_AREA ); + return PJD_ERR_GRID_AREA; } } @@ -223,3 +231,64 @@ int pj_apply_vgridshift( PJ *defn, const char *listname, return 0; } +/**********************************************/ +int proj_vgrid_init(PJ* P, const char *grids) { +/********************************************** + + 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. + +***********************************************/ + + /* 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->vgridlist_geoid == NULL) { + P->vgridlist_geoid = pj_gridlist_from_nadgrids( + P->ctx, + pj_param(P->ctx, P->params, sgrids).s, + &(P->vgridlist_geoid_count) + ); + + if( P->vgridlist_geoid == NULL || P->vgridlist_geoid_count == 0 ) { + pj_dealloc(sgrids); + return 0; + } + } + + if (P->vgridlist_geoid_count == 0) { + proj_errno_set(P, PJD_ERR_FAILED_TO_LOAD_GRID); + } + + pj_dealloc(sgrids); + return P->vgridlist_geoid_count; +} + +/***********************************************/ +double proj_vgrid_value(PJ *P, LP lp){ +/*********************************************** + + Read grid value at position lp in grids loaded + with proj_grid_init. + + Returns the grid value of the given coordinate. + +************************************************/ + + struct CTABLE used_grid; + double value; + memset(&used_grid, 0, sizeof(struct CTABLE)); + + value = pj_read_vgrid_value(P, lp, &(P->vgridlist_geoid_count), P->vgridlist_geoid, &used_grid); + proj_log_trace(P, "proj_vgrid_value: (%f, %f) = %f", lp.lam*RAD_TO_DEG, lp.phi*RAD_TO_DEG, value); + + return value; +} + diff --git a/src/proj_internal.h b/src/proj_internal.h index 5773637c..29e73e70 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -89,7 +89,12 @@ PJ_OBS pj_invobs (PJ_OBS obs, PJ *P); PJ_COORD pj_fwdcoord (PJ_COORD coo, PJ *P); PJ_COORD pj_invcoord (PJ_COORD coo, PJ *P); - +/* Grid functionality */ +int proj_vgrid_init(PJ *P, const char *grids); +int proj_hgrid_init(PJ *P, const char *grids); +double proj_vgrid_value(PJ *P, LP lp); +LP proj_hgrid_value(PJ *P, LP lp); +LP proj_hgrid_apply(PJ *P, LP lp, PJ_DIRECTION direction); /* High level functionality for handling thread contexts */ enum proj_log_level { |