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 68 : static float DoubleToFloatClamp(double val) {
10 68 : if (val >= FLT_MAX) return FLT_MAX;
11 68 : if (val <= -FLT_MAX) return -FLT_MAX;
12 68 : return (float)val;
13 : }
14 : #endif
15 :
16 34 : 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 34 : g2int nbitsd=0,isign;
60 34 : g2int j,iofst,ival1,ival2,minsd,itemp,l,k,n,non=0;
61 34 : g2int *ifld=NULL,*ifldmiss=NULL;
62 34 : 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 34 : g2int errorOccurred = 0;
68 :
69 : //printf('IDRSTMPL: ',(idrstmpl(j),j=1,16)
70 34 : rdieee(idrstmpl+0,&ref,1);
71 : // printf("SAGTref: %f\n",ref);
72 34 : bscale = DoubleToFloatClamp(int_power(2.0,idrstmpl[1]));
73 34 : dscale = DoubleToFloatClamp(int_power(10.0,-idrstmpl[2]));
74 34 : nbitsgref = idrstmpl[3];
75 34 : itype = idrstmpl[4];
76 34 : ngroups = idrstmpl[9];
77 34 : nbitsgwidth = idrstmpl[11];
78 34 : nbitsglen = idrstmpl[15];
79 34 : if (idrsnum == 3)
80 16 : nbitsd=idrstmpl[17]*8;
81 :
82 : // Constant field
83 :
84 34 : 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 26 : if( ngroups < 0 || ngroups - 10 > ndpts / 2 ) {
94 0 : return -1;
95 : }
96 :
97 : /* Early test in particular case for more general test belows */
98 : /* "Test to see if the group widths and lengths are consistent with number of */
99 : /* values, and length of section 7. */
100 26 : if( idrstmpl[12] < 0 || idrstmpl[14] < 0 || idrstmpl[14] > ndpts )
101 0 : return -1;
102 26 : if( nbitsglen == 0 &&
103 0 : ((ngroups > 1 && idrstmpl[12] != (ndpts - idrstmpl[14]) / (ngroups - 1)) ||
104 0 : idrstmpl[12] * (ngroups-1) + idrstmpl[14] != ndpts) )
105 : {
106 0 : return -1;
107 : }
108 :
109 26 : iofst=0;
110 26 : ifld=(g2int *)calloc(ndpts,sizeof(g2int));
111 : //printf("ALLOC ifld: %d %x\n",(int)ndpts,ifld);
112 26 : gref=(g2int *)calloc(ngroups,sizeof(g2int));
113 : //printf("ALLOC gref: %d %x\n",(int)ngroups,gref);
114 26 : gwidth=(g2int *)calloc(ngroups,sizeof(g2int));
115 : //printf("ALLOC gwidth: %d %x\n",(int)ngroups,gwidth);
116 26 : 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 26 : 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 26 : 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 26 : if (idrsnum == 3) {
148 16 : if (nbitsd != 0) {
149 : // one mistake here should be unsigned int
150 16 : gbit(cpack,&ival1,iofst,nbitsd);
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 13 : gbit(cpack,&ival2,iofst,nbitsd);
160 13 : 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 : gbit(cpack,&isign,iofst,1);
168 16 : iofst=iofst+1;
169 16 : gbit(cpack,&minsd,iofst,nbitsd-1);
170 16 : iofst=iofst+nbitsd-1;
171 16 : if (isign == 1 && minsd != INT_MIN) minsd=-minsd;
172 : }
173 : else {
174 0 : ival1=0;
175 0 : ival2=0;
176 0 : minsd=0;
177 : }
178 : //printf("SDu %ld %ld %ld %ld \n",ival1,ival2,minsd,nbitsd);
179 : }
180 : //
181 : // Extract Each Group's reference value
182 : //
183 : //printf("SAG1: %ld %ld %ld \n",nbitsgref,ngroups,iofst);
184 26 : if (nbitsgref != 0) {
185 26 : if( gbits(cpack,cpack_length,gref+0,iofst,nbitsgref,0,ngroups) != 0 )
186 : {
187 0 : free(ifld);
188 0 : free(gwidth);
189 0 : free(gref);
190 0 : return -1;
191 : }
192 26 : itemp=nbitsgref*ngroups;
193 26 : iofst=iofst+itemp;
194 26 : if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
195 : }
196 : #if 0
197 : else {
198 : for (j=0;j<ngroups;j++)
199 : gref[j]=0;
200 : }
201 : #endif
202 :
203 : //
204 : // Extract Each Group's bit width
205 : //
206 : //printf("SAG2: %ld %ld %ld %ld \n",nbitsgwidth,ngroups,iofst,idrstmpl[10]);
207 26 : if (nbitsgwidth != 0) {
208 26 : if( gbits(cpack,cpack_length,gwidth+0,iofst,nbitsgwidth,0,ngroups) != 0 )
209 : {
210 0 : free(ifld);
211 0 : free(gwidth);
212 0 : free(gref);
213 0 : return -1;
214 : }
215 26 : itemp=nbitsgwidth*ngroups;
216 26 : iofst=iofst+itemp;
217 26 : if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
218 : }
219 : #if 0
220 : else {
221 : for (j=0;j<ngroups;j++)
222 : gwidth[j]=0;
223 : }
224 : #endif
225 :
226 105696 : for (j=0;j<ngroups;j++)
227 : {
228 105670 : if( gwidth[j] > INT_MAX - idrstmpl[10] )
229 : {
230 0 : free(ifld);
231 0 : free(gwidth);
232 0 : free(gref);
233 0 : return -1;
234 : }
235 105670 : gwidth[j]=gwidth[j]+idrstmpl[10];
236 : }
237 :
238 : //
239 : // Extract Each Group's length (number of values in each group)
240 : //
241 26 : glen=(g2int *)calloc(ngroups,sizeof(g2int));
242 26 : if( glen == NULL )
243 : {
244 0 : free(ifld);
245 0 : free(gwidth);
246 0 : free(gref);
247 0 : return -1;
248 : }
249 : //printf("ALLOC glen: %d %x\n",(int)ngroups,glen);
250 : //printf("SAG3: %ld %ld %ld %ld %ld \n",nbitsglen,ngroups,iofst,idrstmpl[13],idrstmpl[12]);
251 26 : if (nbitsglen != 0) {
252 26 : if( gbits(cpack,cpack_length,glen,iofst,nbitsglen,0,ngroups) != 0 )
253 : {
254 0 : free(ifld);
255 0 : free(gwidth);
256 0 : free(glen);
257 0 : free(gref);
258 0 : return -1;
259 : }
260 26 : itemp=nbitsglen*ngroups;
261 26 : iofst=iofst+itemp;
262 26 : if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
263 : }
264 : #if 0
265 : else {
266 : for (j=0;j<ngroups;j++)
267 : glen[j]=0;
268 : }
269 : #endif
270 :
271 : // Note: iterate only to ngroups - 1 since we override just after the
272 : // loop. Otherwise the security checks in the loop might trigger, like
273 : // on band 23 of
274 : // https://nomads.ncep.noaa.gov/pub/data/nccf/com/blend/prod/blend.20200605/16/grib2/blend.t16z.master.f001.co.grib2
275 105670 : for (j=0;j<ngroups - 1;j++)
276 : {
277 105644 : if( glen[j] < 0 ||
278 105644 : (idrstmpl[13] != 0 && glen[j] > INT_MAX / idrstmpl[13]) ||
279 105644 : glen[j] * idrstmpl[13] > INT_MAX - idrstmpl[12] )
280 : {
281 0 : free(ifld);
282 0 : free(gwidth);
283 0 : free(glen);
284 0 : free(gref);
285 0 : return -1;
286 : }
287 105644 : glen[j]=(glen[j]*idrstmpl[13])+idrstmpl[12];
288 : }
289 26 : glen[ngroups-1]=idrstmpl[14];
290 : //
291 : // Test to see if the group widths and lengths are consistent with number of
292 : // values, and length of section 7.
293 : //
294 26 : totBit = 0;
295 26 : totLen = 0;
296 105696 : for (j=0;j<ngroups;j++) {
297 : g2int width_mult_len;
298 105670 : if( gwidth[j] < 0 || glen[j] < 0 ||
299 105670 : (gwidth[j] > 0 && glen[j] > INT_MAX / gwidth[j]) )
300 : {
301 0 : errorOccurred = 1;
302 0 : break;
303 : }
304 105670 : width_mult_len = gwidth[j]*glen[j];
305 105670 : if( totBit > INT_MAX - width_mult_len )
306 : {
307 0 : errorOccurred = 1;
308 0 : break;
309 : }
310 105670 : totBit += width_mult_len;
311 105670 : if( totLen > INT_MAX - glen[j] )
312 : {
313 0 : errorOccurred = 1;
314 0 : break;
315 : }
316 105670 : totLen += glen[j];
317 : }
318 26 : if (errorOccurred || totLen != ndpts || totBit / 8. > lensec) {
319 0 : free(ifld);
320 0 : free(gwidth);
321 0 : free(glen);
322 0 : free(gref);
323 0 : return 1;
324 : }
325 : //
326 : // For each group, unpack data values
327 : //
328 26 : if ( idrstmpl[6] == 0 ) { // no missing values
329 9 : n=0;
330 18401 : for (j=0;j<ngroups;j++) {
331 18392 : if (gwidth[j] != 0) {
332 17275 : if( gbits(cpack,cpack_length,ifld+n,iofst,gwidth[j],0,glen[j]) != 0 )
333 : {
334 0 : free(ifld);
335 0 : free(gwidth);
336 0 : free(glen);
337 0 : free(gref);
338 0 : return -1;
339 : }
340 270272 : for (k=0;k<glen[j];k++) {
341 252997 : ifld[n]=ifld[n]+gref[j];
342 252997 : n=n+1;
343 : }
344 : }
345 : else {
346 2321 : for (l=n;l<n+glen[j];l++) ifld[l]=gref[j];
347 1117 : n=n+glen[j];
348 : }
349 18392 : iofst=iofst+(gwidth[j]*glen[j]);
350 : }
351 : }
352 17 : else if ( idrstmpl[6]==1 || idrstmpl[6]==2 ) {
353 : // missing values included
354 17 : ifldmiss=(g2int *)malloc(ndpts*sizeof(g2int));
355 : //printf("ALLOC ifldmiss: %d %x\n",(int)ndpts,ifldmiss);
356 : //for (j=0;j<ndpts;j++) ifldmiss[j]=0;
357 17 : n=0;
358 17 : non=0;
359 87295 : for (j=0;j<ngroups;j++) {
360 : //printf(" SAGNGP %d %d %d %d\n",j,gwidth[j],glen[j],gref[j]);
361 87278 : if (gwidth[j] != 0) {
362 75982 : msng1=(g2int)int_power(2.0,gwidth[j])-1;
363 75982 : msng2=msng1-1;
364 75982 : if( gbits(cpack,cpack_length,ifld+n,iofst,gwidth[j],0,glen[j]) != 0 )
365 : {
366 0 : free(ifld);
367 0 : free(gwidth);
368 0 : free(glen);
369 0 : free(gref);
370 0 : free(ifldmiss);
371 0 : return -1;
372 : }
373 75982 : iofst=iofst+(gwidth[j]*glen[j]);
374 1157860 : for (k=0;k<glen[j];k++) {
375 1081880 : if (ifld[n] == msng1) {
376 86597 : ifldmiss[n]=1;
377 : //ifld[n]=0;
378 : }
379 995279 : else if (idrstmpl[6]==2 && ifld[n]==msng2) {
380 0 : ifldmiss[n]=2;
381 : //ifld[n]=0;
382 : }
383 : else {
384 995279 : ifldmiss[n]=0;
385 995279 : ifld[non++]=ifld[n]+gref[j];
386 : }
387 1081880 : n++;
388 : }
389 : }
390 : else {
391 11296 : msng1=(g2int)int_power(2.0,nbitsgref)-1;
392 11296 : msng2=msng1-1;
393 11296 : if (gref[j] == msng1) {
394 1303530 : for (l=n;l<n+glen[j];l++) ifldmiss[l]=1;
395 : }
396 0 : else if (idrstmpl[6]==2 && gref[j]==msng2) {
397 0 : for (l=n;l<n+glen[j];l++) ifldmiss[l]=2;
398 : }
399 : else {
400 0 : for (l=n;l<n+glen[j];l++) ifldmiss[l]=0;
401 0 : for (l=non;l<non+glen[j];l++) ifld[l]=gref[j];
402 0 : non += glen[j];
403 : }
404 11296 : n=n+glen[j];
405 : }
406 : }
407 : }
408 :
409 26 : if ( gref != 0 ) free(gref);
410 26 : if ( gwidth != 0 ) free(gwidth);
411 26 : if ( glen != 0 ) free(glen);
412 : //
413 : // If using spatial differences, add overall min value, and
414 : // sum up recursively
415 : //
416 : //printf("SAGod: %ld %ld\n",idrsnum,idrstmpl[16]);
417 26 : if (idrsnum == 3) { // spatial differencing
418 16 : if (idrstmpl[16] == 1) { // first order
419 3 : ifld[0]=ival1;
420 3 : if ( idrstmpl[6] == 0 ) itemp=ndpts; // no missing values
421 1 : else itemp=non;
422 19877 : for (n=1;n<itemp;n++) {
423 19874 : if( (minsd > 0 && ifld[n] > INT_MAX - minsd) ||
424 19874 : (minsd < 0 && ifld[n] < INT_MIN - minsd) )
425 : {
426 0 : free(ifldmiss);
427 0 : free(ifld);
428 0 : return -1;
429 : }
430 19874 : ifld[n]=ifld[n]+minsd;
431 19874 : if( (ifld[n-1] > 0 && ifld[n] > INT_MAX - ifld[n-1]) ||
432 19874 : (ifld[n-1] < 0 && ifld[n] < INT_MIN - ifld[n-1]) )
433 : {
434 0 : free(ifldmiss);
435 0 : free(ifld);
436 0 : return -1;
437 : }
438 19874 : ifld[n]=ifld[n]+ifld[n-1];
439 : }
440 : }
441 13 : else if (idrstmpl[16] == 2) { // second order
442 13 : ifld[0]=ival1;
443 13 : ifld[1]=ival2;
444 13 : if ( idrstmpl[6] == 0 ) itemp=ndpts; // no missing values
445 11 : else itemp=non;
446 461235 : for (n=2;n<itemp;n++) {
447 : /* Lazy way of detecting int overflows: operate on double */
448 461222 : double tmp = (double)ifld[n]+(double)minsd+(2.0 * ifld[n-1])-ifld[n-2];
449 461222 : if( tmp > INT_MAX || tmp < INT_MIN )
450 : {
451 0 : free(ifldmiss);
452 0 : free(ifld);
453 0 : return -1;
454 : }
455 461222 : ifld[n]=(int)tmp;
456 : }
457 : }
458 : }
459 : //
460 : // Scale data back to original form
461 : //
462 : //printf("SAGT: %f %f %f\n",ref,bscale,dscale);
463 26 : if ( idrstmpl[6] == 0 ) { // no missing values
464 254210 : for (n=0;n<ndpts;n++) {
465 254201 : fld[n]=(((g2float)ifld[n]*bscale)+ref)*dscale;
466 : }
467 : }
468 17 : else if ( idrstmpl[6]==1 || idrstmpl[6]==2 ) {
469 : // missing values included
470 17 : non=0;
471 2374130 : for (n=0;n<ndpts;n++) {
472 2374110 : if ( ifldmiss[n] == 0 ) {
473 995279 : fld[n]=(((g2float)ifld[non++]*bscale)+ref)*dscale;
474 : //printf(" SAG %d %f %d %f %f %f\n",n,fld[n],ifld[non-1],bscale,ref,dscale);
475 : }
476 1378830 : else if ( ifldmiss[n] == 1 )
477 1378830 : fld[n]=rmiss1;
478 0 : else if ( ifldmiss[n] == 2 )
479 0 : fld[n]=rmiss2;
480 : }
481 17 : if ( ifldmiss != 0 ) free(ifldmiss);
482 : }
483 :
484 26 : if ( ifld != 0 ) free(ifld);
485 :
486 26 : return(0);
487 :
488 : }
|
