source: trunk/WRF.COMMON/WRFV2/external/io_grib1/MEL_grib1/pack_spatial.c

/* 20jun97/atn:  always scaling data up whether dsf is -/+; 
*/
#include <stdio.h> 
#include <stdlib.h>
#include <string.h> 
#include <math.h>
#ifdef XT3_Catamount
#include <features.h>
#undef htonl
#define htonl(x)     swap_byte4(x)
#else
#include <netinet/in.h>
#endif
#include "dprints.h"            /* for dprints */
#include "gribfuncs.h"          /* prototypes */
#include "isdb.h"               /* WORD_BIT_CNT defn */

/*
****************************************************************
* A. FUNCTION:   pack_spatial
*      pack gridded data values into a bitstream 
*
*    INTERFACE:
*      int     pack_spatial (pt_cnt, bit_cnt, pack_null, fbuff, ppbitstream,
*                           dec_scl_fctr, BDSlength, errmsg)
*
*    ARGUMENTS (I=input, O=output, I&O=input and output):
*     (I) long *pt_cnt;                count of points in grid 
*     (I) long *bit_cnt;               count of bits to pack value in.  
*                                      will be calculated if set to zero
*     (I) float *pack_null;            parameter value for null (huge number) 
*   (I&O) float *fbuff;                array containing grid values to pack
*                                     returned scaled up by Decimal Scale Factor
*     (O) unsigned long **ppbitstream; Null upon entry; 
*                                      returned pointing to new Storage 
*                                      holding packed bitstream;
*     (I) short dec_scl_fctr;          decimal scale factor to apply to data
*     (O) long  *BDSlength;            updated with #bytes in packed bitstream
*     (O) char  *errmsg                returned filled if error occurred
*
*    RETURN CODE: 
*      0> success, ppbitstream contains packed values
*   else> error: errmsg holds msg;
****************************************************************
*/
#if PROTOTYPE_NEEDED
int     pack_spatial ( long             *pt_cnt,
                        unsigned short  *bit_cnt,
                        float           *pack_null,
                        float           *fbuff,
                        unsigned long   **ppbitstream,
                        short           dec_scl_fctr,
                        long            *BDSlength,
                        char            *errmsg)

#else
int     pack_spatial (  pt_cnt, bit_cnt, pack_null, fbuff, ppbitstream,
                        dec_scl_fctr, BDSlength, errmsg)
                        long            *pt_cnt;
                        unsigned short  *bit_cnt;
                        float           *pack_null;
                        float           *fbuff;
                        unsigned long   **ppbitstream;
                        short           dec_scl_fctr;
                        long            *BDSlength;
                        char            *errmsg;
