source: trunk/WRF.COMMON/WRFV3/external/io_grib1/MEL_grib1/gribgetbds.c @ 3567

#include <stdio.h>
#include <stdlib.h>
#include <math.h> 
#include "dprints.h"            /* for dprints */
#include "gribfuncs.h"          /* prototypes */
/*
  REVISION/MODIFICATION HISTORY:
       03/07/94 written by Mugur Georgescu CSC, Monterey CA
       02/01/96 modified by Steve Lowe SAIC, Monterey CA
       04/17/96 modified by Alice Nakajima SAIC, Monterey CA
       06/19/96 add hdrprint;/nakajima
* 
* ********************************************************************
* A.  FUNCTION:  gribgetbds 
*       decodes the Binary Data Section of the GRIB message 
*       and filling grib_data float array.
*
*    INTERFACE:
*       int gribgetbds (curr_ptr, deci_scale, bms, gds, 
*                       ppgrib_data, bds_head, errmsg)
*
*    ARGUMENTS (I=input, O=output, I&O=input and output):
*      (I)  char *curr_ptr;
*           points to first Octet of the BDS to be decoded;
*      (I)  short  deci_scale;
*           decimal scaling factor to be applied to data;
*      (I)  BMS_INPUT *bms;
*           points to the decoded internal Bit Map Section Struct 
*      (I/O) grid_desc_sec *gds
*           points to decoded internal grid definition struct
*      (O)  float **ppgrib_data;
*           double pointer to array of float, null upon entry;
*           upon successful exit, holds the unpacked and restored float data
*           in a newly malloced array;
*      (O)  BDS_HEAD_INPUT *bds_head; 
*           points to Binary Data Sect hdr struct, empty upon entry;
*           to be filled with decoded BDS info;
*      (O)  char  *errmsg;
*           Only returned filled when error occurred;
*
*      RETURN CODE:
*        0>   no errors
*        1>   unrecognized packing algorithm
*        2>   number of points does not match bitmap
*        3>   number of points does not match grid size in GDS
*        4>   malloc error
* ********************************************************************
*/

#if PROTOTYPE_NEEDED
int gribgetbds ( char *curr_ptr, short  deci_scale, BMS_INPUT *bms,
        grid_desc_sec *gds, float **ppgrib_data,
        BDS_HEAD_INPUT *bds_head, char  *errmsg)
#else
int gribgetbds (curr_ptr, deci_scale, bms, gds, ppgrib_data, 
                bds_head, errmsg)
                char *curr_ptr; 
                short  deci_scale; 
                BMS_INPUT *bms;
                grid_desc_sec *gds;
                float **ppgrib_data; 
                BDS_HEAD_INPUT *bds_head; 
                char  *errmsg;
