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