#endif
{
    char *func="pack_spatial";
    long ipt;                   /* index over points */
    int null_flag;              /* flag indicating presence of null values */
    int bit1;                   /* starting bit in current word */ 
    int empty;                  /* number of empty bits in word */ 
    int diff;                   /* difference of empty - bit1 */ 
    long max_value;             /* max value storable in bit_cnt bits */
    unsigned long itemp;        /* temporary unsigned integer */
    unsigned long *bstr;        /* pointer running across bitstream */
    int pack_bit_cnt;           /* count of bits to pack parameter values */ 
    int unused_bit_cnt;         /* count of unused bits for i2 words */
    /*long byte4_cnt;           /- count of bytes using i4 words */
    long byte2_cnt;             /* count of bytes using i2 words */
    short scl_fctr;             /* scaling factor for grid values */
    double pow_scl;             /* 2 ** (-scl_fctr) */
    double pwr10toD;            /* 10 ** (D) */
    float reference;            /* reference = minimum value in grid */
    float max_grid;             /* maximum value in grid */
    float ftemp;                /* temporary float containing grid value */
    unsigned long *pBitstream;
    unsigned long grib_local_ibm();
    int wordnum;
    int zero_cnt;
    int prec_too_high = 0;
    unsigned char bdshdr[14];   /* Character array to temporarily hold bds
                                 *   header */
    int hdrwords;

    DPRINT1 ( "Entering %s....\n", func );

/*
*
* A.1       IF (no data in grid) THEN
*              PRINT message
*              RETURN Stat= -1
*           ENDIF
*/
    if (*pt_cnt <= 0) {
        DPRINT2 ("%s; invalid pt_cnt = %d\n", func,*pt_cnt);
        sprintf(errmsg, "%s; invalid pt_cnt = %d\n", func,*pt_cnt);
        return (-1);
    }

/*
*
* A.2       IF (number of bits to pack into is greater than 30) THEN
*               PRINT message
*               RETURN Stat= -1
*           ENDIF
*           SET pack_bit_cnt for local use
*/
    if ( *bit_cnt > 30 ) {
        DPRINT2 ("%s; invalid bit_cnt = %d\n", func,*bit_cnt);
        sprintf(errmsg, "%s; invalid bit_cnt = %d\n", func,*bit_cnt);
        return (-1);
    }
    pack_bit_cnt        = (int) *bit_cnt; 
    DPRINT1 ("  use Pack_bit_cnt= %d\n", pack_bit_cnt);

/*
*               
* A.3      FOR (each data point) DO
*              SCALE all values of data array  !multiply by 10**DSF
*          ENDDO
*/
        pwr10toD=  pow ( 10., (double) dec_scl_fctr );
        for (ipt=0; ipt < *pt_cnt; ipt++)  fbuff[ipt] *= pwr10toD;

        DPRINT2 ("  Decimal Scale Fctr= %d, scale data by 10**dsf "\
        "(Fbuff *= %lf)\n", dec_scl_fctr,  pwr10toD);
/*
*
* A.4       INIT reference, max_grid, null_flag 
*/
    reference   = 1.e30;
    max_grid    = -1.e30;
    null_flag   = 0;

/*
*
* A.5       FOR (each data point) DO
*              IF (value < reference) THEN
*                  SET reference to this value   !smallest value
*              ENDIF
*              IF (value > max_grid AND not a missing value ) THEN
*                  SET max_grid to this value    !largest value
*              ENDIF
*              IF (value >= missing value ) THEN
*                  SET null_flag to 1            !grid contains nulls
*              ENDIF
*           ENDDO
  Find reference (minimum) and maximum values of the grid points 
*/
    for (ipt = 0; ipt < *pt_cnt; ipt++) {
        ftemp   = *(fbuff+ipt);
        if (ftemp < reference) reference = ftemp;    /* REF is SCALED UP */
        if (ftemp > max_grid && ftemp < *pack_null) max_grid = ftemp;
        if (ftemp >= *pack_null) null_flag = 1;
    }

    DPRINT2 ("  Max before taking out Ref =%.4lf\n  Null flag=%d\n",
    max_grid, null_flag);

/*  Compute maximum range of grid (max_grid - reference) */
/*
*
* A.6       IF (max value is same as smallest value AND
*               null_flag is zero) THEN
*               CLEAR pack_bit_cnt  !constant values, no nulls
*               CLEAR max_grid      !set grid range to 0
*/
    if (((max_grid - reference) < 1.0) && null_flag == 0) {   
        pack_bit_cnt = 0;
        max_grid =  0;

/*
* A.6.a     ELSE IF (max value is same as smallest value AND
*                    null_flag is set) THEN
*               SET max_grid to 1      !const values, some nulls
*/
    } else if (((max_grid - reference) < 1.0) && null_flag == 1) {
        max_grid = 1.;

/*
* A.6.b     ELSE IF (max value <= -1.e29 AND null_flag is set) THEN
*               PRINT message
*               RETURN Stat= -1
*/
    } else if (max_grid <= -1.e29 && null_flag == 1) {
        DPRINT1 ("%s; Grid contains all NULLS\n",func);
        sprintf(errmsg, "%s; Grid contains all NULLS\n",func);
        return (-1);

/*
* A.6.c     ELSE IF (max value not equal to reference) THEN
*               SET max_grid (max_grid-reference) !non-constant values w/wo nulls
*/
    } else if (max_grid != reference) {
        max_grid -= reference;

/*
* A.6       ENDIF
*/
    }

/*
* 
* A.7       DEBUG print grid range and reference value
*/
    DPRINT2 ( "  Reference = %f\n  Max_grid (range) = %f\n",  
    reference, max_grid);

/*  Find minimum number of bits to pack data */
/*
*
* A.8.a     IF (grid range is not zero) THEN
*/
    if ( max_grid != 0 )    
    {

/*
*
* A.8.a.1      DEBUG print input bit width
*              IF (input bit_num is zero) THEN
*                 CALCULATE number of bits needed to store grid range
*                 DEBUG print calculated bit count
*              ENDIF
*/
       DPRINT1 ( "  Input bit cnt = %d\n", pack_bit_cnt );
       if ( pack_bit_cnt == 0 )
       {
          pack_bit_cnt = (int)(floor(log10((double)max_grid) / log10(2.)) +1); 
          DPRINT1 ( "  Calculated bit cnt = %d\n", pack_bit_cnt );
       }
       if ( (pack_bit_cnt < 0) || (pack_bit_cnt > 30) )
       {
          DPRINT1 ("%s: Calculated bit count OUT OF RANGE [0 - 30] !!\n", func);
          sprintf (errmsg, "%s: Calculated bit count OUT OF RANGE!! bit_cnt: %d  max: %f\n", func,pack_bit_cnt,max_grid);
          return (-1);
       }
/*
*
* A.8.a.2      CALCULATE various byte counters
*              !itemp: #bits required for header + grid
*              !Byte2_cnt: #bytes rounded up to next 2-byte bdry
*              !Byte4_cnt: #bytes rounded up to next 4-byte bdry
*              !Unused_bit_cnt: #unused bits at end using byte2_cnt
*              DEBUG print expected length and unused bits
*/
       itemp = *pt_cnt * pack_bit_cnt + 11 * BYTE_BIT_CNT;
       byte2_cnt = (long) ceil(((double) itemp / BYTE_BIT_CNT) / 2.) * 2;
       /*byte4_cnt = (long) ceil(((double) itemp / BYTE_BIT_CNT) / 4.) * 4;*/
       unused_bit_cnt = byte2_cnt * BYTE_BIT_CNT - itemp;
       DPRINT1 ( "  Calculated length = %ld bytes (Rnd2)\n", byte2_cnt);
       DPRINT1 ( "  Bitstream padding = %ld bits\n",unused_bit_cnt);

/*
*
* A.8.a.3      CALCULATE maximum storable value 
*              CALCULATE scl_fctr required to fit grid range
*                        into available bit width
*/
       max_value = (long) pow(2., (double) pack_bit_cnt) - 1;
       if (max_value < 2) max_value = 2;
       scl_fctr = -(short) floor(log10((double) (max_value-1) / 
                          (double) max_grid) / log10(2.));
       pow_scl  = pow(2., (double) -scl_fctr); 
       DPRINT1 ( "  Calculated Binary scale = %d\n",scl_fctr);
    }

/*
*
* A.8.b     ELSE       !max_grid = 0, all zero data or constant values
*              SET number of bits to pack to zero
*              SET lengths to 12 bytes
*              SET unused bits to 8 (1 byte of padding)
*              SET scl_fctr to 0
*              DEBUG print constant grid
*           ENDIF
*/
    else
    {
       pack_bit_cnt = 0;
       byte2_cnt = 12;
       /*byte4_cnt = 12;*/
       unused_bit_cnt   = 8;
       scl_fctr = 0;
       DPRINT0 ( "  Constant grid.  Using bit cnt = 0\n");
    }

/*
*
* A.9       MALLOC space for bitstream (Rnd2_cnt)
*           IF (failed) THEN
*              PRINT error mesg
*              RETURN Stat= 999;
*           ENDIF
*/ 
    pBitstream = ( unsigned long * ) malloc ( sizeof( unsigned long ) * 
                                              byte2_cnt );
    if ( !pBitstream )
    {
       DPRINT1 ("%s:  MAlloc failed pBitstream\n", func );
       sprintf(errmsg, "%s:  MAlloc failed pBitstream\n", func );
       return (999);
    }

/*
*
* A.10      SET ptr to bitstream
*           UPDATE bit_cnt for input structure
*/
    *bit_cnt = (unsigned short) pack_bit_cnt;
    DPRINT1 ("  Updated input bit cnt to %d\n", *bit_cnt);
    bstr = pBitstream;
 
/*
*
* A.11      ZERO out entire bitstream
*/
   zero_cnt = ceil(byte2_cnt / (float)sizeof(long)) * sizeof(long);
   memset ((void *)pBitstream, '\0', zero_cnt);

/*
*
* A.12      PUT packing info into first 11 bytes: 
*           NOTE: The Table 11 Flag stored in the first
*                 4 bits of Octet 4 is HARDCODED to 0000.
*                 This implies Simple packing of float
*                 grid point data with no additional flags.
*           Octet 1-3  = Byte2_cnt
*           Octet 4    = Table 11 Flag & unused_bit_cnt
*           Octet 5-6  = Scl_fctr
*           Octet 7-10 = Reference truncated to DSF precision
*           Octet 11   = Pack_bit_cnt
*           Octet 12   = Bitstream starts (bit 25 of word 3)
*/


   set_bytes_u(byte2_cnt, 3, bdshdr);

   itemp = unused_bit_cnt;
   set_bytes_u(itemp, 1, bdshdr+3);

   set_bytes_u(scl_fctr, 2, bdshdr+4);

   DPRINT1 ("  Reference (%f) ", reference);
   reference   =  floor((double) reference + .5);
   DPRINT1 ("truncated to DSF precision= %lf\n", reference);
   itemp        = grib_local_ibm(reference);
   DPRINT1 ("  Reference converted to local IBM format= 0x%x\n", itemp);

   set_bytes_u(itemp, 4, bdshdr+6);

   set_bytes_u(pack_bit_cnt, 1, bdshdr+10);

   bit1 = 25;

   memcpy(bstr,bdshdr,11);

   /* 
    * For non-internet byte order machines (i.e., linux), 
    * We reverse the order of the last byte in the bds header, since
    *   it will be reversed once again below.
    */
   hdrwords = 11/(WORD_BIT_CNT/BYTE_BIT_CNT);
   set_bytes_u(bstr[hdrwords], WORD_BIT_CNT/BYTE_BIT_CNT,
              (char *)(bstr+hdrwords) );

   bstr += hdrwords;


/*
    itemp       = unused_bit_cnt;
    *bstr       = (byte2_cnt << 8) | itemp;
    bstr++;
    *bstr       = scl_fctr;

    DPRINT1 ("  Reference (%f) ", reference);
    reference   =  floor((double) reference + .5);
    DPRINT1 ("truncated to DSF precision= %lf\n", reference);
    itemp       = grib_local_ibm(reference);
    DPRINT1 ("  Reference converted to local IBM format= 0x%x\n", itemp);

    *bstr       = (*bstr << 16) | (itemp >> 16);
    bstr++;
    *bstr       = (itemp << 16) | (pack_bit_cnt << 8);
    bit1        = 25;   */       /* starting bit within current bstr word */

/*
*
* A.13      IF (grid values are not constant) THEN
*/ 
    if (pack_bit_cnt > 0) {


/*
* A.13.1       SET empty value
*/
        empty = WORD_BIT_CNT - pack_bit_cnt + 1;

        for (ipt=0; ipt < 5; ipt++) DPRINT4 (
            " ITEMP= (*(fbuff+ipt) - reference) * pow_scl + .5=\n"\
            "        (%lf -%lf) * %lf + .5 = %lf\n",
            *(fbuff+ipt), reference, pow_scl, 
            (*(fbuff+ipt) - reference) * pow_scl + .5);
           
/*
* A.13.2       FOR (each point in bitstream) DO
*/
        for (ipt = 0; ipt < *pt_cnt; ipt++) {

/*
* A.13.2.1         IF ( data value < pack_null) THEN
*                     SET itemp to (value - reference) * pow_scl + .5;
*                  ELSE
*                     SET itemp to max value;
*                  ENDIF
*/
            if (*(fbuff+ipt) < *pack_null) {
                itemp  = (*(fbuff+ipt) - reference) * pow_scl + .5;
            } else {
                itemp  = max_value; 
                DPRINT1 ("%s: Setting to max_value: Precision may be too high !!\n", func);
                sprintf (errmsg, "%s: Setting grid point to max value, precision may be too high", func);
                /*              return (-1); */
            }

/*
* A.13.2.2         COMPUTE if data point can fit in current word
*/
            diff        = empty - bit1;

/*
* A.13.2.3.a       IF (data point falls within word ) THEN
*                      SHIFT value to the correct bit position 
*                      COMBINE it with data in current word of bitstream
*                      CALCULATE starting bit in curr word for next time
*/
            if (diff > 0)
            {
                *bstr   |= itemp << diff;
                bit1    += pack_bit_cnt;
            }

/*
* A.13.2.3.b       ELSE IF (data point ends at word boundary) THEN
*                      COMBINE value with data in current word of bitstream
*                      SET starting bit to 1
*                      BUMP word counter in bitstream up by 1 word
*/
            else if (diff == 0) 
            {
                *bstr |= itemp;
                bit1   = 1;
                bstr++;
            }

/*
* A.13.2.3.c       ELSE     !point crosses word boundary
*                      STORE "diff" bits of value in current word of bitstream
*                      BUMP word counter in bitstream up by 1 word
*                      STORE remaining bits of value in next word
*                      CALCULATE starting bit in curr word for next time
*                  ENDIF    !word location check
*/
            else                /* pixel crosses word boundary */
            {
                *bstr |= itemp >> -diff;
                bstr++;
                *bstr |= itemp << (WORD_BIT_CNT + diff);
                bit1   = -diff + 1;
            }

/*
* A.13.2        ENDFOR loop over grid points
*/
       }

/*
* A.13       ENDIF (pack_bit_cnt > 0)
*/
    }

/* For little endian machines, swap the bytes in the bstr pointer */
    /*    for (wordnum = 0; */
    for (wordnum = hdrwords;
         wordnum < ceil(byte2_cnt/(float)(WORD_BIT_CNT/BYTE_BIT_CNT)); 
         wordnum++) {
      set_bytes_u(pBitstream[wordnum], WORD_BIT_CNT/BYTE_BIT_CNT, 
                  (char *)(pBitstream+wordnum) );
    }

/*
*
* A.14       ASSIGN bitstream block to ppbitstream pointer
*            SET BDSlength (size rnded to next 2 byte boundary)
*            RETURN Status 0  ! success
*/
    *ppbitstream = pBitstream;
    *BDSlength = (long) byte2_cnt;

    DPRINT1 ("Exiting pack_spatial, BDS len=%ld,  Status=0\n" , *BDSlength);

    if (prec_too_high) {
        fprintf(stderr,"pack_spatial: Warning: Precision for a parameter may be too high in gribmap.txt\n");
    }
    return (0);
/*
* END OF FUNCTION
*
*/
}

