1 #include <float.h>
2 #include <limits.h>
3 #include <stdio.h>
4 #include <stdlib.h>
5 #include "grib2.h"
6
DoubleToFloatClamp(double val)7 static float DoubleToFloatClamp(double val) {
8 if (val >= FLT_MAX) return FLT_MAX;
9 if (val <= -FLT_MAX) return -FLT_MAX;
10 return (float)val;
11 }
12
comunpack(unsigned char * cpack,g2int cpack_length,g2int lensec,g2int idrsnum,g2int * idrstmpl,g2int ndpts,g2float * fld)13 int comunpack(unsigned char *cpack,g2int cpack_length,g2int lensec,g2int idrsnum,g2int *idrstmpl,g2int ndpts,g2float *fld)
14 ////$$$ SUBPROGRAM DOCUMENTATION BLOCK
15 // . . . .
16 // SUBPROGRAM: comunpack
17 // PRGMMR: Gilbert ORG: W/NP11 DATE: 2002-10-29
18 //
19 // ABSTRACT: This subroutine unpacks a data field that was packed using a
20 // complex packing algorithm as defined in the GRIB2 documentation,
21 // using info from the GRIB2 Data Representation Template 5.2 or 5.3.
22 // Supports GRIB2 complex packing templates with or without
23 // spatial differences (i.e. DRTs 5.2 and 5.3).
24 //
25 // PROGRAM HISTORY LOG:
26 // 2002-10-29 Gilbert
27 // 2004-12-16 Gilbert - Added test ( provided by Arthur Taylor/MDL )
28 // to verify that group widths and lengths are
29 // consistent with section length.
30 //
31 // USAGE: int comunpack(unsigned char *cpack,g2int lensec,g2int idrsnum,
32 // g2int *idrstmpl, g2int ndpts,g2float *fld)
33 // INPUT ARGUMENT LIST:
34 // cpack - pointer to the packed data field.
35 // lensec - length of section 7 (used for error checking).
36 // idrsnum - Data Representation Template number 5.N
37 // Must equal 2 or 3.
38 // idrstmpl - pointer to the array of values for Data Representation
39 // Template 5.2 or 5.3
40 // ndpts - The number of data values to unpack
41 //
42 // OUTPUT ARGUMENT LIST:
43 // fld - Contains the unpacked data values. fld must be allocated
44 // with at least ndpts*sizeof(g2float) bytes before
45 // calling this routine.
46 //
47 // REMARKS: None
48 //
49 // ATTRIBUTES:
50 // LANGUAGE: C
51 // MACHINE:
52 //
53 //$$$//
54 {
55
56 g2int nbitsd=0,isign;
57 g2int j,iofst,ival1,ival2,minsd,itemp,l,k,n,non=0;
58 g2int *ifld=NULL,*ifldmiss=NULL;
59 g2int *gref=NULL,*gwidth=NULL,*glen=NULL;
60 g2int itype,ngroups,nbitsgref,nbitsgwidth,nbitsglen;
61 g2int msng1,msng2;
62 g2float ref,bscale,dscale,rmiss1,rmiss2;
63 g2int totBit, totLen;
64 g2int errorOccurred = 0;
65
66 //printf('IDRSTMPL: ',(idrstmpl(j),j=1,16)
67 rdieee(idrstmpl+0,&ref,1);
68 // printf("SAGTref: %f\n",ref);
69 bscale = DoubleToFloatClamp(int_power(2.0,idrstmpl[1]));
70 dscale = DoubleToFloatClamp(int_power(10.0,-idrstmpl[2]));
71 nbitsgref = idrstmpl[3];
72 itype = idrstmpl[4];
73 ngroups = idrstmpl[9];
74 nbitsgwidth = idrstmpl[11];
75 nbitsglen = idrstmpl[15];
76 if (idrsnum == 3)
77 nbitsd=idrstmpl[17]*8;
78
79 // Constant field
80
81 if (ngroups == 0) {
82 for (j=0;j<ndpts;j++) fld[j]=ref;
83 return(0);
84 }
85
86 /* To avoid excessive memory allocations. Not completely sure */
87 /* if this test is appropriate (the 10 and 2 are arbitrary), */
88 /* but it doesn't seem to make sense to have ngroups much larger than */
89 /* ndpts */
90 if( ngroups < 0 || ngroups - 10 > ndpts / 2 ) {
91 return -1;
92 }
93
94 /* Early test in particular case for more general test belows */
95 /* "Test to see if the group widths and lengths are consistent with number of */
96 /* values, and length of section 7. */
97 if( idrstmpl[12] < 0 || idrstmpl[14] < 0 || idrstmpl[14] > ndpts )
98 return -1;
99 if( nbitsglen == 0 &&
100 ((ngroups > 1 && idrstmpl[12] != (ndpts - idrstmpl[14]) / (ngroups - 1)) ||
101 idrstmpl[12] * (ngroups-1) + idrstmpl[14] != ndpts) )
102 {
103 return -1;
104 }
105
106 iofst=0;
107 ifld=(g2int *)calloc(ndpts,sizeof(g2int));
108 //printf("ALLOC ifld: %d %x\n",(int)ndpts,ifld);
109 gref=(g2int *)calloc(ngroups,sizeof(g2int));
110 //printf("ALLOC gref: %d %x\n",(int)ngroups,gref);
111 gwidth=(g2int *)calloc(ngroups,sizeof(g2int));
112 //printf("ALLOC gwidth: %d %x\n",(int)ngroups,gwidth);
113 if( ifld == NULL || gref == NULL || gwidth == NULL )
114 {
115 free(ifld);
116 free(gref);
117 free(gwidth);
118 return -1;
119 }
120 //
121 // Get missing values, if supplied
122 //
123 if ( idrstmpl[6] == 1 ) {
124 if (itype == 0)
125 rdieee(idrstmpl+7,&rmiss1,1);
126 else
127 rmiss1=(g2float)idrstmpl[7];
128 }
129 if ( idrstmpl[6] == 2 ) {
130 if (itype == 0) {
131 rdieee(idrstmpl+7,&rmiss1,1);
132 rdieee(idrstmpl+8,&rmiss2,1);
133 }
134 else {
135 rmiss1=(g2float)idrstmpl[7];
136 rmiss2=(g2float)idrstmpl[8];
137 }
138 }
139
140 //printf("RMISSs: %f %f %f \n",rmiss1,rmiss2,ref);
141 //
142 // Extract spatial differencing values, if using DRS Template 5.3
143 //
144 if (idrsnum == 3) {
145 if (nbitsd != 0) {
146 // one mistake here should be unsigned int
147 gbit(cpack,&ival1,iofst,nbitsd);
148 iofst=iofst+nbitsd;
149 // gbit(cpack,&isign,iofst,1);
150 // iofst=iofst+1
151 // gbit(cpack,&ival1,iofst,nbitsd-1);
152 // iofst=iofst+nbitsd-1;
153 // if (isign == 1) ival1=-ival1;
154 if (idrstmpl[16] == 2) {
155 // one mistake here should be unsigned int
156 gbit(cpack,&ival2,iofst,nbitsd);
157 iofst=iofst+nbitsd;
158 // gbit(cpack,&isign,iofst,1);
159 // iofst=iofst+1;
160 // gbit(cpack,&ival2,iofst,nbitsd-1);
161 // iofst=iofst+nbitsd-1;
162 // if (isign == 1) ival2=-ival2;
163 }
164 gbit(cpack,&isign,iofst,1);
165 iofst=iofst+1;
166 gbit(cpack,&minsd,iofst,nbitsd-1);
167 iofst=iofst+nbitsd-1;
168 if (isign == 1 && minsd != INT_MIN) minsd=-minsd;
169 }
170 else {
171 ival1=0;
172 ival2=0;
173 minsd=0;
174 }
175 //printf("SDu %ld %ld %ld %ld \n",ival1,ival2,minsd,nbitsd);
176 }
177 //
178 // Extract Each Group's reference value
179 //
180 //printf("SAG1: %ld %ld %ld \n",nbitsgref,ngroups,iofst);
181 if (nbitsgref != 0) {
182 if( gbits(cpack,cpack_length,gref+0,iofst,nbitsgref,0,ngroups) != 0 )
183 {
184 free(ifld);
185 free(gwidth);
186 free(gref);
187 return -1;
188 }
189 itemp=nbitsgref*ngroups;
190 iofst=iofst+itemp;
191 if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
192 }
193 #if 0
194 else {
195 for (j=0;j<ngroups;j++)
196 gref[j]=0;
197 }
198 #endif
199
200 //
201 // Extract Each Group's bit width
202 //
203 //printf("SAG2: %ld %ld %ld %ld \n",nbitsgwidth,ngroups,iofst,idrstmpl[10]);
204 if (nbitsgwidth != 0) {
205 if( gbits(cpack,cpack_length,gwidth+0,iofst,nbitsgwidth,0,ngroups) != 0 )
206 {
207 free(ifld);
208 free(gwidth);
209 free(gref);
210 return -1;
211 }
212 itemp=nbitsgwidth*ngroups;
213 iofst=iofst+itemp;
214 if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
215 }
216 #if 0
217 else {
218 for (j=0;j<ngroups;j++)
219 gwidth[j]=0;
220 }
221 #endif
222
223 for (j=0;j<ngroups;j++)
224 {
225 if( gwidth[j] > INT_MAX - idrstmpl[10] )
226 {
227 free(ifld);
228 free(gwidth);
229 free(gref);
230 return -1;
231 }
232 gwidth[j]=gwidth[j]+idrstmpl[10];
233 }
234
235 //
236 // Extract Each Group's length (number of values in each group)
237 //
238 glen=(g2int *)calloc(ngroups,sizeof(g2int));
239 if( glen == NULL )
240 {
241 free(ifld);
242 free(gwidth);
243 free(gref);
244 return -1;
245 }
246 //printf("ALLOC glen: %d %x\n",(int)ngroups,glen);
247 //printf("SAG3: %ld %ld %ld %ld %ld \n",nbitsglen,ngroups,iofst,idrstmpl[13],idrstmpl[12]);
248 if (nbitsglen != 0) {
249 if( gbits(cpack,cpack_length,glen,iofst,nbitsglen,0,ngroups) != 0 )
250 {
251 free(ifld);
252 free(gwidth);
253 free(glen);
254 free(gref);
255 return -1;
256 }
257 itemp=nbitsglen*ngroups;
258 iofst=iofst+itemp;
259 if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
260 }
261 #if 0
262 else {
263 for (j=0;j<ngroups;j++)
264 glen[j]=0;
265 }
266 #endif
267
268 // Note: iterate only to ngroups - 1 since we override just after the
269 // loop. Otherwise the security checks in the loop might trigger, like
270 // on band 23 of
271 // https://nomads.ncep.noaa.gov/pub/data/nccf/com/blend/prod/blend.20200605/16/grib2/blend.t16z.master.f001.co.grib2
272 for (j=0;j<ngroups - 1;j++)
273 {
274 if( glen[j] < 0 ||
275 (idrstmpl[13] != 0 && glen[j] > INT_MAX / idrstmpl[13]) ||
276 glen[j] * idrstmpl[13] > INT_MAX - idrstmpl[12] )
277 {
278 free(ifld);
279 free(gwidth);
280 free(glen);
281 free(gref);
282 return -1;
283 }
284 glen[j]=(glen[j]*idrstmpl[13])+idrstmpl[12];
285 }
286 glen[ngroups-1]=idrstmpl[14];
287 //
288 // Test to see if the group widths and lengths are consistent with number of
289 // values, and length of section 7.
290 //
291 totBit = 0;
292 totLen = 0;
293 for (j=0;j<ngroups;j++) {
294 g2int width_mult_len;
295 if( gwidth[j] < 0 || glen[j] < 0 ||
296 (gwidth[j] > 0 && glen[j] > INT_MAX / gwidth[j]) )
297 {
298 errorOccurred = 1;
299 break;
300 }
301 width_mult_len = gwidth[j]*glen[j];
302 if( totBit > INT_MAX - width_mult_len )
303 {
304 errorOccurred = 1;
305 break;
306 }
307 totBit += width_mult_len;
308 if( totLen > INT_MAX - glen[j] )
309 {
310 errorOccurred = 1;
311 break;
312 }
313 totLen += glen[j];
314 }
315 if (errorOccurred || totLen != ndpts || totBit / 8. > lensec) {
316 free(ifld);
317 free(gwidth);
318 free(glen);
319 free(gref);
320 return 1;
321 }
322 //
323 // For each group, unpack data values
324 //
325 if ( idrstmpl[6] == 0 ) { // no missing values
326 n=0;
327 for (j=0;j<ngroups;j++) {
328 if (gwidth[j] != 0) {
329 if( gbits(cpack,cpack_length,ifld+n,iofst,gwidth[j],0,glen[j]) != 0 )
330 {
331 free(ifld);
332 free(gwidth);
333 free(glen);
334 free(gref);
335 return -1;
336 }
337 for (k=0;k<glen[j];k++) {
338 ifld[n]=ifld[n]+gref[j];
339 n=n+1;
340 }
341 }
342 else {
343 for (l=n;l<n+glen[j];l++) ifld[l]=gref[j];
344 n=n+glen[j];
345 }
346 iofst=iofst+(gwidth[j]*glen[j]);
347 }
348 }
349 else if ( idrstmpl[6]==1 || idrstmpl[6]==2 ) {
350 // missing values included
351 ifldmiss=(g2int *)malloc(ndpts*sizeof(g2int));
352 //printf("ALLOC ifldmiss: %d %x\n",(int)ndpts,ifldmiss);
353 //for (j=0;j<ndpts;j++) ifldmiss[j]=0;
354 n=0;
355 non=0;
356 for (j=0;j<ngroups;j++) {
357 //printf(" SAGNGP %d %d %d %d\n",j,gwidth[j],glen[j],gref[j]);
358 if (gwidth[j] != 0) {
359 msng1=(g2int)int_power(2.0,gwidth[j])-1;
360 msng2=msng1-1;
361 if( gbits(cpack,cpack_length,ifld+n,iofst,gwidth[j],0,glen[j]) != 0 )
362 {
363 free(ifld);
364 free(gwidth);
365 free(glen);
366 free(gref);
367 free(ifldmiss);
368 return -1;
369 }
370 iofst=iofst+(gwidth[j]*glen[j]);
371 for (k=0;k<glen[j];k++) {
372 if (ifld[n] == msng1) {
373 ifldmiss[n]=1;
374 //ifld[n]=0;
375 }
376 else if (idrstmpl[6]==2 && ifld[n]==msng2) {
377 ifldmiss[n]=2;
378 //ifld[n]=0;
379 }
380 else {
381 ifldmiss[n]=0;
382 ifld[non++]=ifld[n]+gref[j];
383 }
384 n++;
385 }
386 }
387 else {
388 msng1=(g2int)int_power(2.0,nbitsgref)-1;
389 msng2=msng1-1;
390 if (gref[j] == msng1) {
391 for (l=n;l<n+glen[j];l++) ifldmiss[l]=1;
392 }
393 else if (idrstmpl[6]==2 && gref[j]==msng2) {
394 for (l=n;l<n+glen[j];l++) ifldmiss[l]=2;
395 }
396 else {
397 for (l=n;l<n+glen[j];l++) ifldmiss[l]=0;
398 for (l=non;l<non+glen[j];l++) ifld[l]=gref[j];
399 non += glen[j];
400 }
401 n=n+glen[j];
402 }
403 }
404 }
405
406 if ( gref != 0 ) free(gref);
407 if ( gwidth != 0 ) free(gwidth);
408 if ( glen != 0 ) free(glen);
409 //
410 // If using spatial differences, add overall min value, and
411 // sum up recursively
412 //
413 //printf("SAGod: %ld %ld\n",idrsnum,idrstmpl[16]);
414 if (idrsnum == 3) { // spatial differencing
415 if (idrstmpl[16] == 1) { // first order
416 ifld[0]=ival1;
417 if ( idrstmpl[6] == 0 ) itemp=ndpts; // no missing values
418 else itemp=non;
419 for (n=1;n<itemp;n++) {
420 if( (minsd > 0 && ifld[n] > INT_MAX - minsd) ||
421 (minsd < 0 && ifld[n] < INT_MIN - minsd) )
422 {
423 free(ifldmiss);
424 free(ifld);
425 return -1;
426 }
427 ifld[n]=ifld[n]+minsd;
428 if( (ifld[n-1] > 0 && ifld[n] > INT_MAX - ifld[n-1]) ||
429 (ifld[n-1] < 0 && ifld[n] < INT_MIN - ifld[n-1]) )
430 {
431 free(ifldmiss);
432 free(ifld);
433 return -1;
434 }
435 ifld[n]=ifld[n]+ifld[n-1];
436 }
437 }
438 else if (idrstmpl[16] == 2) { // second order
439 ifld[0]=ival1;
440 ifld[1]=ival2;
441 if ( idrstmpl[6] == 0 ) itemp=ndpts; // no missing values
442 else itemp=non;
443 for (n=2;n<itemp;n++) {
444 /* Lazy way of detecting int overflows: operate on double */
445 double tmp = (double)ifld[n]+(double)minsd+(2.0 * ifld[n-1])-ifld[n-2];
446 if( tmp > INT_MAX || tmp < INT_MIN )
447 {
448 free(ifldmiss);
449 free(ifld);
450 return -1;
451 }
452 ifld[n]=(int)tmp;
453 }
454 }
455 }
456 //
457 // Scale data back to original form
458 //
459 //printf("SAGT: %f %f %f\n",ref,bscale,dscale);
460 if ( idrstmpl[6] == 0 ) { // no missing values
461 for (n=0;n<ndpts;n++) {
462 fld[n]=(((g2float)ifld[n]*bscale)+ref)*dscale;
463 }
464 }
465 else if ( idrstmpl[6]==1 || idrstmpl[6]==2 ) {
466 // missing values included
467 non=0;
468 for (n=0;n<ndpts;n++) {
469 if ( ifldmiss[n] == 0 ) {
470 fld[n]=(((g2float)ifld[non++]*bscale)+ref)*dscale;
471 //printf(" SAG %d %f %d %f %f %f\n",n,fld[n],ifld[non-1],bscale,ref,dscale);
472 }
473 else if ( ifldmiss[n] == 1 )
474 fld[n]=rmiss1;
475 else if ( ifldmiss[n] == 2 )
476 fld[n]=rmiss2;
477 }
478 if ( ifldmiss != 0 ) free(ifldmiss);
479 }
480
481 if ( ifld != 0 ) free(ifld);
482
483 return(0);
484
485 }
486