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