Line data Source code
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 : */
119 1244 : 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 1244 : 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 1244 : curLen = 8;
142 1244 : if (*buffLen < curLen) {
143 1233 : *buffLen = curLen;
144 1233 : *buff = (char *) realloc ((void *) *buff, *buffLen * sizeof (char));
145 : }
146 1244 : if (VSIFReadL(*buff, sizeof (char), curLen, fp) != curLen) {
147 1 : errSprintf ("ERROR: Couldn't find 'GRIB' or 'TDLP'\n");
148 1 : 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 3325 : while (
159 : #if 0
160 : /* tdlpack is no longer supported by GDAL */
161 : (tdlpMatch != 4) &&
162 : #endif
163 4568 : (gribMatch != 4)) {
164 5407 : for (i = curLen - 8; i + 7 < curLen; i++) {
165 3325 : if ((*buff)[i] == 'G') {
166 1245 : if (((*buff)[i + 1] == 'R') && ((*buff)[i + 2] == 'I') &&
167 1243 : ((*buff)[i + 3] == 'B')) {
168 1243 : if (((*buff)[i + 7] == 1) ||
169 1102 : ((*buff)[i + 7] == 2)) {
170 1243 : gribMatch = 4;
171 1243 : 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 3325 : stillNeed = i - (curLen - 8);
187 : /* Read enough of message to have the first 8 bytes (including ID). */
188 3325 : if (stillNeed != 0) {
189 2082 : curLen += stillNeed;
190 2082 : if ((limit >= 0) && (curLen > (size_t) limit)) {
191 0 : errSprintf ("ERROR: Couldn't find type in %ld bytes\n", limit);
192 0 : *buffLen = curLen - stillNeed;
193 0 : return -1;
194 : }
195 2082 : if (*buffLen < curLen) {
196 22 : myAssert (200 > stillNeed);
197 22 : *buffLen = *buffLen + 200;
198 : /* *buffLen = curLen; */
199 22 : *buff = (char *) realloc ((void *) *buff,
200 22 : *buffLen * sizeof (char));
201 : }
202 2082 : if (VSIFReadL((*buff) + (curLen - stillNeed), sizeof (char), stillNeed, fp) != stillNeed) {
203 0 : errSprintf ("ERROR: Ran out of file reading SECT0\n");
204 0 : *buffLen = curLen;
205 0 : 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 1243 : *buffLen = curLen;
213 :
214 : /* curLen and (*buff) hold 8 bytes of section 0. */
215 1243 : curLen -= 8;
216 1243 : memcpy (&(sect0[0]), (*buff) + curLen, 4);
217 : #ifdef DEBUG
218 : #ifdef LITTLE_ENDIAN
219 1243 : myAssert ((sect0[0] == 1112101447L) || (sect0[0] == 1347175508L));
220 : #else
221 : myAssert ((sect0[0] == 1196575042L) || (sect0[0] == 1413762128L));
222 : #endif
223 : #endif
224 1243 : memcpy (&(sect0[1]), *buff + curLen + 4, 4);
225 : /* Make sure we don't pass back part of "GRIB" in the buffer. */
226 1243 : (*buff)[curLen] = '\0';
227 1243 : *buffLen = curLen;
228 :
229 1243 : 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 1243 : if (word.buffer[3] == 1) {
249 141 : *version = 1;
250 : /* Find out the GRIB Message Length */
251 141 : *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 141 : if (*gribLen < 52) {
255 0 : errSprintf ("GRIB1 length %ld was < 52?\n", *gribLen);
256 0 : return -5;
257 : }
258 1102 : } else if (word.buffer[3] == 2) {
259 1102 : *version = 2;
260 : /* Make sure we still have enough file for the rest of section 0. */
261 1102 : if (VSIFReadL(sect0 + 2, sizeof (sInt4), 2, fp) != 2) {
262 0 : errSprintf ("ERROR: Ran out of file reading SECT0\n");
263 0 : return -2;
264 : }
265 1102 : if (sect0[2] != 0) {
266 0 : errSprintf ("Most significant sInt4 of GRIB length was not 0?\n");
267 0 : 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 0 : return -4;
271 : }
272 : #ifdef LITTLE_ENDIAN
273 1102 : revmemcpy (gribLen, &(sect0[3]), sizeof (sInt4));
274 : #else
275 : *gribLen = sect0[3];
276 : #endif
277 : } else {
278 0 : errSprintf ("ERROR: Not TDLPack, and Grib edition is not 1 or 2\n");
279 0 : return -3;
280 : }
281 1243 : 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 : */
319 0 : int FindGRIBMsg (VSILFILE *fp, int msgNum, sInt4 *offset, int *curMsg)
320 : {
321 : int cnt; /* The current message we are looking at. */
322 0 : 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 0 : cnt = *curMsg + 1;
331 0 : buffLen = 0;
332 0 : while (VSIFReadL(&c, sizeof(char), 1, fp) == 1) {
333 0 : VSIFSeekL(fp, VSIFTellL(fp) - sizeof(char), SEEK_SET);
334 :
335 0 : if (cnt >= msgNum) {
336 : /* 12/1/2004 version 1.63 forgot to free buff */
337 0 : free (buff);
338 0 : *curMsg = cnt;
339 0 : return 0;
340 : }
341 : /* Read section 0 to find gribLen and wmoLen. */
342 0 : if (ReadSECT0 (fp, &buff, &buffLen, GRIB_LIMIT, sect0, &gribLen,
343 0 : &version) < 0) {
344 0 : preErrSprintf ("Inside FindGRIBMsg\n");
345 0 : free (buff);
346 0 : return -1;
347 : }
348 0 : myAssert ((version == 1) || (version == 2) || (version == -1));
349 : /* Continue on to the next grib message. */
350 0 : if ((version == 1) || (version == -1)) {
351 0 : jump = gribLen - 8;
352 : } else {
353 0 : jump = gribLen - 16;
354 : }
355 0 : VSIFSeekL(fp, jump, SEEK_CUR);
356 0 : *offset = *offset + gribLen + buffLen;
357 0 : cnt++;
358 : }
359 0 : free (buff);
360 0 : *curMsg = cnt - 1;
361 : /* Return -2 since we reached the end of file. This may not be an error
362 : * (multiple file option). */
363 0 : 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 : */
403 466 : 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 466 : 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 457 : if (*curTot > gribLen - (6 + 4)) {
414 0 : errSprintf ("ERROR: Ran out of data in Section 2 or 3\n");
415 0 : return -1;
416 : }
417 : /* Handle optional section 2. */
418 457 : if (c_ipack[*curTot + 4] == 2) {
419 421 : MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
420 421 : if( sectLen < 0 || *curTot > INT_MAX - sectLen ) {
421 0 : errSprintf ("ERROR: Invalid sectLen for section 2\n");
422 0 : return -1;
423 : }
424 421 : *curTot = *curTot + sectLen;
425 421 : if (ns[2] < sectLen)
426 421 : ns[2] = sectLen;
427 : /* ERO: check change from + 5 to +6+4 per r39022 */
428 421 : if (*curTot + 6 + 4 > gribLen) {
429 0 : errSprintf ("ERROR: Ran out of data in Section 3\n");
430 0 : return -1;
431 : }
432 : }
433 : /* Handle section 3. */
434 457 : if (c_ipack[*curTot + 4] != 3) {
435 0 : errSprintf ("ERROR: Section 3 labeled as %d\n",
436 0 : c_ipack[*curTot + 4]);
437 0 : return -2;
438 : }
439 457 : MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
440 457 : if( sectLen < 0 || *curTot > INT_MAX - sectLen ) {
441 0 : errSprintf ("ERROR: Invalid sectLen for section 3\n");
442 0 : return -1;
443 : }
444 457 : if (ns[3] < sectLen)
445 457 : ns[3] = sectLen;
446 : /* While we are here, grab the total number of grid points nd2x3. */
447 457 : MEMCPY_BIG (&li_temp, c_ipack + *curTot + 6, 4);
448 457 : if (*nd2x3 < li_temp)
449 457 : *nd2x3 = li_temp;
450 457 : *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 466 : if (*curTot > gribLen - 5) {
460 0 : errSprintf ("ERROR: Ran out of data in Section 4\n");
461 0 : return -1;
462 : }
463 466 : if (c_ipack[*curTot + 4] != 4) {
464 0 : errSprintf ("ERROR: Section 4 labeled as %d\n", c_ipack[*curTot + 4]);
465 0 : return -2;
466 : }
467 466 : MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
468 466 : if( sectLen < 0 || *curTot > INT_MAX - sectLen ) {
469 0 : errSprintf ("ERROR: Invalid sectLen for section 4\n");
470 0 : return -1;
471 : }
472 466 : if (ns[4] < sectLen)
473 457 : ns[4] = sectLen;
474 466 : *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 466 : if (*curTot > gribLen - (9 + 2)) {
484 0 : errSprintf ("ERROR: Ran out of data in Section 5\n");
485 0 : return -1;
486 : }
487 466 : if (c_ipack[*curTot + 4] != 5) {
488 0 : errSprintf ("ERROR: Section 5 labeled as %d\n", c_ipack[*curTot + 4]);
489 0 : return -2;
490 : }
491 466 : MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
492 466 : if( sectLen < 0 || *curTot > INT_MAX - sectLen ) {
493 0 : errSprintf ("ERROR: Invalid sectLen for section 5\n");
494 0 : return -1;
495 : }
496 : /* While we are here, grab the packing method. */
497 466 : MEMCPY_BIG (table50, c_ipack + *curTot + 9, 2);
498 466 : if (ns[5] < sectLen)
499 457 : ns[5] = sectLen;
500 466 : *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 466 : if (*curTot > gribLen - 5) {
509 0 : errSprintf ("ERROR: Ran out of data in Section 6\n");
510 0 : return -1;
511 : }
512 466 : if (c_ipack[*curTot + 4] != 6) {
513 0 : errSprintf ("ERROR: Section 6 labeled as %d\n", c_ipack[*curTot + 4]);
514 0 : return -2;
515 : }
516 466 : MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
517 466 : if( sectLen < 0 || *curTot > INT_MAX - sectLen ) {
518 0 : errSprintf ("ERROR: Invalid sectLen for section 6\n");
519 0 : return -1;
520 : }
521 466 : if (ns[6] < sectLen)
522 457 : ns[6] = sectLen;
523 466 : *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 466 : if (*curTot > gribLen - 5) {
532 0 : errSprintf ("ERROR: Ran out of data in Section 7\n");
533 0 : return -1;
534 : }
535 466 : if (c_ipack[*curTot + 4] != 7) {
536 0 : errSprintf ("ERROR: Section 7 labeled as %d\n", c_ipack[*curTot + 4]);
537 0 : return -2;
538 : }
539 466 : MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
540 466 : if( sectLen < 0 || *curTot > INT_MAX - sectLen ) {
541 0 : errSprintf ("ERROR: Invalid sectLen for section 7\n");
542 0 : return -1;
543 : }
544 466 : if (ns[7] < sectLen)
545 457 : ns[7] = sectLen;
546 466 : *curTot = *curTot + sectLen;
547 : /*
548 : #ifdef DEBUG
549 : printf ("Section len (7=%ld < %ld)\n", sectLen, ns[7]);
550 : #endif
551 : */
552 466 : 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 : */
590 457 : 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 457 : ns[0] = SECT0LEN_WORD * 4;
600 457 : curTot = ns[0];
601 :
602 : /* Figure out the size of section 1. */
603 457 : if (curTot + 5 > gribLen) {
604 0 : errSprintf ("ERROR: Ran out of data in Section 1\n");
605 0 : return -1;
606 : }
607 457 : if (c_ipack[curTot + 4] != 1) {
608 0 : errSprintf ("ERROR: Section 1 labeled as %d\n", c_ipack[curTot + 4]);
609 0 : return -2;
610 : }
611 457 : MEMCPY_BIG (&(ns[1]), c_ipack + curTot, 4);
612 457 : curTot += ns[1];
613 : /*
614 : #ifdef DEBUG
615 : printf ("Section len (0=%ld) (1=%ld)\n", ns[0], ns[1]);
616 : #endif
617 : */
618 457 : sectNum = 2;
619 3199 : for (i = 2; i < 8; i++) {
620 2742 : ns[i] = -1;
621 : }
622 457 : *nd2x3 = -1;
623 9 : do {
624 466 : if ((ans = FindSectLen2to7 (c_ipack, gribLen, ns, sectNum, &curTot,
625 466 : nd2x3, table50)) != 0) {
626 0 : return ans;
627 : }
628 466 : if( curTot + 4 > gribLen ) {
629 0 : errSprintf ("ERROR: Ran out of data in Section 1\n");
630 0 : 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 466 : memcpy (§Len, c_ipack + curTot, 4);
636 466 : if (sectLen == 926365495L) {
637 457 : sectNum = 8;
638 : } else {
639 9 : if (curTot + 4 == gribLen) {
640 0 : errSprintf("ERROR: Ran out of data in Section 1\n");
641 0 : return -1;
642 : }
643 9 : sectNum = c_ipack[curTot + 4];
644 9 : if ((sectNum < 2) || (sectNum > 7)) {
645 0 : errSprintf ("ERROR (FindSectLen): Couldn't find the end of the "
646 : "message\n");
647 0 : errSprintf ("and it doesn't appear to repeat sections.\n");
648 0 : errSprintf ("so it is probably an ASCII / binary bug\n");
649 0 : errSprintf ("Max Sect Lengths: %ld %ld %ld %ld %ld %ld %ld"
650 0 : " %ld\n", ns[0], ns[1], ns[2], ns[3], ns[4], ns[5],
651 0 : ns[6], ns[7]);
652 0 : return -2;
653 : }
654 : }
655 466 : } while (sectNum != 8);
656 457 : 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 : */
693 523 : void IS_Init (IS_dataType *is)
694 : {
695 : int i; /* A simple loop counter. */
696 523 : is->ns[0] = 16;
697 523 : is->ns[1] = 21;
698 523 : is->ns[2] = 7;
699 523 : is->ns[3] = 96;
700 523 : is->ns[4] = 130; /* 60->130 in case there are some S4 time intervals */
701 523 : is->ns[5] = 49;
702 523 : is->ns[6] = 6;
703 523 : is->ns[7] = 8;
704 4707 : for (i = 0; i < 8; i++) {
705 4184 : is->is[i] = (sInt4 *) calloc (is->ns[i], sizeof (sInt4));
706 : }
707 : /* Allocate grid memory. */
708 523 : is->nd2x3 = 0;
709 523 : is->iain = nullptr;
710 523 : is->ib = nullptr;
711 : /* Allocate section 2 int memory. */
712 523 : is->nidat = 0;
713 523 : is->idat = nullptr;
714 : /* Allocate section 2 float memory. */
715 523 : is->nrdat = 0;
716 523 : is->rdat = nullptr;
717 : /* Allocate storage for ipack. */
718 523 : is->ipackLen = 0;
719 523 : is->ipack = nullptr;
720 523 : }
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 : */
747 523 : void IS_Free (IS_dataType *is)
748 : {
749 : int i; /* A simple loop counter. */
750 4707 : for (i = 0; i < 8; i++) {
751 4184 : free (is->is[i]);
752 4184 : is->is[i] = nullptr;
753 4184 : is->ns[i] = 0;
754 : }
755 : /* Free grid memory. */
756 523 : free (is->iain);
757 523 : is->iain = nullptr;
758 523 : free (is->ib);
759 523 : is->ib = nullptr;
760 523 : is->nd2x3 = 0;
761 : /* Free section 2 int memory. */
762 523 : free (is->idat);
763 523 : is->idat = nullptr;
764 523 : is->nidat = 0;
765 : /* Free section 2 float memory. */
766 523 : free (is->rdat);
767 523 : is->rdat = nullptr;
768 523 : is->nrdat = 0;
769 : /* Free storage for ipack. */
770 523 : free (is->ipack);
771 523 : is->ipack = nullptr;
772 523 : is->ipackLen = 0;
773 523 : }
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
882 523 : 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 523 : 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 523 : 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 523 : sInt4 iclean = 0; /* 0 embed the missing values, 1 don't. */
911 : int j; /* Counter used to find the desired subgrid. */
912 523 : sInt4 kfildo = 5; /* FORTRAN Unit number for diagnostic info. Ignored,
913 : * unless library is compiled a particular way. */
914 523 : sInt4 ibitmap = 0; /* 0 means no bitmap returned, otherwise 1. */
915 523 : 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 523 : 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 523 : 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 523 : l3264b = sizeof (sInt4) * 8;
948 523 : buffLen = 0;
949 523 : if (*f_endMsg == 1) {
950 523 : if (ReadSECT0 (fp, &buff, &buffLen, -1, sect0, &gribLen, &version) < 0) {
951 0 : preErrSprintf ("Inside ReadGrib2Record\n");
952 0 : free (buff);
953 0 : return -1;
954 : }
955 523 : meta->GribVersion = version;
956 523 : if (version == 1) {
957 66 : if (ReadGrib1Record (fp, f_unit, Grib_Data, grib_DataLen, meta, IS,
958 66 : sect0, gribLen, majEarth, minEarth) != 0) {
959 0 : preErrSprintf ("Problems with ReadGrib1Record called by "
960 : "ReadGrib2Record\n");
961 0 : free (buff);
962 0 : return -1;
963 : }
964 66 : *f_endMsg = 1;
965 66 : free (buff);
966 66 : 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 457 : if( gribLen > UINT_MAX - 3 )
988 : {
989 0 : errSprintf("Invalid value of gribLen");
990 0 : free (buff);
991 0 : return -1;
992 : }
993 457 : nd5 = (gribLen + 3) / 4;
994 457 : if (nd5 > IS->ipackLen) {
995 457 : if( gribLen > 100 * 1024 * 1024 )
996 : {
997 0 : vsi_l_offset curPos = VSIFTellL(fp);
998 0 : VSIFSeekL(fp, 0, SEEK_END);
999 0 : vsi_l_offset fileSize = VSIFTellL(fp);
1000 0 : VSIFSeekL(fp, curPos, SEEK_SET);
1001 0 : if( fileSize < gribLen )
1002 : {
1003 0 : errSprintf("File too short");
1004 0 : free (buff);
1005 0 : return -1;
1006 : }
1007 : }
1008 457 : sInt4* ipackNew = (sInt4 *) realloc ((void *) (IS->ipack),
1009 457 : nd5 * sizeof (sInt4));
1010 457 : if( ipackNew == nullptr )
1011 : {
1012 0 : errSprintf("Out of memory");
1013 0 : free (buff);
1014 0 : return -1;
1015 : }
1016 457 : IS->ipackLen = nd5;
1017 457 : IS->ipack = ipackNew;
1018 : }
1019 457 : c_ipack = (unsigned char *) IS->ipack;
1020 : /* Init last sInt4 to 0, to make sure that the padded bytes are 0. */
1021 457 : IS->ipack[nd5 - 1] = 0;
1022 : /* Init first 4 sInt4 to sect0. */
1023 457 : memcpy (c_ipack, sect0, SECT0LEN_WORD * 4);
1024 : /* Read in the rest of the message. */
1025 914 : if (VSIFReadL (c_ipack + SECT0LEN_WORD * 4, sizeof (char),
1026 457 : (gribLen - SECT0LEN_WORD * 4), fp) != (gribLen - SECT0LEN_WORD * 4)) {
1027 0 : errSprintf ("GribLen = %ld, SECT0Len_WORD = %d\n", gribLen,
1028 : SECT0LEN_WORD);
1029 0 : errSprintf ("Ran out of file\n");
1030 0 : free (buff);
1031 0 : 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 457 : if (FindSectLen (c_ipack, gribLen, local_ns, &nd2x3, &table50) < 0) {
1040 0 : preErrSprintf ("Inside ReadGrib2Record.. Calling FindSectLen\n");
1041 0 : free (buff);
1042 0 : return -2;
1043 : }
1044 457 : if( nd2x3 < 0 || nd2x3 > INT_MAX / static_cast<int>(sizeof (sInt4)) )
1045 : {
1046 0 : preErrSprintf ("Invalid nd2x3\n");
1047 0 : free (buff);
1048 0 : 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 3656 : for (i = 0; i < 7; i++) {
1054 3199 : if (local_ns[i] > IS->ns[i]) {
1055 10 : IS->ns[i] = local_ns[i];
1056 10 : IS->is[i] = (sInt4 *) realloc ((void *) (IS->is[i]),
1057 10 : 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 457 : if (local_ns[2] == -1) {
1063 36 : nidat = 10;
1064 36 : 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 421 : nidat = (32 * local_ns[2] < SECT2_INIT_SIZE) ? SECT2_INIT_SIZE :
1073 : 32 * local_ns[2];
1074 421 : nrdat = nidat;
1075 : }
1076 457 : if (nidat > IS->nidat) {
1077 457 : IS->idat = (sInt4 *) realloc ((void *) IS->idat,
1078 457 : nidat * sizeof (sInt4));
1079 : // Initialization needed to avoid use of uninitialized memory in
1080 : // ParseSect2_Unknown() on dataset of https://github.com/OSGeo/gdal/issues/5290
1081 457 : if( IS->nidat == 0 )
1082 457 : IS->idat[0] = 0;
1083 457 : IS->nidat = nidat;
1084 : }
1085 457 : if (nrdat > IS->nrdat) {
1086 457 : IS->rdat = (float *) realloc ((void *) IS->rdat,
1087 457 : nrdat * sizeof (float));
1088 : // Initialization needed to avoid use of uninitialized memory in
1089 : // ParseSect2_Unknown() on dataset of https://github.com/OSGeo/gdal/issues/5290
1090 457 : if( IS->nrdat == 0 )
1091 457 : IS->rdat[0] = 0;
1092 457 : IS->nrdat = nrdat;
1093 : }
1094 : /* Make sure we have room for the GRID part of the output. */
1095 457 : if (nd2x3 > IS->nd2x3) {
1096 457 : if( nd2x3 > 100 * 1024 * 1024 )
1097 : {
1098 :
1099 0 : vsi_l_offset curPos = VSIFTellL(fp);
1100 0 : VSIFSeekL(fp, 0, SEEK_END);
1101 0 : vsi_l_offset fileSize = VSIFTellL(fp);
1102 0 : VSIFSeekL(fp, curPos, SEEK_SET);
1103 : // allow a compression ratio of 1:1000
1104 0 : if( (vsi_l_offset)(nd2x3 / 1000) > fileSize )
1105 : {
1106 0 : preErrSprintf ("File too short\n");
1107 0 : free (buff);
1108 0 : return -2;
1109 : }
1110 : }
1111 :
1112 457 : IS->nd2x3 = nd2x3;
1113 457 : if( Grib_Data )
1114 : {
1115 157 : IS->iain = (sInt4 *) realloc ((void *) IS->iain,
1116 157 : IS->nd2x3 * sizeof (sInt4));
1117 : }
1118 457 : IS->ib = (sInt4 *) realloc ((void *) IS->ib,
1119 457 : IS->nd2x3 * sizeof (sInt4));
1120 : }
1121 : /* See note 3) If table50 == 3, unpacker library needs nd5 >= nd2x3. */
1122 457 : if ((table50 == 3) || (table50 == 0)) {
1123 179 : if (nd5 < nd2x3) {
1124 158 : nd5 = nd2x3;
1125 158 : if (nd5 > IS->ipackLen) {
1126 158 : sInt4* ipackNew = (sInt4 *) realloc ((void *) (IS->ipack),
1127 158 : nd5 * sizeof (sInt4));
1128 158 : if( ipackNew == nullptr )
1129 : {
1130 0 : errSprintf("Out of memory");
1131 0 : free (buff);
1132 0 : return -1;
1133 : }
1134 : // zero initialize to avoid later use of undefined values
1135 : // Not sure if it is strictly needed (perhaps some later
1136 : // steps should detect errors / data shortening), but
1137 : // cannot hurt.
1138 : // see https://trac.osgeo.org/gdal/ticket/6967
1139 158 : memset(ipackNew + IS->ipackLen, 0,
1140 158 : (nd5 - IS->ipackLen) * sizeof(sInt4));
1141 158 : IS->ipackLen = nd5;
1142 158 : IS->ipack = ipackNew;
1143 : }
1144 : /* Don't need to do the following, but we do in case code
1145 : * changes. */
1146 158 : c_ipack = (unsigned char *) IS->ipack;
1147 : }
1148 : }
1149 457 : IS->nd5 = nd5;
1150 : /* Unpacker library requires ipack to be MSB. */
1151 : /*
1152 : #ifdef DEBUG
1153 : if (1==1) {
1154 : FILE *fp = fopen ("test.bin", "wb");
1155 : fwrite (IS->ipack, sizeof (sInt4), IS->nd5, fp);
1156 : fclose (fp);
1157 : }
1158 : #endif
1159 : */
1160 : /*
1161 : #ifdef LITTLE_ENDIAN
1162 : memswp (IS->ipack, sizeof (sInt4), IS->nd5);
1163 : #endif
1164 : */
1165 : } else {
1166 0 : c_ipack = (unsigned char *) IS->ipack;
1167 : /* GRIB2 files are in big endian so c_ipack is as well. */
1168 : #ifdef LITTLE_ENDIAN
1169 0 : revmemcpy (&gribLen, &(c_ipack[12]), sizeof (sInt4));
1170 : #else
1171 : memcpy (&gribLen, &(c_ipack[12]), sizeof (sInt4));
1172 : #endif
1173 : }
1174 457 : free (buff);
1175 :
1176 : /* Loop through the grib message looking for the subgNum grid. subgNum
1177 : * goes from 0 to n-1. */
1178 917 : for (j = 0; j <= subgNum; j++) {
1179 460 : if (j == 0) {
1180 457 : inew = 1;
1181 : } else {
1182 3 : inew = 0;
1183 : }
1184 :
1185 : /* Note we are getting data back either as a float or an int, but not
1186 : * both, so we don't need to allocated room for both. */
1187 460 : unpk_g2ncep (&kfildo,
1188 : j == subgNum ? (float *) (IS->iain) : nullptr,
1189 : j == subgNum ? IS->iain : nullptr,
1190 : &(IS->nd2x3),
1191 : IS->idat, &(IS->nidat), IS->rdat, &(IS->nrdat), IS->is[0],
1192 : &(IS->ns[0]), IS->is[1], &(IS->ns[1]), IS->is[2],
1193 : &(IS->ns[2]), IS->is[3], &(IS->ns[3]), IS->is[4],
1194 : &(IS->ns[4]), IS->is[5], &(IS->ns[5]), IS->is[6],
1195 : &(IS->ns[6]), IS->is[7], &(IS->ns[7]), IS->ib, &ibitmap,
1196 : c_ipack, &(IS->nd5), &xmissp, &xmisss, &inew, &iclean,
1197 : &l3264b, f_endMsg, jer, &ndjer, &kjer);
1198 : /*
1199 : unpk_grib2 (&kfildo, (float *) (IS->iain), IS->iain, &(IS->nd2x3),
1200 : IS->idat, &(IS->nidat), IS->rdat, &(IS->nrdat), IS->is[0],
1201 : &(IS->ns[0]), IS->is[1], &(IS->ns[1]), IS->is[2],
1202 : &(IS->ns[2]), IS->is[3], &(IS->ns[3]), IS->is[4],
1203 : &(IS->ns[4]), IS->is[5], &(IS->ns[5]), IS->is[6],
1204 : &(IS->ns[6]), IS->is[7], &(IS->ns[7]), IS->ib, &ibitmap,
1205 : IS->ipack, &(IS->nd5), &xmissp, &xmisss, &inew, &iclean,
1206 : &l3264b, f_endMsg, jer, &ndjer, &kjer);
1207 : */
1208 : /*
1209 : * Check for error messages...
1210 : * If we get an error message, print it, and return.
1211 : */
1212 4140 : for (i = 0; i < (uInt4) kjer; i++) {
1213 3680 : if (jer[ndjer + i] == 0) {
1214 : /* no error. */
1215 0 : } else if (jer[ndjer + i] == 1) {
1216 : /* Warning. */
1217 : #ifdef DEBUG
1218 0 : printf ("Warning: Unpack library warning code (%d %d)\n",
1219 0 : jer[i], jer[ndjer + i]);
1220 : #endif
1221 : } else {
1222 : /* BAD Error. */
1223 0 : errSprintf ("ERROR: Unpack library error code (%ld %ld)\n",
1224 0 : jer[i], jer[ndjer + i]);
1225 0 : return -3;
1226 : }
1227 : }
1228 : }
1229 :
1230 : /* Parse the meta data out. */
1231 457 : if (MetaParse (meta, IS->is[0], IS->ns[0], IS->is[1], IS->ns[1],
1232 : IS->is[2], IS->ns[2], IS->rdat, IS->nrdat, IS->idat,
1233 : IS->nidat, IS->is[3], IS->ns[3], IS->is[4], IS->ns[4],
1234 : IS->is[5], IS->ns[5], gribLen, xmissp, xmisss, simpVer, simpWWA)
1235 457 : != 0) {
1236 : #ifdef DEBUG
1237 : FILE *l_fp;
1238 0 : if ((l_fp = fopen ("dump.is0", "wt")) != nullptr) {
1239 0 : for (i = 0; i < 8; i++) {
1240 0 : fprintf (l_fp, "---Section %d---\n", (int) i);
1241 0 : for (j = 1; j <= IS->ns[i]; j++) {
1242 0 : fprintf (l_fp, "IS%d Item %d = %d\n", (int) i, (int) j, IS->is[i][j - 1]);
1243 : }
1244 : }
1245 0 : fclose (l_fp);
1246 : }
1247 : #endif
1248 0 : preErrSprintf ("Inside ReadGrib2Record.. Problems in MetaParse\n");
1249 0 : return -4;
1250 : }
1251 :
1252 457 : if ((majEarth > 6000) && (majEarth < 7000)) {
1253 0 : if ((minEarth > 6000) && (minEarth < 7000)) {
1254 0 : meta->gds.f_sphere = 0;
1255 0 : meta->gds.majEarth = majEarth;
1256 0 : meta->gds.minEarth = minEarth;
1257 : } else {
1258 0 : meta->gds.f_sphere = 1;
1259 0 : meta->gds.majEarth = majEarth;
1260 0 : meta->gds.minEarth = majEarth;
1261 : }
1262 : }
1263 :
1264 : /* Figure out an equation to pass to ParseGrid to convert the units for
1265 : * this grid. */
1266 : /*
1267 : if (ComputeUnit (meta->pds2.prodType, meta->pds2.sect4.templat,
1268 : meta->pds2.sect4.cat, meta->pds2.sect4.subcat, f_unit,
1269 : &unitM, &unitB, unitName) == 0) {
1270 : */
1271 457 : if (ComputeUnit (meta->convert, meta->unitName, f_unit, &unitM, &unitB,
1272 457 : unitName) == 0) {
1273 61 : unitLen = static_cast<int>(strlen (unitName));
1274 61 : meta->unitName = (char *) realloc ((void *) (meta->unitName),
1275 61 : 1 + unitLen * sizeof (char));
1276 61 : memcpy (meta->unitName, unitName, unitLen);
1277 61 : meta->unitName[unitLen] = '\0';
1278 : }
1279 :
1280 : /* compute the subgrid. */
1281 : /*
1282 : if ((lwlf->lat != -100) && (uprt->lat != -100)) {
1283 : Nx = meta->gds.Nx;
1284 : Ny = meta->gds.Ny;
1285 : if (computeSubGrid (lwlf, &x1, &y1, uprt, &x2, &y2, &(meta->gds),
1286 : &newGds) != 0) {
1287 : preErrSprintf ("ERROR: In compute subgrid.\n");
1288 : return 1;
1289 : }
1290 : // I couldn't decide if I should "permanently" change the GDS or not.
1291 : // when I wrote computeSubGrid. If next line stays, really should
1292 : // rewrite computeSubGrid.
1293 : memcpy (&(meta->gds), &newGds, sizeof (gdsType));
1294 : f_subGrid = 1;
1295 : } else {
1296 : Nx = meta->gds.Nx;
1297 : Ny = meta->gds.Ny;
1298 : x1 = 1;
1299 : x2 = Nx;
1300 : y1 = 1;
1301 : y2 = Ny;
1302 : f_subGrid = 0;
1303 : }
1304 : */
1305 457 : if( meta->gds.Nx == 0 || meta->gds.Nx > INT_MAX ||
1306 457 : meta->gds.Ny == 0 || meta->gds.Ny > INT_MAX )
1307 : {
1308 0 : errSprintf ("Invalid grid dimension.\n");
1309 0 : return -3;
1310 : }
1311 457 : Nx = meta->gds.Nx;
1312 457 : Ny = meta->gds.Ny;
1313 457 : x1 = 1;
1314 457 : x2 = Nx;
1315 457 : y1 = 1;
1316 457 : y2 = Ny;
1317 457 : f_subGrid = 0;
1318 :
1319 : #ifdef deadcode
1320 : /* Figure out if we need iain or ain, and set it to Grib_Data. At the
1321 : * same time handle any bitmaps, and compute some statistics. */
1322 : if ((f_subGrid) && (meta->gds.scan != 64)) {
1323 : errSprintf ("Can not do a subgrid of non scanmode 64 grid yet.\n");
1324 : return -3;
1325 : }
1326 : #endif
1327 :
1328 457 : if (Grib_Data != nullptr && strcmp (meta->element, "Wx") != 0) {
1329 157 : if (strcmp (meta->element, "WWA") != 0) {
1330 157 : ParseGrid (fp, &(meta->gridAttrib), Grib_Data, grib_DataLen, Nx, Ny,
1331 157 : meta->gds.scan, IS->nd2x3, IS->iain, ibitmap, IS->ib, unitM, unitB, 0,
1332 : 0, nullptr, f_subGrid, x1, y1, x2, y2);
1333 : } else {
1334 0 : ParseGrid (fp, &(meta->gridAttrib), Grib_Data, grib_DataLen, Nx, Ny,
1335 0 : meta->gds.scan, IS->nd2x3, IS->iain, ibitmap, IS->ib, unitM, unitB, 1,
1336 : meta->pds2.sect2.hazard.dataLen, meta->pds2.sect2.hazard.f_valid, f_subGrid, x1, y1,
1337 : x2, y2);
1338 : /* compact the table to only those which are actually used. */
1339 0 : cnt = 0;
1340 0 : for (i = 0; i < meta->pds2.sect2.hazard.dataLen; i++) {
1341 0 : if (meta->pds2.sect2.hazard.f_valid[i] == 2) {
1342 0 : meta->pds2.sect2.hazard.haz[i].validIndex = cnt;
1343 0 : cnt++;
1344 0 : } else if (meta->pds2.sect2.hazard.f_valid[i] == 3) {
1345 0 : meta->pds2.sect2.hazard.f_valid[i] = 0;
1346 0 : meta->pds2.sect2.hazard.haz[i].validIndex = cnt;
1347 0 : cnt++;
1348 : } else {
1349 0 : meta->pds2.sect2.hazard.haz[i].validIndex = -1;
1350 : }
1351 : }
1352 : }
1353 300 : } else if( Grib_Data != nullptr ) {
1354 : /* Handle weather grid. ParseGrid looks up the values... If they are
1355 : * "<Invalid>" it sets it to missing (or creates one). If the table
1356 : * entry is used it sets f_valid to 2. */
1357 0 : ParseGrid (fp, &(meta->gridAttrib), Grib_Data, grib_DataLen, Nx, Ny,
1358 0 : meta->gds.scan, IS->nd2x3, IS->iain, ibitmap, IS->ib, unitM, unitB, 1,
1359 : meta->pds2.sect2.wx.dataLen, meta->pds2.sect2.wx.f_valid, f_subGrid, x1, y1,
1360 : x2, y2);
1361 0 : if( *Grib_Data == nullptr )
1362 0 : return -1;
1363 :
1364 : /* compact the table to only those which are actually used. */
1365 0 : cnt = 0;
1366 0 : for (i = 0; i < meta->pds2.sect2.wx.dataLen; i++) {
1367 0 : if (meta->pds2.sect2.wx.f_valid[i] == 2) {
1368 0 : meta->pds2.sect2.wx.ugly[i].validIndex = cnt;
1369 0 : cnt++;
1370 0 : } else if (meta->pds2.sect2.wx.f_valid[i] == 3) {
1371 0 : meta->pds2.sect2.wx.f_valid[i] = 0;
1372 0 : meta->pds2.sect2.wx.ugly[i].validIndex = cnt;
1373 0 : cnt++;
1374 : } else {
1375 0 : meta->pds2.sect2.wx.ugly[i].validIndex = -1;
1376 : }
1377 : }
1378 : } else {
1379 : // Simulate part of ParseGrid() "Resolve bitmap (if there is one) in the data"
1380 : // behavior just to set nodata value.
1381 300 : if (ibitmap && (meta->gridAttrib.f_miss != 1) && (meta->gridAttrib.f_miss != 2)) {
1382 6 : meta->gridAttrib.f_miss = 1;
1383 6 : meta->gridAttrib.missPri = 9999;
1384 : }
1385 : }
1386 :
1387 : /* Figure out some other non-section oriented meta data. */
1388 : /* strftime (meta->refTime, 20, "%Y%m%d%H%M",
1389 : gmtime (&(meta->pds2.refTime)));
1390 : */
1391 457 : Clock_Print (meta->refTime, 20, meta->pds2.refTime, "%Y%m%d%H%M", 0);
1392 : /*
1393 : strftime (meta->validTime, 20, "%Y%m%d%H%M",
1394 : gmtime (&(meta->pds2.sect4.validTime)));
1395 : */
1396 457 : Clock_Print (meta->validTime, 20, meta->pds2.sect4.validTime,
1397 : "%Y%m%d%H%M", 0);
1398 :
1399 457 : const double deltTime = meta->pds2.sect4.validTime - meta->pds2.refTime;
1400 914 : if (deltTime >= std::numeric_limits<sInt4>::max() ||
1401 457 : deltTime <= std::numeric_limits<sInt4>::min()) {
1402 0 : meta->deltTime = 0;
1403 0 : preErrSprintf ("deltTime over range\n");
1404 0 : return -4;
1405 : }
1406 457 : meta->deltTime = static_cast<sInt4>(deltTime);
1407 :
1408 457 : return 0;
1409 : }
1410 :
1411 : #if 0 // unused by GDAL
1412 : int ReadGrib2RecordFast (FILE *fp, sChar f_unit, double **Grib_Data,
1413 : uInt4 *grib_DataLen, grib_MetaData *meta,
1414 : IS_dataType *IS, int subgNum, double majEarth,
1415 : double minEarth, int simpVer, int simpWWA,
1416 : sInt4 *f_endMsg, LatLon *lwlf, LatLon *uprt)
1417 : {
1418 : sInt4 l3264b; /* Number of bits in a sInt4. Needed by FORTRAN
1419 : * unpack library to determine if system has a 4
1420 : * byte_ sInt4 or an 8 byte sInt4. */
1421 : char *buff; /* Holds the info between records. */
1422 : uInt4 buffLen; /* Length of info between records. */
1423 : sInt4 sect0[SECT0LEN_WORD]; /* Holds the current Section 0. */
1424 : uInt4 gribLen; /* Length of the current GRIB message. */
1425 : sInt4 nd5; /* Size of grib message rounded up to the nearest
1426 : * sInt4. */
1427 : unsigned char *c_ipack; /* A char ptr to the message stored in IS->ipack */
1428 : sInt4 local_ns[8]; /* Local copy of section lengths. */
1429 : sInt4 nd2x3; /* Total number of grid points. */
1430 : short int table50; /* Type of packing used. (See code table 5.0)
1431 : * (GS5_SIMPLE==0, GS5_CMPLX==2, GS5_CMPLXSEC==3) */
1432 : sInt4 nidat; /* Size of section 2 if it contains integer data. */
1433 : sInt4 nrdat; /* Size of section 2 if it contains float data. */
1434 : sInt4 inew; /* 1 if this is the first grid we are reading. 0 if
1435 : * this is the second or later grid from the same
1436 : * GRIB message. */
1437 : sInt4 iclean = 0; /* 0 embed the missing values, 1 don't. */
1438 : int j; /* Counter used to find the desired subgrid. */
1439 : sInt4 kfildo = 5; /* FORTRAN Unit number for diagnostic info. Ignored,
1440 : * unless library is compiled a particular way. */
1441 : sInt4 ibitmap; /* 0 means no bitmap returned, otherwise 1. */
1442 : float xmissp; /* The primary missing value. If iclean = 0, this
1443 : * value is embedded in grid, otherwise it is the
1444 : * value returned from the GRIB message. */
1445 : float xmisss; /* The secondary missing value. If iclean = 0, this
1446 : * value is embedded in grid, otherwise it is the
1447 : * value returned from the GRIB message. */
1448 : sInt4 jer[UNPK_NUM_ERRORS * 2]; /* Any Error codes along with their *
1449 : * severity levels generated using the *
1450 : * unpack GRIB2 library. */
1451 : sInt4 ndjer = UNPK_NUM_ERRORS; /* The number of rows in JER( ). */
1452 : sInt4 kjer; /* The actual number of errors returned in JER. */
1453 : size_t i; /* counter as we loop through jer. */
1454 : double unitM, unitB; /* values in y = m x + b used for unit conversion. */
1455 : char unitName[15]; /* Holds the string name of the current unit. */
1456 : int unitLen; /* String length of string name of current unit. */
1457 : int version; /* Which version of GRIB is in this message. */
1458 : gdsType newGds; /* The GDS of the subgrid if needed. */
1459 : int x1, y1; /* The original grid coordinates of the lower left
1460 : * corner of the subgrid. */
1461 : int x2, y2; /* The original grid coordinates of the upper right
1462 : * corner of the subgrid. */
1463 : uChar f_subGrid; /* True if we have a subgrid. */
1464 : sInt4 Nx, Ny; /* original size of the data. */
1465 :
1466 : /*
1467 : * f_endMsg is 1 if in the past we either completed reading a message,
1468 : * or we haven't read any messages. In either case we need to read the
1469 : * next message from file.
1470 : * If f_endMsg is false, then there is more to read from IS->ipack, so we
1471 : * don't want to throw it out, nor have to re-read ipack from disk.
1472 : */
1473 : l3264b = sizeof (sInt4) * 8;
1474 : buff = NULL;
1475 : buffLen = 0;
1476 : if (*f_endMsg == 1) {
1477 : if (ReadSECT0 (fp, &buff, &buffLen, -1, sect0, &gribLen, &version) < 0) {
1478 : preErrSprintf ("Inside ReadGrib2Record\n");
1479 : free (buff);
1480 : return -1;
1481 : }
1482 : meta->GribVersion = version;
1483 : if (version != 2) {
1484 : printf ("Fast parsing doesn't handle this version because ReadGrib1Record/ReadTDLPRecord used Grib_Data[]\n");
1485 : return -1;
1486 : }
1487 :
1488 : /*
1489 : * Make room for entire message, and read it in.
1490 : */
1491 : /* nd5 needs to be gribLen in (sInt4) units rounded up. */
1492 : nd5 = (gribLen + 3) / 4;
1493 : if (nd5 > IS->ipackLen) {
1494 : IS->ipackLen = nd5;
1495 : IS->ipack = (sInt4 *) realloc ((void *) (IS->ipack),
1496 : (IS->ipackLen) * sizeof (sInt4));
1497 : }
1498 : c_ipack = (unsigned char *) IS->ipack;
1499 : /* Init last sInt4 to 0, to make sure that the padded bytes are 0. */
1500 : IS->ipack[nd5 - 1] = 0;
1501 : /* Init first 4 sInt4 to sect0. */
1502 : memcpy (c_ipack, sect0, SECT0LEN_WORD * 4);
1503 : /* Read in the rest of the message. */
1504 : if (fread (c_ipack + SECT0LEN_WORD * 4, sizeof (char),
1505 : (gribLen - SECT0LEN_WORD * 4),
1506 : fp) != (gribLen - SECT0LEN_WORD * 4)) {
1507 : errSprintf ("GribLen = %ld, SECT0Len_WORD = %d\n", gribLen,
1508 : SECT0LEN_WORD);
1509 : errSprintf ("Ran out of file\n");
1510 : free (buff);
1511 : return -1;
1512 : }
1513 :
1514 : /*
1515 : * Make sure the arrays are large enough for call to unpacker library.
1516 : */
1517 : /* FindSectLen Does not want (ipack / c_ipack) word swapped, because
1518 : * that would make it much more confusing to find bytes in c_ipack. */
1519 : if (FindSectLen (c_ipack, gribLen, local_ns, &nd2x3, &table50) < 0) {
1520 : preErrSprintf ("Inside ReadGrib2Record.. Calling FindSectLen\n");
1521 : free (buff);
1522 : return -2;
1523 : }
1524 :
1525 : /* Make sure all 'is' arrays except ns[7] are MAX (IS.ns[] ,
1526 : * local_ns[]). See note 1 for reason to exclude ns[7] from MAX (). */
1527 : for (i = 0; i < 7; i++) {
1528 : if (local_ns[i] > IS->ns[i]) {
1529 : IS->ns[i] = local_ns[i];
1530 : IS->is[i] = (sInt4 *) realloc ((void *) (IS->is[i]),
1531 : IS->ns[i] * sizeof (sInt4));
1532 : }
1533 : }
1534 :
1535 : /* Allocate room for sect 2. If local_ns[2] = -1 there is no sect 2. */
1536 : if (local_ns[2] == -1) {
1537 : nidat = 10;
1538 : nrdat = 10;
1539 : } else {
1540 : /*
1541 : * See note 2) We have a section 2, so use:
1542 : * MAX (32 * local_ns[2],SECT2_INTSIZE)
1543 : * and MAX (32 * local_ns[2],SECT2_FLOATSIZE)
1544 : * for size of section 2 unpacked.
1545 : */
1546 : nidat = (32 * local_ns[2] < SECT2_INIT_SIZE) ? SECT2_INIT_SIZE :
1547 : 32 * local_ns[2];
1548 : nrdat = nidat;
1549 : }
1550 : if (nidat > IS->nidat) {
1551 : IS->nidat = nidat;
1552 : IS->idat = (sInt4 *) realloc ((void *) IS->idat,
1553 : IS->nidat * sizeof (sInt4));
1554 : }
1555 : if (nrdat > IS->nrdat) {
1556 : IS->nrdat = nrdat;
1557 : IS->rdat = (float *) realloc ((void *) IS->rdat,
1558 : IS->nrdat * sizeof (float));
1559 : }
1560 : /* Make sure we have room for the GRID part of the output. */
1561 : if (nd2x3 > IS->nd2x3) {
1562 : IS->nd2x3 = nd2x3;
1563 : IS->iain = (sInt4 *) realloc ((void *) IS->iain,
1564 : IS->nd2x3 * sizeof (sInt4));
1565 : IS->ib = (sInt4 *) realloc ((void *) IS->ib,
1566 : IS->nd2x3 * sizeof (sInt4));
1567 : }
1568 : /* See note 3) If table50 == 3, unpacker library needs nd5 >= nd2x3. */
1569 : if ((table50 == 3) || (table50 == 0)) {
1570 : if (nd5 < nd2x3) {
1571 : nd5 = nd2x3;
1572 : if (nd5 > IS->ipackLen) {
1573 : IS->ipackLen = nd5;
1574 : IS->ipack = (sInt4 *) realloc ((void *) (IS->ipack),
1575 : IS->ipackLen * sizeof (sInt4));
1576 : }
1577 : /* Don't need to do the following, but we do in case code
1578 : * changes. */
1579 : c_ipack = (unsigned char *) IS->ipack;
1580 : }
1581 : }
1582 : IS->nd5 = nd5;
1583 : /* Unpacker library requires ipack to be MSB. */
1584 : /*
1585 : #ifdef DEBUG
1586 : if (1==1) {
1587 : FILE *fp = fopen ("test.bin", "wb");
1588 : fwrite (IS->ipack, sizeof (sInt4), IS->nd5, fp);
1589 : fclose (fp);
1590 : }
1591 : #endif
1592 : */
1593 : } else {
1594 : c_ipack = (unsigned char *) IS->ipack;
1595 : /* GRIB2 files are in big endian so c_ipack is as well. */
1596 : #ifdef LITTLE_ENDIAN
1597 : revmemcpy (&gribLen, &(c_ipack[12]), sizeof (sInt4));
1598 : #else
1599 : memcpy (&gribLen, &(c_ipack[12]), sizeof (sInt4));
1600 : #endif
1601 : }
1602 : free (buff);
1603 :
1604 : /* Loop through the grib message looking for the subgNum grid. subgNum
1605 : * goes from 0 to n-1. */
1606 : for (j = 0; j <= subgNum; j++) {
1607 : if (j == 0) {
1608 : inew = 1;
1609 : } else {
1610 : inew = 0;
1611 : }
1612 :
1613 : /* Note we are getting data back either as a float or an int, but not
1614 : * both, so we don't need to allocated room for both. */
1615 : /*
1616 : unpk_grib2 (&kfildo, (float *) (IS->iain), IS->iain, &(IS->nd2x3),
1617 : IS->idat, &(IS->nidat), IS->rdat, &(IS->nrdat), IS->is[0],
1618 : &(IS->ns[0]), IS->is[1], &(IS->ns[1]), IS->is[2],
1619 : &(IS->ns[2]), IS->is[3], &(IS->ns[3]), IS->is[4],
1620 : &(IS->ns[4]), IS->is[5], &(IS->ns[5]), IS->is[6],
1621 : &(IS->ns[6]), IS->is[7], &(IS->ns[7]), IS->ib, &ibitmap,
1622 : IS->ipack, &(IS->nd5), &xmissp, &xmisss, &inew, &iclean,
1623 : &l3264b, f_endMsg, jer, &ndjer, &kjer);
1624 : */
1625 : c_ipack = (unsigned char *)IS->ipack;
1626 : unpk_g2ncep(&kfildo, (float *) (IS->iain), IS->iain, &(IS->nd2x3),
1627 : IS->idat, &(IS->nidat), IS->rdat, &(IS->nrdat), IS->is[0],
1628 : &(IS->ns[0]), IS->is[1], &(IS->ns[1]), IS->is[2],
1629 : &(IS->ns[2]), IS->is[3], &(IS->ns[3]), IS->is[4],
1630 : &(IS->ns[4]), IS->is[5], &(IS->ns[5]), IS->is[6],
1631 : &(IS->ns[6]), IS->is[7], &(IS->ns[7]), IS->ib, &ibitmap,
1632 : c_ipack, &(IS->nd5), &xmissp, &xmisss, &inew, &iclean,
1633 : &l3264b, f_endMsg, jer, &ndjer, &kjer);
1634 :
1635 :
1636 : /*
1637 : * Check for error messages...
1638 : * If we get an error message, print it, and return.
1639 : */
1640 : for (i = 0; i < (uInt4) kjer; i++) {
1641 : if (jer[ndjer + i] == 0) {
1642 : /* no error. */
1643 : } else if (jer[ndjer + i] == 1) {
1644 : /* Warning. */
1645 : #ifdef DEBUG
1646 : printf ("Warning: Unpack library warning code (%ld %ld)\n",
1647 : jer[i], jer[ndjer + i]);
1648 : #endif
1649 : } else {
1650 : /* BAD Error. */
1651 : errSprintf ("ERROR: Unpack library error code (%ld %ld)\n",
1652 : jer[i], jer[ndjer + i]);
1653 : return -3;
1654 : }
1655 : }
1656 : }
1657 :
1658 : /* Parse the meta data out. */
1659 : if (MetaParse (meta, IS->is[0], IS->ns[0], IS->is[1], IS->ns[1],
1660 : IS->is[2], IS->ns[2], IS->rdat, IS->nrdat, IS->idat,
1661 : IS->nidat, IS->is[3], IS->ns[3], IS->is[4], IS->ns[4],
1662 : IS->is[5], IS->ns[5], gribLen, xmissp, xmisss, simpVer, simpWWA)
1663 : != 0) {
1664 : #ifdef DEBUG
1665 : FILE *fp;
1666 : if ((fp = fopen ("dump.is0", "wt")) != NULL) {
1667 : for (i = 0; i < 8; i++) {
1668 : fprintf (fp, "---Section %d---\n", i);
1669 : for (j = 1; j <= IS->ns[i]; j++) {
1670 : fprintf (fp, "IS%d Item %d = %ld\n", i, j, IS->is[i][j - 1]);
1671 : }
1672 : }
1673 : fclose (fp);
1674 : }
1675 : #endif
1676 : preErrSprintf ("Inside ReadGrib2Record.. Problems in MetaParse\n");
1677 : return -4;
1678 : }
1679 :
1680 : if ((majEarth > 6000) && (majEarth < 7000)) {
1681 : if ((minEarth > 6000) && (minEarth < 7000)) {
1682 : meta->gds.f_sphere = 0;
1683 : meta->gds.majEarth = majEarth;
1684 : meta->gds.minEarth = minEarth;
1685 : } else {
1686 : meta->gds.f_sphere = 1;
1687 : meta->gds.majEarth = majEarth;
1688 : meta->gds.minEarth = majEarth;
1689 : }
1690 : }
1691 :
1692 : /* Figure out an equation to pass to ParseGrid to convert the units for
1693 : * this grid. */
1694 : /*
1695 : if (ComputeUnit (meta->pds2.prodType, meta->pds2.sect4.templat,
1696 : meta->pds2.sect4.cat, meta->pds2.sect4.subcat, f_unit,
1697 : &unitM, &unitB, unitName) == 0) {
1698 : */
1699 : if (ComputeUnit (meta->convert, meta->unitName, f_unit, &unitM, &unitB,
1700 : unitName) == 0) {
1701 : unitLen = strlen (unitName);
1702 : meta->unitName = (char *) realloc ((void *) (meta->unitName),
1703 : 1 + unitLen * sizeof (char));
1704 : strncpy (meta->unitName, unitName, unitLen);
1705 : meta->unitName[unitLen] = '\0';
1706 : }
1707 :
1708 : /* compute the subgrid. */
1709 : if ((lwlf->lat != -100) && (uprt->lat != -100)) {
1710 : Nx = meta->gds.Nx;
1711 : Ny = meta->gds.Ny;
1712 : if (computeSubGrid (lwlf, &x1, &y1, uprt, &x2, &y2, &(meta->gds),
1713 : &newGds) != 0) {
1714 : preErrSprintf ("ERROR: In compute subgrid.\n");
1715 : return 1;
1716 : }
1717 : /* I couldn't decide if I should "permanently" change the GDS or not.
1718 : * when I wrote computeSubGrid. If next line stays, really should
1719 : * rewrite computeSubGrid. */
1720 : memcpy (&(meta->gds), &newGds, sizeof (gdsType));
1721 : f_subGrid = 1;
1722 : } else {
1723 : Nx = meta->gds.Nx;
1724 : Ny = meta->gds.Ny;
1725 : x1 = 1;
1726 : x2 = Nx;
1727 : y1 = 1;
1728 : y2 = Ny;
1729 : f_subGrid = 0;
1730 : }
1731 :
1732 : if ((f_subGrid) && (meta->gds.scan != 64)) {
1733 : errSprintf ("Can not do a subgrid of non scanmode 64 grid yet.\n");
1734 : return -3;
1735 : }
1736 :
1737 : /* Figure out some other non-section oriented meta data. */
1738 : /* strftime (meta->refTime, 20, "%Y%m%d%H%M",
1739 : gmtime (&(meta->pds2.refTime)));
1740 : */
1741 : Clock_Print (meta->refTime, 20, meta->pds2.refTime, "%Y%m%d%H%M", 0);
1742 : /*
1743 : strftime (meta->validTime, 20, "%Y%m%d%H%M",
1744 : gmtime (&(meta->pds2.sect4.validTime)));
1745 : */
1746 : Clock_Print (meta->validTime, 20, meta->pds2.sect4.validTime,
1747 : "%Y%m%d%H%M", 0);
1748 :
1749 : meta->deltTime = (sInt4) (meta->pds2.sect4.validTime - meta->pds2.refTime);
1750 :
1751 : return 0;
1752 : }
1753 : #endif
|
