pj_gridinfo.cpp

projapi是一个关于GIS行业投影转换的程序库

#include "stdafx.h"

#define PJ_LIB__

#include "projects.h"
#include <string.h>
#ifndef UNICODE
#include <errno.h>
#endif
#include <assert.h>
#include <math.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 )

{
    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->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 );
    
    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 )
    {
        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.                                                    */
/*      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( (int)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 )
    {
        float	*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 = (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 )
        {
            pj_errno = -38;
            return 0;
        }
        
        for( row = 0; row < gi->ct->lim.phi; row++ )
        {
            int	    i;
            FLP     *cvs;
            float   *diff_seconds;

            if( (int)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, 4, 
                            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 );