implement heirarchical NTv2 gridinfos

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

3 files changed, 151 insertions, 23 deletions

diff --git a/src/pj_apply_gridshift.c b/src/pj_apply_gridshift.c
index 287ffb07..dbfddda7 100644
--- a/src/pj_apply_gridshift.c
+++ b/src/pj_apply_gridshift.c
@@ -31,6 +31,9 @@
  ******************************************************************************
  *
  * $Log$
+ * Revision 1.6  2003/03/17 18:56:34  warmerda
+ * implement heirarchical NTv2 gridinfos
+ *
  * Revision 1.5  2003/03/15 06:02:02  warmerda
  * preliminary NTv2 support, major restructure of datum shifting
  *
@@ -66,6 +69,7 @@ int pj_apply_gridshift( const char *nadgrids, int inverse,
     int grid_count = 0;
     PJ_GRIDINFO   **tables = pj_gridlist_from_nadgrids( nadgrids, &grid_count);
     int  i;
+    int debug_flag = getenv( "PROJ_DEBUG" ) != NULL;
 
     if( tables == NULL || grid_count == 0 )
         return pj_errno;
@@ -82,7 +86,8 @@ int pj_apply_gridshift( const char *nadgrids, int inverse,
         /* keep trying till we find a table that works */
         for( itable = 0; itable < grid_count; itable++ )
         {
-            struct CTABLE *ct = tables[itable]->ct;
+            PJ_GRIDINFO *gi = tables[itable];
+            struct CTABLE *ct = gi->ct;
             
             /* skip tables that don't match our point at all.  */
             if( ct->ll.phi > input.phi || ct->ll.lam > input.lam
@@ -90,22 +95,58 @@ int pj_apply_gridshift( const char *nadgrids, int inverse,
                 || ct->ll.lam + ct->lim.lam * 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*ct1->del.phi < input.phi
+                        || ct1->ll.lam + ct1->lim.lam*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( tables[itable]->ct->cvs == NULL 
-                && !pj_gridinfo_load( tables[itable] ) )
+            if( ct->cvs == NULL && !pj_gridinfo_load( gi ) )
             {
                 pj_errno = -38;
                 return pj_errno;
             }
             
-            output = nad_cvt( input, inverse, tables[itable]->ct );
+            output = nad_cvt( input, inverse, ct );
             if( output.lam != HUGE_VAL )
+            {
+                if( debug_flag )
+                {
+                    static int debug_count = 0;
+
+                    if( debug_count < 20 )
+                        fprintf( stderr,
+                                 "pj_apply_gridshift(): used %s\n",
+                                 ct->id );
+                }
+
                 break;
+            }
         }
 
         if( output.lam == HUGE_VAL )
         {
-            if( getenv( "PROJ_DEBUG" ) != NULL )
+            if( debug_flag )
             {
                 fprintf( stderr, 
                          "pj_apply_gridshift(): failed to find a grid shift table for\n"
diff --git a/src/pj_gridinfo.c b/src/pj_gridinfo.c
index 5db29ab2..8f73e447 100644
--- a/src/pj_gridinfo.c
+++ b/src/pj_gridinfo.c
@@ -29,6 +29,9 @@
  ******************************************************************************
  *
  * $Log$
+ * Revision 1.2  2003/03/17 18:56:34  warmerda
+ * implement heirarchical NTv2 gridinfos
+ *
  * Revision 1.1  2003/03/15 06:01:18  warmerda
  * New
  *
@@ -89,6 +92,17 @@ void pj_gridinfo_free( PJ_GRIDINFO *gi )
     if( gi == NULL )
         return;
 
+    if( gi->child != NULL )
+    {
+        PJ_GRIDINFO *child, *next;
+
+        for( child = gi->child; child != NULL; child=next)
+        {
+            next=child->next;
+            pj_gridinfo_free( child );
+        }
+    }
+
     if( gi->ct != NULL )
         nad_free( gi->ct );
     
@@ -118,7 +132,24 @@ int pj_gridinfo_load( PJ_GRIDINFO *gi )
 /*      no real reason not to support delayed loading for it as well.   */
 /* -------------------------------------------------------------------- */
     if( strcmp(gi->format,"ctable") == 0 )
-        return 0;
+    {
+        FILE *fid;
+        int result;
+
+        fid = pj_open_lib( gi->filename, "rb" );
+        
+        if( fid == NULL )
+        {
+            pj_errno = -38;
+            return 0;
+        }
+
+        result = nad_ctable_load( gi->ct, fid );
+
+        fclose( fid );
+
+        return result;
+    }
 
 /* -------------------------------------------------------------------- */
 /*      NTv1 format.                                                    */
@@ -198,10 +229,15 @@ int pj_gridinfo_load( PJ_GRIDINFO *gi )
 /* -------------------------------------------------------------------- */
     else if( strcmp(gi->format,"ntv2") == 0 )
     {
-        double	*row_buf;
+        float	*row_buf;
         int	row;
         FILE *fid;
 
+        if( getenv("PROJ_DEBUG") != NULL )
+        {
+            fprintf( stderr, "NTv2 - loading grid %s\n", gi->ct->id );
+        }
+
         fid = pj_open_lib( gi->filename, "rb" );
         
         if( fid == NULL )
@@ -212,7 +248,7 @@ int pj_gridinfo_load( PJ_GRIDINFO *gi )
 
         fseek( fid, gi->grid_offset, SEEK_SET );
 
-        row_buf = (double *) pj_malloc(gi->ct->lim.lam * sizeof(double) * 4);
+        row_buf = (float *) pj_malloc(gi->ct->lim.lam * sizeof(float) * 4);
         gi->ct->cvs = (FLP *) pj_malloc(gi->ct->lim.lam*gi->ct->lim.phi*sizeof(FLP));
         if( row_buf == NULL || gi->ct->cvs == NULL )
         {
@@ -237,7 +273,7 @@ int pj_gridinfo_load( PJ_GRIDINFO *gi )
             }
 
             if( !IS_LSB )
-                swap_words( (unsigned char *) row_buf, 8, 
+                swap_words( (unsigned char *) row_buf, 4, 
                             gi->ct->lim.lam*4 );
 
             /* convert seconds to radians */
@@ -375,13 +411,19 @@ static int pj_gridinfo_init_ntv2( FILE *fid, PJ_GRIDINFO *gilist )
         ct->lim.lam = (int) (fabs(ur.lam-ct->ll.lam)/ct->del.lam + 0.5) + 1;
         ct->lim.phi = (int) (fabs(ur.phi-ct->ll.phi)/ct->del.phi + 0.5) + 1;
 
+        if( getenv("PROJ_DEBUG") != NULL )
+            fprintf( stderr, 
+                     "NTv2 %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)\n",
+                     ct->id, 
+                     ct->lim.lam, ct->lim.phi,
+                     ct->ll.lam/3600.0, ct->ll.phi/3600.0,
+                     ur.lam/3600.0, ur.phi/3600.0 );
+
         ct->ll.lam *= DEG_TO_RAD/3600.0;
         ct->ll.phi *= DEG_TO_RAD/3600.0;
         ct->del.lam *= DEG_TO_RAD/3600.0;
         ct->del.phi *= DEG_TO_RAD/3600.0;
 
-        printf( "NTv2 Bounds (%s): LL=(%.12g,%.12g)\n",
-                ct->id, ct->ll.lam, ct->ll.phi );
         memcpy( &gs_count, header + 8 + 16*10, 4 );
         if( gs_count != ct->lim.lam * ct->lim.phi )
         {
@@ -397,24 +439,18 @@ static int pj_gridinfo_init_ntv2( FILE *fid, PJ_GRIDINFO *gilist )
 
 /* -------------------------------------------------------------------- */
 /*      Create a new gridinfo for this if we aren't processing the      */
-/*      1st subfile, and add it to the list.                            */
+/*      1st subfile, and initialize our grid info.                      */
 /* -------------------------------------------------------------------- */
         if( subfile == 0 )
             gi = gilist;
         else
         {
-            PJ_GRIDINFO *lnk;
-
             gi = (PJ_GRIDINFO *) pj_malloc(sizeof(PJ_GRIDINFO));
             memset( gi, 0, sizeof(PJ_GRIDINFO) );
     
             gi->gridname = strdup( gilist->gridname );
             gi->filename = strdup( gilist->filename );
             gi->next = NULL;
-
-            for( lnk = gilist; lnk->next != NULL; lnk = lnk->next ) {}
-
-            lnk->next = gi;
         }
 
         gi->ct = ct;
@@ -422,6 +458,49 @@ static int pj_gridinfo_init_ntv2( FILE *fid, PJ_GRIDINFO *gilist )
         gi->grid_offset = ftell( fid );
 
 /* -------------------------------------------------------------------- */
+/*      Attach to the correct list or sublist.                          */
+/* -------------------------------------------------------------------- */
+        if( strncmp(header+24,"NONE",4) == 0 )
+        {
+            if( gi != gilist )
+            {
+                PJ_GRIDINFO *lnk;
+
+                for( lnk = gilist; lnk->next != NULL; lnk = lnk->next ) {}
+                lnk->next = gi;
+            }
+        }
+
+        else
+        {
+            PJ_GRIDINFO *lnk;
+            PJ_GRIDINFO *gp = gilist;
+            
+            while( gp != NULL && strncmp(gp->ct->id,header+24,8) != 0 )
+                gp = gp->next;
+
+            if( gp == NULL )
+            {
+                if( getenv("PROJ_DEBUG") != NULL )
+                    fprintf( stderr, "pj_gridinfo_init_ntv2(): "
+                             "failed to find parent %8.8s for %.\n", 
+                             header+24, gi->ct->id );
+
+                for( lnk = gp; lnk->next != NULL; lnk = lnk->next ) {}
+                lnk->next = gi;
+            }
+            else if( gp->child == NULL )
+            {
+                gp->child = gi;
+            }
+            else
+            {
+                for( lnk = gp->child; lnk->next != NULL; lnk = lnk->next ) {}
+                lnk->next = gi;
+            }
+        }
+
+/* -------------------------------------------------------------------- */
 /*      Seek past the data.                                             */
 /* -------------------------------------------------------------------- */
         fseek( fid, gs_count * 16, SEEK_CUR );
@@ -498,15 +577,18 @@ static int pj_gridinfo_init_ntv1( FILE * fid, PJ_GRIDINFO *gi )
     ct->lim.lam = (int) (fabs(ur.lam-ct->ll.lam)/ct->del.lam + 0.5) + 1;
     ct->lim.phi = (int) (fabs(ur.phi-ct->ll.phi)/ct->del.phi + 0.5) + 1;
 
+    if( getenv("PROJ_DEBUG") != NULL )
+        fprintf( stderr, 
+                 "NTv1 %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)\n",
+                 ct->lim.lam, ct->lim.phi,
+                 ct->ll.lam, ct->ll.phi, ur.lam, ur.phi );
+
     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;
 
-    printf( "NTv1 Bounds: LL=(%.12g,%.12g)\n",
-            ct->ll.lam, ct->ll.phi );
-
     gi->ct = ct;
     gi->grid_offset = ftell( fid );
     gi->format = "ntv1";
@@ -588,7 +670,7 @@ PJ_GRIDINFO *pj_gridinfo_init( const char *gridname )
     
     else
     {
-        struct CTABLE *ct = nad_load_ctable( fp );
+        struct CTABLE *ct = nad_ctable_init( fp );
 
         gilist->format = "ctable";
         gilist->ct = ct;
diff --git a/src/projects.h b/src/projects.h
index f417e4a8..c9bd1672 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -28,6 +28,9 @@
  ******************************************************************************
  *
  * $Log$
+ * Revision 1.17  2003/03/17 18:56:34  warmerda
+ * implement heirarchical NTv2 gridinfos
+ *
  * Revision 1.16  2003/03/15 06:02:02  warmerda
  * preliminary NTv2 support, major restructure of datum shifting
  *
@@ -333,6 +336,7 @@ typedef struct _pj_gi {
     struct CTABLE *ct;
 
     struct _pj_gi *next;
+    struct _pj_gi *child;
 } PJ_GRIDINFO;
 
 /* procedure prototypes */
@@ -390,7 +394,8 @@ int bch2bps(projUV, projUV, projUV **, int, int);
 LP nad_intr(LP, struct CTABLE *);
 LP nad_cvt(LP, int, struct CTABLE *);
 struct CTABLE *nad_init(char *);
-struct CTABLE *nad_load_ctable( FILE * fid );
+struct CTABLE *nad_ctable_init( FILE * fid );
+int nad_ctable_load( struct CTABLE *, FILE * fid );
 void nad_free(struct CTABLE *);
 
 /* higher level handling of datum grid shift files */