Line data Source code
1 : /*****************************************************************************
2 : * degrib1.c
3 : *
4 : * DESCRIPTION
5 : * This file contains the main driver routines to unpack GRIB 1 files.
6 : *
7 : * HISTORY
8 : * 4/2003 Arthur Taylor (MDL / RSIS): Created.
9 : *
10 : * NOTES
11 : * GRIB 1 files are assumed to be big endian.
12 : *****************************************************************************
13 : */
14 :
15 : #include <assert.h>
16 : #include <stdio.h>
17 : #include <string.h>
18 : #include <stdlib.h>
19 : #include <math.h>
20 :
21 : #include <limits>
22 :
23 : #include "degrib2.h"
24 : #include "myerror.h"
25 : #include "myassert.h"
26 : #include "tendian.h"
27 : #include "scan.h"
28 : #include "degrib1.h"
29 : #include "metaname.h"
30 : #include "clock.h"
31 : #include "cpl_error.h"
32 : #include "cpl_string.h"
33 :
34 : /* default missing data value (see: bitmap GRIB1: sect3 and sect4) */
35 : /* UNDEFINED is default, UNDEFINED_PRIM is desired choice. */
36 : #define UNDEFINED 9.999e20
37 : #define UNDEFINED_PRIM 9999
38 :
39 : #define GRIB_UNSIGN_INT3(a,b,c) ((a<<16)+(b<<8)+c)
40 : #define GRIB_UNSIGN_INT2(a,b) ((a<<8)+b)
41 : #define GRIB_SIGN_INT3(a,b,c) ((1-(int) ((unsigned) (a & 0x80) >> 6)) * (int) (((a & 127) << 16)+(b<<8)+c))
42 : #define GRIB_SIGN_INT2(a,b) ((1-(int) ((unsigned) (a & 0x80) >> 6)) * (int) (((a & 0x7f) << 8) + b))
43 :
44 : /* various centers */
45 : #define NMC 7
46 : #define US_OTHER 9
47 : #define CPTEC 46
48 : /* Canada Center */
49 : #define CMC 54
50 : #define AFWA 57
51 : #define DWD 78
52 : #define ECMWF 98
53 : #define ATHENS 96
54 : #define NORWAY 88
55 :
56 : /* various subcenters */
57 : #define SUBCENTER_MDL 14
58 : #define SUBCENTER_TDL 11
59 :
60 : /* The idea of rean or opn is to give a warning about default choice of
61 : which table to use. */
62 : #define DEF_NCEP_TABLE rean_nowarn
63 : enum Def_NCEP_Table { rean, opn, rean_nowarn, opn_nowarn };
64 :
65 : #if 0 // moved by GDAL in degrib1.h
66 :
67 : /* For an update to these tables see:
68 : * http://www.nco.ncep.noaa.gov/pmb/docs/on388/table2.html
69 : */
70 :
71 : extern GRIB1ParmTable parm_table_ncep_opn[256];
72 : extern GRIB1ParmTable parm_table_ncep_reanal[256];
73 : extern GRIB1ParmTable parm_table_ncep_tdl[256];
74 : extern GRIB1ParmTable parm_table_ncep_mdl[256];
75 : extern GRIB1ParmTable parm_table_omb[256];
76 : extern GRIB1ParmTable parm_table_nceptab_129[256];
77 : extern GRIB1ParmTable parm_table_nceptab_130[256];
78 : extern GRIB1ParmTable parm_table_nceptab_131[256];
79 : extern GRIB1ParmTable parm_table_nceptab_133[256];
80 : extern GRIB1ParmTable parm_table_nceptab_140[256];
81 : extern GRIB1ParmTable parm_table_nceptab_141[256];
82 :
83 : extern GRIB1ParmTable parm_table_nohrsc[256];
84 :
85 : extern GRIB1ParmTable parm_table_cptec_254[256];
86 :
87 : extern GRIB1ParmTable parm_table_afwa_000[256];
88 : extern GRIB1ParmTable parm_table_afwa_001[256];
89 : extern GRIB1ParmTable parm_table_afwa_002[256];
90 : extern GRIB1ParmTable parm_table_afwa_003[256];
91 : extern GRIB1ParmTable parm_table_afwa_010[256];
92 : extern GRIB1ParmTable parm_table_afwa_011[256];
93 :
94 : extern GRIB1ParmTable parm_table_dwd_002[256];
95 : extern GRIB1ParmTable parm_table_dwd_201[256];
96 : extern GRIB1ParmTable parm_table_dwd_202[256];
97 : extern GRIB1ParmTable parm_table_dwd_203[256];
98 :
99 : extern GRIB1ParmTable parm_table_ecmwf_128[256];
100 : extern GRIB1ParmTable parm_table_ecmwf_129[256];
101 : extern GRIB1ParmTable parm_table_ecmwf_130[256];
102 : extern GRIB1ParmTable parm_table_ecmwf_131[256];
103 : extern GRIB1ParmTable parm_table_ecmwf_140[256];
104 : extern GRIB1ParmTable parm_table_ecmwf_150[256];
105 : extern GRIB1ParmTable parm_table_ecmwf_160[256];
106 : extern GRIB1ParmTable parm_table_ecmwf_170[256];
107 : extern GRIB1ParmTable parm_table_ecmwf_180[256];
108 :
109 : extern GRIB1ParmTable parm_table_athens[256];
110 :
111 : extern GRIB1ParmTable parm_table_norway128[256];
112 :
113 : extern GRIB1ParmTable parm_table_cmc[256];
114 :
115 : extern GRIB1ParmTable parm_table_undefined[256];
116 :
117 : #endif
118 :
119 : /*****************************************************************************
120 : * Choose_ParmTable() --
121 : *
122 : * Arthur Taylor / MDL
123 : *
124 : * PURPOSE
125 : * Chooses the correct Parameter table depending on what is in the GRIB1
126 : * message's "Product Definition Section".
127 : *
128 : * ARGUMENTS
129 : * pdsMeta = The filled out pdsMeta data structure to base choice on. (Input)
130 : * center = The Center that created the data (Input)
131 : * subcenter = The Sub Center that created the data (Input)
132 : *
133 : * FILES/DATABASES: None
134 : *
135 : * RETURNS: ParmTable (appropriate parameter table.)
136 : *
137 : * HISTORY
138 : * <unknown> : wgrib library : cnames.c
139 : * 4/2003 Arthur Taylor (MDL/RSIS): Modified
140 : * 10/2005 AAT: Adjusted to take center, subcenter
141 : *
142 : * NOTES
143 : *****************************************************************************
144 : */
145 132 : static const GRIB1ParmTable *Choose_ParmTable (pdsG1Type *pdsMeta,
146 : unsigned short int center,
147 : unsigned short int subcenter)
148 : {
149 : int process; /* The process ID from the GRIB1 message. */
150 :
151 132 : switch (center) {
152 128 : case NMC:
153 128 : if (pdsMeta->mstrVersion <= 3) {
154 109 : switch (subcenter) {
155 0 : case 1:
156 0 : return &parm_table_ncep_reanal[0];
157 0 : case SUBCENTER_TDL:
158 0 : return &parm_table_ncep_tdl[0];
159 0 : case SUBCENTER_MDL:
160 0 : return &parm_table_ncep_mdl[0];
161 : }
162 : }
163 : /* figure out if NCEP opn or reanalysis */
164 128 : switch (pdsMeta->mstrVersion) {
165 74 : case 0:
166 74 : return &parm_table_ncep_opn[0];
167 35 : case 1:
168 : case 2:
169 35 : process = pdsMeta->genProcess;
170 35 : if ((subcenter != 0) || ((process != 80) && (process != 180))) {
171 35 : return &parm_table_ncep_opn[0];
172 : }
173 : #if 0
174 : /* At this point could be either the opn or reanalysis table */
175 : switch (DEF_NCEP_TABLE) {
176 : case opn_nowarn:
177 : return &parm_table_ncep_opn[0];
178 : case rean_nowarn:
179 : return &parm_table_ncep_reanal[0];
180 : }
181 : break;
182 : #else
183 : // ERO: this is the non convoluted version of the above code
184 0 : return &parm_table_ncep_reanal[0];
185 : #endif
186 0 : case 3:
187 0 : return &parm_table_ncep_opn[0];
188 0 : case 128:
189 0 : return &parm_table_omb[0];
190 19 : case 129:
191 19 : return &parm_table_nceptab_129[0];
192 0 : case 130:
193 0 : return &parm_table_nceptab_130[0];
194 0 : case 131:
195 0 : return &parm_table_nceptab_131[0];
196 0 : case 133:
197 0 : return &parm_table_nceptab_133[0];
198 0 : case 140:
199 0 : return &parm_table_nceptab_140[0];
200 0 : case 141:
201 0 : return &parm_table_nceptab_141[0];
202 : }
203 0 : break;
204 0 : case AFWA:
205 0 : switch (subcenter) {
206 0 : case 0:
207 0 : return &parm_table_afwa_000[0];
208 0 : case 1:
209 : case 4:
210 0 : return &parm_table_afwa_001[0];
211 0 : case 2:
212 0 : return &parm_table_afwa_002[0];
213 0 : case 3:
214 0 : return &parm_table_afwa_003[0];
215 0 : case 10:
216 0 : return &parm_table_afwa_010[0];
217 0 : case 11:
218 0 : return &parm_table_afwa_011[0];
219 : /* case 5:*/
220 : /* Didn't have a table 5. */
221 : }
222 0 : break;
223 0 : case ECMWF:
224 0 : switch (pdsMeta->mstrVersion) {
225 0 : case 128:
226 0 : return &parm_table_ecmwf_128[0];
227 0 : case 129:
228 0 : return &parm_table_ecmwf_129[0];
229 0 : case 130:
230 0 : return &parm_table_ecmwf_130[0];
231 0 : case 131:
232 0 : return &parm_table_ecmwf_131[0];
233 0 : case 140:
234 0 : return &parm_table_ecmwf_140[0];
235 0 : case 150:
236 0 : return &parm_table_ecmwf_150[0];
237 0 : case 160:
238 0 : return &parm_table_ecmwf_160[0];
239 0 : case 170:
240 0 : return &parm_table_ecmwf_170[0];
241 0 : case 180:
242 0 : return &parm_table_ecmwf_180[0];
243 0 : case 228:
244 0 : return &parm_table_ecmwf_228[0];
245 : }
246 0 : break;
247 0 : case DWD:
248 0 : switch (pdsMeta->mstrVersion) {
249 0 : case 2:
250 0 : return &parm_table_dwd_002[0];
251 0 : case 201:
252 0 : return &parm_table_dwd_201[0];
253 0 : case 202:
254 0 : return &parm_table_dwd_202[0];
255 0 : case 203:
256 0 : return &parm_table_dwd_203[0];
257 : }
258 0 : break;
259 0 : case CPTEC:
260 0 : switch (pdsMeta->mstrVersion) {
261 0 : case 254:
262 0 : return &parm_table_cptec_254[0];
263 : }
264 0 : break;
265 0 : case US_OTHER:
266 0 : switch (subcenter) {
267 0 : case 163:
268 0 : return &parm_table_nohrsc[0];
269 : /* Based on 11/7/2006 email with Rob Doornbos, mimic'd what wgrib
270 : * did which was to use parm_table_ncep_opn. */
271 0 : case 161:
272 0 : return &parm_table_ncep_opn[0];
273 : }
274 0 : break;
275 0 : case ATHENS:
276 0 : return &parm_table_athens[0];
277 : break;
278 0 : case NORWAY:
279 0 : if (pdsMeta->mstrVersion == 128) {
280 0 : return &parm_table_norway128[0];
281 : }
282 0 : break;
283 0 : case CMC:
284 0 : return &parm_table_cmc[0];
285 : break;
286 : }
287 4 : if (pdsMeta->mstrVersion > 3) {
288 0 : CPLError( CE_Warning, CPLE_AppDefined, "GRIB: Don't understand the parameter table, since center %d-%d used\n"
289 : "parameter table version %d instead of 3 (international exchange).\n"
290 : "Using default for now (which might lead to erroneous interpretation), but please email arthur.taylor@noaa.gov\n"
291 : "about adding this table to his 'degrib1.c' and 'grib1tab.c' files.",
292 0 : center, subcenter, pdsMeta->mstrVersion);
293 : }
294 4 : if (pdsMeta->cat > 127) {
295 2 : CPLError(CE_Warning, CPLE_AppDefined, "GRIB: Parameter %d is > 127, so it falls in the local use section of\n"
296 : "the parameter table (and is undefined on the international table.\n"
297 : "Using default for now(which might lead to erroneous interpretation), but please email arthur.taylor@noaa.gov\n"
298 : "about adding this table to his 'degrib1.c' and 'grib1tab.c' files.",
299 2 : pdsMeta->cat);
300 : }
301 4 : return &parm_table_undefined[0];
302 : }
303 :
304 : /*****************************************************************************
305 : * GRIB1_Table2LookUp() --
306 : *
307 : * Arthur Taylor / MDL
308 : *
309 : * PURPOSE
310 : * Returns the variable name (type of data) and comment (longer form of the
311 : * name) for the data that is in the GRIB1 message.
312 : *
313 : * ARGUMENTS
314 : * name = A pointer to the resulting short name. (Output)
315 : * comment = A pointer to the resulting long name. (Output)
316 : * pdsMeta = The filled out pdsMeta data structure to base choice on. (Input)
317 : * center = The Center that created the data (Input)
318 : * subcenter = The Sub Center that created the data (Input)
319 : *
320 : * FILES/DATABASES: None
321 : *
322 : * RETURNS: void
323 : *
324 : * HISTORY
325 : * 4/2003 Arthur Taylor (MDL/RSIS): Created
326 : * 10/2005 AAT: Adjusted to take center, subcenter
327 : *
328 : * NOTES
329 : *****************************************************************************
330 : */
331 132 : static void GRIB1_Table2LookUp (pdsG1Type *pdsMeta, const char **name,
332 : const char **comment, const char **unit,
333 : int *convert,
334 : unsigned short int center,
335 : unsigned short int subcenter)
336 : {
337 : const GRIB1ParmTable *table; /* The parameter table chosen by the pdsMeta data */
338 :
339 132 : table = Choose_ParmTable (pdsMeta, center, subcenter);
340 132 : if ((center == NMC) && (pdsMeta->mstrVersion == 129)
341 19 : && (pdsMeta->cat == 180)) {
342 0 : if (pdsMeta->timeRange == 3) {
343 0 : *name = "AVGOZCON";
344 0 : *comment = "Average Ozone Concentration";
345 0 : *unit = "PPB";
346 0 : *convert = UC_NONE;
347 0 : return;
348 : }
349 : }
350 132 : *name = table[pdsMeta->cat].name;
351 132 : if( strcmp(*name, CPLSPrintf("var%d", pdsMeta->cat)) == 0 )
352 : {
353 2 : if( center == ECMWF )
354 0 : *name = CPLSPrintf("var%d of table %d of center ECMWF", pdsMeta->cat, pdsMeta->mstrVersion);
355 : else
356 2 : *name = CPLSPrintf("var%d of table %d of center %d", pdsMeta->cat, pdsMeta->mstrVersion, center);
357 : }
358 132 : *comment = table[pdsMeta->cat].comment;
359 132 : *unit = table[pdsMeta->cat].unit;
360 132 : *convert = table[pdsMeta->cat].convert;
361 : /* printf ("%s %s %s\n", *name, *comment, *unit);*/
362 : }
363 :
364 : /* Similar to metaname.c :: ParseLevelName() */
365 132 : static void GRIB1_Table3LookUp (pdsG1Type *pdsMeta, char **shortLevelName,
366 : char **longLevelName)
367 : {
368 132 : uChar type = pdsMeta->levelType;
369 : uChar level1, level2;
370 :
371 132 : free (*shortLevelName);
372 132 : *shortLevelName = nullptr;
373 132 : free (*longLevelName);
374 132 : *longLevelName = nullptr;
375 : /* Find out if val is a 2 part value or not */
376 132 : if (GRIB1Surface[type].f_twoPart) {
377 0 : level1 = (pdsMeta->levelVal >> 8);
378 0 : level2 = (pdsMeta->levelVal & 0xff);
379 0 : reallocSprintf (shortLevelName, "%d-%d-%s", level1, level2,
380 0 : GRIB1Surface[type].name);
381 0 : reallocSprintf (longLevelName, "%d-%d[%s] %s (%s)", level1, level2,
382 0 : GRIB1Surface[type].unit, GRIB1Surface[type].name,
383 0 : GRIB1Surface[type].comment);
384 : } else {
385 132 : reallocSprintf (shortLevelName, "%d-%s", pdsMeta->levelVal,
386 132 : GRIB1Surface[type].name);
387 132 : reallocSprintf (longLevelName, "%d[%s] %s (%s)", pdsMeta->levelVal,
388 132 : GRIB1Surface[type].unit, GRIB1Surface[type].name,
389 132 : GRIB1Surface[type].comment);
390 : }
391 132 : }
392 :
393 : /*****************************************************************************
394 : * fval_360() --
395 : *
396 : * Albion Taylor / ARL
397 : *
398 : * PURPOSE
399 : * Converts an IBM360 floating point number to an IEEE floating point
400 : * number. The IBM floating point spec represents the fraction as the last
401 : * 3 bytes of the number, with 0xffffff being just shy of 1.0. The first byte
402 : * leads with a sign bit, and the last seven bits represent the powers of 16
403 : * (not 2), with a bias of 0x40 giving 16^0.
404 : *
405 : * ARGUMENTS
406 : * aval = A sInt4 containing the original IBM 360 number. (Input)
407 : *
408 : * FILES/DATABASES: None
409 : *
410 : * RETURNS: double = the value that aval represents.
411 : *
412 : * HISTORY
413 : * <unknown> Albion Taylor (ARL): Created
414 : * 4/2003 Arthur Taylor (MDL/RSIS): Cleaned up.
415 : * 5/2003 AAT: some kind of Bug due to optimizations...
416 : * -1055916032 => 0 instead of -1
417 : *
418 : * NOTES
419 : *****************************************************************************
420 : */
421 67 : static double fval_360 (uInt4 aval)
422 : {
423 : short int ptr[4];
424 : #ifdef LITTLE_ENDIAN
425 67 : ptr[3] = ((((aval >> 24) & 0x7f) << 2) + (0x3ff - 0x100)) << 4;
426 67 : ptr[2] = 0;
427 67 : ptr[1] = 0;
428 67 : ptr[0] = 0;
429 : #else
430 : ptr[0] = ((((aval >> 24) & 0x7f) << 2) + (0x3ff - 0x100)) << 4;
431 : ptr[1] = 0;
432 : ptr[2] = 0;
433 : ptr[3] = 0;
434 : #endif
435 : double pow16;
436 67 : memcpy(&pow16, ptr, 8);
437 67 : return ((aval & 0x80000000) ? -pow16 : pow16) *
438 67 : (aval & 0xffffff) / ((double) 0x1000000);
439 : }
440 :
441 : /*****************************************************************************
442 : * ReadGrib1Sect1() --
443 : *
444 : * Arthur Taylor / MDL
445 : *
446 : * PURPOSE
447 : * Parses the GRIB1 "Product Definition Section" or section 1, filling out
448 : * the pdsMeta data structure.
449 : *
450 : * ARGUMENTS
451 : * pds = The compressed part of the message dealing with "PDS". (Input)
452 : * pdsLen = Size of pds in bytes. (Input)
453 : * gribLen = The total length of the GRIB1 message. (Input)
454 : * curLoc = Current location in the GRIB1 message. (Output)
455 : * pdsMeta = The filled out pdsMeta data structure. (Output)
456 : * f_gds = boolean if there is a Grid Definition Section. (Output)
457 : * gridID = The Grid ID. (Output)
458 : * f_bms = boolean if there is a Bitmap Section. (Output)
459 : * DSF = Decimal Scale Factor for unpacking the data. (Output)
460 : * center = The Center that created the data (Output)
461 : * subcenter = The Sub Center that created the data (Output)
462 : *
463 : * FILES/DATABASES: None
464 : *
465 : * RETURNS: int (could use errSprintf())
466 : * 0 = OK
467 : * -1 = gribLen is too small.
468 : *
469 : * HISTORY
470 : * 4/2003 Arthur Taylor (MDL/RSIS): Created.
471 : * 5/2004 AAT: Paid attention to table 5 (Time range indicator) and which of
472 : * P1 and P2 are the valid times.
473 : * 10/2005 AAT: Adjusted to take center, subcenter
474 : *
475 : * NOTES
476 : *****************************************************************************
477 : */
478 132 : static int ReadGrib1Sect1 (uChar *pds, uInt4 pdsLen, uInt4 gribLen, uInt4 *curLoc,
479 : pdsG1Type *pdsMeta, char *f_gds, uChar *gridID,
480 : char *f_bms, short int *DSF,
481 : unsigned short int *center,
482 : unsigned short int *subcenter)
483 : {
484 : uInt4 sectLen; /* Length in bytes of the current section. */
485 : int year; /* The year of the GRIB1 Message. */
486 : double P1_DeltaTime; /* Used to parse the time for P1 */
487 : double P2_DeltaTime; /* Used to parse the time for P2 */
488 : uInt4 uli_temp;
489 : #ifdef DEBUG
490 : /*
491 : int i;
492 : */
493 : #endif
494 : /* We will read the first required 28 bytes */
495 132 : if( pdsLen < 28 )
496 0 : return -1;
497 :
498 132 : sectLen = GRIB_UNSIGN_INT3 (*pds, pds[1], pds[2]);
499 132 : if( sectLen > pdsLen )
500 0 : return -1;
501 : #ifdef DEBUG
502 : /*
503 : printf ("Section 1 length = %ld\n", sectLen);
504 : for (i = 0; i < sectLen; i++) {
505 : printf ("Sect1: item %d = %d\n", i + 1, pds[i]);
506 : }
507 : printf ("Century is item 25\n");
508 : */
509 : #endif
510 132 : *curLoc += sectLen;
511 132 : if (*curLoc > gribLen) {
512 0 : errSprintf ("Ran out of data in PDS (GRIB 1 Section 1)\n");
513 0 : return -1;
514 : }
515 132 : pds += 3;
516 132 : pdsMeta->mstrVersion = *(pds++);
517 132 : *center = *(pds++);
518 132 : pdsMeta->genProcess = *(pds++);
519 132 : *gridID = *(pds++);
520 132 : *f_gds = GRIB2BIT_1 & *pds;
521 132 : *f_bms = GRIB2BIT_2 & *pds;
522 132 : pds++;
523 132 : pdsMeta->cat = *(pds++);
524 132 : pdsMeta->levelType = *(pds++);
525 132 : pdsMeta->levelVal = GRIB_UNSIGN_INT2 (*pds, pds[1]);
526 132 : pds += 2;
527 132 : if (*pds == 0) {
528 : /* The 12 is because we have increased pds by 12. (but 25 is in
529 : * reference of 1..25, so we need another -1) */
530 0 : year = (pds[25 - 13] * 100);
531 : } else {
532 : /* The 12 is because we have increased pds by 12. (but 25 is in
533 : * reference of 1..25, so we need another -1) */
534 132 : year = *pds + ((pds[25 - 13] - 1) * 100);
535 : }
536 :
537 132 : if (ParseTime (&(pdsMeta->refTime), year, pds[1], pds[2], pds[3], pds[4],
538 132 : 0) != 0) {
539 0 : preErrSprintf ("Error In call to ParseTime\n");
540 0 : errSprintf ("(Probably a corrupt file)\n");
541 0 : return -1;
542 : }
543 132 : pds += 5;
544 132 : pdsMeta->timeRange = pds[3];
545 132 : if (ParseSect4Time2secV1 (pds[1], *pds, &P1_DeltaTime) == 0) {
546 132 : pdsMeta->P1 = pdsMeta->refTime + P1_DeltaTime;
547 : } else {
548 0 : pdsMeta->P1 = pdsMeta->refTime;
549 0 : printf ("Warning! : Can't figure out time unit of %u\n", *pds);
550 : }
551 132 : if (ParseSect4Time2secV1 (pds[2], *pds, &P2_DeltaTime) == 0) {
552 132 : pdsMeta->P2 = pdsMeta->refTime + P2_DeltaTime;
553 : } else {
554 0 : pdsMeta->P2 = pdsMeta->refTime;
555 0 : printf ("Warning! : Can't figure out time unit of %u\n", *pds);
556 : }
557 : /* The following is based on Table 5. */
558 : /* Note: For ensemble forecasts, 119 has meaning. */
559 132 : switch (pdsMeta->timeRange) {
560 107 : case 0:
561 : case 1:
562 : case 113:
563 : case 114:
564 : case 115:
565 : case 116:
566 : case 117:
567 : case 118:
568 : case 123:
569 : case 124:
570 107 : pdsMeta->validTime = pdsMeta->P1;
571 107 : break;
572 0 : case 2:
573 : /* Puzzling case. */
574 0 : pdsMeta->validTime = pdsMeta->P2;
575 0 : break;
576 0 : case 3:
577 : case 4:
578 : case 5:
579 : case 51:
580 0 : pdsMeta->validTime = pdsMeta->P2;
581 0 : break;
582 25 : case 10:
583 25 : if (ParseSect4Time2secV1 (GRIB_UNSIGN_INT2 (pds[1], pds[2]), *pds,
584 25 : &P1_DeltaTime) == 0) {
585 25 : pdsMeta->P2 = pdsMeta->P1 = pdsMeta->refTime + P1_DeltaTime;
586 : } else {
587 0 : pdsMeta->P2 = pdsMeta->P1 = pdsMeta->refTime;
588 0 : printf ("Warning! : Can't figure out time unit of %u\n", *pds);
589 : }
590 25 : pdsMeta->validTime = pdsMeta->P1;
591 25 : break;
592 0 : default:
593 0 : pdsMeta->validTime = pdsMeta->P1;
594 : }
595 132 : pds += 4;
596 132 : pdsMeta->Average = GRIB_UNSIGN_INT2 (*pds, pds[1]);
597 132 : pds += 2;
598 132 : pdsMeta->numberMissing = *(pds++);
599 : /* Skip over centry of reference time. */
600 132 : pds++;
601 132 : *subcenter = *(pds++);
602 132 : *DSF = GRIB_SIGN_INT2 (*pds, pds[1]);
603 132 : pds += 2;
604 132 : pdsMeta->f_hasEns = 0;
605 132 : pdsMeta->f_hasProb = 0;
606 132 : pdsMeta->f_hasCluster = 0;
607 132 : if (sectLen < 41) {
608 132 : return 0;
609 : }
610 : /* Following is based on:
611 : * http://www.emc.ncep.noaa.gov/gmb/ens/info/ens_grib.html */
612 0 : if ((*center == NMC) && (*subcenter == 2)) {
613 0 : if (sectLen < 45) {
614 0 : printf ("Warning! Problems with Ensemble section\n");
615 0 : return 0;
616 : }
617 0 : pdsMeta->f_hasEns = 1;
618 0 : pdsMeta->ens.BitFlag = *(pds++);
619 : /* octet21 = pdsMeta->timeRange; = 119 has meaning now */
620 0 : pds += 11;
621 0 : pdsMeta->ens.Application = *(pds++);
622 0 : pdsMeta->ens.Type = *(pds++);
623 0 : pdsMeta->ens.Number = *(pds++);
624 0 : pdsMeta->ens.ProdID = *(pds++);
625 0 : pdsMeta->ens.Smooth = *(pds++);
626 0 : if ((pdsMeta->cat == 191) || (pdsMeta->cat == 192) ||
627 0 : (pdsMeta->cat == 193)) {
628 0 : if (sectLen < 60) {
629 0 : printf ("Warning! Problems with Ensemble Probability section\n");
630 0 : return 0;
631 : }
632 0 : pdsMeta->f_hasProb = 1;
633 0 : pdsMeta->prob.Cat = pdsMeta->cat;
634 0 : pdsMeta->cat = *(pds++);
635 0 : pdsMeta->prob.Type = *(pds++);
636 0 : MEMCPY_BIG (&uli_temp, pds, sizeof (sInt4));
637 0 : pdsMeta->prob.lower = fval_360 (uli_temp);
638 0 : pds += 4;
639 0 : MEMCPY_BIG (&uli_temp, pds, sizeof (sInt4));
640 0 : pdsMeta->prob.upper = fval_360 (uli_temp);
641 0 : pds += 4;
642 0 : pds += 4;
643 : }
644 0 : if ((pdsMeta->ens.Type == 4) || (pdsMeta->ens.Type == 5)) {
645 : /* 87 ... 100 was reserved, but may not be encoded */
646 0 : if ((sectLen < 100) && (sectLen != 86)) {
647 0 : printf ("Warning! Problems with Ensemble Clustering section\n");
648 0 : printf ("Section length == %u\n", sectLen);
649 0 : return 0;
650 : }
651 0 : if (pdsMeta->f_hasProb == 0) {
652 0 : pds += 14;
653 : }
654 0 : pdsMeta->f_hasCluster = 1;
655 0 : pdsMeta->cluster.ensSize = *(pds++);
656 0 : pdsMeta->cluster.clusterSize = *(pds++);
657 0 : pdsMeta->cluster.Num = *(pds++);
658 0 : pdsMeta->cluster.Method = *(pds++);
659 0 : pdsMeta->cluster.NorLat = GRIB_UNSIGN_INT3 (*pds, pds[1], pds[2]);
660 0 : pdsMeta->cluster.NorLat = pdsMeta->cluster.NorLat / 1000.;
661 0 : pds += 3;
662 0 : pdsMeta->cluster.SouLat = GRIB_UNSIGN_INT3 (*pds, pds[1], pds[2]);
663 0 : pdsMeta->cluster.SouLat = pdsMeta->cluster.SouLat / 1000.;
664 0 : pds += 3;
665 0 : pdsMeta->cluster.EasLon = GRIB_UNSIGN_INT3 (*pds, pds[1], pds[2]);
666 0 : pdsMeta->cluster.EasLon = pdsMeta->cluster.EasLon / 1000.;
667 0 : pds += 3;
668 0 : pdsMeta->cluster.WesLon = GRIB_UNSIGN_INT3 (*pds, pds[1], pds[2]);
669 0 : pdsMeta->cluster.WesLon = pdsMeta->cluster.WesLon / 1000.;
670 0 : pds += 3;
671 0 : memcpy (pdsMeta->cluster.Member, pds, 10);
672 0 : pdsMeta->cluster.Member[10] = '\0';
673 : }
674 : /* Following based on:
675 : * http://www.ecmwf.int/publications/manuals/libraries/gribex/
676 : * localGRIBUsage.html */
677 0 : } else if (*center == ECMWF) {
678 0 : if (sectLen < 45) {
679 0 : printf ("Warning! Problems with ECMWF PDS extension\n");
680 0 : return 0;
681 : }
682 : /*
683 : sInt4 i_temp;
684 : pds += 12;
685 : i_temp = GRIB_SIGN_INT2 (pds[3], pds[4]);
686 : printf ("ID %d Class %d Type %d Stream %d", pds[0], pds[1], pds[2],
687 : i_temp);
688 : pds += 5;
689 : printf (" Ver %c%c%c%c, ", pds[0], pds[1], pds[2], pds[3]);
690 : pds += 4;
691 : printf ("Octet-50 %d, Octet-51 %d SectLen %d\n", pds[0], pds[1],
692 : sectLen);
693 : */
694 : } else {
695 0 : printf ("Un-handled possible ensemble section center %u "
696 0 : "subcenter %u\n", *center, *subcenter);
697 : }
698 0 : return 0;
699 : }
700 :
701 : /*****************************************************************************
702 : * Grib1_Inventory() --
703 : *
704 : * Arthur Taylor / MDL
705 : *
706 : * PURPOSE
707 : * Parses the GRIB1 "Product Definition Section" for enough information to
708 : * fill out the inventory data structure so we can do a simple inventory on
709 : * the file in a similar way to how we did it for GRIB2.
710 : *
711 : * ARGUMENTS
712 : * fp = An opened GRIB2 file already at the correct message. (Input)
713 : * gribLen = The total length of the GRIB1 message. (Input)
714 : * inv = The inventory data structure that we need to fill. (Output)
715 : *
716 : * FILES/DATABASES: None
717 : *
718 : * RETURNS: int (could use errSprintf())
719 : * 0 = OK
720 : * -1 = gribLen is too small.
721 : *
722 : * HISTORY
723 : * 4/2003 Arthur Taylor (MDL/RSIS): Created.
724 : *
725 : * NOTES
726 : *****************************************************************************
727 : */
728 66 : int GRIB1_Inventory (VSILFILE *fp, uInt4 gribLen, inventoryType *inv)
729 : {
730 : uChar temp[3]; /* Used to determine the section length. */
731 : uInt4 sectLen; /* Length in bytes of the current section. */
732 : uChar *pds; /* The part of the message dealing with the PDS. */
733 : pdsG1Type pdsMeta; /* The pds parsed into a usable data structure. */
734 : char f_gds; /* flag if there is a gds section. */
735 : char f_bms; /* flag if there is a bms section. */
736 : short int DSF; /* Decimal Scale Factor for unpacking the data. */
737 : uChar gridID; /* Which GDS specs to use. */
738 : const char *varName; /* The name of the data stored in the grid. */
739 : const char *varComment; /* Extra comments about the data stored in grid. */
740 : const char *varUnit; /* Holds the name of the unit [K] [%] .. etc */
741 : int convert; /* Conversion method for this variable's unit. */
742 : uInt4 curLoc; /* Where we are in the current GRIB message. */
743 : unsigned short int center; /* The Center that created the data */
744 : unsigned short int subcenter; /* The Sub Center that created the data */
745 :
746 66 : curLoc = 8;
747 66 : if (VSIFReadL(temp, sizeof (char), 3, fp) != 3) {
748 0 : errSprintf ("Ran out of file.\n");
749 0 : return -1;
750 : }
751 66 : sectLen = GRIB_UNSIGN_INT3 (*temp, temp[1], temp[2]);
752 66 : if (curLoc + sectLen > gribLen) {
753 0 : errSprintf ("Ran out of data in PDS (GRIB1_Inventory)\n");
754 0 : return -1;
755 : }
756 66 : if( sectLen < 3 )
757 : {
758 0 : errSprintf ("Invalid sectLen.\n");
759 0 : return -1;
760 : }
761 66 : pds = (uChar *) malloc (sectLen * sizeof (uChar));
762 66 : if( pds == nullptr )
763 : {
764 0 : errSprintf ("Ran out of memory.\n");
765 0 : return -1;
766 : }
767 66 : *pds = *temp;
768 66 : pds[1] = temp[1];
769 66 : pds[2] = temp[2];
770 66 : if (VSIFReadL(pds + 3, sizeof (char), sectLen - 3, fp) + 3 != sectLen) {
771 0 : errSprintf ("Ran out of file.\n");
772 0 : free (pds);
773 0 : return -1;
774 : }
775 :
776 66 : if (ReadGrib1Sect1 (pds, sectLen, gribLen, &curLoc, &pdsMeta, &f_gds, &gridID,
777 66 : &f_bms, &DSF, ¢er, &subcenter) != 0) {
778 0 : preErrSprintf ("Inside GRIB1_Inventory\n");
779 0 : free (pds);
780 0 : return -1;
781 : }
782 66 : free (pds);
783 66 : inv->refTime = pdsMeta.refTime;
784 66 : inv->validTime = pdsMeta.validTime;
785 66 : inv->foreSec = inv->validTime - inv->refTime;
786 66 : GRIB1_Table2LookUp (&(pdsMeta), &varName, &varComment, &varUnit,
787 : &convert, center, subcenter);
788 66 : inv->element = (char *) malloc ((1 + strlen (varName)) * sizeof (char));
789 66 : strcpy (inv->element, varName);
790 66 : inv->unitName = (char *) malloc ((1 + 2 + strlen (varUnit)) *
791 : sizeof (char));
792 66 : snprintf (inv->unitName, (1 + 2 + strlen (varUnit)) *
793 : sizeof (char), "[%s]", varUnit);
794 66 : inv->comment = (char *) malloc ((1 + strlen (varComment) +
795 66 : strlen (varUnit) + 2 + 1) *
796 : sizeof (char));
797 66 : snprintf (inv->comment, (1 + strlen (varComment) +
798 66 : strlen (varUnit) + 2 + 1) *
799 : sizeof (char),
800 : "%s [%s]", varComment, varUnit);
801 :
802 66 : GRIB1_Table3LookUp (&(pdsMeta), &(inv->shortFstLevel),
803 : &(inv->longFstLevel));
804 :
805 : /* Get to the end of the GRIB1 message. */
806 : /* (inventory.c : GRIB2Inventory), is responsible for this. */
807 : /* fseek (fp, gribLen - sectLen, SEEK_CUR); */
808 66 : return 0;
809 : }
810 :
811 0 : int GRIB1_RefTime (VSILFILE *fp, uInt4 gribLen, double *refTime)
812 : {
813 : uChar temp[3]; /* Used to determine the section length. */
814 : uInt4 sectLen; /* Length in bytes of the current section. */
815 : uChar *pds; /* The part of the message dealing with the PDS. */
816 : pdsG1Type pdsMeta; /* The pds parsed into a usable data structure. */
817 : char f_gds; /* flag if there is a gds section. */
818 : char f_bms; /* flag if there is a bms section. */
819 : short int DSF; /* Decimal Scale Factor for unpacking the data. */
820 : uChar gridID; /* Which GDS specs to use. */
821 : uInt4 curLoc; /* Where we are in the current GRIB message. */
822 : unsigned short int center; /* The Center that created the data */
823 : unsigned short int subcenter; /* The Sub Center that created the data */
824 :
825 0 : curLoc = 8;
826 0 : if (VSIFReadL (temp, sizeof (char), 3, fp) != 3) {
827 0 : errSprintf ("Ran out of file.\n");
828 0 : return -1;
829 : }
830 0 : sectLen = GRIB_UNSIGN_INT3 (*temp, temp[1], temp[2]);
831 0 : if (curLoc + sectLen > gribLen) {
832 0 : errSprintf ("Ran out of data in PDS (GRIB1_Inventory)\n");
833 0 : return -1;
834 : }
835 0 : pds = (uChar *) malloc (sectLen * sizeof (uChar));
836 0 : *pds = *temp;
837 0 : pds[1] = temp[1];
838 0 : pds[2] = temp[2];
839 0 : if (VSIFReadL (pds + 3, sizeof (char), sectLen - 3, fp) + 3 != sectLen) {
840 0 : errSprintf ("Ran out of file.\n");
841 0 : free (pds);
842 0 : return -1;
843 : }
844 :
845 0 : if (ReadGrib1Sect1 (pds, sectLen, gribLen, &curLoc, &pdsMeta, &f_gds, &gridID,
846 0 : &f_bms, &DSF, ¢er, &subcenter) != 0) {
847 0 : preErrSprintf ("Inside GRIB1_Inventory\n");
848 0 : free (pds);
849 0 : return -1;
850 : }
851 0 : free (pds);
852 :
853 0 : *refTime = pdsMeta.refTime;
854 :
855 : /* Get to the end of the GRIB1 message. */
856 : /* (inventory.c : GRIB2Inventory), is responsible for this. */
857 : /* fseek (fp, gribLen - sectLen, SEEK_CUR); */
858 0 : return 0;
859 : }
860 :
861 : /*****************************************************************************
862 : * ReadGrib1Sect2() --
863 : *
864 : * Arthur Taylor / MDL
865 : *
866 : * PURPOSE
867 : * Parses the GRIB1 "Grid Definition Section" or section 2, filling out
868 : * the gdsMeta data structure.
869 : *
870 : * ARGUMENTS
871 : * gds = The compressed part of the message dealing with "GDS". (Input)
872 : * gribLen = The total length of the GRIB1 message. (Input)
873 : * curLoc = Current location in the GRIB1 message. (Output)
874 : * gdsMeta = The filled out gdsMeta data structure. (Output)
875 : *
876 : * FILES/DATABASES: None
877 : *
878 : * RETURNS: int (could use errSprintf())
879 : * 0 = OK
880 : * -1 = gribLen is too small.
881 : * -2 = unexpected values in gds.
882 : *
883 : * HISTORY
884 : * 4/2003 Arthur Taylor (MDL/RSIS): Created.
885 : * 12/2003 AAT: adas data encoder seems to have # of vertical data = 1, but
886 : * parameters of vertical data = 255, which doesn't make sense.
887 : * Changed the error from "fatal" to a warning in debug mode.
888 : * 6/2004 AAT: Modified to allow "extended" lat/lon grids (i.e. stretched or
889 : * stretched and rotated).
890 : *
891 : * NOTES
892 : *****************************************************************************
893 : */
894 66 : static int ReadGrib1Sect2 (uChar *gds, uInt4 gribLen, uInt4 *curLoc,
895 : gdsType *gdsMeta)
896 : {
897 : uInt4 sectLen; /* Length in bytes of the current section. */
898 : int gridType; /* Which type of grid. (see enumerated types). */
899 66 : double unit = 1e-3; /* Used for converting to the correct unit. */
900 : uInt4 uli_temp; /* Used for reading a GRIB1 float. */
901 : int i;
902 : int f_allZero; /* Used to find out if the "lat/lon" extension part
903 : * is all 0 hence missing. */
904 : int f_allOne; /* Used to find out if the "lat/lon" extension part
905 : * is all 1 hence missing. */
906 :
907 66 : if( gribLen < *curLoc || gribLen - *curLoc < 6 ) {
908 0 : errSprintf ("Ran out of data in GDS (GRIB 1 Section 2)\n");
909 0 : return -1;
910 : }
911 66 : sectLen = GRIB_UNSIGN_INT3 (*gds, gds[1], gds[2]);
912 : #ifdef DEBUG
913 : /*
914 : printf ("Section 2 length = %ld\n", sectLen);
915 : */
916 : #endif
917 66 : *curLoc += sectLen;
918 66 : if (*curLoc > gribLen) {
919 0 : errSprintf ("Ran out of data in GDS (GRIB 1 Section 2)\n");
920 0 : return -1;
921 : }
922 66 : gds += 3;
923 : /*
924 : #ifdef DEBUG
925 : if ((*gds != 0) || (gds[1] != 255)) {
926 : printf ("GRIB1 GDS: Expect (NV = 0) != %d, (PV = 255) != %d\n",
927 : *gds, gds[1]);
928 : errSprintf ("SectLen == %ld\n", sectLen);
929 : errSprintf ("GridType == %d\n", gds[2]);
930 : }
931 : #endif
932 : */
933 : #ifdef DEBUG
934 : /*if (gds[1] != 255) {
935 : printf ("\n\tCaution: GRIB1 GDS: FOR ALL NWS products, PV should be "
936 : "255 rather than %u\n", gds[1]);
937 : }*/
938 : #endif
939 66 : if ((gds[1] != 255) && (gds[1] > 6)) {
940 : //errSprintf ("GRIB1 GDS: Expect PV = 255 != %d\n", gds[1]);
941 : //return -2;
942 : }
943 66 : gds += 2;
944 66 : gridType = *(gds++);
945 66 : switch (gridType) {
946 65 : case GB1S2_LATLON: // Latitude/Longitude Grid
947 : case GB1S2_GAUSSIAN_LATLON: // Gaussian Latitude/Longitude
948 : case GB1S2_ROTATED: // Rotated Latitude/Longitude
949 65 : if( gridType == GB1S2_ROTATED && (sectLen < 42)) {
950 0 : errSprintf ("For Rotated LatLon GDS, should have at least 42 bytes "
951 : "of data\n");
952 0 : return -1;
953 : }
954 65 : if ((sectLen < 32)) {
955 0 : errSprintf ("For LatLon GDS, should have at least 32 bytes "
956 : "of data\n");
957 0 : return -1;
958 : }
959 : switch(gridType) {
960 0 : case GB1S2_GAUSSIAN_LATLON:
961 0 : gdsMeta->projType = GS3_GAUSSIAN_LATLON;
962 0 : break;
963 1 : case GB1S2_ROTATED:
964 1 : gdsMeta->projType = GS3_ROTATED_LATLON;
965 1 : break;
966 64 : default:
967 64 : gdsMeta->projType = GS3_LATLON;
968 64 : break;
969 : }
970 65 : gdsMeta->orientLon = 0;
971 65 : gdsMeta->meshLat = 0;
972 65 : gdsMeta->scaleLat1 = 0;
973 65 : gdsMeta->scaleLat2 = 0;
974 65 : gdsMeta->southLat = 0;
975 65 : gdsMeta->southLon = 0;
976 65 : gdsMeta->center = 0;
977 :
978 65 : gdsMeta->Nx = GRIB_UNSIGN_INT2 (*gds, gds[1]);
979 65 : if( gdsMeta->Nx == 65535 )
980 : {
981 : /* https://rda.ucar.edu/docs/formats/grib/gribdoc/llgrid.html */
982 0 : errSprintf ("Quasi rectangular grid with varying number of grids points per row are not supported\n");
983 0 : return -1;
984 : }
985 65 : gds += 2;
986 65 : gdsMeta->Ny = GRIB_UNSIGN_INT2 (*gds, gds[1]);
987 65 : gds += 2;
988 65 : gdsMeta->lat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
989 65 : gds += 3;
990 65 : gdsMeta->lon1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
991 65 : gds += 3;
992 :
993 65 : gdsMeta->resFlag = *(gds++);
994 65 : if (gdsMeta->resFlag & 0x40) {
995 0 : gdsMeta->f_sphere = 0;
996 0 : gdsMeta->majEarth = 6378.160;
997 0 : gdsMeta->minEarth = 6356.775;
998 : } else {
999 65 : gdsMeta->f_sphere = 1;
1000 65 : gdsMeta->majEarth = 6367.47;
1001 65 : gdsMeta->minEarth = 6367.47;
1002 : }
1003 :
1004 65 : gdsMeta->lat2 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1005 65 : gds += 3;
1006 65 : gdsMeta->lon2 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1007 65 : gds += 3;
1008 65 : gdsMeta->Dx = GRIB_UNSIGN_INT2 (*gds, gds[1]) * unit;
1009 65 : gds += 2;
1010 65 : if (gridType == GB1S2_GAUSSIAN_LATLON) {
1011 0 : int np = GRIB_UNSIGN_INT2 (*gds, gds[1]); /* parallels between a pole and the equator */
1012 0 : if( np == 0 )
1013 : {
1014 0 : errSprintf ("Invalid Gaussian LatLon\n" );
1015 0 : return -1;
1016 : }
1017 0 : gdsMeta->Dy = 90.0 / np;
1018 : } else
1019 65 : gdsMeta->Dy = GRIB_UNSIGN_INT2 (*gds, gds[1]) * unit;
1020 65 : gds += 2;
1021 65 : gdsMeta->scan = *gds;
1022 65 : gdsMeta->f_typeLatLon = 0;
1023 : #ifdef DEBUG
1024 : /*
1025 : printf ("sectLen %ld\n", sectLen);
1026 : */
1027 : #endif
1028 65 : if (gridType == GB1S2_ROTATED && sectLen >= 42) {
1029 : /* Check if all 0's or all 1's, which means f_typeLatLon == 0 */
1030 1 : f_allZero = 1;
1031 1 : f_allOne = 1;
1032 11 : for (i = 0; i < 10; i++) {
1033 10 : if (gds[i] != 0)
1034 6 : f_allZero = 0;
1035 10 : if (gds[i] != 255)
1036 10 : f_allOne = 0;
1037 : }
1038 1 : if (!f_allZero && !f_allOne) {
1039 1 : gdsMeta->f_typeLatLon = 3;
1040 1 : gds += 5;
1041 1 : gdsMeta->southLat = (GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) *
1042 : unit);
1043 1 : gds += 3;
1044 1 : gdsMeta->southLon = (GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) *
1045 : unit);
1046 1 : gds += 3;
1047 1 : MEMCPY_BIG (&uli_temp, gds, sizeof (sInt4));
1048 1 : gdsMeta->angleRotate = fval_360 (uli_temp);
1049 : }
1050 : }
1051 : #if 0
1052 : else if (gridType == GB1S2_STRETCHED && sectLen >= 42) {
1053 : gds += 5;
1054 : /* Check if all 0's or all 1's, which means f_typeLatLon == 0 */
1055 : f_allZero = 1;
1056 : f_allOne = 1;
1057 : for (i = 0; i < 20; i++) {
1058 : if (gds[i] != 0)
1059 : f_allZero = 0;
1060 : if (gds[i] != 255)
1061 : f_allOne = 0;
1062 : }
1063 : if (!f_allZero && !f_allOne) {
1064 : gdsMeta->f_typeLatLon = 1;
1065 : gdsMeta->poleLat = (GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) *
1066 : unit);
1067 : gds += 3;
1068 : gdsMeta->poleLon = (GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) *
1069 : unit);
1070 : gds += 3;
1071 : MEMCPY_BIG (&uli_temp, gds, sizeof (sInt4));
1072 : gdsMeta->stretchFactor = fval_360 (uli_temp);
1073 : }
1074 : else if (gridType == GB1S2_ROTATED_STRETCHED && sectLen >= 52) {
1075 : gds += 5;
1076 : /* Check if all 0's or all 1's, which means f_typeLatLon == 0 */
1077 : f_allZero = 1;
1078 : f_allOne = 1;
1079 : for (i = 0; i < 20; i++) {
1080 : if (gds[i] != 0)
1081 : f_allZero = 0;
1082 : if (gds[i] != 255)
1083 : f_allOne = 0;
1084 : }
1085 : if (!f_allZero && !f_allOne) {
1086 : gdsMeta->f_typeLatLon = 2;
1087 : gdsMeta->southLat = (GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) *
1088 : unit);
1089 : gds += 3;
1090 : gdsMeta->southLon = (GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) *
1091 : unit);
1092 : gds += 3;
1093 : MEMCPY_BIG (&uli_temp, gds, sizeof (sInt4));
1094 : gdsMeta->angleRotate = fval_360 (uli_temp);
1095 : gds += 4;
1096 : gdsMeta->poleLat = (GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) *
1097 : unit);
1098 : gds += 3;
1099 : gdsMeta->poleLon = (GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) *
1100 : unit);
1101 : gds += 3;
1102 : MEMCPY_BIG (&uli_temp, gds, sizeof (sInt4));
1103 : gdsMeta->stretchFactor = fval_360 (uli_temp);
1104 : }
1105 : #endif
1106 :
1107 65 : break;
1108 :
1109 1 : case GB1S2_POLAR:
1110 1 : if (sectLen < 32) {
1111 0 : errSprintf ("For Polar GDS, should have 32 bytes of data\n");
1112 0 : return -1;
1113 : }
1114 1 : gdsMeta->projType = GS3_POLAR;
1115 1 : gdsMeta->lat2 = 0;
1116 1 : gdsMeta->lon2 = 0;
1117 1 : gdsMeta->southLat = 0;
1118 1 : gdsMeta->southLon = 0;
1119 :
1120 1 : gdsMeta->Nx = GRIB_UNSIGN_INT2 (*gds, gds[1]);
1121 1 : gds += 2;
1122 1 : gdsMeta->Ny = GRIB_UNSIGN_INT2 (*gds, gds[1]);
1123 1 : gds += 2;
1124 1 : gdsMeta->lat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1125 1 : gds += 3;
1126 1 : gdsMeta->lon1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1127 1 : gds += 3;
1128 :
1129 1 : gdsMeta->resFlag = *(gds++);
1130 1 : if (gdsMeta->resFlag & 0x40) {
1131 0 : gdsMeta->f_sphere = 0;
1132 0 : gdsMeta->majEarth = 6378.160;
1133 0 : gdsMeta->minEarth = 6356.775;
1134 : } else {
1135 1 : gdsMeta->f_sphere = 1;
1136 1 : gdsMeta->majEarth = 6367.47;
1137 1 : gdsMeta->minEarth = 6367.47;
1138 : }
1139 :
1140 1 : gdsMeta->orientLon = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1141 1 : gds += 3;
1142 1 : gdsMeta->Dx = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]);
1143 1 : gds += 3;
1144 1 : gdsMeta->Dy = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]);
1145 1 : gds += 3;
1146 1 : gdsMeta->meshLat = 60; /* Depends on hemisphere. */
1147 1 : gdsMeta->center = *(gds++);
1148 1 : if (gdsMeta->center & GRIB2BIT_1) {
1149 : /* South polar stereographic. */
1150 1 : gdsMeta->scaleLat1 = gdsMeta->scaleLat2 = -90;
1151 1 : gdsMeta->meshLat = -gdsMeta->meshLat;
1152 : } else {
1153 : /* North polar stereographic. */
1154 0 : gdsMeta->scaleLat1 = gdsMeta->scaleLat2 = 90;
1155 : }
1156 1 : gdsMeta->scan = *gds;
1157 1 : break;
1158 :
1159 0 : case GB1S2_LAMBERT:
1160 0 : if (sectLen < 42) {
1161 0 : errSprintf ("For Lambert GDS, should have 42 bytes of data\n");
1162 0 : return -1;
1163 : }
1164 0 : gdsMeta->projType = GS3_LAMBERT;
1165 0 : gdsMeta->lat2 = 0;
1166 0 : gdsMeta->lon2 = 0;
1167 :
1168 0 : gdsMeta->Nx = GRIB_UNSIGN_INT2 (*gds, gds[1]);
1169 0 : gds += 2;
1170 0 : gdsMeta->Ny = GRIB_UNSIGN_INT2 (*gds, gds[1]);
1171 0 : gds += 2;
1172 0 : gdsMeta->lat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1173 0 : gds += 3;
1174 0 : gdsMeta->lon1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1175 0 : gds += 3;
1176 :
1177 0 : gdsMeta->resFlag = *(gds++);
1178 0 : if (gdsMeta->resFlag & 0x40) {
1179 0 : gdsMeta->f_sphere = 0;
1180 0 : gdsMeta->majEarth = 6378.160;
1181 0 : gdsMeta->minEarth = 6356.775;
1182 : } else {
1183 0 : gdsMeta->f_sphere = 1;
1184 0 : gdsMeta->majEarth = 6367.47;
1185 0 : gdsMeta->minEarth = 6367.47;
1186 : }
1187 :
1188 0 : gdsMeta->orientLon = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1189 0 : gds += 3;
1190 0 : gdsMeta->Dx = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]);
1191 0 : gds += 3;
1192 0 : gdsMeta->Dy = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]);
1193 0 : gds += 3;
1194 0 : gdsMeta->center = *(gds++);
1195 0 : gdsMeta->scan = *(gds++);
1196 0 : gdsMeta->scaleLat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1197 0 : gds += 3;
1198 0 : gdsMeta->scaleLat2 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1199 0 : gds += 3;
1200 0 : gdsMeta->meshLat = gdsMeta->scaleLat1;
1201 0 : gdsMeta->southLat = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1202 0 : gds += 3;
1203 0 : gdsMeta->southLon = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1204 0 : break;
1205 :
1206 0 : case GB1S2_MERCATOR:
1207 0 : if (sectLen < 42) {
1208 0 : errSprintf ("For Mercator GDS, should have 42 bytes of data\n");
1209 0 : return -1;
1210 : }
1211 0 : gdsMeta->projType = GS3_MERCATOR;
1212 0 : gdsMeta->southLat = 0;
1213 0 : gdsMeta->southLon = 0;
1214 0 : gdsMeta->orientLon = 0;
1215 0 : gdsMeta->center = 0;
1216 :
1217 0 : gdsMeta->Nx = GRIB_UNSIGN_INT2 (*gds, gds[1]);
1218 0 : gds += 2;
1219 0 : gdsMeta->Ny = GRIB_UNSIGN_INT2 (*gds, gds[1]);
1220 0 : gds += 2;
1221 0 : gdsMeta->lat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1222 0 : gds += 3;
1223 0 : gdsMeta->lon1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1224 0 : gds += 3;
1225 :
1226 0 : gdsMeta->resFlag = *(gds++);
1227 0 : if (gdsMeta->resFlag & 0x40) {
1228 0 : gdsMeta->f_sphere = 0;
1229 0 : gdsMeta->majEarth = 6378.160;
1230 0 : gdsMeta->minEarth = 6356.775;
1231 : } else {
1232 0 : gdsMeta->f_sphere = 1;
1233 0 : gdsMeta->majEarth = 6367.47;
1234 0 : gdsMeta->minEarth = 6367.47;
1235 : }
1236 :
1237 0 : gdsMeta->lat2 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1238 0 : gds += 3;
1239 0 : gdsMeta->lon2 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1240 0 : gds += 3;
1241 0 : gdsMeta->scaleLat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1242 0 : gds += 3;
1243 0 : gdsMeta->scaleLat2 = gdsMeta->scaleLat1;
1244 0 : gdsMeta->meshLat = gdsMeta->scaleLat1;
1245 : /* Reserved set to 0. */
1246 0 : gds++;
1247 0 : gdsMeta->scan = *(gds++);
1248 0 : gdsMeta->Dx = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]);
1249 0 : gds += 3;
1250 0 : gdsMeta->Dy = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]);
1251 0 : break;
1252 :
1253 0 : default:
1254 0 : errSprintf ("Grid projection number is %d\n", gridType);
1255 0 : errSprintf ("Don't know how to handle this grid projection.\n");
1256 0 : return -2;
1257 : }
1258 66 : gdsMeta->numPts = gdsMeta->Nx * gdsMeta->Ny;
1259 : #ifdef DEBUG
1260 : /*
1261 : printf ("NumPts = %ld\n", gdsMeta->numPts);
1262 : printf ("Nx = %ld, Ny = %ld\n", gdsMeta->Nx, gdsMeta->Ny);
1263 : */
1264 : #endif
1265 66 : return 0;
1266 : }
1267 :
1268 : /*****************************************************************************
1269 : * ReadGrib1Sect3() --
1270 : *
1271 : * Arthur Taylor / MDL
1272 : *
1273 : * PURPOSE
1274 : * Parses the GRIB1 "Bit Map Section" or section 3, filling out the bitmap
1275 : * as needed.
1276 : *
1277 : * ARGUMENTS
1278 : * bms = The compressed part of the message dealing with "BMS". (Input)
1279 : * gribLen = The total length of the GRIB1 message. (Input)
1280 : * curLoc = Current location in the GRIB1 message. (Output)
1281 : * pBitmap = Pointer to the extracted bitmap. (Output)
1282 : * NxNy = The total size of the grid. (Input)
1283 : *
1284 : * FILES/DATABASES: None
1285 : *
1286 : * RETURNS: int (could use errSprintf())
1287 : * 0 = OK
1288 : * -1 = gribLen is too small.
1289 : * -2 = unexpected values in bms.
1290 : *
1291 : * HISTORY
1292 : * 4/2003 Arthur Taylor (MDL/RSIS): Created.
1293 : *
1294 : * NOTES
1295 : *****************************************************************************
1296 : */
1297 51 : static int ReadGrib1Sect3 (uChar *bms, uInt4 gribLen, uInt4 *curLoc,
1298 : uChar **pBitmap, uInt4 NxNy)
1299 : {
1300 : uInt4 sectLen; /* Length in bytes of the current section. */
1301 : short int numeric; /* Determine if this is a predefined bitmap */
1302 : uChar bits; /* Used to locate which bit we are currently using. */
1303 : uInt4 i; /* Helps traverse the bitmap. */
1304 :
1305 51 : *pBitmap = nullptr;
1306 :
1307 51 : uInt4 bmsRemainingSize = gribLen - *curLoc;
1308 51 : if( bmsRemainingSize < 6 )
1309 : {
1310 0 : errSprintf ("Ran out of data in BMS (GRIB 1 Section 3)\n");
1311 0 : return -1;
1312 : }
1313 51 : sectLen = GRIB_UNSIGN_INT3 (*bms, bms[1], bms[2]);
1314 : #ifdef DEBUG
1315 : /*
1316 : printf ("Section 3 length = %ld\n", sectLen);
1317 : */
1318 : #endif
1319 51 : *curLoc += sectLen;
1320 51 : if (*curLoc > gribLen) {
1321 0 : errSprintf ("Ran out of data in BMS (GRIB 1 Section 3)\n");
1322 0 : return -1;
1323 : }
1324 51 : bms += 3;
1325 : /* Assert: *bms currently points to number of unused bits at end of BMS. */
1326 51 : if (NxNy + *bms + 6 * 8 != sectLen * 8) {
1327 0 : errSprintf ("NxNy + # of unused bits != # of available bits\n");
1328 : // commented out to avoid unsigned integer overflow
1329 : // (sInt4) (NxNy + *bms), (sInt4) ((sectLen - 6) * 8));
1330 0 : return -2;
1331 : }
1332 51 : bms++;
1333 : /* Assert: Non-zero "numeric" means predefined bitmap. */
1334 51 : numeric = GRIB_UNSIGN_INT2 (*bms, bms[1]);
1335 51 : bms += 2;
1336 51 : if (numeric != 0) {
1337 0 : errSprintf ("Don't handle predefined bitmaps yet.\n");
1338 0 : return -2;
1339 : }
1340 51 : bmsRemainingSize -= 6;
1341 51 : if( bmsRemainingSize < (NxNy+7) / 8 )
1342 : {
1343 0 : errSprintf ("Ran out of data in BMS (GRIB 1 Section 3)\n");
1344 0 : return -1;
1345 : }
1346 51 : *pBitmap = (uChar*) malloc(NxNy);
1347 51 : auto bitmap = *pBitmap;
1348 51 : if( bitmap== nullptr )
1349 : {
1350 0 : errSprintf ("Ran out of memory in allocating bitmap (GRIB 1 Section 3)\n");
1351 0 : return -1;
1352 : }
1353 51 : bits = 0x80;
1354 156819 : for (i = 0; i < NxNy; i++) {
1355 156768 : *(bitmap++) = (*bms) & bits;
1356 156768 : bits = bits >> 1;
1357 156768 : if (bits == 0) {
1358 19566 : bms++;
1359 19566 : bits = 0x80;
1360 : }
1361 : }
1362 51 : return 0;
1363 : }
1364 :
1365 : #ifdef USE_UNPACKCMPLX
1366 : static int UnpackCmplx (uChar *bds, CPL_UNUSED uInt4 gribLen, CPL_UNUSED uInt4 *curLoc,
1367 : CPL_UNUSED short int DSF, CPL_UNUSED double *data, CPL_UNUSED grib_MetaData *meta,
1368 : CPL_UNUSED char f_bms, CPL_UNUSED uChar *bitmap, CPL_UNUSED double unitM,
1369 : CPL_UNUSED double unitB, CPL_UNUSED short int ESF, CPL_UNUSED double refVal,
1370 : uChar numBits, uChar f_octet14)
1371 : {
1372 : uInt4 secLen;
1373 : int N1;
1374 : int N2;
1375 : int P1;
1376 : int P2;
1377 : uChar octet14;
1378 : uChar f_maxtrixValues;
1379 : uChar f_secBitmap = 0;
1380 : uChar f_secValDiffWid;
1381 : int i;
1382 : uInt4 uli_temp; /* Used to store sInt4s (temporarily) */
1383 : uChar bufLoc; /* Keeps track of where to start getting more data
1384 : * out of the packed data stream. */
1385 : size_t numUsed; /* How many bytes were used in a given call to
1386 : * memBitRead. */
1387 : uChar *width;
1388 :
1389 : secLen = 11;
1390 : N1 = GRIB_UNSIGN_INT2 (bds[0], bds[1]);
1391 : octet14 = bds[2];
1392 : printf ("octet14, %d\n", octet14);
1393 : if (f_octet14) {
1394 : f_maxtrixValues = octet14 & GRIB2BIT_2;
1395 : f_secBitmap = octet14 & GRIB2BIT_3;
1396 : f_secValDiffWid = octet14 & GRIB2BIT_4;
1397 : printf ("f_matrixValues, f_secBitmap, f_secValeDiffWid %d %d %d\n",
1398 : f_maxtrixValues, f_secBitmap, f_secValDiffWid);
1399 : }
1400 : N2 = GRIB_UNSIGN_INT2 (bds[3], bds[4]);
1401 : P1 = GRIB_UNSIGN_INT2 (bds[5], bds[6]);
1402 : P2 = GRIB_UNSIGN_INT2 (bds[7], bds[8]);
1403 : printf ("N1 N2 P1 P2 : %d %d %d %d\n", N1, N2, P1, P2);
1404 : printf ("Reserved %u\n", bds[9]);
1405 : bds += 10;
1406 : secLen += 10;
1407 :
1408 : width = (uChar *) malloc (P1 * sizeof (uChar));
1409 :
1410 : for (i = 0; i < P1; i++) {
1411 : width[i] = *bds;
1412 : printf ("(Width %d %u)\n", i, width[i]);
1413 : bds++;
1414 : secLen++;
1415 : }
1416 : if (f_secBitmap) {
1417 : bufLoc = 8;
1418 : for (i = 0; i < P2; i++) {
1419 : memBitRead (&uli_temp, sizeof (sInt4), bds, 1, &bufLoc, &numUsed);
1420 : printf ("(%d %u) ", i, uli_temp);
1421 : if (numUsed != 0) {
1422 : printf ("\n");
1423 : bds += numUsed;
1424 : secLen++;
1425 : }
1426 : }
1427 : if (bufLoc != 8) {
1428 : bds++;
1429 : secLen++;
1430 : }
1431 : printf ("Observed Sec Len %u\n", secLen);
1432 : } else {
1433 : /* Jump over widths and secondary bitmap */
1434 : bds += (N1 - 21);
1435 : secLen += (N1 - 21);
1436 : }
1437 :
1438 : bufLoc = 8;
1439 : for (i = 0; i < P1; i++) {
1440 : memBitRead (&uli_temp, sizeof (sInt4), bds, numBits, &bufLoc, &numUsed);
1441 : printf ("(%d %u) (numUsed %ld numBits %d)", i, uli_temp,
1442 : (long) numUsed, numBits);
1443 : if (numUsed != 0) {
1444 : printf ("\n");
1445 : bds += numUsed;
1446 : secLen++;
1447 : }
1448 : }
1449 : if (bufLoc != 8) {
1450 : // cppcheck-suppress uselessAssignmentPtrArg
1451 : bds++;
1452 : secLen++;
1453 : }
1454 :
1455 : printf ("Observed Sec Len %u\n", secLen);
1456 : printf ("N2 = %d\n", N2);
1457 :
1458 : errSprintf ("Don't know how to handle Complex GRIB1 packing yet.\n");
1459 : free (width);
1460 : return -2;
1461 :
1462 : }
1463 : #endif /* USE_UNPACKCMPLX */
1464 :
1465 : /*****************************************************************************
1466 : * ReadGrib1Sect4() --
1467 : *
1468 : * Arthur Taylor / MDL
1469 : *
1470 : * PURPOSE
1471 : * Unpacks the "Binary Data Section" or section 4.
1472 : *
1473 : * ARGUMENTS
1474 : * bds = The compressed part of the message dealing with "BDS". (Input)
1475 : * gribLen = The total length of the GRIB1 message. (Input)
1476 : * curLoc = Current location in the GRIB1 message. (Input/Output)
1477 : * DSF = Decimal Scale Factor for unpacking the data. (Input)
1478 : * data = The extracted grid. (Output)
1479 : * meta = The meta data associated with the grid (Input/Output)
1480 : * f_bms = True if bitmap is to be used. (Input)
1481 : * bitmap = 0 if missing data, 1 if valid data. (Input)
1482 : * unitM = The M unit conversion value in equation y = Mx + B. (Input)
1483 : * unitB = The B unit conversion value in equation y = Mx + B. (Input)
1484 : *
1485 : * FILES/DATABASES: None
1486 : *
1487 : * RETURNS: int (could use errSprintf())
1488 : * 0 = OK
1489 : * -1 = gribLen is too small.
1490 : * -2 = unexpected values in bds.
1491 : *
1492 : * HISTORY
1493 : * 4/2003 Arthur Taylor (MDL/RSIS): Created
1494 : * 3/2004 AAT: Switched {# Pts * (# Bits in a Group) +
1495 : * # of unused bits != # of available bits} to a warning from an
1496 : * error.
1497 : *
1498 : * NOTES
1499 : * 1) See metaparse.c : ParseGrid()
1500 : * 2) Currently, only handles "Simple pack".
1501 : *****************************************************************************
1502 : */
1503 66 : static int ReadGrib1Sect4 (uChar *bds, uInt4 gribLen, uInt4 *curLoc,
1504 : short int DSF, double *data, grib_MetaData *meta,
1505 : char f_bms, uChar *bitmap, double unitM,
1506 : double unitB)
1507 : {
1508 : uInt4 sectLen; /* Length in bytes of the current section. */
1509 : short int ESF; /* Power of 2 scaling factor. */
1510 : uInt4 uli_temp; /* Used to store sInt4s (temporarily) */
1511 : double refVal; /* The reference value for the grid, also the minimum
1512 : * value. */
1513 : uChar numBits; /* # of bits for a single element of data. */
1514 : uChar numUnusedBit; /* # of extra bits at end of record. */
1515 : uChar f_spherHarm; /* Flag if data contains Spherical Harmonics. */
1516 : uChar f_cmplxPack; /* Flag if complex packing was used. */
1517 : #ifdef USE_UNPACKCMPLX
1518 : uChar f_octet14; /* Flag if octet 14 was used. */
1519 : #endif
1520 : uChar bufLoc; /* Keeps track of where to start getting more data
1521 : * out of the packed data stream. */
1522 : uChar f_convert; /* Determine if scan mode implies that we have to do
1523 : * manipulation as we read the grid to get desired
1524 : * internal scan mode. */
1525 : uInt4 i; /* Used to traverse the grid. */
1526 : size_t numUsed; /* How many bytes were used in a given call to
1527 : * memBitRead. */
1528 : double d_temp; /* Holds the extracted data until we put it in data */
1529 : sInt4 newIndex; /* Where to put the answer (primarily if f_convert) */
1530 : sInt4 x; /* Used to help compute newIndex , if f_convert. */
1531 : sInt4 y; /* Used to help compute newIndex , if f_convert. */
1532 : double resetPrim; /* If possible, used to reset the primary missing
1533 : * value from 9.999e20 to a reasonable # (9999) */
1534 :
1535 66 : if (meta->gds.Nx * meta->gds.Ny != meta->gds.numPts) {
1536 0 : errSprintf ("(Nx * Ny != numPts) ?? in BDS (GRIB 1 Section 4)\n");
1537 0 : return -2;
1538 : }
1539 66 : if( *curLoc >= gribLen )
1540 0 : return -1;
1541 :
1542 66 : uInt4 bdsRemainingSize = gribLen - *curLoc;
1543 66 : if( bdsRemainingSize < 3 )
1544 0 : return -1;
1545 :
1546 66 : sectLen = GRIB_UNSIGN_INT3 (*bds, bds[1], bds[2]);
1547 : #ifdef DEBUG
1548 : /*
1549 : printf ("Section 4 length = %ld\n", sectLen);
1550 : */
1551 : #endif
1552 66 : *curLoc += sectLen;
1553 66 : if (*curLoc > gribLen) {
1554 0 : errSprintf ("Ran out of data in BDS (GRIB 1 Section 4)\n");
1555 0 : return -1;
1556 : }
1557 66 : bds += 3;
1558 66 : bdsRemainingSize -= 3;
1559 :
1560 : /* Assert: bds now points to the main pack flag. */
1561 66 : if( bdsRemainingSize < 1 )
1562 0 : return -1;
1563 66 : f_spherHarm = (*bds) & GRIB2BIT_1;
1564 66 : f_cmplxPack = (*bds) & GRIB2BIT_2;
1565 66 : meta->gridAttrib.fieldType = (*bds) & GRIB2BIT_3;
1566 : #ifdef USE_UNPACKCMPLX
1567 : f_octet14 = (*bds) & GRIB2BIT_4;
1568 : #endif
1569 :
1570 66 : numUnusedBit = (*bds) & 0x0f;
1571 : #ifdef DEBUG
1572 : /*
1573 : printf ("bds byte flag = %d\n", *bds);
1574 : printf ("Number of unused bits = %d\n", numUnusedBit);
1575 : */
1576 : #endif
1577 66 : if (f_spherHarm) {
1578 0 : errSprintf ("Don't know how to handle Spherical Harmonics yet.\n");
1579 0 : return -2;
1580 : }
1581 : /*
1582 : if (f_octet14) {
1583 : errSprintf ("Don't know how to handle Octet 14 data yet.\n");
1584 : errSprintf ("bds byte flag = %d\n", *bds);
1585 : errSprintf ("bds byte: %d %d %d %d\n", f_spherHarm, f_cmplxPack,
1586 : meta->gridAttrib.fieldType, f_octet14);
1587 : return -2;
1588 : }
1589 : */
1590 66 : if (f_cmplxPack) {
1591 0 : meta->gridAttrib.packType = 2;
1592 : } else {
1593 66 : meta->gridAttrib.packType = 0;
1594 : }
1595 66 : bds++;
1596 66 : bdsRemainingSize --;
1597 :
1598 : /* Assert: bds now points to E (power of 2 scaling factor). */
1599 66 : if( bdsRemainingSize < 2 )
1600 0 : return -1;
1601 66 : ESF = GRIB_SIGN_INT2 (*bds, bds[1]);
1602 66 : bds += 2;
1603 66 : bdsRemainingSize -= 2;
1604 :
1605 66 : if( bdsRemainingSize < 4 )
1606 0 : return -1;
1607 66 : MEMCPY_BIG (&uli_temp, bds, sizeof (sInt4));
1608 66 : refVal = fval_360 (uli_temp);
1609 66 : bds += 4;
1610 66 : bdsRemainingSize -= 4;
1611 :
1612 : /* Assert: bds is now the number of bits in a group. */
1613 66 : if( bdsRemainingSize < 1 )
1614 0 : return -1;
1615 66 : numBits = *bds;
1616 : /*
1617 : #ifdef DEBUG
1618 : printf ("refValue %f numBits %d\n", refVal, numBits);
1619 : printf ("ESF %d DSF %d\n", ESF, DSF);
1620 : #endif
1621 : */
1622 66 : if (f_cmplxPack) {
1623 : #ifdef USE_UNPACKCMPLX
1624 : bds++;
1625 : bdsRemainingSize --;
1626 : return UnpackCmplx (bds, gribLen, curLoc, DSF, data, meta, f_bms,
1627 : bitmap, unitM, unitB, ESF, refVal, numBits,
1628 : f_octet14);
1629 : #else
1630 0 : errSprintf ("Don't know how to handle Complex GRIB1 packing yet.\n");
1631 0 : return -2;
1632 : #endif
1633 : }
1634 :
1635 66 : if (!f_bms && (
1636 15 : sectLen < 11 ||
1637 15 : (numBits > 0 && meta->gds.numPts > UINT_MAX / numBits) ||
1638 15 : (meta->gds.numPts * numBits > UINT_MAX - numUnusedBit))) {
1639 0 : printf ("numPts * (numBits in a Group) + # of unused bits != "
1640 : "# of available bits\n");
1641 : }
1642 66 : else if (!f_bms &&
1643 15 : (meta->gds.numPts * numBits + numUnusedBit) != (sectLen - 11) * 8) {
1644 0 : printf ("numPts * (numBits in a Group) + # of unused bits %d != "
1645 : "# of available bits %d\n",
1646 0 : (sInt4) (meta->gds.numPts * numBits + numUnusedBit),
1647 0 : (sInt4) ((sectLen - 11) * 8));
1648 : /*
1649 : errSprintf ("numPts * (numBits in a Group) + # of unused bits %ld != "
1650 : "# of available bits %ld\n",
1651 : (sInt4) (meta->gds.numPts * numBits + numUnusedBit),
1652 : (sInt4) ((sectLen - 11) * 8));
1653 : return -2;
1654 : */
1655 : }
1656 66 : if (numBits > 32) {
1657 0 : errSprintf ("The number of bits per number is larger than 32?\n");
1658 0 : return -2;
1659 : }
1660 66 : bds++;
1661 66 : bdsRemainingSize -= 1;
1662 :
1663 : /* Convert Units. */
1664 : {
1665 66 : double pow_10_DSF = pow (10.0, DSF);
1666 66 : if( pow_10_DSF == 0.0 ) {
1667 0 : errSprintf ("pow_10_DSF == 0.0\n");
1668 0 : return -2;
1669 : }
1670 66 : if (unitM == -10) {
1671 0 : meta->gridAttrib.min = pow (10.0, (refVal * pow (2.0, ESF) /
1672 : pow_10_DSF));
1673 : } else {
1674 : /* meta->gridAttrib.min = unitM * (refVal / pow (10.0, DSF)) + unitB; */
1675 66 : meta->gridAttrib.min = unitM * (refVal * pow (2.0, ESF) /
1676 66 : pow_10_DSF) + unitB;
1677 : }
1678 : }
1679 :
1680 66 : meta->gridAttrib.max = meta->gridAttrib.min;
1681 66 : meta->gridAttrib.f_maxmin = 1;
1682 66 : meta->gridAttrib.numMiss = 0;
1683 66 : if (refVal >= std::numeric_limits<float>::max() || CPLIsNan(refVal)) {
1684 0 : meta->gridAttrib.refVal = std::numeric_limits<float>::max();
1685 66 : } else if (refVal <= -std::numeric_limits<float>::max()) {
1686 0 : meta->gridAttrib.refVal = -std::numeric_limits<float>::max();
1687 : } else {
1688 66 : meta->gridAttrib.refVal = static_cast<float>(refVal);
1689 : }
1690 66 : meta->gridAttrib.ESF = ESF;
1691 66 : meta->gridAttrib.DSF = DSF;
1692 66 : bufLoc = 8;
1693 : /* Internally we use scan = 0100. Scan is usually 0100 but if need be, we
1694 : * can convert it. */
1695 66 : f_convert = ((meta->gds.scan & 0xe0) != 0x40);
1696 :
1697 66 : if (f_bms) {
1698 51 : meta->gridAttrib.f_miss = 1;
1699 51 : meta->gridAttrib.missPri = UNDEFINED;
1700 : }
1701 : else
1702 : {
1703 15 : meta->gridAttrib.f_miss = 0;
1704 : }
1705 :
1706 66 : if (f_bms) {
1707 : /*
1708 : #ifdef DEBUG
1709 : printf ("There is a bitmap?\n");
1710 : #endif
1711 : */
1712 : /* Start unpacking the data, assuming there is a bitmap. */
1713 156819 : for (i = 0; i < meta->gds.numPts; i++) {
1714 : /* Find the destination index. */
1715 156768 : if (f_convert) {
1716 : /* ScanIndex2XY returns value as if scan was 0100 */
1717 12936 : ScanIndex2XY (i, &x, &y, meta->gds.scan, meta->gds.Nx,
1718 12936 : meta->gds.Ny);
1719 12936 : newIndex = (x - 1) + (y - 1) * meta->gds.Nx;
1720 : } else {
1721 143832 : newIndex = i;
1722 : }
1723 : /* A 0 in bitmap means no data. A 1 in bitmap means data. */
1724 : // cppcheck-suppress nullPointer
1725 156768 : if (!bitmap[i]) {
1726 60042 : meta->gridAttrib.numMiss++;
1727 60042 : if( data ) data[newIndex] = UNDEFINED;
1728 : } else {
1729 96726 : if (numBits != 0) {
1730 96726 : if( ((int)bdsRemainingSize - 1) * 8 + bufLoc < (int)numBits )
1731 0 : return -1;
1732 96726 : memBitRead (&uli_temp, sizeof (sInt4), bds, numBits,
1733 : &bufLoc, &numUsed);
1734 96726 : assert( numUsed <= bdsRemainingSize );
1735 96726 : bds += numUsed;
1736 96726 : bdsRemainingSize -= (uInt4)numUsed;
1737 96726 : d_temp = (refVal + (uli_temp * pow (2.0, ESF))) / pow (10.0, DSF);
1738 : /* Convert Units. */
1739 96726 : if (unitM == -10) {
1740 0 : d_temp = pow (10.0, d_temp);
1741 : } else {
1742 96726 : d_temp = unitM * d_temp + unitB;
1743 : }
1744 96726 : if (meta->gridAttrib.max < d_temp) {
1745 645 : meta->gridAttrib.max = d_temp;
1746 : }
1747 96726 : if( data ) data[newIndex] = d_temp;
1748 : } else {
1749 : /* Assert: d_temp = unitM * refVal / pow (10.0,DSF) + unitB. */
1750 : /* Assert: min = unitM * refVal / pow (10.0, DSF) + unitB. */
1751 0 : if( data ) data[newIndex] = meta->gridAttrib.min;
1752 : }
1753 : }
1754 : }
1755 : /* Reset the missing value to UNDEFINED_PRIM if possible. If not
1756 : * possible, make sure UNDEFINED is outside the range. If UNDEFINED
1757 : * is_ in the range, choose max + 1 for missing. */
1758 51 : resetPrim = 0;
1759 51 : if ((meta->gridAttrib.max < UNDEFINED_PRIM) ||
1760 3 : (meta->gridAttrib.min > UNDEFINED_PRIM)) {
1761 48 : resetPrim = UNDEFINED_PRIM;
1762 3 : } else if ((meta->gridAttrib.max >= UNDEFINED) &&
1763 0 : (meta->gridAttrib.min <= UNDEFINED)) {
1764 0 : resetPrim = meta->gridAttrib.max + 1;
1765 : }
1766 51 : if (resetPrim != 0) {
1767 48 : meta->gridAttrib.missPri = resetPrim;
1768 : }
1769 51 : if (data != nullptr && resetPrim != 0) {
1770 43044 : for (i = 0; i < meta->gds.numPts; i++) {
1771 : /* Find the destination index. */
1772 43026 : if (f_convert) {
1773 : /* ScanIndex2XY returns value as if scan was 0100 */
1774 6006 : ScanIndex2XY (i, &x, &y, meta->gds.scan, meta->gds.Nx,
1775 6006 : meta->gds.Ny);
1776 6006 : newIndex = (x - 1) + (y - 1) * meta->gds.Nx;
1777 : } else {
1778 37020 : newIndex = i;
1779 : }
1780 : // cppcheck-suppress nullPointer
1781 43026 : if (!bitmap[i]) {
1782 12911 : data[newIndex] = resetPrim;
1783 : }
1784 : }
1785 : }
1786 :
1787 : } else {
1788 :
1789 15 : if( !data )
1790 8 : return 0;
1791 :
1792 : #ifdef DEBUG
1793 : /*
1794 : printf ("There is no bitmap?\n");
1795 : */
1796 : #endif
1797 :
1798 : /* Start unpacking the data, assuming there is NO bitmap. */
1799 4123 : for (i = 0; i < meta->gds.numPts; i++) {
1800 4116 : if (numBits != 0) {
1801 : /* Find the destination index. */
1802 4116 : if (f_convert) {
1803 : /* ScanIndex2XY returns value as if scan was 0100 */
1804 4116 : ScanIndex2XY (i, &x, &y, meta->gds.scan, meta->gds.Nx,
1805 4116 : meta->gds.Ny);
1806 4116 : newIndex = (x - 1) + (y - 1) * meta->gds.Nx;
1807 : } else {
1808 0 : newIndex = i;
1809 : }
1810 :
1811 4116 : if( ((int)bdsRemainingSize - 1) * 8 + bufLoc < (int)numBits )
1812 0 : return -1;
1813 4116 : memBitRead (&uli_temp, sizeof (sInt4), bds, numBits, &bufLoc,
1814 : &numUsed);
1815 4116 : assert( numUsed <= bdsRemainingSize );
1816 4116 : bds += numUsed;
1817 4116 : bdsRemainingSize -= (uInt4)numUsed;
1818 4116 : d_temp = (refVal + (uli_temp * pow (2.0, ESF))) / pow (10.0, DSF);
1819 :
1820 : #ifdef DEBUG
1821 : /*
1822 : if (i == 1) {
1823 : printf ("refVal %f, uli_temp %ld, ans %f\n", refVal, uli_temp,
1824 : d_temp);
1825 : printf ("numBits %d, bufLoc %d, numUsed %d\n", numBits,
1826 : bufLoc, numUsed);
1827 : }
1828 : */
1829 : #endif
1830 :
1831 : /* Convert Units. */
1832 4116 : if (unitM == -10) {
1833 0 : d_temp = pow (10.0, d_temp);
1834 : } else {
1835 4116 : d_temp = unitM * d_temp + unitB;
1836 : }
1837 4116 : if (meta->gridAttrib.max < d_temp) {
1838 247 : meta->gridAttrib.max = d_temp;
1839 : }
1840 4116 : data[newIndex] = d_temp;
1841 : } else {
1842 : /* Assert: whole array = unitM * refVal + unitB. */
1843 : /* Assert: *min = unitM * refVal + unitB. */
1844 0 : data[i] = meta->gridAttrib.min;
1845 : }
1846 : }
1847 : }
1848 58 : return 0;
1849 : }
1850 :
1851 : /*****************************************************************************
1852 : * ReadGrib1Record() --
1853 : *
1854 : * Arthur Taylor / MDL
1855 : *
1856 : * PURPOSE
1857 : * Reads in a GRIB1 message, and parses the data into various data
1858 : * structures, for use with other code.
1859 : *
1860 : * ARGUMENTS
1861 : * fp = An opened GRIB2 file already at the correct message. (Input)
1862 : * f_unit = 0 use GRIB2 units, 1 use English, 2 use metric. (Input)
1863 : * Grib_Data = The read in GRIB2 grid. (Output)
1864 : * grib_DataLen = Size of Grib_Data. (Output)
1865 : * meta = A filled in meta structure (Output)
1866 : * IS = The structure containing all the arrays that the
1867 : * unpacker uses (Output)
1868 : * sect0 = Already read in section 0 data. (Input)
1869 : * gribLen = Length of the GRIB1 message. (Input)
1870 : * majEarth = Used to override the GRIB major axis of earth. (Input)
1871 : * minEarth = Used to override the GRIB minor axis of earth. (Input)
1872 : *
1873 : * FILES/DATABASES:
1874 : * An already opened file pointing to the desired GRIB1 message.
1875 : *
1876 : * RETURNS: int (could use errSprintf())
1877 : * 0 = OK
1878 : * -1 = Problems reading in the PDS.
1879 : * -2 = Problems reading in the GDS.
1880 : * -3 = Problems reading in the BMS.
1881 : * -4 = Problems reading in the BDS.
1882 : * -5 = Problems reading the closing section.
1883 : *
1884 : * HISTORY
1885 : * 4/2003 Arthur Taylor (MDL/RSIS): Created
1886 : * 5/2003 AAT: Was not updating offset. It should be updated by
1887 : * calling routine anyways, so I got rid of the parameter.
1888 : * 7/2003 AAT: Allowed user to override the radius of earth.
1889 : * 8/2003 AAT: Found a memory Leak (Had been setting unitName to NULL).
1890 : * 2/2004 AAT: Added maj/min earth override.
1891 : * 3/2004 AAT: Added ability to change units.
1892 : *
1893 : * NOTES
1894 : * 1) Could also compare GDS with the one specified by gridID
1895 : * 2) Could add gridID support.
1896 : * 3) Should add unitM / unitB support.
1897 : *****************************************************************************
1898 : */
1899 66 : int ReadGrib1Record (VSILFILE *fp, sChar f_unit, double **Grib_Data,
1900 : uInt4 *grib_DataLen, grib_MetaData *meta,
1901 : IS_dataType *IS, sInt4 sect0[SECT0LEN_WORD],
1902 : uInt4 gribLen, double majEarth, double minEarth)
1903 : {
1904 : sInt4 nd5; /* Size of grib message rounded up to the nearest *
1905 : * sInt4. */
1906 : uChar *c_ipack; /* A char ptr to the message stored in IS->ipack */
1907 : uInt4 curLoc; /* Current location in the GRIB message. */
1908 : char f_gds; /* flag if there is a gds section. */
1909 : char f_bms; /* flag if there is a bms section. */
1910 66 : double *grib_Data = nullptr; /* A pointer to Grib_Data for ease of manipulation. */
1911 66 : uChar *bitmap = nullptr; /* A char field (0=noData, 1=data) set up in BMS. */
1912 : short int DSF; /* Decimal Scale Factor for unpacking the data. */
1913 66 : double unitM = 1; /* M in y = Mx + B, for unit conversion. */
1914 66 : double unitB = 0; /* B in y = Mx + B, for unit conversion. */
1915 : uChar gridID; /* Which GDS specs to use. */
1916 : const char *varName; /* The name of the data stored in the grid. */
1917 : const char *varComment; /* Extra comments about the data stored in grid. */
1918 : const char *varUnit; /* Holds the name of the unit [K] [%] .. etc */
1919 : sInt4 li_temp; /* Used to make sure section 5 is 7777. */
1920 : char unitName[15]; /* Holds the string name of the current unit. */
1921 : int unitLen; /* String length of string name of current unit. */
1922 :
1923 : /* Make room for entire message, and read it in. */
1924 : /* nd5 needs to be gribLen in (sInt4) units rounded up. */
1925 66 : nd5 = (gribLen + 3) / 4;
1926 66 : if (nd5 > IS->ipackLen) {
1927 66 : if( gribLen > 100 * 1024 * 1024 )
1928 : {
1929 0 : vsi_l_offset curPos = VSIFTellL(fp);
1930 0 : VSIFSeekL(fp, 0, SEEK_END);
1931 0 : vsi_l_offset fileSize = VSIFTellL(fp);
1932 0 : VSIFSeekL(fp, curPos, SEEK_SET);
1933 0 : if( fileSize < gribLen )
1934 : {
1935 0 : errSprintf("File too short");
1936 0 : return -1;
1937 : }
1938 : }
1939 66 : sInt4* newipack = (sInt4 *) realloc ((void *) (IS->ipack),
1940 66 : nd5 * sizeof (sInt4));
1941 66 : if( newipack == nullptr )
1942 : {
1943 0 : errSprintf ("Out of memory\n");
1944 0 : return -1;
1945 : }
1946 66 : IS->ipackLen = nd5;
1947 66 : IS->ipack = newipack;
1948 : }
1949 66 : c_ipack = (uChar *) IS->ipack;
1950 : /* Init last sInt4 to 0, to make sure that the padded bytes are 0. */
1951 66 : IS->ipack[nd5 - 1] = 0;
1952 : /* Init first 2 sInt4 to sect0. */
1953 66 : memcpy (c_ipack, sect0, SECT0LEN_WORD * 2);
1954 : /* Read in the rest of the message. */
1955 132 : if (VSIFReadL (c_ipack + SECT0LEN_WORD * 2, sizeof (char),
1956 66 : (gribLen - SECT0LEN_WORD * 2), fp) + SECT0LEN_WORD * 2 != gribLen) {
1957 0 : errSprintf ("Ran out of file\n");
1958 0 : return -1;
1959 : }
1960 :
1961 : /* Preceding was in degrib2, next part is specific to GRIB1. */
1962 66 : curLoc = 8;
1963 66 : if (ReadGrib1Sect1 (c_ipack + curLoc, gribLen - curLoc, gribLen, &curLoc, &(meta->pds1),
1964 : &f_gds, &gridID, &f_bms, &DSF, &(meta->center),
1965 66 : &(meta->subcenter)) != 0) {
1966 0 : preErrSprintf ("Inside ReadGrib1Record\n");
1967 0 : return -1;
1968 : }
1969 :
1970 : /* Get the Grid Definition Section. */
1971 66 : if (f_gds) {
1972 66 : if (ReadGrib1Sect2 (c_ipack + curLoc, gribLen, &curLoc,
1973 66 : &(meta->gds)) != 0) {
1974 0 : preErrSprintf ("Inside ReadGrib1Record\n");
1975 0 : return -2;
1976 : }
1977 : /* Could also compare GDS with the one specified by gridID? */
1978 : } else {
1979 0 : errSprintf ("Don't know how to handle a gridID lookup yet.\n");
1980 0 : return -2;
1981 : }
1982 66 : meta->pds1.gridID = gridID;
1983 : /* Allow data originating from NCEP to be 6371.2 by default. */
1984 66 : if (meta->center == NMC) {
1985 64 : if (meta->gds.majEarth == 6367.47) {
1986 64 : meta->gds.f_sphere = 1;
1987 64 : meta->gds.majEarth = 6371.2;
1988 64 : meta->gds.minEarth = 6371.2;
1989 : }
1990 : }
1991 66 : if ((majEarth > 6300) && (majEarth < 6400)) {
1992 0 : if ((minEarth > 6300) && (minEarth < 6400)) {
1993 0 : meta->gds.f_sphere = 0;
1994 0 : meta->gds.majEarth = majEarth;
1995 0 : meta->gds.minEarth = minEarth;
1996 0 : if (majEarth == minEarth) {
1997 0 : meta->gds.f_sphere = 1;
1998 : }
1999 : } else {
2000 0 : meta->gds.f_sphere = 1;
2001 0 : meta->gds.majEarth = majEarth;
2002 0 : meta->gds.minEarth = majEarth;
2003 : }
2004 : }
2005 :
2006 66 : if( Grib_Data )
2007 : {
2008 : /* Allocate memory for the grid. */
2009 26 : if (meta->gds.numPts > *grib_DataLen) {
2010 26 : if( meta->gds.numPts > 100 * 1024 * 1024 )
2011 : {
2012 0 : vsi_l_offset curPos = VSIFTellL(fp);
2013 0 : VSIFSeekL(fp, 0, SEEK_END);
2014 0 : vsi_l_offset fileSize = VSIFTellL(fp);
2015 0 : VSIFSeekL(fp, curPos, SEEK_SET);
2016 : // allow a compression ratio of 1:1000
2017 0 : if( meta->gds.numPts / 1000 > (uInt4)fileSize )
2018 : {
2019 0 : errSprintf ("ERROR: File too short\n");
2020 0 : *grib_DataLen = 0;
2021 0 : *Grib_Data = nullptr;
2022 0 : return -2;
2023 : }
2024 : }
2025 :
2026 : #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
2027 : if( meta->gds.numPts > static_cast<size_t>(INT_MAX) / sizeof (double) )
2028 : {
2029 : errSprintf ("Memory allocation failed due to being bigger than 2 GB in fuzzing mode");
2030 : *grib_DataLen = 0;
2031 : *Grib_Data = nullptr;
2032 : return -2;
2033 : }
2034 : #endif
2035 :
2036 26 : *grib_DataLen = meta->gds.numPts;
2037 26 : *Grib_Data = (double *) realloc ((void *) (*Grib_Data),
2038 26 : (*grib_DataLen) * sizeof (double));
2039 26 : if (!(*Grib_Data))
2040 : {
2041 0 : *grib_DataLen = 0;
2042 0 : return -1;
2043 : }
2044 : }
2045 26 : grib_Data = *Grib_Data;
2046 : }
2047 :
2048 : /* Get the Bit Map Section. */
2049 66 : if (f_bms) {
2050 51 : if (ReadGrib1Sect3 (c_ipack + curLoc, gribLen, &curLoc, &bitmap,
2051 51 : meta->gds.numPts) != 0) {
2052 0 : free (bitmap);
2053 0 : preErrSprintf ("Inside ReadGrib1Record\n");
2054 0 : return -3;
2055 : }
2056 : }
2057 :
2058 : /* Figure out some basic stuff about the grid. */
2059 : /* Following is similar to metaparse.c : ParseElemName */
2060 66 : GRIB1_Table2LookUp (&(meta->pds1), &varName, &varComment, &varUnit,
2061 66 : &(meta->convert), meta->center, meta->subcenter);
2062 66 : meta->element = (char *) realloc ((void *) (meta->element),
2063 66 : (1 + strlen (varName)) * sizeof (char));
2064 66 : strcpy (meta->element, varName);
2065 66 : meta->unitName = (char *) realloc ((void *) (meta->unitName),
2066 66 : (1 + 2 + strlen (varUnit)) *
2067 : sizeof (char));
2068 66 : snprintf (meta->unitName,
2069 66 : (1 + 2 + strlen (varUnit)) *
2070 : sizeof (char),
2071 : "[%s]", varUnit);
2072 66 : meta->comment = (char *) realloc ((void *) (meta->comment),
2073 66 : (1 + strlen (varComment) +
2074 66 : strlen (varUnit)
2075 : + 2 + 1) * sizeof (char));
2076 66 : snprintf (meta->comment,
2077 66 : (1 + strlen (varComment) +
2078 66 : strlen (varUnit)
2079 : + 2 + 1) * sizeof (char),
2080 : "%s [%s]", varComment, varUnit);
2081 :
2082 66 : if (ComputeUnit (meta->convert, meta->unitName, f_unit, &unitM, &unitB,
2083 66 : unitName) == 0) {
2084 0 : unitLen = static_cast<int>(strlen (unitName));
2085 0 : meta->unitName = (char *) realloc ((void *) (meta->unitName),
2086 0 : 1 + unitLen * sizeof (char));
2087 0 : memcpy (meta->unitName, unitName, unitLen);
2088 0 : meta->unitName[unitLen] = '\0';
2089 : }
2090 :
2091 : /* Read the GRID. */
2092 66 : if (ReadGrib1Sect4 (c_ipack + curLoc, gribLen, &curLoc, DSF, grib_Data,
2093 66 : meta, f_bms, bitmap, unitM, unitB) != 0) {
2094 0 : free (bitmap);
2095 0 : preErrSprintf ("Inside ReadGrib1Record\n");
2096 0 : return -4;
2097 : }
2098 66 : if (f_bms) {
2099 51 : free (bitmap);
2100 : }
2101 :
2102 66 : GRIB1_Table3LookUp (&(meta->pds1), &(meta->shortFstLevel),
2103 : &(meta->longFstLevel));
2104 : /* printf ("%s .. %s\n", meta->shortFstLevel, meta->longFstLevel);*/
2105 :
2106 : /*
2107 : strftime (meta->refTime, 20, "%Y%m%d%H%M",
2108 : gmtime (&(meta->pds1.refTime)));
2109 : */
2110 66 : Clock_Print (meta->refTime, 20, meta->pds1.refTime, "%Y%m%d%H%M", 0);
2111 :
2112 : /*
2113 : strftime (meta->validTime, 20, "%Y%m%d%H%M",
2114 : gmtime (&(meta->pds1.validTime)));
2115 : */
2116 66 : Clock_Print (meta->validTime, 20, meta->pds1.validTime, "%Y%m%d%H%M", 0);
2117 :
2118 66 : double deltaTime = meta->pds1.validTime - meta->pds1.refTime;
2119 66 : if (deltaTime >= std::numeric_limits<sInt4>::max()) {
2120 0 : printf ("Clamped deltaTime. Was %lf\n", deltaTime);
2121 0 : deltaTime = std::numeric_limits<sInt4>::max();
2122 : }
2123 66 : if (deltaTime <= std::numeric_limits<sInt4>::min()) {
2124 0 : printf ("Clamped deltaTime. Was %lf\n", deltaTime);
2125 0 : deltaTime = std::numeric_limits<sInt4>::min();
2126 : }
2127 66 : meta->deltTime = static_cast<sInt4>(deltaTime);
2128 :
2129 : /* Read section 5. If it is "7777" == 926365495 we are done. */
2130 66 : if (curLoc == gribLen) {
2131 0 : printf ("Warning: either gribLen did not account for section 5, or "
2132 : "section 5 is missing\n");
2133 0 : return 0;
2134 : }
2135 66 : if (curLoc + 4 != gribLen) {
2136 0 : errSprintf ("Invalid number of bytes for the end of the message.\n");
2137 0 : return -5;
2138 : }
2139 66 : memcpy (&li_temp, c_ipack + curLoc, 4);
2140 66 : if (li_temp != 926365495L) {
2141 0 : errSprintf ("Did not find the end of the message.\n");
2142 0 : return -5;
2143 : }
2144 :
2145 66 : return 0;
2146 : }
2147 :
2148 : /*****************************************************************************
2149 : * main() --
2150 : *
2151 : * Arthur Taylor / MDL
2152 : *
2153 : * PURPOSE
2154 : * To test the capabilities of this module, and give an example as to how
2155 : * ReadGrib1Record expects to be called.
2156 : *
2157 : * ARGUMENTS
2158 : * argc = The number of arguments on the command line. (Input)
2159 : * argv = The arguments on the command line. (Input)
2160 : *
2161 : * FILES/DATABASES:
2162 : * A GRIB1 file.
2163 : *
2164 : * RETURNS: int
2165 : * 0 = OK
2166 : *
2167 : * HISTORY
2168 : * 4/2003 Arthur Taylor (MDL/RSIS): Created
2169 : * 6/2003 Matthew T. Kallio (matt@wunderground.com):
2170 : * "wmo" dimension increased to WMO_HEADER_LEN + 1 (for '\0' char)
2171 : *
2172 : * NOTES
2173 : *****************************************************************************
2174 : */
2175 : #ifdef DEBUG_DEGRIB1
2176 : int main (int argc, char **argv)
2177 : {
2178 : VSILFILE * grib_fp; /* The opened grib2 file for input. */
2179 : char *buff = nullptr;
2180 : uInt4 buffLen = 0;
2181 : sInt4 sect0[SECT0LEN_WORD] = { 0 }; /* Holds the current Section 0. */
2182 : uInt4 gribLen; /* Length of the current GRIB message. */
2183 : char *msg;
2184 : int version;
2185 : sChar f_unit = 0;
2186 : double *grib_Data;
2187 : uInt4 grib_DataLen;
2188 : grib_MetaData meta;
2189 :
2190 : double majEarth = 0.0; // -radEarth if < 6000 ignore, otherwise use this
2191 : // to override the radEarth in the GRIB1 or GRIB2
2192 : // message. Needed because NCEP uses 6371.2 but
2193 : // GRIB1 could only state 6367.47.
2194 : double minEarth = 0.0; // -minEarth if < 6000 ignore, otherwise use this
2195 : // to override the minEarth in the GRIB1 or GRIB2
2196 : // message.
2197 :
2198 :
2199 : IS_dataType is; /* Un-parsed meta data for this GRIB2 message. As
2200 : * well as some memory used by the unpacker. */
2201 :
2202 : if ((grib_fp = VSIFOpenL (argv[1], "rb")) == NULL) {
2203 : printf ("Problems opening %s for read\n", argv[1]);
2204 : return 1;
2205 : }
2206 : IS_Init (&is);
2207 : MetaInit (&meta);
2208 :
2209 : if (ReadSECT0 (grib_fp, &buff, &buffLen, -1, sect0, &gribLen, &version) < 0) {
2210 : VSIFCloseL(grib_fp);
2211 : free(buff);
2212 : msg = errSprintf (NULL);
2213 : printf ("%s\n", msg);
2214 : return -1;
2215 : }
2216 : free(buff);
2217 :
2218 : grib_DataLen = 0;
2219 : grib_Data = NULL;
2220 : if (version == 1) {
2221 : meta.GribVersion = version;
2222 : ReadGrib1Record (grib_fp, f_unit, &grib_Data, &grib_DataLen, &meta,
2223 : &is, sect0, gribLen, majEarth, minEarth);
2224 : }
2225 :
2226 : MetaFree (&meta);
2227 : IS_Free (&is);
2228 : free (grib_Data);
2229 : VSIFCloseL(grib_fp);
2230 : return 0;
2231 : }
2232 : #endif
|
