1 #include <stdlib.h>
2 #include <math.h>
3 #include <limits.h>
4 #include <float.h>
5 #include "grib2.h"
6
misspack(g2float * fld,g2int ndpts,g2int idrsnum,g2int * idrstmpl,unsigned char * cpack,g2int * lcpack)7 void misspack(g2float *fld,g2int ndpts,g2int idrsnum,g2int *idrstmpl,
8 unsigned char *cpack, g2int *lcpack)
9 //$$$ SUBPROGRAM DOCUMENTATION BLOCK
10 // . . . .
11 // SUBPROGRAM: misspack
12 // PRGMMR: Gilbert ORG: W/NP11 DATE: 2000-06-21
13 //
14 // ABSTRACT: This subroutine packs up a data field using a complex
15 // packing algorithm as defined in the GRIB2 documentation. It
16 // supports GRIB2 complex packing templates with or without
17 // spatial differences (i.e. DRTs 5.2 and 5.3).
18 // It also fills in GRIB2 Data Representation Template 5.2 or 5.3
19 // with the appropriate values.
20 // This version assumes that Missing Value Management is being used and that
21 // 1 or 2 missing values appear in the data.
22 //
23 // PROGRAM HISTORY LOG:
24 // 2000-06-21 Gilbert
25 //
26 // USAGE: misspack(g2float *fld,g2int ndpts,g2int idrsnum,g2int *idrstmpl,
27 // unsigned char *cpack, g2int *lcpack)
28 // INPUT ARGUMENT LIST:
29 // fld[] - Contains the data values to pack
30 // ndpts - The number of data values in array fld[]
31 // idrsnum - Data Representation Template number 5.N
32 // Must equal 2 or 3.
33 // idrstmpl - Contains the array of values for Data Representation
34 // Template 5.2 or 5.3
35 // [0] = Reference value - ignored on input
36 // [1] = Binary Scale Factor
37 // [2] = Decimal Scale Factor
38 // .
39 // .
40 // [6] = Missing value management
41 // [7] = Primary missing value
42 // [8] = Secondary missing value
43 // .
44 // .
45 // [16] = Order of Spatial Differencing ( 1 or 2 )
46 // .
47 // .
48 //
49 // OUTPUT ARGUMENT LIST:
50 // idrstmpl - Contains the array of values for Data Representation
51 // Template 5.3
52 // [0] = Reference value - set by misspack routine.
53 // [1] = Binary Scale Factor - unchanged from input
54 // [2] = Decimal Scale Factor - unchanged from input
55 // .
56 // .
57 // cpack - The packed data field (character*1 array)
58 // *lcpack - length of packed field cpack(). (or -1 in case of error)
59 //
60 // REMARKS: None
61 //
62 // ATTRIBUTES:
63 // LANGUAGE: C
64 // MACHINE:
65 //
66 //$$$
67 {
68
69 g2int *ifld, *ifldmiss, *jfld;
70 g2int *jmin, *jmax, *lbit;
71 const g2int zero=0;
72 g2int *gref, *gwidth, *glen;
73 g2int glength, grpwidth;
74 g2int i, n, iofst, ival1, ival2, isd, minsd, nbitsd = 0;
75 g2int nbitsgref, left, iwmax, ngwidthref, nbitsgwidth, ilmax;
76 g2int nglenref, nglenlast, nbitsglen /* , ij */;
77 g2int j, missopt, nonmiss, itemp, maxorig, nbitorig, miss1, miss2;
78 g2int ngroups, ng, num0, num1, num2;
79 g2int imax, lg, mtemp, ier = 0, igmax;
80 g2int kfildo, minpk, inc, maxgrps, ibit, jbit, kbit, novref, lbitref;
81 g2float rmissp, rmisss, bscale, dscale, rmin, rmax, temp;
82 const g2int simple_alg = 0;
83 const g2float alog2=0.69314718f; // ln(2.0)
84 const g2int one=1;
85
86 bscale=(float)int_power(2.0,-idrstmpl[1]);
87 dscale=(float)int_power(10.0,idrstmpl[2]);
88 missopt=idrstmpl[6];
89 if ( missopt != 1 && missopt != 2 ) {
90 printf("misspack: Unrecognized option.\n");
91 *lcpack=-1;
92 return;
93 }
94 else { // Get missing values
95 rdieee(idrstmpl+7,&rmissp,1);
96 if (missopt == 2) rdieee(idrstmpl+8,&rmisss,1);
97 }
98 //
99 // Find min value of non-missing values in the data,
100 // AND set up missing value mapping of the field.
101 //
102 ifldmiss = calloc(ndpts,sizeof(g2int));
103 if( ifldmiss == NULL )
104 {
105 *lcpack = -1;
106 return;
107 }
108 rmin=FLT_MAX;
109 rmax=-FLT_MAX;
110 if ( missopt == 1 ) { // Primary missing value only
111 for ( j=0; j<ndpts; j++) {
112 if (fld[j] == rmissp) {
113 ifldmiss[j]=1;
114 }
115 else {
116 ifldmiss[j]=0;
117 if (fld[j] < rmin) rmin=fld[j];
118 if (fld[j] > rmax) rmax=fld[j];
119 }
120 }
121 }
122 if ( missopt == 2 ) { // Primary and secondary missing values
123 for ( j=0; j<ndpts; j++ ) {
124 if (fld[j] == rmissp) {
125 ifldmiss[j]=1;
126 }
127 else if (fld[j] == rmisss) {
128 ifldmiss[j]=2;
129 }
130 else {
131 ifldmiss[j]=0;
132 if (fld[j] < rmin) rmin=fld[j];
133 if (fld[j] > rmax) rmax=fld[j];
134 }
135 }
136 }
137
138 if( !(floor(rmin*dscale) >= -FLT_MAX && floor(rmin*dscale) <= FLT_MAX) )
139 {
140 fprintf(stderr,
141 "Scaled min value not representable on IEEE754 "
142 "single precision float\n");
143 *lcpack = -1;
144 free(ifldmiss);
145 return;
146 }
147 if (idrstmpl[1] == 0)
148 {
149 double max_diff = RINT(rmax*dscale) - RINT(rmin*dscale);
150 if( !(max_diff >= 0 && max_diff <= INT_MAX) )
151 {
152 fprintf(stderr,
153 "Difference of scaled max value with scaled min value "
154 "not representable on int \n");
155 *lcpack = -1;
156 free(ifldmiss);
157 return;
158 }
159 }
160 else
161 {
162 double max_diff = RINT(rmax*dscale - rmin*dscale) * bscale;
163 if( !(max_diff >= 0 && max_diff <= INT_MAX) )
164 {
165 fprintf(stderr,
166 "Difference of scaled max value with scaled min value "
167 "not representable on int \n");
168 *lcpack = -1;
169 free(ifldmiss);
170 return;
171 }
172 }
173 //
174 // Allocate work arrays:
175 // Note: -ifldmiss[j],j=0,ndpts-1 is a map of original field indicating
176 // which of the original data values
177 // are primary missing (1), secondary missing (2) or non-missing (0).
178 // -jfld[j],j=0,nonmiss-1 is a subarray of just the non-missing values
179 // from the original field.
180 //
181 //if (rmin != rmax) {
182 iofst=0;
183 ifld = calloc(ndpts,sizeof(g2int));
184 jfld = calloc(ndpts,sizeof(g2int));
185 gref = calloc(ndpts,sizeof(g2int));
186 gwidth = calloc(ndpts,sizeof(g2int));
187 glen = calloc(ndpts,sizeof(g2int));
188 if( ifld == NULL || jfld == NULL || gref == NULL || gwidth == NULL ||
189 glen == NULL )
190 {
191 free(ifld);
192 free(jfld);
193 free(ifldmiss);
194 free(gref);
195 free(gwidth);
196 free(glen);
197 *lcpack = -1;
198 return;
199 }
200 //
201 // Scale original data
202 //
203 nonmiss=0;
204 if (idrstmpl[1] == 0) { // No binary scaling
205 rmin=(g2float)RINT(rmin*dscale);
206 for ( j=0; j<ndpts; j++) {
207 if (ifldmiss[j] == 0) {
208 jfld[nonmiss]=(g2int)(RINT(fld[j]*dscale)-rmin);
209 nonmiss++;
210 }
211 }
212 }
213 else { // Use binary scaling factor
214 rmin=rmin*dscale;
215 //rmax=rmax*dscale;
216 for ( j=0; j<ndpts; j++ ) {
217 if (ifldmiss[j] == 0) {
218 jfld[nonmiss]=(g2int)RINT(((fld[j]*dscale)-rmin)*bscale);
219 nonmiss++;
220 }
221 }
222 }
223 //
224 // Calculate Spatial differences, if using DRS Template 5.3
225 //
226 if (idrsnum == 3) { // spatial differences
227 if (idrstmpl[16]!=1 && idrstmpl[16]!=2) idrstmpl[16]=2;
228 if (idrstmpl[16] == 1) { // first order
229 ival1=jfld[0];
230 for ( j=nonmiss-1; j>0; j--)
231 jfld[j]=jfld[j]-jfld[j-1];
232 jfld[0]=0;
233 }
234 else if (idrstmpl[16] == 2) { // second order
235 ival1=jfld[0];
236 ival2=jfld[1];
237 for ( j=nonmiss-1; j>1; j--)
238 jfld[j]=jfld[j]-(2*jfld[j-1])+jfld[j-2];
239 jfld[0]=0;
240 jfld[1]=0;
241 }
242 //
243 // subtract min value from spatial diff field
244 //
245 isd=idrstmpl[16];
246 minsd=jfld[isd];
247 for ( j=isd; j<nonmiss; j++ ) if ( jfld[j] < minsd ) minsd=jfld[j];
248 for ( j=isd; j<nonmiss; j++ ) jfld[j]=jfld[j]-minsd;
249 //
250 // find num of bits need to store minsd and add 1 extra bit
251 // to indicate sign
252 //
253 temp=(float)(log((double)(abs(minsd)+1))/alog2);
254 nbitsd=(g2int)ceil(temp)+1;
255 //
256 // find num of bits need to store ifld[0] ( and ifld[1]
257 // if using 2nd order differencing )
258 //
259 maxorig=ival1;
260 if (idrstmpl[16]==2 && ival2>ival1) maxorig=ival2;
261 temp=(float)(log((double)(maxorig+1))/alog2);
262 nbitorig=(g2int)ceil(temp)+1;
263 if (nbitorig > nbitsd) nbitsd=nbitorig;
264 // increase number of bits to even multiple of 8 ( octet )
265 if ( (nbitsd%8) != 0) nbitsd=nbitsd+(8-(nbitsd%8));
266 //
267 // Store extra spatial differencing info into the packed
268 // data section.
269 //
270 if (nbitsd != 0) {
271 // pack first original value
272 if (ival1 >= 0) {
273 sbit(cpack,&ival1,iofst,nbitsd);
274 iofst=iofst+nbitsd;
275 }
276 else {
277 sbit(cpack,&one,iofst,1);
278 iofst=iofst+1;
279 itemp=abs(ival1);
280 sbit(cpack,&itemp,iofst,nbitsd-1);
281 iofst=iofst+nbitsd-1;
282 }
283 if (idrstmpl[16] == 2) {
284 // pack second original value
285 if (ival2 >= 0) {
286 sbit(cpack,&ival2,iofst,nbitsd);
287 iofst=iofst+nbitsd;
288 }
289 else {
290 sbit(cpack,&one,iofst,1);
291 iofst=iofst+1;
292 itemp=abs(ival2);
293 sbit(cpack,&itemp,iofst,nbitsd-1);
294 iofst=iofst+nbitsd-1;
295 }
296 }
297 // pack overall min of spatial differences
298 if (minsd >= 0) {
299 sbit(cpack,&minsd,iofst,nbitsd);
300 iofst=iofst+nbitsd;
301 }
302 else {
303 sbit(cpack,&one,iofst,1);
304 iofst=iofst+1;
305 itemp=abs(minsd);
306 sbit(cpack,&itemp,iofst,nbitsd-1);
307 iofst=iofst+nbitsd-1;
308 }
309 }
310 //print *,'SDp ',ival1,ival2,minsd,nbitsd
311 } // end of spatial diff section
312 //
313 // Expand non-missing data values to original grid.
314 //
315 miss1=jfld[0];
316 for ( j=0; j<nonmiss; j++) if (jfld[j] < miss1) miss1 = jfld[j];
317 if( miss1 <= INT_MIN + 1 )
318 {
319 // E. Rouault: no idea if this is correct, but avoids integer
320 // wrap over
321 miss1++;
322 miss2 = miss1 + 1;
323 }
324 else
325 {
326 miss1--;
327 miss2=miss1-1;
328 }
329 n=0;
330 for ( j=0; j<ndpts; j++) {
331 if ( ifldmiss[j] == 0 ) {
332 ifld[j]=jfld[n];
333 n++;
334 }
335 else if ( ifldmiss[j] == 1 ) {
336 ifld[j]=miss1;
337 }
338 else if ( ifldmiss[j] == 2 ) {
339 ifld[j]=miss2;
340 }
341 }
342 //
343 // Determine Groups to be used.
344 //
345 if ( simple_alg == 1 ) {
346 // set group length to 10 : calculate number of groups
347 // and length of last group
348 ngroups=ndpts/10;
349 for (j=0;j<ngroups;j++) glen[j]=10;
350 itemp=ndpts%10;
351 if (itemp != 0) {
352 ngroups++;
353 glen[ngroups-1]=itemp;
354 }
355 }
356 else {
357 // Use Dr. Glahn's algorithm for determining grouping.
358 //
359 kfildo=6;
360 minpk=10;
361 inc=1;
362 maxgrps=(ndpts/minpk)+1;
363 jmin = calloc(maxgrps,sizeof(g2int));
364 jmax = calloc(maxgrps,sizeof(g2int));
365 lbit = calloc(maxgrps,sizeof(g2int));
366 pack_gp(&kfildo,ifld,&ndpts,&missopt,&minpk,&inc,&miss1,&miss2,
367 jmin,jmax,lbit,glen,&maxgrps,&ngroups,&ibit,&jbit,
368 &kbit,&novref,&lbitref,&ier);
369 //printf("SAGier = %d %d %d %d %d %d\n",ier,ibit,jbit,kbit,novref,lbitref);
370 for ( ng=0; ng<ngroups; ng++) glen[ng]=glen[ng]+novref;
371 free(jmin);
372 free(jmax);
373 free(lbit);
374 if( ier != 0 )
375 {
376 free(ifld);
377 free(jfld);
378 free(ifldmiss);
379 free(gref);
380 free(gwidth);
381 free(glen);
382 *lcpack = -1;
383 return;
384 }
385 }
386 //
387 // For each group, find the group's reference value (min)
388 // and the number of bits needed to hold the remaining values
389 //
390 n=0;
391 for ( ng=0; ng<ngroups; ng++) {
392 // how many of each type?
393 num0=num1=num2=0;
394 for (j=n; j<n+glen[ng]; j++) {
395 if (ifldmiss[j] == 0 ) num0++;
396 if (ifldmiss[j] == 1 ) num1++;
397 if (ifldmiss[j] == 2 ) num2++;
398 }
399 if ( num0 == 0 ) { // all missing values
400 if ( num1 == 0 ) { // all secondary missing
401 gref[ng]=-2;
402 gwidth[ng]=0;
403 }
404 else if ( num2 == 0 ) { // all primary missing
405 gref[ng]=-1;
406 gwidth[ng]=0;
407 }
408 else { // both primary and secondary
409 gref[ng]=0;
410 gwidth[ng]=1;
411 }
412 }
413 else { // contains some non-missing data
414 // find max and min values of group
415 gref[ng]=2147483647;
416 imax=-2147483647;
417 j=n;
418 for ( lg=0; lg<glen[ng]; lg++ ) {
419 if ( ifldmiss[j] == 0 ) {
420 if (ifld[j] < gref[ng]) gref[ng]=ifld[j];
421 if (ifld[j] > imax) imax=ifld[j];
422 }
423 j++;
424 }
425 if (missopt == 1) imax=imax+1;
426 if (missopt == 2) imax=imax+2;
427 // calc num of bits needed to hold data
428 if ( gref[ng] != imax ) {
429 temp=(float)(log((double)(imax-gref[ng]+1))/alog2);
430 gwidth[ng]=(g2int)ceil(temp);
431 }
432 else {
433 gwidth[ng]=0;
434 }
435 }
436 // Subtract min from data
437 j=n;
438 mtemp=(g2int)int_power(2.,gwidth[ng]);
439 for ( lg=0; lg<glen[ng]; lg++ ) {
440 if (ifldmiss[j] == 0) // non-missing
441 ifld[j]=ifld[j]-gref[ng];
442 else if (ifldmiss[j] == 1) // primary missing
443 ifld[j]=mtemp-1;
444 else if (ifldmiss[j] == 2) // secondary missing
445 ifld[j]=mtemp-2;
446
447 j++;
448 }
449 // increment fld array counter
450 n=n+glen[ng];
451 }
452 //
453 // Find max of the group references and calc num of bits needed
454 // to pack each groups reference value, then
455 // pack up group reference values
456 //
457 //printf(" GREFS: ");
458 //for (j=0;j<ngroups;j++) printf(" %d",gref[j]); printf("\n");
459 igmax=gref[0];
460 for (j=1;j<ngroups;j++) if (gref[j] > igmax) igmax=gref[j];
461 if (missopt == 1) igmax=igmax+1;
462 if (missopt == 2) igmax=igmax+2;
463 if (igmax != 0) {
464 temp=(float)(log((double)(igmax+1))/alog2);
465 nbitsgref=(g2int)ceil(temp);
466 if( nbitsgref < 0 || nbitsgref >= 31 )
467 {
468 free(ifld);
469 free(jfld);
470 free(ifldmiss);
471 free(gref);
472 free(gwidth);
473 free(glen);
474 *lcpack = -1;
475 return;
476 }
477 // reset the ref values of any "missing only" groups.
478 mtemp=(g2int)int_power(2.,nbitsgref);
479 for ( j=0; j<ngroups; j++ ) {
480 if (gref[j] == -1) gref[j]=mtemp-1;
481 if (gref[j] == -2) gref[j]=mtemp-2;
482 }
483 sbits(cpack,gref,iofst,nbitsgref,0,ngroups);
484 itemp=nbitsgref*ngroups;
485 iofst=iofst+itemp;
486 // Pad last octet with Zeros, if necessary,
487 if ( (itemp%8) != 0) {
488 left=8-(itemp%8);
489 sbit(cpack,&zero,iofst,left);
490 iofst=iofst+left;
491 }
492 }
493 else {
494 nbitsgref=0;
495 }
496 //
497 // Find max/min of the group widths and calc num of bits needed
498 // to pack each groups width value, then
499 // pack up group width values
500 //
501 //write(77,*)'GWIDTHS: ',(gwidth(j),j=1,ngroups)
502 iwmax=gwidth[0];
503 ngwidthref=gwidth[0];
504 for (j=1;j<ngroups;j++) {
505 if (gwidth[j] > iwmax) iwmax=gwidth[j];
506 if (gwidth[j] < ngwidthref) ngwidthref=gwidth[j];
507 }
508 if (iwmax != ngwidthref) {
509 temp=(float)(log((double)(iwmax-ngwidthref+1))/alog2);
510 nbitsgwidth=(g2int)ceil(temp);
511 for ( i=0; i<ngroups; i++) gwidth[i]=gwidth[i]-ngwidthref;
512 sbits(cpack,gwidth,iofst,nbitsgwidth,0,ngroups);
513 itemp=nbitsgwidth*ngroups;
514 iofst=iofst+itemp;
515 // Pad last octet with Zeros, if necessary,
516 if ( (itemp%8) != 0) {
517 left=8-(itemp%8);
518 sbit(cpack,&zero,iofst,left);
519 iofst=iofst+left;
520 }
521 }
522 else {
523 nbitsgwidth=0;
524 for (i=0;i<ngroups;i++) gwidth[i]=0;
525 }
526 //
527 // Find max/min of the group lengths and calc num of bits needed
528 // to pack each groups length value, then
529 // pack up group length values
530 //
531 //printf(" GLENS: ");
532 //for (j=0;j<ngroups;j++) printf(" %d",glen[j]); printf("\n");
533 ilmax=glen[0];
534 nglenref=glen[0];
535 for (j=1;j<ngroups-1;j++) {
536 if (glen[j] > ilmax) ilmax=glen[j];
537 if (glen[j] < nglenref) nglenref=glen[j];
538 }
539 nglenlast=glen[ngroups-1];
540 if (ilmax != nglenref) {
541 temp=(float)(log((double)(ilmax-nglenref+1))/alog2);
542 nbitsglen=(g2int)ceil(temp);
543 for ( i=0; i<ngroups-1; i++) glen[i]=glen[i]-nglenref;
544 sbits(cpack,glen,iofst,nbitsglen,0,ngroups);
545 itemp=nbitsglen*ngroups;
546 iofst=iofst+itemp;
547 // Pad last octet with Zeros, if necessary,
548 if ( (itemp%8) != 0) {
549 left=8-(itemp%8);
550 sbit(cpack,&zero,iofst,left);
551 iofst=iofst+left;
552 }
553 }
554 else {
555 nbitsglen=0;
556 for (i=0;i<ngroups;i++) glen[i]=0;
557 }
558 //
559 // For each group, pack data values
560 //
561 //write(77,*)'IFLDS: ',(ifld(j),j=1,ndpts)
562 n=0;
563 // ij=0;
564 for ( ng=0; ng<ngroups; ng++) {
565 glength=glen[ng]+nglenref;
566 if (ng == (ngroups-1) ) glength=nglenlast;
567 grpwidth=gwidth[ng]+ngwidthref;
568 //write(77,*)'NGP ',ng,grpwidth,glength,gref(ng)
569 if ( grpwidth != 0 ) {
570 sbits(cpack,ifld+n,iofst,grpwidth,0,glength);
571 iofst=iofst+(grpwidth*glength);
572 }
573 // do kk=1,glength
574 // ij=ij+1
575 //write(77,*)'SAG ',ij,fld(ij),ifld(ij),gref(ng),bscale,rmin,dscale
576 // enddo
577 n=n+glength;
578 }
579 // Pad last octet with Zeros, if necessary,
580 if ( (iofst%8) != 0) {
581 left=8-(iofst%8);
582 sbit(cpack,&zero,iofst,left);
583 iofst=iofst+left;
584 }
585 *lcpack=iofst/8;
586 //
587 free(ifld);
588 free(jfld);
589 free(ifldmiss);
590 free(gref);
591 free(gwidth);
592 free(glen);
593 //}
594 //else { // Constant field ( max = min )
595 // nbits=0;
596 // *lcpack=0;
597 // nbitsgref=0;
598 // ngroups=0;
599 //}
600
601 //
602 // Fill in ref value and number of bits in Template 5.2
603 //
604 mkieee(&rmin,idrstmpl+0,1); // ensure reference value is IEEE format
605 idrstmpl[3]=nbitsgref;
606 idrstmpl[4]=0; // original data were reals
607 idrstmpl[5]=1; // general group splitting
608 idrstmpl[9]=ngroups; // Number of groups
609 idrstmpl[10]=ngwidthref; // reference for group widths
610 idrstmpl[11]=nbitsgwidth; // num bits used for group widths
611 idrstmpl[12]=nglenref; // Reference for group lengths
612 idrstmpl[13]=1; // length increment for group lengths
613 idrstmpl[14]=nglenlast; // True length of last group
614 idrstmpl[15]=nbitsglen; // num bits used for group lengths
615 if (idrsnum == 3) {
616 idrstmpl[17]=nbitsd/8; // num bits used for extra spatial
617 // differencing values
618 }
619
620 }
621