Line data Source code
1 : #include <stdlib.h>
2 : #include <math.h>
3 : #include "grib2.h"
4 :
5 :
6 9 : void compack(g2float *fld,g2int ndpts,g2int idrsnum,g2int *idrstmpl,
7 : unsigned char *cpack,g2int *lcpack)
8 : //$$$ SUBPROGRAM DOCUMENTATION BLOCK
9 : // . . . .
10 : // SUBPROGRAM: compack
11 : // PRGMMR: Gilbert ORG: W/NP11 DATE: 2002-11-07
12 : //
13 : // ABSTRACT: This subroutine packs up a data field using a complex
14 : // packing algorithm as defined in the GRIB2 documentation. It
15 : // supports GRIB2 complex packing templates with or without
16 : // spatial differences (i.e. DRTs 5.2 and 5.3).
17 : // It also fills in GRIB2 Data Representation Template 5.2 or 5.3
18 : // with the appropriate values.
19 : //
20 : // PROGRAM HISTORY LOG:
21 : // 2002-11-07 Gilbert
22 : //
23 : // USAGE: void compack(g2float *fld,g2int ndpts,g2int idrsnum,
24 : // g2int *idrstmpl,unsigned char *cpack,g2int *lcpack)
25 : //
26 : // INPUT ARGUMENTS:
27 : // fld[] - Contains the data values to pack
28 : // ndpts - The number of data values in array fld[]
29 : // idrsnum - Data Representation Template number 5.N
30 : // Must equal 2 or 3.
31 : // idrstmpl - Contains the array of values for Data Representation
32 : // Template 5.2 or 5.3
33 : // [0] = Reference value - ignored on input
34 : // [1] = Binary Scale Factor
35 : // [2] = Decimal Scale Factor
36 : // .
37 : // .
38 : // [6] = Missing value management
39 : // [7] = Primary missing value
40 : // [8] = Secondary missing value
41 : // .
42 : // .
43 : // [16] = Order of Spatial Differencing ( 1 or 2 )
44 : // .
45 : // .
46 : //
47 : // OUTPUT ARGUMENTS:
48 : // idrstmpl - Contains the array of values for Data Representation
49 : // Template 5.3
50 : // [0] = Reference value - set by compack routine.
51 : // [1] = Binary Scale Factor - unchanged from input
52 : // [2] = Decimal Scale Factor - unchanged from input
53 : // .
54 : // .
55 : // cpack - The packed data field
56 : // lcpack - length of packed field cpack (or -1 in case of error)
57 : //
58 : // REMARKS: None
59 : //
60 : // ATTRIBUTES:
61 : // LANGUAGE: C
62 : // MACHINE:
63 : //
64 : //$$$
65 : {
66 :
67 9 : const g2int zero=0;
68 : g2int *ifld,*gref,*glen,*gwidth;
69 : g2int *jmin, *jmax, *lbit;
70 : g2int i,j,n, /* nbits, */ imin,imax,left;
71 9 : g2int isd,itemp,ilmax,ngwidthref=0,nbitsgwidth=0;
72 9 : g2int nglenref=0,nglenlast=0,iofst,ival1,ival2=0;
73 9 : g2int minsd,nbitsd=0,maxorig,nbitorig,ngroups;
74 : g2int lg,ng,igmax,iwmax,nbitsgref;
75 9 : g2int glength,grpwidth,nbitsglen=0;
76 : g2int kfildo, minpk, inc, maxgrps, ibit, jbit, kbit, novref, lbitref;
77 9 : g2int missopt, miss1, miss2, ier = 0;
78 : g2float bscale,dscale,rmax,rmin,temp;
79 9 : const g2int simple_alg = 0;
80 9 : const g2float alog2=0.69314718f; // ln(2.0)
81 9 : const g2int one=1;
82 :
83 9 : bscale=(float)int_power(2.0,-idrstmpl[1]);
84 9 : dscale=(float)int_power(10.0,idrstmpl[2]);
85 : //
86 : // Find max and min values in the data
87 : //
88 9 : rmax=fld[0];
89 9 : rmin=fld[0];
90 3204 : for (j=1;j<ndpts;j++) {
91 3195 : if (fld[j] > rmax) rmax=fld[j];
92 3195 : if (fld[j] < rmin) rmin=fld[j];
93 : }
94 :
95 : //
96 : // If max and min values are not equal, pack up field.
97 : // If they are equal, we have a constant field, and the reference
98 : // value (rmin) is the value for each point in the field and
99 : // set nbits to 0.
100 : //
101 9 : if (rmin != rmax) {
102 9 : iofst=0;
103 9 : ifld=calloc(ndpts,sizeof(g2int));
104 9 : gref=calloc(ndpts,sizeof(g2int));
105 9 : gwidth=calloc(ndpts,sizeof(g2int));
106 9 : glen=calloc(ndpts,sizeof(g2int));
107 9 : if( ifld == NULL || gref == NULL || gwidth == NULL || glen == NULL )
108 : {
109 0 : free(ifld);
110 0 : free(gref);
111 0 : free(gwidth);
112 0 : free(glen);
113 0 : *lcpack = -1;
114 1 : return;
115 : }
116 : //
117 : // Scale original data
118 : //
119 9 : if (idrstmpl[1] == 0) { // No binary scaling
120 7 : imin=(g2int)RINT(rmin*dscale);
121 : //imax=(g2int)rint(rmax*dscale);
122 7 : rmin=(g2float)imin;
123 2807 : for (j=0;j<ndpts;j++)
124 2800 : ifld[j]=(g2int)RINT(fld[j]*dscale)-imin;
125 : }
126 : else { // Use binary scaling factor
127 2 : rmin=rmin*dscale;
128 : //rmax=rmax*dscale;
129 406 : for (j=0;j<ndpts;j++)
130 404 : ifld[j]=(g2int)RINT(((fld[j]*dscale)-rmin)*bscale);
131 : }
132 : //
133 : // Calculate spatial differences, if using DRS Template 5.3.
134 : //
135 9 : if (idrsnum == 3) { // spatial differences
136 3 : if (idrstmpl[16]!=1 && idrstmpl[16]!=2) idrstmpl[16]=1;
137 3 : if (idrstmpl[16] == 1) { // first order
138 2 : ival1=ifld[0];
139 800 : for (j=ndpts-1;j>0;j--)
140 798 : ifld[j]=ifld[j]-ifld[j-1];
141 2 : ifld[0]=0;
142 : }
143 1 : else if (idrstmpl[16] == 2) { // second order
144 1 : ival1=ifld[0];
145 1 : ival2=ifld[1];
146 399 : for (j=ndpts-1;j>1;j--)
147 398 : ifld[j]=ifld[j]-(2*ifld[j-1])+ifld[j-2];
148 1 : ifld[0]=0;
149 1 : ifld[1]=0;
150 : }
151 : //
152 : // subtract min value from spatial diff field
153 : //
154 3 : isd=idrstmpl[16];
155 3 : minsd=ifld[isd];
156 1199 : for (j=isd;j<ndpts;j++) if ( ifld[j] < minsd ) minsd=ifld[j];
157 1199 : for (j=isd;j<ndpts;j++) ifld[j]=ifld[j]-minsd;
158 : //
159 : // find num of bits need to store minsd and add 1 extra bit
160 : // to indicate sign
161 : //
162 3 : temp=(float)(log((double)(abs(minsd)+1))/alog2);
163 3 : nbitsd=(g2int)ceil(temp)+1;
164 : //
165 : // find num of bits need to store ifld[0] ( and ifld[1]
166 : // if using 2nd order differencing )
167 : //
168 3 : maxorig=ival1;
169 3 : if (idrstmpl[16]==2 && ival2>ival1) maxorig=ival2;
170 3 : temp=(float)(log((double)(maxorig+1))/alog2);
171 3 : nbitorig=(g2int)ceil(temp)+1;
172 3 : if (nbitorig > nbitsd) nbitsd=nbitorig;
173 : // increase number of bits to even multiple of 8 ( octet )
174 3 : if ( (nbitsd%8) != 0) nbitsd=nbitsd+(8-(nbitsd%8));
175 : //
176 : // Store extra spatial differencing info into the packed
177 : // data section.
178 : //
179 3 : if (nbitsd != 0) {
180 : // pack first original value
181 3 : if (ival1 >= 0) {
182 3 : sbit(cpack,&ival1,iofst,nbitsd);
183 3 : iofst=iofst+nbitsd;
184 : }
185 : else {
186 0 : sbit(cpack,&one,iofst,1);
187 0 : iofst=iofst+1;
188 0 : itemp=abs(ival1);
189 0 : sbit(cpack,&itemp,iofst,nbitsd-1);
190 0 : iofst=iofst+nbitsd-1;
191 : }
192 3 : if (idrstmpl[16] == 2) {
193 : // pack second original value
194 1 : if (ival2 >= 0) {
195 1 : sbit(cpack,&ival2,iofst,nbitsd);
196 1 : iofst=iofst+nbitsd;
197 : }
198 : else {
199 0 : sbit(cpack,&one,iofst,1);
200 0 : iofst=iofst+1;
201 0 : itemp=abs(ival2);
202 0 : sbit(cpack,&itemp,iofst,nbitsd-1);
203 0 : iofst=iofst+nbitsd-1;
204 : }
205 : }
206 : // pack overall min of spatial differences
207 3 : if (minsd >= 0) {
208 0 : sbit(cpack,&minsd,iofst,nbitsd);
209 0 : iofst=iofst+nbitsd;
210 : }
211 : else {
212 3 : sbit(cpack,&one,iofst,1);
213 3 : iofst=iofst+1;
214 3 : itemp=abs(minsd);
215 3 : sbit(cpack,&itemp,iofst,nbitsd-1);
216 3 : iofst=iofst+nbitsd-1;
217 : }
218 : }
219 : //printf("SDp %ld %ld %ld %ld\n",ival1,ival2,minsd,nbitsd);
220 : } // end of spatial diff section
221 : //
222 : // Determine Groups to be used.
223 : //
224 9 : if ( simple_alg == 1 ) {
225 : // set group length to 10; calculate number of groups
226 : // and length of last group
227 0 : ngroups=ndpts/10;
228 0 : for (j=0;j<ngroups;j++) glen[j]=10;
229 0 : itemp=ndpts%10;
230 0 : if (itemp != 0) {
231 0 : ngroups=ngroups+1;
232 0 : glen[ngroups-1]=itemp;
233 : }
234 : }
235 : else {
236 : // Use Dr. Glahn's algorithm for determining grouping.
237 : //
238 9 : kfildo=6;
239 9 : minpk=10;
240 9 : inc=1;
241 9 : maxgrps=(ndpts/minpk)+1;
242 9 : jmin = calloc(maxgrps,sizeof(g2int));
243 9 : jmax = calloc(maxgrps,sizeof(g2int));
244 9 : lbit = calloc(maxgrps,sizeof(g2int));
245 9 : if( jmin == NULL || jmax == NULL || lbit == NULL )
246 : {
247 0 : free(jmin);
248 0 : free(jmax);
249 0 : free(lbit);
250 :
251 0 : free(ifld);
252 0 : free(gref);
253 0 : free(gwidth);
254 0 : free(glen);
255 0 : *lcpack = -1;
256 0 : return;
257 : }
258 9 : missopt=0;
259 9 : pack_gp(&kfildo,ifld,&ndpts,&missopt,&minpk,&inc,&miss1,&miss2,
260 : jmin,jmax,lbit,glen,&maxgrps,&ngroups,&ibit,&jbit,
261 : &kbit,&novref,&lbitref,&ier);
262 : //print *,'SAGier = ',ier,ibit,jbit,kbit,novref,lbitref
263 225 : for ( ng=0; ng<ngroups; ng++) glen[ng]=glen[ng]+novref;
264 9 : free(jmin);
265 9 : free(jmax);
266 9 : free(lbit);
267 9 : if( ier != 0 )
268 : {
269 1 : free(ifld);
270 1 : free(gref);
271 1 : free(gwidth);
272 1 : free(glen);
273 1 : *lcpack = -1;
274 1 : return;
275 : }
276 : }
277 : //
278 : // For each group, find the group's reference value
279 : // and the number of bits needed to hold the remaining values
280 : //
281 8 : n=0;
282 223 : for (ng=0;ng<ngroups;ng++) {
283 : // find max and min values of group
284 215 : gref[ng]=ifld[n];
285 215 : imax=ifld[n];
286 215 : j=n+1;
287 3200 : for (lg=1;lg<glen[ng];lg++) {
288 2985 : if (ifld[j] < gref[ng]) gref[ng]=ifld[j];
289 2985 : if (ifld[j] > imax) imax=ifld[j];
290 2985 : j++;
291 : }
292 : // calc num of bits needed to hold data
293 215 : if ( gref[ng] != imax ) {
294 198 : temp=(float)(log((double)(imax-gref[ng]+1))/alog2);
295 198 : gwidth[ng]=(g2int)ceil(temp);
296 : }
297 : else
298 17 : gwidth[ng]=0;
299 : // Subtract min from data
300 215 : j=n;
301 3415 : for (lg=0;lg<glen[ng];lg++) {
302 3200 : ifld[j]=ifld[j]-gref[ng];
303 3200 : j++;
304 : }
305 : // increment fld array counter
306 215 : n=n+glen[ng];
307 : }
308 : //
309 : // Find max of the group references and calc num of bits needed
310 : // to pack each groups reference value, then
311 : // pack up group reference values
312 : //
313 8 : igmax=gref[0];
314 215 : for (j=1;j<ngroups;j++) if (gref[j] > igmax) igmax=gref[j];
315 8 : if (igmax != 0) {
316 8 : temp=(float)(log((double)(igmax+1))/alog2);
317 8 : nbitsgref=(g2int)ceil(temp);
318 8 : sbits(cpack,gref,iofst,nbitsgref,0,ngroups);
319 8 : itemp=nbitsgref*ngroups;
320 8 : iofst=iofst+itemp;
321 : // Pad last octet with Zeros, if necessary,
322 8 : if ( (itemp%8) != 0) {
323 2 : left=8-(itemp%8);
324 2 : sbit(cpack,&zero,iofst,left);
325 2 : iofst=iofst+left;
326 : }
327 : }
328 : else
329 0 : nbitsgref=0;
330 : //
331 : // Find max/min of the group widths and calc num of bits needed
332 : // to pack each groups width value, then
333 : // pack up group width values
334 : //
335 8 : iwmax=gwidth[0];
336 8 : ngwidthref=gwidth[0];
337 215 : for (j=1;j<ngroups;j++) {
338 207 : if (gwidth[j] > iwmax) iwmax=gwidth[j];
339 207 : if (gwidth[j] < ngwidthref) ngwidthref=gwidth[j];
340 : }
341 8 : if (iwmax != ngwidthref) {
342 8 : temp=(float)(log((double)(iwmax-ngwidthref+1))/alog2);
343 8 : nbitsgwidth=(g2int)ceil(temp);
344 223 : for (i=0;i<ngroups;i++)
345 215 : gwidth[i]=gwidth[i]-ngwidthref;
346 8 : sbits(cpack,gwidth,iofst,nbitsgwidth,0,ngroups);
347 8 : itemp=nbitsgwidth*ngroups;
348 8 : iofst=iofst+itemp;
349 : // Pad last octet with Zeros, if necessary,
350 8 : if ( (itemp%8) != 0) {
351 7 : left=8-(itemp%8);
352 7 : sbit(cpack,&zero,iofst,left);
353 7 : iofst=iofst+left;
354 : }
355 : }
356 : else {
357 0 : nbitsgwidth=0;
358 0 : for (i=0;i<ngroups;i++) gwidth[i]=0;
359 : }
360 : //
361 : // Find max/min of the group lengths and calc num of bits needed
362 : // to pack each groups length value, then
363 : // pack up group length values
364 : //
365 : //write(77,*)'GLENS: ',(glen(j),j=1,ngroups)
366 8 : ilmax=glen[0];
367 8 : nglenref=glen[0];
368 207 : for (j=1;j<ngroups-1;j++) {
369 199 : if (glen[j] > ilmax) ilmax=glen[j];
370 199 : if (glen[j] < nglenref) nglenref=glen[j];
371 : }
372 8 : nglenlast=glen[ngroups-1];
373 8 : if (ilmax != nglenref) {
374 8 : temp=(float)(log((double)(ilmax-nglenref+1))/alog2);
375 8 : nbitsglen=(g2int)ceil(temp);
376 215 : for (i=0;i<ngroups-1;i++) glen[i]=glen[i]-nglenref;
377 8 : sbits(cpack,glen,iofst,nbitsglen,0,ngroups);
378 8 : itemp=nbitsglen*ngroups;
379 8 : iofst=iofst+itemp;
380 : // Pad last octet with Zeros, if necessary,
381 8 : if ( (itemp%8) != 0) {
382 7 : left=8-(itemp%8);
383 7 : sbit(cpack,&zero,iofst,left);
384 7 : iofst=iofst+left;
385 : }
386 : }
387 : else {
388 0 : nbitsglen=0;
389 0 : for (i=0;i<ngroups;i++) glen[i]=0;
390 : }
391 : //
392 : // For each group, pack data values
393 : //
394 8 : n=0;
395 223 : for (ng=0;ng<ngroups;ng++) {
396 215 : glength=glen[ng]+nglenref;
397 215 : if (ng == (ngroups-1) ) glength=nglenlast;
398 215 : grpwidth=gwidth[ng]+ngwidthref;
399 215 : if ( grpwidth != 0 ) {
400 198 : sbits(cpack,ifld+n,iofst,grpwidth,0,glength);
401 198 : iofst=iofst+(grpwidth*glength);
402 : }
403 215 : n=n+glength;
404 : }
405 : // Pad last octet with Zeros, if necessary,
406 8 : if ( (iofst%8) != 0) {
407 7 : left=8-(iofst%8);
408 7 : sbit(cpack,&zero,iofst,left);
409 7 : iofst=iofst+left;
410 : }
411 8 : *lcpack=iofst/8;
412 : //
413 8 : free(ifld);
414 8 : free(gref);
415 8 : free(gwidth);
416 8 : free(glen);
417 : }
418 : else { // Constant field ( max = min )
419 : /* nbits=0; */
420 0 : *lcpack=0;
421 0 : nbitsgref=0;
422 0 : ngroups=0;
423 : }
424 :
425 : //
426 : // Fill in ref value and number of bits in Template 5.2
427 : //
428 8 : mkieee(&rmin,idrstmpl+0,1); // ensure reference value is IEEE format
429 8 : idrstmpl[3]=nbitsgref;
430 8 : idrstmpl[4]=0; // original data were reals
431 8 : idrstmpl[5]=1; // general group splitting
432 8 : idrstmpl[6]=0; // No internal missing values
433 8 : idrstmpl[7]=0; // Primary missing value
434 8 : idrstmpl[8]=0; // secondary missing value
435 8 : idrstmpl[9]=ngroups; // Number of groups
436 8 : idrstmpl[10]=ngwidthref; // reference for group widths
437 8 : idrstmpl[11]=nbitsgwidth; // num bits used for group widths
438 8 : idrstmpl[12]=nglenref; // Reference for group lengths
439 8 : idrstmpl[13]=1; // length increment for group lengths
440 8 : idrstmpl[14]=nglenlast; // True length of last group
441 8 : idrstmpl[15]=nbitsglen; // num bits used for group lengths
442 8 : if (idrsnum == 3) {
443 3 : idrstmpl[17]=nbitsd/8; // num bits used for extra spatial
444 : // differencing values
445 : }
446 :
447 : }
|
