Line data Source code
1 : #include <stdlib.h>
2 : #include <math.h>
3 : #include <limits.h>
4 : #include <float.h>
5 : #include "grib2.h"
6 :
7 :
8 65 : void simpack(g2float *fld,g2int ndpts,g2int *idrstmpl,unsigned char *cpack,g2int *lcpack)
9 : //$$$ SUBPROGRAM DOCUMENTATION BLOCK
10 : // . . . .
11 : // SUBPROGRAM: simpack
12 : // PRGMMR: Gilbert ORG: W/NP11 DATE: 2002-11-06
13 : //
14 : // ABSTRACT: This subroutine packs up a data field using the simple
15 : // packing algorithm as defined in the GRIB2 documentation. It
16 : // also fills in GRIB2 Data Representation Template 5.0 with the
17 : // appropriate values.
18 : //
19 : // PROGRAM HISTORY LOG:
20 : // 2002-11-06 Gilbert
21 : //
22 : // USAGE: CALL simpack(fld,ndpts,idrstmpl,cpack,lcpack)
23 : // INPUT ARGUMENT LIST:
24 : // fld[] - Contains the data values to pack
25 : // ndpts - The number of data values in array fld[]
26 : // idrstmpl - Contains the array of values for Data Representation
27 : // Template 5.0
28 : // [0] = Reference value - ignored on input
29 : // [1] = Binary Scale Factor
30 : // [2] = Decimal Scale Factor
31 : // [3] = Number of bits used to pack data, if value is
32 : // > 0 and <= 31.
33 : // If this input value is 0 or outside above range
34 : // then the num of bits is calculated based on given
35 : // data and scale factors.
36 : // [4] = Original field type - currently ignored on input
37 : // Data values assumed to be reals.
38 : //
39 : // OUTPUT ARGUMENT LIST:
40 : // idrstmpl - Contains the array of values for Data Representation
41 : // Template 5.0
42 : // [0] = Reference value - set by simpack routine.
43 : // [1] = Binary Scale Factor - unchanged from input
44 : // [2] = Decimal Scale Factor - unchanged from input
45 : // [3] = Number of bits used to pack data, unchanged from
46 : // input if value is between 0 and 31.
47 : // If this input value is 0 or outside above range
48 : // then the num of bits is calculated based on given
49 : // data and scale factors.
50 : // [4] = Original field type - currently set = 0 on output.
51 : // Data values assumed to be reals.
52 : // cpack - The packed data field
53 : // lcpack - length of packed field starting at cpack.
54 : //
55 : // REMARKS: None
56 : //
57 : // ATTRIBUTES:
58 : // LANGUAGE: C
59 : // MACHINE:
60 : //
61 : //$$$
62 : {
63 :
64 65 : const g2int zero=0;
65 : g2int *ifld;
66 : g2int j,nbits,maxdif,nbittot,left;
67 : g2float bscale,dscale,rmax,rmin,temp, ref, rmin_dscaled, rmax_dscaled;
68 : double maxnum;
69 65 : const g2float alog2=0.69314718f; // ln(2.0)
70 :
71 65 : bscale=(float)int_power(2.0,-idrstmpl[1]);
72 65 : dscale=(float)int_power(10.0,idrstmpl[2]);
73 65 : if (idrstmpl[3] <= 0 || idrstmpl[3] > 31)
74 55 : nbits=0;
75 : else
76 10 : nbits=idrstmpl[3];
77 :
78 65 : if( dscale == 0.0 )
79 : {
80 0 : fprintf(stderr, "Invalid dscale == 0 value\n");
81 0 : *lcpack = -1;
82 0 : return;
83 : }
84 : //
85 : // Find max and min values in the data
86 : //
87 65 : rmax=fld[0];
88 65 : rmin=fld[0];
89 60456 : for (j=1;j<ndpts;j++) {
90 60391 : if (fld[j] > rmax) rmax=fld[j];
91 60391 : if (fld[j] < rmin) rmin=fld[j];
92 : }
93 65 : if( !(floor((double)rmin*dscale) >= -FLT_MAX && floor((double)rmin*dscale) <= FLT_MAX) )
94 : {
95 0 : fprintf(stderr,
96 : "Scaled min value not representable on IEEE754 "
97 : "single precision float\n");
98 0 : *lcpack = -1;
99 0 : return;
100 : }
101 65 : if( !(floor((double)rmax*dscale) >= -FLT_MAX && floor((double)rmax*dscale) <= FLT_MAX) )
102 : {
103 0 : fprintf(stderr,
104 : "Scaled max value not representable on IEEE754 "
105 : "single precision float\n");
106 0 : *lcpack = -1;
107 0 : return;
108 : }
109 65 : rmin_dscaled = rmin*dscale;
110 65 : rmax_dscaled = rmax*dscale;
111 :
112 65 : ifld=calloc(ndpts,sizeof(g2int));
113 65 : if( ifld == NULL )
114 : {
115 0 : fprintf(stderr, "Cannot allocate ifld in simpack()\n");
116 0 : *lcpack = -1;
117 0 : return;
118 : }
119 :
120 : //
121 : // If max and min values are not equal, pack up field.
122 : // If they are equal, we have a constant field, and the reference
123 : // value (rmin) is the value for each point in the field and
124 : // set nbits to 0.
125 : //
126 65 : if ( (rmax_dscaled - rmin_dscaled >= 1) ||
127 2 : (rmin != rmax && nbits!=0 && idrstmpl[1]==0) ) {
128 34 : int done = 0;
129 : //
130 : // Determine which algorithm to use based on user-supplied
131 : // binary scale factor and number of bits.
132 : //
133 34 : if (nbits==0 && idrstmpl[1]==0) {
134 : //
135 : // No binary scaling and calculate minimum number of
136 : // bits in which the data will fit.
137 : //
138 24 : if( dscale != 1.0 )
139 : {
140 1 : rmin_dscaled = (float)floor(rmin_dscaled);
141 : }
142 24 : if( (double)(rmax_dscaled - rmin_dscaled) > (double)INT_MAX )
143 : {
144 1 : nbits = 31;
145 : }
146 : else
147 : {
148 23 : temp=(float)(log(ceil(rmax_dscaled - rmin_dscaled))/alog2);
149 23 : nbits=(g2int)ceil(temp);
150 : // scale data
151 23 : if( nbits > 31 )
152 : {
153 0 : nbits = 31;
154 : }
155 : else
156 : {
157 23 : done = 1;
158 9499 : for(j=0;j<ndpts;j++)
159 9476 : ifld[j]=(g2int)RINT(fld[j]*dscale -rmin_dscaled);
160 23 : ref = rmin_dscaled;
161 : }
162 : }
163 : }
164 :
165 34 : if (!done && nbits!=0 && idrstmpl[1]==0) {
166 : //
167 : // Use number of bits specified by user and
168 : // adjust binary scaling factor to accommodate data.
169 : //
170 11 : if( dscale != 1.0 )
171 : {
172 4 : rmin_dscaled = (float)floor(rmin_dscaled);
173 : }
174 11 : maxnum=int_power(2.0,nbits)-1;
175 11 : temp=(float)(log(maxnum/(rmax_dscaled-rmin_dscaled))/alog2);
176 11 : idrstmpl[1]=(g2int)ceil(-1.0*temp);
177 11 : bscale=(float)int_power(2.0,-idrstmpl[1]);
178 : // scale data
179 48089 : for (j=0;j<ndpts;j++)
180 48078 : ifld[j]=(g2int)RINT(((fld[j]*dscale)-rmin_dscaled)*bscale);
181 11 : ref=rmin_dscaled;
182 : }
183 23 : else if (nbits==0 && idrstmpl[1]!=0) {
184 : //
185 : // Use binary scaling factor and calculate minimum number of
186 : // bits in which the data will fit.
187 : //
188 0 : maxdif=(g2int)RINT((rmax_dscaled-rmin_dscaled)*bscale);
189 0 : temp=(float)(log((double)(maxdif+1))/alog2);
190 0 : nbits=(g2int)ceil(temp);
191 : // scale data
192 0 : for (j=0;j<ndpts;j++)
193 0 : ifld[j]=(g2int)RINT(((fld[j]*dscale)-rmin_dscaled)*bscale);
194 0 : ref=rmin_dscaled;
195 : }
196 23 : else if (nbits!=0 && idrstmpl[1]!=0) {
197 : //
198 : // Use binary scaling factor and use minimum number of
199 : // bits specified by user. Dangerous - may loose
200 : // information if binary scale factor and nbits not set
201 : // properly by user.
202 : //
203 : // scale data
204 0 : for (j=0;j<ndpts;j++)
205 0 : ifld[j]=(g2int)RINT(((fld[j]*dscale)-rmin_dscaled)*bscale);
206 0 : ref=rmin_dscaled;
207 : }
208 : //
209 : // Pack data, Pad last octet with Zeros, if necessary,
210 : // and calculate the length of the packed data in bytes
211 : //
212 34 : sbits(cpack,ifld+0,0,nbits,0,ndpts);
213 34 : nbittot=nbits*ndpts;
214 34 : left=8-(nbittot%8);
215 34 : if (left != 8) {
216 3 : sbit(cpack,&zero,nbittot,left); // Pad with zeros to fill Octet
217 3 : nbittot=nbittot+left;
218 : }
219 34 : *lcpack=nbittot/8;
220 : }
221 : else {
222 : /* Force E and D to 0 to avoid compatibility issues */
223 31 : idrstmpl[1]=0;
224 31 : idrstmpl[2]=0;
225 31 : if( dscale != 1.0 )
226 : {
227 0 : ref = (float)floor((double)rmin * dscale) / dscale;
228 : }
229 : else
230 : {
231 31 : ref = rmin;
232 : }
233 31 : nbits=0;
234 31 : *lcpack=0;
235 : }
236 :
237 : //
238 : // Fill in ref value and number of bits in Template 5.0
239 : //
240 : //printf("SAGmkieee %f\n",ref);
241 65 : mkieee(&ref,idrstmpl+0,1); // ensure reference value is IEEE format
242 : //printf("SAGmkieee %ld\n",idrstmpl[0]);
243 65 : idrstmpl[3]=nbits;
244 65 : idrstmpl[4]=0; // original data were reals
245 :
246 65 : free(ifld);
247 : }
|
