proj/src/pj_gridinfo.c

   1 /******************************************************************************
   2  * $Id$
   3  *
   4  * Project:  PROJ.4
   5  * Purpose:  Functions for handling individual PJ_GRIDINFO's.  Includes
   6  *           loaders for all formats but CTABLE (in nad_init.c).
   7  * Author:   Frank Warmerdam, warmerdam@pobox.com
   8  *
   9  ******************************************************************************
  10  * Copyright (c) 2000, Frank Warmerdam <warmerdam@pobox.com>
  11  *
  12  * Permission is hereby granted, free of charge, to any person obtaining a
  13  * copy of this software and associated documentation files (the "Software"),
  14  * to deal in the Software without restriction, including without limitation
  15  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  16  * and/or sell copies of the Software, and to permit persons to whom the
  17  * Software is furnished to do so, subject to the following conditions:
  18  *
  19  * The above copyright notice and this permission notice shall be included
  20  * in all copies or substantial portions of the Software.
  21  *
  22  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  23  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  24  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  25  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  26  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  27  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  28  * DEALINGS IN THE SOFTWARE.
  29  *****************************************************************************/
  30
  31 #define PJ_LIB__
  32
  33 #include <projects.h>
  34 #include <string.h>
  35 #include <math.h>
  36 #include <errno.h>
  37
  38 #ifdef _WIN32_WCE
  39 /* assert.h includes all Windows API headers and causes 'LP' name clash.
  40  * Here assert we disable assert() for Windows CE.
  41  * TODO - mloskot: re-implement porting friendly assert
  42  */
  43 # define assert(exp)    ((void)0)
  44 #else
  45 # include <assert.h>
  46 #endif /* _WIN32_WCE */
  47
  48 /************************************************************************/
  49 /*                             swap_words()                             */
  50 /*                                                                      */
  51 /*      Convert the byte order of the given word(s) in place.           */
  52 /************************************************************************/
  53
  54 static int  byte_order_test = 1;
  55 #define IS_LSB  (((unsigned char *) (&byte_order_test))[0] == 1)
  56
  57 static void swap_words( unsigned char *data, int word_size, int word_count )
  58
  59 {
  60     int word;
  61
  62     for( word = 0; word < word_count; word++ )
  63     {
  64         int     i;
  65
  66         for( i = 0; i < word_size/2; i++ )
  67         {
  68             int t;
  69
  70             t = data[i];
  71             data[i] = data[word_size-i-1];
  72             data[word_size-i-1] = t;
  73         }
  74
  75         data += word_size;
  76     }
  77 }
  78
  79 /************************************************************************/
  80 /*                          pj_gridinfo_free()                          */
  81 /************************************************************************/
  82
  83 void pj_gridinfo_free( projCtx ctx, PJ_GRIDINFO *gi )
  84
  85 {
  86     if( gi == NULL )
  87         return;
  88
  89     if( gi->child != NULL )
  90     {
  91         PJ_GRIDINFO *child, *next;
  92
  93         for( child = gi->child; child != NULL; child=next)
  94         {
  95             next=child->next;
  96             pj_gridinfo_free( ctx, child );
  97         }
  98     }
  99
 100     if( gi->ct != NULL )
 101         nad_free( gi->ct );
 102
 103     free( gi->gridname );
 104     if( gi->filename != NULL )
 105         free( gi->filename );
 106
 107     pj_dalloc( gi );
 108 }
 109
 110 /************************************************************************/
 111 /*                          pj_gridinfo_load()                          */
 112 /*                                                                      */
 113 /*      This function is intended to implement delayed loading of       */
 114 /*      the data contents of a grid file.  The header and related       */
 115 /*      stuff are loaded by pj_gridinfo_init().                         */
 116 /************************************************************************/
 117
 118 int pj_gridinfo_load( projCtx ctx, PJ_GRIDINFO *gi )
 119
 120 {
 121     struct CTABLE ct_tmp;
 122
 123     if( gi == NULL || gi->ct == NULL )
 124         return 0;
 125
 126     pj_acquire_lock();
 127     if( gi->ct->cvs != NULL )
 128     {
 129         pj_release_lock();
 130         return 1;
 131     }
 132
 133     memcpy(&ct_tmp, gi->ct, sizeof(struct CTABLE));
 134
 135 /* -------------------------------------------------------------------- */
 136 /*      Original platform specific CTable format.                       */
 137 /* -------------------------------------------------------------------- */
 138     if( strcmp(gi->format,"ctable") == 0 )
 139     {
 140         PAFile fid;
 141         int result;
 142
 143         fid = pj_open_lib( ctx, gi->filename, "rb" );
 144
 145         if( fid == NULL )
 146         {
 147             pj_ctx_set_errno( ctx, -38 );
 148             pj_release_lock();
 149             return 0;
 150         }
 151
 152         result = nad_ctable_load( ctx, &ct_tmp, fid );
 153
 154         pj_ctx_fclose( ctx, fid );
 155
 156         gi->ct->cvs = ct_tmp.cvs;
 157         pj_release_lock();
 158
 159         return result;
 160     }
 161
 162 /* -------------------------------------------------------------------- */
 163 /*      CTable2 format.                                                 */
 164 /* -------------------------------------------------------------------- */
 165     else if( strcmp(gi->format,"ctable2") == 0 )
 166     {
 167         PAFile fid;
 168         int result;
 169
 170         fid = pj_open_lib( ctx, gi->filename, "rb" );
 171
 172         if( fid == NULL )
 173         {
 174             pj_ctx_set_errno( ctx, -38 );
 175             pj_release_lock();
 176             return 0;
 177         }
 178
 179         result = nad_ctable2_load( ctx, &ct_tmp, fid );
 180
 181         pj_ctx_fclose( ctx, fid );
 182
 183         gi->ct->cvs = ct_tmp.cvs;
 184
 185         pj_release_lock();
 186         return result;
 187     }
 188
 189 /* -------------------------------------------------------------------- */
 190 /*      NTv1 format.                                                    */
 191 /*      We process one line at a time.  Note that the array storage     */
 192 /*      direction (e-w) is different in the NTv1 file and what          */
 193 /*      the CTABLE is supposed to have.  The phi/lam are also           */
 194 /*      reversed, and we have to be aware of byte swapping.             */
 195 /* -------------------------------------------------------------------- */
 196     else if( strcmp(gi->format,"ntv1") == 0 )
 197     {
 198         double  *row_buf;
 199         int     row;
 200         PAFile fid;
 201
 202         fid = pj_open_lib( ctx, gi->filename, "rb" );
 203
 204         if( fid == NULL )
 205         {
 206             pj_ctx_set_errno( ctx, -38 );
 207             pj_release_lock();
 208             return 0;
 209         }
 210
 211         pj_ctx_fseek( ctx, fid, gi->grid_offset, SEEK_SET );
 212
 213         row_buf = (double *) pj_malloc(gi->ct->lim.lam * sizeof(double) * 2);
 214         ct_tmp.cvs = (FLP *) pj_malloc(gi->ct->lim.lam*gi->ct->lim.phi*sizeof(FLP));
 215         if( row_buf == NULL || ct_tmp.cvs == NULL )
 216         {
 217             pj_ctx_set_errno( ctx, -38 );
 218             pj_release_lock();
 219             return 0;
 220         }
 221
 222         for( row = 0; row < gi->ct->lim.phi; row++ )
 223         {
 224             int     i;
 225             FLP     *cvs;
 226             double  *diff_seconds;
 227
 228             if( pj_ctx_fread( ctx, row_buf,
 229                               sizeof(double), gi->ct->lim.lam * 2, fid )
 230                 != 2 * gi->ct->lim.lam )
 231             {
 232                 pj_dalloc( row_buf );
 233                 pj_dalloc( ct_tmp.cvs );
 234                 pj_ctx_set_errno( ctx, -38 );
 235                 return 0;
 236             }
 237
 238             if( IS_LSB )
 239                 swap_words( (unsigned char *) row_buf, 8, gi->ct->lim.lam*2 );
 240
 241             /* convert seconds to radians */
 242             diff_seconds = row_buf;
 243
 244             for( i = 0; i < gi->ct->lim.lam; i++ )
 245             {
 246                 cvs = ct_tmp.cvs + (row) * gi->ct->lim.lam
 247                     + (gi->ct->lim.lam - i - 1);
 248
 249                 cvs->phi = *(diff_seconds++) * ((PI/180.0) / 3600.0);
 250                 cvs->lam = *(diff_seconds++) * ((PI/180.0) / 3600.0);
 251             }
 252         }
 253
 254         pj_dalloc( row_buf );
 255
 256         pj_ctx_fclose( ctx, fid );
 257
 258         gi->ct->cvs = ct_tmp.cvs;
 259         pj_release_lock();
 260
 261         return 1;
 262     }
 263
 264 /* -------------------------------------------------------------------- */
 265 /*      NTv2 format.                                                    */
 266 /*      We process one line at a time.  Note that the array storage     */
 267 /*      direction (e-w) is different in the NTv2 file and what          */
 268 /*      the CTABLE is supposed to have.  The phi/lam are also           */
 269 /*      reversed, and we have to be aware of byte swapping.             */
 270 /* -------------------------------------------------------------------- */
 271     else if( strcmp(gi->format,"ntv2") == 0 )
 272     {
 273         float   *row_buf;
 274         int     row;
 275         PAFile fid;
 276
 277         pj_log( ctx, PJ_LOG_DEBUG_MINOR,
 278                 "NTv2 - loading grid %s", gi->ct->id );
 279
 280         fid = pj_open_lib( ctx, gi->filename, "rb" );
 281
 282         if( fid == NULL )
 283         {
 284             pj_ctx_set_errno( ctx, -38 );
 285             pj_release_lock();
 286             return 0;
 287         }
 288
 289         pj_ctx_fseek( ctx, fid, gi->grid_offset, SEEK_SET );
 290
 291         row_buf = (float *) pj_malloc(gi->ct->lim.lam * sizeof(float) * 4);
 292         ct_tmp.cvs = (FLP *) pj_malloc(gi->ct->lim.lam*gi->ct->lim.phi*sizeof(FLP));
 293         if( row_buf == NULL || ct_tmp.cvs == NULL )
 294         {
 295             pj_ctx_set_errno( ctx, -38 );
 296             pj_release_lock();
 297             return 0;
 298         }
 299
 300         for( row = 0; row < gi->ct->lim.phi; row++ )
 301         {
 302             int     i;
 303             FLP     *cvs;
 304             float   *diff_seconds;
 305
 306             if( pj_ctx_fread( ctx, row_buf, sizeof(float),
 307                               gi->ct->lim.lam*4, fid )
 308                 != 4 * gi->ct->lim.lam )
 309             {
 310                 pj_dalloc( row_buf );
 311                 pj_dalloc( ct_tmp.cvs );
 312                 pj_ctx_set_errno( ctx, -38 );
 313                 pj_release_lock();
 314                 return 0;
 315             }
 316
 317             if( !IS_LSB )
 318                 swap_words( (unsigned char *) row_buf, 4,
 319                             gi->ct->lim.lam*4 );
 320
 321             /* convert seconds to radians */
 322             diff_seconds = row_buf;
 323
 324             for( i = 0; i < gi->ct->lim.lam; i++ )
 325             {
 326                 cvs = ct_tmp.cvs + (row) * gi->ct->lim.lam
 327                     + (gi->ct->lim.lam - i - 1);
 328
 329                 cvs->phi = *(diff_seconds++) * ((PI/180.0) / 3600.0);
 330                 cvs->lam = *(diff_seconds++) * ((PI/180.0) / 3600.0);
 331                 diff_seconds += 2; /* skip accuracy values */
 332             }
 333         }
 334
 335         pj_dalloc( row_buf );
 336
 337         pj_ctx_fclose( ctx, fid );
 338
 339         gi->ct->cvs = ct_tmp.cvs;
 340
 341         pj_release_lock();
 342         return 1;
 343     }
 344
 345 /* -------------------------------------------------------------------- */
 346 /*      GTX format.                                                     */
 347 /* -------------------------------------------------------------------- */
 348     else if( strcmp(gi->format,"gtx") == 0 )
 349     {
 350         int   words = gi->ct->lim.lam * gi->ct->lim.phi;
 351         PAFile fid;
 352
 353         fid = pj_open_lib( ctx, gi->filename, "rb" );
 354
 355         if( fid == NULL )
 356         {
 357             pj_ctx_set_errno( ctx, -38 );
 358             pj_release_lock();
 359             return 0;
 360         }
 361
 362         pj_ctx_fseek( ctx, fid, gi->grid_offset, SEEK_SET );
 363
 364         ct_tmp.cvs = (FLP *) pj_malloc(words*sizeof(float));
 365         if( ct_tmp.cvs == NULL )
 366         {
 367             pj_ctx_set_errno( ctx, -38 );
 368             pj_release_lock();
 369             return 0;
 370         }
 371
 372         if( pj_ctx_fread( ctx, ct_tmp.cvs, sizeof(float), words, fid )
 373             != words )
 374         {
 375             pj_dalloc( ct_tmp.cvs );
 376             pj_release_lock();
 377             return 0;
 378         }
 379
 380         if( IS_LSB )
 381             swap_words( (unsigned char *) ct_tmp.cvs, 4, words );
 382
 383         pj_ctx_fclose( ctx, fid );
 384         gi->ct->cvs = ct_tmp.cvs;
 385         pj_release_lock();
 386         return 1;
 387     }
 388
 389     else
 390     {
 391         pj_release_lock();
 392         return 0;
 393     }
 394 }
 395
 396 /************************************************************************/
 397 /*                       pj_gridinfo_init_ntv2()                        */
 398 /*                                                                      */
 399 /*      Seek a parent grid file by name from a grid list                */
 400 /************************************************************************/
 401
 402 static PJ_GRIDINFO* pj_gridinfo_parent( PJ_GRIDINFO *gilist,
 403         const char *name, int length )
 404 {
 405     while( gilist )
 406     {
 407         if( strncmp(gilist->ct->id,name,length) == 0 ) return gilist;
 408         if( gilist->child )
 409         {
 410             PJ_GRIDINFO *parent=pj_gridinfo_parent( gilist->child, name, length );
 411             if( parent ) return parent;
 412         }
 413         gilist=gilist->next;
 414     }
 415     return gilist;
 416 }
 417
 418 /************************************************************************/
 419 /*                       pj_gridinfo_init_ntv2()                        */
 420 /*                                                                      */
 421 /*      Load a ntv2 (.gsb) file.                                        */
 422 /************************************************************************/
 423
 424 static int pj_gridinfo_init_ntv2( projCtx ctx, PAFile fid, PJ_GRIDINFO *gilist )
 425 {
 426     unsigned char header[11*16];
 427     int num_subfiles, subfile;
 428
 429     assert( sizeof(int) == 4 );
 430     assert( sizeof(double) == 8 );
 431     if( sizeof(int) != 4 || sizeof(double) != 8 )
 432     {
 433         pj_log( ctx, PJ_LOG_ERROR,
 434              "basic types of inappropraiate size in pj_gridinfo_init_ntv2()" );
 435         pj_ctx_set_errno( ctx, -38 );
 436         return 0;
 437     }
 438
 439 /* -------------------------------------------------------------------- */
 440 /*      Read the overview header.                                       */
 441 /* -------------------------------------------------------------------- */
 442     if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 )
 443     {
 444         pj_ctx_set_errno( ctx, -38 );
 445         return 0;
 446     }
 447
 448 /* -------------------------------------------------------------------- */
 449 /*      Byte swap interesting fields if needed.                         */
 450 /* -------------------------------------------------------------------- */
 451     if( !IS_LSB )
 452     {
 453         swap_words( header+8, 4, 1 );
 454         swap_words( header+8+16, 4, 1 );
 455         swap_words( header+8+32, 4, 1 );
 456         swap_words( header+8+7*16, 8, 1 );
 457         swap_words( header+8+8*16, 8, 1 );
 458         swap_words( header+8+9*16, 8, 1 );
 459         swap_words( header+8+10*16, 8, 1 );
 460     }
 461
 462 /* -------------------------------------------------------------------- */
 463 /*      Get the subfile count out ... all we really use for now.        */
 464 /* -------------------------------------------------------------------- */
 465     memcpy( &num_subfiles, header+8+32, 4 );
 466
 467 /* ==================================================================== */
 468 /*      Step through the subfiles, creating a PJ_GRIDINFO for each.     */
 469 /* ==================================================================== */
 470     for( subfile = 0; subfile < num_subfiles; subfile++ )
 471     {
 472         struct CTABLE *ct;
 473         LP ur;
 474         int gs_count;
 475         PJ_GRIDINFO *gi;
 476
 477 /* -------------------------------------------------------------------- */
 478 /*      Read header.                                                    */
 479 /* -------------------------------------------------------------------- */
 480         if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 )
 481         {
 482             pj_ctx_set_errno( ctx, -38 );
 483             return 0;
 484         }
 485
 486         if( strncmp((const char *) header,"SUB_NAME",8) != 0 )
 487         {
 488             pj_ctx_set_errno( ctx, -38 );
 489             return 0;
 490         }
 491
 492 /* -------------------------------------------------------------------- */
 493 /*      Byte swap interesting fields if needed.                         */
 494 /* -------------------------------------------------------------------- */
 495         if( !IS_LSB )
 496         {
 497             swap_words( header+8+16*4, 8, 1 );
 498             swap_words( header+8+16*5, 8, 1 );
 499             swap_words( header+8+16*6, 8, 1 );
 500             swap_words( header+8+16*7, 8, 1 );
 501             swap_words( header+8+16*8, 8, 1 );
 502             swap_words( header+8+16*9, 8, 1 );
 503             swap_words( header+8+16*10, 4, 1 );
 504         }
 505
 506 /* -------------------------------------------------------------------- */
 507 /*      Initialize a corresponding "ct" structure.                      */
 508 /* -------------------------------------------------------------------- */
 509         ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE));
 510         strncpy( ct->id, (const char *) header + 8, 8 );
 511         ct->id[8] = '\0';
 512
 513         ct->ll.lam = - *((double *) (header+7*16+8)); /* W_LONG */
 514         ct->ll.phi = *((double *) (header+4*16+8));   /* S_LAT */
 515
 516         ur.lam = - *((double *) (header+6*16+8));     /* E_LONG */
 517         ur.phi = *((double *) (header+5*16+8));       /* N_LAT */
 518
 519         ct->del.lam = *((double *) (header+9*16+8));
 520         ct->del.phi = *((double *) (header+8*16+8));
 521
 522         ct->lim.lam = (int) (fabs(ur.lam-ct->ll.lam)/ct->del.lam + 0.5) + 1;
 523         ct->lim.phi = (int) (fabs(ur.phi-ct->ll.phi)/ct->del.phi + 0.5) + 1;
 524
 525         pj_log( ctx, PJ_LOG_DEBUG_MINOR,
 526                 "NTv2 %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)\n",
 527                 ct->id,
 528                 ct->lim.lam, ct->lim.phi,
 529                 ct->ll.lam/3600.0, ct->ll.phi/3600.0,
 530                 ur.lam/3600.0, ur.phi/3600.0 );
 531
 532         ct->ll.lam *= DEG_TO_RAD/3600.0;
 533         ct->ll.phi *= DEG_TO_RAD/3600.0;
 534         ct->del.lam *= DEG_TO_RAD/3600.0;
 535         ct->del.phi *= DEG_TO_RAD/3600.0;
 536
 537         memcpy( &gs_count, header + 8 + 16*10, 4 );
 538         if( gs_count != ct->lim.lam * ct->lim.phi )
 539         {
 540             pj_log( ctx, PJ_LOG_ERROR,
 541                     "GS_COUNT(%d) does not match expected cells (%dx%d=%d)\n",
 542                     gs_count, ct->lim.lam, ct->lim.phi,
 543                     ct->lim.lam * ct->lim.phi );
 544             pj_ctx_set_errno( ctx, -38 );
 545             return 0;
 546         }
 547
 548         ct->cvs = NULL;
 549
 550 /* -------------------------------------------------------------------- */
 551 /*      Create a new gridinfo for this if we aren't processing the      */
 552 /*      1st subfile, and initialize our grid info.                      */
 553 /* -------------------------------------------------------------------- */
 554         if( subfile == 0 )
 555             gi = gilist;
 556         else
 557         {
 558             gi = (PJ_GRIDINFO *) pj_malloc(sizeof(PJ_GRIDINFO));
 559             memset( gi, 0, sizeof(PJ_GRIDINFO) );
 560
 561             gi->gridname = strdup( gilist->gridname );
 562             gi->filename = strdup( gilist->filename );
 563             gi->next = NULL;
 564         }
 565
 566         gi->ct = ct;
 567         gi->format = "ntv2";
 568         gi->grid_offset = pj_ctx_ftell( ctx, fid );
 569
 570 /* -------------------------------------------------------------------- */
 571 /*      Attach to the correct list or sublist.                          */
 572 /* -------------------------------------------------------------------- */
 573         if( strncmp((const char *)header+24,"NONE",4) == 0 )
 574         {
 575             if( gi != gilist )
 576             {
 577                 PJ_GRIDINFO *lnk;
 578
 579                 for( lnk = gilist; lnk->next != NULL; lnk = lnk->next ) {}
 580                 lnk->next = gi;
 581             }
 582         }
 583
 584         else
 585         {
 586             PJ_GRIDINFO *lnk;
 587             PJ_GRIDINFO *gp = pj_gridinfo_parent(gilist,
 588                                                  (const char*)header+24,8);
 589
 590             if( gp == NULL )
 591             {
 592                 pj_log( ctx, PJ_LOG_ERROR,
 593                         "pj_gridinfo_init_ntv2(): "
 594                         "failed to find parent %8.8s for %s.\n",
 595                         (const char *) header+24, gi->ct->id );
 596
 597                 for( lnk = gilist; lnk->next != NULL; lnk = lnk->next ) {}
 598                 lnk->next = gi;
 599             }
 600             else
 601             {
 602                 if( gp->child == NULL )
 603                 {
 604                     gp->child = gi;
 605                 }
 606                 else
 607                 {
 608                     for( lnk = gp->child; lnk->next != NULL; lnk = lnk->next ) {}
 609                     lnk->next = gi;
 610                 }
 611             }
 612         }
 613
 614 /* -------------------------------------------------------------------- */
 615 /*      Seek past the data.                                             */
 616 /* -------------------------------------------------------------------- */
 617         pj_ctx_fseek( ctx, fid, gs_count * 16, SEEK_CUR );
 618     }
 619
 620     return 1;
 621 }
 622
 623 /************************************************************************/
 624 /*                       pj_gridinfo_init_ntv1()                        */
 625 /*                                                                      */
 626 /*      Load an NTv1 style Canadian grid shift file.                    */
 627 /************************************************************************/
 628
 629 static int pj_gridinfo_init_ntv1( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi )
 630
 631 {
 632     unsigned char header[176];
 633     struct CTABLE *ct;
 634     LP          ur;
 635
 636     assert( sizeof(int) == 4 );
 637     assert( sizeof(double) == 8 );
 638     if( sizeof(int) != 4 || sizeof(double) != 8 )
 639     {
 640         pj_log( ctx, PJ_LOG_ERROR,
 641                  "basic types of inappropraiate size in nad_load_ntv1()" );
 642         pj_ctx_set_errno( ctx, -38 );
 643         return 0;
 644     }
 645
 646 /* -------------------------------------------------------------------- */
 647 /*      Read the header.                                                */
 648 /* -------------------------------------------------------------------- */
 649     if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 )
 650     {
 651         pj_ctx_set_errno( ctx, -38 );
 652         return 0;
 653     }
 654
 655 /* -------------------------------------------------------------------- */
 656 /*      Regularize fields of interest.                                  */
 657 /* -------------------------------------------------------------------- */
 658     if( IS_LSB )
 659     {
 660         swap_words( header+8, 4, 1 );
 661         swap_words( header+24, 8, 1 );
 662         swap_words( header+40, 8, 1 );
 663         swap_words( header+56, 8, 1 );
 664         swap_words( header+72, 8, 1 );
 665         swap_words( header+88, 8, 1 );
 666         swap_words( header+104, 8, 1 );
 667     }
 668
 669     if( *((int *) (header+8)) != 12 )
 670     {
 671         pj_log( ctx, PJ_LOG_ERROR,
 672                 "NTv1 grid shift file has wrong record count, corrupt?" );
 673         pj_ctx_set_errno( ctx, -38 );
 674         return 0;
 675     }
 676
 677 /* -------------------------------------------------------------------- */
 678 /*      Fill in CTABLE structure.                                       */
 679 /* -------------------------------------------------------------------- */
 680     ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE));
 681     strcpy( ct->id, "NTv1 Grid Shift File" );
 682
 683     ct->ll.lam = - *((double *) (header+72));
 684     ct->ll.phi = *((double *) (header+24));
 685     ur.lam = - *((double *) (header+56));
 686     ur.phi = *((double *) (header+40));
 687     ct->del.lam = *((double *) (header+104));
 688     ct->del.phi = *((double *) (header+88));
 689     ct->lim.lam = (int) (fabs(ur.lam-ct->ll.lam)/ct->del.lam + 0.5) + 1;
 690     ct->lim.phi = (int) (fabs(ur.phi-ct->ll.phi)/ct->del.phi + 0.5) + 1;
 691
 692     pj_log( ctx, PJ_LOG_DEBUG_MINOR,
 693             "NTv1 %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)",
 694             ct->lim.lam, ct->lim.phi,
 695             ct->ll.lam, ct->ll.phi, ur.lam, ur.phi );
 696
 697     ct->ll.lam *= DEG_TO_RAD;
 698     ct->ll.phi *= DEG_TO_RAD;
 699     ct->del.lam *= DEG_TO_RAD;
 700     ct->del.phi *= DEG_TO_RAD;
 701     ct->cvs = NULL;
 702
 703     gi->ct = ct;
 704     gi->grid_offset = pj_ctx_ftell( ctx, fid );
 705     gi->format = "ntv1";
 706
 707     return 1;
 708 }
 709
 710 /************************************************************************/
 711 /*                       pj_gridinfo_init_gtx()                         */
 712 /*                                                                      */
 713 /*      Load a NOAA .gtx vertical datum shift file.                     */
 714 /************************************************************************/
 715
 716 static int pj_gridinfo_init_gtx( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi )
 717
 718 {
 719     unsigned char header[40];
 720     struct CTABLE *ct;
 721     double      xorigin,yorigin,xstep,ystep;
 722     int         rows, columns;
 723
 724     assert( sizeof(int) == 4 );
 725     assert( sizeof(double) == 8 );
 726     if( sizeof(int) != 4 || sizeof(double) != 8 )
 727     {
 728         pj_log( ctx, PJ_LOG_ERROR,
 729                 "basic types of inappropraiate size in nad_load_gtx()" );
 730         pj_ctx_set_errno( ctx, -38 );
 731         return 0;
 732     }
 733
 734 /* -------------------------------------------------------------------- */
 735 /*      Read the header.                                                */
 736 /* -------------------------------------------------------------------- */
 737     if( pj_ctx_fread( ctx, header, sizeof(header), 1, fid ) != 1 )
 738     {
 739         pj_ctx_set_errno( ctx, -38 );
 740         return 0;
 741     }
 742
 743 /* -------------------------------------------------------------------- */
 744 /*      Regularize fields of interest and extract.                      */
 745 /* -------------------------------------------------------------------- */
 746     if( IS_LSB )
 747     {
 748         swap_words( header+0, 8, 4 );
 749         swap_words( header+32, 4, 2 );
 750     }
 751
 752     memcpy( &yorigin, header+0, 8 );
 753     memcpy( &xorigin, header+8, 8 );
 754     memcpy( &ystep, header+16, 8 );
 755     memcpy( &xstep, header+24, 8 );
 756
 757     memcpy( &rows, header+32, 4 );
 758     memcpy( &columns, header+36, 4 );
 759
 760     if( xorigin < -360 || xorigin > 360
 761         || yorigin < -90 || yorigin > 90 )
 762     {
 763         pj_log( ctx, PJ_LOG_ERROR,
 764                 "gtx file header has invalid extents, corrupt?");
 765         pj_ctx_set_errno( ctx, -38 );
 766         return 0;
 767     }
 768
 769 /* -------------------------------------------------------------------- */
 770 /*      Fill in CTABLE structure.                                       */
 771 /* -------------------------------------------------------------------- */
 772     ct = (struct CTABLE *) pj_malloc(sizeof(struct CTABLE));
 773     strcpy( ct->id, "GTX Vertical Grid Shift File" );
 774
 775     ct->ll.lam = xorigin;
 776     ct->ll.phi = yorigin;
 777     ct->del.lam = xstep;
 778     ct->del.phi = ystep;
 779     ct->lim.lam = columns;
 780     ct->lim.phi = rows;
 781
 782     /* some GTX files come in 0-360 and we shift them back into the
 783        expected -180 to 180 range if possible.  This does not solve
 784        problems with grids spanning the dateline. */
 785     if( ct->ll.lam >= 180.0 )
 786         ct->ll.lam -= 360.0;
 787
 788     if( ct->ll.lam >= 0.0 && ct->ll.lam + ct->del.lam * ct->lim.lam > 180.0 )
 789     {
 790         pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
 791                 "This GTX spans the dateline!  This will cause problems." );
 792     }
 793
 794     pj_log( ctx, PJ_LOG_DEBUG_MINOR,
 795             "GTX %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)",
 796             ct->lim.lam, ct->lim.phi,
 797             ct->ll.lam, ct->ll.phi,
 798             ct->ll.lam + (columns-1)*xstep, ct->ll.phi + (rows-1)*ystep);
 799
 800     ct->ll.lam *= DEG_TO_RAD;
 801     ct->ll.phi *= DEG_TO_RAD;
 802     ct->del.lam *= DEG_TO_RAD;
 803     ct->del.phi *= DEG_TO_RAD;
 804     ct->cvs = NULL;
 805
 806     gi->ct = ct;
 807     gi->grid_offset = 40;
 808     gi->format = "gtx";
 809
 810     return 1;
 811 }
 812
 813 /************************************************************************/
 814 /*                          pj_gridinfo_init()                          */
 815 /*                                                                      */
 816 /*      Open and parse header details from a datum gridshift file       */
 817 /*      returning a list of PJ_GRIDINFOs for the grids in that          */
 818 /*      file.  This superceeds use of nad_init() for modern             */
 819 /*      applications.                                                   */
 820 /************************************************************************/
 821
 822 PJ_GRIDINFO *pj_gridinfo_init( projCtx ctx, const char *gridname )
 823
 824 {
 825     char        fname[MAX_PATH_FILENAME+1];
 826     PJ_GRIDINFO *gilist;
 827     PAFile      fp;
 828     char        header[160];
 829
 830     errno = pj_errno = 0;
 831     ctx->last_errno = 0;
 832
 833 /* -------------------------------------------------------------------- */
 834 /*      Initialize a GRIDINFO with stub info we would use if it         */
 835 /*      cannot be loaded.                                               */
 836 /* -------------------------------------------------------------------- */
 837     gilist = (PJ_GRIDINFO *) pj_malloc(sizeof(PJ_GRIDINFO));
 838     memset( gilist, 0, sizeof(PJ_GRIDINFO) );
 839
 840     gilist->gridname = strdup( gridname );
 841     gilist->filename = NULL;
 842     gilist->format = "missing";
 843     gilist->grid_offset = 0;
 844     gilist->ct = NULL;
 845     gilist->next = NULL;
 846
 847 /* -------------------------------------------------------------------- */
 848 /*      Open the file using the usual search rules.                     */
 849 /* -------------------------------------------------------------------- */
 850     strcpy(fname, gridname);
 851     if (!(fp = pj_open_lib(ctx, fname, "rb"))) {
 852         ctx->last_errno = 0; /* don't treat as a persistent error */
 853         return gilist;
 854     }
 855
 856     gilist->filename = strdup(fname);
 857
 858 /* -------------------------------------------------------------------- */
 859 /*      Load a header, to determine the file type.                      */
 860 /* -------------------------------------------------------------------- */
 861     if( pj_ctx_fread( ctx, header, sizeof(header), 1, fp ) != 1 )
 862     {
 863         pj_ctx_fclose( ctx, fp );
 864         pj_ctx_set_errno( ctx, -38 );
 865         return gilist;
 866     }
 867
 868     pj_ctx_fseek( ctx, fp, SEEK_SET, 0 );
 869
 870 /* -------------------------------------------------------------------- */
 871 /*      Determine file type.                                            */
 872 /* -------------------------------------------------------------------- */
 873     if( strncmp(header + 0, "HEADER", 6) == 0
 874         && strncmp(header + 96, "W GRID", 6) == 0
 875         && strncmp(header + 144, "TO      NAD83   ", 16) == 0 )
 876     {
 877         pj_gridinfo_init_ntv1( ctx, fp, gilist );
 878     }
 879
 880     else if( strncmp(header + 0, "NUM_OREC", 8) == 0
 881              && strncmp(header + 48, "GS_TYPE", 7) == 0 )
 882     {
 883         pj_gridinfo_init_ntv2( ctx, fp, gilist );
 884     }
 885
 886     else if( strlen(gridname) > 4
 887              && (strcmp(gridname+strlen(gridname)-3,"gtx") == 0
 888                  || strcmp(gridname+strlen(gridname)-3,"GTX") == 0) )
 889     {
 890         pj_gridinfo_init_gtx( ctx, fp, gilist );
 891     }
 892
 893     else if( strncmp(header + 0,"CTABLE V2",9) == 0 )
 894     {
 895         struct CTABLE *ct = nad_ctable2_init( ctx, fp );
 896
 897         gilist->format = "ctable2";
 898         gilist->ct = ct;
 899
 900         pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
 901                 "Ctable2 %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)\n",
 902                 ct->id,
 903                 ct->lim.lam, ct->lim.phi,
 904                 ct->ll.lam * RAD_TO_DEG, ct->ll.phi * RAD_TO_DEG,
 905                 (ct->ll.lam + (ct->lim.lam-1)*ct->del.lam) * RAD_TO_DEG,
 906                 (ct->ll.phi + (ct->lim.phi-1)*ct->del.phi) * RAD_TO_DEG );
 907     }
 908
 909     else
 910     {
 911         struct CTABLE *ct = nad_ctable_init( ctx, fp );
 912         if (ct == NULL)
 913         {
 914             pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
 915                     "CTABLE ct is NULL.");
 916         } else
 917         {
 918             gilist->format = "ctable";
 919             gilist->ct = ct;
 920
 921             pj_log( ctx, PJ_LOG_DEBUG_MAJOR,
 922                     "Ctable %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)\n",
 923                     ct->id,
 924                     ct->lim.lam, ct->lim.phi,
 925                     ct->ll.lam * RAD_TO_DEG, ct->ll.phi * RAD_TO_DEG,
 926                     (ct->ll.lam + (ct->lim.lam-1)*ct->del.lam) * RAD_TO_DEG,
 927                     (ct->ll.phi + (ct->lim.phi-1)*ct->del.phi) * RAD_TO_DEG );
 928         }
 929     }
 930
 931     pj_ctx_fclose(ctx, fp);
 932
 933     return gilist;
 934 }