pj_gridinfo.c

开源投影系统 Cartographic Projections library originally written by Gerald Evenden then of the USGS. The

        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;

        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;

        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 initialize our grid info.                      */
/* -------------------------------------------------------------------- */
        if( subfile == 0 )
            gi = gilist;
        else
        {
            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;
        }

        gi->ct = ct;
        gi->format = "ntv2";
        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 );
    }

    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;

    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;

    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[160];

    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_ctable_init( fp );

        gilist->format = "ctable";
        gilist->ct = ct;

        if( getenv("PROJ_DEBUG") != NULL )
            fprintf( stderr, 
                     "Ctable %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)\n",
                     ct->id, 
                     ct->lim.lam, ct->lim.phi,
                     ct->ll.lam * RAD_TO_DEG, ct->ll.phi * RAD_TO_DEG,
                     (ct->ll.lam + ct->lim.lam*ct->del.lam) * RAD_TO_DEG, 
                     (ct->ll.phi + ct->lim.phi*ct->del.lam) * RAD_TO_DEG );
    }

    fclose(fp);

    return gilist;
}