1 /*****************************************************************************
2 * degrib2.c
3 *
4 * DESCRIPTION
5 * This file contains the main driver routines to call the unpack grib2
6 * library functions. It also contains the code needed to figure out the
7 * dimensions of the arrays before calling 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 <limits>
23
24 #include "myassert.h"
25 #include "myerror.h"
26 #include "tendian.h"
27 #include "meta.h"
28 #include "metaname.h"
29 //#include "write.h"
30 #include "degrib2.h"
31 #include "degrib1.h"
32 #if 0
33 /* tdlpack is no longer supported by GDAL */
34 #include "tdlpack.h"
35 #endif
36 #include "grib2api.h"
37 //#include "mymapf.h"
38 #include "clock.h"
39
40 #define GRIB_UNSIGN_INT3(a,b,c) ((a<<16)+(b<<8)+c)
41
42 /*****************************************************************************
43 * ReadSect0() -- Review 12/2002
44 *
45 * Arthur Taylor / MDL
46 *
47 * PURPOSE
48 * Looks for the next GRIB message, by looking for the keyword "GRIB". It
49 * expects the message in "expect" bytes from the start, but could find the
50 * message in "expect2" bytes or 0 bytes from the start. Returns -1 if it
51 * can't find "GRIB", 1 if "GRIB" is not 0, "expect", or "expect2" bytes from
52 * the start.
53 * It stores the bytes it reads (a max of "expect") up to but not including
54 * the 'G' in "GRIB" in wmo.
55 *
56 * After it finds section 0, it then parses the 16 bytes that make up
57 * section 0 so that it can return the length of the entire GRIB message.
58 *
59 * When done, it sets fp to point to the end of Sect0.
60 *
61 * The reason for this procedure is so that we can read in the size of the
62 * grib message, and thus allocate enough memory to read the message in before
63 * making it Big endian, and passing it to the library for unpacking.
64 *
65 * ARGUMENTS
66 * fp = A pointer to an opened file in which to read.
67 * When done, this points to the start of section 1. (Input/Output)
68 * buff = The data between messages. (Input/Output)
69 * buffLen = The length of buff (Output)
70 * limit = How many bytes to read before giving up and stating it is not
71 * a proper message. (-1 means no limit). (Input)
72 * sect0 = The read in Section 0 (as seen on disk). (Output)
73 * gribLen = Length of this GRIB message. (Output)
74 * version = 1 if GRIB1 message, 2 if GRIB2 message, -1 if TDLP message.
75 * (Output)
76 * expect = The expected number of bytes to find "GRIB" in. (Input)
77 * expect2 = The second possible number of bytes to find "GRIB" in. (Input)
78 * wmo = Assumed allocated to be at least size "expect".
79 * Holds the bytes before the first "GRIB" message.
80 * expect should be > expect2, but is up to caller (Output)
81 * wmoLen = Length of wmo (total number of bytes read - SECT0LEN_WORD * 4).
82 * (Output)
83 *
84 * FILES/DATABASES:
85 * An already opened "GRIB2" File
86 *
87 * RETURNS: int (could use errSprintf())
88 * 1 = Length of wmo was != 0 and was != expect
89 * 0 = OK
90 * -1 = Couldn't find "GRIB" part of message.
91 * -2 = Ran out of file while reading this section.
92 * -3 = Grib version was not 1 or 2.
93 * -4 = Most significant sInt4 of GRIB length was not 0
94 * -5 = Grib message length was <= 16 (can't be smaller than just sect 0)
95 *
96 * HISTORY
97 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
98 * 11/2002 AAT: Combined with ReadWMOHeader
99 * 12/2002 (TK,AC,TB,&MS): Code Review.
100 * 1/2003 AAT: Bug found. wmo access out of bounds of expect when setting
101 * the /0 element, if wmoLen > expect.
102 * 4/2003 AAT: Added ability to handle GRIB version 1.
103 * 5/2003 AAT: Added limit option.
104 * 8/2003 AAT: Removed dependence on offset, and fileLen.
105 * 10/2004 AAT: Modified to allow for TDLP files
106 *
107 * NOTES
108 * 1a) 1196575042L == ASCII representation of "GRIB" (GRIB in MSB)
109 * 1b) 1112101447L == ASCII representation of "BIRG" (GRIB in LSB)
110 * 1c) 1413762128L == ASCII representation of "TDLP" (TDLP in MSB)
111 * 1d) 1347175508L == ASCII representation of "PLDT" (TDLP in LSB)
112 * 2) Takes advantage of the wordType to check that the edition is correct.
113 * 3) May want to return prodType.
114 * 4) WMO_HEADER_ORIG_LEN was added for backward compatibility... should be
115 * removed when we no longer use old format. (say in a year from 11/2002)
116 *
117 *****************************************************************************
118 */
ReadSECT0(VSILFILE * fp,char ** buff,uInt4 * buffLen,sInt4 limit,sInt4 sect0[SECT0LEN_WORD],uInt4 * gribLen,int * version)119 int ReadSECT0 (VSILFILE *fp, char **buff, uInt4 *buffLen, sInt4 limit,
120 sInt4 sect0[SECT0LEN_WORD], uInt4 *gribLen, int *version)
121 {
122 typedef union {
123 sInt4 li;
124 unsigned char buffer[4];
125 } wordType;
126
127 uChar gribMatch = 0; /* Counts how many letters in GRIB we've matched. */
128 #if 0
129 /* tdlpack is no longer supported by GDAL */
130 uChar tdlpMatch = 0; /* Counts how many letters in TDLP we've matched. */
131 #endif
132 wordType word; /* Used to check that the edition is correct. */
133 uInt4 curLen; /* Where we currently are in buff. */
134 uInt4 i; /* Used to loop over the first few char's */
135 uInt4 stillNeed; /* Number of bytes still needed to get 1st 8 bytes of
136 * message into memory. */
137
138 /* Get first 8 bytes. If GRIB we don't care. If TDLP, this is the length
139 * of record. Read at least 1 record (length + 2 * 8) + 8 (next record
140 * length) + 8 bytes before giving up. */
141 curLen = 8;
142 if (*buffLen < curLen) {
143 *buffLen = curLen;
144 *buff = (char *) realloc ((void *) *buff, *buffLen * sizeof (char));
145 }
146 if (VSIFReadL(*buff, sizeof (char), curLen, fp) != curLen) {
147 errSprintf ("ERROR: Couldn't find 'GRIB' or 'TDLP'\n");
148 return -1;
149 }
150 /*
151 Can't do the following because we don't know if the file is a GRIB file or
152 not, or if it was a FORTRAN file.
153 if (limit > 0) {
154 MEMCPY_BIG (&recLen, *buff, 4);
155 limit = (limit > recLen + 32) ? limit : recLen + 32;
156 }
157 */
158 while (
159 #if 0
160 /* tdlpack is no longer supported by GDAL */
161 (tdlpMatch != 4) &&
162 #endif
163 (gribMatch != 4)) {
164 for (i = curLen - 8; i + 7 < curLen; i++) {
165 if ((*buff)[i] == 'G') {
166 if (((*buff)[i + 1] == 'R') && ((*buff)[i + 2] == 'I') &&
167 ((*buff)[i + 3] == 'B')) {
168 if (((*buff)[i + 7] == 1) ||
169 ((*buff)[i + 7] == 2)) {
170 gribMatch = 4;
171 break;
172 }
173 }
174 }
175 #if 0
176 /* tdlpack is no longer supported by GDAL */
177 else if ((*buff)[i] == 'T') {
178 if (((*buff)[i + 1] == 'D') && ((*buff)[i + 2] == 'L') &&
179 ((*buff)[i + 3] == 'P')) {
180 tdlpMatch = 4;
181 break;
182 }
183 }
184 #endif
185 }
186 stillNeed = i - (curLen - 8);
187 /* Read enough of message to have the first 8 bytes (including ID). */
188 if (stillNeed != 0) {
189 curLen += stillNeed;
190 if ((limit >= 0) && (curLen > (size_t) limit)) {
191 errSprintf ("ERROR: Couldn't find type in %ld bytes\n", limit);
192 *buffLen = curLen - stillNeed;
193 return -1;
194 }
195 if (*buffLen < curLen) {
196 myAssert (200 > stillNeed);
197 *buffLen = *buffLen + 200;
198 /* *buffLen = curLen; */
199 *buff = (char *) realloc ((void *) *buff,
200 *buffLen * sizeof (char));
201 }
202 if (VSIFReadL((*buff) + (curLen - stillNeed), sizeof (char), stillNeed, fp) != stillNeed) {
203 errSprintf ("ERROR: Ran out of file reading SECT0\n");
204 *buffLen = curLen;
205 return -1;
206 }
207 }
208 }
209 /* Following is needed because we are increasing buffLen at a rate of
210 * 200 (to save reallocs), so it may not actually line up with the length
211 * of buffer. curLen should always be the length of buffer. */
212 *buffLen = curLen;
213
214 /* curLen and (*buff) hold 8 bytes of section 0. */
215 curLen -= 8;
216 memcpy (&(sect0[0]), (*buff) + curLen, 4);
217 #ifdef DEBUG
218 #ifdef LITTLE_ENDIAN
219 myAssert ((sect0[0] == 1112101447L) || (sect0[0] == 1347175508L));
220 #else
221 myAssert ((sect0[0] == 1196575042L) || (sect0[0] == 1413762128L));
222 #endif
223 #endif
224 memcpy (&(sect0[1]), *buff + curLen + 4, 4);
225 /* Make sure we don't pass back part of "GRIB" in the buffer. */
226 (*buff)[curLen] = '\0';
227 *buffLen = curLen;
228
229 word.li = sect0[1];
230 #if 0
231 /* tdlpack is no longer supported by GDAL */
232 if (tdlpMatch == 4) {
233 if (word.buffer[3] != 0) {
234 errSprintf ("ERROR: unexpected version of TDLP in SECT0\n");
235 return -2;
236 }
237 *version = -1;
238 /* Find out the GRIB Message Length */
239 *gribLen = GRIB_UNSIGN_INT3 (word.buffer[0], word.buffer[1],
240 word.buffer[2]);
241 /* Min message size: GRIB1=52, TDLP=59, GRIB2=86. */
242 if (*gribLen < 59) {
243 errSprintf ("TDLP length %ld was < 59?\n", *gribLen);
244 return -5;
245 }
246 } else
247 #endif
248 if (word.buffer[3] == 1) {
249 *version = 1;
250 /* Find out the GRIB Message Length */
251 *gribLen = GRIB_UNSIGN_INT3 (word.buffer[0], word.buffer[1],
252 word.buffer[2]);
253 /* Min message size: GRIB1=52, TDLP=59, GRIB2=86. */
254 if (*gribLen < 52) {
255 errSprintf ("GRIB1 length %ld was < 52?\n", *gribLen);
256 return -5;
257 }
258 } else if (word.buffer[3] == 2) {
259 *version = 2;
260 /* Make sure we still have enough file for the rest of section 0. */
261 if (VSIFReadL(sect0 + 2, sizeof (sInt4), 2, fp) != 2) {
262 errSprintf ("ERROR: Ran out of file reading SECT0\n");
263 return -2;
264 }
265 if (sect0[2] != 0) {
266 errSprintf ("Most significant sInt4 of GRIB length was not 0?\n");
267 errSprintf ("This is either an error, or we have a single GRIB "
268 "message which is larger than 2^31 = 2,147,283,648 "
269 "bytes.\n");
270 return -4;
271 }
272 #ifdef LITTLE_ENDIAN
273 revmemcpy (gribLen, &(sect0[3]), sizeof (sInt4));
274 #else
275 *gribLen = sect0[3];
276 #endif
277 } else {
278 errSprintf ("ERROR: Not TDLPack, and Grib edition is not 1 or 2\n");
279 return -3;
280 }
281 return 0;
282 }
283
284 /*****************************************************************************
285 * FindGRIBMsg() -- Review 12/2002
286 *
287 * Arthur Taylor / MDL
288 *
289 * PURPOSE
290 * Jumps through a GRIB2 file looking for a specific message. Currently
291 * that message is determined by msgNum which is in the range of 1..n.
292 * In the future we may be searching based on projection or date.
293 *
294 * ARGUMENTS
295 * fp = The current GRIB2 file to look through. (Input)
296 * msgNum = Which message to look for. (Input)
297 * offset = Where in the file the message starts (this is before the
298 * wmo ASCII part if there is one.) (Output)
299 * curMsg = The current # of messages we have looked through. (In/Out)
300 *
301 * FILES/DATABASES:
302 * An already opened "GRIB2" File
303 *
304 * RETURNS: int (could use errSprintf())
305 * 0 = OK
306 * -1 = Problems reading Section 0.
307 * -2 = Ran out of file.
308 *
309 * HISTORY
310 * 11/2002 Arthur Taylor (MDL/RSIS): Created.
311 * 12/2002 (TK,AC,TB,&MS): Code Review.
312 * 6/2003 Matthew T. Kallio (matt@wunderground.com):
313 * "wmo" dimension increased to WMO_HEADER_LEN + 1 (for '\0' char)
314 * 8/2003 AAT: Removed dependence on offset and fileLen.
315 *
316 * NOTES
317 *****************************************************************************
318 */
FindGRIBMsg(VSILFILE * fp,int msgNum,sInt4 * offset,int * curMsg)319 int FindGRIBMsg (VSILFILE *fp, int msgNum, sInt4 *offset, int *curMsg)
320 {
321 int cnt; /* The current message we are looking at. */
322 char *buff = nullptr; /* Holds the info between records. */
323 uInt4 buffLen; /* Length of info between records. */
324 sInt4 sect0[SECT0LEN_WORD]; /* Holds the current Section 0. */
325 uInt4 gribLen; /* Length of the current GRIB message. */
326 int version; /* Which version of GRIB is in this message. */
327 char c; /* Determine if end of the file without fileLen. */
328 sInt4 jump; /* How far to jump to get to past GRIB message. */
329
330 cnt = *curMsg + 1;
331 buffLen = 0;
332 while (VSIFReadL(&c, sizeof(char), 1, fp) == 1) {
333 VSIFSeekL(fp, VSIFTellL(fp) - sizeof(char), SEEK_SET);
334
335 if (cnt >= msgNum) {
336 /* 12/1/2004 version 1.63 forgot to free buff */
337 free (buff);
338 *curMsg = cnt;
339 return 0;
340 }
341 /* Read section 0 to find gribLen and wmoLen. */
342 if (ReadSECT0 (fp, &buff, &buffLen, GRIB_LIMIT, sect0, &gribLen,
343 &version) < 0) {
344 preErrSprintf ("Inside FindGRIBMsg\n");
345 free (buff);
346 return -1;
347 }
348 myAssert ((version == 1) || (version == 2) || (version == -1));
349 /* Continue on to the next grib message. */
350 if ((version == 1) || (version == -1)) {
351 jump = gribLen - 8;
352 } else {
353 jump = gribLen - 16;
354 }
355 VSIFSeekL(fp, jump, SEEK_CUR);
356 *offset = *offset + gribLen + buffLen;
357 cnt++;
358 }
359 free (buff);
360 *curMsg = cnt - 1;
361 /* Return -2 since we reached the end of file. This may not be an error
362 * (multiple file option). */
363 return -2;
364 /*
365 errSprintf ("ERROR: Ran out of file looking for msgNum %d.\n", msgNum);
366 errSprintf (" Current msgNum %d\n", cnt);
367 */
368 }
369
370 /*****************************************************************************
371 * FindSectLen2to7() --
372 *
373 * Arthur Taylor / MDL
374 *
375 * PURPOSE
376 * Looks through a GRIB message and finds out the maximum size of each
377 * section. Simpler if there is only one grid in the message.
378 *
379 * ARGUMENTS
380 * c_ipack = The current GRIB2 message. (Input)
381 * gribLen = Length of c_ipack. (Input)
382 * ns = Array of section lengths. (Output)
383 * sectNum = Which section to start with. (Input)
384 * curTot = on going total read from c_ipack. (Input)
385 * nd2x3 = Total number of grid points (Output)
386 * table50 = Type of packing used. (See code table 5.0) (GS5_SIMPLE = 0,
387 * GS5_CMPLX = 2, GS5_CMPLXSEC = 3) (Output)
388 *
389 * FILES/DATABASES: None
390 *
391 * RETURNS: int (could use errSprintf())
392 * 0 = OK
393 * -1 = Ran of data in a section
394 * -2 = Section not properly labeled.
395 *
396 * HISTORY
397 * 3/2003 AAT: Created
398 *
399 * NOTES
400 * 1) Assumes that the pack method of multiple grids are the same.
401 *****************************************************************************
402 */
FindSectLen2to7(unsigned char * c_ipack,sInt4 gribLen,sInt4 ns[8],char sectNum,sInt4 * curTot,sInt4 * nd2x3,short int * table50)403 static int FindSectLen2to7 (unsigned char *c_ipack, sInt4 gribLen, sInt4 ns[8],
404 char sectNum, sInt4 *curTot, sInt4 *nd2x3,
405 short int *table50)
406 {
407 sInt4 sectLen; /* The length of the current section. */
408 sInt4 li_temp; /* A temporary holder of sInt4s. */
409
410 if ((sectNum == 2) || (sectNum == 3)) {
411 /* Figure out the size of section 2 and 3. */
412 /* ERO: check change from + 5 to +6+4 per r39022 */
413 if (*curTot > gribLen - (6 + 4)) {
414 errSprintf ("ERROR: Ran out of data in Section 2 or 3\n");
415 return -1;
416 }
417 /* Handle optional section 2. */
418 if (c_ipack[*curTot + 4] == 2) {
419 MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
420 if( sectLen < 0 || *curTot > INT_MAX - sectLen ) {
421 errSprintf ("ERROR: Invalid sectLen for section 2\n");
422 return -1;
423 }
424 *curTot = *curTot + sectLen;
425 if (ns[2] < sectLen)
426 ns[2] = sectLen;
427 /* ERO: check change from + 5 to +6+4 per r39022 */
428 if (*curTot + 6 + 4 > gribLen) {
429 errSprintf ("ERROR: Ran out of data in Section 3\n");
430 return -1;
431 }
432 }
433 /* Handle section 3. */
434 if (c_ipack[*curTot + 4] != 3) {
435 errSprintf ("ERROR: Section 3 labeled as %d\n",
436 c_ipack[*curTot + 4]);
437 return -2;
438 }
439 MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
440 if( sectLen < 0 || *curTot > INT_MAX - sectLen ) {
441 errSprintf ("ERROR: Invalid sectLen for section 3\n");
442 return -1;
443 }
444 if (ns[3] < sectLen)
445 ns[3] = sectLen;
446 /* While we are here, grab the total number of grid points nd2x3. */
447 MEMCPY_BIG (&li_temp, c_ipack + *curTot + 6, 4);
448 if (*nd2x3 < li_temp)
449 *nd2x3 = li_temp;
450 *curTot = *curTot + sectLen;
451 }
452 /*
453 #ifdef DEBUG
454 printf ("Section len (2=%ld) (3=%ld)\n", ns[2], ns[3]);
455 #endif
456 */
457
458 /* Figure out the size of section 4. */
459 if (*curTot > gribLen - 5) {
460 errSprintf ("ERROR: Ran out of data in Section 4\n");
461 return -1;
462 }
463 if (c_ipack[*curTot + 4] != 4) {
464 errSprintf ("ERROR: Section 4 labeled as %d\n", c_ipack[*curTot + 4]);
465 return -2;
466 }
467 MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
468 if( sectLen < 0 || *curTot > INT_MAX - sectLen ) {
469 errSprintf ("ERROR: Invalid sectLen for section 4\n");
470 return -1;
471 }
472 if (ns[4] < sectLen)
473 ns[4] = sectLen;
474 *curTot = *curTot + sectLen;
475 /*
476 #ifdef DEBUG
477 printf ("Section len (4=%ld < %ld)\n", sectLen, ns[4]);
478 #endif
479 */
480
481 /* Figure out the size of section 5. */
482 /* ERO: check change from + 5 to +9+2 per r39127 */
483 if (*curTot > gribLen - (9 + 2)) {
484 errSprintf ("ERROR: Ran out of data in Section 5\n");
485 return -1;
486 }
487 if (c_ipack[*curTot + 4] != 5) {
488 errSprintf ("ERROR: Section 5 labeled as %d\n", c_ipack[*curTot + 4]);
489 return -2;
490 }
491 MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
492 if( sectLen < 0 || *curTot > INT_MAX - sectLen ) {
493 errSprintf ("ERROR: Invalid sectLen for section 5\n");
494 return -1;
495 }
496 /* While we are here, grab the packing method. */
497 MEMCPY_BIG (table50, c_ipack + *curTot + 9, 2);
498 if (ns[5] < sectLen)
499 ns[5] = sectLen;
500 *curTot = *curTot + sectLen;
501 /*
502 #ifdef DEBUG
503 printf ("Section len (5=%ld < %ld)\n", sectLen, ns[5]);
504 #endif
505 */
506
507 /* Figure out the size of section 6. */
508 if (*curTot > gribLen - 5) {
509 errSprintf ("ERROR: Ran out of data in Section 6\n");
510 return -1;
511 }
512 if (c_ipack[*curTot + 4] != 6) {
513 errSprintf ("ERROR: Section 6 labeled as %d\n", c_ipack[*curTot + 4]);
514 return -2;
515 }
516 MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
517 if( sectLen < 0 || *curTot > INT_MAX - sectLen ) {
518 errSprintf ("ERROR: Invalid sectLen for section 6\n");
519 return -1;
520 }
521 if (ns[6] < sectLen)
522 ns[6] = sectLen;
523 *curTot = *curTot + sectLen;
524 /*
525 #ifdef DEBUG
526 printf ("Section len (6=%ld < %ld)\n", sectLen, ns[6]);
527 #endif
528 */
529
530 /* Figure out the size of section 7. */
531 if (*curTot > gribLen - 5) {
532 errSprintf ("ERROR: Ran out of data in Section 7\n");
533 return -1;
534 }
535 if (c_ipack[*curTot + 4] != 7) {
536 errSprintf ("ERROR: Section 7 labeled as %d\n", c_ipack[*curTot + 4]);
537 return -2;
538 }
539 MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
540 if( sectLen < 0 || *curTot > INT_MAX - sectLen ) {
541 errSprintf ("ERROR: Invalid sectLen for section 7\n");
542 return -1;
543 }
544 if (ns[7] < sectLen)
545 ns[7] = sectLen;
546 *curTot = *curTot + sectLen;
547 /*
548 #ifdef DEBUG
549 printf ("Section len (7=%ld < %ld)\n", sectLen, ns[7]);
550 #endif
551 */
552 return 0;
553 }
554
555 /*****************************************************************************
556 * FindSectLen() -- Review 12/2002
557 *
558 * Arthur Taylor / MDL
559 *
560 * PURPOSE
561 * Looks through a GRIB message and finds out how big each section is.
562 *
563 * ARGUMENTS
564 * c_ipack = The current GRIB2 message. (Input)
565 * gribLen = Length of c_ipack. (Input)
566 * ns = Array of section lengths. (Output)
567 * nd2x3 = Total number of grid points (Output)
568 * table50 = Type of packing used. (See code table 5.0) (GS5_SIMPLE = 0,
569 * GS5_CMPLX = 2, GS5_CMPLXSEC = 3) (Output)
570 *
571 * FILES/DATABASES: None
572 *
573 * RETURNS: int (could use errSprintf())
574 * 0 = OK
575 * -1 = Ran of data in a section
576 * -2 = Section not properly labeled.
577 *
578 * HISTORY
579 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
580 * 11/2002 AAT: Updated.
581 * 12/2002 (TK,AC,TB,&MS): Code Review.
582 * 3/2003 AAT: Made it handle multiple grids in the same GRIB2 message.
583 * 5/2003 AAT: Bug: Initialized size of section 2..6 to -1, instead
584 * of 2..7.
585 *
586 * NOTES
587 * 1) Assumes that the pack method of multiple grids are the same.
588 *****************************************************************************
589 */
FindSectLen(unsigned char * c_ipack,sInt4 gribLen,sInt4 ns[8],sInt4 * nd2x3,short int * table50)590 static int FindSectLen (unsigned char *c_ipack, sInt4 gribLen, sInt4 ns[8],
591 sInt4 *nd2x3, short int *table50)
592 {
593 sInt4 curTot; /* Where we are in the current GRIB message. */
594 char sectNum; /* Which section we are working with. */
595 int ans; /* The return error code of FindSectLen2to7. */
596 sInt4 sectLen; /* The length of the current section. */
597 int i; /* counter as we init ns[]. */
598
599 ns[0] = SECT0LEN_WORD * 4;
600 curTot = ns[0];
601
602 /* Figure out the size of section 1. */
603 if (curTot + 5 > gribLen) {
604 errSprintf ("ERROR: Ran out of data in Section 1\n");
605 return -1;
606 }
607 if (c_ipack[curTot + 4] != 1) {
608 errSprintf ("ERROR: Section 1 labeled as %d\n", c_ipack[curTot + 4]);
609 return -2;
610 }
611 MEMCPY_BIG (&(ns[1]), c_ipack + curTot, 4);
612 curTot += ns[1];
613 /*
614 #ifdef DEBUG
615 printf ("Section len (0=%ld) (1=%ld)\n", ns[0], ns[1]);
616 #endif
617 */
618 sectNum = 2;
619 for (i = 2; i < 8; i++) {
620 ns[i] = -1;
621 }
622 *nd2x3 = -1;
623 do {
624 if ((ans = FindSectLen2to7 (c_ipack, gribLen, ns, sectNum, &curTot,
625 nd2x3, table50)) != 0) {
626 return ans;
627 }
628 if( curTot + 4 > gribLen ) {
629 errSprintf ("ERROR: Ran out of data in Section 1\n");
630 return -1;
631 }
632 /* Try to read section 8. If it is "7777" == 926365495 regardless of
633 * endian'ness then we have a simple message, otherwise it is complex,
634 * and we need to read more. */
635 memcpy (§Len, c_ipack + curTot, 4);
636 if (sectLen == 926365495L) {
637 sectNum = 8;
638 } else {
639 if (curTot + 4 == gribLen) {
640 errSprintf("ERROR: Ran out of data in Section 1\n");
641 return -1;
642 }
643 sectNum = c_ipack[curTot + 4];
644 if ((sectNum < 2) || (sectNum > 7)) {
645 errSprintf ("ERROR (FindSectLen): Couldn't find the end of the "
646 "message\n");
647 errSprintf ("and it doesn't appear to repeat sections.\n");
648 errSprintf ("so it is probably an ASCII / binary bug\n");
649 errSprintf ("Max Sect Lengths: %ld %ld %ld %ld %ld %ld %ld"
650 " %ld\n", ns[0], ns[1], ns[2], ns[3], ns[4], ns[5],
651 ns[6], ns[7]);
652 return -2;
653 }
654 }
655 } while (sectNum != 8);
656 return 0;
657 }
658
659 /*****************************************************************************
660 * IS_Init() -- Review 12/2002
661 *
662 * Arthur Taylor / MDL
663 *
664 * PURPOSE
665 * Initialize the IS data structure. The IS_dataType is used to organize
666 * and allocate all the arrays that the unpack library uses.
667 * This makes an initial guess for the size of the arrays, and we use
668 * realloc to increase the size if needed. The write up: "UNPK_GRIB2
669 * 3/15/02" by Bryon Lawrence, Bob Glahn, David Rudack suggested these numbers
670 *
671 * ARGUMENTS
672 * is = The data structure to initialize. (Output)
673 *
674 * FILES/DATABASES: None
675 *
676 * RETURNS: void
677 *
678 * HISTORY
679 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
680 * 11/2002 AAT : Updated.
681 * 12/2002 (TK,AC,TB,&MS): Code Review.
682 *
683 * NOTES
684 * 1) Numbers not found in document were discussed with Bob Glahn on 8/29/2002
685 * 2) Possible exceptions:
686 * template 3.120 could need ns[3] = 1600
687 * template 4.30 could need a different ns4.
688 * 3) sizeof(float) == sizeof(sInt4), and unpacker does not use both ain
689 * and iain, so it is possible to have ain and iain point to the same
690 * memory. Not sure if this is safe, so I haven't done it.
691 *****************************************************************************
692 */
IS_Init(IS_dataType * is)693 void IS_Init (IS_dataType *is)
694 {
695 int i; /* A simple loop counter. */
696 is->ns[0] = 16;
697 is->ns[1] = 21;
698 is->ns[2] = 7;
699 is->ns[3] = 96;
700 is->ns[4] = 130; /* 60->130 in case there are some S4 time intervals */
701 is->ns[5] = 49;
702 is->ns[6] = 6;
703 is->ns[7] = 8;
704 for (i = 0; i < 8; i++) {
705 is->is[i] = (sInt4 *) calloc (is->ns[i], sizeof (sInt4));
706 }
707 /* Allocate grid memory. */
708 is->nd2x3 = 0;
709 is->iain = nullptr;
710 is->ib = nullptr;
711 /* Allocate section 2 int memory. */
712 is->nidat = 0;
713 is->idat = nullptr;
714 /* Allocate section 2 float memory. */
715 is->nrdat = 0;
716 is->rdat = nullptr;
717 /* Allocate storage for ipack. */
718 is->ipackLen = 0;
719 is->ipack = nullptr;
720 }
721
722 /*****************************************************************************
723 * IS_Free() -- Review 12/2002
724 *
725 * Arthur Taylor / MDL
726 *
727 * PURPOSE
728 * Free the memory allocated in the IS data structure.
729 * The IS_dataType is used to organize and allocate all the arrays that the
730 * unpack library uses.
731 *
732 * ARGUMENTS
733 * is = The data structure to Free. (Output)
734 *
735 * FILES/DATABASES: None
736 *
737 * RETURNS: void
738 *
739 * HISTORY
740 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
741 * 11/2002 AAT : Updated.
742 * 12/2002 (TK,AC,TB,&MS): Code Review.
743 *
744 * NOTES
745 *****************************************************************************
746 */
IS_Free(IS_dataType * is)747 void IS_Free (IS_dataType *is)
748 {
749 int i; /* A simple loop counter. */
750 for (i = 0; i < 8; i++) {
751 free (is->is[i]);
752 is->is[i] = nullptr;
753 is->ns[i] = 0;
754 }
755 /* Free grid memory. */
756 free (is->iain);
757 is->iain = nullptr;
758 free (is->ib);
759 is->ib = nullptr;
760 is->nd2x3 = 0;
761 /* Free section 2 int memory. */
762 free (is->idat);
763 is->idat = nullptr;
764 is->nidat = 0;
765 /* Free section 2 float memory. */
766 free (is->rdat);
767 is->rdat = nullptr;
768 is->nrdat = 0;
769 /* Free storage for ipack. */
770 free (is->ipack);
771 is->ipack = nullptr;
772 is->ipackLen = 0;
773 }
774
775 /*****************************************************************************
776 * ReadGrib2Record() -- Review 12/2002
777 *
778 * Arthur Taylor / MDL
779 *
780 * PURPOSE
781 * Reads a GRIB2 message from a file which is already opened and is pointing
782 * at the correct message. It reads in the message storing the results in
783 * Grib_Data which is of size grib_DataLen. If needed, it increases
784 * grib_DataLen enough to fit the current message's grid. It converts (if
785 * appropriate) the data in Grib_Data to the units specified in f_unit.
786 *
787 * In addition it updates offset, and stores the meta data returned by the
788 * unpacker library in both IS, and (after parsing it) in meta.
789 *
790 * Note: It expects meta and IS to already be initialized through calls to
791 * MetaInit(&meta) and IS_Init(&is) respectively.
792 *
793 * ARGUMENTS
794 * fp = An opened GRIB2 file already at the correct message. (Input)
795 * fileLen = Length of the opened file. (Input)
796 * f_unit = 0 use GRIB2 units, 1 use English, 2 use metric. (Input)
797 * Grib_Data = The read in GRIB2 grid. (Output)
798 * grib_DataLen = Size of Grib_Data. (Output)
799 * meta = A filled in meta structure (Output)
800 * IS = The structure containing all the arrays that the
801 * unpacker uses (Output)
802 * subgNum = Which subgrid in the GRIB2 message is of interest.
803 * (0 = first grid), if it can't find message subgNum,
804 * returns -5, and an error message (Input)
805 * majEarth = Used to override the GRIB major axis of earth. (Input)
806 * minEarth = Used to override the GRIB minor axis of earth. (Input)
807 * simpVer = The version of the simple weather code to use when parsing
808 * the WxString. (Input)
809 * f_endMsg = 1 means we finished reading the previous GRIB message, or
810 * there was no previous GRIB message. 0 means that we need
811 * to read a subgrid of the previous message. (Input/Output)
812 * lwlt, uprt = If the lat is not -100, then lwlt, and uprt define a
813 * subgrid that the user is interested in. Get the map
814 * projection out of the GRIB2 message, and do everything
815 * on the subgrid. (if lwlt, and uprt are not "correct", the
816 * lat/lons may get swapped) (Input/Output)
817 *
818 * FILES/DATABASES:
819 * An already opened "GRIB2" File
820 *
821 * RETURNS: int (could use errSprintf())
822 * 0 = OK
823 * -1 = Problems in section 0
824 * -2 = Problems figuring out the Section Lengths.
825 * -3 = Error returned by unpack library.
826 * -4 = Problems parsing the Meta Data.
827 *
828 * HISTORY
829 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
830 * 11/2002 AAT: Updated.
831 * 12/2002 (TK,AC,TB,&MS): Code Review.
832 * 1/2003 AAT: It was not error coded 208, but rather 202 to look for.
833 * 3/2003 AAT: Modified handling of section 2 stuff (no loop)
834 * 3/2003 AAT: Added ability to handle multiple grids in same message.
835 * 4/2003 AAT: Added ability to call GRIB1 decoder for GRIB1 messages.
836 * 5/2003 AAT: Update the offset for ReadGrib1.
837 * 6/2003 Matthew T. Kallio (matt@wunderground.com):
838 * "wmo" dimension increased to WMO_HEADER_LEN + 1 (for '\0' char)
839 * 7/2003 AAT: switched to checking against element name for Wx instead
840 * of pds2.sect2.ptrType == GS2_WXTYPE
841 * 7/2003 AAT: Allowed user to override the radius of earth.
842 * 8/2003 AAT: Removed dependence on fileLen and offset.
843 * 2/2004 AAT: Added "f_endMsg" logic.
844 * 2/2004 AAT: Added subgrid potential.
845 * 2/2004 AAT: Added maj/min Earth override.
846 *
847 * NOTES
848 * 1) Reason ns[7] is not MAX (IS.ns[], local_ns[]) is because local_ns[7]
849 * is size of the packed message, but ns[7] refers to the returned meta
850 * data which unpacker library found in section 7, which is a lot smaller.
851 * 2) Problem: MDL's sect2 is packed and we have no idea how large it is
852 * when unpacked. So we allocate room for 4000 sInt4s and 500 floats.
853 * We then check 'jer' for error "202", if we find it we double the size
854 * and call the unpacker again.
855 * 3/26/2003: Changed this to be: try once with size
856 * = max (32 * packed size, 4000)
857 * Should be fewer calls (more memory intensive) same result, since we had
858 * been doubling it 5 times.
859 * 3) For Complex second order packing (GS5_CMPLXSEC) the unpacker needs nd5
860 * (which is size of message) to be >= nd2x3 (size of grid).
861 * 3a) Appears to also need this if simple packing, and has a bitmap.
862 * 4) inew = 1: Currently we only expect 1 grid in 1 GRIB message, although
863 * the library does allow for multiple grids in a GRIB message.
864 * 5) iclean = 1: This only maters if there is bitmap data, otherwise it is
865 * ignored. For bitmap data, if == 0, it embeds the given values for
866 * xmissp, and xmisss. We don't embed because we don't know what to set
867 * xmissp or xmisss to. Instead after we know the range, we choose a value
868 * and walk through the bitmap setting grib_Data appropriately.
869 * 5a) iclean = 0; This is because we do want the missing values embedded.
870 * that is we want the missing values to be place holders.
871 * 6) f_endMsg is true if in the past we either completed reading a message,
872 * or we haven't read any messages. In either case we need to read the
873 * next message from file. If f_endMsg is false, then there is more to read
874 * from IS->ipack, so we don't want to throw it out, nor have to re-read
875 * ipack from disk.
876 *
877 * Question: Should we double ns[2] when we double nrdat, and nidat?
878 *****************************************************************************
879 */
880 #define SECT2_INIT_SIZE 4000
881 #define UNPK_NUM_ERRORS 22
ReadGrib2Record(VSILFILE * fp,sChar f_unit,double ** Grib_Data,uInt4 * grib_DataLen,grib_MetaData * meta,IS_dataType * IS,int subgNum,double majEarth,double minEarth,int simpVer,int simpWWA,sInt4 * f_endMsg,CPL_UNUSED LatLon * lwlf,CPL_UNUSED LatLon * uprt)882 int ReadGrib2Record (VSILFILE *fp, sChar f_unit, double **Grib_Data,
883 uInt4 *grib_DataLen, grib_MetaData *meta,
884 IS_dataType *IS, int subgNum, double majEarth,
885 double minEarth, int simpVer, int simpWWA,
886 sInt4 *f_endMsg,
887 CPL_UNUSED LatLon *lwlf,
888 CPL_UNUSED LatLon *uprt)
889 {
890 sInt4 l3264b; /* Number of bits in a sInt4. Needed by FORTRAN
891 * unpack library to determine if system has a 4
892 * byte_ sInt4 or an 8 byte sInt4. */
893 char *buff = nullptr; /* Holds the info between records. */
894 uInt4 buffLen; /* Length of info between records. */
895 sInt4 sect0[SECT0LEN_WORD]; /* Holds the current Section 0. */
896 uInt4 gribLen; /* Length of the current GRIB message. */
897 sInt4 nd5; /* Size of grib message rounded up to the nearest
898 * sInt4. */
899 /* A char ptr to the message stored in IS->ipack */
900 unsigned char *c_ipack = nullptr;
901 sInt4 local_ns[8]; /* Local copy of section lengths. */
902 sInt4 nd2x3; /* Total number of grid points. */
903 short int table50; /* Type of packing used. (See code table 5.0)
904 * (GS5_SIMPLE==0, GS5_CMPLX==2, GS5_CMPLXSEC==3) */
905 sInt4 nidat; /* Size of section 2 if it contains integer data. */
906 sInt4 nrdat; /* Size of section 2 if it contains float data. */
907 sInt4 inew; /* 1 if this is the first grid we are reading. 0 if
908 * this is the second or later grid from the same
909 * GRIB message. */
910 sInt4 iclean = 0; /* 0 embed the missing values, 1 don't. */
911 int j; /* Counter used to find the desired subgrid. */
912 sInt4 kfildo = 5; /* FORTRAN Unit number for diagnostic info. Ignored,
913 * unless library is compiled a particular way. */
914 sInt4 ibitmap = 0; /* 0 means no bitmap returned, otherwise 1. */
915 float xmissp = 0.0f; /* The primary missing value. If iclean = 0, this
916 * value is embedded in grid, otherwise it is the
917 * value returned from the GRIB message. */
918 float xmisss = 0.0f; /* The secondary missing value. If iclean = 0, this
919 * value is embedded in grid, otherwise it is the
920 * value returned from the GRIB message. */
921 sInt4 jer[UNPK_NUM_ERRORS * 2]; /* Any Error codes along with their *
922 * severity levels generated using the *
923 * unpack GRIB2 library. */
924 sInt4 ndjer = UNPK_NUM_ERRORS; /* The number of rows in JER( ). */
925 sInt4 kjer; /* The actual number of errors returned in JER. */
926 size_t i; /* counter as we loop through jer. */
927 double unitM, unitB; /* values in y = m x + b used for unit conversion. */
928 char unitName[15]; /* Holds the string name of the current unit. */
929 int unitLen; /* String length of string name of current unit. */
930 int version; /* Which version of GRIB is in this message. */
931 sInt4 cnt; /* Used to help compact the weather table. */
932 //gdsType newGds; /* The GDS of the subgrid if needed. */
933 int x1, y1; /* The original grid coordinates of the lower left
934 * corner of the subgrid. */
935 int x2, y2; /* The original grid coordinates of the upper right
936 * corner of the subgrid. */
937 uChar f_subGrid; /* True if we have a subgrid. */
938 sInt4 Nx, Ny; /* original size of the data. */
939
940 /*
941 * f_endMsg is 1 if in the past we either completed reading a message,
942 * or we haven't read any messages. In either case we need to read the
943 * next message from file.
944 * If f_endMsg is false, then there is more to read from IS->ipack, so we
945 * don't want to throw it out, nor have to re-read ipack from disk.
946 */
947 l3264b = sizeof (sInt4) * 8;
948 buffLen = 0;
949 if (*f_endMsg == 1) {
950 if (ReadSECT0 (fp, &buff, &buffLen, -1, sect0, &gribLen, &version) < 0) {
951 preErrSprintf ("Inside ReadGrib2Record\n");
952 free (buff);
953 return -1;
954 }
955 meta->GribVersion = version;
956 if (version == 1) {
957 if (ReadGrib1Record (fp, f_unit, Grib_Data, grib_DataLen, meta, IS,
958 sect0, gribLen, majEarth, minEarth) != 0) {
959 preErrSprintf ("Problems with ReadGrib1Record called by "
960 "ReadGrib2Record\n");
961 free (buff);
962 return -1;
963 }
964 *f_endMsg = 1;
965 free (buff);
966 return 0;
967 }
968 #if 0
969 /* tdlpack is no longer supported by GDAL */
970 else if (version == -1) {
971 if (ReadTDLPRecord (fp, Grib_Data, grib_DataLen, meta, IS,
972 sect0, gribLen, majEarth, minEarth) != 0) {
973 preErrSprintf ("Problems with ReadGrib1Record called by "
974 "ReadGrib2Record\n");
975 free (buff);
976 return -1;
977 }
978 free (buff);
979 return 0;
980 }
981 #endif
982
983 /*
984 * Make room for entire message, and read it in.
985 */
986 /* nd5 needs to be gribLen in (sInt4) units rounded up. */
987 if( gribLen > UINT_MAX - 3 )
988 {
989 errSprintf("Invalid value of gribLen");
990 free (buff);
991 return -1;
992 }
993 nd5 = (gribLen + 3) / 4;
994 if (nd5 > IS->ipackLen) {
995 if( gribLen > 100 * 1024 * 1024 )
996 {
997 vsi_l_offset curPos = VSIFTellL(fp);
998 VSIFSeekL(fp, 0, SEEK_END);
999 vsi_l_offset fileSize = VSIFTellL(fp);
1000 VSIFSeekL(fp, curPos, SEEK_SET);
1001 if( fileSize < gribLen )
1002 {
1003 errSprintf("File too short");
1004 free (buff);
1005 return -1;
1006 }
1007 }
1008 sInt4* ipackNew = (sInt4 *) realloc ((void *) (IS->ipack),
1009 nd5 * sizeof (sInt4));
1010 if( ipackNew == nullptr )
1011 {
1012 errSprintf("Out of memory");
1013 free (buff);
1014 return -1;
1015 }
1016 IS->ipackLen = nd5;
1017 IS->ipack = ipackNew;
1018 }
1019 c_ipack = (unsigned char *) IS->ipack;
1020 /* Init last sInt4 to 0, to make sure that the padded bytes are 0. */
1021 IS->ipack[nd5 - 1] = 0;
1022 /* Init first 4 sInt4 to sect0. */
1023 memcpy (c_ipack, sect0, SECT0LEN_WORD * 4);
1024 /* Read in the rest of the message. */
1025 if (VSIFReadL (c_ipack + SECT0LEN_WORD * 4, sizeof (char),
1026 (gribLen - SECT0LEN_WORD * 4), fp) != (gribLen - SECT0LEN_WORD * 4)) {
1027 errSprintf ("GribLen = %ld, SECT0Len_WORD = %d\n", gribLen,
1028 SECT0LEN_WORD);
1029 errSprintf ("Ran out of file\n");
1030 free (buff);
1031 return -1;
1032 }
1033
1034 /*
1035 * Make sure the arrays are large enough for call to unpacker library.
1036 */
1037 /* FindSectLen Does not want (ipack / c_ipack) word swapped, because
1038 * that would make it much more confusing to find bytes in c_ipack. */
1039 if (FindSectLen (c_ipack, gribLen, local_ns, &nd2x3, &table50) < 0) {
1040 preErrSprintf ("Inside ReadGrib2Record.. Calling FindSectLen\n");
1041 free (buff);
1042 return -2;
1043 }
1044 if( nd2x3 < 0 || nd2x3 > INT_MAX / static_cast<int>(sizeof (sInt4)) )
1045 {
1046 preErrSprintf ("Invalid nd2x3\n");
1047 free (buff);
1048 return -2;
1049 }
1050
1051 /* Make sure all 'is' arrays except ns[7] are MAX (IS.ns[] ,
1052 * local_ns[]). See note 1 for reason to exclude ns[7] from MAX (). */
1053 for (i = 0; i < 7; i++) {
1054 if (local_ns[i] > IS->ns[i]) {
1055 IS->ns[i] = local_ns[i];
1056 IS->is[i] = (sInt4 *) realloc ((void *) (IS->is[i]),
1057 IS->ns[i] * sizeof (sInt4));
1058 }
1059 }
1060
1061 /* Allocate room for sect 2. If local_ns[2] = -1 there is no sect 2. */
1062 if (local_ns[2] == -1) {
1063 nidat = 10;
1064 nrdat = 10;
1065 } else {
1066 /*
1067 * See note 2) We have a section 2, so use:
1068 * MAX (32 * local_ns[2],SECT2_INTSIZE)
1069 * and MAX (32 * local_ns[2],SECT2_FLOATSIZE)
1070 * for size of section 2 unpacked.
1071 */
1072 nidat = (32 * local_ns[2] < SECT2_INIT_SIZE) ? SECT2_INIT_SIZE :
1073 32 * local_ns[2];
1074 nrdat = nidat;
1075 }
1076 if (nidat > IS->nidat) {
1077 IS->nidat = nidat;
1078 IS->idat = (sInt4 *) realloc ((void *) IS->idat,
1079 IS->nidat * sizeof (sInt4));
1080 }
1081 if (nrdat > IS->nrdat) {
1082 IS->nrdat = nrdat;
1083 IS->rdat = (float *) realloc ((void *) IS->rdat,
1084 IS->nrdat * sizeof (float));
1085 }
1086 /* Make sure we have room for the GRID part of the output. */
1087 if (nd2x3 > IS->nd2x3) {
1088 if( nd2x3 > 100 * 1024 * 1024 )
1089 {
1090
1091 vsi_l_offset curPos = VSIFTellL(fp);
1092 VSIFSeekL(fp, 0, SEEK_END);
1093 vsi_l_offset fileSize = VSIFTellL(fp);
1094 VSIFSeekL(fp, curPos, SEEK_SET);
1095 // allow a compression ratio of 1:1000
1096 if( (vsi_l_offset)(nd2x3 / 1000) > fileSize )
1097 {
1098 preErrSprintf ("File too short\n");
1099 free (buff);
1100 return -2;
1101 }
1102 }
1103
1104 IS->nd2x3 = nd2x3;
1105 if( Grib_Data )
1106 {
1107 IS->iain = (sInt4 *) realloc ((void *) IS->iain,
1108 IS->nd2x3 * sizeof (sInt4));
1109 }
1110 IS->ib = (sInt4 *) realloc ((void *) IS->ib,
1111 IS->nd2x3 * sizeof (sInt4));
1112 }
1113 /* See note 3) If table50 == 3, unpacker library needs nd5 >= nd2x3. */
1114 if ((table50 == 3) || (table50 == 0)) {
1115 if (nd5 < nd2x3) {
1116 nd5 = nd2x3;
1117 if (nd5 > IS->ipackLen) {
1118 sInt4* ipackNew = (sInt4 *) realloc ((void *) (IS->ipack),
1119 nd5 * sizeof (sInt4));
1120 if( ipackNew == nullptr )
1121 {
1122 errSprintf("Out of memory");
1123 free (buff);
1124 return -1;
1125 }
1126 // zero initialize to avoid later use of undefined values
1127 // Not sure if it is strictly needed (perhaps some later
1128 // steps should detect errors / data shortening), but
1129 // cannot hurt.
1130 // see https://trac.osgeo.org/gdal/ticket/6967
1131 memset(ipackNew + IS->ipackLen, 0,
1132 (nd5 - IS->ipackLen) * sizeof(sInt4));
1133 IS->ipackLen = nd5;
1134 IS->ipack = ipackNew;
1135 }
1136 /* Don't need to do the following, but we do in case code
1137 * changes. */
1138 c_ipack = (unsigned char *) IS->ipack;
1139 }
1140 }
1141 IS->nd5 = nd5;
1142 /* Unpacker library requires ipack to be MSB. */
1143 /*
1144 #ifdef DEBUG
1145 if (1==1) {
1146 FILE *fp = fopen ("test.bin", "wb");
1147 fwrite (IS->ipack, sizeof (sInt4), IS->nd5, fp);
1148 fclose (fp);
1149 }
1150 #endif
1151 */
1152 /*
1153 #ifdef LITTLE_ENDIAN
1154 memswp (IS->ipack, sizeof (sInt4), IS->nd5);
1155 #endif
1156 */
1157 } else {
1158 c_ipack = (unsigned char *) IS->ipack;
1159 /* GRIB2 files are in big endian so c_ipack is as well. */
1160 #ifdef LITTLE_ENDIAN
1161 revmemcpy (&gribLen, &(c_ipack[12]), sizeof (sInt4));
1162 #else
1163 memcpy (&gribLen, &(c_ipack[12]), sizeof (sInt4));
1164 #endif
1165 }
1166 free (buff);
1167
1168 /* Loop through the grib message looking for the subgNum grid. subgNum
1169 * goes from 0 to n-1. */
1170 for (j = 0; j <= subgNum; j++) {
1171 if (j == 0) {
1172 inew = 1;
1173 } else {
1174 inew = 0;
1175 }
1176
1177 /* Note we are getting data back either as a float or an int, but not
1178 * both, so we don't need to allocated room for both. */
1179 unpk_g2ncep (&kfildo,
1180 j == subgNum ? (float *) (IS->iain) : nullptr,
1181 j == subgNum ? IS->iain : nullptr,
1182 &(IS->nd2x3),
1183 IS->idat, &(IS->nidat), IS->rdat, &(IS->nrdat), IS->is[0],
1184 &(IS->ns[0]), IS->is[1], &(IS->ns[1]), IS->is[2],
1185 &(IS->ns[2]), IS->is[3], &(IS->ns[3]), IS->is[4],
1186 &(IS->ns[4]), IS->is[5], &(IS->ns[5]), IS->is[6],
1187 &(IS->ns[6]), IS->is[7], &(IS->ns[7]), IS->ib, &ibitmap,
1188 c_ipack, &(IS->nd5), &xmissp, &xmisss, &inew, &iclean,
1189 &l3264b, f_endMsg, jer, &ndjer, &kjer);
1190 /*
1191 unpk_grib2 (&kfildo, (float *) (IS->iain), IS->iain, &(IS->nd2x3),
1192 IS->idat, &(IS->nidat), IS->rdat, &(IS->nrdat), IS->is[0],
1193 &(IS->ns[0]), IS->is[1], &(IS->ns[1]), IS->is[2],
1194 &(IS->ns[2]), IS->is[3], &(IS->ns[3]), IS->is[4],
1195 &(IS->ns[4]), IS->is[5], &(IS->ns[5]), IS->is[6],
1196 &(IS->ns[6]), IS->is[7], &(IS->ns[7]), IS->ib, &ibitmap,
1197 IS->ipack, &(IS->nd5), &xmissp, &xmisss, &inew, &iclean,
1198 &l3264b, f_endMsg, jer, &ndjer, &kjer);
1199 */
1200 /*
1201 * Check for error messages...
1202 * If we get an error message, print it, and return.
1203 */
1204 for (i = 0; i < (uInt4) kjer; i++) {
1205 if (jer[ndjer + i] == 0) {
1206 /* no error. */
1207 } else if (jer[ndjer + i] == 1) {
1208 /* Warning. */
1209 #ifdef DEBUG
1210 printf ("Warning: Unpack library warning code (%d %d)\n",
1211 jer[i], jer[ndjer + i]);
1212 #endif
1213 } else {
1214 /* BAD Error. */
1215 errSprintf ("ERROR: Unpack library error code (%ld %ld)\n",
1216 jer[i], jer[ndjer + i]);
1217 return -3;
1218 }
1219 }
1220 }
1221
1222 /* Parse the meta data out. */
1223 if (MetaParse (meta, IS->is[0], IS->ns[0], IS->is[1], IS->ns[1],
1224 IS->is[2], IS->ns[2], IS->rdat, IS->nrdat, IS->idat,
1225 IS->nidat, IS->is[3], IS->ns[3], IS->is[4], IS->ns[4],
1226 IS->is[5], IS->ns[5], gribLen, xmissp, xmisss, simpVer, simpWWA)
1227 != 0) {
1228 #ifdef DEBUG
1229 FILE *l_fp;
1230 if ((l_fp = fopen ("dump.is0", "wt")) != nullptr) {
1231 for (i = 0; i < 8; i++) {
1232 fprintf (l_fp, "---Section %d---\n", (int) i);
1233 for (j = 1; j <= IS->ns[i]; j++) {
1234 fprintf (l_fp, "IS%d Item %d = %d\n", (int) i, (int) j, IS->is[i][j - 1]);
1235 }
1236 }
1237 fclose (l_fp);
1238 }
1239 #endif
1240 preErrSprintf ("Inside ReadGrib2Record.. Problems in MetaParse\n");
1241 return -4;
1242 }
1243
1244 if ((majEarth > 6000) && (majEarth < 7000)) {
1245 if ((minEarth > 6000) && (minEarth < 7000)) {
1246 meta->gds.f_sphere = 0;
1247 meta->gds.majEarth = majEarth;
1248 meta->gds.minEarth = minEarth;
1249 } else {
1250 meta->gds.f_sphere = 1;
1251 meta->gds.majEarth = majEarth;
1252 meta->gds.minEarth = majEarth;
1253 }
1254 }
1255
1256 /* Figure out an equation to pass to ParseGrid to convert the units for
1257 * this grid. */
1258 /*
1259 if (ComputeUnit (meta->pds2.prodType, meta->pds2.sect4.templat,
1260 meta->pds2.sect4.cat, meta->pds2.sect4.subcat, f_unit,
1261 &unitM, &unitB, unitName) == 0) {
1262 */
1263 if (ComputeUnit (meta->convert, meta->unitName, f_unit, &unitM, &unitB,
1264 unitName) == 0) {
1265 unitLen = static_cast<int>(strlen (unitName));
1266 meta->unitName = (char *) realloc ((void *) (meta->unitName),
1267 1 + unitLen * sizeof (char));
1268 memcpy (meta->unitName, unitName, unitLen);
1269 meta->unitName[unitLen] = '\0';
1270 }
1271
1272 /* compute the subgrid. */
1273 /*
1274 if ((lwlf->lat != -100) && (uprt->lat != -100)) {
1275 Nx = meta->gds.Nx;
1276 Ny = meta->gds.Ny;
1277 if (computeSubGrid (lwlf, &x1, &y1, uprt, &x2, &y2, &(meta->gds),
1278 &newGds) != 0) {
1279 preErrSprintf ("ERROR: In compute subgrid.\n");
1280 return 1;
1281 }
1282 // I couldn't decide if I should "permanently" change the GDS or not.
1283 // when I wrote computeSubGrid. If next line stays, really should
1284 // rewrite computeSubGrid.
1285 memcpy (&(meta->gds), &newGds, sizeof (gdsType));
1286 f_subGrid = 1;
1287 } else {
1288 Nx = meta->gds.Nx;
1289 Ny = meta->gds.Ny;
1290 x1 = 1;
1291 x2 = Nx;
1292 y1 = 1;
1293 y2 = Ny;
1294 f_subGrid = 0;
1295 }
1296 */
1297 if( meta->gds.Nx == 0 || meta->gds.Nx > INT_MAX ||
1298 meta->gds.Ny == 0 || meta->gds.Ny > INT_MAX )
1299 {
1300 errSprintf ("Invalid grid dimension.\n");
1301 return -3;
1302 }
1303 Nx = meta->gds.Nx;
1304 Ny = meta->gds.Ny;
1305 x1 = 1;
1306 x2 = Nx;
1307 y1 = 1;
1308 y2 = Ny;
1309 f_subGrid = 0;
1310
1311 #ifdef deadcode
1312 /* Figure out if we need iain or ain, and set it to Grib_Data. At the
1313 * same time handle any bitmaps, and compute some statistics. */
1314 if ((f_subGrid) && (meta->gds.scan != 64)) {
1315 errSprintf ("Can not do a subgrid of non scanmode 64 grid yet.\n");
1316 return -3;
1317 }
1318 #endif
1319
1320 if (Grib_Data != nullptr && strcmp (meta->element, "Wx") != 0) {
1321 if (strcmp (meta->element, "WWA") != 0) {
1322 ParseGrid (fp, &(meta->gridAttrib), Grib_Data, grib_DataLen, Nx, Ny,
1323 meta->gds.scan, IS->nd2x3, IS->iain, ibitmap, IS->ib, unitM, unitB, 0,
1324 0, nullptr, f_subGrid, x1, y1, x2, y2);
1325 } else {
1326 ParseGrid (fp, &(meta->gridAttrib), Grib_Data, grib_DataLen, Nx, Ny,
1327 meta->gds.scan, IS->nd2x3, IS->iain, ibitmap, IS->ib, unitM, unitB, 1,
1328 meta->pds2.sect2.hazard.dataLen, meta->pds2.sect2.hazard.f_valid, f_subGrid, x1, y1,
1329 x2, y2);
1330 /* compact the table to only those which are actually used. */
1331 cnt = 0;
1332 for (i = 0; i < meta->pds2.sect2.hazard.dataLen; i++) {
1333 if (meta->pds2.sect2.hazard.f_valid[i] == 2) {
1334 meta->pds2.sect2.hazard.haz[i].validIndex = cnt;
1335 cnt++;
1336 } else if (meta->pds2.sect2.hazard.f_valid[i] == 3) {
1337 meta->pds2.sect2.hazard.f_valid[i] = 0;
1338 meta->pds2.sect2.hazard.haz[i].validIndex = cnt;
1339 cnt++;
1340 } else {
1341 meta->pds2.sect2.hazard.haz[i].validIndex = -1;
1342 }
1343 }
1344 }
1345 } else if( Grib_Data != nullptr ) {
1346 /* Handle weather grid. ParseGrid looks up the values... If they are
1347 * "<Invalid>" it sets it to missing (or creates one). If the table
1348 * entry is used it sets f_valid to 2. */
1349 ParseGrid (fp, &(meta->gridAttrib), Grib_Data, grib_DataLen, Nx, Ny,
1350 meta->gds.scan, IS->nd2x3, IS->iain, ibitmap, IS->ib, unitM, unitB, 1,
1351 meta->pds2.sect2.wx.dataLen, meta->pds2.sect2.wx.f_valid, f_subGrid, x1, y1,
1352 x2, y2);
1353 if( *Grib_Data == nullptr )
1354 return -1;
1355
1356 /* compact the table to only those which are actually used. */
1357 cnt = 0;
1358 for (i = 0; i < meta->pds2.sect2.wx.dataLen; i++) {
1359 if (meta->pds2.sect2.wx.f_valid[i] == 2) {
1360 meta->pds2.sect2.wx.ugly[i].validIndex = cnt;
1361 cnt++;
1362 } else if (meta->pds2.sect2.wx.f_valid[i] == 3) {
1363 meta->pds2.sect2.wx.f_valid[i] = 0;
1364 meta->pds2.sect2.wx.ugly[i].validIndex = cnt;
1365 cnt++;
1366 } else {
1367 meta->pds2.sect2.wx.ugly[i].validIndex = -1;
1368 }
1369 }
1370 }
1371
1372 /* Figure out some other non-section oriented meta data. */
1373 /* strftime (meta->refTime, 20, "%Y%m%d%H%M",
1374 gmtime (&(meta->pds2.refTime)));
1375 */
1376 Clock_Print (meta->refTime, 20, meta->pds2.refTime, "%Y%m%d%H%M", 0);
1377 /*
1378 strftime (meta->validTime, 20, "%Y%m%d%H%M",
1379 gmtime (&(meta->pds2.sect4.validTime)));
1380 */
1381 Clock_Print (meta->validTime, 20, meta->pds2.sect4.validTime,
1382 "%Y%m%d%H%M", 0);
1383
1384 const double deltTime = meta->pds2.sect4.validTime - meta->pds2.refTime;
1385 if (deltTime >= std::numeric_limits<sInt4>::max() ||
1386 deltTime <= std::numeric_limits<sInt4>::min()) {
1387 meta->deltTime = 0;
1388 preErrSprintf ("deltTime over range\n");
1389 return -4;
1390 }
1391 meta->deltTime = static_cast<sInt4>(deltTime);
1392
1393 return 0;
1394 }
1395
1396 #if 0 // unused by GDAL
1397 int ReadGrib2RecordFast (FILE *fp, sChar f_unit, double **Grib_Data,
1398 uInt4 *grib_DataLen, grib_MetaData *meta,
1399 IS_dataType *IS, int subgNum, double majEarth,
1400 double minEarth, int simpVer, int simpWWA,
1401 sInt4 *f_endMsg, LatLon *lwlf, LatLon *uprt)
1402 {
1403 sInt4 l3264b; /* Number of bits in a sInt4. Needed by FORTRAN
1404 * unpack library to determine if system has a 4
1405 * byte_ sInt4 or an 8 byte sInt4. */
1406 char *buff; /* Holds the info between records. */
1407 uInt4 buffLen; /* Length of info between records. */
1408 sInt4 sect0[SECT0LEN_WORD]; /* Holds the current Section 0. */
1409 uInt4 gribLen; /* Length of the current GRIB message. */
1410 sInt4 nd5; /* Size of grib message rounded up to the nearest
1411 * sInt4. */
1412 unsigned char *c_ipack; /* A char ptr to the message stored in IS->ipack */
1413 sInt4 local_ns[8]; /* Local copy of section lengths. */
1414 sInt4 nd2x3; /* Total number of grid points. */
1415 short int table50; /* Type of packing used. (See code table 5.0)
1416 * (GS5_SIMPLE==0, GS5_CMPLX==2, GS5_CMPLXSEC==3) */
1417 sInt4 nidat; /* Size of section 2 if it contains integer data. */
1418 sInt4 nrdat; /* Size of section 2 if it contains float data. */
1419 sInt4 inew; /* 1 if this is the first grid we are reading. 0 if
1420 * this is the second or later grid from the same
1421 * GRIB message. */
1422 sInt4 iclean = 0; /* 0 embed the missing values, 1 don't. */
1423 int j; /* Counter used to find the desired subgrid. */
1424 sInt4 kfildo = 5; /* FORTRAN Unit number for diagnostic info. Ignored,
1425 * unless library is compiled a particular way. */
1426 sInt4 ibitmap; /* 0 means no bitmap returned, otherwise 1. */
1427 float xmissp; /* The primary missing value. If iclean = 0, this
1428 * value is embedded in grid, otherwise it is the
1429 * value returned from the GRIB message. */
1430 float xmisss; /* The secondary missing value. If iclean = 0, this
1431 * value is embedded in grid, otherwise it is the
1432 * value returned from the GRIB message. */
1433 sInt4 jer[UNPK_NUM_ERRORS * 2]; /* Any Error codes along with their *
1434 * severity levels generated using the *
1435 * unpack GRIB2 library. */
1436 sInt4 ndjer = UNPK_NUM_ERRORS; /* The number of rows in JER( ). */
1437 sInt4 kjer; /* The actual number of errors returned in JER. */
1438 size_t i; /* counter as we loop through jer. */
1439 double unitM, unitB; /* values in y = m x + b used for unit conversion. */
1440 char unitName[15]; /* Holds the string name of the current unit. */
1441 int unitLen; /* String length of string name of current unit. */
1442 int version; /* Which version of GRIB is in this message. */
1443 gdsType newGds; /* The GDS of the subgrid if needed. */
1444 int x1, y1; /* The original grid coordinates of the lower left
1445 * corner of the subgrid. */
1446 int x2, y2; /* The original grid coordinates of the upper right
1447 * corner of the subgrid. */
1448 uChar f_subGrid; /* True if we have a subgrid. */
1449 sInt4 Nx, Ny; /* original size of the data. */
1450
1451 /*
1452 * f_endMsg is 1 if in the past we either completed reading a message,
1453 * or we haven't read any messages. In either case we need to read the
1454 * next message from file.
1455 * If f_endMsg is false, then there is more to read from IS->ipack, so we
1456 * don't want to throw it out, nor have to re-read ipack from disk.
1457 */
1458 l3264b = sizeof (sInt4) * 8;
1459 buff = NULL;
1460 buffLen = 0;
1461 if (*f_endMsg == 1) {
1462 if (ReadSECT0 (fp, &buff, &buffLen, -1, sect0, &gribLen, &version) < 0) {
1463 preErrSprintf ("Inside ReadGrib2Record\n");
1464 free (buff);
1465 return -1;
1466 }
1467 meta->GribVersion = version;
1468 if (version != 2) {
1469 printf ("Fast parsing doesn't handle this version because ReadGrib1Record/ReadTDLPRecord used Grib_Data[]\n");
1470 return -1;
1471 }
1472
1473 /*
1474 * Make room for entire message, and read it in.
1475 */
1476 /* nd5 needs to be gribLen in (sInt4) units rounded up. */
1477 nd5 = (gribLen + 3) / 4;
1478 if (nd5 > IS->ipackLen) {
1479 IS->ipackLen = nd5;
1480 IS->ipack = (sInt4 *) realloc ((void *) (IS->ipack),
1481 (IS->ipackLen) * sizeof (sInt4));
1482 }
1483 c_ipack = (unsigned char *) IS->ipack;
1484 /* Init last sInt4 to 0, to make sure that the padded bytes are 0. */
1485 IS->ipack[nd5 - 1] = 0;
1486 /* Init first 4 sInt4 to sect0. */
1487 memcpy (c_ipack, sect0, SECT0LEN_WORD * 4);
1488 /* Read in the rest of the message. */
1489 if (fread (c_ipack + SECT0LEN_WORD * 4, sizeof (char),
1490 (gribLen - SECT0LEN_WORD * 4),
1491 fp) != (gribLen - SECT0LEN_WORD * 4)) {
1492 errSprintf ("GribLen = %ld, SECT0Len_WORD = %d\n", gribLen,
1493 SECT0LEN_WORD);
1494 errSprintf ("Ran out of file\n");
1495 free (buff);
1496 return -1;
1497 }
1498
1499 /*
1500 * Make sure the arrays are large enough for call to unpacker library.
1501 */
1502 /* FindSectLen Does not want (ipack / c_ipack) word swapped, because
1503 * that would make it much more confusing to find bytes in c_ipack. */
1504 if (FindSectLen (c_ipack, gribLen, local_ns, &nd2x3, &table50) < 0) {
1505 preErrSprintf ("Inside ReadGrib2Record.. Calling FindSectLen\n");
1506 free (buff);
1507 return -2;
1508 }
1509
1510 /* Make sure all 'is' arrays except ns[7] are MAX (IS.ns[] ,
1511 * local_ns[]). See note 1 for reason to exclude ns[7] from MAX (). */
1512 for (i = 0; i < 7; i++) {
1513 if (local_ns[i] > IS->ns[i]) {
1514 IS->ns[i] = local_ns[i];
1515 IS->is[i] = (sInt4 *) realloc ((void *) (IS->is[i]),
1516 IS->ns[i] * sizeof (sInt4));
1517 }
1518 }
1519
1520 /* Allocate room for sect 2. If local_ns[2] = -1 there is no sect 2. */
1521 if (local_ns[2] == -1) {
1522 nidat = 10;
1523 nrdat = 10;
1524 } else {
1525 /*
1526 * See note 2) We have a section 2, so use:
1527 * MAX (32 * local_ns[2],SECT2_INTSIZE)
1528 * and MAX (32 * local_ns[2],SECT2_FLOATSIZE)
1529 * for size of section 2 unpacked.
1530 */
1531 nidat = (32 * local_ns[2] < SECT2_INIT_SIZE) ? SECT2_INIT_SIZE :
1532 32 * local_ns[2];
1533 nrdat = nidat;
1534 }
1535 if (nidat > IS->nidat) {
1536 IS->nidat = nidat;
1537 IS->idat = (sInt4 *) realloc ((void *) IS->idat,
1538 IS->nidat * sizeof (sInt4));
1539 }
1540 if (nrdat > IS->nrdat) {
1541 IS->nrdat = nrdat;
1542 IS->rdat = (float *) realloc ((void *) IS->rdat,
1543 IS->nrdat * sizeof (float));
1544 }
1545 /* Make sure we have room for the GRID part of the output. */
1546 if (nd2x3 > IS->nd2x3) {
1547 IS->nd2x3 = nd2x3;
1548 IS->iain = (sInt4 *) realloc ((void *) IS->iain,
1549 IS->nd2x3 * sizeof (sInt4));
1550 IS->ib = (sInt4 *) realloc ((void *) IS->ib,
1551 IS->nd2x3 * sizeof (sInt4));
1552 }
1553 /* See note 3) If table50 == 3, unpacker library needs nd5 >= nd2x3. */
1554 if ((table50 == 3) || (table50 == 0)) {
1555 if (nd5 < nd2x3) {
1556 nd5 = nd2x3;
1557 if (nd5 > IS->ipackLen) {
1558 IS->ipackLen = nd5;
1559 IS->ipack = (sInt4 *) realloc ((void *) (IS->ipack),
1560 IS->ipackLen * sizeof (sInt4));
1561 }
1562 /* Don't need to do the following, but we do in case code
1563 * changes. */
1564 c_ipack = (unsigned char *) IS->ipack;
1565 }
1566 }
1567 IS->nd5 = nd5;
1568 /* Unpacker library requires ipack to be MSB. */
1569 /*
1570 #ifdef DEBUG
1571 if (1==1) {
1572 FILE *fp = fopen ("test.bin", "wb");
1573 fwrite (IS->ipack, sizeof (sInt4), IS->nd5, fp);
1574 fclose (fp);
1575 }
1576 #endif
1577 */
1578 } else {
1579 c_ipack = (unsigned char *) IS->ipack;
1580 /* GRIB2 files are in big endian so c_ipack is as well. */
1581 #ifdef LITTLE_ENDIAN
1582 revmemcpy (&gribLen, &(c_ipack[12]), sizeof (sInt4));
1583 #else
1584 memcpy (&gribLen, &(c_ipack[12]), sizeof (sInt4));
1585 #endif
1586 }
1587 free (buff);
1588
1589 /* Loop through the grib message looking for the subgNum grid. subgNum
1590 * goes from 0 to n-1. */
1591 for (j = 0; j <= subgNum; j++) {
1592 if (j == 0) {
1593 inew = 1;
1594 } else {
1595 inew = 0;
1596 }
1597
1598 /* Note we are getting data back either as a float or an int, but not
1599 * both, so we don't need to allocated room for both. */
1600 /*
1601 unpk_grib2 (&kfildo, (float *) (IS->iain), IS->iain, &(IS->nd2x3),
1602 IS->idat, &(IS->nidat), IS->rdat, &(IS->nrdat), IS->is[0],
1603 &(IS->ns[0]), IS->is[1], &(IS->ns[1]), IS->is[2],
1604 &(IS->ns[2]), IS->is[3], &(IS->ns[3]), IS->is[4],
1605 &(IS->ns[4]), IS->is[5], &(IS->ns[5]), IS->is[6],
1606 &(IS->ns[6]), IS->is[7], &(IS->ns[7]), IS->ib, &ibitmap,
1607 IS->ipack, &(IS->nd5), &xmissp, &xmisss, &inew, &iclean,
1608 &l3264b, f_endMsg, jer, &ndjer, &kjer);
1609 */
1610 c_ipack = (unsigned char *)IS->ipack;
1611 unpk_g2ncep(&kfildo, (float *) (IS->iain), IS->iain, &(IS->nd2x3),
1612 IS->idat, &(IS->nidat), IS->rdat, &(IS->nrdat), IS->is[0],
1613 &(IS->ns[0]), IS->is[1], &(IS->ns[1]), IS->is[2],
1614 &(IS->ns[2]), IS->is[3], &(IS->ns[3]), IS->is[4],
1615 &(IS->ns[4]), IS->is[5], &(IS->ns[5]), IS->is[6],
1616 &(IS->ns[6]), IS->is[7], &(IS->ns[7]), IS->ib, &ibitmap,
1617 c_ipack, &(IS->nd5), &xmissp, &xmisss, &inew, &iclean,
1618 &l3264b, f_endMsg, jer, &ndjer, &kjer);
1619
1620
1621 /*
1622 * Check for error messages...
1623 * If we get an error message, print it, and return.
1624 */
1625 for (i = 0; i < (uInt4) kjer; i++) {
1626 if (jer[ndjer + i] == 0) {
1627 /* no error. */
1628 } else if (jer[ndjer + i] == 1) {
1629 /* Warning. */
1630 #ifdef DEBUG
1631 printf ("Warning: Unpack library warning code (%ld %ld)\n",
1632 jer[i], jer[ndjer + i]);
1633 #endif
1634 } else {
1635 /* BAD Error. */
1636 errSprintf ("ERROR: Unpack library error code (%ld %ld)\n",
1637 jer[i], jer[ndjer + i]);
1638 return -3;
1639 }
1640 }
1641 }
1642
1643 /* Parse the meta data out. */
1644 if (MetaParse (meta, IS->is[0], IS->ns[0], IS->is[1], IS->ns[1],
1645 IS->is[2], IS->ns[2], IS->rdat, IS->nrdat, IS->idat,
1646 IS->nidat, IS->is[3], IS->ns[3], IS->is[4], IS->ns[4],
1647 IS->is[5], IS->ns[5], gribLen, xmissp, xmisss, simpVer, simpWWA)
1648 != 0) {
1649 #ifdef DEBUG
1650 FILE *fp;
1651 if ((fp = fopen ("dump.is0", "wt")) != NULL) {
1652 for (i = 0; i < 8; i++) {
1653 fprintf (fp, "---Section %d---\n", i);
1654 for (j = 1; j <= IS->ns[i]; j++) {
1655 fprintf (fp, "IS%d Item %d = %ld\n", i, j, IS->is[i][j - 1]);
1656 }
1657 }
1658 fclose (fp);
1659 }
1660 #endif
1661 preErrSprintf ("Inside ReadGrib2Record.. Problems in MetaParse\n");
1662 return -4;
1663 }
1664
1665 if ((majEarth > 6000) && (majEarth < 7000)) {
1666 if ((minEarth > 6000) && (minEarth < 7000)) {
1667 meta->gds.f_sphere = 0;
1668 meta->gds.majEarth = majEarth;
1669 meta->gds.minEarth = minEarth;
1670 } else {
1671 meta->gds.f_sphere = 1;
1672 meta->gds.majEarth = majEarth;
1673 meta->gds.minEarth = majEarth;
1674 }
1675 }
1676
1677 /* Figure out an equation to pass to ParseGrid to convert the units for
1678 * this grid. */
1679 /*
1680 if (ComputeUnit (meta->pds2.prodType, meta->pds2.sect4.templat,
1681 meta->pds2.sect4.cat, meta->pds2.sect4.subcat, f_unit,
1682 &unitM, &unitB, unitName) == 0) {
1683 */
1684 if (ComputeUnit (meta->convert, meta->unitName, f_unit, &unitM, &unitB,
1685 unitName) == 0) {
1686 unitLen = strlen (unitName);
1687 meta->unitName = (char *) realloc ((void *) (meta->unitName),
1688 1 + unitLen * sizeof (char));
1689 strncpy (meta->unitName, unitName, unitLen);
1690 meta->unitName[unitLen] = '\0';
1691 }
1692
1693 /* compute the subgrid. */
1694 if ((lwlf->lat != -100) && (uprt->lat != -100)) {
1695 Nx = meta->gds.Nx;
1696 Ny = meta->gds.Ny;
1697 if (computeSubGrid (lwlf, &x1, &y1, uprt, &x2, &y2, &(meta->gds),
1698 &newGds) != 0) {
1699 preErrSprintf ("ERROR: In compute subgrid.\n");
1700 return 1;
1701 }
1702 /* I couldn't decide if I should "permanently" change the GDS or not.
1703 * when I wrote computeSubGrid. If next line stays, really should
1704 * rewrite computeSubGrid. */
1705 memcpy (&(meta->gds), &newGds, sizeof (gdsType));
1706 f_subGrid = 1;
1707 } else {
1708 Nx = meta->gds.Nx;
1709 Ny = meta->gds.Ny;
1710 x1 = 1;
1711 x2 = Nx;
1712 y1 = 1;
1713 y2 = Ny;
1714 f_subGrid = 0;
1715 }
1716
1717 if ((f_subGrid) && (meta->gds.scan != 64)) {
1718 errSprintf ("Can not do a subgrid of non scanmode 64 grid yet.\n");
1719 return -3;
1720 }
1721
1722 /* Figure out some other non-section oriented meta data. */
1723 /* strftime (meta->refTime, 20, "%Y%m%d%H%M",
1724 gmtime (&(meta->pds2.refTime)));
1725 */
1726 Clock_Print (meta->refTime, 20, meta->pds2.refTime, "%Y%m%d%H%M", 0);
1727 /*
1728 strftime (meta->validTime, 20, "%Y%m%d%H%M",
1729 gmtime (&(meta->pds2.sect4.validTime)));
1730 */
1731 Clock_Print (meta->validTime, 20, meta->pds2.sect4.validTime,
1732 "%Y%m%d%H%M", 0);
1733
1734 meta->deltTime = (sInt4) (meta->pds2.sect4.validTime - meta->pds2.refTime);
1735
1736 return 0;
1737 }
1738 #endif
1739