Line data Source code
1 : /*****************************************************************************
2 : * inventory.c
3 : *
4 : * DESCRIPTION
5 : * This file contains the code needed to do a quick inventory of the GRIB2
6 : * file. The intent is to enable one to figure out which message in a GRIB
7 : * file one is after without needing to call the FORTRAN library.
8 : *
9 : * HISTORY
10 : * 9/2002 Arthur Taylor (MDL / RSIS): Created.
11 : * 12/2002 Tim Kempisty, Ana Canizares, Tim Boyer, & Marc Saccucci
12 : * (TK,AC,TB,&MS): Code Review 1.
13 : *
14 : * NOTES
15 : *****************************************************************************
16 : */
17 : #include <stdio.h>
18 : #include <string.h>
19 : #include <stdlib.h>
20 : #include <math.h>
21 :
22 : #include <algorithm>
23 : #include <cmath>
24 : #include <limits>
25 :
26 : #include "clock.h"
27 : #include "tendian.h"
28 : #include "degrib2.h"
29 : #include "degrib1.h"
30 : #if 0
31 : /* tdlpack is no longer supported by GDAL */
32 : #include "tdlpack.h"
33 : #endif
34 : #include "myerror.h"
35 : #include "myutil.h"
36 : #include "myassert.h"
37 : #include "inventory.h"
38 : #include "metaname.h"
39 :
40 : #include "cpl_port.h"
41 : #include "cpl_error.h"
42 :
43 : #define SECT0LEN_BYTE 16
44 :
45 0 : static sInt4 DoubleToSInt4Clamp(double val) {
46 0 : if (val >= INT_MAX) return INT_MAX;
47 0 : if (val <= INT_MIN) return INT_MIN;
48 0 : if (std::isnan(val)) return 0;
49 0 : return (sInt4)val;
50 : }
51 :
52 : typedef union {
53 : sInt4 li;
54 : char buffer[4];
55 : } wordType;
56 :
57 : /*****************************************************************************
58 : * GRIB2InventoryFree() -- Review 12/2002
59 : *
60 : * Arthur Taylor / MDL
61 : *
62 : * PURPOSE
63 : * Free's any memory that was allocated for the inventory of a single grib
64 : * message
65 : *
66 : * ARGUMENTS
67 : * inv = Pointer to the inventory of a single grib message. (Input/Output)
68 : *
69 : * FILES/DATABASES: None
70 : *
71 : * RETURNS: void
72 : *
73 : * HISTORY
74 : * 9/2002 Arthur Taylor (MDL/RSIS): Created.
75 : * 12/2002 (TK,AC,TB,&MS): Code Review.
76 : * 7/2003 AAT: memwatch detected unfreed inv->unitName
77 : *
78 : * NOTES
79 : *****************************************************************************
80 : */
81 422 : void GRIB2InventoryFree (inventoryType *inv)
82 : {
83 422 : free (inv->element);
84 422 : inv->element = nullptr;
85 422 : free (inv->comment);
86 422 : inv->comment = nullptr;
87 422 : free (inv->unitName);
88 422 : inv->unitName = nullptr;
89 422 : free (inv->shortFstLevel);
90 422 : inv->shortFstLevel = nullptr;
91 422 : free (inv->longFstLevel);
92 422 : inv->longFstLevel = nullptr;
93 422 : }
94 :
95 : /*****************************************************************************
96 : * GRIB2InventoryPrint() -- Review 12/2002
97 : *
98 : * Arthur Taylor / MDL
99 : *
100 : * PURPOSE
101 : * Prints to standard out, an inventory of the file, assuming one has an
102 : * array of inventories of single GRIB messages.
103 : *
104 : * ARGUMENTS
105 : * Inv = Pointer to an Array of inventories to print. (Input)
106 : * LenInv = Length of the Array Inv (Input)
107 : *
108 : * FILES/DATABASES: None
109 : *
110 : * RETURNS: void
111 : *
112 : * HISTORY
113 : * 9/2002 Arthur Taylor (MDL/RSIS): Created.
114 : * 12/2002 (TK,AC,TB,&MS): Code Review.
115 : * 1/2004 AAT: Added short form of First level to print out.
116 : * 3/2004 AAT: Switched from "#, Byte, ..." to "MsgNum, Byte, ..."
117 : *
118 : * NOTES
119 : *****************************************************************************
120 : */
121 0 : void GRIB2InventoryPrint (inventoryType *Inv, uInt4 LenInv)
122 : {
123 : uInt4 i; /* Counter of which inventory we are printing. */
124 : double delta; /* Difference between valid and reference time. */
125 : char refTime[25]; /* Used to store the formatted reference time. */
126 : char validTime[25]; /* Used to store the formatted valid time. */
127 :
128 0 : printf ("MsgNum, Byte, GRIB-Version, elem, level, reference(UTC),"
129 : " valid(UTC), Proj(hr)\n");
130 0 : fflush (stdout);
131 0 : for (i = 0; i < LenInv; i++) {
132 : /* strftime (refTime, 25, "%m/%d/%Y %H:%M", gmtime (&(Inv[i].refTime)));*/
133 0 : Clock_Print (refTime, 25, Inv[i].refTime, "%m/%d/%Y %H:%M", 0);
134 : /* strftime (validTime, 25, "%m/%d/%Y %H:%M",
135 : gmtime (&(Inv[i].validTime)));*/
136 0 : Clock_Print (validTime, 25, Inv[i].validTime, "%m/%d/%Y %H:%M", 0);
137 0 : delta = (Inv[i].validTime - Inv[i].refTime) / 3600.;
138 0 : delta = myRound (delta, 2);
139 0 : if (Inv[i].comment == nullptr) {
140 0 : printf ("%u.%u, " CPL_FRMT_GUIB ", %d, %s, %s, %s, %s, %.2f\n",
141 0 : Inv[i].msgNum, Inv[i].subgNum, Inv[i].start,
142 0 : Inv[i].GribVersion, Inv[i].element, Inv[i].shortFstLevel,
143 : refTime, validTime, delta);
144 0 : fflush (stdout);
145 : } else {
146 0 : printf ("%u.%u, " CPL_FRMT_GUIB ", %d, %s=\"%s\", %s, %s, %s, %.2f\n",
147 0 : Inv[i].msgNum, Inv[i].subgNum, Inv[i].start,
148 0 : Inv[i].GribVersion, Inv[i].element, Inv[i].comment,
149 0 : Inv[i].shortFstLevel, refTime, validTime, delta);
150 0 : fflush (stdout);
151 : }
152 : }
153 0 : }
154 :
155 : /*****************************************************************************
156 : * InventoryParseTime() -- Review 12/2002
157 : *
158 : * Arthur Taylor / MDL
159 : *
160 : * PURPOSE
161 : * To parse the time data from a grib2 char array to a time_t in UTC
162 : * seconds from the epoch. This is very similar to metaparse.c:ParseTime
163 : * except using char * instead of sInt4
164 : *
165 : * ARGUMENTS
166 : * is = The char array to read the time info from. (Input)
167 : * AnsTime = The time_t value to fill with the resulting time. (Output)
168 : *
169 : * FILES/DATABASES: None
170 : *
171 : * RETURNS: void
172 : *
173 : * HISTORY
174 : * 11/2002 Arthur Taylor (MDL/RSIS): Created.
175 : * 12/2002 (TK,AC,TB,&MS): Code Review.
176 : *
177 : * NOTES
178 : * 1) Couldn't use the default time_zone variable (concern over portability
179 : * issues), so we print the hours, and compare them to the hours we had
180 : * intended. Then subtract the difference from the AnsTime.
181 : * 2) Similar to metaparse.c:ParseTime except using char * instead of sInt4
182 : *****************************************************************************
183 : */
184 392 : static int InventoryParseTime (char *is, double *AnsTime)
185 : {
186 : /* struct tm time; *//* A temporary variable to put the time info into. */
187 : /* char buffer[10]; *//* Used when printing the AnsTime's Hr. */
188 : /* int timeZone; *//* The adjustment in Hr needed to get the right UTC * time. */
189 : short int si_temp; /* Temporarily stores the year as a short int to fix
190 : * possible endian problems. */
191 :
192 : /* memset (&time, 0, sizeof (struct tm));*/
193 392 : MEMCPY_BIG (&si_temp, is + 0, sizeof (short int));
194 392 : if ((si_temp < 1900) || (si_temp > 2100)) {
195 0 : return -1;
196 : }
197 392 : if ((is[2] > 12) || (is[3] == 0) || (is[3] > 31) || (is[4] > 24) ||
198 392 : (is[5] > 60) || (is[6] > 61)) {
199 0 : return -1;
200 : }
201 392 : Clock_ScanDate (AnsTime, si_temp, is[2], is[3]);
202 392 : *AnsTime += is[4] * 3600. + is[5] * 60. + is[6];
203 : /*
204 : time.tm_year = si_temp - 1900;
205 : time.tm_mon = is[2] - 1;
206 : time.tm_mday = is[3];
207 : time.tm_hour = is[4];
208 : time.tm_min = is[5];
209 : time.tm_sec = is[6];
210 : *AnsTime = mktime (&time) - (Clock_GetTimeZone () * 3600);
211 : */
212 : /* Cheap method of getting global time_zone variable. */
213 : /*
214 : strftime (buffer, 10, "%H", gmtime (AnsTime));
215 : timeZone = atoi (buffer) - is[4];
216 : if (timeZone < 0) {
217 : timeZone += 24;
218 : }
219 : *AnsTime = *AnsTime - (timeZone * 3600);
220 : */
221 392 : return 0;
222 : }
223 :
224 : /*****************************************************************************
225 : * GRIB2SectToBuffer() -- Review 12/2002
226 : *
227 : * Arthur Taylor / MDL
228 : *
229 : * PURPOSE
230 : * To read in a GRIB2 section into a buffer. Reallocates space for the
231 : * section if buffLen < secLen. Reads in secLen and checks that the section
232 : * is valid, and the file is large enough to hold the entire section.
233 : *
234 : * ARGUMENTS
235 : * fp = Opened file pointing to the section in question. (Input/Output)
236 : * gribLen = The total length of the grib message. (Input)
237 : * sect = Which section we think we are reading.
238 : * If it is -1, then set it to the section the file says we are
239 : * reading (useful for optional sect 2)) (Input/Output).
240 : * secLen = The length of this section (Output)
241 : * buffLen = Allocated length of buff (Input/Output)
242 : * buff = Stores the section (Output)
243 : *
244 : * FILES/DATABASES:
245 : * An already opened GRIB2 file pointer, already at section in question.
246 : *
247 : * RETURNS: int (could use errSprintf())
248 : * 0 = Ok.
249 : * -1 = Ran out of file.
250 : * -2 = Section was miss-labeled.
251 : *
252 : * HISTORY
253 : * 11/2002 Arthur Taylor (MDL/RSIS): Created.
254 : * 12/2002 (TK,AC,TB,&MS): Code Review.
255 : * 8/2003 AAT: Removed dependence on curTot
256 : *
257 : * NOTES
258 : * May want to put this in degrib2.c
259 : *****************************************************************************
260 : */
261 704 : static int GRIB2SectToBuffer (VSILFILE *fp,
262 : uInt4 gribLen,
263 : sChar *sect,
264 : uInt4 *secLen, uInt4 *buffLen, char **buff)
265 : {
266 704 : char *buffer = *buff; /* Local ptr to buff to reduce ptr confusion. */
267 :
268 704 : if (FREAD_BIG (secLen, sizeof (sInt4), 1, fp) != 1) {
269 1 : if (*sect != -1) {
270 1 : errSprintf ("ERROR: Ran out of file in Section %d\n", *sect);
271 : } else {
272 0 : errSprintf ("ERROR: Ran out of file in GRIB2SectToBuffer\n");
273 : }
274 1 : return -1;
275 : }
276 703 : if( *secLen <= sizeof(sInt4) || *secLen > gribLen )
277 : {
278 0 : errSprintf ("ERROR: Wrong secLen in GRIB2SectToBuffer\n");
279 0 : return -1;
280 : }
281 703 : if (*buffLen < *secLen) {
282 590 : if( *secLen > 100 * 1024 * 1024 )
283 : {
284 0 : vsi_l_offset curPos = VSIFTellL(fp);
285 0 : VSIFSeekL(fp, 0, SEEK_END);
286 0 : vsi_l_offset fileSize = VSIFTellL(fp);
287 0 : VSIFSeekL(fp, curPos, SEEK_SET);
288 0 : if( *secLen > fileSize )
289 : {
290 0 : errSprintf ("ERROR: File too short\n");
291 0 : return -1;
292 : }
293 : }
294 590 : char* buffnew = (char *) realloc ((void *) *buff, *secLen * sizeof (char));
295 590 : if( buffnew == nullptr )
296 : {
297 0 : errSprintf ("ERROR: Ran out of memory in GRIB2SectToBuffer\n");
298 0 : return -1;
299 : }
300 590 : *buffLen = *secLen;
301 590 : *buff = buffnew;
302 590 : buffer = *buff;
303 : }
304 :
305 703 : if (VSIFReadL(buffer, sizeof (char), *secLen - sizeof (sInt4), fp) !=
306 703 : *secLen - sizeof (sInt4)) {
307 1 : if (*sect != -1) {
308 1 : errSprintf ("ERROR: Ran out of file in Section %d\n", *sect);
309 : } else {
310 0 : errSprintf ("ERROR: Ran out of file in GRIB2SectToBuffer\n");
311 : }
312 1 : return -1;
313 : }
314 702 : if (*sect == -1) {
315 0 : *sect = buffer[0];
316 702 : } else if (buffer[0] != *sect) {
317 0 : errSprintf ("ERROR: Section %d mislabeled\n", *sect);
318 0 : return -2;
319 : }
320 702 : return 0;
321 : }
322 :
323 : /*****************************************************************************
324 : * GRIB2SectJump() --
325 : *
326 : * Arthur Taylor / MDL
327 : *
328 : * PURPOSE
329 : * To jump past a GRIB2 section. Reads in secLen and checks that the
330 : * section is valid.
331 : *
332 : * ARGUMENTS
333 : * fp = Opened file pointing to the section in question. (Input/Output)
334 : * gribLen = The total length of the grib message. (Input)
335 : * sect = Which section we think we are reading.
336 : * If it is -1, then set it to the section the file says we are
337 : * reading (useful for optional sect 2)) (Input/Output).
338 : * secLen = The length of this section (Output)
339 : *
340 : * FILES/DATABASES:
341 : * An already opened GRIB2 file pointer, already at section in question.
342 : *
343 : * RETURNS: int (could use errSprintf())
344 : * 0 = Ok.
345 : * -1 = Ran out of file.
346 : * -2 = Section was miss-labeled.
347 : *
348 : * HISTORY
349 : * 3/2003 Arthur Taylor (MDL/RSIS): Created.
350 : * 8/2003 AAT: Removed dependence on curTot, which was used to compute if
351 : * the file should be large enough for the fseek, but didn't check
352 : * if it actually was.
353 : *
354 : * NOTES
355 : * May want to put this in degrib2.c
356 : *****************************************************************************
357 : */
358 1707 : static int GRIB2SectJump (VSILFILE *fp,
359 : CPL_UNUSED sInt4 gribLen, sChar *sect, uInt4 *secLen)
360 : {
361 : char sectNum; /* Validates that we are on the correct section. */
362 :
363 1707 : if (FREAD_BIG (secLen, sizeof (sInt4), 1, fp) != 1) {
364 1 : if (*sect != -1) {
365 1 : errSprintf ("ERROR: Ran out of file in Section %d\n", *sect);
366 : } else {
367 0 : errSprintf ("ERROR: Ran out of file in GRIB2SectSkip\n");
368 : }
369 1 : return -1;
370 : }
371 1706 : if (*secLen < 5 || VSIFReadL (§Num, sizeof (char), 1, fp) != 1) {
372 0 : if (*sect != -1) {
373 0 : errSprintf ("ERROR: Ran out of file in Section %d\n", *sect);
374 : } else {
375 0 : errSprintf ("ERROR: Ran out of file in GRIB2SectSkip\n");
376 : }
377 0 : return -1;
378 : }
379 1706 : if (*sect == -1) {
380 352 : *sect = sectNum;
381 1354 : } else if (sectNum != *sect) {
382 0 : errSprintf ("ERROR: Section %d mislabeled\n", *sect);
383 0 : return -2;
384 : }
385 : /* Since fseek does not give an error if we jump outside the file, we test
386 : * it by using fgetc / ungetc. */
387 1706 : VSIFSeekL(fp, *secLen - 5, SEEK_CUR);
388 1706 : if (VSIFReadL (§Num, sizeof (char), 1, fp) != 1) {
389 5 : errSprintf ("ERROR: Ran out of file in Section %d\n", *sect);
390 5 : return -1;
391 : } else {
392 1701 : VSIFSeekL(fp, VSIFTellL(fp) - sizeof(char), SEEK_SET);
393 : }
394 1701 : return 0;
395 : }
396 :
397 : /*****************************************************************************
398 : * GRIB2Inventory2to7() --
399 : *
400 : * Arthur Taylor / MDL
401 : *
402 : * PURPOSE
403 : * Inventories sections 3 to 7, filling out the inv record with the data in
404 : * section 4. (Note: No Call to FORTRAN routines here).
405 : *
406 : * ARGUMENTS
407 : * sectNum = Which section we are currently reading. (Input)
408 : * fp = An opened file pointer to the file to the inventory of (In/Out)
409 : * gribLen = The total length of the grib message. (Input)
410 : * buffLen = length of buffer. (Input)
411 : * buffer = Holds a given section. (Input)
412 : * inv = The current inventory record to fill out. (Output)
413 : * prodType = The GRIB2 type of product: 0 is meteo product, 1 is hydro,
414 : * 2 is land, 3 is space, 10 is oceanographic. (Input)
415 : * center = Who produced it (Input)
416 : * subcenter = A sub group of center that actually produced it (Input)
417 : *
418 : * FILES/DATABASES:
419 : *
420 : * RETURNS: int (could use errSprintf())
421 : * 0 = "Ok"
422 : * -5 = Problems Reading in section 2 or 3
423 : * -6 = Problems Reading in section 3
424 : * -7 = Problems Reading in section 4
425 : * -8 = Problems Parsing section 4.
426 : * -9 = Problems Reading in section 5
427 : * -10 = Problems Reading in section 6
428 : * -11 = Problems Reading in section 7
429 : *
430 : * HISTORY
431 : * 3/2003 Arthur Taylor (MDL/RSIS): Created.
432 : * 4/2003 AAT: Modified to not have prodType, cat, subcat, templat in
433 : * inventoryType structure.
434 : * 8/2003 AAT: curTot no longer serves a purpose.
435 : * 1/2004 AAT: Added center/subcenter.
436 : *
437 : * NOTES
438 : *****************************************************************************
439 : */
440 354 : static int GRIB2Inventory2to7 (sChar sectNum, VSILFILE *fp, sInt4 gribLen,
441 : uInt4 *buffLen, char **buffer,
442 : inventoryType *inv, uChar prodType,
443 : unsigned short int center,
444 : unsigned short int subcenter,
445 : uChar mstrVersion)
446 : {
447 : uInt4 secLen; /* The length of the current section. */
448 : sInt4 foreTime; /* forecast time (NDFD treats as "projection") */
449 : uChar foreTimeUnit; /* The time unit of the "forecast time". */
450 : /* char *element; *//* Holds the name of the current variable. */
451 : /* char *comment; *//* Holds more comments about current variable. */
452 : /* char *unitName; *//* Holds the name of the unit [K] [%] .. etc */
453 : int convert; /* Enum type of unit conversions (metaname.c),
454 : * Conversion method for this variable's unit. */
455 : uChar cat; /* General category of Meteo Product. */
456 : unsigned short int templat; /* The section 4 template number. */
457 : uChar subcat; /* Specific subcategory of Product. */
458 : uChar genProcess; /* What type of generate process (Analysis,
459 : Forecast, Probability Forecast, etc). */
460 : uChar fstSurfType; /* Type of the first fixed surface. */
461 : double fstSurfValue; /* Value of first fixed surface. */
462 : sInt4 value; /* The scaled value from GRIB2 file. */
463 : sChar factor; /* The scaled factor from GRIB2 file */
464 : sChar scale; /* Surface scale as opposed to probability factor. */
465 : uChar sndSurfType; /* Type of the second fixed surface. */
466 : double sndSurfValue; /* Value of second fixed surface. */
467 : sChar f_sndValue; /* flag if SndValue is valid. */
468 : sChar f_fstValue; /* flag if FstValue is valid. */
469 : uChar timeRangeUnit;
470 354 : uChar statProcessID = 255;
471 : sInt4 lenTime; /* Used by parseTime to tell difference between 8hr
472 : * average and 1hr average ozone. */
473 : uChar genID; /* The Generating process ID (used for GFS MOS) */
474 : uChar probType; /* The probability type */
475 : double lowerProb; /* The lower limit on probability forecast if
476 : * template 4.5 or 4.9 */
477 : double upperProb; /* The upper limit on probability forecast if
478 : * template 4.5 or 4.9 */
479 : uChar timeIncrType;
480 354 : sChar percentile = 0;
481 :
482 354 : if ((sectNum == 2) || (sectNum == 3)) {
483 : /* Jump past section (2 or 3). */
484 352 : sectNum = -1;
485 352 : if (GRIB2SectJump (fp, gribLen, §Num, &secLen) != 0) {
486 0 : errSprintf ("ERROR: Problems Jumping past section 2 || 3\n");
487 0 : return -6;
488 : }
489 352 : if ((sectNum != 2) && (sectNum != 3)) {
490 0 : errSprintf ("ERROR: Section 2 or 3 mislabeled\n");
491 0 : return -5;
492 352 : } else if (sectNum == 2) {
493 : /* Jump past section 3. */
494 312 : sectNum = 3;
495 312 : if (GRIB2SectJump (fp, gribLen, §Num, &secLen) != 0) {
496 4 : errSprintf ("ERROR: Problems Jumping past section 3\n");
497 4 : return -6;
498 : }
499 : }
500 : }
501 : /* Read section 4 into buffer. */
502 350 : sectNum = 4;
503 350 : if (GRIB2SectToBuffer (fp, gribLen, §Num, &secLen, buffLen,
504 350 : buffer) != 0) {
505 0 : errSprintf ("ERROR: Problems with section 4\n");
506 0 : return -7;
507 : }
508 : /*
509 : enum { GS4_ANALYSIS, GS4_ENSEMBLE, GS4_DERIVED, GS4_PROBABIL_PNT = 5,
510 : GS4_PERCENT_PNT, GS4_STATISTIC = 8, GS4_PROBABIL_TIME = 9,
511 : GS4_PERCENT_TIME = 10, GS4_RADAR = 20, GS4_SATELLITE = 30
512 : };
513 : */
514 350 : if( secLen < 11 )
515 0 : return -8;
516 :
517 : /* Parse the interesting data out of sect 4. */
518 350 : MEMCPY_BIG (&templat, *buffer + 8 - 5, sizeof (short int));
519 350 : if ((templat != GS4_ANALYSIS) && (templat != GS4_ENSEMBLE)
520 63 : && (templat != GS4_DERIVED)
521 63 : && (templat != GS4_PROBABIL_PNT) && (templat != GS4_PERCENT_PNT)
522 63 : && (templat != GS4_ERROR)
523 63 : && (templat != GS4_STATISTIC)
524 25 : && (templat != GS4_PROBABIL_TIME) && (templat != GS4_PERCENT_TIME)
525 25 : && (templat != GS4_ENSEMBLE_STAT)
526 25 : && (templat != GS4_STATISTIC_SPATIAL_AREA)
527 24 : && (templat != GS4_RADAR) && (templat != GS4_SATELLITE)
528 24 : && (templat != GS4_SATELLITE_SYNTHETIC)
529 17 : && (templat != GS4_DERIVED_INTERVAL)
530 15 : && (templat != GS4_ANALYSIS_CHEMICAL)
531 6 : && (templat != GS4_OPTICAL_PROPERTIES_AEROSOL) ) {
532 5 : errSprintf ("This was only designed for templates 0, 1, 2, 5, 6, 7, 8, 9, "
533 : "10, 11, 12, 15, 20, 30, 32, 40, 48. Template found = %d\n", templat);
534 :
535 5 : inv->validTime = 0;
536 5 : inv->foreSec = 0;
537 5 : reallocSprintf (&(inv->element), "unknown");
538 5 : reallocSprintf (&(inv->comment), "unknown");
539 5 : reallocSprintf (&(inv->unitName), "unknown");
540 5 : reallocSprintf (&(inv->shortFstLevel), "unknown");
541 5 : reallocSprintf (&(inv->longFstLevel), "unknown");
542 :
543 : /* Jump past section 5. */
544 5 : sectNum = 5;
545 5 : if (GRIB2SectJump (fp, gribLen, §Num, &secLen) != 0) {
546 0 : errSprintf ("ERROR: Problems Jumping past section 5\n");
547 0 : return -9;
548 : }
549 : /* Jump past section 6. */
550 5 : sectNum = 6;
551 5 : if (GRIB2SectJump (fp, gribLen, §Num, &secLen) != 0) {
552 0 : errSprintf ("ERROR: Problems Jumping past section 6\n");
553 0 : return -10;
554 : }
555 : /* Jump past section 7. */
556 5 : sectNum = 7;
557 5 : if (GRIB2SectJump (fp, gribLen, §Num, &secLen) != 0) {
558 0 : errSprintf ("ERROR: Problems Jumping past section 7\n");
559 0 : return -11;
560 : }
561 5 : return 1;
562 : }
563 :
564 345 : unsigned nOffset = 0;
565 345 : if( templat == GS4_ANALYSIS_CHEMICAL ) {
566 9 : nOffset = 16 - 14;
567 : }
568 336 : else if( templat == GS4_OPTICAL_PROPERTIES_AEROSOL )
569 : {
570 1 : nOffset = 38 - 14;
571 : }
572 :
573 345 : if( secLen < nOffset + 19 - 5 + 4 )
574 1 : return -8;
575 :
576 344 : cat = (*buffer)[10 - 5];
577 344 : subcat = (*buffer)[11 - 5];
578 344 : genProcess = (*buffer)[nOffset + 12 - 5];
579 344 : genID = 0;
580 344 : probType = 0;
581 344 : lowerProb = 0;
582 344 : upperProb = 0;
583 344 : if ((templat == GS4_RADAR) || (templat == GS4_SATELLITE) ||
584 344 : (templat == 254)) {
585 0 : inv->foreSec = 0;
586 0 : inv->validTime = inv->refTime;
587 0 : timeIncrType = 255;
588 0 : timeRangeUnit = 255;
589 0 : lenTime = 0;
590 : } else {
591 344 : genID = (*buffer)[nOffset + 14 - 5];
592 : /* Compute forecast time. */
593 344 : foreTimeUnit = (*buffer)[nOffset + 18 - 5];
594 344 : MEMCPY_BIG (&foreTime, *buffer + nOffset + 19 - 5, sizeof (sInt4));
595 344 : if (ParseSect4Time2sec (inv->refTime, foreTime, foreTimeUnit, &(inv->foreSec)) != 0) {
596 0 : errSprintf ("unable to convert TimeUnit: %d \n", foreTimeUnit);
597 0 : return -8;
598 : }
599 : /* Compute valid time. */
600 344 : inv->validTime = inv->refTime + inv->foreSec;
601 344 : timeIncrType = 255;
602 344 : timeRangeUnit = 1;
603 344 : lenTime = static_cast<sInt4>(
604 344 : std::max(static_cast<double>(std::numeric_limits<sInt4>::min()),
605 1032 : std::min(static_cast<double>(std::numeric_limits<sInt4>::max()),
606 688 : inv->foreSec / 3600.0)));
607 344 : switch (templat) {
608 0 : case GS4_PROBABIL_PNT: /* 4.5 */
609 0 : if( secLen < 44 - 5 + 4)
610 0 : return -8;
611 0 : probType = (*buffer)[37 - 5];
612 0 : factor = sbit_2Comp_oneByte((sChar) (*buffer)[38 - 5]);
613 0 : MEMCPY_BIG (&value, *buffer + 39 - 5, sizeof (sInt4));
614 0 : value = sbit_2Comp_fourByte(value);
615 0 : lowerProb = value * pow (10.0, -1 * factor);
616 0 : factor = sbit_2Comp_oneByte((sChar) (*buffer)[43 - 5]);
617 0 : MEMCPY_BIG (&value, *buffer + 44 - 5, sizeof (sInt4));
618 0 : value = sbit_2Comp_fourByte(value);
619 0 : upperProb = value * pow (10.0, -1 * factor);
620 0 : break;
621 :
622 0 : case GS4_PERCENT_PNT: /* 4.6 */
623 0 : if( secLen < 35 - 5 + 1)
624 0 : return -8;
625 0 : percentile = (*buffer)[35 - 5];
626 0 : break;
627 :
628 2 : case GS4_DERIVED_INTERVAL: /* 4.12 */
629 2 : if( secLen < 52 - 5 + 4)
630 0 : return -8;
631 2 : if (InventoryParseTime (*buffer + 37 - 5, &(inv->validTime)) != 0) {
632 0 : printf ("Warning: Investigate Template 4.12 bytes 37-43\n");
633 0 : inv->validTime = inv->refTime + inv->foreSec;
634 : }
635 2 : timeIncrType = (*buffer)[50 - 5];
636 2 : timeRangeUnit = (*buffer)[51 - 5];
637 2 : MEMCPY_BIG (&lenTime, *buffer + 52 - 5, sizeof (sInt4));
638 : /* If lenTime == missing (2^32 -1) we might do something, but not with 255.*/
639 : /*
640 : if (lenTime == 255) {
641 : lenTime = (inv->validTime -
642 : (inv->refTime + inv->foreSec)) / 3600;
643 : }
644 : */
645 2 : break;
646 :
647 0 : case GS4_PERCENT_TIME: /* 4.10 */
648 0 : if( secLen < 51 - 5 + 4)
649 0 : return -8;
650 0 : percentile = (*buffer)[35 - 5];
651 0 : if (InventoryParseTime (*buffer + 36 - 5, &(inv->validTime)) != 0) {
652 0 : printf ("Warning: Investigate Template 4.10 bytes 36-42\n");
653 0 : inv->validTime = inv->refTime + inv->foreSec;
654 : }
655 0 : timeIncrType = (*buffer)[49 - 5];
656 0 : timeRangeUnit = (*buffer)[50 - 5];
657 0 : MEMCPY_BIG (&lenTime, *buffer + 51 - 5, sizeof (sInt4));
658 : /* If lenTime == missing (2^32 -1) we might do something, but not with 255.*/
659 : /*
660 : if (lenTime == 255) {
661 : lenTime = (inv->validTime -
662 : (inv->refTime + inv->foreSec)) / 3600;
663 : }
664 : */
665 0 : break;
666 38 : case GS4_STATISTIC: /* 4.8 */
667 38 : if( secLen < 50 - 5 + 4)
668 0 : return -8;
669 38 : if (InventoryParseTime (*buffer + 35 - 5, &(inv->validTime)) != 0) {
670 0 : printf ("Warning: Investigate Template 4.8 bytes 35-41\n");
671 0 : inv->validTime = inv->refTime + inv->foreSec;
672 : }
673 38 : statProcessID = (*buffer)[47 -5];
674 38 : timeIncrType = (*buffer)[48 - 5];
675 38 : timeRangeUnit = (*buffer)[49 - 5];
676 38 : MEMCPY_BIG (&lenTime, *buffer + 50 - 5, sizeof (sInt4));
677 : /* If lenTime == missing (2^32 -1) we might do something, but not with 255.*/
678 : /*
679 : if (lenTime == 255) {
680 : lenTime = (inv->validTime -
681 : (inv->refTime + inv->foreSec)) / 3600;
682 : }
683 : */
684 38 : break;
685 0 : case GS4_ENSEMBLE_STAT: /* 4.11 */
686 0 : if( secLen < 53 - 5 + 4)
687 0 : return -8;
688 0 : if (InventoryParseTime (*buffer + 38 - 5, &(inv->validTime)) != 0) {
689 0 : printf ("Warning: Investigate Template 4.11 bytes 38-44\n");
690 0 : inv->validTime = inv->refTime + inv->foreSec;
691 : }
692 0 : timeIncrType = (*buffer)[51 - 5];
693 0 : timeRangeUnit = (*buffer)[52 - 5];
694 0 : MEMCPY_BIG (&lenTime, *buffer + 53 - 5, sizeof (sInt4));
695 : /* If lenTime == missing (2^32 -1) we might do something, but not with 255.*/
696 : /*
697 : if (lenTime == 255) {
698 : lenTime = (inv->validTime -
699 : (inv->refTime + inv->foreSec)) / 3600;
700 : }
701 : */
702 0 : break;
703 0 : case GS4_PROBABIL_TIME: /* 4.9 */
704 0 : if( secLen < 63 - 5 + 4)
705 0 : return -8;
706 0 : probType = (*buffer)[37 - 5];
707 0 : factor = sbit_2Comp_oneByte((sChar) (*buffer)[38 - 5]);
708 0 : MEMCPY_BIG (&value, *buffer + 39 - 5, sizeof (sInt4));
709 0 : value = sbit_2Comp_fourByte(value);
710 0 : lowerProb = value * pow (10.0, -1 * factor);
711 0 : factor = sbit_2Comp_oneByte((sChar) (*buffer)[43 - 5]);
712 0 : MEMCPY_BIG (&value, *buffer + 44 - 5, sizeof (sInt4));
713 0 : value = sbit_2Comp_fourByte(value);
714 0 : upperProb = value * pow (10.0, -1 * factor);
715 :
716 0 : if (InventoryParseTime (*buffer + 48 - 5, &(inv->validTime)) != 0) {
717 0 : printf ("Warning: Investigate Template 4.9 bytes 48-54\n");
718 0 : inv->validTime = inv->refTime + inv->foreSec;
719 : }
720 0 : timeIncrType = (*buffer)[61 - 5];
721 0 : timeRangeUnit = (*buffer)[62 - 5];
722 0 : MEMCPY_BIG (&lenTime, *buffer + 63 - 5, sizeof (sInt4));
723 : /* If lenTime == missing (2^32 -1) we might do something, but not with 255.*/
724 : /*
725 : if (lenTime == 255) {
726 : lenTime = (inv->validTime -
727 : (inv->refTime + inv->foreSec)) / 3600;
728 : }
729 : */
730 0 : break;
731 : }
732 : }
733 :
734 344 : uChar derivedFcst = (uChar)-1; // not determined
735 344 : switch (templat) {
736 2 : case GS4_DERIVED:
737 : case GS4_DERIVED_CLUSTER_RECTANGULAR_AREA:
738 : case GS4_DERIVED_CLUSTER_CIRCULAR_AREA:
739 : case GS4_DERIVED_INTERVAL:
740 : case GS4_DERIVED_INTERVAL_CLUSTER_RECTANGULAR_AREA:
741 : case GS4_DERIVED_INTERVAL_CLUSTER_CIRCULAR_AREA:
742 2 : if( secLen >= 35 ) {
743 2 : derivedFcst = (uChar) (*buffer)[35 - 5];
744 : }
745 2 : break;
746 342 : default:
747 342 : break;
748 : }
749 :
750 344 : if (timeRangeUnit == 255) {
751 0 : timeRangeUnit = 1;
752 0 : lenTime = DoubleToSInt4Clamp(
753 0 : (inv->validTime - inv->foreSec - inv->refTime) / 3600.0);
754 : }
755 : /* myAssert (timeRangeUnit == 1);*/
756 : /* Try to convert lenTime to hourly. */
757 344 : if (timeRangeUnit == 0) {
758 0 : lenTime = (sInt4) (lenTime / 60.);
759 0 : timeRangeUnit = 1;
760 344 : } else if (timeRangeUnit == 1) {
761 0 : } else if (timeRangeUnit == 2) {
762 0 : if( lenTime < INT_MIN / 24 || lenTime > INT_MAX / 24 )
763 0 : return -8;
764 0 : lenTime = lenTime * 24;
765 0 : timeRangeUnit = 1;
766 0 : } else if (timeRangeUnit == 10) {
767 0 : if( lenTime < INT_MIN / 3 || lenTime > INT_MAX / 3 )
768 0 : return -8;
769 0 : lenTime = lenTime * 3;
770 0 : timeRangeUnit = 1;
771 0 : } else if (timeRangeUnit == 11) {
772 0 : if( lenTime < INT_MIN / 6 || lenTime > INT_MAX / 6 )
773 0 : return -8;
774 0 : lenTime = lenTime * 6;
775 0 : timeRangeUnit = 1;
776 0 : } else if (timeRangeUnit == 12) {
777 0 : if( lenTime < INT_MIN / 12 || lenTime > INT_MAX / 12 )
778 0 : return -8;
779 0 : lenTime = lenTime * 12;
780 0 : timeRangeUnit = 1;
781 0 : } else if (timeRangeUnit == 13) {
782 0 : lenTime = (sInt4) (lenTime / 3600.);
783 0 : timeRangeUnit = 1;
784 0 : } else if (timeRangeUnit == 3) { /* month */
785 : /* Actually use the timeRangeUnit == 3 */
786 : /*
787 : lenTime = (inv->validTime - Clock_AddMonthYear (inv->validTime, -1 * lenTime, 0)) / 3600.;
788 : timeRangeUnit = 1;
789 : */
790 0 : } else if (timeRangeUnit == 4) { /* month */
791 : /* Actually use the timeRangeUnit == 4 */
792 : /*
793 : lenTime = (inv->validTime - Clock_AddMonthYear (inv->validTime, 0, -1 * lenTime)) / 3600.;
794 : timeRangeUnit = 1;
795 : */
796 0 : } else if (timeRangeUnit == 5) { /* decade */
797 0 : if( lenTime < INT_MIN / 10 || lenTime > INT_MAX / 10 )
798 0 : return -8;
799 0 : lenTime = lenTime * 10;
800 0 : timeRangeUnit = 4;
801 : /*
802 : lenTime = (inv->validTime - Clock_AddMonthYear (inv->validTime, 0, -10 * lenTime)) / 3600.;
803 : timeRangeUnit = 1;
804 : */
805 0 : } else if (timeRangeUnit == 6) { /* normal */
806 0 : if( lenTime < INT_MIN / 30 || lenTime > INT_MAX / 30 )
807 0 : return -8;
808 0 : lenTime = lenTime * 30;
809 0 : timeRangeUnit = 4;
810 : /*
811 : lenTime = (inv->validTime - Clock_AddMonthYear (inv->validTime, 0, -30 * lenTime)) / 3600.;
812 : timeRangeUnit = 1;
813 : */
814 0 : } else if (timeRangeUnit == 7) { /* century */
815 0 : if( lenTime < INT_MIN / 100 || lenTime > INT_MAX / 100 )
816 0 : return -8;
817 0 : lenTime = lenTime * 100;
818 0 : timeRangeUnit = 4;
819 : /*
820 : lenTime = (inv->validTime - Clock_AddMonthYear (inv->validTime, 0, -100 * lenTime)) / 3600.;
821 : timeRangeUnit = 1;
822 : */
823 : } else {
824 0 : printf ("Can't handle this timeRangeUnit\n");
825 : //myAssert (timeRangeUnit == 1);
826 0 : return -8;
827 : }
828 344 : if (lenTime == GRIB2MISSING_s4) {
829 0 : lenTime = 0;
830 : }
831 :
832 344 : if ((templat == GS4_RADAR) || (templat == GS4_SATELLITE)
833 344 : || (templat == GS4_SATELLITE_SYNTHETIC)
834 337 : || (templat == 254) || (templat == 1000) || (templat == 1001)
835 337 : || (templat == 1002)) {
836 7 : fstSurfValue = 0;
837 7 : f_fstValue = 0;
838 7 : fstSurfType = 0;
839 7 : sndSurfValue = 0;
840 7 : f_sndValue = 0;
841 : } else {
842 337 : if( secLen < nOffset + 31 - 5 + 4)
843 0 : return -8;
844 337 : fstSurfType = (*buffer)[nOffset + 23 - 5];
845 337 : unsigned char u_scale = ((unsigned char*)(*buffer))[nOffset + 24 - 5];
846 337 : scale = (u_scale & 0x80) ? -(u_scale & 0x7f) : u_scale;
847 : unsigned int u_value;
848 337 : MEMCPY_BIG (&u_value, *buffer + nOffset + 25 - 5, sizeof (u_value));
849 337 : value = (u_value & 0x80000000U) ? -static_cast<int>(u_value & 0x7fffffff) : static_cast<int>(u_value);
850 337 : if ((value == GRIB2MISSING_s4) || (scale == GRIB2MISSING_s1) ||
851 : (fstSurfType == GRIB2MISSING_u1)) {
852 12 : fstSurfValue = 0;
853 12 : f_fstValue = 1;
854 : } else {
855 325 : fstSurfValue = value * pow (10.0, (int) (-1 * scale));
856 325 : f_fstValue = 1;
857 : }
858 337 : sndSurfType = (*buffer)[nOffset + 29 - 5];
859 337 : u_scale = ((unsigned char*)(*buffer))[nOffset + 30 - 5];
860 337 : scale = (u_scale & 0x80) ? -(u_scale & 0x7f) : u_scale;
861 337 : MEMCPY_BIG (&u_value, *buffer + nOffset + 31 - 5, sizeof (u_value));
862 337 : value = (u_value & 0x80000000U) ? -static_cast<int>(u_value & 0x7fffffff) : static_cast<int>(u_value);
863 337 : if ((value == GRIB2MISSING_s4) || (scale == GRIB2MISSING_s1) ||
864 : (sndSurfType == GRIB2MISSING_u1)) {
865 337 : sndSurfValue = 0;
866 337 : f_sndValue = 0;
867 : } else {
868 0 : sndSurfValue = value * pow (10.0, -1 * scale);
869 0 : f_sndValue = 1;
870 : }
871 : }
872 :
873 : /* Find out what the name of this variable is. */
874 344 : ParseElemName (mstrVersion, center, subcenter, prodType, templat, cat, subcat,
875 : lenTime, timeRangeUnit, statProcessID, timeIncrType, genID, probType, lowerProb,
876 : upperProb,
877 : derivedFcst,
878 : &(inv->element), &(inv->comment),
879 : &(inv->unitName), &convert, percentile, genProcess,
880 : f_fstValue, fstSurfValue, f_sndValue, sndSurfValue);
881 :
882 344 : if (! f_fstValue) {
883 7 : reallocSprintf (&(inv->shortFstLevel), "0 undefined");
884 7 : reallocSprintf (&(inv->longFstLevel), "0.000[-] undefined ()");
885 : } else {
886 337 : ParseLevelName (center, subcenter, fstSurfType, fstSurfValue,
887 : f_sndValue, sndSurfValue, &(inv->shortFstLevel),
888 : &(inv->longFstLevel));
889 : }
890 :
891 : /* Jump past section 5. */
892 344 : sectNum = 5;
893 344 : if (GRIB2SectJump (fp, gribLen, §Num, &secLen) != 0) {
894 2 : errSprintf ("ERROR: Problems Jumping past section 5\n");
895 2 : return -9;
896 : }
897 : /* Jump past section 6. */
898 342 : sectNum = 6;
899 342 : if (GRIB2SectJump (fp, gribLen, §Num, &secLen) != 0) {
900 0 : errSprintf ("ERROR: Problems Jumping past section 6\n");
901 0 : return -10;
902 : }
903 : /* Jump past section 7. */
904 342 : sectNum = 7;
905 342 : if (GRIB2SectJump (fp, gribLen, §Num, &secLen) != 0) {
906 0 : errSprintf ("ERROR: Problems Jumping past section 7\n");
907 0 : return -11;
908 : }
909 342 : return 0;
910 : }
911 :
912 : /*****************************************************************************
913 : * GRIB2Inventory() -- Review 12/2002
914 : *
915 : * Arthur Taylor / MDL
916 : *
917 : * PURPOSE
918 : * Fills out an inventory structure for each GRIB message in a GRIB file,
919 : * without calling the FORTRAN routines to unpack the message. It returns
920 : * the number of messages it found, or a negative number signifying an error.
921 : *
922 : * ARGUMENTS
923 : * filename = File to do the inventory of. (Input)
924 : * Inv = The resultant array of inventories. (Output)
925 : * LenInv = Length of the Array Inv (Output)
926 : * numMsg = # of messages to inventory (0 = all, 1 = just first) (In)
927 : * msgNum = MsgNum to start with, MsgNum of last message (Input/Output)
928 : *
929 : * FILES/DATABASES:
930 : * Opens a GRIB2 file for reading given its filename.
931 : *
932 : * RETURNS: int (could use errSprintf())
933 : * +# = number of GRIB2 messages in the file.
934 : * -1 = Problems opening file for read.
935 : * -2 = Problems in section 0
936 : * -3 = Ran out of file.
937 : * -4 = Problems Reading in section 1
938 : * -5 = Problems Reading in section 2 or 3
939 : * -6 = Problems Reading in section 3
940 : * -7 = Problems Reading in section 4
941 : * -8 = Problems Parsing section 4.
942 : * -9 = Problems Reading in section 5
943 : * -10 = Problems Reading in section 6
944 : * -11 = Problems Reading in section 7
945 : * -12 = Problems inventory'ing a GRIB1 record
946 : * -13 = Problems inventory'ing a TDLP record
947 : *
948 : * HISTORY
949 : * 9/2002 Arthur Taylor (MDL/RSIS): Created.
950 : * 11/2002 AAT: Revised.
951 : * 12/2002 (TK,AC,TB,&MS): Code Review.
952 : * 3/2003 AAT: Corrected some satellite type mistakes.
953 : * 3/2003 AAT: Implemented multiple grid inventories in the same GRIB2
954 : * message.
955 : * 4/2003 AAT: Started adding GRIB1 support
956 : * 6/2003 Matthew T. Kallio (matt at wunderground.com):
957 : * "wmo" dimension increased to WMO_HEADER_LEN + 1 (for '\0' char)
958 : * 7/2003 AAT: Added numMsg so we can quickly find the reference time for
959 : * a file by inventorying just the first message.
960 : * 8/2003 AAT: Adjusted use of GRIB_LIMIT to only affect the first message
961 : * after we know we have a GRIB file, we don't want "trailing" bytes
962 : * to break the program.
963 : * 8/2003 AAT: switched fileLen to only be computed for an error message.
964 : * 8/2003 AAT: curTot no longer serves a purpose.
965 : * 5/2004 AAT: Added a check for section number 2..8 for the repeated
966 : * section (otherwise error)
967 : * 10/2004 AAT: Added ability to inventory TDLP records.
968 : *
969 : * NOTES
970 : *****************************************************************************
971 : */
972 310 : int GRIB2Inventory (VSILFILE *fp, inventoryType **Inv, uInt4 *LenInv,
973 : int numMsg, int *MsgNum)
974 : {
975 310 : vsi_l_offset offset = 0; /* Where we are in the file. */
976 : sInt4 msgNum; /* Which GRIB2 message we are on. */
977 : uInt4 gribLen; /* Length of the current GRIB message. */
978 : uInt4 secLen; /* Length of current section. */
979 : sChar sectNum; /* Which section we are reading. */
980 : char *buff; /* Holds the info between records. */
981 : uInt4 buffLen; /* Length of info between records. */
982 : sInt4 sect0[SECT0LEN_WORD]; /* Holds the current Section 0. */
983 310 : char *buffer = nullptr; /* Holds a given section. */
984 310 : uInt4 bufferLen = 0; /* Size of buffer. */
985 : inventoryType *inv; /* Local ptr to Inv to reduce ptr confusion. */
986 : inventoryType *lastInv; /* Used to point to last inventory record when
987 : * there are multiple grids in the same message. */
988 : wordType word; /* Used to parse the prodType out of Sect 0. */
989 : int ans; /* The return error code of ReadSect0. */
990 : #if 0
991 : char *msg; /* Used to pop messages off the error Stack. */
992 : #endif
993 : int version; /* Which version of GRIB is in this message. */
994 : uChar prodType; /* Which GRIB2 type of product, 0 is meteo, 1 is
995 : * hydro, 2 is land, 3 is space, 10 is oceanographic.
996 : */
997 : int grib_limit; /* How many bytes to look for before the first "GRIB"
998 : * in the file. If not found, is not a GRIB file. */
999 : char c; /* Determine if end of the file without fileLen. */
1000 : #ifdef DEBUG
1001 : vsi_l_offset fileLen; /* Length of the GRIB2 file. */
1002 : #endif
1003 : unsigned short int center, subcenter; /* Who produced it. */
1004 : uChar mstrVersion; /* The master table version (is it 255?) */
1005 : // char *ptr; /* used to find the file extension. */
1006 :
1007 310 : grib_limit = GRIB_LIMIT;
1008 : /*
1009 : if (filename != NULL) {
1010 : //if ((fp = fopen (filename, "rb")) == NULL) {
1011 : // errSprintf ("ERROR: Problems opening %s for read.", filename);
1012 : // return -1;
1013 : //}
1014 : //fp = FileDataSource(filename);
1015 : ptr = strrchr (filename, '.');
1016 : if (ptr != NULL) {
1017 : if (strcmp (ptr, ".tar") == 0) {
1018 : grib_limit = 5000;
1019 : }
1020 : }
1021 : } else {
1022 : //fp = stdin; // TODO!!
1023 : }
1024 : */
1025 310 : msgNum = *MsgNum;
1026 :
1027 310 : buff = nullptr;
1028 310 : buffLen = 0;
1029 :
1030 721 : while (VSIFReadL(&c, sizeof(char), 1, fp) == 1) {
1031 421 : VSIFSeekL(fp, VSIFTellL(fp) - sizeof(char), SEEK_SET);
1032 : // ungetc (c, fp);
1033 : /* msgNum++ done first so any error messages range from 1..n, instead
1034 : * of 0.. n-1. Note msgNum should end up as n not (n-1) */
1035 421 : msgNum++;
1036 : /* Used when testing inventory of large TDLPack files. */
1037 : /*
1038 : #ifdef DEBUG
1039 : myAssert (msgNum < 32500L);
1040 : if (msgNum % 10 == 0) {
1041 : printf ("%ld :: %f\n", msgNum, clock () / (double) CLOCKS_PER_SEC);
1042 : }
1043 : #endif
1044 : */
1045 : /* Allow 2nd, 3rd, etc messages to have no limit to finding "GRIB". */
1046 421 : if (msgNum > 1) {
1047 111 : grib_limit = -1;
1048 : }
1049 : /* Read in the wmo header and sect0. */
1050 421 : if (ReadSECT0 (fp, &buff, &buffLen, grib_limit, sect0, &gribLen,
1051 421 : &version) < 0) {
1052 1 : if (msgNum == 1) {
1053 : /* Handle case where we couldn't find 'GRIB' in the message. */
1054 0 : preErrSprintf ("Inside GRIB2Inventory, Message # %d\n", msgNum);
1055 0 : free (buffer);
1056 0 : free (buff);
1057 : //fclose (fp);
1058 0 : return -2;
1059 : } else {
1060 : /* Handle case where there are trailing bytes. */
1061 : #if 0
1062 : msg = errSprintf (nullptr);
1063 : printf ("Warning: Inside GRIB2Inventory, Message # %d\n",
1064 : msgNum);
1065 : printf ("%s", msg);
1066 : free (msg);
1067 : #ifdef DEBUG
1068 : /* find out how big the file is. */
1069 : VSIFSeekL (fp, 0L, SEEK_END);
1070 : fileLen = VSIFTellL(fp);
1071 : /* fseek (fp, 0L, SEEK_SET); */
1072 : printf ("There were " CPL_FRMT_GUIB " trailing bytes in the file.\n",
1073 : fileLen - offset);
1074 : #endif
1075 : #else
1076 : #ifdef DEBUG
1077 : /* find out how big the file is. */
1078 1 : VSIFSeekL (fp, 0L, SEEK_END);
1079 1 : fileLen = VSIFTellL(fp);
1080 1 : CPLDebug("GRIB", CPL_FRMT_GUIB " trailing byte(s) found after message #%d", fileLen - offset, msgNum);
1081 : #else
1082 : CPLDebug("GRIB", "Trailing byte(s) found after message #%d", msgNum);
1083 : #endif
1084 : #endif
1085 1 : free (buffer);
1086 1 : free (buff);
1087 : //fclose (fp);
1088 1 : msgNum --;
1089 1 : *MsgNum = msgNum;
1090 1 : return msgNum;
1091 : }
1092 : }
1093 :
1094 : /* Make room for this GRIB message in the inventory list. */
1095 420 : *LenInv = *LenInv + 1;
1096 420 : *Inv = (inventoryType *) realloc ((void *) *Inv,
1097 420 : *LenInv * sizeof (inventoryType));
1098 420 : inv = *Inv + (*LenInv - 1);
1099 :
1100 : /* Start parsing the message. */
1101 420 : inv->GribVersion = version;
1102 420 : inv->msgNum = msgNum;
1103 420 : inv->subgNum = 0;
1104 420 : inv->start = offset;
1105 420 : inv->element = nullptr;
1106 420 : inv->comment = nullptr;
1107 420 : inv->unitName = nullptr;
1108 420 : inv->shortFstLevel = nullptr;
1109 420 : inv->longFstLevel = nullptr;
1110 :
1111 420 : if (version == 1) {
1112 66 : if (GRIB1_Inventory (fp, gribLen, inv) != 0) {
1113 0 : preErrSprintf ("Inside GRIB2Inventory \n");
1114 0 : free (buffer);
1115 0 : free (buff);
1116 : //fclose (fp);
1117 0 : return -12;
1118 : }
1119 : #if 0
1120 : /* tdlpack is no longer supported by GDAL */
1121 : } else if (version == -1) {
1122 : if (TDLP_Inventory (fp, gribLen, inv) != 0) {
1123 : preErrSprintf ("Inside GRIB2Inventory \n");
1124 : free (buffer);
1125 : free (buff);
1126 : //fclose (fp);
1127 : return -13;
1128 : }
1129 : #endif
1130 : } else {
1131 354 : word.li = sect0[1];
1132 354 : prodType = word.buffer[2];
1133 :
1134 : /* Read section 1 into buffer. */
1135 354 : sectNum = 1;
1136 354 : if (GRIB2SectToBuffer (fp, gribLen, §Num, &secLen, &bufferLen,
1137 354 : &buffer) != 0) {
1138 2 : errSprintf ("ERROR: Problems with section 1\n");
1139 2 : free (buffer);
1140 2 : free (buff);
1141 : //fclose (fp);
1142 2 : return -4;
1143 : }
1144 : /* Parse the interesting data out of sect 1. */
1145 352 : if( bufferLen < 13 - 5 + 7 )
1146 : {
1147 0 : errSprintf ("ERROR: Problems with section 1\n");
1148 0 : free (buffer);
1149 0 : free (buff);
1150 0 : return -4;
1151 : }
1152 : /* InventoryParseTime reads 7 bytes */
1153 352 : if( InventoryParseTime (buffer + 13 - 5, &(inv->refTime)) < 0 )
1154 : {
1155 0 : errSprintf ("ERROR: Problems with section 1: invalid refTime\n");
1156 0 : free (buffer);
1157 0 : free (buff);
1158 0 : return -4;
1159 : }
1160 352 : MEMCPY_BIG (¢er, buffer + 6 - 5, sizeof (short int));
1161 352 : MEMCPY_BIG (&subcenter, buffer + 8 - 5, sizeof (short int));
1162 352 : MEMCPY_BIG (&mstrVersion, buffer + 10 - 5, sizeof (uChar));
1163 :
1164 352 : sectNum = 2;
1165 2 : do {
1166 : /* Look at sections 2 to 7 */
1167 354 : if ((ans = GRIB2Inventory2to7 (sectNum, fp, gribLen, &bufferLen,
1168 : &buffer, inv, prodType, center,
1169 354 : subcenter, mstrVersion)) < 0) {
1170 : //fclose (fp);
1171 7 : free (buffer);
1172 7 : free (buff);
1173 7 : return ans;
1174 : }
1175 : /* Try to read section 8. If it is "7777" = 926365495 regardless
1176 : * of endian'ness then we have a simple message, otherwise it is
1177 : * complex, and we need to read more. */
1178 347 : if (FREAD_BIG (&secLen, sizeof (sInt4), 1, fp) != 1) {
1179 0 : errSprintf ("ERROR: Ran out of file looking for Sect 8.\n");
1180 0 : free (buffer);
1181 0 : free (buff);
1182 : // fclose (fp);
1183 0 : return -4;
1184 : }
1185 347 : if (secLen == 926365495L) {
1186 345 : sectNum = 8;
1187 : } else {
1188 2 : if (VSIFReadL (§Num, sizeof (char), 1, fp) != 1) {
1189 0 : errSprintf ("ERROR: Ran out of file looking for "
1190 : "subMessage.\n");
1191 0 : free (buffer);
1192 0 : free (buff);
1193 : //fclose (fp);
1194 0 : return -4;
1195 : }
1196 2 : if ((sectNum < 2) || (sectNum > 7)) {
1197 0 : errSprintf ("ERROR (GRIB2Inventory): Couldn't find the end"
1198 : " of message\n");
1199 0 : errSprintf ("and it doesn't appear to repeat sections.\n");
1200 0 : errSprintf ("so it is probably an ASCII / binary bug\n");
1201 0 : free (buffer);
1202 0 : free (buff);
1203 : //fclose (fp);
1204 0 : return -4;
1205 : }
1206 2 : VSIFSeekL(fp, VSIFTellL(fp) - 5, SEEK_SET);
1207 : /* Make room for the next part of this GRIB message in the
1208 : * inventory list. This is for when we have sub-grids. */
1209 2 : *LenInv = *LenInv + 1;
1210 2 : *Inv = (inventoryType *) realloc ((void *) *Inv,
1211 2 : *LenInv *
1212 : sizeof (inventoryType));
1213 2 : inv = *Inv + (*LenInv - 1);
1214 2 : lastInv = *Inv + (*LenInv - 2);
1215 :
1216 2 : inv->GribVersion = version;
1217 2 : inv->msgNum = msgNum;
1218 2 : inv->subgNum = lastInv->subgNum + 1;
1219 2 : inv->start = offset;
1220 2 : inv->element = nullptr;
1221 2 : inv->comment = nullptr;
1222 2 : inv->unitName = nullptr;
1223 2 : inv->shortFstLevel = nullptr;
1224 2 : inv->longFstLevel = nullptr;
1225 :
1226 2 : word.li = sect0[1];
1227 2 : prodType = word.buffer[2];
1228 2 : inv->refTime = lastInv->refTime;
1229 : }
1230 347 : } while (sectNum != 8);
1231 : }
1232 :
1233 : /* added to inventory either first msgNum messages, or all messages */
1234 411 : if (numMsg == msgNum) {
1235 0 : break;
1236 : }
1237 : {
1238 : uInt4 increment;
1239 : /* Continue on to the next GRIB2 message. */
1240 : #if 0
1241 : /* tdlpack is no longer supported by GDAL */
1242 : if (version == -1) {
1243 : /* TDLPack uses 4 bytes for FORTRAN record size, then another 8
1244 : * bytes for the size of the record (so FORTRAN can see it), then
1245 : * the data rounded up to an 8 byte boundary, then a trailing 4
1246 : * bytes for a final FORTRAN record size. However it only stores
1247 : * in_ the gribLen the non-rounded amount, so we need to take care
1248 : * of the rounding, and the trailing 4 bytes here. */
1249 : increment = buffLen + ((uInt4) ceil (gribLen / 8.0)) * 8 + 4;
1250 : } else
1251 : #endif
1252 : {
1253 411 : if( buffLen > UINT_MAX - gribLen )
1254 0 : break;
1255 411 : increment = buffLen + gribLen;
1256 : }
1257 411 : if( /* increment < buffLen || */ increment > (VSI_L_OFFSET_MAX - offset) )
1258 0 : break;
1259 411 : offset += increment;
1260 : }
1261 411 : VSIFSeekL(fp, offset, SEEK_SET);
1262 : }
1263 300 : free (buffer);
1264 300 : free (buff);
1265 : //fclose (fp);
1266 300 : *MsgNum = msgNum;
1267 300 : return msgNum;
1268 : }
1269 :
1270 0 : int GRIB2RefTime (const char *filename, double *refTime)
1271 : {
1272 0 : VSILFILE * fp = nullptr;
1273 0 : vsi_l_offset offset = 0; /* Where we are in the file. */
1274 : sInt4 msgNum; /* Which GRIB2 message we are on. */
1275 : uInt4 gribLen; /* Length of the current GRIB message. */
1276 : uInt4 secLen; /* Length of current section. */
1277 : sChar sectNum; /* Which section we are reading. */
1278 : char *buff; /* Holds the info between records. */
1279 : uInt4 buffLen; /* Length of info between records. */
1280 : sInt4 sect0[SECT0LEN_WORD]; /* Holds the current Section 0. */
1281 0 : char *buffer = nullptr; /* Holds a given section. */
1282 0 : uInt4 bufferLen = 0; /* Size of buffer. */
1283 : /* wordType word; */ /* Used to parse the prodType out of Sect 0. */
1284 : int ans; /* The return error code of ReadSect0. */
1285 : char *msg; /* Used to pop messages off the error Stack. */
1286 : int version; /* Which version of GRIB is in this message. */
1287 : /* uChar prodType; */ /* Which GRIB2 type of product, 0 is meteo, 1 is
1288 : * hydro, 2 is land, 3 is space, 10 is oceanographic.
1289 : */
1290 : int grib_limit; /* How many bytes to look for before the first "GRIB"
1291 : * in the file. If not found, is not a GRIB file. */
1292 : char c; /* Determine if end of the file without fileLen. */
1293 : #ifdef DEBUG
1294 : vsi_l_offset fileLen; /* Length of the GRIB2 file. */
1295 : #endif
1296 : const char *ptr; /* used to find the file extension. */
1297 : double refTime1;
1298 :
1299 0 : grib_limit = GRIB_LIMIT;
1300 : /*if (filename != NULL)*/ {
1301 0 : if ((fp = VSIFOpenL(filename, "rb")) == nullptr) {
1302 0 : return -1;
1303 : }
1304 0 : ptr = strrchr (filename, '.');
1305 0 : if (ptr != nullptr) {
1306 0 : if (strcmp (ptr, ".tar") == 0) {
1307 0 : grib_limit = 5000;
1308 : }
1309 : }
1310 : }// else {
1311 : // fp = stdin; // TODO!!
1312 : //}
1313 0 : msgNum = 0;
1314 :
1315 0 : buff = nullptr;
1316 0 : buffLen = 0;
1317 0 : while (VSIFReadL(&c, sizeof(char), 1, fp) == 1) {
1318 0 : VSIFSeekL(fp, VSIFTellL(fp) - sizeof(char), SEEK_SET);
1319 :
1320 : /* msgNum++ done first so any error messages range from 1..n, instead
1321 : * of 0.. n-1. Note msgNum should end up as n not (n-1) */
1322 0 : msgNum++;
1323 : /* Make it so the second, third, etc messages have no limit to finding
1324 : * the "GRIB" keyword. */
1325 0 : if (msgNum > 1) {
1326 0 : grib_limit = -1;
1327 : }
1328 : /* Read in the wmo header and sect0. */
1329 0 : if ((ans = ReadSECT0 (fp, &buff, &buffLen, grib_limit, sect0, &gribLen,
1330 0 : &version)) < 0) {
1331 0 : if (msgNum == 1) {
1332 : /* Handle case where we couldn't find 'GRIB' in the message. */
1333 0 : preErrSprintf ("Inside GRIB2RefTime, Message # %d\n", msgNum);
1334 0 : free (buffer);
1335 0 : free (buff);
1336 : //fclose (fp);
1337 0 : return -2;
1338 : } else {
1339 : /* Handle case where there are trailing bytes. */
1340 0 : msg = errSprintf (nullptr);
1341 0 : printf ("Warning: Inside GRIB2RefTime, Message # %d\n", msgNum);
1342 0 : printf ("%s", msg);
1343 0 : free (msg);
1344 : #ifdef DEBUG
1345 : /* find out how big the file is. */
1346 0 : VSIFSeekL(fp, 0L, SEEK_END);
1347 0 : fileLen = VSIFTellL(fp);
1348 : /* fseek (fp, 0L, SEEK_SET); */
1349 0 : printf ("There were " CPL_FRMT_GUIB " trailing bytes in the file.\n",
1350 : fileLen - offset);
1351 : #endif
1352 0 : free (buffer);
1353 0 : free (buff);
1354 : //fclose (fp);
1355 0 : return msgNum;
1356 : }
1357 : }
1358 :
1359 0 : if (version == 1) {
1360 0 : if (GRIB1_RefTime (fp, gribLen, &(refTime1)) != 0) {
1361 0 : preErrSprintf ("Inside GRIB1_RefTime\n");
1362 0 : free (buffer);
1363 0 : free (buff);
1364 : //fclose (fp);
1365 0 : return -12;
1366 : }
1367 : #if 0
1368 : /* tdlpack is no longer supported by GDAL */
1369 : } else if (version == -1) {
1370 : if (TDLP_RefTime (fp, gribLen, &(refTime1)) != 0) {
1371 : preErrSprintf ("Inside TDLP_RefTime\n");
1372 : free (buffer);
1373 : free (buff);
1374 : //fclose (fp);
1375 : return -13;
1376 : }
1377 : #endif
1378 : } else {
1379 : /* word.li = sect0[1]; */
1380 : /* prodType = word.buffer[2]; */
1381 :
1382 : /* Read section 1 into buffer. */
1383 0 : sectNum = 1;
1384 0 : if (GRIB2SectToBuffer (fp, gribLen, §Num, &secLen, &bufferLen,
1385 0 : &buffer) != 0) {
1386 0 : errSprintf ("ERROR: Problems with section 1\n");
1387 0 : free (buffer);
1388 : //fclose (fp);
1389 0 : return -4;
1390 : }
1391 : /* Parse the interesting data out of sect 1. */
1392 0 : if( InventoryParseTime (buffer + 13 - 5, &(refTime1)) < 0 )
1393 0 : refTime1 = 0.0;
1394 : }
1395 0 : if (msgNum == 1) {
1396 0 : *refTime = refTime1;
1397 : } else {
1398 0 : if (*refTime > refTime1) {
1399 0 : *refTime = refTime1;
1400 : }
1401 : }
1402 :
1403 : /* Continue on to the next GRIB2 message. */
1404 0 : offset += gribLen + buffLen;
1405 0 : VSIFSeekL (fp, offset, SEEK_SET);
1406 : }
1407 0 : free (buffer);
1408 0 : free (buff);
1409 : //fclose (fp);
1410 0 : return 0;
1411 : }
|