#endif
{
char *func="gribgetbds";
char *in = curr_ptr;      /* pointer to beginning of BDS */
long length;              /* size of the Binary Data Section */
long scale;               /* scaling factor */
float ref_val;            /* reference value (minimum value) */
unsigned long something;  /* generic value from message */
long data_width;          /* number of bits that data occupies */
int halfBYTE4;            /* the first 4 bits in 4-th byte */
int status=0;
int sign;                 /* sign + or - */
float dscale;             /* 10 to the decimal scaling power */
float bscale;             /* 2 to the binary scaling power   */
unsigned long skip=0;              /* number of bits to be skipped */
long c_ristic;            /* characteristic for float representation */
long mantissa;            /* mantissa for float representation */
long numpts;              /* number of bits left at end of bitstream */
unsigned long data_pts;   /* number of data points in bitstream */
unsigned long num_calc;   /* temp work var */
float *grib_data=0;       /* local work array for grid data */
float fdata=0;            /* data value stored in reference */
int     i,j;         /* array counter */
int xsize, ysize; 
float *outdata;

  DPRINT1 ("Entering %s()\n", func);
/*
*
* A.1       FUNCTION gbyte !get bds length
*/

  gbyte(in,(unsigned long *) &length,&skip,24);
  DPRINT0 ("bds_head->length\n");
  bds_head->length = (unsigned long) length;

/*
*
* A.2       FUNCTION gbyte !get BDS flag 
*/
  gbyte(in, (unsigned long *) &halfBYTE4, &skip, 4);
  DPRINT0 ("bds_head->usBDS_flag\n");
  bds_head->usBDS_flag = (short) halfBYTE4;  /* get BDS Flag (Table 11) */

/*
*
* A.3       IF (unsupported packing algorithm)  THEN
*               RETURN 1
*           ENDIF
*/
  /* need to check on packing algorithm */
  if (halfBYTE4)  /* unrecognized packing algorithm */
    {
     DPRINT1 ("%s:  error, unrecognized packing algorithm\n", func);
     sprintf(errmsg, "%s:  unrecognized packing algorithm\n", func);
     status= (1);
     goto BYE;
    }

/*
*
* A.4       FUNCTION gbyte !get number of unused bits
*/
  gbyte(in,(unsigned long *) &numpts,&skip,4); /* get #bits at end of BDS */
  DPRINT0 ("numpts\n");

/*
*
* A.5       FUNCTION gbyte !get Binary Scale Factor
*/
  gbyte(in,&something,&skip,16);
  DPRINT0 ("Sign & bds_head->Bin_sc_fctr\n");
  sign = (int)(something >> 15) & 1;  /* get sign for scale */
  scale = (int)(something) & 32767;   /* get scale */
  if(sign)                            /* scale negative */
     scale = -scale;                  /* multiply scale by -1 */
  bds_head->Bin_sc_fctr = (int) scale;  /* get binary scale factor */
  DPRINT1 ("Binary Scale Factor = %d\n", scale);

/*
*
* A.6       CALCULATE Reference value from IBM representation
*             !FUNCTION gbyte !get the sign of reference
*             !FUNCTION gbyte !get charateristic
*             !FUNCTION gbyte !get the mantissa
*/
  gbyte(in,&something,&skip,8);
  DPRINT0 ("Sign & Reference)\n");
  sign = (int)(something >> 7) & 1; /* get the sign for reference value */

  skip -= 7;
  gbyte(in,(unsigned long *)&c_ristic,&skip,7); /*characteristic for the float*/
  DPRINT0 ("c_ristic\n");

  gbyte(in,(unsigned long*)&mantissa,&skip,24); /*mantissa for the float */
  DPRINT0 ("mantissa\n");
  c_ristic -= 64;                   /* subtract 64 from characteristic */
  ref_val = (float) mantissa * (float)(pow(16.0,(double)c_ristic)) * 
            (pow(2.0,-24.0));
  if(sign)                          /* negative reference value */
     ref_val = -ref_val;            /* multiply ref_val by -1 */
  bds_head->fReference = (float)ref_val;
  DPRINT1 ("Reference = %f\n", ref_val);

/*
*
* A.7       FUNCTION gbyte !get data width
*/
  gbyte(in, (unsigned long*) &data_width,&skip,8);       /* get data width */
  DPRINT0 ("bds_head->usBit_pack_num\n");
  bds_head->usBit_pack_num = (short)data_width;

/*
*
* A.8       SET Binary and Decimal Scale Factors
*/
  /* set binary scale */
  bscale = (float)pow (2.0,(double) scale);

  /* set decimal scale */
  dscale = (float)pow (10.0, (double) deci_scale);

  DPRINT2 ("Scaled-up BSF= (2**%d) = %f\n", scale, bscale);
  DPRINT2 ("Scaled-up DSF= (10**%d) = %f\n", deci_scale,dscale);

/*
*
* A.9       IF (data_width is zero) THEN
*               ! grid contains a constant value
*               IF expected grid size is invalid THEN
*                  FORCE grid size of 1
*               ENDIF
*               ALLOCATE array for grid of expected grid size
*               FILL grid with the constant value
*               RETURN 0 !success
*           ENDIF
*/
  if (!data_width)   /* grid contains constant value, success, all done */
  {
/* Used to send back array of 1 element:
#     bds_head->ulGrid_size = 1;
#     fdata = (float) (ref_val / dscale);
#     *ppgrib_data = (float *) malloc(sizeof(float));
#     **ppgrib_data = fdata;
*/

     fdata = (float) (ref_val / dscale);
     if (bds_head->ulGrid_size <= 0) {
        fprintf(stdout,
        "WARNING:  gribgetbds detects bad ulGrid_size (%ld); "\
        "Set to 1 to hold constant value %lf\n", bds_head->ulGrid_size, fdata);
        bds_head->ulGrid_size = 1;
        }
    
     grib_data =(float *) malloc(bds_head->ulGrid_size * sizeof(float));
     if (!grib_data) {
        sprintf(errmsg,
        "%s: failed to create array[%d] for grid to hold Constant data", 
        func, bds_head->ulGrid_size);
        goto BYE;
        }

     *ppgrib_data = grib_data;          /* store address */
     for (i=0; i < bds_head->ulGrid_size ; ) 
        grib_data[i++] = fdata;         /* fill grid with constant val */

     DPRINT3("%s:  grid[%ld] contains convant value %lf\n",
        func, bds_head->ulGrid_size, fdata);

     status = (0);  /* no errors */
     goto BYE;
  }

  /* fill the data array with values from message */
  /*     - Assume that GDS may not be included so that
   *         the number of grid points may not be defined.       
   *     - Compute space to malloc based on BDS length,
   *         data_width, and numpts.
   *     - if grid_size from GDS is zero, use
   *         computed number of points.
   */

/*
*
* A.10      CALCULATE number of data points actually in BDS
*/
  num_calc = ((length - 11)*8 - numpts) / data_width;

  /* Check the number of points computed against info in the BMS
     or GDS, if they are available */

/*
*
* A.11      IF (BMS is present and has included bitmap) THEN
*               IF (#calculated not same as #bits set in BMS) THEN
*                   RETURN 2
*               ENDIF
*/

  if (bms->uslength > 6)
  {
      if (bms->ulbits_set != num_calc) {
        DPRINT3 ("%s:  BMS present, #datapts calculated (%d) " \
        "not same as BMS's set bits (%d)\n",func, num_calc, bms->ulbits_set);

        sprintf(errmsg,"%s:  BMS present, #datapts calculated (%d) " \
        "not same as BMS's set bits (%d)\n",func, num_calc, bms->ulbits_set);
        status= (2); goto BYE;
        }
  }
/*
* A.11.1    ELSE  !no bms
*               IF (GDS is present AND
*                   #calculated not same as GDS's grid size)
*               THEN
*                   RETURN 3
*               ENDIF
*/
  
  else
    
  {    
      if (bds_head->ulGrid_size && bds_head->ulGrid_size != num_calc) {
        DPRINT0("Averting failure of grid size test\n");
        /*
        DPRINT3 ( "%s:  GDS present, #datapts calculated (%d) " \
           "not same as GDS's grid size (%d)\n",
           func,num_calc,bds_head->ulGrid_size);

           sprintf(errmsg,"%s:  GDS present, #datapts calculated (%d) " \
           "not same as GDS's grid size (%d)\n",
           func,num_calc,bds_head->ulGrid_size);
           status=(3); goto BYE;
        */
         }          
  }
    
  
/*
* A.11      ENDIF (BMS present)
*/

  /* Only reach this point if number of points in BDS matches info
     in BMS or GDS.  This number is unchecked if no BMS or GDS. */

  /* Make sure number of points in BDS is value used for extracting */
/*
*
* A.12      SET #datapoints
*/
  bds_head->ulGrid_size = num_calc;
  data_pts= num_calc;

/*
*
* A.13      ALLOCATE storage for float array size
*           IF (error) THEN
*               RETURN 4
*           ENDIF
*/
  grib_data =(float *) malloc(data_pts * sizeof(float));
  if (grib_data==NULL) { 
        DPRINT1 ("%s: failed to malloc Grib_Data\n",func);
        sprintf(errmsg,"%s: failed to malloc Grib_Data\n",func);
        status=(4); goto BYE; }

/*
*
* A.14      SET data array pointer to local data array
*/
  *ppgrib_data = grib_data;

/*
*
* A.15      FOR (each data point) DO
*               FUNCTION gbyte_quiet   !get data_width bits
*               INCREMENT skip by data_width
*               COMPUTE and STORE value in float array
*           ENDDO
*/

  DPRINT3 ( "Restore float data by = (float)(%f + X * %f)) / %f;\n",
  ref_val, bscale, dscale);

  for(i=0;i < data_pts ;i++)  
     {
     gbyte_quiet(in,&something,&skip,data_width);
     grib_data[i]= (float)(ref_val + (something * bscale))/dscale;
     }

/*
 * Unthin grid if it is thinned
 */
  if (gds->head.thin != NULL) {
    if ((gds->head.usData_type == LATLON_PRJ) || 
        (gds->head.usData_type == GAUSS_PRJ) ||
        (gds->head.usData_type == ROT_LATLON_PRJ) ||
        (gds->head.usData_type == ROT_GAUSS_PRJ) ||
        (gds->head.usData_type == STR_LATLON_PRJ) ||
        (gds->head.usData_type == STR_GAUSS_PRJ) ||
        (gds->head.usData_type == STR_ROT_LATLON_PRJ) ||
        (gds->head.usData_type == STR_ROT_GAUSS_PRJ)) {
      ysize = gds->llg.usNj;
    } else if (gds->head.usData_type == MERC_PRJ) {
      ysize = gds->merc.rows;
    } else if (gds->head.usData_type == POLAR_PRJ) {
      ysize = gds->pol.usNy;
    } else if ((gds->head.usData_type == LAMB_PRJ) ||
               (gds->head.usData_type == ALBERS_PRJ) ||
               (gds->head.usData_type == OBLIQ_LAMB_PRJ)) {
      ysize = gds->lam.iNy;
    } else {
       DPRINT2 ("%s:  unknown datatype=%d\n",func, gds->head.usData_type);
       sprintf(errmsg,"%s:  unknown datatype=%d\n",func, gds->head.usData_type);
       status=1;         /* set status to failure */
    }

    xsize = 0;
    for (j = 0; j<ysize; j++) {
      if (gds->head.thin[j] > xsize) {
        xsize = gds->head.thin[j];
      }
    }
    outdata = (float *)malloc(ysize*xsize*sizeof(float));
    grib_unthin(grib_data,outdata,gds->head.thin,ysize,
          &xsize);
    free(grib_data);
    grib_data = (float *)malloc(sizeof(float)*ysize*xsize);
    *ppgrib_data = grib_data;
    memcpy(grib_data,outdata,sizeof(float)*ysize*xsize);
    free(outdata);
    free(gds->head.thin);
    gds->head.thin = NULL;
    if ((gds->head.usData_type == LATLON_PRJ) || 
        (gds->head.usData_type == GAUSS_PRJ) ||
        (gds->head.usData_type == ROT_LATLON_PRJ) ||
        (gds->head.usData_type == ROT_GAUSS_PRJ) ||
        (gds->head.usData_type == STR_LATLON_PRJ) ||
        (gds->head.usData_type == STR_GAUSS_PRJ) ||
        (gds->head.usData_type == STR_ROT_LATLON_PRJ) ||
        (gds->head.usData_type == STR_ROT_GAUSS_PRJ)) {
      gds->llg.usNi = xsize;
      gds->llg.iDi = abs(gds->llg.lLat2 - gds->llg.lLat1)/(xsize-1);
    } else if (gds->head.usData_type == MERC_PRJ) {
      gds->merc.cols = xsize;
    } else if (gds->head.usData_type == POLAR_PRJ) {
      gds->pol.usNx = xsize;
    } else if ((gds->head.usData_type == LAMB_PRJ) ||
               (gds->head.usData_type == ALBERS_PRJ) ||
               (gds->head.usData_type == OBLIQ_LAMB_PRJ)) {
      gds->lam.iNx = xsize;
    } else {
       DPRINT2 ("%s:  unknown datatype=%d\n",func, gds->head.usData_type);
       sprintf(errmsg,"%s:  unknown datatype=%d\n",func, gds->head.usData_type);
       status=1;         /* set status to failure */
    }
  }

  DPRINT0 ("Sample of first 30 unpacked & restored datapoints=\n");
  for (i=0; i < 30; i+=5)
     DPRINT6 ("%03d:  %f %f %f %f %f\n",
        i, grib_data[i], grib_data[i+1],grib_data[i+2],
        grib_data[i+3], grib_data[i+4] );
  
BYE:
/*
*
* A.16      RETURN Status;
*/
  DPRINT2 ("Exiting %s, status=%d\n", func, status);
  return(status); 
/* 
* END OF FUNCTION
*
*/
}