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