1 #include <stdlib.h>
2 #include <math.h>
3 #include "grib2.h"
4
5
compack(g2float * fld,g2int ndpts,g2int idrsnum,g2int * idrstmpl,unsigned char * cpack,g2int * lcpack)6 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 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 g2int isd,itemp,ilmax,ngwidthref=0,nbitsgwidth=0;
72 g2int nglenref=0,nglenlast=0,iofst,ival1,ival2=0;
73 g2int minsd,nbitsd=0,maxorig,nbitorig,ngroups;
74 g2int lg,ng,igmax,iwmax,nbitsgref;
75 g2int glength,grpwidth,nbitsglen=0;
76 g2int kfildo, minpk, inc, maxgrps, ibit, jbit, kbit, novref, lbitref;
77 g2int missopt, miss1, miss2, ier = 0;
78 g2float bscale,dscale,rmax,rmin,temp;
79 const g2int simple_alg = 0;
80 const g2float alog2=0.69314718f; // ln(2.0)
81 const g2int one=1;
82
83 bscale=(float)int_power(2.0,-idrstmpl[1]);
84 dscale=(float)int_power(10.0,idrstmpl[2]);
85 //
86 // Find max and min values in the data
87 //
88 rmax=fld[0];
89 rmin=fld[0];
90 for (j=1;j<ndpts;j++) {
91 if (fld[j] > rmax) rmax=fld[j];
92 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 if (rmin != rmax) {
102 iofst=0;
103 ifld=calloc(ndpts,sizeof(g2int));
104 gref=calloc(ndpts,sizeof(g2int));
105 gwidth=calloc(ndpts,sizeof(g2int));
106 glen=calloc(ndpts,sizeof(g2int));
107 if( ifld == NULL || gref == NULL || gwidth == NULL || glen == NULL )
108 {
109 free(ifld);
110 free(gref);
111 free(gwidth);
112 free(glen);
113 *lcpack = -1;
114 return;
115 }
116 //
117 // Scale original data
118 //
119 if (idrstmpl[1] == 0) { // No binary scaling
120 imin=(g2int)RINT(rmin*dscale);
121 //imax=(g2int)rint(rmax*dscale);
122 rmin=(g2float)imin;
123 for (j=0;j<ndpts;j++)
124 ifld[j]=(g2int)RINT(fld[j]*dscale)-imin;
125 }
126 else { // Use binary scaling factor
127 rmin=rmin*dscale;
128 //rmax=rmax*dscale;
129 for (j=0;j<ndpts;j++)
130 ifld[j]=(g2int)RINT(((fld[j]*dscale)-rmin)*bscale);
131 }
132 //
133 // Calculate spatial differences, if using DRS Template 5.3.
134 //
135 if (idrsnum == 3) { // spatial differences
136 if (idrstmpl[16]!=1 && idrstmpl[16]!=2) idrstmpl[16]=1;
137 if (idrstmpl[16] == 1) { // first order
138 ival1=ifld[0];
139 for (j=ndpts-1;j>0;j--)
140 ifld[j]=ifld[j]-ifld[j-1];
141 ifld[0]=0;
142 }
143 else if (idrstmpl[16] == 2) { // second order
144 ival1=ifld[0];
145 ival2=ifld[1];
146 for (j=ndpts-1;j>1;j--)
147 ifld[j]=ifld[j]-(2*ifld[j-1])+ifld[j-2];
148 ifld[0]=0;
149 ifld[1]=0;
150 }
151 //
152 // subtract min value from spatial diff field
153 //
154 isd=idrstmpl[16];
155 minsd=ifld[isd];
156 for (j=isd;j<ndpts;j++) if ( ifld[j] < minsd ) minsd=ifld[j];
157 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 temp=(float)(log((double)(abs(minsd)+1))/alog2);
163 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 maxorig=ival1;
169 if (idrstmpl[16]==2 && ival2>ival1) maxorig=ival2;
170 temp=(float)(log((double)(maxorig+1))/alog2);
171 nbitorig=(g2int)ceil(temp)+1;
172 if (nbitorig > nbitsd) nbitsd=nbitorig;
173 // increase number of bits to even multiple of 8 ( octet )
174 if ( (nbitsd%8) != 0) nbitsd=nbitsd+(8-(nbitsd%8));
175 //
176 // Store extra spatial differencing info into the packed
177 // data section.
178 //
179 if (nbitsd != 0) {
180 // pack first original value
181 if (ival1 >= 0) {
182 sbit(cpack,&ival1,iofst,nbitsd);
183 iofst=iofst+nbitsd;
184 }
185 else {
186 sbit(cpack,&one,iofst,1);
187 iofst=iofst+1;
188 itemp=abs(ival1);
189 sbit(cpack,&itemp,iofst,nbitsd-1);
190 iofst=iofst+nbitsd-1;
191 }
192 if (idrstmpl[16] == 2) {
193 // pack second original value
194 if (ival2 >= 0) {
195 sbit(cpack,&ival2,iofst,nbitsd);
196 iofst=iofst+nbitsd;
197 }
198 else {
199 sbit(cpack,&one,iofst,1);
200 iofst=iofst+1;
201 itemp=abs(ival2);
202 sbit(cpack,&itemp,iofst,nbitsd-1);
203 iofst=iofst+nbitsd-1;
204 }
205 }
206 // pack overall min of spatial differences
207 if (minsd >= 0) {
208 sbit(cpack,&minsd,iofst,nbitsd);
209 iofst=iofst+nbitsd;
210 }
211 else {
212 sbit(cpack,&one,iofst,1);
213 iofst=iofst+1;
214 itemp=abs(minsd);
215 sbit(cpack,&itemp,iofst,nbitsd-1);
216 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 if ( simple_alg == 1 ) {
225 // set group length to 10; calculate number of groups
226 // and length of last group
227 ngroups=ndpts/10;
228 for (j=0;j<ngroups;j++) glen[j]=10;
229 itemp=ndpts%10;
230 if (itemp != 0) {
231 ngroups=ngroups+1;
232 glen[ngroups-1]=itemp;
233 }
234 }
235 else {
236 // Use Dr. Glahn's algorithm for determining grouping.
237 //
238 kfildo=6;
239 minpk=10;
240 inc=1;
241 maxgrps=(ndpts/minpk)+1;
242 jmin = calloc(maxgrps,sizeof(g2int));
243 jmax = calloc(maxgrps,sizeof(g2int));
244 lbit = calloc(maxgrps,sizeof(g2int));
245 if( jmin == NULL || jmax == NULL || lbit == NULL )
246 {
247 free(jmin);
248 free(jmax);
249 free(lbit);
250
251 free(ifld);
252 free(gref);
253 free(gwidth);
254 free(glen);
255 *lcpack = -1;
256 return;
257 }
258 missopt=0;
259 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 for ( ng=0; ng<ngroups; ng++) glen[ng]=glen[ng]+novref;
264 free(jmin);
265 free(jmax);
266 free(lbit);
267 if( ier != 0 )
268 {
269 free(ifld);
270 free(gref);
271 free(gwidth);
272 free(glen);
273 *lcpack = -1;
274 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 n=0;
282 for (ng=0;ng<ngroups;ng++) {
283 // find max and min values of group
284 gref[ng]=ifld[n];
285 imax=ifld[n];
286 j=n+1;
287 for (lg=1;lg<glen[ng];lg++) {
288 if (ifld[j] < gref[ng]) gref[ng]=ifld[j];
289 if (ifld[j] > imax) imax=ifld[j];
290 j++;
291 }
292 // calc num of bits needed to hold data
293 if ( gref[ng] != imax ) {
294 temp=(float)(log((double)(imax-gref[ng]+1))/alog2);
295 gwidth[ng]=(g2int)ceil(temp);
296 }
297 else
298 gwidth[ng]=0;
299 // Subtract min from data
300 j=n;
301 for (lg=0;lg<glen[ng];lg++) {
302 ifld[j]=ifld[j]-gref[ng];
303 j++;
304 }
305 // increment fld array counter
306 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 igmax=gref[0];
314 for (j=1;j<ngroups;j++) if (gref[j] > igmax) igmax=gref[j];
315 if (igmax != 0) {
316 temp=(float)(log((double)(igmax+1))/alog2);
317 nbitsgref=(g2int)ceil(temp);
318 sbits(cpack,gref,iofst,nbitsgref,0,ngroups);
319 itemp=nbitsgref*ngroups;
320 iofst=iofst+itemp;
321 // Pad last octet with Zeros, if necessary,
322 if ( (itemp%8) != 0) {
323 left=8-(itemp%8);
324 sbit(cpack,&zero,iofst,left);
325 iofst=iofst+left;
326 }
327 }
328 else
329 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 iwmax=gwidth[0];
336 ngwidthref=gwidth[0];
337 for (j=1;j<ngroups;j++) {
338 if (gwidth[j] > iwmax) iwmax=gwidth[j];
339 if (gwidth[j] < ngwidthref) ngwidthref=gwidth[j];
340 }
341 if (iwmax != ngwidthref) {
342 temp=(float)(log((double)(iwmax-ngwidthref+1))/alog2);
343 nbitsgwidth=(g2int)ceil(temp);
344 for (i=0;i<ngroups;i++)
345 gwidth[i]=gwidth[i]-ngwidthref;
346 sbits(cpack,gwidth,iofst,nbitsgwidth,0,ngroups);
347 itemp=nbitsgwidth*ngroups;
348 iofst=iofst+itemp;
349 // Pad last octet with Zeros, if necessary,
350 if ( (itemp%8) != 0) {
351 left=8-(itemp%8);
352 sbit(cpack,&zero,iofst,left);
353 iofst=iofst+left;
354 }
355 }
356 else {
357 nbitsgwidth=0;
358 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 ilmax=glen[0];
367 nglenref=glen[0];
368 for (j=1;j<ngroups-1;j++) {
369 if (glen[j] > ilmax) ilmax=glen[j];
370 if (glen[j] < nglenref) nglenref=glen[j];
371 }
372 nglenlast=glen[ngroups-1];
373 if (ilmax != nglenref) {
374 temp=(float)(log((double)(ilmax-nglenref+1))/alog2);
375 nbitsglen=(g2int)ceil(temp);
376 for (i=0;i<ngroups-1;i++) glen[i]=glen[i]-nglenref;
377 sbits(cpack,glen,iofst,nbitsglen,0,ngroups);
378 itemp=nbitsglen*ngroups;
379 iofst=iofst+itemp;
380 // Pad last octet with Zeros, if necessary,
381 if ( (itemp%8) != 0) {
382 left=8-(itemp%8);
383 sbit(cpack,&zero,iofst,left);
384 iofst=iofst+left;
385 }
386 }
387 else {
388 nbitsglen=0;
389 for (i=0;i<ngroups;i++) glen[i]=0;
390 }
391 //
392 // For each group, pack data values
393 //
394 n=0;
395 for (ng=0;ng<ngroups;ng++) {
396 glength=glen[ng]+nglenref;
397 if (ng == (ngroups-1) ) glength=nglenlast;
398 grpwidth=gwidth[ng]+ngwidthref;
399 if ( grpwidth != 0 ) {
400 sbits(cpack,ifld+n,iofst,grpwidth,0,glength);
401 iofst=iofst+(grpwidth*glength);
402 }
403 n=n+glength;
404 }
405 // Pad last octet with Zeros, if necessary,
406 if ( (iofst%8) != 0) {
407 left=8-(iofst%8);
408 sbit(cpack,&zero,iofst,left);
409 iofst=iofst+left;
410 }
411 *lcpack=iofst/8;
412 //
413 free(ifld);
414 free(gref);
415 free(gwidth);
416 free(glen);
417 }
418 else { // Constant field ( max = min )
419 /* nbits=0; */
420 *lcpack=0;
421 nbitsgref=0;
422 ngroups=0;
423 }
424
425 //
426 // Fill in ref value and number of bits in Template 5.2
427 //
428 mkieee(&rmin,idrstmpl+0,1); // ensure reference value is IEEE format
429 idrstmpl[3]=nbitsgref;
430 idrstmpl[4]=0; // original data were reals
431 idrstmpl[5]=1; // general group splitting
432 idrstmpl[6]=0; // No internal missing values
433 idrstmpl[7]=0; // Primary missing value
434 idrstmpl[8]=0; // secondary missing value
435 idrstmpl[9]=ngroups; // Number of groups
436 idrstmpl[10]=ngwidthref; // reference for group widths
437 idrstmpl[11]=nbitsgwidth; // num bits used for group widths
438 idrstmpl[12]=nglenref; // Reference for group lengths
439 idrstmpl[13]=1; // length increment for group lengths
440 idrstmpl[14]=nglenlast; // True length of last group
441 idrstmpl[15]=nbitsglen; // num bits used for group lengths
442 if (idrsnum == 3) {
443 idrstmpl[17]=nbitsd/8; // num bits used for extra spatial
444 // differencing values
445 }
446
447 }
448