geogrid/src/read_geogrid.c

   1 /* File: read_geogrid.c
   2
   3    Sample subroutine to read an array from the geogrid binary format.
   4
   5    Notes: Depending on the compiler and compiler flags, the name of
   6    the read_geogrid() routine may need to be adjusted with respect
   7    to the number of trailing underscores when calling from Fortran.
   8
   9    Michael G. Duda, NCAR/MMM
  10 */
  11
  12 #include <stdlib.h>
  13 #include <stdio.h>
  14 #include <string.h>
  15
  16 #ifdef _UNDERSCORE
  17 #define read_geogrid read_geogrid_
  18 #endif
  19 #ifdef _DOUBLEUNDERSCORE
  20 #define read_geogrid read_geogrid__
  21 #endif
  22
  23 #define BIG_ENDIAN    0
  24 #define LITTLE_ENDIAN 1
  25
  26 int read_geogrid(
  27       char * fname,            /* The name of the file to read from */
  28       int * len,               /* The length of the filename */
  29       float * rarray,          /* The array to be filled */
  30       int * nx,                /* x-dimension of the array */
  31       int * ny,                /* y-dimension of the array */
  32       int * nz,                /* z-dimension of the array */
  33       int * isigned,           /* 0=unsigned data, 1=signed data */
  34       int * endian,            /* 0=big endian, 1=little endian */
  35       float * scalefactor,     /* value to multiply array elements by before truncation to integers */
  36       int * wordsize,          /* number of bytes to use for each array element */
  37       int * status)
  38 {
  39    int i, ival, cnt, narray;
  40    int A2, B2;
  41    int A3, B3, C3;
  42    int A4, B4, C4, D4;
  43    unsigned char * c;
  44    char local_fname[1024];
  45    FILE * bfile;
  46
  47    *status = 0;
  48
  49    narray = (*nx) * (*ny) * (*nz);
  50
  51    /* Make a null-terminated local copy of the filename */
  52    strncpy(local_fname,fname,*len);
  53    local_fname[*len]='\0';
  54
  55    /* Attempt to open file for reading */
  56    if (!(bfile = fopen(local_fname,"rb")))
  57    {
  58       *status = 1;
  59       return 1;
  60    }
  61
  62    /* Allocate memory to hold bytes from file and read data */
  63    c = (unsigned char *)malloc(sizeof(unsigned char)*(*wordsize) * narray);
  64    cnt = fread((void *)c, sizeof(unsigned char), narray*(*wordsize), bfile);
  65
  66    fclose(bfile);
  67
  68    if (cnt == 0)
  69    {
  70       *status = 1;
  71       return 1;
  72    }
  73
  74    /*
  75       Set up byte offsets for each wordsize depending on byte order.
  76       A, B, C, D give the offsets of the LSB through MSB (i.e., for
  77       word ABCD, A=MSB, D=LSB) in the array from the beginning of a word
  78    */
  79    if (*endian == BIG_ENDIAN) {
  80       A2 = 0; B2 = 1;
  81       A3 = 0; B3 = 1; C3 = 2;
  82       A4 = 0; B4 = 1; C4 = 2; D4 = 3;
  83    }
  84    else {
  85       B2 = 0; A2 = 1;
  86       C3 = 0; B3 = 1; A3 = 2;
  87       D4 = 0; C4 = 1; B4 = 2; A4 = 3;
  88    }
  89
  90    /* Convert words from native byte order */
  91    switch(*wordsize) {
  92       case 1:
  93          for(i=0; i<narray; i++)
  94          {
  95             ival = (int)(c[i]);
  96             if ((*isigned) && (ival > (1 << 7))) ival -= (1 << 8);
  97             rarray[i] = (float)ival;
  98          }
  99          break;
 100
 101       case 2:
 102          for(i=0; i<narray; i++)
 103          {
 104             ival = (int)((c[2*i+A2]<<8) | (c[2*i+B2]));
 105             if ((*isigned) && (ival > (1 << 15))) ival -= (1 << 16);
 106             rarray[i] = (float)ival;
 107          }
 108          break;
 109
 110       case 3:
 111          for(i=0; i<narray; i++)
 112          {
 113             ival = (int)((c[3*i+A3]<<16) | (c[3*i+B3]<<8) | c[3*i+C3]);
 114             if ((*isigned) * (ival > (1 << 23))) ival -= (1 << 24);
 115             rarray[i] = (float)ival;
 116          }
 117          break;
 118
 119       case 4:
 120          for(i=0; i<narray; i++)
 121          {
 122             ival = (int)((c[4*i+A4]<<24) | (c[4*i+B4]<<16) | (c[4*i+C4]<<8) | c[4*i+D4]);
 123             if ((*isigned) && (ival > (1 << 31))) ival -= (1 << 32);
 124             rarray[i] = (float)ival;
 125          }
 126          break;
 127    }
 128
 129    free(c);
 130
 131    /* Scale real-valued array by scalefactor */
 132    if (*scalefactor != 1.0)
 133    {
 134       for (i=0; i<narray; i++)
 135          rarray[i] = rarray[i] * (*scalefactor);
 136    }
 137
 138    return 0;
 139 }