1 /******************************************************************************
2 * Project: PROJ.4
3 * Purpose: Apply vertical datum shifts based on grid shift files, normally
4 * geoid grids mapping WGS84 to NAVD88 or something similar.
5 * Author: Frank Warmerdam, warmerdam@pobox.com
6 *
7 ******************************************************************************
8 * Copyright (c) 2010, Frank Warmerdam <warmerdam@pobox.com>
9 *
10 * Permission is hereby granted, free of charge, to any person obtaining a
11 * copy of this software and associated documentation files (the "Software"),
12 * to deal in the Software without restriction, including without limitation
13 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14 * and/or sell copies of the Software, and to permit persons to whom the
15 * Software is furnished to do so, subject to the following conditions:
16 *
17 * The above copyright notice and this permission notice shall be included
18 * in all copies or substantial portions of the Software.
19 *
20 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 * DEALINGS IN THE SOFTWARE.
27 *****************************************************************************/
28
29 #define PJ_LIB__
30
31 #include <projects.h>
32 #include <string.h>
33 #include <math.h>
34
35 /************************************************************************/
36 /* pj_apply_vgridshift() */
37 /* */
38 /* This implmentation takes uses the gridlist from a coordinate */
39 /* system definition. If the gridlist has not yet been */
40 /* populated in the coordinate system definition we set it up */
41 /* now. */
42 /************************************************************************/
43
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)44 int pj_apply_vgridshift( PJ *defn, const char *listname,
45 PJ_GRIDINFO ***gridlist_p,
46 int *gridlist_count_p,
47 int inverse,
48 long point_count, int point_offset,
49 double *x, double *y, double *z )
50
51 {
52 int i;
53 static int debug_count = 0;
54 PJ_GRIDINFO **tables;
55
56 if( *gridlist_p == NULL )
57 {
58 *gridlist_p =
59 pj_gridlist_from_nadgrids( pj_get_ctx(defn),
60 pj_param(defn->ctx,defn->params,listname).s,
61 gridlist_count_p );
62
63 if( *gridlist_p == NULL || *gridlist_count_p == 0 )
64 return defn->ctx->last_errno;
65 }
66
67 if( *gridlist_count_p == 0 )
68 {
69 pj_ctx_set_errno( defn->ctx, -38);
70 return -38;
71 }
72
73 tables = *gridlist_p;
74 defn->ctx->last_errno = 0;
75
76 for( i = 0; i < point_count; i++ )
77 {
78 long io = i * point_offset;
79 LP input;
80 int itable;
81 double value = HUGE_VAL;
82
83 input.phi = y[io];
84 input.lam = x[io];
85
86 /* keep trying till we find a table that works */
87 for( itable = 0; itable < *gridlist_count_p; itable++ )
88 {
89 PJ_GRIDINFO *gi = tables[itable];
90 struct CTABLE *ct = gi->ct;
91 double grid_x, grid_y;
92 int grid_ix, grid_iy;
93 float *cvs;
94
95 /* skip tables that don't match our point at all. */
96 if( ct->ll.phi > input.phi || ct->ll.lam > input.lam
97 || ct->ll.phi + (ct->lim.phi-1) * ct->del.phi < input.phi
98 || ct->ll.lam + (ct->lim.lam-1) * ct->del.lam < input.lam )
99 continue;
100
101 /* If we have child nodes, check to see if any of them apply. */
102 while( gi->child != NULL )
103 {
104 PJ_GRIDINFO *child;
105
106 for( child = gi->child; child != NULL; child = child->next )
107 {
108 struct CTABLE *ct1 = child->ct;
109
110 if( ct1->ll.phi > input.phi || ct1->ll.lam > input.lam
111 || ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi < input.phi
112 || ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam < input.lam)
113 continue;
114
115 break;
116 }
117
118 /* we didn't find a more refined child node to use, so go with current grid */
119 if( child == NULL )
120 {
121 break;
122 }
123
124 /* Otherwise let's try for childrens children .. */
125 gi = child;
126 ct = child->ct;
127 }
128
129 /* load the grid shift info if we don't have it. */
130 if( ct->cvs == NULL && !pj_gridinfo_load( pj_get_ctx(defn), gi ) )
131 {
132 pj_ctx_set_errno( defn->ctx, -38 );
133 return -38;
134 }
135
136 /* Interpolation a location within the grid */
137 grid_x = (input.lam - ct->ll.lam) / ct->del.lam;
138 grid_y = (input.phi - ct->ll.phi) / ct->del.phi;
139 grid_ix = (int) floor(grid_x);
140 grid_iy = (int) floor(grid_y);
141 grid_x -= grid_ix;
142 grid_y -= grid_iy;
143
144 cvs = (float *) ct->cvs;
145 value = cvs[grid_ix + grid_iy * ct->lim.lam]
146 * (1.0-grid_x) * (1.0-grid_y)
147 + cvs[grid_ix + 1 + grid_iy * ct->lim.lam]
148 * (grid_x) * (1.0-grid_y)
149 + cvs[grid_ix + (grid_iy+1) * ct->lim.lam]
150 * (1.0-grid_x) * (grid_y)
151 + cvs[grid_ix + 1 + (grid_iy+1) * ct->lim.lam]
152 * (grid_x) * (grid_y);
153
154 /* nodata? */
155 /* GTX official nodata value if -88.88880f, but some grids also */
156 /* use other big values for nodata (e.g naptrans2008.gtx has */
157 /* nodata values like -2147479936), so test them too */
158 if( value > 1000 || value < -1000 || value == -88.88880f )
159 value = HUGE_VAL;
160 else
161 {
162 if( inverse )
163 z[io] -= value;
164 else
165 z[io] += value;
166 }
167
168 if( value != HUGE_VAL )
169 {
170 if( debug_count++ < 20 )
171 pj_log( defn->ctx, PJ_LOG_DEBUG_MINOR,
172 "pj_apply_gridshift(): used %s",
173 ct->id );
174 break;
175 }
176 }
177
178 if( value == HUGE_VAL )
179 {
180 char gridlist[3000];
181
182 pj_log( defn->ctx, PJ_LOG_DEBUG_MAJOR,
183 "pj_apply_vgridshift(): failed to find a grid shift table for\n"
184 " location (%.7fdW,%.7fdN)",
185 x[io] * RAD_TO_DEG,
186 y[io] * RAD_TO_DEG );
187
188 gridlist[0] = '\0';
189 for( itable = 0; itable < *gridlist_count_p; itable++ )
190 {
191 PJ_GRIDINFO *gi = tables[itable];
192 if( strlen(gridlist) + strlen(gi->gridname) > sizeof(gridlist)-100 )
193 {
194 strcat( gridlist, "..." );
195 break;
196 }
197
198 if( itable == 0 )
199 sprintf( gridlist, " tried: %s", gi->gridname );
200 else
201 sprintf( gridlist+strlen(gridlist), ",%s", gi->gridname );
202 }
203 pj_log( defn->ctx, PJ_LOG_DEBUG_MAJOR,
204 "%s", gridlist );
205
206 pj_ctx_set_errno( defn->ctx, PJD_ERR_GRID_AREA );
207 return PJD_ERR_GRID_AREA;
208 }
209 }
210
211 return 0;
212 }
213
214