New

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

2 files changed, 841 insertions, 0 deletions

diff --git a/src/pj_gridinfo.c b/src/pj_gridinfo.c
new file mode 100644
index 00000000..5db29ab2
--- /dev/null
+++ b/src/pj_gridinfo.c
@@ -0,0 +1,600 @@
+/******************************************************************************
+ * $Id$
+ *
+ * Project:  PROJ.4
+ * Purpose:  Functions for handling individual PJ_GRIDINFO's.  Includes
+ *           loaders for all formats but CTABLE (in nad_init.c).
+ * Author:   Frank Warmerdam, warmerdam@pobox.com
+ *
+ ******************************************************************************
+ * Copyright (c) 2000, 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.
+ ******************************************************************************
+ *
+ * $Log$
+ * Revision 1.1  2003/03/15 06:01:18  warmerda
+ * New
+ *
+ */
+
+#define PJ_LIB__
+
+#include <projects.h>
+#include <string.h>
+#include <math.h>
+#include <errno.h>
+#include <assert.h>
+
+/************************************************************************/
+/*                             swap_words()                             */
+/*                                                                      */
+/*      Convert the byte order of the given word(s) in place.           */
+/************************************************************************/
+
+static int  byte_order_test = 1;
+#define IS_LSB	(((unsigned char *) (&byte_order_test))[0] == 1)
+
+static void swap_words( unsigned char *data, int word_size, int word_count )
+
+{
+    /* We only need to do work on LSB machines.  Perhaps we should 
+       convert the data files into LSB order to cut workload! */
+
+    if( IS_LSB )
+    {
+        int	word;
+
+        for( word = 0; word < word_count; word++ )
+        {
+            int	i;
+
+            for( i = 0; i < word_size/2; i++ )
+            {
+                int	t;
+
+                t = data[i];
+                data[i] = data[word_size-i-1];
+                data[word_size-i-1] = t;
+            }
+
+            data += word_size;
+        }
+    }
+}
+
+/************************************************************************/
+/*                          pj_gridinfo_free()                          */
+/************************************************************************/
+
+void pj_gridinfo_free( PJ_GRIDINFO *gi )
+
+{
+    if( gi == NULL )
+        return;
+
+    if( gi->ct != NULL )
+        nad_free( gi->ct );
+    
+    free( gi->gridname );
+    if( gi->filename != NULL )
+        free( gi->filename );
+
+    pj_dalloc( gi );
+}
+
+/************************************************************************/
+/*                          pj_gridinfo_load()                          */
+/*                                                                      */
+/*      This function is intended to implement delayed loading of       */
+/*      the data contents of a grid file.  The header and related       */
+/*      stuff are loaded by pj_gridinfo_init().                         */
+/************************************************************************/
+
+int pj_gridinfo_load( PJ_GRIDINFO *gi )
+
+{
+    if( gi == NULL || gi->ct == NULL )
+        return 0;
+
+/* -------------------------------------------------------------------- */
+/*      ctable is currently loaded on initialization though there is    */
+/*      no real reason not to support delayed loading for it as well.   */
+/* -------------------------------------------------------------------- */
+    if( strcmp(gi->format,"ctable") == 0 )
+        return 0;
+
+/* -------------------------------------------------------------------- */
+/*      NTv1 format.                                                    */
+/*      We process one line at a time.  Note that the array storage     */
+/*      direction (e-w) is different in the NTv1 file and what          */
+/*      the CTABLE is supposed to have.  The phi/lam are also           */
+/*      reversed, and we have to be aware of byte swapping.             */
+/* -------------------------------------------------------------------- */
+    else if( strcmp(gi->format,"ntv1") == 0 )
+    {
+        double	*row_buf;
+        int	row;
+        FILE *fid;
+
+        fid = pj_open_lib( gi->filename, "rb" );
+        
+        if( fid == NULL )
+        {
+            pj_errno = -38;
+            return 0;
+        }
+
+        fseek( fid, gi->grid_offset, SEEK_SET );
+
+        row_buf = (double *) pj_malloc(gi->ct->lim.lam * sizeof(double) * 2);
+        gi->ct->cvs = (FLP *) pj_malloc(gi->ct->lim.lam*gi->ct->lim.phi*sizeof(FLP));
+        if( row_buf == NULL || gi->ct->cvs == NULL )
+        {
+            pj_errno = -38;
+            return 0;
+        }
+        
+        for( row = 0; row < gi->ct->lim.phi; row++ )
+        {
+            int	    i;
+            FLP     *cvs;
+            double  *diff_seconds;
+
+            if( fread( row_buf, sizeof(double), gi->ct->lim.lam * 2, fid ) 
+                != 2 * gi->ct->lim.lam )
+            {
+                pj_dalloc( row_buf );
+                pj_dalloc( gi->ct->cvs );
+                pj_errno = -38;
+                return 0;
+            }
+
+            if( IS_LSB )
+                swap_words( (unsigned char *) row_buf, 8, gi->ct->lim.lam*2 );
+
+            /* convert seconds to radians */
+            diff_seconds = row_buf;
+
+            for( i = 0; i < gi->ct->lim.lam; i++ )
+            {
+                cvs = gi->ct->cvs + (row) * gi->ct->lim.lam
+                    + (gi->ct->lim.lam - i - 1);
+
+                cvs->phi = *(diff_seconds++) * ((PI/180.0) / 3600.0);
+                cvs->lam = *(diff_seconds++) * ((PI/180.0) / 3600.0);
+            }
+        }
+
+        pj_dalloc( row_buf );
+
+        fclose( fid );
+
+        return 1;
+    }
+
+/* -------------------------------------------------------------------- */
+/*      NTv2 format.                                                    */
+/*      We process one line at a time.  Note that the array storage     */
+/*      direction (e-w) is different in the NTv2 file and what          */
+/*      the CTABLE is supposed to have.  The phi/lam are also           */
+/*      reversed, and we have to be aware of byte swapping.             */
+/* -------------------------------------------------------------------- */
+    else if( strcmp(gi->format,"ntv2") == 0 )
+    {
+        double	*row_buf;
+        int	row;
+        FILE *fid;
+
+        fid = pj_open_lib( gi->filename, "rb" );
+        
+        if( fid == NULL )
+        {
+            pj_errno = -38;
+            return 0;
+        }
+
+        fseek( fid, gi->grid_offset, SEEK_SET );
+
+        row_buf = (double *) pj_malloc(gi->ct->lim.lam * sizeof(double) * 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 )
+        {
+            pj_errno = -38;
+            return 0;
+        }
+        
+        for( row = 0; row < gi->ct->lim.phi; row++ )
+        {
+            int	    i;
+            FLP     *cvs;
+            float   *diff_seconds;
+
+            if( fread( row_buf, sizeof(float), gi->ct->lim.lam*4, fid ) 
+                != 4 * gi->ct->lim.lam )
+            {
+                pj_dalloc( row_buf );
+                pj_dalloc( gi->ct->cvs );
+                gi->ct->cvs = NULL;
+                pj_errno = -38;
+                return 0;
+            }
+
+            if( !IS_LSB )
+                swap_words( (unsigned char *) row_buf, 8, 
+                            gi->ct->lim.lam*4 );
+
+            /* convert seconds to radians */
+            diff_seconds = row_buf;
+
+            for( i = 0; i < gi->ct->lim.lam; i++ )
+            {
+                cvs = gi->ct->cvs + (row) * gi->ct->lim.lam
+                    + (gi->ct->lim.lam - i - 1);
+
+                cvs->phi = *(diff_seconds++) * ((PI/180.0) / 3600.0);
+                cvs->lam = *(diff_seconds++) * ((PI/180.0) / 3600.0);
+                diff_seconds += 2; /* skip accuracy values */
+            }
+        }
+
+        pj_dalloc( row_buf );
+
+        fclose( fid );
+
+        return 1;
+    }
+
+    else
+    {
+        return 0;
+    }
+}
+
+/************************************************************************/
+/*                       pj_gridinfo_init_ntv2()                        */
+/*                                                                      */
+/*      Load a ntv2 (.gsb) file.                                        */
+/************************************************************************/
+
+static int pj_gridinfo_init_ntv2( FILE *fid, PJ_GRIDINFO *gilist )
+
+{
+    unsigned char header[11*16];
+    int num_subfiles, subfile;
+
+    assert( sizeof(int) == 4 );
+    assert( sizeof(double) == 8 );
+    if( sizeof(int) != 4 || sizeof(double) != 8 )
+    {
+        fprintf( stderr, 
+                 "basic types of inappropraiate size in pj_gridinfo_init_ntv2()\n" );
+        pj_errno = -38;
+        return 0;
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Read the overview header.                                       */
+/* -------------------------------------------------------------------- */
+    if( fread( header, sizeof(header), 1, fid ) != 1 )
+    {
+        pj_errno = -38;
+        return 0;
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Byte swap interesting fields if needed.                         */
+/* -------------------------------------------------------------------- */
+    if( !IS_LSB )
+    {
+        swap_words( header+8, 4, 1 );
+        swap_words( header+8+16, 4, 1 );
+        swap_words( header+8+32, 4, 1 );
+        swap_words( header+8+7*16, 8, 1 );
+        swap_words( header+8+8*16, 8, 1 );
+        swap_words( header+8+9*16, 8, 1 );
+        swap_words( header+8+10*16, 8, 1 );
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Get the subfile count out ... all we really use for now.        */
+/* -------------------------------------------------------------------- */
+    memcpy( &num_subfiles, header+8+32, 4 );
+
+/* ==================================================================== */
+/*      Step through the subfiles, creating a PJ_GRIDINFO for each.     */
+/* ==================================================================== */
+    for( subfile = 0; subfile < num_subfiles; subfile++ )
+    {
+        struct CTABLE *ct;
+        LP ur;
+        int gs_count;
+        PJ_GRIDINFO *gi;
+
+/* -------------------------------------------------------------------- */
+/*      Read header.                                                    */
+/* -------------------------------------------------------------------- */
+        if( fread( header, sizeof(header), 1, fid ) != 1 )
+        {
+            pj_errno = -38;
+            return 0;
+        }
+
+        if( strncmp(header,"SUB_NAME",8) != 0 )
+        {
+            pj_errno = -38;
+            return 0;
+        }
+        
+/* -------------------------------------------------------------------- */
+/*      Byte swap interesting fields if needed.                         */
+/* -------------------------------------------------------------------- */
+        if( !IS_LSB )
+        {
+            swap_words( header+8+16*4, 8, 1 );
+            swap_words( header+8+16*5, 8, 1 );
+            swap_words( header+8+16*6, 8, 1 );
+            swap_words( header+8+16*7, 8, 1 );
+            swap_words( header+8+16*8, 8, 1 );
+            swap_words( header+8+16*9, 8, 1 );
+            swap_words( header+8+16*10, 4, 1 );
+        }
+        
+/* -------------------------------------------------------------------- */
+/*      Initialize a corresponding "ct" structure.                      */
+/* -------------------------------------------------------------------- */
+        ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE));
+        strncpy( ct->id, header + 8, 8 );
+        ct->id[8] = '\0';
+
+        ct->ll.lam = - *((double *) (header+7*16+8)); /* W_LONG */
+        ct->ll.phi = *((double *) (header+4*16+8));   /* S_LAT */
+
+        ur.lam = - *((double *) (header+6*16+8));     /* E_LONG */
+        ur.phi = *((double *) (header+5*16+8));       /* N_LAT */
+
+        ct->del.lam = *((double *) (header+9*16+8));
+        ct->del.phi = *((double *) (header+8*16+8));
+
+        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;
+
+        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 )
+        {
+            fprintf( stderr, 
+                     "GS_COUNT(%d) does not match expected cells (%dx%d=%d)\n",
+                     gs_count, ct->lim.lam, ct->lim.phi, 
+                     ct->lim.lam * ct->lim.phi );
+            pj_errno = -38;
+            return 0;
+        }
+
+        ct->cvs = NULL;
+
+/* -------------------------------------------------------------------- */
+/*      Create a new gridinfo for this if we aren't processing the      */
+/*      1st subfile, and add it to the list.                            */
+/* -------------------------------------------------------------------- */
+        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;
+        gi->format = "ntv2";
+        gi->grid_offset = ftell( fid );
+
+/* -------------------------------------------------------------------- */
+/*      Seek past the data.                                             */
+/* -------------------------------------------------------------------- */
+        fseek( fid, gs_count * 16, SEEK_CUR );
+    }
+
+    return 1;
+}
+
+/************************************************************************/
+/*                       pj_gridinfo_init_ntv1()                        */
+/*                                                                      */
+/*      Load an NTv1 style Canadian grid shift file.                    */
+/************************************************************************/
+
+static int pj_gridinfo_init_ntv1( FILE * fid, PJ_GRIDINFO *gi )
+
+{
+    unsigned char header[176];
+    struct CTABLE *ct;
+    LP		ur;
+    
+    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_ntv1()\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.                                  */
+/* -------------------------------------------------------------------- */
+    if( IS_LSB )
+    {
+        swap_words( header+8, 4, 1 );
+        swap_words( header+24, 8, 1 );
+        swap_words( header+40, 8, 1 );
+        swap_words( header+56, 8, 1 );
+        swap_words( header+72, 8, 1 );
+        swap_words( header+88, 8, 1 );
+        swap_words( header+104, 8, 1 );
+    }
+
+    if( *((int *) (header+8)) != 12 )
+    {
+        pj_errno = -38;
+        printf("NTv1 grid shift file has wrong record count, corrupt?\n");
+        return 0;
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Fill in CTABLE structure.                                       */
+/* -------------------------------------------------------------------- */
+    ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE));
+    strcpy( ct->id, "NTv1 Grid Shift File" );
+
+    ct->ll.lam = - *((double *) (header+72));
+    ct->ll.phi = *((double *) (header+24));
+    ur.lam = - *((double *) (header+56));
+    ur.phi = *((double *) (header+40));
+    ct->del.lam = *((double *) (header+104));
+    ct->del.phi = *((double *) (header+88));
+    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;
+
+    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";
+
+    return 1;
+}
+
+/************************************************************************/
+/*                          pj_gridinfo_init()                          */
+/*                                                                      */
+/*      Open and parse header details from a datum gridshift file       */
+/*      returning a list of PJ_GRIDINFOs for the grids in that          */
+/*      file.  This superceeds use of nad_init() for modern             */
+/*      applications.                                                   */
+/************************************************************************/
+
+PJ_GRIDINFO *pj_gridinfo_init( const char *gridname )
+
+{
+    char 	fname[MAX_PATH_FILENAME+1];
+    PJ_GRIDINFO *gilist;
+    FILE 	*fp;
+    char	header[512];
+
+    errno = pj_errno = 0;
+
+/* -------------------------------------------------------------------- */
+/*      Initialize a GRIDINFO with stub info we would use if it         */
+/*      cannot be loaded.                                               */
+/* -------------------------------------------------------------------- */
+    gilist = (PJ_GRIDINFO *) pj_malloc(sizeof(PJ_GRIDINFO));
+    memset( gilist, 0, sizeof(PJ_GRIDINFO) );
+    
+    gilist->gridname = strdup( gridname );
+    gilist->filename = NULL;
+    gilist->format = "missing";
+    gilist->grid_offset = 0;
+    gilist->ct = NULL;
+    gilist->next = NULL;
+
+/* -------------------------------------------------------------------- */
+/*      Open the file using the usual search rules.                     */
+/* -------------------------------------------------------------------- */
+    strcpy(fname, gridname);
+    if (!(fp = pj_open_lib(fname, "rb"))) {
+        pj_errno = errno;
+        return gilist;
+    }
+
+    gilist->filename = strdup(fname);
+    
+/* -------------------------------------------------------------------- */
+/*      Load a header, to determine the file type.                      */
+/* -------------------------------------------------------------------- */
+    if( fread( header, sizeof(header), 1, fp ) != 1 )
+    {
+        fclose( fp );
+        pj_errno = -38;
+        return gilist;
+    }
+
+    fseek( fp, SEEK_SET, 0 );
+
+/* -------------------------------------------------------------------- */
+/*      Determine file type.                                            */
+/* -------------------------------------------------------------------- */
+    if( strncmp(header + 0, "HEADER", 6) == 0 
+        && strncmp(header + 96, "W GRID", 6) == 0 
+        && strncmp(header + 144, "TO      NAD83   ", 16) == 0 )
+    {
+        pj_gridinfo_init_ntv1( fp, gilist );
+    }
+    
+    else if( strncmp(header + 0, "NUM_OREC", 8) == 0 
+             && strncmp(header + 48, "GS_TYPE", 7) == 0 )
+    {
+        pj_gridinfo_init_ntv2( fp, gilist );
+    }
+    
+    else
+    {
+        struct CTABLE *ct = nad_load_ctable( fp );
+
+        gilist->format = "ctable";
+        gilist->ct = ct;
+    }
+
+    fclose(fp);
+
+    return gilist;
+}
diff --git a/src/pj_gridlist.c b/src/pj_gridlist.c
new file mode 100644
index 00000000..4e5a1583
--- /dev/null
+++ b/src/pj_gridlist.c
@@ -0,0 +1,241 @@
+/******************************************************************************
+ * $Id$
+ *
+ * Project:  PROJ.4
+ * Purpose:  Code to manage the list of currently loaded (cached) PJ_GRIDINFOs
+ *           See pj_gridinfo.c for details of loading individual grids.
+ * Author:   Frank Warmerdam, warmerdam@pobox.com
+ *
+ ******************************************************************************
+ * Copyright (c) 2000, 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.
+ ******************************************************************************
+ *
+ * $Log$
+ * Revision 1.1  2003/03/15 06:01:18  warmerda
+ * New
+ *
+ */
+
+#define PJ_LIB__
+
+#include <projects.h>
+#include <string.h>
+#include <math.h>
+#include <assert.h>
+
+static PJ_GRIDINFO *grid_list = NULL;
+
+/* used only by pj_load_nadgrids() and pj_deallocate_grids() */
+
+static int           last_nadgrids_max = 0;
+static int           last_nadgrids_count = 0;
+static PJ_GRIDINFO **last_nadgrids_list = NULL;
+static char         *last_nadgrids = NULL;
+
+/************************************************************************/
+/*                        pj_deallocate_grids()                         */
+/*                                                                      */
+/*      Deallocate all loaded grids.                                    */
+/************************************************************************/
+
+void pj_deallocate_grids()
+
+{
+    while( grid_list != NULL )
+    {
+        PJ_GRIDINFO *item = grid_list;
+        grid_list = grid_list->next;
+        item->next = NULL;
+
+        pj_gridinfo_free( item );
+    }
+
+    if( last_nadgrids != NULL )
+    {
+        pj_dalloc( last_nadgrids );
+        last_nadgrids = NULL;
+
+        pj_dalloc( last_nadgrids_list );
+        last_nadgrids_list = NULL;
+
+        last_nadgrids_count = 0;
+        last_nadgrids_max = 0;
+    }
+}
+
+/************************************************************************/
+/*                       pj_gridlist_merge_grid()                       */
+/*                                                                      */
+/*      Find/load the named gridfile and merge it into the              */
+/*      last_nadgrids_list.                                             */
+/************************************************************************/
+
+static int pj_gridlist_merge_gridfile( const char *gridname )
+
+{
+    int i, got_match=0;
+    PJ_GRIDINFO *this_grid, *tail = NULL;
+
+/* -------------------------------------------------------------------- */
+/*      Try to find in the existing list of loaded grids.  Add all      */
+/*      matching grids as with NTv2 we can get many grids from one      */
+/*      file (one shared gridname).                                     */
+/* -------------------------------------------------------------------- */
+    for( this_grid = grid_list; this_grid != NULL; this_grid = this_grid->next)
+    {
+        if( strcmp(this_grid->gridname,gridname) == 0 )
+        {
+            got_match = 1;
+
+            /* dont add to the list if it is invalid. */
+            if( this_grid->ct == NULL )
+                return 0;
+
+            /* do we need to grow the list? */
+            if( last_nadgrids_count >= last_nadgrids_max - 2 )
+            {
+                PJ_GRIDINFO **new_list;
+                int new_max = last_nadgrids_max + 20;
+
+                new_list = (PJ_GRIDINFO **) pj_malloc(sizeof(void*) * new_max);
+                if( last_nadgrids_list != NULL )
+                {
+                    memcpy( new_list, last_nadgrids_list, 
+                            sizeof(void*) * last_nadgrids_max );
+                    pj_dalloc( last_nadgrids_list );
+                }
+
+                last_nadgrids_list = new_list;
+                last_nadgrids_max = new_max;
+            }
+
+            /* add to the list */
+            last_nadgrids_list[last_nadgrids_count++] = this_grid;
+            last_nadgrids_list[last_nadgrids_count] = NULL;
+        }
+
+        tail = this_grid;
+    }
+
+    if( got_match )
+        return 1;
+
+/* -------------------------------------------------------------------- */
+/*      Try to load the named grid.                                     */
+/* -------------------------------------------------------------------- */
+    this_grid = pj_gridinfo_init( gridname );
+
+    if( this_grid == NULL )
+    {
+        /* we should get at least a stub grid with a missing "ct" member */
+        assert( FALSE );
+        return 0;
+    }
+    
+    if( tail != NULL )
+        tail->next = this_grid;
+    else
+        grid_list = this_grid;
+
+/* -------------------------------------------------------------------- */
+/*      Recurse to add the grid now that it is loaded.                  */
+/* -------------------------------------------------------------------- */
+    return pj_gridlist_merge_gridfile( gridname );
+}
+
+/************************************************************************/
+/*                     pj_gridlist_from_nadgrids()                      */
+/*                                                                      */
+/*      This functions loads the list of grids corresponding to a       */
+/*      particular nadgrids string into a list, and returns it.  The    */
+/*      list is kept around till a request is made with a different     */
+/*      string in order to cut down on the string parsing cost, and     */
+/*      the cost of building the list of tables each time.              */
+/************************************************************************/
+
+PJ_GRIDINFO **pj_gridlist_from_nadgrids( const char *nadgrids, int *grid_count)
+
+{
+    const char *s;
+
+    pj_errno = 0;
+    *grid_count = 0;
+
+    if( last_nadgrids != NULL 
+        && strcmp(nadgrids,last_nadgrids) == 0 )
+    {
+        *grid_count = last_nadgrids_count;
+        return last_nadgrids_list;
+    }
+
+/* -------------------------------------------------------------------- */
+/*      Free old one, if any, and make space for new list.              */
+/* -------------------------------------------------------------------- */
+    if( last_nadgrids != NULL )
+    {
+        pj_dalloc(last_nadgrids);
+    }
+    
+    last_nadgrids = (char *) pj_malloc(strlen(nadgrids)+1);
+    strcpy( last_nadgrids, nadgrids );
+
+    last_nadgrids_count = 0;
+
+/* -------------------------------------------------------------------- */
+/*      Loop processing names out of nadgrids one at a time.            */
+/* -------------------------------------------------------------------- */
+    for( s = nadgrids; *s != '\0'; )
+    {
+        int   end_char;
+        char  name[128];
+
+        for( end_char = 0; 
+             s[end_char] != '\0' && s[end_char] != ','; 
+             end_char++ ) {}
+
+        if( end_char > sizeof(name) )
+        {
+            pj_errno = -38;
+            return NULL;
+        }
+        
+        strncpy( name, s, end_char );
+        name[end_char] = '\0';
+
+        s += end_char;
+        if( *s == ',' )
+            s++;
+
+        if( !pj_gridlist_merge_gridfile( name ) )
+        {
+            pj_errno = -38;
+            return NULL;
+        }
+    }
+
+    if( last_nadgrids_count > 0 )
+    {
+        *grid_count = last_nadgrids_count;
+        return last_nadgrids_list;
+    }
+    else
+        return NULL;
+}