/*
*
**************************************************************
* B. FUNCTION  grib_local_ibm
*      convert local_float from local floating point to
*      IBM floating point stored in a 4-byte integer.
*
*    INTERFACE:
*      unsigned long grib_local_ibm (local_float)
*
*    ARGUMENTS (I=input, O=output, I&O=input and output):
*      (I)  double local_float      float value in local format
*
*    RETURNS:
*       the actual IBM floating point value
**************************************************************
*
*/
#if PROTOTYPE_NEEDED
unsigned long grib_local_ibm (double local_float)
#else
unsigned long grib_local_ibm (local_float)
double local_float;
#endif
{
    long a, b; 
    unsigned long ibm_float;
/*
*
* B.1.a     IF (local float value is zero) THEN
*               SET the ibm float to zero too
*/
    if (local_float == 0.) { 
        ibm_float = 0;
    } else {
/*
* B.1.b     ELSE
*              CONVERT to IBM floating point
*              ! IBM floating point is stored in 4 bytes as:  
*              ! saaaaaaa bbbbbbbb bbbbbbbb bbbbbbbb
*              ! where s is sign bit, 0 -> positive, 1 -> negative
*              !       a is 7-bit characteristic
*              !       b is 24-bit fraction
*              ! s, a and b are obtained from local_float (local 32-bit float) as
*              !       s = sign(local_float)
*              !       a = ceil(log10(local_float) / log10(16.)) + 64 
*              !       b = local_float / 16**(a-64) * 2**24
* B.1.b     ENDIF
*/
        a = ceil(log10(fabs(local_float)) / log10(16.)) + 64;
        /* Added by Todd Hutchinson, 8/13/99 */
        /* This fixes a problem when local_float == 256, etc. */
        if ( fmod((log10(fabs(local_float))/log10(16.)),1.) == 0) {
          a++;
        }
/*      Local_float == +/-1. is a special case because of log function */
        if ( (local_float == 1.) || (local_float == -1.)) a = 65;
        b = (long) (fabs(local_float) * pow(16.,(double) (70 - a)) +.5)
            & 0x00ffffff; 
        ibm_float = (((local_float > 0.) ? 0 : 1) << 31) | (a << 24) | b;
    }
/*
*
* B.2       RETURN the ibm float value
*/
    return ibm_float;
/*
*
* END OF FUNCTION
*
*/ }

