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