preliminary addition of vertical datum shifting capability

git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@1839 4e78687f-474d-0410-85f9-8d5e500ac6b2

1 files changed, 200 insertions, 0 deletions

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;
+}
+