/wrfv2_fire/external/io_grib1/MEL_grib1/pack_spatial.c
C | 652 lines | 257 code | 56 blank | 339 comment | 52 complexity | 25355c1a41c79f4e2f25f1c4a3f7ad45 MD5 | raw file
Possible License(s): AGPL-1.0
- /* 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)
- #elif defined(_WIN32)
- #include <Winsock2.h>
- #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
- *
- */
- }