1 : /*****************************************************************************
2 : * degrib1.c
3 : *
4 : * DESCRIPTION
5 : * This file contains the main driver routines to unpack GRIB 1 files.
6 : *
7 : * HISTORY
8 : * 4/2003 Arthur Taylor (MDL / RSIS): Created.
9 : *
10 : * NOTES
11 : * GRIB 1 files are assumed to be big endian.
12 : *****************************************************************************
13 : */
14 :
15 : #include <stdio.h>
16 : #include <string.h>
17 : #include <stdlib.h>
18 : #include <math.h>
19 : #include "degrib2.h"
20 : #include "myerror.h"
21 : #include "myassert.h"
22 : #include "memendian.h"
23 : #include "scan.h"
24 : #include "degrib1.h"
25 : #include "metaname.h"
26 : #include "clock.h"
27 : #include "cpl_error.h"
28 :
29 : /* default missing data value (see: bitmap GRIB1: sect3 and sect4) */
30 : /* UNDEFINED is default, UNDEFINED_PRIM is desired choice. */
31 : #define UNDEFINED 9.999e20
32 : #define UNDEFINED_PRIM 9999
33 :
34 : #define GRIB_UNSIGN_INT3(a,b,c) ((a<<16)+(b<<8)+c)
35 : #define GRIB_UNSIGN_INT2(a,b) ((a<<8)+b)
36 : #define GRIB_SIGN_INT3(a,b,c) ((1-(int) ((unsigned) (a & 0x80) >> 6)) * (int) (((a & 127) << 16)+(b<<8)+c))
37 : #define GRIB_SIGN_INT2(a,b) ((1-(int) ((unsigned) (a & 0x80) >> 6)) * (int) (((a & 0x7f) << 8) + b))
38 :
39 : /* various centers */
40 : #define NMC 7
41 : #define US_OTHER 9
42 : #define CPTEC 46
43 : /* Canada Center */
44 : #define CMC 54
45 : #define AFWA 57
46 : #define DWD 78
47 : #define ECMWF 98
48 : #define ATHENS 96
49 :
50 : /* various subcenters */
51 : #define SUBCENTER_MDL 14
52 : #define SUBCENTER_TDL 11
53 :
54 : /* The idea of rean or opn is to give a warning about default choice of
55 : which table to use. */
56 : #define DEF_NCEP_TABLE rean_nowarn
57 : enum Def_NCEP_Table { rean, opn, rean_nowarn, opn_nowarn };
58 :
59 : extern GRIB1ParmTable parm_table_ncep_opn[256];
60 : extern GRIB1ParmTable parm_table_ncep_reanal[256];
61 : extern GRIB1ParmTable parm_table_ncep_tdl[256];
62 : extern GRIB1ParmTable parm_table_ncep_mdl[256];
63 : extern GRIB1ParmTable parm_table_omb[256];
64 : extern GRIB1ParmTable parm_table_nceptab_129[256];
65 : extern GRIB1ParmTable parm_table_nceptab_130[256];
66 : extern GRIB1ParmTable parm_table_nceptab_131[256];
67 :
68 : extern GRIB1ParmTable parm_table_nohrsc[256];
69 :
70 : extern GRIB1ParmTable parm_table_cptec_254[256];
71 :
72 : extern GRIB1ParmTable parm_table_afwa_000[256];
73 : extern GRIB1ParmTable parm_table_afwa_001[256];
74 : extern GRIB1ParmTable parm_table_afwa_002[256];
75 : extern GRIB1ParmTable parm_table_afwa_003[256];
76 : extern GRIB1ParmTable parm_table_afwa_010[256];
77 : extern GRIB1ParmTable parm_table_afwa_011[256];
78 :
79 : extern GRIB1ParmTable parm_table_dwd_002[256];
80 : extern GRIB1ParmTable parm_table_dwd_201[256];
81 : extern GRIB1ParmTable parm_table_dwd_202[256];
82 : extern GRIB1ParmTable parm_table_dwd_203[256];
83 :
84 : extern GRIB1ParmTable parm_table_ecmwf_128[256];
85 : extern GRIB1ParmTable parm_table_ecmwf_129[256];
86 : extern GRIB1ParmTable parm_table_ecmwf_130[256];
87 : extern GRIB1ParmTable parm_table_ecmwf_131[256];
88 : extern GRIB1ParmTable parm_table_ecmwf_140[256];
89 : extern GRIB1ParmTable parm_table_ecmwf_150[256];
90 : extern GRIB1ParmTable parm_table_ecmwf_160[256];
91 : extern GRIB1ParmTable parm_table_ecmwf_170[256];
92 : extern GRIB1ParmTable parm_table_ecmwf_180[256];
93 :
94 : extern GRIB1ParmTable parm_table_athens[256];
95 :
96 : extern GRIB1ParmTable parm_table_cmc[256];
97 :
98 : extern GRIB1ParmTable parm_table_undefined[256];
99 :
100 : /*****************************************************************************
101 : * Choose_ParmTable() --
102 : *
103 : * Arthur Taylor / MDL
104 : *
105 : * PURPOSE
106 : * Chooses the correct Parameter table depending on what is in the GRIB1
107 : * message's "Product Definition Section".
108 : *
109 : * ARGUMENTS
110 : * pdsMeta = The filled out pdsMeta data structure to base choice on. (Input)
111 : * center = The Center that created the data (Input)
112 : * subcenter = The Sub Center that created the data (Input)
113 : *
114 : * FILES/DATABASES: None
115 : *
116 : * RETURNS: ParmTable (appropriate parameter table.)
117 : *
118 : * HISTORY
119 : * <unknown> : wgrib library : cnames.c
120 : * 4/2003 Arthur Taylor (MDL/RSIS): Modified
121 : * 10/2005 AAT: Adjusted to take center, subcenter
122 : *
123 : * NOTES
124 : *****************************************************************************
125 : */
126 20 : static GRIB1ParmTable *Choose_ParmTable (pdsG1Type *pdsMeta,
127 : unsigned short int center,
128 : unsigned short int subcenter)
129 : {
130 : int process; /* The process ID from the GRIB1 message. */
131 :
132 20 : switch (center) {
133 : case NMC:
134 20 : if (pdsMeta->mstrVersion <= 3) {
135 16 : switch (subcenter) {
136 : case 1:
137 0 : return &parm_table_ncep_reanal[0];
138 : case SUBCENTER_TDL:
139 0 : return &parm_table_ncep_tdl[0];
140 : case SUBCENTER_MDL:
141 0 : return &parm_table_ncep_mdl[0];
142 : }
143 : }
144 : /* figure out if NCEP opn or reanalysis */
145 20 : switch (pdsMeta->mstrVersion) {
146 : case 0:
147 10 : return &parm_table_ncep_opn[0];
148 : case 1:
149 : case 2:
150 6 : process = pdsMeta->genProcess;
151 6 : if ((subcenter != 0) || ((process != 80) && (process != 180))) {
152 6 : return &parm_table_ncep_opn[0];
153 : }
154 : /* At this point could be either the opn or reanalysis table */
155 : switch (DEF_NCEP_TABLE) {
156 : case opn_nowarn:
157 : return &parm_table_ncep_opn[0];
158 : case rean_nowarn:
159 0 : return &parm_table_ncep_reanal[0];
160 : }
161 : break;
162 : case 3:
163 0 : return &parm_table_ncep_opn[0];
164 : case 128:
165 0 : return &parm_table_omb[0];
166 : case 129:
167 4 : return &parm_table_nceptab_129[0];
168 : case 130:
169 0 : return &parm_table_nceptab_130[0];
170 : case 131:
171 0 : return &parm_table_nceptab_131[0];
172 : }
173 0 : break;
174 : case AFWA:
175 0 : switch (subcenter) {
176 : case 0:
177 0 : return &parm_table_afwa_000[0];
178 : case 1:
179 : case 4:
180 0 : return &parm_table_afwa_001[0];
181 : case 2:
182 0 : return &parm_table_afwa_002[0];
183 : case 3:
184 0 : return &parm_table_afwa_003[0];
185 : case 10:
186 0 : return &parm_table_afwa_010[0];
187 : case 11:
188 0 : return &parm_table_afwa_011[0];
189 : /* case 5:*/
190 : /* Didn't have a table 5. */
191 : }
192 0 : break;
193 : case ECMWF:
194 0 : switch (pdsMeta->mstrVersion) {
195 : case 128:
196 0 : return &parm_table_ecmwf_128[0];
197 : case 129:
198 0 : return &parm_table_ecmwf_129[0];
199 : case 130:
200 0 : return &parm_table_ecmwf_130[0];
201 : case 131:
202 0 : return &parm_table_ecmwf_131[0];
203 : case 140:
204 0 : return &parm_table_ecmwf_140[0];
205 : case 150:
206 0 : return &parm_table_ecmwf_150[0];
207 : case 160:
208 0 : return &parm_table_ecmwf_160[0];
209 : case 170:
210 0 : return &parm_table_ecmwf_170[0];
211 : case 180:
212 0 : return &parm_table_ecmwf_180[0];
213 : }
214 0 : break;
215 : case DWD:
216 0 : switch (pdsMeta->mstrVersion) {
217 : case 2:
218 0 : return &parm_table_dwd_002[0];
219 : case 201:
220 0 : return &parm_table_dwd_201[0];
221 : case 202:
222 0 : return &parm_table_dwd_202[0];
223 : case 203:
224 0 : return &parm_table_dwd_203[0];
225 : }
226 0 : break;
227 : case CPTEC:
228 0 : switch (pdsMeta->mstrVersion) {
229 : case 254:
230 0 : return &parm_table_cptec_254[0];
231 : }
232 0 : break;
233 : case US_OTHER:
234 0 : switch (subcenter) {
235 : case 163:
236 0 : return &parm_table_nohrsc[0];
237 : }
238 0 : break;
239 : case ATHENS:
240 0 : return &parm_table_athens[0];
241 : break;
242 : case CMC:
243 0 : return &parm_table_cmc[0];
244 : break;
245 : }
246 0 : if ((pdsMeta->mstrVersion > 3) || (pdsMeta->cat > 127)) {
247 : CPLDebug (
248 : "GRIB",
249 : "Undefined parameter table (center %d-%d table %d).",
250 0 : center, subcenter, pdsMeta->mstrVersion);
251 : }
252 0 : return &parm_table_undefined[0];
253 : }
254 :
255 : /*****************************************************************************
256 : * GRIB1_Table2LookUp() --
257 : *
258 : * Arthur Taylor / MDL
259 : *
260 : * PURPOSE
261 : * Returns the variable name (type of data) and comment (longer form of the
262 : * name) for the data that is in the GRIB1 message.
263 : *
264 : * ARGUMENTS
265 : * name = A pointer to the resulting short name. (Output)
266 : * comment = A pointer to the resulting long name. (Output)
267 : * pdsMeta = The filled out pdsMeta data structure to base choice on. (Input)
268 : * center = The Center that created the data (Input)
269 : * subcenter = The Sub Center that created the data (Input)
270 : *
271 : * FILES/DATABASES: None
272 : *
273 : * RETURNS: void
274 : *
275 : * HISTORY
276 : * 4/2003 Arthur Taylor (MDL/RSIS): Created
277 : * 10/2005 AAT: Adjusted to take center, subcenter
278 : *
279 : * NOTES
280 : *****************************************************************************
281 : */
282 20 : static void GRIB1_Table2LookUp (pdsG1Type *pdsMeta, const char **name,
283 : const char **comment, const char **unit,
284 : int *convert,
285 : unsigned short int center,
286 : unsigned short int subcenter)
287 : {
288 : GRIB1ParmTable *table; /* The paramter table choosen by the pdsMeta data */
289 :
290 20 : table = Choose_ParmTable (pdsMeta, center, subcenter);
291 20 : if ((center == NMC) && (pdsMeta->mstrVersion == 129)
292 : && (pdsMeta->cat == 180)) {
293 0 : if (pdsMeta->timeRange == 3) {
294 0 : *name = "AVGOZCON";
295 0 : *comment = "Average Ozone Concentration";
296 0 : *unit = "PPB";
297 0 : *convert = UC_NONE;
298 0 : return;
299 : }
300 : }
301 20 : *name = table[pdsMeta->cat].name;
302 20 : *comment = table[pdsMeta->cat].comment;
303 20 : *unit = table[pdsMeta->cat].unit;
304 20 : *convert = table[pdsMeta->cat].convert;
305 : /* printf ("%s %s %s\n", *name, *comment, *unit);*/
306 : }
307 :
308 : extern GRIB1SurfTable GRIB1Surface[256];
309 :
310 : /* Similar to metaname.c :: ParseLevelName() */
311 20 : static void GRIB1_Table3LookUp (pdsG1Type *pdsMeta, char **shortLevelName,
312 : char **longLevelName)
313 : {
314 20 : uChar type = pdsMeta->levelType;
315 : uChar level1, level2;
316 :
317 20 : free (*shortLevelName);
318 20 : *shortLevelName = NULL;
319 20 : free (*longLevelName);
320 20 : *longLevelName = NULL;
321 : /* Find out if val is a 2 part value or not */
322 20 : if (GRIB1Surface[type].f_twoPart) {
323 0 : level1 = (pdsMeta->levelVal >> 8);
324 0 : level2 = (pdsMeta->levelVal & 0xff);
325 : reallocSprintf (shortLevelName, "%d-%d-%s", level1, level2,
326 0 : GRIB1Surface[type].name);
327 : reallocSprintf (longLevelName, "%d-%d[%s] %s (%s)", level1, level2,
328 : GRIB1Surface[type].unit, GRIB1Surface[type].name,
329 0 : GRIB1Surface[type].comment);
330 : } else {
331 : reallocSprintf (shortLevelName, "%d-%s", pdsMeta->levelVal,
332 20 : GRIB1Surface[type].name);
333 : reallocSprintf (longLevelName, "%d[%s] %s (%s)", pdsMeta->levelVal,
334 : GRIB1Surface[type].unit, GRIB1Surface[type].name,
335 20 : GRIB1Surface[type].comment);
336 : }
337 20 : }
338 :
339 : /*****************************************************************************
340 : * fval_360() --
341 : *
342 : * Albion Taylor / ARL
343 : *
344 : * PURPOSE
345 : * Converts an IBM360 floating point number to an IEEE floating point
346 : * number. The IBM floating point spec represents the fraction as the last
347 : * 3 bytes of the number, with 0xffffff being just shy of 1.0. The first byte
348 : * leads with a sign bit, and the last seven bits represent the powers of 16
349 : * (not 2), with a bias of 0x40 giving 16^0.
350 : *
351 : * ARGUMENTS
352 : * aval = A sInt4 containing the original IBM 360 number. (Input)
353 : *
354 : * FILES/DATABASES: None
355 : *
356 : * RETURNS: double = the value that aval represents.
357 : *
358 : * HISTORY
359 : * <unknown> Albion Taylor (ARL): Created
360 : * 4/2003 Arthur Taylor (MDL/RSIS): Cleaned up.
361 : * 5/2003 AAT: some kind of Bug due to optimizations...
362 : * -1055916032 => 0 instead of -1
363 : *
364 : * NOTES
365 : *****************************************************************************
366 : */
367 4 : static double fval_360 (uInt4 aval)
368 : {
369 : double pow16;
370 : /* short int *ptr = (short int *) (&pow16);*/
371 4 : void *voidPtr = &pow16;
372 4 : short int *ptr = (short int *) voidPtr;
373 :
374 : #ifdef LITTLE_ENDIAN
375 4 : ptr[3] = ((((aval >> 24) & 0x7f) << 2) + (0x3ff - 0x100)) << 4;
376 4 : ptr[2] = 0;
377 4 : ptr[1] = 0;
378 4 : ptr[0] = 0;
379 : #else
380 : ptr[0] = ((((aval >> 24) & 0x7f) << 2) + (0x3ff - 0x100)) << 4;
381 : ptr[1] = 0;
382 : ptr[2] = 0;
383 : ptr[3] = 0;
384 : #endif
385 : return ((aval & 0x80000000) ? -pow16 : pow16) *
386 4 : (aval & 0xffffff) / ((double) 0x1000000);
387 : }
388 :
389 : /*****************************************************************************
390 : * ReadGrib1Sect1() --
391 : *
392 : * Arthur Taylor / MDL
393 : *
394 : * PURPOSE
395 : * Parses the GRIB1 "Product Definition Section" or section 1, filling out
396 : * the pdsMeta data structure.
397 : *
398 : * ARGUMENTS
399 : * pds = The compressed part of the message dealing with "PDS". (Input)
400 : * gribLen = The total length of the GRIB1 message. (Input)
401 : * curLoc = Current location in the GRIB1 message. (Output)
402 : * pdsMeta = The filled out pdsMeta data structure. (Output)
403 : * f_gds = boolean if there is a Grid Definition Section. (Output)
404 : * gridID = The Grid ID. (Output)
405 : * f_bms = boolean if there is a Bitmap Section. (Output)
406 : * DSF = Decimal Scale Factor for unpacking the data. (Output)
407 : * center = The Center that created the data (Output)
408 : * subcenter = The Sub Center that created the data (Output)
409 : *
410 : * FILES/DATABASES: None
411 : *
412 : * RETURNS: int (could use errSprintf())
413 : * 0 = OK
414 : * -1 = gribLen is too small.
415 : *
416 : * HISTORY
417 : * 4/2003 Arthur Taylor (MDL/RSIS): Created.
418 : * 5/2004 AAT: Paid attention to table 5 (Time range indicator) and which of
419 : * P1 and P2 are the valid times.
420 : * 10/2005 AAT: Adjusted to take center, subcenter
421 : *
422 : * NOTES
423 : *****************************************************************************
424 : */
425 20 : static int ReadGrib1Sect1 (uChar *pds, uInt4 gribLen, uInt4 *curLoc,
426 : pdsG1Type *pdsMeta, char *f_gds, uChar *gridID,
427 : char *f_bms, short int *DSF,
428 : unsigned short int *center,
429 : unsigned short int *subcenter)
430 : {
431 : sInt4 sectLen; /* Length in bytes of the current section. */
432 : int year; /* The year of the GRIB1 Message. */
433 : double P1_DeltaTime; /* Used to parse the time for P1 */
434 : double P2_DeltaTime; /* Used to parse the time for P2 */
435 : uInt4 uli_temp;
436 : #ifdef DEBUG
437 : /*
438 : int i;
439 : */
440 : #endif
441 :
442 20 : sectLen = GRIB_UNSIGN_INT3 (*pds, pds[1], pds[2]);
443 : #ifdef DEBUG
444 : /*
445 : printf ("Section 1 length = %ld\n", sectLen);
446 : for (i = 0; i < sectLen; i++) {
447 : printf ("Sect1: item %d = %d\n", i + 1, pds[i]);
448 : }
449 : printf ("Century is item 25\n");
450 : */
451 : #endif
452 20 : *curLoc += sectLen;
453 20 : if (*curLoc > gribLen) {
454 0 : errSprintf ("Ran out of data in PDS (GRIB 1 Section 1)\n");
455 0 : return -1;
456 : }
457 20 : pds += 3;
458 20 : pdsMeta->mstrVersion = *(pds++);
459 20 : *center = *(pds++);
460 20 : pdsMeta->genProcess = *(pds++);
461 20 : *gridID = *(pds++);
462 20 : *f_gds = GRIB2BIT_1 & *pds;
463 20 : *f_bms = GRIB2BIT_2 & *pds;
464 20 : pds++;
465 20 : pdsMeta->cat = *(pds++);
466 20 : pdsMeta->levelType = *(pds++);
467 20 : pdsMeta->levelVal = GRIB_UNSIGN_INT2 (*pds, pds[1]);
468 20 : pds += 2;
469 20 : if (*pds == 0) {
470 : /* The 12 is because we have increased pds by 12. (but 25 is in
471 : * reference of 1..25, so we need another -1) */
472 0 : year = (pds[25 - 13] * 100);
473 : } else {
474 : /* The 12 is because we have increased pds by 12. (but 25 is in
475 : * reference of 1..25, so we need another -1) */
476 20 : year = *pds + ((pds[25 - 13] - 1) * 100);
477 :
478 : /* It seems like some old files (such as spring/I000176.grb)
479 : do not have a century byte, and assum 19xx. */
480 : // if( (year < 1900 || year > 2100) && *pds >= 0 && *pds < 100 )
481 : // year = *pds + 1900;
482 : }
483 :
484 20 : if (ParseTime (&(pdsMeta->refTime), year, pds[1], pds[2], pds[3], pds[4],
485 : 0) != 0) {
486 0 : preErrSprintf ("Error In call to ParseTime\n");
487 0 : errSprintf ("(Probably a corrupt file)\n");
488 0 : return -1;
489 : }
490 20 : pds += 5;
491 20 : pdsMeta->timeRange = pds[3];
492 20 : if (ParseSect4Time2secV1 (pds[1], *pds, &P1_DeltaTime) == 0) {
493 20 : pdsMeta->P1 = pdsMeta->refTime + P1_DeltaTime;
494 : } else {
495 0 : pdsMeta->P1 = pdsMeta->refTime;
496 0 : printf ("Warning! : Can't figure out time unit of %d\n", *pds);
497 : }
498 20 : if (ParseSect4Time2secV1 (pds[2], *pds, &P2_DeltaTime) == 0) {
499 20 : pdsMeta->P2 = pdsMeta->refTime + P2_DeltaTime;
500 : } else {
501 0 : pdsMeta->P2 = pdsMeta->refTime;
502 0 : printf ("Warning! : Can't figure out time unit of %d\n", *pds);
503 : }
504 : /* The following is based on Table 5. */
505 : /* Note: For ensemble forecasts, 119 has meaning. */
506 20 : switch (pdsMeta->timeRange) {
507 : case 0:
508 : case 1:
509 : case 113:
510 : case 114:
511 : case 115:
512 : case 116:
513 : case 117:
514 : case 118:
515 : case 123:
516 : case 124:
517 16 : pdsMeta->validTime = pdsMeta->P1;
518 16 : break;
519 : case 2:
520 : /* Puzzling case. */
521 0 : pdsMeta->validTime = pdsMeta->P2;
522 0 : break;
523 : case 3:
524 : case 4:
525 : case 5:
526 : case 51:
527 0 : pdsMeta->validTime = pdsMeta->P2;
528 0 : break;
529 : case 10:
530 4 : if (ParseSect4Time2secV1 (GRIB_UNSIGN_INT2 (pds[1], pds[2]), *pds,
531 : &P1_DeltaTime) == 0) {
532 4 : pdsMeta->P2 = pdsMeta->P1 = pdsMeta->refTime + P1_DeltaTime;
533 : } else {
534 0 : pdsMeta->P2 = pdsMeta->P1 = pdsMeta->refTime;
535 0 : printf ("Warning! : Can't figure out time unit of %d\n", *pds);
536 : }
537 4 : pdsMeta->validTime = pdsMeta->P1;
538 4 : break;
539 : default:
540 0 : pdsMeta->validTime = pdsMeta->P1;
541 : }
542 20 : pds += 4;
543 20 : pdsMeta->Average = GRIB_UNSIGN_INT2 (*pds, pds[1]);
544 20 : pds += 2;
545 20 : pdsMeta->numberMissing = *(pds++);
546 : /* Skip over centry of reference time. */
547 20 : pds++;
548 20 : *subcenter = *(pds++);
549 20 : *DSF = GRIB_SIGN_INT2 (*pds, pds[1]);
550 20 : pds += 2;
551 20 : pdsMeta->f_hasEns = 0;
552 20 : pdsMeta->f_hasProb = 0;
553 20 : pdsMeta->f_hasCluster = 0;
554 20 : if (sectLen < 41) {
555 20 : return 0;
556 : }
557 : /* Following is based on:
558 : * http://www.emc.ncep.noaa.gov/gmb/ens/info/ens_grib.html */
559 0 : if ((*center == NMC) && (*subcenter == 2)) {
560 0 : if (sectLen < 45) {
561 0 : printf ("Warning! Problems with Ensemble section\n");
562 0 : return 0;
563 : }
564 0 : pdsMeta->f_hasEns = 1;
565 0 : pdsMeta->ens.BitFlag = *(pds++);
566 : /* octet21 = pdsMeta->timeRange; = 119 has meaning now */
567 0 : pds += 11;
568 0 : pdsMeta->ens.Application = *(pds++);
569 0 : pdsMeta->ens.Type = *(pds++);
570 0 : pdsMeta->ens.Number = *(pds++);
571 0 : pdsMeta->ens.ProdID = *(pds++);
572 0 : pdsMeta->ens.Smooth = *(pds++);
573 0 : if ((pdsMeta->cat == 191) || (pdsMeta->cat == 192) ||
574 : (pdsMeta->cat == 193)) {
575 0 : if (sectLen < 60) {
576 0 : printf ("Warning! Problems with Ensemble Probability section\n");
577 0 : return 0;
578 : }
579 0 : pdsMeta->f_hasProb = 1;
580 0 : pdsMeta->prob.Cat = pdsMeta->cat;
581 0 : pdsMeta->cat = *(pds++);
582 0 : pdsMeta->prob.Type = *(pds++);
583 0 : MEMCPY_BIG (&uli_temp, pds, sizeof (sInt4));
584 0 : pdsMeta->prob.lower = fval_360 (uli_temp);
585 0 : pds += 4;
586 0 : MEMCPY_BIG (&uli_temp, pds, sizeof (sInt4));
587 0 : pdsMeta->prob.upper = fval_360 (uli_temp);
588 0 : pds += 4;
589 0 : pds += 4;
590 : }
591 0 : if ((pdsMeta->ens.Type == 4) || (pdsMeta->ens.Type == 5)) {
592 : /* 87 ... 100 was reserved, but may not be encoded */
593 0 : if ((sectLen < 100) && (sectLen != 86)) {
594 0 : printf ("Warning! Problems with Ensemble Clustering section\n");
595 0 : printf ("Section length == %d\n", sectLen);
596 0 : return 0;
597 : }
598 0 : if (pdsMeta->f_hasProb == 0) {
599 0 : pds += 14;
600 : }
601 0 : pdsMeta->f_hasCluster = 1;
602 0 : pdsMeta->cluster.ensSize = *(pds++);
603 0 : pdsMeta->cluster.clusterSize = *(pds++);
604 0 : pdsMeta->cluster.Num = *(pds++);
605 0 : pdsMeta->cluster.Method = *(pds++);
606 0 : pdsMeta->cluster.NorLat = GRIB_UNSIGN_INT3 (*pds, pds[1], pds[2]);
607 0 : pdsMeta->cluster.NorLat = pdsMeta->cluster.NorLat / 1000.;
608 0 : pds += 3;
609 0 : pdsMeta->cluster.SouLat = GRIB_UNSIGN_INT3 (*pds, pds[1], pds[2]);
610 0 : pdsMeta->cluster.SouLat = pdsMeta->cluster.SouLat / 1000.;
611 0 : pds += 3;
612 0 : pdsMeta->cluster.EasLon = GRIB_UNSIGN_INT3 (*pds, pds[1], pds[2]);
613 0 : pdsMeta->cluster.EasLon = pdsMeta->cluster.EasLon / 1000.;
614 0 : pds += 3;
615 0 : pdsMeta->cluster.WesLon = GRIB_UNSIGN_INT3 (*pds, pds[1], pds[2]);
616 0 : pdsMeta->cluster.WesLon = pdsMeta->cluster.WesLon / 1000.;
617 0 : pds += 3;
618 0 : memcpy (pdsMeta->cluster.Member, pds, 10);
619 0 : pdsMeta->cluster.Member[10] = '\0';
620 : }
621 : /* Following based on:
622 : * http://www.ecmwf.int/publications/manuals/libraries/gribex/
623 : * localGRIBUsage.html */
624 0 : } else if (*center == ECMWF) {
625 0 : if (sectLen < 45) {
626 0 : printf ("Warning! Problems with ECMWF PDS extension\n");
627 0 : return 0;
628 : }
629 : /*
630 : sInt4 i_temp;
631 : pds += 12;
632 : i_temp = GRIB_SIGN_INT2 (pds[3], pds[4]);
633 : printf ("ID %d Class %d Type %d Stream %d", pds[0], pds[1], pds[2],
634 : i_temp);
635 : pds += 5;
636 : printf (" Ver %c%c%c%c, ", pds[0], pds[1], pds[2], pds[3]);
637 : pds += 4;
638 : printf ("Octet-50 %d, Octet-51 %d SectLen %d\n", pds[0], pds[1],
639 : sectLen);
640 : */
641 : } else {
642 : printf ("Un-handled possible ensemble section center %d "
643 0 : "subcenter %d\n", *center, *subcenter);
644 : }
645 0 : return 0;
646 : }
647 :
648 : /*****************************************************************************
649 : * Grib1_Inventory() --
650 : *
651 : * Arthur Taylor / MDL
652 : *
653 : * PURPOSE
654 : * Parses the GRIB1 "Product Definition Section" for enough information to
655 : * fill out the inventory data structure so we can do a simple inventory on
656 : * the file in a similar way to how we did it for GRIB2.
657 : *
658 : * ARGUMENTS
659 : * fp = An opened GRIB2 file already at the correct message. (Input)
660 : * gribLen = The total length of the GRIB1 message. (Input)
661 : * inv = The inventory data structure that we need to fill. (Output)
662 : *
663 : * FILES/DATABASES: None
664 : *
665 : * RETURNS: int (could use errSprintf())
666 : * 0 = OK
667 : * -1 = gribLen is too small.
668 : *
669 : * HISTORY
670 : * 4/2003 Arthur Taylor (MDL/RSIS): Created.
671 : *
672 : * NOTES
673 : *****************************************************************************
674 : */
675 16 : int GRIB1_Inventory (DataSource &fp, uInt4 gribLen, inventoryType *inv)
676 : {
677 : char temp[3]; /* Used to determine the section length. */
678 : uInt4 sectLen; /* Length in bytes of the current section. */
679 : uChar *pds; /* The part of the message dealing with the PDS. */
680 : pdsG1Type pdsMeta; /* The pds parsed into a usable data structure. */
681 : char f_gds; /* flag if there is a gds section. */
682 : char f_bms; /* flag if there is a bms section. */
683 : short int DSF; /* Decimal Scale Factor for unpacking the data. */
684 : uChar gridID; /* Which GDS specs to use. */
685 : const char *varName; /* The name of the data stored in the grid. */
686 : const char *varComment; /* Extra comments about the data stored in grid. */
687 : const char *varUnit; /* Holds the name of the unit [K] [%] .. etc */
688 : int convert; /* Conversion method for this variable's unit. */
689 : uInt4 curLoc; /* Where we are in the current GRIB message. */
690 : unsigned short int center; /* The Center that created the data */
691 : unsigned short int subcenter; /* The Sub Center that created the data */
692 :
693 16 : curLoc = 8;
694 16 : if (fp.DataSourceFread(temp, sizeof (char), 3) != 3) {
695 0 : errSprintf ("Ran out of file.\n");
696 0 : return -1;
697 : }
698 16 : sectLen = GRIB_UNSIGN_INT3 (*temp, temp[1], temp[2]);
699 16 : if (curLoc + sectLen > gribLen) {
700 0 : errSprintf ("Ran out of data in PDS (GRIB1_Inventory)\n");
701 0 : return -1;
702 : }
703 16 : pds = (uChar *) malloc (sectLen * sizeof (uChar));
704 16 : *pds = *temp;
705 16 : pds[1] = temp[1];
706 16 : pds[2] = temp[2];
707 16 : if (fp.DataSourceFread(pds + 3, sizeof (char), sectLen - 3) + 3 != sectLen) {
708 0 : errSprintf ("Ran out of file.\n");
709 0 : free (pds);
710 0 : return -1;
711 : }
712 :
713 16 : if (ReadGrib1Sect1 (pds, gribLen, &curLoc, &pdsMeta, &f_gds, &gridID,
714 : &f_bms, &DSF, ¢er, &subcenter) != 0) {
715 0 : preErrSprintf ("Inside GRIB1_Inventory\n");
716 0 : free (pds);
717 0 : return -1;
718 : }
719 16 : free (pds);
720 16 : inv->refTime = pdsMeta.refTime;
721 16 : inv->validTime = pdsMeta.validTime;
722 16 : inv->foreSec = inv->validTime - inv->refTime;
723 : GRIB1_Table2LookUp (&(pdsMeta), &varName, &varComment, &varUnit,
724 16 : &convert, center, subcenter);
725 16 : inv->element = (char *) malloc ((1 + strlen (varName)) * sizeof (char));
726 16 : strcpy (inv->element, varName);
727 : inv->unitName = (char *) malloc ((1 + 2 + strlen (varUnit)) *
728 16 : sizeof (char));
729 16 : sprintf (inv->unitName, "[%s]", varUnit);
730 : inv->comment = (char *) malloc ((1 + strlen (varComment) +
731 : strlen (varUnit) + 2 + 1) *
732 16 : sizeof (char));
733 16 : sprintf (inv->comment, "%s [%s]", varComment, varUnit);
734 :
735 : GRIB1_Table3LookUp (&(pdsMeta), &(inv->shortFstLevel),
736 16 : &(inv->longFstLevel));
737 :
738 : /* Get to the end of the GRIB1 message. */
739 : /* (inventory.c : GRIB2Inventory), is responsible for this. */
740 : /* fseek (fp, gribLen - sectLen, SEEK_CUR); */
741 16 : return 0;
742 : }
743 :
744 0 : int GRIB1_RefTime (DataSource &fp, uInt4 gribLen, double *refTime)
745 : {
746 : char temp[3]; /* Used to determine the section length. */
747 : uInt4 sectLen; /* Length in bytes of the current section. */
748 : uChar *pds; /* The part of the message dealing with the PDS. */
749 : pdsG1Type pdsMeta; /* The pds parsed into a usable data structure. */
750 : char f_gds; /* flag if there is a gds section. */
751 : char f_bms; /* flag if there is a bms section. */
752 : short int DSF; /* Decimal Scale Factor for unpacking the data. */
753 : uChar gridID; /* Which GDS specs to use. */
754 : uInt4 curLoc; /* Where we are in the current GRIB message. */
755 : unsigned short int center; /* The Center that created the data */
756 : unsigned short int subcenter; /* The Sub Center that created the data */
757 :
758 0 : curLoc = 8;
759 0 : if (fp.DataSourceFread (temp, sizeof (char), 3) != 3) {
760 0 : errSprintf ("Ran out of file.\n");
761 0 : return -1;
762 : }
763 0 : sectLen = GRIB_UNSIGN_INT3 (*temp, temp[1], temp[2]);
764 0 : if (curLoc + sectLen > gribLen) {
765 0 : errSprintf ("Ran out of data in PDS (GRIB1_Inventory)\n");
766 0 : return -1;
767 : }
768 0 : pds = (uChar *) malloc (sectLen * sizeof (uChar));
769 0 : *pds = *temp;
770 0 : pds[1] = temp[1];
771 0 : pds[2] = temp[2];
772 0 : if (fp.DataSourceFread (pds + 3, sizeof (char), sectLen - 3) + 3 != sectLen) {
773 0 : errSprintf ("Ran out of file.\n");
774 0 : free (pds);
775 0 : return -1;
776 : }
777 :
778 0 : if (ReadGrib1Sect1 (pds, gribLen, &curLoc, &pdsMeta, &f_gds, &gridID,
779 : &f_bms, &DSF, ¢er, &subcenter) != 0) {
780 0 : preErrSprintf ("Inside GRIB1_Inventory\n");
781 0 : free (pds);
782 0 : return -1;
783 : }
784 0 : free (pds);
785 :
786 0 : *refTime = pdsMeta.refTime;
787 :
788 : /* Get to the end of the GRIB1 message. */
789 : /* (inventory.c : GRIB2Inventory), is responsible for this. */
790 : /* fseek (fp, gribLen - sectLen, SEEK_CUR); */
791 0 : return 0;
792 : }
793 :
794 : /*****************************************************************************
795 : * ReadGrib1Sect2() --
796 : *
797 : * Arthur Taylor / MDL
798 : *
799 : * PURPOSE
800 : * Parses the GRIB1 "Grid Definition Section" or section 2, filling out
801 : * the gdsMeta data structure.
802 : *
803 : * ARGUMENTS
804 : * gds = The compressed part of the message dealing with "GDS". (Input)
805 : * gribLen = The total length of the GRIB1 message. (Input)
806 : * curLoc = Current location in the GRIB1 message. (Output)
807 : * gdsMeta = The filled out gdsMeta data structure. (Output)
808 : *
809 : * FILES/DATABASES: None
810 : *
811 : * RETURNS: int (could use errSprintf())
812 : * 0 = OK
813 : * -1 = gribLen is too small.
814 : * -2 = unexpected values in gds.
815 : *
816 : * HISTORY
817 : * 4/2003 Arthur Taylor (MDL/RSIS): Created.
818 : * 12/2003 AAT: adas data encoder seems to have # of vertical data = 1, but
819 : * parameters of vertical data = 255, which doesn't make sense.
820 : * Changed the error from "fatal" to a warning in debug mode.
821 : * 6/2004 AAT: Modified to allow "extended" lat/lon grids (ie stretched or
822 : * stretched and rotated).
823 : *
824 : * NOTES
825 : *****************************************************************************
826 : */
827 4 : static int ReadGrib1Sect2 (uChar *gds, uInt4 gribLen, uInt4 *curLoc,
828 : gdsType *gdsMeta)
829 : {
830 : uInt4 sectLen; /* Length in bytes of the current section. */
831 : int gridType; /* Which type of grid. (see enumerated types). */
832 4 : double unit = 1e-3; /* Used for converting to the correct unit. */
833 : uInt4 uli_temp; /* Used for reading a GRIB1 float. */
834 : int i;
835 : int f_allZero; /* Used to find out if the "lat/lon" extension part
836 : * is all 0 hence missing. */
837 : int f_allOne; /* Used to find out if the "lat/lon" extension part
838 : * is all 1 hence missing. */
839 :
840 4 : sectLen = GRIB_UNSIGN_INT3 (*gds, gds[1], gds[2]);
841 : #ifdef DEBUG
842 : /*
843 : printf ("Section 2 length = %ld\n", sectLen);
844 : */
845 : #endif
846 4 : *curLoc += sectLen;
847 4 : if (*curLoc > gribLen) {
848 0 : errSprintf ("Ran out of data in GDS (GRIB 1 Section 2)\n");
849 0 : return -1;
850 : }
851 4 : gds += 3;
852 : /*
853 : #ifdef DEBUG
854 : if ((*gds != 0) || (gds[1] != 255)) {
855 : printf ("GRIB1 GDS: Expect (NV = 0) != %d, (PV = 255) != %d\n",
856 : *gds, gds[1]);
857 : errSprintf ("SectLen == %ld\n", sectLen);
858 : errSprintf ("GridType == %d\n", gds[2]);
859 : }
860 : #endif
861 : */
862 : #ifdef DEBUG
863 : if (gds[1] != 255) {
864 : printf ("\n\tCaution: GRIB1 GDS: FOR ALL NWS products, PV should be "
865 : "255 rather than %d\n", gds[1]);
866 : }
867 : #endif
868 4 : if ((gds[1] != 255) && (gds[1] > 6)) {
869 0 : errSprintf ("GRIB1 GDS: Expect PV = 255 != %d\n", gds[1]);
870 0 : return -2;
871 : }
872 4 : gds += 2;
873 4 : gridType = *(gds++);
874 4 : switch (gridType) {
875 : case GB1S2_LATLON: // Latitude/Longitude Grid
876 : case GB1S2_GAUSSIAN_LATLON: // Gaussian Latitude/Longitude
877 : case GB1S2_ROTATED_LATLON: // Rotated Latitude/Longitude
878 4 : if ((sectLen != 32) && (sectLen != 42) && (sectLen != 52)) {
879 : errSprintf ("For LatLon GDS, should have 32 or 42 or 52 bytes "
880 0 : "of data\n");
881 0 : return -1;
882 : }
883 4 : switch(gridType) {
884 : case GB1S2_GAUSSIAN_LATLON:
885 0 : gdsMeta->projType = GS3_GAUSSIAN_LATLON;
886 0 : break;
887 : case GB1S2_ROTATED_LATLON:
888 0 : gdsMeta->projType = GS3_ROTATED_LATLON;
889 0 : break;
890 : default:
891 4 : gdsMeta->projType = GS3_LATLON;
892 : break;
893 : }
894 4 : gdsMeta->orientLon = 0;
895 4 : gdsMeta->meshLat = 0;
896 4 : gdsMeta->scaleLat1 = 0;
897 4 : gdsMeta->scaleLat2 = 0;
898 4 : gdsMeta->southLat = 0;
899 4 : gdsMeta->southLon = 0;
900 4 : gdsMeta->center = 0;
901 :
902 4 : gdsMeta->Nx = GRIB_UNSIGN_INT2 (*gds, gds[1]);
903 4 : gds += 2;
904 4 : gdsMeta->Ny = GRIB_UNSIGN_INT2 (*gds, gds[1]);
905 4 : gds += 2;
906 4 : gdsMeta->lat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
907 4 : gds += 3;
908 4 : gdsMeta->lon1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
909 4 : gds += 3;
910 :
911 4 : gdsMeta->resFlag = *(gds++);
912 4 : if (gdsMeta->resFlag & 0x40) {
913 0 : gdsMeta->f_sphere = 0;
914 0 : gdsMeta->majEarth = 6378.160;
915 0 : gdsMeta->minEarth = 6356.775;
916 : } else {
917 4 : gdsMeta->f_sphere = 1;
918 4 : gdsMeta->majEarth = 6367.47;
919 4 : gdsMeta->minEarth = 6367.47;
920 : }
921 :
922 4 : gdsMeta->lat2 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
923 4 : gds += 3;
924 4 : gdsMeta->lon2 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
925 4 : gds += 3;
926 4 : gdsMeta->Dx = GRIB_UNSIGN_INT2 (*gds, gds[1]) * unit;
927 4 : gds += 2;
928 4 : if (gridType == GB1S2_GAUSSIAN_LATLON) {
929 0 : int np = GRIB_UNSIGN_INT2 (*gds, gds[1]); /* parallels between a pole and the equator */
930 0 : gdsMeta->Dy = 90.0 / np;
931 : } else
932 4 : gdsMeta->Dy = GRIB_UNSIGN_INT2 (*gds, gds[1]) * unit;
933 4 : gds += 2;
934 4 : gdsMeta->scan = *gds;
935 4 : gdsMeta->f_typeLatLon = 0;
936 : #ifdef DEBUG
937 : /*
938 : printf ("sectLen %ld\n", sectLen);
939 : */
940 : #endif
941 4 : if (sectLen == 42) {
942 : /* Check if all 0's or all 1's, which means f_typeLatLon == 0 */
943 0 : f_allZero = 1;
944 0 : f_allOne = 1;
945 0 : for (i = 0; i < 10; i++) {
946 0 : if (gds[i] != 0)
947 0 : f_allZero = 0;
948 0 : if (gds[i] != 255)
949 0 : f_allOne = 0;
950 : }
951 0 : if (!f_allZero && !f_allOne) {
952 0 : gdsMeta->f_typeLatLon = 1;
953 0 : gds += 5;
954 0 : gdsMeta->poleLat = (GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) *
955 0 : unit);
956 0 : gds += 3;
957 0 : gdsMeta->poleLon = (GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) *
958 0 : unit);
959 0 : gds += 3;
960 0 : MEMCPY_BIG (&uli_temp, gds, sizeof (sInt4));
961 0 : gdsMeta->stretchFactor = fval_360 (uli_temp);
962 : }
963 4 : } else if (sectLen == 52) {
964 0 : gds += 5;
965 : /* Check if all 0's or all 1's, which means f_typeLatLon == 0 */
966 0 : f_allZero = 1;
967 0 : f_allOne = 1;
968 0 : for (i = 0; i < 20; i++) {
969 0 : if (gds[i] != 0)
970 0 : f_allZero = 0;
971 0 : if (gds[i] != 255)
972 0 : f_allOne = 0;
973 : }
974 0 : if (!f_allZero && !f_allOne) {
975 0 : gdsMeta->f_typeLatLon = 2;
976 0 : gdsMeta->southLat = (GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) *
977 0 : unit);
978 0 : gds += 3;
979 0 : gdsMeta->southLon = (GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) *
980 0 : unit);
981 0 : gds += 3;
982 0 : MEMCPY_BIG (&uli_temp, gds, sizeof (sInt4));
983 0 : gdsMeta->angleRotate = fval_360 (uli_temp);
984 0 : gds += 4;
985 0 : gdsMeta->poleLat = (GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) *
986 0 : unit);
987 0 : gds += 3;
988 0 : gdsMeta->poleLon = (GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) *
989 0 : unit);
990 0 : gds += 3;
991 0 : MEMCPY_BIG (&uli_temp, gds, sizeof (sInt4));
992 0 : gdsMeta->stretchFactor = fval_360 (uli_temp);
993 : }
994 : #ifdef DEBUG
995 : /*
996 : if (gdsMeta->lon2 == 360.25)
997 : gdsMeta->lon2 = 359.75;
998 : */
999 : /*
1000 : printf ("south %f %f rotate %f pole %f %f stretch %f\n",
1001 : gdsMeta->southLat, gdsMeta->southLon,
1002 : gdsMeta->angleRotate, gdsMeta->poleLat, gdsMeta->poleLon,
1003 : gdsMeta->stretchFactor);
1004 : printf ("lat/lon type %d \n", gdsMeta->f_typeLatLon);
1005 : */
1006 : #endif
1007 : }
1008 4 : break;
1009 :
1010 : case GB1S2_POLAR:
1011 0 : if (sectLen != 32) {
1012 0 : errSprintf ("For Polar GDS, should have 32 bytes of data\n");
1013 0 : return -1;
1014 : }
1015 0 : gdsMeta->projType = GS3_POLAR;
1016 0 : gdsMeta->lat2 = 0;
1017 0 : gdsMeta->lon2 = 0;
1018 0 : gdsMeta->southLat = 0;
1019 0 : gdsMeta->southLon = 0;
1020 :
1021 0 : gdsMeta->Nx = GRIB_UNSIGN_INT2 (*gds, gds[1]);
1022 0 : gds += 2;
1023 0 : gdsMeta->Ny = GRIB_UNSIGN_INT2 (*gds, gds[1]);
1024 0 : gds += 2;
1025 0 : gdsMeta->lat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1026 0 : gds += 3;
1027 0 : gdsMeta->lon1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1028 0 : gds += 3;
1029 :
1030 0 : gdsMeta->resFlag = *(gds++);
1031 0 : if (gdsMeta->resFlag & 0x40) {
1032 0 : gdsMeta->f_sphere = 0;
1033 0 : gdsMeta->majEarth = 6378.160;
1034 0 : gdsMeta->minEarth = 6356.775;
1035 : } else {
1036 0 : gdsMeta->f_sphere = 1;
1037 0 : gdsMeta->majEarth = 6367.47;
1038 0 : gdsMeta->minEarth = 6367.47;
1039 : }
1040 :
1041 0 : gdsMeta->orientLon = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1042 0 : gds += 3;
1043 0 : gdsMeta->Dx = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]);
1044 0 : gds += 3;
1045 0 : gdsMeta->Dy = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]);
1046 0 : gds += 3;
1047 0 : gdsMeta->meshLat = 60; /* Depends on hemisphere. */
1048 0 : gdsMeta->center = *(gds++);
1049 0 : if (gdsMeta->center & GRIB2BIT_1) {
1050 : /* South polar stereographic. */
1051 0 : gdsMeta->scaleLat1 = gdsMeta->scaleLat2 = -90;
1052 : } else {
1053 : /* North polar stereographic. */
1054 0 : gdsMeta->scaleLat1 = gdsMeta->scaleLat2 = 90;
1055 : }
1056 0 : gdsMeta->scan = *gds;
1057 0 : break;
1058 :
1059 : case GB1S2_LAMBERT:
1060 0 : if (sectLen != 42) {
1061 0 : errSprintf ("For Lambert GDS, should have 42 bytes of data\n");
1062 0 : return -1;
1063 : }
1064 0 : gdsMeta->projType = GS3_LAMBERT;
1065 0 : gdsMeta->lat2 = 0;
1066 0 : gdsMeta->lon2 = 0;
1067 :
1068 0 : gdsMeta->Nx = GRIB_UNSIGN_INT2 (*gds, gds[1]);
1069 0 : gds += 2;
1070 0 : gdsMeta->Ny = GRIB_UNSIGN_INT2 (*gds, gds[1]);
1071 0 : gds += 2;
1072 0 : gdsMeta->lat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1073 0 : gds += 3;
1074 0 : gdsMeta->lon1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1075 0 : gds += 3;
1076 :
1077 0 : gdsMeta->resFlag = *(gds++);
1078 0 : if (gdsMeta->resFlag & 0x40) {
1079 0 : gdsMeta->f_sphere = 0;
1080 0 : gdsMeta->majEarth = 6378.160;
1081 0 : gdsMeta->minEarth = 6356.775;
1082 : } else {
1083 0 : gdsMeta->f_sphere = 1;
1084 0 : gdsMeta->majEarth = 6367.47;
1085 0 : gdsMeta->minEarth = 6367.47;
1086 : }
1087 :
1088 0 : gdsMeta->orientLon = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1089 0 : gds += 3;
1090 0 : gdsMeta->Dx = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]);
1091 0 : gds += 3;
1092 0 : gdsMeta->Dy = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]);
1093 0 : gds += 3;
1094 0 : gdsMeta->center = *(gds++);
1095 0 : gdsMeta->scan = *(gds++);
1096 0 : gdsMeta->scaleLat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1097 0 : gds += 3;
1098 0 : gdsMeta->scaleLat2 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1099 0 : gds += 3;
1100 0 : gdsMeta->meshLat = gdsMeta->scaleLat1;
1101 0 : gdsMeta->southLat = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1102 0 : gds += 3;
1103 0 : gdsMeta->southLon = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1104 0 : break;
1105 :
1106 : case GB1S2_MERCATOR:
1107 0 : if (sectLen != 42) {
1108 0 : errSprintf ("For Mercator GDS, should have 42 bytes of data\n");
1109 0 : return -1;
1110 : }
1111 0 : gdsMeta->projType = GS3_MERCATOR;
1112 0 : gdsMeta->southLat = 0;
1113 0 : gdsMeta->southLon = 0;
1114 0 : gdsMeta->orientLon = 0;
1115 0 : gdsMeta->center = 0;
1116 :
1117 0 : gdsMeta->Nx = GRIB_UNSIGN_INT2 (*gds, gds[1]);
1118 0 : gds += 2;
1119 0 : gdsMeta->Ny = GRIB_UNSIGN_INT2 (*gds, gds[1]);
1120 0 : gds += 2;
1121 0 : gdsMeta->lat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1122 0 : gds += 3;
1123 0 : gdsMeta->lon1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1124 0 : gds += 3;
1125 :
1126 0 : gdsMeta->resFlag = *(gds++);
1127 0 : if (gdsMeta->resFlag & 0x40) {
1128 0 : gdsMeta->f_sphere = 0;
1129 0 : gdsMeta->majEarth = 6378.160;
1130 0 : gdsMeta->minEarth = 6356.775;
1131 : } else {
1132 0 : gdsMeta->f_sphere = 1;
1133 0 : gdsMeta->majEarth = 6367.47;
1134 0 : gdsMeta->minEarth = 6367.47;
1135 : }
1136 :
1137 0 : gdsMeta->lat2 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1138 0 : gds += 3;
1139 0 : gdsMeta->lon2 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1140 0 : gds += 3;
1141 0 : gdsMeta->scaleLat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) * unit;
1142 0 : gds += 3;
1143 0 : gdsMeta->scaleLat2 = gdsMeta->scaleLat1;
1144 0 : gdsMeta->meshLat = gdsMeta->scaleLat1;
1145 : /* Reserved set to 0. */
1146 0 : gds++;
1147 0 : gdsMeta->scan = *(gds++);
1148 0 : gdsMeta->Dx = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]);
1149 0 : gds += 3;
1150 0 : gdsMeta->Dy = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]);
1151 0 : break;
1152 :
1153 : default:
1154 0 : errSprintf ("Grid projection number is %d\n", gridType);
1155 0 : errSprintf ("Don't know how to handle this grid projection.\n");
1156 0 : return -2;
1157 : }
1158 4 : gdsMeta->numPts = gdsMeta->Nx * gdsMeta->Ny;
1159 : #ifdef DEBUG
1160 : /*
1161 : printf ("NumPts = %ld\n", gdsMeta->numPts);
1162 : printf ("Nx = %ld, Ny = %ld\n", gdsMeta->Nx, gdsMeta->Ny);
1163 : */
1164 : #endif
1165 4 : return 0;
1166 : }
1167 :
1168 : /*****************************************************************************
1169 : * ReadGrib1Sect3() --
1170 : *
1171 : * Arthur Taylor / MDL
1172 : *
1173 : * PURPOSE
1174 : * Parses the GRIB1 "Bit Map Section" or section 3, filling out the bitmap
1175 : * as needed.
1176 : *
1177 : * ARGUMENTS
1178 : * bms = The compressed part of the message dealing with "BMS". (Input)
1179 : * gribLen = The total length of the GRIB1 message. (Input)
1180 : * curLoc = Current location in the GRIB1 message. (Output)
1181 : * bitmap = The extracted bitmap. (Output)
1182 : * NxNy = The total size of the grid. (Input)
1183 : *
1184 : * FILES/DATABASES: None
1185 : *
1186 : * RETURNS: int (could use errSprintf())
1187 : * 0 = OK
1188 : * -1 = gribLen is too small.
1189 : * -2 = unexpected values in bms.
1190 : *
1191 : * HISTORY
1192 : * 4/2003 Arthur Taylor (MDL/RSIS): Created.
1193 : *
1194 : * NOTES
1195 : *****************************************************************************
1196 : */
1197 3 : static int ReadGrib1Sect3 (uChar *bms, uInt4 gribLen, uInt4 *curLoc,
1198 : uChar *bitmap, uInt4 NxNy)
1199 : {
1200 : uInt4 sectLen; /* Length in bytes of the current section. */
1201 : short int numeric; /* Determine if this is a predefined bitmap */
1202 : uChar bits; /* Used to locate which bit we are currently using. */
1203 : uInt4 i; /* Helps traverse the bitmap. */
1204 :
1205 3 : sectLen = GRIB_UNSIGN_INT3 (*bms, bms[1], bms[2]);
1206 : #ifdef DEBUG
1207 : /*
1208 : printf ("Section 3 length = %ld\n", sectLen);
1209 : */
1210 : #endif
1211 3 : *curLoc += sectLen;
1212 3 : if (*curLoc > gribLen) {
1213 0 : errSprintf ("Ran out of data in BMS (GRIB 1 Section 3)\n");
1214 0 : return -1;
1215 : }
1216 3 : bms += 3;
1217 : /* Assert: *bms currently points to number of unused bits at end of BMS. */
1218 3 : if (NxNy + *bms + 6 * 8 != sectLen * 8) {
1219 : errSprintf ("NxNy + # of unused bits %ld != # of available bits %ld\n",
1220 0 : (sInt4) (NxNy + *bms), (sInt4) ((sectLen - 6) * 8));
1221 0 : return -2;
1222 : }
1223 3 : bms++;
1224 : /* Assert: Non-zero "numeric" means predefined bitmap. */
1225 3 : numeric = GRIB_UNSIGN_INT2 (*bms, bms[1]);
1226 3 : bms += 2;
1227 3 : if (numeric != 0) {
1228 0 : errSprintf ("Don't handle predefined bitmaps yet.\n");
1229 0 : return -2;
1230 : }
1231 3 : bits = 0x80;
1232 17805 : for (i = 0; i < NxNy; i++) {
1233 17802 : *(bitmap++) = (*bms) & bits;
1234 17802 : bits = bits >> 1;
1235 17802 : if (bits == 0) {
1236 2224 : bms++;
1237 2224 : bits = 0x80;
1238 : }
1239 : }
1240 3 : return 0;
1241 : }
1242 :
1243 0 : static int UnpackCmplx (uChar *bds, uInt4 gribLen, uInt4 *curLoc,
1244 : short int DSF, double *data, grib_MetaData *meta,
1245 : char f_bms, uChar *bitmap, double unitM,
1246 : double unitB, short int ESF, double refVal,
1247 : uChar numBits, uChar f_octet14)
1248 : {
1249 : uInt4 secLen;
1250 : int N1;
1251 : int N2;
1252 : int P1;
1253 : int P2;
1254 : uChar octet14;
1255 : uChar f_maxtrixValues;
1256 0 : uChar f_secBitmap = 0;
1257 : uChar f_secValDiffWid;
1258 : int i;
1259 : uInt4 uli_temp; /* Used to store sInt4s (temporarily) */
1260 : uChar bufLoc; /* Keeps track of where to start getting more data
1261 : * out of the packed data stream. */
1262 : size_t numUsed; /* How many bytes were used in a given call to
1263 : * memBitRead. */
1264 : uChar *width;
1265 :
1266 0 : secLen = 11;
1267 0 : N1 = GRIB_UNSIGN_INT2 (bds[0], bds[1]);
1268 0 : octet14 = bds[2];
1269 0 : printf ("octet14, %d\n", octet14);
1270 0 : if (f_octet14) {
1271 0 : f_maxtrixValues = octet14 & GRIB2BIT_2;
1272 0 : f_secBitmap = octet14 & GRIB2BIT_3;
1273 0 : f_secValDiffWid = octet14 & GRIB2BIT_4;
1274 : printf ("f_matrixValues, f_secBitmap, f_secValeDiffWid %d %d %d\n",
1275 0 : f_maxtrixValues, f_secBitmap, f_secValDiffWid);
1276 : }
1277 0 : N2 = GRIB_UNSIGN_INT2 (bds[3], bds[4]);
1278 0 : P1 = GRIB_UNSIGN_INT2 (bds[5], bds[6]);
1279 0 : P2 = GRIB_UNSIGN_INT2 (bds[7], bds[8]);
1280 0 : printf ("N1 N2 P1 P2 : %d %d %d %d\n", N1, N2, P1, P2);
1281 0 : printf ("Reserved %d\n", bds[9]);
1282 0 : bds += 10;
1283 0 : secLen += 10;
1284 :
1285 0 : width = (uChar *) malloc (P1 * sizeof (uChar));
1286 :
1287 0 : for (i = 0; i < P1; i++) {
1288 0 : width[i] = *bds;
1289 0 : printf ("(Width %d %d)\n", i, width[i]);
1290 0 : bds++;
1291 0 : secLen++;
1292 : }
1293 0 : if (f_secBitmap) {
1294 0 : bufLoc = 8;
1295 0 : for (i = 0; i < P2; i++) {
1296 0 : memBitRead (&uli_temp, sizeof (sInt4), bds, 1, &bufLoc, &numUsed);
1297 0 : printf ("(%d %d) ", i, uli_temp);
1298 0 : if (numUsed != 0) {
1299 0 : printf ("\n");
1300 0 : bds += numUsed;
1301 0 : secLen++;
1302 : }
1303 : }
1304 0 : if (bufLoc != 8) {
1305 0 : bds++;
1306 0 : secLen++;
1307 : }
1308 0 : printf ("Observed Sec Len %d\n", secLen);
1309 : } else {
1310 : /* Jump over widths and secondary bitmap */
1311 0 : bds += (N1 - 21);
1312 0 : secLen += (N1 - 21);
1313 : }
1314 :
1315 0 : bufLoc = 8;
1316 0 : for (i = 0; i < P1; i++) {
1317 0 : memBitRead (&uli_temp, sizeof (sInt4), bds, numBits, &bufLoc, &numUsed);
1318 : printf ("(%d %d) (numUsed %ld numBits %d)", i, uli_temp,
1319 0 : (long) numUsed, numBits);
1320 0 : if (numUsed != 0) {
1321 0 : printf ("\n");
1322 0 : bds += numUsed;
1323 0 : secLen++;
1324 : }
1325 : }
1326 0 : if (bufLoc != 8) {
1327 0 : bds++;
1328 0 : secLen++;
1329 : }
1330 :
1331 0 : printf ("Observed Sec Len %d\n", secLen);
1332 0 : printf ("N2 = %d\n", N2);
1333 :
1334 0 : errSprintf ("Don't know how to handle Complex GRIB1 packing yet.\n");
1335 0 : free (width);
1336 0 : return -2;
1337 :
1338 : }
1339 :
1340 : /*****************************************************************************
1341 : * ReadGrib1Sect4() --
1342 : *
1343 : * Arthur Taylor / MDL
1344 : *
1345 : * PURPOSE
1346 : * Unpacks the "Binary Data Section" or section 4.
1347 : *
1348 : * ARGUMENTS
1349 : * bds = The compressed part of the message dealing with "BDS". (Input)
1350 : * gribLen = The total length of the GRIB1 message. (Input)
1351 : * curLoc = Current location in the GRIB1 message. (Output)
1352 : * DSF = Decimal Scale Factor for unpacking the data. (Input)
1353 : * data = The extracted grid. (Output)
1354 : * meta = The meta data associated with the grid (Input/Output)
1355 : * f_bms = True if bitmap is to be used. (Input)
1356 : * bitmap = 0 if missing data, 1 if valid data. (Input)
1357 : * unitM = The M unit conversion value in equation y = Mx + B. (Input)
1358 : * unitB = The B unit conversion value in equation y = Mx + B. (Input)
1359 : *
1360 : * FILES/DATABASES: None
1361 : *
1362 : * RETURNS: int (could use errSprintf())
1363 : * 0 = OK
1364 : * -1 = gribLen is too small.
1365 : * -2 = unexpected values in bds.
1366 : *
1367 : * HISTORY
1368 : * 4/2003 Arthur Taylor (MDL/RSIS): Created
1369 : * 3/2004 AAT: Switched {# Pts * (# Bits in a Group) +
1370 : * # of unused bits != # of available bits} to a warning from an
1371 : * error.
1372 : *
1373 : * NOTES
1374 : * 1) See metaparse.c : ParseGrid()
1375 : * 2) Currently, only handles "Simple pack".
1376 : *****************************************************************************
1377 : */
1378 4 : static int ReadGrib1Sect4 (uChar *bds, uInt4 gribLen, uInt4 *curLoc,
1379 : short int DSF, double *data, grib_MetaData *meta,
1380 : char f_bms, uChar *bitmap, double unitM,
1381 : double unitB)
1382 : {
1383 : uInt4 sectLen; /* Length in bytes of the current section. */
1384 : short int ESF; /* Power of 2 scaling factor. */
1385 : uInt4 uli_temp; /* Used to store sInt4s (temporarily) */
1386 : double refVal; /* The refrence value for the grid, also the minimum
1387 : * value. */
1388 : uChar numBits; /* # of bits for a single element of data. */
1389 : uChar numUnusedBit; /* # of extra bits at end of record. */
1390 : uChar f_spherHarm; /* Flag if data contains Spherical Harmonics. */
1391 : uChar f_cmplxPack; /* Flag if complex packing was used. */
1392 : uChar f_octet14; /* Flag if octet 14 was used. */
1393 : uChar bufLoc; /* Keeps track of where to start getting more data
1394 : * out of the packed data stream. */
1395 : uChar f_convert; /* Determine if scan mode implies that we have to do
1396 : * manipulation as we read the grid to get desired
1397 : * internal scan mode. */
1398 : uInt4 i; /* Used to traverse the grid. */
1399 : size_t numUsed; /* How many bytes were used in a given call to
1400 : * memBitRead. */
1401 : double d_temp; /* Holds the extracted data until we put it in data */
1402 : sInt4 newIndex; /* Where to put the answer (primarily if f_convert) */
1403 : sInt4 x; /* Used to help compute newIndex , if f_convert. */
1404 : sInt4 y; /* Used to help compute newIndex , if f_convert. */
1405 : double resetPrim; /* If possible, used to reset the primary missing
1406 : * value from 9.999e20 to a reasonable # (9999) */
1407 :
1408 4 : if (meta->gds.Nx * meta->gds.Ny != meta->gds.numPts) {
1409 0 : errSprintf ("(Nx * Ny != numPts) ?? in BDS (GRIB 1 Section 4)\n");
1410 0 : return -2;
1411 : }
1412 4 : sectLen = GRIB_UNSIGN_INT3 (*bds, bds[1], bds[2]);
1413 : #ifdef DEBUG
1414 : /*
1415 : printf ("Section 4 length = %ld\n", sectLen);
1416 : */
1417 : #endif
1418 4 : *curLoc += sectLen;
1419 4 : if (*curLoc > gribLen) {
1420 0 : errSprintf ("Ran out of data in BDS (GRIB 1 Section 4)\n");
1421 0 : return -1;
1422 : }
1423 4 : bds += 3;
1424 :
1425 : /* Assert: bds now points to the main pack flag. */
1426 4 : f_spherHarm = (*bds) & GRIB2BIT_1;
1427 4 : f_cmplxPack = (*bds) & GRIB2BIT_2;
1428 4 : meta->gridAttrib.fieldType = (*bds) & GRIB2BIT_3;
1429 4 : f_octet14 = (*bds) & GRIB2BIT_4;
1430 :
1431 4 : numUnusedBit = (*bds) & 0x0f;
1432 : #ifdef DEBUG
1433 : /*
1434 : printf ("bds byte flag = %d\n", *bds);
1435 : printf ("Number of unused bits = %d\n", numUnusedBit);
1436 : */
1437 : #endif
1438 4 : if (f_spherHarm) {
1439 0 : errSprintf ("Don't know how to handle Spherical Harmonics yet.\n");
1440 0 : return -2;
1441 : }
1442 : /*
1443 : if (f_octet14) {
1444 : errSprintf ("Don't know how to handle Octet 14 data yet.\n");
1445 : errSprintf ("bds byte flag = %d\n", *bds);
1446 : errSprintf ("bds byte: %d %d %d %d\n", f_spherHarm, f_cmplxPack,
1447 : meta->gridAttrib.fieldType, f_octet14);
1448 : return -2;
1449 : }
1450 : */
1451 4 : if (f_cmplxPack) {
1452 0 : meta->gridAttrib.packType = 2;
1453 : } else {
1454 4 : meta->gridAttrib.packType = 0;
1455 : }
1456 4 : bds++;
1457 :
1458 : /* Assert: bds now points to E (power of 2 scaling factor). */
1459 4 : ESF = GRIB_SIGN_INT2 (*bds, bds[1]);
1460 4 : bds += 2;
1461 4 : MEMCPY_BIG (&uli_temp, bds, sizeof (sInt4));
1462 4 : refVal = fval_360 (uli_temp);
1463 4 : bds += 4;
1464 :
1465 : /* Assert: bds is now the number of bits in a group. */
1466 4 : numBits = *bds;
1467 : /*
1468 : #ifdef DEBUG
1469 : printf ("refValue %f numBits %d\n", refVal, numBits);
1470 : printf ("ESF %d DSF %d\n", ESF, DSF);
1471 : #endif
1472 : */
1473 4 : if (f_cmplxPack) {
1474 0 : bds++;
1475 : #ifdef DEBUG
1476 : return UnpackCmplx (bds, gribLen, curLoc, DSF, data, meta, f_bms,
1477 : bitmap, unitM, unitB, ESF, refVal, numBits,
1478 : f_octet14);
1479 : #else
1480 0 : errSprintf ("Don't know how to handle Complex GRIB1 packing yet.\n");
1481 0 : return -2;
1482 : #endif
1483 : }
1484 :
1485 4 : if (!f_bms && (meta->gds.numPts * numBits + numUnusedBit) !=
1486 : (sectLen - 11) * 8) {
1487 : printf ("numPts * (numBits in a Group) + # of unused bits %d != "
1488 : "# of available bits %d\n",
1489 : (sInt4) (meta->gds.numPts * numBits + numUnusedBit),
1490 0 : (sInt4) ((sectLen - 11) * 8));
1491 : /*
1492 : errSprintf ("numPts * (numBits in a Group) + # of unused bits %ld != "
1493 : "# of available bits %ld\n",
1494 : (sInt4) (meta->gds.numPts * numBits + numUnusedBit),
1495 : (sInt4) ((sectLen - 11) * 8));
1496 : return -2;
1497 : */
1498 : }
1499 4 : if (numBits > 32) {
1500 0 : errSprintf ("The number of bits per number is larger than 32?\n");
1501 0 : return -2;
1502 : }
1503 4 : bds++;
1504 :
1505 : /* Convert Units. */
1506 4 : if (unitM == -10) {
1507 : meta->gridAttrib.min = pow (10.0, (refVal * pow (2.0, ESF) /
1508 0 : pow (10.0, DSF)));
1509 : } else {
1510 : /* meta->gridAttrib.min = unitM * (refVal / pow (10.0, DSF)) + unitB; */
1511 : meta->gridAttrib.min = unitM * (refVal * pow (2.0, ESF) /
1512 4 : pow (10.0, DSF)) + unitB;
1513 : }
1514 4 : meta->gridAttrib.max = meta->gridAttrib.min;
1515 4 : meta->gridAttrib.f_maxmin = 1;
1516 4 : meta->gridAttrib.numMiss = 0;
1517 4 : meta->gridAttrib.refVal = refVal;
1518 4 : meta->gridAttrib.ESF = ESF;
1519 4 : meta->gridAttrib.DSF = DSF;
1520 4 : bufLoc = 8;
1521 : /* Internally we use scan = 0100. Scan is usually 0100 but if need be, we
1522 : * can convert it. */
1523 4 : f_convert = ((meta->gds.scan & 0xe0) != 0x40);
1524 :
1525 4 : if (f_bms) {
1526 : /*
1527 : #ifdef DEBUG
1528 : printf ("There is a bitmap?\n");
1529 : #endif
1530 : */
1531 : /* Start unpacking the data, assuming there is a bitmap. */
1532 3 : meta->gridAttrib.f_miss = 1;
1533 3 : meta->gridAttrib.missPri = UNDEFINED;
1534 17805 : for (i = 0; i < meta->gds.numPts; i++) {
1535 : /* Find the destination index. */
1536 17802 : if (f_convert) {
1537 : /* ScanIndex2XY returns value as if scan was 0100 */
1538 : ScanIndex2XY (i, &x, &y, meta->gds.scan, meta->gds.Nx,
1539 0 : meta->gds.Ny);
1540 0 : newIndex = (x - 1) + (y - 1) * meta->gds.Nx;
1541 : } else {
1542 17802 : newIndex = i;
1543 : }
1544 : /* A 0 in bitmap means no data. A 1 in bitmap means data. */
1545 17802 : if (!bitmap[i]) {
1546 7542 : meta->gridAttrib.numMiss++;
1547 7542 : data[newIndex] = UNDEFINED;
1548 : } else {
1549 10260 : if (numBits != 0) {
1550 : memBitRead (&uli_temp, sizeof (sInt4), bds, numBits,
1551 10260 : &bufLoc, &numUsed);
1552 10260 : bds += numUsed;
1553 10260 : d_temp = (refVal + (uli_temp * pow (2.0, ESF))) / pow (10.0, DSF);
1554 : /* Convert Units. */
1555 10260 : if (unitM == -10) {
1556 0 : d_temp = pow (10.0, d_temp);
1557 : } else {
1558 10260 : d_temp = unitM * d_temp + unitB;
1559 : }
1560 10260 : if (meta->gridAttrib.max < d_temp) {
1561 16 : meta->gridAttrib.max = d_temp;
1562 : }
1563 10260 : data[newIndex] = d_temp;
1564 : } else {
1565 : /* Assert: d_temp = unitM * refVal / pow (10.0,DSF) + unitB. */
1566 : /* Assert: min = unitM * refVal / pow (10.0, DSF) + unitB. */
1567 0 : data[newIndex] = meta->gridAttrib.min;
1568 : }
1569 : }
1570 : }
1571 : /* Reset the missing value to UNDEFINED_PRIM if possible. If not
1572 : * possible, make sure UNDEFINED is outside the range. If UNDEFINED
1573 : * is_ in the range, choose max + 1 for missing. */
1574 3 : resetPrim = 0;
1575 5 : if ((meta->gridAttrib.max < UNDEFINED_PRIM) ||
1576 : (meta->gridAttrib.min > UNDEFINED_PRIM)) {
1577 2 : resetPrim = UNDEFINED_PRIM;
1578 1 : } else if ((meta->gridAttrib.max >= UNDEFINED) &&
1579 : (meta->gridAttrib.min <= UNDEFINED)) {
1580 0 : resetPrim = meta->gridAttrib.max + 1;
1581 : }
1582 3 : if (resetPrim != 0) {
1583 2 : meta->gridAttrib.missPri = resetPrim;
1584 12920 : for (i = 0; i < meta->gds.numPts; i++) {
1585 : /* Find the destination index. */
1586 12918 : if (f_convert) {
1587 : /* ScanIndex2XY returns value as if scan was 0100 */
1588 : ScanIndex2XY (i, &x, &y, meta->gds.scan, meta->gds.Nx,
1589 0 : meta->gds.Ny);
1590 0 : newIndex = (x - 1) + (y - 1) * meta->gds.Nx;
1591 : } else {
1592 12918 : newIndex = i;
1593 : }
1594 12918 : if (!bitmap[i]) {
1595 4852 : data[newIndex] = resetPrim;
1596 : }
1597 : }
1598 : }
1599 :
1600 : } else {
1601 :
1602 : #ifdef DEBUG
1603 : /*
1604 : printf ("There is no bitmap?\n");
1605 : */
1606 : #endif
1607 :
1608 : /* Start unpacking the data, assuming there is NO bitmap. */
1609 1 : meta->gridAttrib.f_miss = 0;
1610 589 : for (i = 0; i < meta->gds.numPts; i++) {
1611 588 : if (numBits != 0) {
1612 : /* Find the destination index. */
1613 588 : if (f_convert) {
1614 : /* ScanIndex2XY returns value as if scan was 0100 */
1615 : ScanIndex2XY (i, &x, &y, meta->gds.scan, meta->gds.Nx,
1616 588 : meta->gds.Ny);
1617 588 : newIndex = (x - 1) + (y - 1) * meta->gds.Nx;
1618 : } else {
1619 0 : newIndex = i;
1620 : }
1621 :
1622 : memBitRead (&uli_temp, sizeof (sInt4), bds, numBits, &bufLoc,
1623 588 : &numUsed);
1624 588 : bds += numUsed;
1625 588 : d_temp = (refVal + (uli_temp * pow (2.0, ESF))) / pow (10.0, DSF);
1626 :
1627 : #ifdef DEBUG
1628 : /*
1629 : if (i == 1) {
1630 : printf ("refVal %f, uli_temp %ld, ans %f\n", refVal, uli_temp,
1631 : d_temp);
1632 : printf ("numBits %d, bufLoc %d, numUsed %d\n", numBits,
1633 : bufLoc, numUsed);
1634 : }
1635 : */
1636 : #endif
1637 :
1638 : /* Convert Units. */
1639 588 : if (unitM == -10) {
1640 0 : d_temp = pow (10.0, d_temp);
1641 : } else {
1642 588 : d_temp = unitM * d_temp + unitB;
1643 : }
1644 588 : if (meta->gridAttrib.max < d_temp) {
1645 19 : meta->gridAttrib.max = d_temp;
1646 : }
1647 588 : data[newIndex] = d_temp;
1648 : } else {
1649 : /* Assert: whole array = unitM * refVal + unitB. */
1650 : /* Assert: *min = unitM * refVal + unitB. */
1651 0 : data[i] = meta->gridAttrib.min;
1652 : }
1653 : }
1654 : }
1655 4 : return 0;
1656 : }
1657 :
1658 : /*****************************************************************************
1659 : * ReadGrib1Record() --
1660 : *
1661 : * Arthur Taylor / MDL
1662 : *
1663 : * PURPOSE
1664 : * Reads in a GRIB1 message, and parses the data into various data
1665 : * structures, for use with other code.
1666 : *
1667 : * ARGUMENTS
1668 : * fp = An opened GRIB2 file already at the correct message. (Input)
1669 : * f_unit = 0 use GRIB2 units, 1 use English, 2 use metric. (Input)
1670 : * Grib_Data = The read in GRIB2 grid. (Output)
1671 : * grib_DataLen = Size of Grib_Data. (Output)
1672 : * meta = A filled in meta structure (Output)
1673 : * IS = The structure containing all the arrays that the
1674 : * unpacker uses (Output)
1675 : * sect0 = Already read in section 0 data. (Input)
1676 : * gribLen = Length of the GRIB1 message. (Input)
1677 : * majEarth = Used to override the GRIB major axis of earth. (Input)
1678 : * minEarth = Used to override the GRIB minor axis of earth. (Input)
1679 : *
1680 : * FILES/DATABASES:
1681 : * An already opened file pointing to the desired GRIB1 message.
1682 : *
1683 : * RETURNS: int (could use errSprintf())
1684 : * 0 = OK
1685 : * -1 = Problems reading in the PDS.
1686 : * -2 = Problems reading in the GDS.
1687 : * -3 = Problems reading in the BMS.
1688 : * -4 = Problems reading in the BDS.
1689 : * -5 = Problems reading the closing section.
1690 : *
1691 : * HISTORY
1692 : * 4/2003 Arthur Taylor (MDL/RSIS): Created
1693 : * 5/2003 AAT: Was not updating offset. It should be updated by
1694 : * calling routine anyways, so I got rid of the parameter.
1695 : * 7/2003 AAT: Allowed user to override the radius of earth.
1696 : * 8/2003 AAT: Found a memory Leak (Had been setting unitName to NULL).
1697 : * 2/2004 AAT: Added maj/min earth override.
1698 : * 3/2004 AAT: Added ability to change units.
1699 : *
1700 : * NOTES
1701 : * 1) Could also compare GDS with the one specified by gridID
1702 : * 2) Could add gridID support.
1703 : * 3) Should add unitM / unitB support.
1704 : *****************************************************************************
1705 : */
1706 4 : int ReadGrib1Record (DataSource &fp, sChar f_unit, double **Grib_Data,
1707 : uInt4 *grib_DataLen, grib_MetaData *meta,
1708 : IS_dataType *IS, sInt4 sect0[SECT0LEN_WORD],
1709 : uInt4 gribLen, double majEarth, double minEarth)
1710 : {
1711 : sInt4 nd5; /* Size of grib message rounded up to the nearest *
1712 : * sInt4. */
1713 : uChar *c_ipack; /* A char ptr to the message stored in IS->ipack */
1714 : uInt4 curLoc; /* Current location in the GRIB message. */
1715 : char f_gds; /* flag if there is a gds section. */
1716 : char f_bms; /* flag if there is a bms section. */
1717 : double *grib_Data; /* A pointer to Grib_Data for ease of manipulation. */
1718 4 : uChar *bitmap = NULL; /* A char field (0=noData, 1=data) set up in BMS. */
1719 : short int DSF; /* Decimal Scale Factor for unpacking the data. */
1720 4 : double unitM = 1; /* M in y = Mx + B, for unit conversion. */
1721 4 : double unitB = 0; /* B in y = Mx + B, for unit conversion. */
1722 : uChar gridID; /* Which GDS specs to use. */
1723 : const char *varName; /* The name of the data stored in the grid. */
1724 : const char *varComment; /* Extra comments about the data stored in grid. */
1725 : const char *varUnit; /* Holds the name of the unit [K] [%] .. etc */
1726 : sInt4 li_temp; /* Used to make sure section 5 is 7777. */
1727 : char unitName[15]; /* Holds the string name of the current unit. */
1728 : int unitLen; /* String length of string name of current unit. */
1729 :
1730 : /* Make room for entire message, and read it in. */
1731 : /* nd5 needs to be gribLen in (sInt4) units rounded up. */
1732 4 : nd5 = (gribLen + 3) / 4;
1733 4 : if (nd5 > IS->ipackLen) {
1734 4 : IS->ipackLen = nd5;
1735 : IS->ipack = (sInt4 *) realloc ((void *) (IS->ipack),
1736 4 : (IS->ipackLen) * sizeof (sInt4));
1737 : }
1738 4 : c_ipack = (uChar *) IS->ipack;
1739 : /* Init last sInt4 to 0, to make sure that the padded bytes are 0. */
1740 4 : IS->ipack[nd5 - 1] = 0;
1741 : /* Init first 2 sInt4 to sect0. */
1742 4 : memcpy (c_ipack, sect0, SECT0LEN_WORD * 2);
1743 : /* Read in the rest of the message. */
1744 4 : if (fp.DataSourceFread (c_ipack + SECT0LEN_WORD * 2, sizeof (char),
1745 4 : (gribLen - SECT0LEN_WORD * 2)) + SECT0LEN_WORD * 2 != gribLen) {
1746 0 : errSprintf ("Ran out of file\n");
1747 0 : return -1;
1748 : }
1749 :
1750 : /* Preceeding was in degrib2, next part is specific to GRIB1. */
1751 4 : curLoc = 8;
1752 4 : if (ReadGrib1Sect1 (c_ipack + curLoc, gribLen, &curLoc, &(meta->pds1),
1753 : &f_gds, &gridID, &f_bms, &DSF, &(meta->center),
1754 : &(meta->subcenter)) != 0) {
1755 0 : preErrSprintf ("Inside ReadGrib1Record\n");
1756 0 : return -1;
1757 : }
1758 :
1759 : /* Get the Grid Definition Section. */
1760 4 : if (f_gds) {
1761 4 : if (ReadGrib1Sect2 (c_ipack + curLoc, gribLen, &curLoc,
1762 : &(meta->gds)) != 0) {
1763 0 : preErrSprintf ("Inside ReadGrib1Record\n");
1764 0 : return -2;
1765 : }
1766 : /* Could also compare GDS with the one specified by gridID? */
1767 : } else {
1768 0 : errSprintf ("Don't know how to handle a gridID lookup yet.\n");
1769 0 : return -2;
1770 : }
1771 4 : meta->pds1.gridID = gridID;
1772 : /* Allow data originating from NCEP to be 6371.2 by default. */
1773 4 : if ((meta->center == NMC)) {
1774 4 : if (meta->gds.majEarth == 6367.47) {
1775 4 : meta->gds.f_sphere = 1;
1776 4 : meta->gds.majEarth = 6371.2;
1777 4 : meta->gds.minEarth = 6371.2;
1778 : }
1779 : }
1780 4 : if ((majEarth > 6300) && (majEarth < 6400)) {
1781 0 : if ((minEarth > 6300) && (minEarth < 6400)) {
1782 0 : meta->gds.f_sphere = 0;
1783 0 : meta->gds.majEarth = majEarth;
1784 0 : meta->gds.minEarth = minEarth;
1785 0 : if (majEarth == minEarth) {
1786 0 : meta->gds.f_sphere = 1;
1787 : }
1788 : } else {
1789 0 : meta->gds.f_sphere = 1;
1790 0 : meta->gds.majEarth = majEarth;
1791 0 : meta->gds.minEarth = majEarth;
1792 : }
1793 : }
1794 :
1795 : /* Allocate memory for the grid. */
1796 4 : if (meta->gds.numPts > *grib_DataLen) {
1797 4 : *grib_DataLen = meta->gds.numPts;
1798 : *Grib_Data = (double *) realloc ((void *) (*Grib_Data),
1799 4 : (*grib_DataLen) * sizeof (double));
1800 : }
1801 4 : grib_Data = *Grib_Data;
1802 :
1803 : /* Get the Bit Map Section. */
1804 4 : if (f_bms) {
1805 3 : bitmap = (uChar *) malloc (meta->gds.numPts * sizeof (char));
1806 3 : if (ReadGrib1Sect3 (c_ipack + curLoc, gribLen, &curLoc, bitmap,
1807 : meta->gds.numPts) != 0) {
1808 0 : free (bitmap);
1809 0 : preErrSprintf ("Inside ReadGrib1Record\n");
1810 0 : return -3;
1811 : }
1812 : }
1813 :
1814 : /* Figure out some basic stuff about the grid. */
1815 : /* Following is similar to metaparse.c : ParseElemName */
1816 : GRIB1_Table2LookUp (&(meta->pds1), &varName, &varComment, &varUnit,
1817 4 : &(meta->convert), meta->center, meta->subcenter);
1818 : meta->element = (char *) realloc ((void *) (meta->element),
1819 4 : (1 + strlen (varName)) * sizeof (char));
1820 4 : strcpy (meta->element, varName);
1821 : meta->unitName = (char *) realloc ((void *) (meta->unitName),
1822 : (1 + 2 + strlen (varUnit)) *
1823 4 : sizeof (char));
1824 4 : sprintf (meta->unitName, "[%s]", varUnit);
1825 : meta->comment = (char *) realloc ((void *) (meta->comment),
1826 : (1 + strlen (varComment) +
1827 : strlen (varUnit)
1828 4 : + 2 + 1) * sizeof (char));
1829 4 : sprintf (meta->comment, "%s [%s]", varComment, varUnit);
1830 :
1831 4 : if (ComputeUnit (meta->convert, meta->unitName, f_unit, &unitM, &unitB,
1832 : unitName) == 0) {
1833 0 : unitLen = strlen (unitName);
1834 : meta->unitName = (char *) realloc ((void *) (meta->unitName),
1835 0 : 1 + unitLen * sizeof (char));
1836 0 : strncpy (meta->unitName, unitName, unitLen);
1837 0 : meta->unitName[unitLen] = '\0';
1838 : }
1839 :
1840 : /* Read the GRID. */
1841 4 : if (ReadGrib1Sect4 (c_ipack + curLoc, gribLen, &curLoc, DSF, grib_Data,
1842 : meta, f_bms, bitmap, unitM, unitB) != 0) {
1843 0 : free (bitmap);
1844 0 : preErrSprintf ("Inside ReadGrib1Record\n");
1845 0 : return -4;
1846 : }
1847 4 : if (f_bms) {
1848 3 : free (bitmap);
1849 : }
1850 :
1851 : GRIB1_Table3LookUp (&(meta->pds1), &(meta->shortFstLevel),
1852 4 : &(meta->longFstLevel));
1853 : /* printf ("%s .. %s\n", meta->shortFstLevel, meta->longFstLevel);*/
1854 :
1855 : /*
1856 : strftime (meta->refTime, 20, "%Y%m%d%H%M",
1857 : gmtime (&(meta->pds1.refTime)));
1858 : */
1859 4 : Clock_Print (meta->refTime, 20, meta->pds1.refTime, "%Y%m%d%H%M", 0);
1860 :
1861 : /*
1862 : strftime (meta->validTime, 20, "%Y%m%d%H%M",
1863 : gmtime (&(meta->pds1.validTime)));
1864 : */
1865 4 : Clock_Print (meta->validTime, 20, meta->pds1.validTime, "%Y%m%d%H%M", 0);
1866 :
1867 4 : meta->deltTime = (sInt4) (meta->pds1.validTime - meta->pds1.refTime);
1868 :
1869 : /* Read section 5. If it is "7777" == 926365495 we are done. */
1870 4 : if (curLoc == gribLen) {
1871 : printf ("Warning: either gribLen did not account for section 5, or "
1872 0 : "section 5 is missing\n");
1873 0 : return 0;
1874 : }
1875 4 : if (curLoc + 4 > gribLen) {
1876 0 : errSprintf ("Ran out of bytes looking for the end of the message.\n");
1877 0 : return -5;
1878 : }
1879 : myAssert (curLoc + 4 == gribLen);
1880 4 : memcpy (&li_temp, c_ipack + curLoc, 4);
1881 4 : if (li_temp != 926365495L) {
1882 0 : errSprintf ("Did not find the end of the message.\n");
1883 0 : return -5;
1884 : }
1885 :
1886 4 : return 0;
1887 : }
1888 :
1889 : /*****************************************************************************
1890 : * main() --
1891 : *
1892 : * Arthur Taylor / MDL
1893 : *
1894 : * PURPOSE
1895 : * To test the capabilities of this module, and give an example as to how
1896 : * ReadGrib1Record expects to be called.
1897 : *
1898 : * ARGUMENTS
1899 : * argc = The number of arguments on the command line. (Input)
1900 : * argv = The arguments on the command line. (Input)
1901 : *
1902 : * FILES/DATABASES:
1903 : * A GRIB1 file.
1904 : *
1905 : * RETURNS: int
1906 : * 0 = OK
1907 : *
1908 : * HISTORY
1909 : * 4/2003 Arthur Taylor (MDL/RSIS): Created
1910 : * 6/2003 Matthew T. Kallio (matt@wunderground.com):
1911 : * "wmo" dimension increased to WMO_HEADER_LEN + 1 (for '\0' char)
1912 : *
1913 : * NOTES
1914 : *****************************************************************************
1915 : */
1916 : #ifdef DEBUG_DEGRIB1
1917 : int main (int argc, char **argv)
1918 : {
1919 : DataSource grib_fp; /* The opened grib2 file for input. */
1920 : sInt4 offset; /* Where we currently are in grib_fp. */
1921 : sInt4 sect0[SECT0LEN_WORD]; /* Holds the current Section 0. */
1922 : char wmo[WMO_HEADER_LEN + 1]; /* Holds the current wmo message. */
1923 : sInt4 gribLen; /* Length of the current GRIB message. */
1924 : sInt4 wmoLen; /* Length of current wmo Message. */
1925 : char *msg;
1926 : int version;
1927 : sChar f_unit = 0;
1928 : double *grib_Data;
1929 : sInt4 grib_DataLen;
1930 : grib_MetaData meta;
1931 : IS_dataType is; /* Un-parsed meta data for this GRIB2 message. As
1932 : * well as some memory used by the unpacker. */
1933 :
1934 : //if ((grib_fp = fopen (argv[1], "rb")) == NULL) {
1935 : // printf ("Problems opening %s for read\n", argv[1]);
1936 : // return 1;
1937 : //}
1938 : grib_fp = FileDataSource(argv[1]);
1939 : IS_Init (&is);
1940 : MetaInit (&meta);
1941 :
1942 : offset = 0;
1943 : if (ReadSECT0 (grib_fp, offset, WMO_HEADER_LEN, WMO_SECOND_LEN, wmo,
1944 : sect0, &gribLen, &wmoLen, &version) < 0) {
1945 : msg = errSprintf (NULL);
1946 : printf ("%s\n", msg);
1947 : return -1;
1948 : }
1949 : grib_DataLen = 0;
1950 : grib_Data = NULL;
1951 : if (version == 1) {
1952 : meta.GribVersion = version;
1953 : ReadGrib1Record (grib_fp, f_unit, &grib_Data, &grib_DataLen, &meta,
1954 : &is, sect0, gribLen);
1955 : offset = offset + gribLen + wmoLen;
1956 : }
1957 :
1958 : MetaFree (&meta);
1959 : IS_Free (&is);
1960 : free (grib_Data);
1961 : //fclose (grib_fp);
1962 : return 0;
1963 : }
1964 : #endif
|
