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 : #include "myassert.h"
22 : #include "myerror.h"
23 : #include "memendian.h"
24 : #include "meta.h"
25 : #include "metaname.h"
26 : //#include "write.h"
27 : #include "degrib2.h"
28 : #include "degrib1.h"
29 : #include "tdlpack.h"
30 : #include "grib2api.h"
31 : //#include "mymapf.h"
32 : #include "clock.h"
33 :
34 : #define GRIB_UNSIGN_INT3(a,b,c) ((a<<16)+(b<<8)+c)
35 :
36 : /*****************************************************************************
37 : * ReadSect0() -- Review 12/2002
38 : *
39 : * Arthur Taylor / MDL
40 : *
41 : * PURPOSE
42 : * Looks for the next GRIB message, by looking for the keyword "GRIB". It
43 : * expects the message in "expect" bytes from the start, but could find the
44 : * message in "expect2" bytes or 0 bytes from the start. Returns -1 if it
45 : * can't find "GRIB", 1 if "GRIB" is not 0, "expect", or "expect2" bytes from
46 : * the start.
47 : * It stores the bytes it reads (a max of "expect") upto but not including
48 : * the 'G' in "GRIB" in wmo.
49 : *
50 : * After it finds section 0, it then parses the 16 bytes that make up
51 : * section 0 so that it can return the length of the entire GRIB message.
52 : *
53 : * When done, it sets fp to point to the end of Sect0.
54 : *
55 : * The reason for this procedure is so that we can read in the size of the
56 : * grib message, and thus allocate enough memory to read the message in before
57 : * making it Big endian, and passing it to the library for unpacking.
58 : *
59 : * ARGUMENTS
60 : * fp = A pointer to an opened file in which to read.
61 : * When done, this points to the start of section 1. (Input/Output)
62 : * buff = The data between messages. (Input/Output)
63 : * buffLen = The length of buff (Output)
64 : * limit = How many bytes to read before giving up and stating it is not
65 : * a proper message. (-1 means no limit). (Input)
66 : * sect0 = The read in Section 0 (as seen on disk). (Output)
67 : * gribLen = Length of this GRIB message. (Output)
68 : * version = 1 if GRIB1 message, 2 if GRIB2 message, -1 if TDLP message.
69 : * (Output)
70 : * expect = The expected number of bytes to find "GRIB" in. (Input)
71 : * expect2 = The second possible number of bytes to find "GRIB" in. (Input)
72 : * wmo = Assumed allocated to be at least size "expect".
73 : * Holds the bytes before the first "GRIB" message.
74 : * expect should be > expect2, but is up to caller (Output)
75 : * wmoLen = Length of wmo (total number of bytes read - SECT0LEN_WORD * 4).
76 : * (Output)
77 : *
78 : * FILES/DATABASES:
79 : * An already opened "GRIB2" File
80 : *
81 : * RETURNS: int (could use errSprintf())
82 : * 1 = Length of wmo was != 0 and was != expect
83 : * 0 = OK
84 : * -1 = Couldn't find "GRIB" part of message.
85 : * -2 = Ran out of file while reading this section.
86 : * -3 = Grib version was not 1 or 2.
87 : * -4 = Most significant sInt4 of GRIB length was not 0
88 : * -5 = Grib message length was <= 16 (can't be smaller than just sect 0)
89 : *
90 : * HISTORY
91 : * 9/2002 Arthur Taylor (MDL/RSIS): Created.
92 : * 11/2002 AAT: Combined with ReadWMOHeader
93 : * 12/2002 (TK,AC,TB,&MS): Code Review.
94 : * 1/2003 AAT: Bug found. wmo access out of bounds of expect when setting
95 : * the /0 element, if wmoLen > expect.
96 : * 4/2003 AAT: Added ability to handle GRIB version 1.
97 : * 5/2003 AAT: Added limit option.
98 : * 8/2003 AAT: Removed dependence on offset, and fileLen.
99 : * 10/2004 AAT: Modified to allow for TDLP files
100 : *
101 : * NOTES
102 : * 1a) 1196575042L == ASCII representation of "GRIB" (GRIB in MSB)
103 : * 1b) 1112101447L == ASCII representation of "BIRG" (GRIB in LSB)
104 : * 1c) 1413762128L == ASCII representation of "TDLP" (TDLP in MSB)
105 : * 1d) 1347175508L == ASCII representation of "PLDT" (TDLP in LSB)
106 : * 2) Takes advantage of the wordType to check that the edition is correct.
107 : * 3) May want to return prodType.
108 : * 4) WMO_HEADER_ORIG_LEN was added for backward compatibility... should be
109 : * removed when we no longer use old format. (say in a year from 11/2002)
110 : *
111 : *****************************************************************************
112 : */
113 27 : int ReadSECT0 (DataSource &fp, char **buff, uInt4 *buffLen, sInt4 limit,
114 : sInt4 sect0[SECT0LEN_WORD], uInt4 *gribLen, int *version)
115 : {
116 : typedef union {
117 : sInt4 li;
118 : unsigned char buffer[4];
119 : } wordType;
120 :
121 27 : uChar gribMatch = 0; /* Counts how many letters in GRIB we've matched. */
122 27 : uChar tdlpMatch = 0; /* Counts how many letters in TDLP we've matched. */
123 : wordType word; /* Used to check that the edition is correct. */
124 : uInt4 curLen; /* Where we currently are in buff. */
125 : uInt4 i; /* Used to loop over the first few char's */
126 : uInt4 stillNeed; /* Number of bytes still needed to get 1st 8 bytes of
127 : * message into memory. */
128 :
129 : /* Get first 8 bytes. If GRIB we don't care. If TDLP, this is the length
130 : * of record. Read at least 1 record (length + 2 * 8) + 8 (next record
131 : * length) + 8 bytes before giving up. */
132 27 : curLen = 8;
133 27 : if (*buffLen < curLen) {
134 26 : *buffLen = curLen;
135 26 : *buff = (char *) realloc ((void *) *buff, *buffLen * sizeof (char));
136 : }
137 27 : if (fp.DataSourceFread(*buff, sizeof (char), curLen) != curLen) {
138 0 : errSprintf ("ERROR: Couldn't find 'GRIB' or 'TDLP'\n");
139 0 : return -1;
140 : }
141 : /*
142 : Can't do the following because we don't know if the file is a GRIB file or
143 : not, or if it was a FORTRAN file.
144 : if (limit > 0) {
145 : MEMCPY_BIG (&recLen, *buff, 4);
146 : limit = (limit > recLen + 32) ? limit : recLen + 32;
147 : }
148 : */
149 145 : while ((tdlpMatch != 4) && (gribMatch != 4)) {
150 411 : for (i = curLen - 8; i + 3 < curLen; i++) {
151 347 : if ((*buff)[i] == 'G') {
152 54 : if (((*buff)[i + 1] == 'R') && ((*buff)[i + 2] == 'I') &&
153 27 : ((*buff)[i + 3] == 'B')) {
154 27 : gribMatch = 4;
155 27 : break;
156 : }
157 320 : } else if ((*buff)[i] == 'T') {
158 0 : if (((*buff)[i + 1] == 'D') && ((*buff)[i + 2] == 'L') &&
159 0 : ((*buff)[i + 3] == 'P')) {
160 0 : tdlpMatch = 4;
161 0 : break;
162 : }
163 : }
164 : }
165 91 : stillNeed = i - (curLen - 8);
166 : /* Read enough of message to have the first 8 bytes (including ID). */
167 91 : if (stillNeed != 0) {
168 64 : curLen += stillNeed;
169 64 : if ((limit >= 0) && (curLen > (size_t) limit)) {
170 0 : errSprintf ("ERROR: Couldn't find type in %ld bytes\n", limit);
171 0 : return -1;
172 : }
173 64 : if (*buffLen < curLen) {
174 56 : *buffLen = curLen;
175 : *buff = (char *) realloc ((void *) *buff,
176 56 : *buffLen * sizeof (char));
177 : }
178 64 : if (fp.DataSourceFread((*buff) + (curLen - stillNeed), sizeof (char), stillNeed) != stillNeed) {
179 0 : errSprintf ("ERROR: Ran out of file reading SECT0\n");
180 0 : return -1;
181 : }
182 : }
183 : }
184 :
185 : /* curLen and (*buff) hold 8 bytes of section 0. */
186 27 : curLen -= 8;
187 27 : memcpy (&(sect0[0]), (*buff) + curLen, 4);
188 : #ifdef DEBUG
189 : #ifdef LITTLE_ENDIAN
190 : myAssert ((sect0[0] == 1112101447L) || (sect0[0] == 1347175508L));
191 : #else
192 : myAssert ((sect0[0] == 1196575042L) || (sect0[0] == 1413762128L));
193 : #endif
194 : #endif
195 27 : memcpy (&(sect0[1]), *buff + curLen + 4, 4);
196 : /* Make sure we don't pass back part of "GRIB" in the buffer. */
197 27 : (*buff)[curLen] = '\0';
198 27 : *buffLen = curLen;
199 :
200 27 : word.li = sect0[1];
201 27 : if (tdlpMatch == 4) {
202 0 : if (word.buffer[3] != 0) {
203 0 : errSprintf ("ERROR: unexpected version of TDLP in SECT0\n");
204 0 : return -2;
205 : }
206 0 : *version = -1;
207 : /* Find out the GRIB Message Length */
208 0 : *gribLen = GRIB_UNSIGN_INT3 (word.buffer[0], word.buffer[1],
209 0 : word.buffer[2]);
210 : /* Min message size: GRIB1=52, TDLP=59, GRIB2=86. */
211 0 : if (*gribLen < 59) {
212 0 : errSprintf ("TDLP length %ld was < 59?\n", *gribLen);
213 0 : return -5;
214 : }
215 27 : } else if (word.buffer[3] == 1) {
216 22 : *version = 1;
217 : /* Find out the GRIB Message Length */
218 66 : *gribLen = GRIB_UNSIGN_INT3 (word.buffer[0], word.buffer[1],
219 66 : word.buffer[2]);
220 : /* Min message size: GRIB1=52, TDLP=59, GRIB2=86. */
221 22 : if (*gribLen < 52) {
222 0 : errSprintf ("GRIB1 length %ld was < 52?\n", *gribLen);
223 0 : return -5;
224 : }
225 5 : } else if (word.buffer[3] == 2) {
226 5 : *version = 2;
227 : /* Make sure we still have enough file for the rest of section 0. */
228 5 : if (fp.DataSourceFread(sect0 + 2, sizeof (sInt4), 2) != 2) {
229 0 : errSprintf ("ERROR: Ran out of file reading SECT0\n");
230 0 : return -2;
231 : }
232 5 : if (sect0[2] != 0) {
233 0 : errSprintf ("Most significant sInt4 of GRIB length was not 0?\n");
234 : errSprintf ("This is either an error, or we have a single GRIB "
235 : "message which is larger than 2^31 = 2,147,283,648 "
236 0 : "bytes.\n");
237 0 : return -4;
238 : }
239 : #ifdef LITTLE_ENDIAN
240 5 : revmemcpy (gribLen, &(sect0[3]), sizeof (sInt4));
241 : #else
242 : *gribLen = sect0[3];
243 : #endif
244 : } else {
245 0 : errSprintf ("ERROR: Not TDLPack, and Grib edition is not 1 or 2\n");
246 0 : return -3;
247 : }
248 27 : return 0;
249 : }
250 :
251 : /*****************************************************************************
252 : * FindGRIBMsg() -- Review 12/2002
253 : *
254 : * Arthur Taylor / MDL
255 : *
256 : * PURPOSE
257 : * Jumps through a GRIB2 file looking for a specific message. Currently
258 : * that message is determined by msgNum which is in the range of 1..n.
259 : * In the future we may be searching based on projection or date.
260 : *
261 : * ARGUMENTS
262 : * fp = The current GRIB2 file to look through. (Input)
263 : * msgNum = Which message to look for. (Input)
264 : * offset = Where in the file the message starts (this is before the
265 : * wmo ASCII part if there is one.) (Output)
266 : * curMsg = The current # of messages we have looked through. (In/Out)
267 : *
268 : * FILES/DATABASES:
269 : * An already opened "GRIB2" File
270 : *
271 : * RETURNS: int (could use errSprintf())
272 : * 0 = OK
273 : * -1 = Problems reading Section 0.
274 : * -2 = Ran out of file.
275 : *
276 : * HISTORY
277 : * 11/2002 Arthur Taylor (MDL/RSIS): Created.
278 : * 12/2002 (TK,AC,TB,&MS): Code Review.
279 : * 6/2003 Matthew T. Kallio (matt@wunderground.com):
280 : * "wmo" dimension increased to WMO_HEADER_LEN + 1 (for '\0' char)
281 : * 8/2003 AAT: Removed dependence on offset and fileLen.
282 : *
283 : * NOTES
284 : *****************************************************************************
285 : */
286 0 : int FindGRIBMsg (DataSource &fp, int msgNum, sInt4 *offset, int *curMsg)
287 : {
288 : int cnt; /* The current message we are looking at. */
289 : char *buff; /* Holds the info between records. */
290 : uInt4 buffLen; /* Length of info between records. */
291 : sInt4 sect0[SECT0LEN_WORD]; /* Holds the current Section 0. */
292 : uInt4 gribLen; /* Length of the current GRIB message. */
293 : int version; /* Which version of GRIB is in this message. */
294 : int c; /* Determine if end of the file without fileLen. */
295 : sInt4 jump; /* How far to jump to get to past GRIB message. */
296 :
297 0 : cnt = *curMsg + 1;
298 0 : buff = NULL;
299 0 : buffLen = 0;
300 0 : while ((c = fp.DataSourceFgetc()) != EOF) {
301 0 : fp.DataSourceUngetc(c);
302 0 : if (cnt >= msgNum) {
303 : /* 12/1/2004 version 1.63 forgot to free buff */
304 0 : free (buff);
305 0 : *curMsg = cnt;
306 0 : return 0;
307 : }
308 : /* Read section 0 to find gribLen and wmoLen. */
309 0 : if (ReadSECT0 (fp, &buff, &buffLen, GRIB_LIMIT, sect0, &gribLen,
310 : &version) < 0) {
311 0 : preErrSprintf ("Inside FindGRIBMsg\n");
312 0 : free (buff);
313 0 : return -1;
314 : }
315 : myAssert ((version == 1) || (version == 2) || (version == -1));
316 : /* Continue on to the next grib message. */
317 0 : if ((version == 1) || (version == -1)) {
318 0 : jump = gribLen - 8;
319 : } else {
320 0 : jump = gribLen - 16;
321 : }
322 0 : fp.DataSourceFseek(jump, SEEK_CUR);
323 0 : *offset = *offset + gribLen + buffLen;
324 0 : cnt++;
325 : }
326 0 : free (buff);
327 0 : *curMsg = cnt - 1;
328 : /* Return -2 since we reached the end of file. This may not be an error
329 : * (multiple file option). */
330 0 : return -2;
331 : /*
332 : errSprintf ("ERROR: Ran out of file looking for msgNum %d.\n", msgNum);
333 : errSprintf (" Current msgNum %d\n", cnt);
334 : */
335 : }
336 :
337 : /*****************************************************************************
338 : * FindSectLen2to7() --
339 : *
340 : * Arthur Taylor / MDL
341 : *
342 : * PURPOSE
343 : * Looks through a GRIB message and finds out the maximum size of each
344 : * section. Simpler if there is only one grid in the message.
345 : *
346 : * ARGUMENTS
347 : * c_ipack = The current GRIB2 message. (Input)
348 : * gribLen = Length of c_ipack. (Input)
349 : * ns = Array of section lengths. (Output)
350 : * sectNum = Which section to start with. (Input)
351 : * curTot = on going total read from c_ipack. (Input)
352 : * nd2x3 = Total number of grid points (Output)
353 : * table50 = Type of packing used. (See code table 5.0) (GS5_SIMPLE = 0,
354 : * GS5_CMPLX = 2, GS5_CMPLXSEC = 3) (Output)
355 : *
356 : * FILES/DATABASES: None
357 : *
358 : * RETURNS: int (could use errSprintf())
359 : * 0 = OK
360 : * -1 = Ran of data in a section
361 : * -2 = Section not properly labeled.
362 : *
363 : * HISTORY
364 : * 3/2003 AAT: Created
365 : *
366 : * NOTES
367 : * 1) Assumes that the pack method of multiple grids are the same.
368 : *****************************************************************************
369 : */
370 2 : static int FindSectLen2to7 (char *c_ipack, sInt4 gribLen, sInt4 ns[8],
371 : char sectNum, sInt4 *curTot, sInt4 *nd2x3,
372 : short int *table50)
373 : {
374 : sInt4 sectLen; /* The length of the current section. */
375 : sInt4 li_temp; /* A temporary holder of sInt4s. */
376 :
377 2 : if ((sectNum == 2) || (sectNum == 3)) {
378 : /* Figure out the size of section 2 and 3. */
379 2 : if (*curTot + 5 > gribLen) {
380 0 : errSprintf ("ERROR: Ran out of data in Section 2 or 3\n");
381 0 : return -1;
382 : }
383 : /* Handle optional section 2. */
384 2 : if (c_ipack[*curTot + 4] == 2) {
385 0 : MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
386 0 : *curTot = *curTot + sectLen;
387 0 : if (ns[2] < sectLen)
388 0 : ns[2] = sectLen;
389 0 : if (*curTot + 5 > gribLen) {
390 0 : errSprintf ("ERROR: Ran out of data in Section 3\n");
391 0 : return -1;
392 : }
393 : }
394 : /* Handle section 3. */
395 2 : if (c_ipack[*curTot + 4] != 3) {
396 : errSprintf ("ERROR: Section 3 labeled as %d\n",
397 0 : c_ipack[*curTot + 4]);
398 0 : return -2;
399 : }
400 2 : MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
401 2 : if (ns[3] < sectLen)
402 2 : ns[3] = sectLen;
403 : /* While we are here, grab the total number of grid points nd2x3. */
404 2 : MEMCPY_BIG (&li_temp, c_ipack + *curTot + 6, 4);
405 2 : if (*nd2x3 < li_temp)
406 2 : *nd2x3 = li_temp;
407 2 : *curTot = *curTot + sectLen;
408 : }
409 : /*
410 : #ifdef DEBUG
411 : printf ("Section len (2=%ld) (3=%ld)\n", ns[2], ns[3]);
412 : #endif
413 : */
414 :
415 : /* Figure out the size of section 4. */
416 2 : if (*curTot + 5 > gribLen) {
417 0 : errSprintf ("ERROR: Ran out of data in Section 4\n");
418 0 : return -1;
419 : }
420 2 : if (c_ipack[*curTot + 4] != 4) {
421 0 : errSprintf ("ERROR: Section 4 labeled as %d\n", c_ipack[*curTot + 4]);
422 0 : return -2;
423 : }
424 2 : MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
425 2 : if (ns[4] < sectLen)
426 2 : ns[4] = sectLen;
427 2 : *curTot = *curTot + sectLen;
428 : /*
429 : #ifdef DEBUG
430 : printf ("Section len (4=%ld < %ld)\n", sectLen, ns[4]);
431 : #endif
432 : */
433 :
434 : /* Figure out the size of section 5. */
435 2 : if (*curTot + 5 > gribLen) {
436 0 : errSprintf ("ERROR: Ran out of data in Section 5\n");
437 0 : return -1;
438 : }
439 2 : if (c_ipack[*curTot + 4] != 5) {
440 0 : errSprintf ("ERROR: Section 5 labeled as %d\n", c_ipack[*curTot + 4]);
441 0 : return -2;
442 : }
443 2 : MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
444 : /* While we are here, grab the packing method. */
445 2 : MEMCPY_BIG (table50, c_ipack + *curTot + 9, 2);
446 2 : if (ns[5] < sectLen)
447 2 : ns[5] = sectLen;
448 2 : *curTot = *curTot + sectLen;
449 : /*
450 : #ifdef DEBUG
451 : printf ("Section len (5=%ld < %ld)\n", sectLen, ns[5]);
452 : #endif
453 : */
454 :
455 : /* Figure out the size of section 6. */
456 2 : if (*curTot + 5 > gribLen) {
457 0 : errSprintf ("ERROR: Ran out of data in Section 6\n");
458 0 : return -1;
459 : }
460 2 : if (c_ipack[*curTot + 4] != 6) {
461 0 : errSprintf ("ERROR: Section 6 labeled as %d\n", c_ipack[*curTot + 4]);
462 0 : return -2;
463 : }
464 2 : MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
465 2 : if (ns[6] < sectLen)
466 2 : ns[6] = sectLen;
467 2 : *curTot = *curTot + sectLen;
468 : /*
469 : #ifdef DEBUG
470 : printf ("Section len (6=%ld < %ld)\n", sectLen, ns[6]);
471 : #endif
472 : */
473 :
474 : /* Figure out the size of section 7. */
475 2 : if (*curTot + 5 > gribLen) {
476 0 : errSprintf ("ERROR: Ran out of data in Section 7\n");
477 0 : return -1;
478 : }
479 2 : if (c_ipack[*curTot + 4] != 7) {
480 0 : errSprintf ("ERROR: Section 7 labeled as %d\n", c_ipack[*curTot + 4]);
481 0 : return -2;
482 : }
483 2 : MEMCPY_BIG (§Len, c_ipack + *curTot, 4);
484 2 : if (ns[7] < sectLen)
485 2 : ns[7] = sectLen;
486 2 : *curTot = *curTot + sectLen;
487 : /*
488 : #ifdef DEBUG
489 : printf ("Section len (7=%ld < %ld)\n", sectLen, ns[7]);
490 : #endif
491 : */
492 2 : return 0;
493 : }
494 :
495 : /*****************************************************************************
496 : * FindSectLen() -- Review 12/2002
497 : *
498 : * Arthur Taylor / MDL
499 : *
500 : * PURPOSE
501 : * Looks through a GRIB message and finds out how big each section is.
502 : *
503 : * ARGUMENTS
504 : * c_ipack = The current GRIB2 message. (Input)
505 : * gribLen = Length of c_ipack. (Input)
506 : * ns = Array of section lengths. (Output)
507 : * nd2x3 = Total number of grid points (Output)
508 : * table50 = Type of packing used. (See code table 5.0) (GS5_SIMPLE = 0,
509 : * GS5_CMPLX = 2, GS5_CMPLXSEC = 3) (Output)
510 : *
511 : * FILES/DATABASES: None
512 : *
513 : * RETURNS: int (could use errSprintf())
514 : * 0 = OK
515 : * -1 = Ran of data in a section
516 : * -2 = Section not properly labeled.
517 : *
518 : * HISTORY
519 : * 9/2002 Arthur Taylor (MDL/RSIS): Created.
520 : * 11/2002 AAT: Updated.
521 : * 12/2002 (TK,AC,TB,&MS): Code Review.
522 : * 3/2003 AAT: Made it handle multiple grids in the same GRIB2 message.
523 : * 5/2003 AAT: Bug: Initialized size of section 2..6 to -1, instead
524 : * of 2..7.
525 : *
526 : * NOTES
527 : * 1) Assumes that the pack method of multiple grids are the same.
528 : *****************************************************************************
529 : */
530 2 : static int FindSectLen (char *c_ipack, sInt4 gribLen, sInt4 ns[8],
531 : sInt4 *nd2x3, short int *table50)
532 : {
533 : sInt4 curTot; /* Where we are in the current GRIB message. */
534 : char sectNum; /* Which section we are working with. */
535 : int ans; /* The return error code of FindSectLen2to7. */
536 : sInt4 sectLen; /* The length of the current section. */
537 : int i; /* counter as we init ns[]. */
538 :
539 2 : ns[0] = SECT0LEN_WORD * 4;
540 2 : curTot = ns[0];
541 :
542 : /* Figure out the size of section 1. */
543 2 : if (curTot + 5 > gribLen) {
544 0 : errSprintf ("ERROR: Ran out of data in Section 1\n");
545 0 : return -1;
546 : }
547 2 : if (c_ipack[curTot + 4] != 1) {
548 0 : errSprintf ("ERROR: Section 1 labeled as %d\n", c_ipack[curTot + 4]);
549 0 : return -2;
550 : }
551 2 : MEMCPY_BIG (&(ns[1]), c_ipack + curTot, 4);
552 2 : curTot += ns[1];
553 : /*
554 : #ifdef DEBUG
555 : printf ("Section len (0=%ld) (1=%ld)\n", ns[0], ns[1]);
556 : #endif
557 : */
558 2 : sectNum = 2;
559 14 : for (i = 2; i < 8; i++) {
560 12 : ns[i] = -1;
561 : }
562 2 : *nd2x3 = -1;
563 2 : do {
564 2 : if ((ans = FindSectLen2to7 (c_ipack, gribLen, ns, sectNum, &curTot,
565 : nd2x3, table50)) != 0) {
566 0 : return ans;
567 : }
568 : /* Try to read section 8. If it is "7777" == 926365495 regardless of
569 : * endian'ness then we have a simple message, otherwise it is complex,
570 : * and we need to read more. */
571 2 : memcpy (§Len, c_ipack + curTot, 4);
572 2 : if (sectLen == 926365495L) {
573 2 : sectNum = 8;
574 : } else {
575 0 : sectNum = c_ipack[curTot + 4];
576 0 : if ((sectNum < 2) || (sectNum > 7)) {
577 : errSprintf ("ERROR (FindSectLen): Couldn't find the end of the "
578 0 : "message\n");
579 0 : errSprintf ("and it doesn't appear to repeat sections.\n");
580 0 : errSprintf ("so it is probably an ASCII / binary bug\n");
581 : errSprintf ("Max Sect Lengths: %ld %ld %ld %ld %ld %ld %ld"
582 : " %ld\n", ns[0], ns[1], ns[2], ns[3], ns[4], ns[5],
583 0 : ns[6], ns[7]);
584 0 : return -2;
585 : }
586 : }
587 : } while (sectNum != 8);
588 2 : return 0;
589 : }
590 :
591 : /*****************************************************************************
592 : * IS_Init() -- Review 12/2002
593 : *
594 : * Arthur Taylor / MDL
595 : *
596 : * PURPOSE
597 : * Initialize the IS data structure. The IS_dataType is used to organize
598 : * and allocate all the arrays that the unpack library uses.
599 : * This makes an initial guess for the size of the arrays, and we use
600 : * realloc to increase the size if needed. The write up: "UNPK_GRIB2
601 : * 3/15/02" by Bryon Lawrence, Bob Glahn, David Rudack suggested these numbers
602 : *
603 : * ARGUMENTS
604 : * is = The data structure to initialize. (Output)
605 : *
606 : * FILES/DATABASES: None
607 : *
608 : * RETURNS: void
609 : *
610 : * HISTORY
611 : * 9/2002 Arthur Taylor (MDL/RSIS): Created.
612 : * 11/2002 AAT : Updated.
613 : * 12/2002 (TK,AC,TB,&MS): Code Review.
614 : *
615 : * NOTES
616 : * 1) Numbers not found in document were discused with Bob Glahn on 8/29/2002
617 : * 2) Possible exceptions:
618 : * template 3.120 could need ns[3] = 1600
619 : * template 4.30 could need a different ns4.
620 : * 3) sizeof(float) == sizeof(sInt4), and unpacker does not use both ain
621 : * and iain, so it is possible to have ain and iain point to the same
622 : * memory. Not sure if this is safe, so I haven't done it.
623 : *****************************************************************************
624 : */
625 6 : void IS_Init (IS_dataType *is)
626 : {
627 : int i; /* A simple loop counter. */
628 6 : is->ns[0] = 16;
629 6 : is->ns[1] = 21;
630 6 : is->ns[2] = 7;
631 6 : is->ns[3] = 96;
632 6 : is->ns[4] = 130; /* 60->130 in case there are some S4 time intervals */
633 6 : is->ns[5] = 49;
634 6 : is->ns[6] = 6;
635 6 : is->ns[7] = 8;
636 54 : for (i = 0; i < 8; i++) {
637 48 : is->is[i] = (sInt4 *) calloc (is->ns[i], sizeof (sInt4));
638 : }
639 : /* Allocate grid memory. */
640 6 : is->nd2x3 = 0;
641 6 : is->iain = NULL;
642 6 : is->ib = NULL;
643 : /* Allocate section 2 int memory. */
644 6 : is->nidat = 0;
645 6 : is->idat = NULL;
646 : /* Allocate section 2 float memory. */
647 6 : is->nrdat = 0;
648 6 : is->rdat = NULL;
649 : /* Allocate storage for ipack. */
650 6 : is->ipackLen = 0;
651 6 : is->ipack = NULL;
652 6 : }
653 :
654 : /*****************************************************************************
655 : * IS_Free() -- Review 12/2002
656 : *
657 : * Arthur Taylor / MDL
658 : *
659 : * PURPOSE
660 : * Free the memory allocated in the IS data structure.
661 : * The IS_dataType is used to organize and allocate all the arrays that the
662 : * unpack library uses.
663 : *
664 : * ARGUMENTS
665 : * is = The data structure to Free. (Output)
666 : *
667 : * FILES/DATABASES: None
668 : *
669 : * RETURNS: void
670 : *
671 : * HISTORY
672 : * 9/2002 Arthur Taylor (MDL/RSIS): Created.
673 : * 11/2002 AAT : Updated.
674 : * 12/2002 (TK,AC,TB,&MS): Code Review.
675 : *
676 : * NOTES
677 : *****************************************************************************
678 : */
679 6 : void IS_Free (IS_dataType *is)
680 : {
681 : int i; /* A simple loop counter. */
682 54 : for (i = 0; i < 8; i++) {
683 48 : free (is->is[i]);
684 48 : is->is[i] = NULL;
685 48 : is->ns[i] = 0;
686 : }
687 : /* Free grid memory. */
688 6 : free (is->iain);
689 6 : is->iain = NULL;
690 6 : free (is->ib);
691 6 : is->ib = NULL;
692 6 : is->nd2x3 = 0;
693 : /* Free section 2 int memory. */
694 6 : free (is->idat);
695 6 : is->idat = NULL;
696 6 : is->nidat = 0;
697 : /* Free section 2 float memory. */
698 6 : free (is->rdat);
699 6 : is->rdat = NULL;
700 6 : is->nrdat = 0;
701 : /* Free storage for ipack. */
702 6 : free (is->ipack);
703 6 : is->ipack = NULL;
704 6 : is->ipackLen = 0;
705 6 : }
706 :
707 : /*****************************************************************************
708 : * ReadGrib2Record() -- Review 12/2002
709 : *
710 : * Arthur Taylor / MDL
711 : *
712 : * PURPOSE
713 : * Reads a GRIB2 message from a file which is already opened and is pointing
714 : * at the correct message. It reads in the message storing the results in
715 : * Grib_Data which is of size grib_DataLen. If needed, it increases
716 : * grib_DataLen enough to fit the current message's grid. It converts (if
717 : * appropriate) the data in Grib_Data to the units specified in f_unit.
718 : *
719 : * In addition it updates offset, and stores the meta data returned by the
720 : * unpacker library in both IS, and (after parsing it) in meta.
721 : *
722 : * Note: It expects meta and IS to already be initialized through calls to
723 : * MetaInit(&meta) and IS_Init(&is) respectively.
724 : *
725 : * ARGUMENTS
726 : * fp = An opened GRIB2 file already at the correct message. (Input)
727 : * fileLen = Length of the opened file. (Input)
728 : * f_unit = 0 use GRIB2 units, 1 use English, 2 use metric. (Input)
729 : * Grib_Data = The read in GRIB2 grid. (Output)
730 : * grib_DataLen = Size of Grib_Data. (Output)
731 : * meta = A filled in meta structure (Output)
732 : * IS = The structure containing all the arrays that the
733 : * unpacker uses (Output)
734 : * subgNum = Which subgrid in the GRIB2 message is of interest.
735 : * (0 = first grid), if it can't find message subgNum,
736 : * returns -5, and an error message (Input)
737 : * majEarth = Used to override the GRIB major axis of earth. (Input)
738 : * minEarth = Used to override the GRIB minor axis of earth. (Input)
739 : * simpVer = The version of the simple weather code to use when parsing
740 : * the WxString. (Input)
741 : * f_endMsg = 1 means we finished reading the previous GRIB message, or
742 : * there was no previous GRIB message. 0 means that we need
743 : * to read a subgrid of the previous message. (Input/Output)
744 : * lwlt, uprt = If the lat is not -100, then lwlt, and uprt define a
745 : * subgrid that the user is interested in. Get the map
746 : * projection out of the GRIB2 message, and do everything
747 : * on the subgrid. (if lwlt, and uprt are not "correct", the
748 : * lat/lons may get swapped) (Input/Output)
749 : *
750 : * FILES/DATABASES:
751 : * An already opened "GRIB2" File
752 : *
753 : * RETURNS: int (could use errSprintf())
754 : * 0 = OK
755 : * -1 = Problems in section 0
756 : * -2 = Problems figuring out the Section Lengths.
757 : * -3 = Error returned by unpack library.
758 : * -4 = Problems parsing the Meta Data.
759 : *
760 : * HISTORY
761 : * 9/2002 Arthur Taylor (MDL/RSIS): Created.
762 : * 11/2002 AAT: Updated.
763 : * 12/2002 (TK,AC,TB,&MS): Code Review.
764 : * 1/2003 AAT: It wasn't error coded 208, but rather 202 to look for.
765 : * 3/2003 AAT: Modified handling of section 2 stuff (no loop)
766 : * 3/2003 AAT: Added ability to handle multiple grids in same message.
767 : * 4/2003 AAT: Added ability to call GRIB1 decoder for GRIB1 messages.
768 : * 5/2003 AAT: Update the offset for ReadGrib1.
769 : * 6/2003 Matthew T. Kallio (matt@wunderground.com):
770 : * "wmo" dimension increased to WMO_HEADER_LEN + 1 (for '\0' char)
771 : * 7/2003 AAT: switched to checking against element name for Wx instead
772 : * of pds2.sect2.ptrType == GS2_WXTYPE
773 : * 7/2003 AAT: Allowed user to override the radius of earth.
774 : * 8/2003 AAT: Removed dependence on fileLen and offset.
775 : * 2/2004 AAT: Added "f_endMsg" logic.
776 : * 2/2004 AAT: Added subgrid potential.
777 : * 2/2004 AAT: Added maj/min Earth override.
778 : *
779 : * NOTES
780 : * 1) Reason ns[7] is not MAX (IS.ns[], local_ns[]) is because local_ns[7]
781 : * is size of the packed message, but ns[7] refers to the returned meta
782 : * data which unpacker library found in section 7, which is a lot smaller.
783 : * 2) Problem: MDL's sect2 is packed and we have no idea how large it is
784 : * when unpacked. So we allocate room for 4000 sInt4s and 500 floats.
785 : * We then check 'jer' for error "202", if we find it we double the size
786 : * and call the unpacker again.
787 : * 3/26/2003: Changed this to be: try once with size
788 : * = max (32 * packed size, 4000)
789 : * Should be fewer calls (more memory intensive) same result, since we had
790 : * been doubling it 5 times.
791 : * 3) For Complex second order packing (GS5_CMPLXSEC) the unpacker needs nd5
792 : * (which is size of message) to be >= nd2x3 (size of grid).
793 : * 3a) Appears to also need this if simple packing, and has a bitmap.
794 : * 4) inew = 1: Currently we only expect 1 grid in 1 GRIB message, although
795 : * the library does allow for multiple grids in a GRIB message.
796 : * 5) iclean = 1: This only maters if there is bitmap data, otherwise it is
797 : * ignored. For bitmap data, if == 0, it embeds the given values for
798 : * xmissp, and xmisss. We don't embed because we don't know what to set
799 : * xmissp or xmisss to. Instead after we know the range, we choose a value
800 : * and walk through the bitmap setting grib_Data appropriately.
801 : * 5a) iclean = 0; This is because we do want the missing values embeded.
802 : * that is we want the missing values to be place holders.
803 : * 6) f_endMsg is true if in the past we either completed reading a message,
804 : * or we haven't read any messages. In either case we need to read the
805 : * next message from file. If f_endMsg is false, then there is more to read
806 : * from IS->ipack, so we don't want to throw it out, nor have to re-read
807 : * ipack from disk.
808 : *
809 : * Question: Should we double ns[2] when we double nrdat, and nidat?
810 : *****************************************************************************
811 : */
812 : #define SECT2_INIT_SIZE 4000
813 : #define UNPK_NUM_ERRORS 22
814 6 : int ReadGrib2Record (DataSource &fp, sChar f_unit, double **Grib_Data,
815 : uInt4 *grib_DataLen, grib_MetaData *meta,
816 : IS_dataType *IS, int subgNum, double majEarth,
817 : double minEarth, int simpVer, sInt4 *f_endMsg,
818 : LatLon *lwlf, LatLon *uprt)
819 : {
820 : sInt4 l3264b; /* Number of bits in a sInt4. Needed by FORTRAN
821 : * unpack library to determine if system has a 4
822 : * byte_ sInt4 or an 8 byte sInt4. */
823 : char *buff; /* Holds the info between records. */
824 : uInt4 buffLen; /* Length of info between records. */
825 : sInt4 sect0[SECT0LEN_WORD]; /* Holds the current Section 0. */
826 : uInt4 gribLen; /* Length of the current GRIB message. */
827 : sInt4 nd5; /* Size of grib message rounded up to the nearest
828 : * sInt4. */
829 : char *c_ipack; /* A char ptr to the message stored in IS->ipack */
830 : sInt4 local_ns[8]; /* Local copy of section lengths. */
831 : sInt4 nd2x3; /* Total number of grid points. */
832 : short int table50; /* Type of packing used. (See code table 5.0)
833 : * (GS5_SIMPLE==0, GS5_CMPLX==2, GS5_CMPLXSEC==3) */
834 : sInt4 nidat; /* Size of section 2 if it contains integer data. */
835 : sInt4 nrdat; /* Size of section 2 if it contains float data. */
836 : sInt4 inew; /* 1 if this is the first grid we are reading. 0 if
837 : * this is the second or later grid from the same
838 : * GRIB message. */
839 6 : sInt4 iclean = 0; /* 0 embed the missing values, 1 don't. */
840 : int j; /* Counter used to find the desired subgrid. */
841 6 : sInt4 kfildo = 5; /* FORTRAN Unit number for diagnostic info. Ignored,
842 : * unless library is compiled a particular way. */
843 : sInt4 ibitmap; /* 0 means no bitmap returned, otherwise 1. */
844 : float xmissp; /* The primary missing value. If iclean = 0, this
845 : * value is embeded in grid, otherwise it is the
846 : * value returned from the GRIB message. */
847 : float xmisss; /* The secondary missing value. If iclean = 0, this
848 : * value is embeded in grid, otherwise it is the
849 : * value returned from the GRIB message. */
850 : sInt4 jer[UNPK_NUM_ERRORS * 2]; /* Any Error codes along with their *
851 : * severity levels generated using the *
852 : * unpack GRIB2 library. */
853 6 : sInt4 ndjer = UNPK_NUM_ERRORS; /* The number of rows in JER( ). */
854 : sInt4 kjer; /* The actual number of errors returned in JER. */
855 : size_t i; /* counter as we loop through jer. */
856 : double unitM, unitB; /* values in y = m x + b used for unit conversion. */
857 : char unitName[15]; /* Holds the string name of the current unit. */
858 : int unitLen; /* String length of string name of current unit. */
859 : int version; /* Which version of GRIB is in this message. */
860 : sInt4 cnt; /* Used to help compact the weather table. */
861 : int x1, y1; /* The original grid coordinates of the lower left
862 : * corner of the subgrid. */
863 : int x2, y2; /* The original grid coordinates of the upper right
864 : * corner of the subgrid. */
865 : uChar f_subGrid; /* True if we have a subgrid. */
866 : sInt4 Nx, Ny; /* original size of the data. */
867 :
868 : /*
869 : * f_endMsg is 1 if in the past we either completed reading a message,
870 : * or we haven't read any messages. In either case we need to read the
871 : * next message from file.
872 : * If f_endMsg is false, then there is more to read from IS->ipack, so we
873 : * don't want to throw it out, nor have to re-read ipack from disk.
874 : */
875 6 : l3264b = sizeof (sInt4) * 8;
876 6 : buff = NULL;
877 6 : buffLen = 0;
878 6 : if (*f_endMsg == 1) {
879 6 : if (ReadSECT0 (fp, &buff, &buffLen, -1, sect0, &gribLen, &version) < 0) {
880 0 : preErrSprintf ("Inside ReadGrib2Record\n");
881 0 : free (buff);
882 0 : return -1;
883 : }
884 6 : meta->GribVersion = version;
885 6 : if (version == 1) {
886 4 : if (ReadGrib1Record (fp, f_unit, Grib_Data, grib_DataLen, meta, IS,
887 : sect0, gribLen, majEarth, minEarth) != 0) {
888 : preErrSprintf ("Problems with ReadGrib1Record called by "
889 0 : "ReadGrib2Record\n");
890 0 : free (buff);
891 0 : return -1;
892 : }
893 4 : *f_endMsg = 1;
894 4 : free (buff);
895 4 : return 0;
896 2 : } else if (version == -1) {
897 0 : if (ReadTDLPRecord (fp, Grib_Data, grib_DataLen, meta, IS,
898 : sect0, gribLen, majEarth, minEarth) != 0) {
899 : preErrSprintf ("Problems with ReadGrib1Record called by "
900 0 : "ReadGrib2Record\n");
901 0 : free (buff);
902 0 : return -1;
903 : }
904 0 : free (buff);
905 0 : return 0;
906 : }
907 :
908 : /*
909 : * Make room for entire message, and read it in.
910 : */
911 : /* nd5 needs to be gribLen in (sInt4) units rounded up. */
912 2 : nd5 = (gribLen + 3) / 4;
913 2 : if (nd5 > IS->ipackLen) {
914 2 : IS->ipackLen = nd5;
915 : IS->ipack = (sInt4 *) realloc ((void *) (IS->ipack),
916 2 : (IS->ipackLen) * sizeof (sInt4));
917 : }
918 2 : c_ipack = (char *) IS->ipack;
919 : /* Init last sInt4 to 0, to make sure that the padded bytes are 0. */
920 2 : IS->ipack[nd5 - 1] = 0;
921 : /* Init first 4 sInt4 to sect0. */
922 2 : memcpy (c_ipack, sect0, SECT0LEN_WORD * 4);
923 : /* Read in the rest of the message. */
924 2 : if (fp.DataSourceFread (c_ipack + SECT0LEN_WORD * 4, sizeof (char),
925 2 : (gribLen - SECT0LEN_WORD * 4)) != (gribLen - SECT0LEN_WORD * 4)) {
926 : errSprintf ("GribLen = %ld, SECT0Len_WORD = %d\n", gribLen,
927 0 : SECT0LEN_WORD);
928 0 : errSprintf ("Ran out of file\n");
929 0 : free (buff);
930 0 : return -1;
931 : }
932 :
933 : /*
934 : * Make sure the arrays are large enough for call to unpacker library.
935 : */
936 : /* FindSectLen Does not want (ipack / c_ipack) word swapped, because
937 : * that would make it much more confusing to find bytes in c_ipack. */
938 2 : if (FindSectLen (c_ipack, gribLen, local_ns, &nd2x3, &table50) < 0) {
939 0 : preErrSprintf ("Inside ReadGrib2Record.. Calling FindSectLen\n");
940 0 : free (buff);
941 0 : return -2;
942 : }
943 :
944 : /* Make sure all 'is' arrays except ns[7] are MAX (IS.ns[] ,
945 : * local_ns[]). See note 1 for reason to exclude ns[7] from MAX (). */
946 16 : for (i = 0; i < 7; i++) {
947 14 : if (local_ns[i] > IS->ns[i]) {
948 0 : IS->ns[i] = local_ns[i];
949 0 : IS->is[i] = (sInt4 *) realloc ((void *) (IS->is[i]),
950 0 : IS->ns[i] * sizeof (sInt4));
951 : }
952 : }
953 :
954 : /* Allocate room for sect 2. If local_ns[2] = -1 there is no sect 2. */
955 2 : if (local_ns[2] == -1) {
956 2 : nidat = 10;
957 2 : nrdat = 10;
958 : } else {
959 : /*
960 : * See note 2) We have a section 2, so use:
961 : * MAX (32 * local_ns[2],SECT2_INTSIZE)
962 : * and MAX (32 * local_ns[2],SECT2_FLOATSIZE)
963 : * for size of section 2 unpacked.
964 : */
965 0 : nidat = (32 * local_ns[2] < SECT2_INIT_SIZE) ? SECT2_INIT_SIZE :
966 0 : 32 * local_ns[2];
967 0 : nrdat = nidat;
968 : }
969 2 : if (nidat > IS->nidat) {
970 2 : IS->nidat = nidat;
971 : IS->idat = (sInt4 *) realloc ((void *) IS->idat,
972 2 : IS->nidat * sizeof (sInt4));
973 : }
974 2 : if (nrdat > IS->nrdat) {
975 2 : IS->nrdat = nrdat;
976 : IS->rdat = (float *) realloc ((void *) IS->rdat,
977 2 : IS->nrdat * sizeof (float));
978 : }
979 : /* Make sure we have room for the GRID part of the output. */
980 2 : if (nd2x3 > IS->nd2x3) {
981 2 : IS->nd2x3 = nd2x3;
982 : IS->iain = (sInt4 *) realloc ((void *) IS->iain,
983 2 : IS->nd2x3 * sizeof (sInt4));
984 : IS->ib = (sInt4 *) realloc ((void *) IS->ib,
985 2 : IS->nd2x3 * sizeof (sInt4));
986 : }
987 : /* See note 3) If table50 == 3, unpacker library needs nd5 >= nd2x3. */
988 2 : if ((table50 == 3) || (table50 == 0)) {
989 2 : if (nd5 < nd2x3) {
990 2 : nd5 = nd2x3;
991 2 : if (nd5 > IS->ipackLen) {
992 2 : IS->ipackLen = nd5;
993 : IS->ipack = (sInt4 *) realloc ((void *) (IS->ipack),
994 2 : IS->ipackLen * sizeof (sInt4));
995 : }
996 : /* Don't need to do the following, but we do in case code
997 : * changes. */
998 2 : c_ipack = (char *) IS->ipack;
999 : }
1000 : }
1001 2 : IS->nd5 = nd5;
1002 : /* Unpacker library requires ipack to be MSB. */
1003 : /*
1004 : #ifdef DEBUG
1005 : if (1==1) {
1006 : FILE *fp = fopen ("test.bin", "wb");
1007 : fwrite (IS->ipack, sizeof (sInt4), IS->nd5, fp);
1008 : fclose (fp);
1009 : }
1010 : #endif
1011 : */
1012 : #ifdef LITTLE_ENDIAN
1013 2 : memswp (IS->ipack, sizeof (sInt4), IS->nd5);
1014 : #endif
1015 : } else {
1016 0 : gribLen = IS->ipack[3];
1017 : }
1018 2 : free (buff);
1019 :
1020 : /* Loop through the grib message looking for the subgNum grid. subgNum
1021 : * goes from 0 to n-1. */
1022 4 : for (j = 0; j <= subgNum; j++) {
1023 2 : if (j == 0) {
1024 2 : inew = 1;
1025 : } else {
1026 0 : inew = 0;
1027 : }
1028 :
1029 : /* Note we are getting data back either as a float or an int, but not
1030 : * both, so we don't need to allocated room for both. */
1031 : unpk_grib2 (&kfildo, (float *) (IS->iain), IS->iain, &(IS->nd2x3),
1032 : IS->idat, &(IS->nidat), IS->rdat, &(IS->nrdat), IS->is[0],
1033 : &(IS->ns[0]), IS->is[1], &(IS->ns[1]), IS->is[2],
1034 : &(IS->ns[2]), IS->is[3], &(IS->ns[3]), IS->is[4],
1035 : &(IS->ns[4]), IS->is[5], &(IS->ns[5]), IS->is[6],
1036 : &(IS->ns[6]), IS->is[7], &(IS->ns[7]), IS->ib, &ibitmap,
1037 : IS->ipack, &(IS->nd5), &xmissp, &xmisss, &inew, &iclean,
1038 2 : &l3264b, f_endMsg, jer, &ndjer, &kjer);
1039 : /*
1040 : * Check for error messages...
1041 : * If we get an error message, print it, and return.
1042 : */
1043 18 : for (i = 0; i < (uInt4) kjer; i++) {
1044 16 : if (jer[ndjer + i] == 0) {
1045 : /* no error. */
1046 0 : } else if (jer[ndjer + i] == 1) {
1047 : /* Warning. */
1048 : #ifdef DEBUG
1049 : printf ("Warning: Unpack library warning code (%d %d)\n",
1050 : jer[i], jer[ndjer + i]);
1051 : #endif
1052 : } else {
1053 : /* BAD Error. */
1054 : errSprintf ("ERROR: Unpack library error code (%ld %ld)\n",
1055 0 : jer[i], jer[ndjer + i]);
1056 0 : return -3;
1057 : }
1058 : }
1059 : }
1060 :
1061 : /* Parse the meta data out. */
1062 2 : if (MetaParse (meta, IS->is[0], IS->ns[0], IS->is[1], IS->ns[1],
1063 : IS->is[2], IS->ns[2], IS->rdat, IS->nrdat, IS->idat,
1064 : IS->nidat, IS->is[3], IS->ns[3], IS->is[4], IS->ns[4],
1065 : IS->is[5], IS->ns[5], gribLen, xmissp, xmisss, simpVer)
1066 : != 0) {
1067 : #ifdef DEBUG
1068 : FILE *fp;
1069 : if ((fp = fopen ("dump.is0", "wt")) != NULL) {
1070 : for (i = 0; i < 8; i++) {
1071 : fprintf (fp, "---Section %d---\n", (int) i);
1072 : for (j = 1; j <= IS->ns[i]; j++) {
1073 : fprintf (fp, "IS%d Item %d = %d\n", (int) i, (int) j, IS->is[i][j - 1]);
1074 : }
1075 : }
1076 : fclose (fp);
1077 : }
1078 : #endif
1079 0 : preErrSprintf ("Inside ReadGrib2Record.. Problems in MetaParse\n");
1080 0 : return -4;
1081 : }
1082 :
1083 2 : if ((majEarth > 6000) && (majEarth < 7000)) {
1084 0 : if ((minEarth > 6000) && (minEarth < 7000)) {
1085 0 : meta->gds.f_sphere = 0;
1086 0 : meta->gds.majEarth = majEarth;
1087 0 : meta->gds.minEarth = minEarth;
1088 : } else {
1089 0 : meta->gds.f_sphere = 1;
1090 0 : meta->gds.majEarth = majEarth;
1091 0 : meta->gds.minEarth = majEarth;
1092 : }
1093 : }
1094 :
1095 : /* Figure out an equation to pass to ParseGrid to convert the units for
1096 : * this grid. */
1097 : /*
1098 : if (ComputeUnit (meta->pds2.prodType, meta->pds2.sect4.templat,
1099 : meta->pds2.sect4.cat, meta->pds2.sect4.subcat, f_unit,
1100 : &unitM, &unitB, unitName) == 0) {
1101 : */
1102 2 : if (ComputeUnit (meta->convert, meta->unitName, f_unit, &unitM, &unitB,
1103 : unitName) == 0) {
1104 2 : unitLen = strlen (unitName);
1105 : meta->unitName = (char *) realloc ((void *) (meta->unitName),
1106 2 : 1 + unitLen * sizeof (char));
1107 2 : strncpy (meta->unitName, unitName, unitLen);
1108 2 : meta->unitName[unitLen] = '\0';
1109 : }
1110 :
1111 : /* compute the subgrid. */
1112 : /*
1113 : if ((lwlf->lat != -100) && (uprt->lat != -100)) {
1114 : Nx = meta->gds.Nx;
1115 : Ny = meta->gds.Ny;
1116 : if (computeSubGrid (lwlf, &x1, &y1, uprt, &x2, &y2, &(meta->gds),
1117 : &newGds) != 0) {
1118 : preErrSprintf ("ERROR: In compute subgrid.\n");
1119 : return 1;
1120 : }
1121 : // I couldn't decide if I should "permanently" change the GDS or not.
1122 : // when I wrote computeSubGrid. If next line stays, really should
1123 : // rewrite computeSubGrid.
1124 : memcpy (&(meta->gds), &newGds, sizeof (gdsType));
1125 : f_subGrid = 1;
1126 : } else {
1127 : Nx = meta->gds.Nx;
1128 : Ny = meta->gds.Ny;
1129 : x1 = 1;
1130 : x2 = Nx;
1131 : y1 = 1;
1132 : y2 = Ny;
1133 : f_subGrid = 0;
1134 : }
1135 : */
1136 2 : Nx = meta->gds.Nx;
1137 2 : Ny = meta->gds.Ny;
1138 2 : x1 = 1;
1139 2 : x2 = Nx;
1140 2 : y1 = 1;
1141 2 : y2 = Ny;
1142 2 : f_subGrid = 0;
1143 :
1144 : /* Figure out if we need iain or ain, and set it to Grib_Data. At the
1145 : * same time handle any bitmaps, and compute some statistics. */
1146 2 : if ((f_subGrid) && (meta->gds.scan != 64)) {
1147 0 : errSprintf ("Can not do a subgrid of non scanmode 64 grid yet.\n");
1148 0 : return -3;
1149 : }
1150 :
1151 2 : if (strcmp (meta->element, "Wx") != 0) {
1152 : ParseGrid (&(meta->gridAttrib), Grib_Data, grib_DataLen, Nx, Ny,
1153 : meta->gds.scan, IS->iain, ibitmap, IS->ib, unitM, unitB, 0,
1154 2 : NULL, f_subGrid, x1, y1, x2, y2);
1155 : } else {
1156 : /* Handle weather grid. ParseGrid looks up the values... If they are
1157 : * "<Invalid>" it sets it to missing (or creates one). If the table
1158 : * entry is used it sets f_valid to 2. */
1159 : ParseGrid (&(meta->gridAttrib), Grib_Data, grib_DataLen, Nx, Ny,
1160 : meta->gds.scan, IS->iain, ibitmap, IS->ib, unitM, unitB, 1,
1161 : (sect2_WxType *) &(meta->pds2.sect2.wx), f_subGrid, x1, y1,
1162 0 : x2, y2);
1163 :
1164 : /* compact the table to only those which are actually used. */
1165 0 : cnt = 0;
1166 0 : for (i = 0; i < meta->pds2.sect2.wx.dataLen; i++) {
1167 0 : if (meta->pds2.sect2.wx.ugly[i].f_valid == 2) {
1168 0 : meta->pds2.sect2.wx.ugly[i].validIndex = cnt;
1169 0 : cnt++;
1170 0 : } else if (meta->pds2.sect2.wx.ugly[i].f_valid == 3) {
1171 0 : meta->pds2.sect2.wx.ugly[i].f_valid = 0;
1172 0 : meta->pds2.sect2.wx.ugly[i].validIndex = cnt;
1173 0 : cnt++;
1174 : } else {
1175 0 : meta->pds2.sect2.wx.ugly[i].validIndex = -1;
1176 : }
1177 : }
1178 : }
1179 :
1180 : /* Figure out some other non-section oriented meta data. */
1181 : /* strftime (meta->refTime, 20, "%Y%m%d%H%M",
1182 : gmtime (&(meta->pds2.refTime)));
1183 : */
1184 2 : Clock_Print (meta->refTime, 20, meta->pds2.refTime, "%Y%m%d%H%M", 0);
1185 : /*
1186 : strftime (meta->validTime, 20, "%Y%m%d%H%M",
1187 : gmtime (&(meta->pds2.sect4.validTime)));
1188 : */
1189 : Clock_Print (meta->validTime, 20, meta->pds2.sect4.validTime,
1190 2 : "%Y%m%d%H%M", 0);
1191 :
1192 2 : meta->deltTime = (sInt4) (meta->pds2.sect4.validTime - meta->pds2.refTime);
1193 :
1194 2 : return 0;
1195 : }
|