/*
*
**************************************************************
* C. FUNCTION:   grib_ibm_local
*      convert local_float from IBM floating point to
*      local floating point.
*
*    INTERFACE:
*      float grib_ibm_local(ibm_float)
*
*    ARGUMENTS (I=input, O=output, I&O=input and output):
*      (I)  double local_float      float value in local format
*
*     RETURNS:
*        the actual local floating point
**************************************************************
*/
#if PROTOTYPE_NEEDED
float grib_ibm_local( unsigned long ibm_float)
#else
float grib_ibm_local( ibm_float)
unsigned long ibm_float;
#endif
{
/*
  Convert ibm_float from IBM floating point stored in a 4-byte integer 
  to local floating point.
*
* C.1       DETERMINE local floating point
*           ! IBM floating point is stored in 4 bytes as:  
*           ! saaaaaaa bbbbbbbb bbbbbbbb bbbbbbbb
*           ! where s is sign bit, 0 -> positive, 1 -> negative
*           !       a is 7-bit characteristic
*           !       b is 24-bit fraction
*           ! local_float (local 32-bit float) is recovered from 
*           ! s, a and b as
*           ! local_float = (-1)**s * 2**(-24) * b * 16**(a-64)
*/

    long a, b;
    float local_float;

    a = (ibm_float >> 24) & 0x0000007f;
    b = ibm_float & 0x00ffffff;
    local_float = (float) b * pow(16., (double) (a - 70));
    if (ibm_float >> 31) local_float = -local_float;
/*
*
* C.2       RETURN floating point
*/
    return local_float;
/*
* END OF FUNCTION
*
*/ 
}