/* 20jun97/atn: always scaling data up whether dsf is -/+; */ #include #include #include #include #ifdef XT3_Catamount #include #undef htonl #define htonl(x) swap_byte4(x) #elif defined(_WIN32) #include #else #include #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 * */ }