1 /*****************************************************************************
2 * metaparse.c
3 *
4 * DESCRIPTION
5 * This file contains the code necessary to initialize the meta data
6 * structure, and parse the meta data that comes out of the GRIB2 decoder.
7 *
8 * HISTORY
9 * 9/2002 Arthur Taylor (MDL / RSIS): Created.
10 *
11 * NOTES
12 * 1) Need to add support for GS3_ORTHOGRAPHIC = 90,
13 * GS3_EQUATOR_EQUIDIST = 110, GS3_AZIMUTH_RANGE = 120
14 * 2) Need to add support for GS4_RADAR = 20
15 *****************************************************************************
16 */
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <string.h>
20 #include <math.h>
21 #include <limits>
22 #include "clock.h"
23 #include "meta.h"
24 #include "metaname.h"
25 #include "myassert.h"
26 #include "myerror.h"
27 #include "scan.h"
28 #include "weather.h"
29 #include "hazard.h"
30 #include "tendian.h"
31 #include "myutil.h"
32
33 #include "cpl_string.h"
34
35 static void debug_printf(const char* fmt, ... ) CPL_PRINT_FUNC_FORMAT (1, 2);
36
debug_printf(const char * fmt,...)37 static void debug_printf(const char* fmt, ... )
38 {
39 va_list args;
40
41 va_start( args, fmt );
42 CPLDebug("GRIB", "%s", CPLString().vPrintf(fmt, args ).c_str() );
43 va_end( args );
44 }
45
46 #undef printf
47 #define printf debug_printf
48
49
50 /*****************************************************************************
51 * MetaInit() --
52 *
53 * Arthur Taylor / MDL
54 *
55 * PURPOSE
56 * To initialize a grib_metaData structure.
57 *
58 * ARGUMENTS
59 * meta = The structure to fill. (Output)
60 *
61 * FILES/DATABASES: None
62 *
63 * RETURNS: void
64 *
65 * HISTORY
66 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
67 *
68 * NOTES
69 *****************************************************************************
70 */
MetaInit(grib_MetaData * meta)71 void MetaInit (grib_MetaData *meta)
72 {
73 meta->element = nullptr;
74 meta->comment = nullptr;
75 meta->unitName = nullptr;
76 meta->convert = 0;
77 meta->shortFstLevel = nullptr;
78 meta->longFstLevel = nullptr;
79 meta->pds2.sect2.ptrType = GS2_NONE;
80
81 meta->pds2.sect2.wx.data = nullptr;
82 meta->pds2.sect2.wx.dataLen = 0;
83 meta->pds2.sect2.wx.maxLen = 0;
84 meta->pds2.sect2.wx.ugly = nullptr;
85 meta->pds2.sect2.unknown.data = nullptr;
86 meta->pds2.sect2.unknown.dataLen = 0;
87 meta->pds2.sect2.hazard.data = nullptr;
88 meta->pds2.sect2.hazard.dataLen = 0;
89 meta->pds2.sect2.hazard.maxLen = 0;
90 meta->pds2.sect2.hazard.haz = nullptr;
91
92 meta->pds2.sect4.numInterval = 0;
93 meta->pds2.sect4.Interval = nullptr;
94 meta->pds2.sect4.numBands = 0;
95 meta->pds2.sect4.bands = nullptr;
96 return;
97 }
98
99 /*****************************************************************************
100 * MetaSect2Free() --
101 *
102 * Arthur Taylor / MDL
103 *
104 * PURPOSE
105 * To free the section 2 data in the grib_metaData structure.
106 *
107 * ARGUMENTS
108 * meta = The structure to free. (Input/Output)
109 *
110 * FILES/DATABASES: None
111 *
112 * RETURNS: void
113 *
114 * HISTORY
115 * 2/2003 Arthur Taylor (MDL/RSIS): Created.
116 * 3/2003 AAT: Cleaned up declaration of variable: WxType.
117 *
118 * NOTES
119 *****************************************************************************
120 */
MetaSect2Free(grib_MetaData * meta)121 void MetaSect2Free (grib_MetaData *meta)
122 {
123 size_t i; /* Counter for use when freeing Wx data. */
124
125 if (meta->pds2.sect2.ptrType == GS2_WXTYPE) {
126 for (i = 0; i < meta->pds2.sect2.wx.dataLen; i++) {
127 free (meta->pds2.sect2.wx.data[i]);
128 FreeUglyString (&(meta->pds2.sect2.wx.ugly[i]));
129 }
130 free (meta->pds2.sect2.wx.ugly);
131 meta->pds2.sect2.wx.ugly = nullptr;
132 free (meta->pds2.sect2.wx.data);
133 meta->pds2.sect2.wx.data = nullptr;
134 free (meta->pds2.sect2.wx.f_valid);
135 meta->pds2.sect2.wx.f_valid = nullptr;
136 meta->pds2.sect2.wx.dataLen = 0;
137 meta->pds2.sect2.wx.maxLen = 0;
138 } else if (meta->pds2.sect2.ptrType == GS2_HAZARD) {
139 for (i = 0; i < meta->pds2.sect2.hazard.dataLen; i++) {
140 free (meta->pds2.sect2.hazard.data[i]);
141 FreeHazardString (&(meta->pds2.sect2.hazard.haz[i]));
142 }
143 free (meta->pds2.sect2.hazard.haz);
144 meta->pds2.sect2.hazard.haz = nullptr;
145 free (meta->pds2.sect2.hazard.data);
146 meta->pds2.sect2.hazard.data = nullptr;
147 free (meta->pds2.sect2.hazard.f_valid);
148 meta->pds2.sect2.hazard.f_valid = nullptr;
149 meta->pds2.sect2.hazard.dataLen = 0;
150 meta->pds2.sect2.hazard.maxLen = 0;
151 } else {
152 free (meta->pds2.sect2.unknown.data);
153 meta->pds2.sect2.unknown.data = nullptr;
154 meta->pds2.sect2.unknown.dataLen = 0;
155 }
156 meta->pds2.sect2.ptrType = GS2_NONE;
157 }
158
159 /*****************************************************************************
160 * MetaFree() --
161 *
162 * Arthur Taylor / MDL
163 *
164 * PURPOSE
165 * To free a grib_metaData structure.
166 *
167 * ARGUMENTS
168 * meta = The structure to free. (Output)
169 *
170 * FILES/DATABASES: None
171 *
172 * RETURNS: void
173 *
174 * HISTORY
175 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
176 *
177 * NOTES
178 *****************************************************************************
179 */
MetaFree(grib_MetaData * meta)180 void MetaFree (grib_MetaData *meta)
181 {
182 free (meta->pds2.sect4.bands);
183 meta->pds2.sect4.bands = nullptr;
184 meta->pds2.sect4.numBands = 0;
185 free (meta->pds2.sect4.Interval);
186 meta->pds2.sect4.Interval = nullptr;
187 meta->pds2.sect4.numInterval = 0;
188 MetaSect2Free (meta);
189 free (meta->unitName);
190 meta->unitName = nullptr;
191 meta->convert = 0;
192 free (meta->comment);
193 meta->comment = nullptr;
194 free (meta->element);
195 meta->element = nullptr;
196 free (meta->shortFstLevel);
197 meta->shortFstLevel = nullptr;
198 free (meta->longFstLevel);
199 meta->longFstLevel = nullptr;
200 }
201
202 /*****************************************************************************
203 * ParseTime() --
204 *
205 * Arthur Taylor / MDL
206 *
207 * PURPOSE
208 * To parse the time data from the grib2 integer array to a time_t in
209 * UTC seconds from the epoch.
210 *
211 * ARGUMENTS
212 * AnsTime = The time_t value to fill with the resulting time. (Output)
213 * year = The year to parse. (Input)
214 * mon = The month to parse. (Input)
215 * day = The day to parse. (Input)
216 * hour = The hour to parse. (Input)
217 * min = The minute to parse. (Input)
218 * sec = The second to parse. (Input)
219 *
220 * FILES/DATABASES: None
221 *
222 * RETURNS: int (could use errSprintf())
223 * 0 = OK
224 * -1 = gribLen is too small.
225 *
226 * HISTORY
227 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
228 * 4/2003 AAT: Modified to use year/mon/day/hour/min/sec instead of an
229 * integer array.
230 * 2/2004 AAT: Added error checks (because of corrupt GRIB1 files)
231 *
232 * NOTES
233 * 1) Couldn't use the default time_zone variable (concern over portability
234 * issues), so we print the hours, and compare them to the hours we had
235 * intended. Then subtract the difference from the AnsTime.
236 * 2) Need error check for times outside of 1902..2037.
237 *****************************************************************************
238 */
ParseTime(double * AnsTime,int year,uChar mon,uChar day,uChar hour,uChar min,uChar sec)239 int ParseTime (double *AnsTime, int year, uChar mon, uChar day, uChar hour,
240 uChar min, uChar sec)
241 {
242 /* struct tm time; *//* A temporary variable to put the time info into. */
243 /* char buffer[10]; *//* Used when printing the AnsTime's Hr. */
244 /* int timeZone; *//* The adjustment in Hr needed to get the right UTC * time. */
245
246 if ((year < 1900) || (year > 2100)) {
247 errSprintf ("ParseTime:: year %d is invalid\n", year);
248 /* return -1; */
249 year += 2000;
250 }
251 /* sec is allowed to be 61 for leap seconds. */
252 if ((mon > 12) || (day == 0) || (day > 31) || (hour > 24) || (min > 60) ||
253 (sec > 61)) {
254 errSprintf ("ParseTime:: Problems with %d/%d %d:%d:%d\n", mon, day,
255 hour, min, sec);
256 return -1;
257 }
258 Clock_ScanDate (AnsTime, year, mon, day);
259 *AnsTime += hour * 3600. + min * 60. + sec;
260 /* *AnsTime -= Clock_GetTimeZone() * 3600;*/
261
262 /*
263 memset (&time, 0, sizeof (struct tm));
264 time.tm_year = year - 1900;
265 time.tm_mon = mon - 1;
266 time.tm_mday = day;
267 time.tm_hour = hour;
268 time.tm_min = min;
269 time.tm_sec = sec;
270 printf ("%ld\n", mktime (&time));
271 *AnsTime = mktime (&time) - (Clock_GetTimeZone () * 3600);
272 */
273 /* Cheap method of getting global time_zone variable. */
274 /*
275 strftime (buffer, 10, "%H", gmtime (AnsTime));
276 timeZone = atoi (buffer) - hour;
277 if (timeZone < 0) {
278 timeZone += 24;
279 }
280 *AnsTime = *AnsTime - (timeZone * 3600);
281 */
282 return 0;
283 }
284
285 /*****************************************************************************
286 * ParseSect0() --
287 *
288 * Arthur Taylor / MDL
289 *
290 * PURPOSE
291 * To verify and parse section 0 data.
292 *
293 * ARGUMENTS
294 * is0 = The unpacked section 0 array. (Input)
295 * ns0 = The size of section 0. (Input)
296 * grib_len = The length of the entire grib message. (Input)
297 * meta = The structure to fill. (Output)
298 *
299 * FILES/DATABASES: None
300 *
301 * RETURNS: int (could use errSprintf())
302 * 0 = OK
303 * -1 = ns0 is too small.
304 * -2 = unexpected values in is0.
305 *
306 * HISTORY
307 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
308 *
309 * NOTES
310 * 1) 1196575042L == ASCII representation of "GRIB"
311 *****************************************************************************
312 */
ParseSect0(sInt4 * is0,sInt4 ns0,sInt4 grib_len,grib_MetaData * meta)313 static int ParseSect0 (sInt4 *is0, sInt4 ns0, sInt4 grib_len,
314 grib_MetaData *meta)
315 {
316 if (ns0 < 9) {
317 return -1;
318 }
319 if ((is0[0] != 1196575042L) || (is0[7] != 2) || (is0[8] != grib_len)) {
320 errSprintf ("ERROR IS0 has unexpected values: %ld %ld %ld\n",
321 is0[0], is0[7], is0[8]);
322 errSprintf ("Should be %ld %d %ld\n", 1196575042L, 2, grib_len);
323 return -2;
324 }
325 meta->pds2.prodType = (uChar) is0[6];
326 return 0;
327 }
328
329 /*****************************************************************************
330 * ParseSect1() --
331 *
332 * Arthur Taylor / MDL
333 *
334 * PURPOSE
335 * To verify and parse section 1 data.
336 *
337 * ARGUMENTS
338 * is1 = The unpacked section 1 array. (Input)
339 * ns1 = The size of section 1. (Input)
340 * meta = The structure to fill. (Output)
341 *
342 * FILES/DATABASES: None
343 *
344 * RETURNS: int (could use errSprintf())
345 * 0 = OK
346 * -1 = ns1 is too small.
347 * -2 = unexpected values in is1.
348 *
349 * HISTORY
350 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
351 *
352 * NOTES
353 *****************************************************************************
354 */
ParseSect1(sInt4 * is1,sInt4 ns1,grib_MetaData * meta)355 static int ParseSect1 (sInt4 *is1, sInt4 ns1, grib_MetaData *meta)
356 {
357 if (ns1 < 21) {
358 return -1;
359 }
360 if (is1[4] != 1) {
361 errSprintf ("ERROR IS1 not labeled correctly. %ld\n", is1[4]);
362 return -2;
363 }
364 meta->center = (unsigned short int) is1[5];
365 meta->subcenter = (unsigned short int) is1[7];
366 meta->pds2.mstrVersion = (uChar) is1[9];
367 meta->pds2.lclVersion = (uChar) is1[10];
368 if (((meta->pds2.mstrVersion < 1) || (meta->pds2.mstrVersion > 5)) ||
369 (meta->pds2.lclVersion > 1)) {
370 if (meta->pds2.mstrVersion == 0) {
371 printf ("Warning: Master table version == 0, was experimental\n"
372 "I don't have a copy, and don't know where to get one\n"
373 "Use meta data at your own risk.\n");
374 } else if (meta->pds2.mstrVersion != 255) {
375 printf ("Warning: use meta data at your own risk.\n");
376 printf ("Supported master table versions: (1,2,3,4,5) yours is %u... ",
377 meta->pds2.mstrVersion);
378 printf ("Supported local table version supported (0,1) yours is %u...\n",
379 meta->pds2.lclVersion);
380 }
381 }
382 meta->pds2.sigTime = (uChar) is1[11];
383 if (ParseTime (&(meta->pds2.refTime), is1[12], is1[14], is1[15], is1[16],
384 is1[17], is1[18]) != 0) {
385 preErrSprintf ("Error in call to ParseTime from ParseSect1 (GRIB2)");
386 return -2;
387 }
388 meta->pds2.operStatus = (uChar) is1[19];
389 meta->pds2.dataType = (uChar) is1[20];
390 return 0;
391 }
392
393 /*****************************************************************************
394 * ParseSect2_Wx() --
395 *
396 * Arthur Taylor / MDL
397 *
398 * PURPOSE
399 * To verify and parse section 2 data when we know the variable is of type
400 * Wx (Weather).
401 *
402 * ARGUMENTS
403 * rdat = The float data in section 2. (Input)
404 * nrdat = Length of rdat. (Input)
405 * idat = The integer data in section 2. (Input)
406 * nidat = Length of idat. (Input)
407 * Wx = The weather structure to fill. (Output)
408 * simpVer = The version of the simple weather code to use when parsing the
409 * WxString. (Input)
410 *
411 * FILES/DATABASES: None
412 *
413 * RETURNS: int (could use errSprintf())
414 * 0 = OK
415 * -1 = nrdat or nidat is too small.
416 * -2 = unexpected values in rdat.
417 *
418 * HISTORY
419 * 2/2003 Arthur Taylor (MDL/RSIS): Created.
420 * 5/2003 AAT: Stopped messing around with the way buffer and data[i]
421 * were allocated. It was confusing the free routine.
422 * 5/2003 AAT: Added maxLen to Wx structure.
423 * 6/2003 AAT: Revisited after Matt (matt@wunderground.com) informed me of
424 * memory problems.
425 * 1) I had a memory leak caused by a buffer+= buffLen
426 * 2) buffLen could have increased out of bounds of buffer.
427 * 8/2003 AAT: Found an invalid "assertion" when dealing with non-NULL
428 * terminated weather groups.
429 *
430 * NOTES
431 * 1) May want to rewrite so that we don't need 'meta->sect2NumGroups'
432 *****************************************************************************
433 */
ParseSect2_Wx(float * rdat,sInt4 nrdat,sInt4 * idat,uInt4 nidat,sect2_WxType * Wx,int simpVer)434 static int ParseSect2_Wx (float *rdat, sInt4 nrdat, sInt4 *idat,
435 uInt4 nidat, sect2_WxType *Wx, int simpVer)
436 {
437 size_t loc; /* Where we currently are in idat. */
438 size_t groupLen; /* Length of current group in idat. */
439 size_t j; /* Counter over the length of the current group. */
440 char *buffer; /* Used to store the current "ugly" string. */
441 int buffLen; /* Length of current "ugly" string. */
442 int len; /* length of current English phrases during creation
443 * of the maxEng[] data. */
444 int i; /* assists in traversing the maxEng[] array. */
445
446 if (nrdat < 1) {
447 return -1;
448 }
449
450 if (rdat[0] != 0) {
451 errSprintf ("ERROR: Expected rdat to be empty when dealing with "
452 "section 2 Weather data\n");
453 return -2;
454 }
455 Wx->dataLen = 0;
456 Wx->data = nullptr;
457 Wx->maxLen = 0;
458 for (i = 0; i < NUM_UGLY_WORD; i++) {
459 Wx->maxEng[i] = 0;
460 }
461
462 loc = 0;
463 if (nidat == 0) {
464 errSprintf ("ERROR: Ran out of idat data\n");
465 return -1;
466 }
467 groupLen = idat[loc++];
468
469 loc++; /* Skip the decimal scale factor data. */
470 /* Note: This also assures that buffLen stays <= nidat. */
471 if (loc + groupLen >= nidat) {
472 errSprintf ("ERROR: Ran out of idat data\n");
473 return -1;
474 }
475
476 buffLen = 0;
477 buffer = (char *) malloc ((nidat + 1) * sizeof (char));
478 while (groupLen > 0) {
479 for (j = 0; j < groupLen; j++) {
480 buffer[buffLen] = (char) idat[loc];
481 buffLen++;
482 loc++;
483 if (buffer[buffLen - 1] == '\0') {
484 Wx->dataLen++;
485 Wx->data = (char **) realloc ((void *) Wx->data,
486 Wx->dataLen * sizeof (char *));
487 /* This is done after the realloc, just to make sure we have
488 * enough memory allocated. */
489 /* Assert: buffLen is 1 more than strlen(buffer). */
490 Wx->data[Wx->dataLen - 1] = (char *)
491 malloc (buffLen * sizeof (char));
492 strcpy (Wx->data[Wx->dataLen - 1], buffer);
493 if (Wx->maxLen < buffLen) {
494 Wx->maxLen = buffLen;
495 }
496 buffLen = 0;
497 }
498 }
499 if (loc >= nidat) {
500 groupLen = 0;
501 } else {
502 groupLen = idat[loc];
503 loc++;
504 if (groupLen != 0) {
505 loc++; /* Skip the decimal scale factor data. */
506 /* Note: This also assures that buffLen stays <= nidat. */
507 if (loc + groupLen >= nidat) {
508 errSprintf ("ERROR: Ran out of idat data\n");
509 free (buffer);
510 return -1;
511 }
512 }
513 }
514 }
515 if (buffLen != 0) {
516 buffer[buffLen] = '\0';
517 Wx->dataLen++;
518 Wx->data = (char **) realloc ((void *) Wx->data,
519 Wx->dataLen * sizeof (char *));
520 /* Assert: buffLen is 1 more than strlen(buffer). -- FALSE -- */
521 buffLen = static_cast<int>(strlen (buffer)) + 1;
522
523 Wx->data[Wx->dataLen - 1] = (char *) malloc (buffLen * sizeof (char));
524 if (Wx->maxLen < buffLen) {
525 Wx->maxLen = buffLen;
526 }
527 strcpy (Wx->data[Wx->dataLen - 1], buffer);
528 }
529 free (buffer);
530 Wx->ugly = (UglyStringType *) malloc (Wx->dataLen *
531 sizeof (UglyStringType));
532 Wx->f_valid = (uChar *) malloc (Wx->dataLen * sizeof (uChar));
533 for (j = 0; j < Wx->dataLen; j++) {
534 if (ParseUglyString (&(Wx->ugly[j]), Wx->data[j], simpVer) == 0) {
535 Wx->f_valid[j] = 1;
536 } else {
537 Wx->f_valid[j] = 0;
538 }
539 }
540 /* We want to know how many bytes we need for each English phrase column,
541 * so we walk through each column calculating that value. */
542 for (i = 0; i < NUM_UGLY_WORD; i++) {
543 /* Assert: Already initialized Wx->maxEng[i]. */
544 for (j = 0; j < Wx->dataLen; j++) {
545 if (Wx->ugly[j].english[i] != nullptr) {
546 len = static_cast<int>(strlen (Wx->ugly[j].english[i]));
547 if (len > Wx->maxEng[i]) {
548 Wx->maxEng[i] = len;
549 }
550 }
551 }
552 }
553 return 0;
554 }
555
ParseSect2_Hazard(float * rdat,sInt4 nrdat,sInt4 * idat,uInt4 nidat,sect2_HazardType * Hazard,int simpWWA)556 static int ParseSect2_Hazard (float *rdat, sInt4 nrdat, sInt4 *idat,
557 uInt4 nidat, sect2_HazardType *Hazard, int simpWWA)
558 {
559 size_t loc; /* Where we currently are in idat. */
560 size_t groupLen; /* Length of current group in idat. */
561 size_t j; /* Counter over the length of the current group. */
562 int len; /* length of current english phrases during creation
563 * of the maxEng[] data. */
564 int i; /* assists in traversing the maxEng[] array. */
565 char *buffer; /* Used to store the current Hazard string. */
566 int buffLen; /* Length of current Hazard string. */
567 /*
568 int k;
569 */
570
571 if (nrdat < 1) {
572 return -1;
573 }
574
575 if (rdat[0] != 0) {
576 errSprintf ("ERROR: Expected rdat to be empty when dealing with "
577 "section 2 Weather data\n");
578 return -2;
579 }
580 Hazard->dataLen = 0;
581 Hazard->data = nullptr;
582 Hazard->maxLen = 0;
583 for (j = 0; j < NUM_HAZARD_WORD; j++) {
584 Hazard->maxEng[j] = 0;
585 }
586
587 loc = 0;
588 if (nidat == 0) {
589 errSprintf ("ERROR: Ran out of idat data\n");
590 return -1;
591 }
592 groupLen = idat[loc++];
593
594 loc++; /* Skip the decimal scale factor data. */
595 /* Note: This also assures that buffLen stays <= nidat. */
596 if (loc + groupLen >= nidat) {
597 errSprintf ("ERROR: Ran out of idat data\n");
598 return -1;
599 }
600
601 buffLen = 0;
602 buffer = (char *) malloc ((nidat + 1) * sizeof (char));
603 while (groupLen > 0) {
604 for (j = 0; j < groupLen; j++) {
605 buffer[buffLen] = (char) idat[loc];
606 buffLen++;
607 loc++;
608 if (buffer[buffLen - 1] == '\0') {
609 Hazard->dataLen++;
610 Hazard->data = (char **) realloc ((void *) Hazard->data,
611 Hazard->dataLen * sizeof (char *));
612 /* This is done after the realloc, just to make sure we have
613 * enough memory allocated. */
614 /* Assert: buffLen is 1 more than strlen(buffer). */
615 Hazard->data[Hazard->dataLen - 1] = (char *)
616 malloc (buffLen * sizeof (char));
617 strcpy (Hazard->data[Hazard->dataLen - 1], buffer);
618 if (Hazard->maxLen < buffLen) {
619 Hazard->maxLen = buffLen;
620 }
621 buffLen = 0;
622 }
623 }
624 if (loc >= nidat) {
625 groupLen = 0;
626 } else {
627 groupLen = idat[loc];
628 loc++;
629 if (groupLen != 0) {
630 loc++; /* Skip the decimal scale factor data. */
631 /* Note: This also assures that buffLen stays <= nidat. */
632 if (loc + groupLen >= nidat) {
633 errSprintf ("ERROR: Ran out of idat data\n");
634 free (buffer);
635 return -1;
636 }
637 }
638 }
639 }
640 if (buffLen != 0) {
641 buffer[buffLen] = '\0';
642 Hazard->dataLen++;
643 Hazard->data = (char **) realloc ((void *) Hazard->data,
644 Hazard->dataLen * sizeof (char *));
645 /* Assert: buffLen is 1 more than strlen(buffer). -- FALSE -- */
646 buffLen = static_cast<int>(strlen (buffer)) + 1;
647
648 Hazard->data[Hazard->dataLen - 1] = (char *) malloc (buffLen * sizeof (char));
649 if (Hazard->maxLen < buffLen) {
650 Hazard->maxLen = buffLen;
651 }
652 strcpy (Hazard->data[Hazard->dataLen - 1], buffer);
653 }
654 free (buffer);
655 Hazard->haz = (HazardStringType *) malloc (Hazard->dataLen *
656 sizeof (HazardStringType));
657 Hazard->f_valid = (uChar *) malloc (Hazard->dataLen * sizeof (uChar));
658 for (j = 0; j < Hazard->dataLen; j++) {
659 ParseHazardString (&(Hazard->haz[j]), Hazard->data[j], simpWWA);
660 Hazard->f_valid[j] = 1;
661 /*
662 printf ("%d : %d : %s", j, Hazard->haz[j].numValid, Hazard->data[j]);
663 for (k = 0; k < Hazard->haz[j].numValid; k++) {
664 printf (": %s", Hazard->haz[j].english[k]);
665 }
666 printf ("\n");
667 */
668 }
669 /* We want to know how many bytes we need for each english phrase column,
670 * so we walk through each column calculating that value. */
671 for (i = 0; i < NUM_HAZARD_WORD; i++) {
672 /* Assert: Already initialized Hazard->maxEng[i]. */
673 for (j = 0; j < Hazard->dataLen; j++) {
674 if (Hazard->haz[j].english[i] != nullptr) {
675 len = static_cast<int>(strlen (Hazard->haz[j].english[i]));
676 if (len > Hazard->maxEng[i]) {
677 Hazard->maxEng[i] = len;
678 }
679 }
680 }
681 }
682 return 0;
683 }
684
685 /*****************************************************************************
686 * ParseSect2_Unknown() --
687 *
688 * Arthur Taylor / MDL
689 *
690 * PURPOSE
691 * To verify and parse section 2 data when we don't know anything more
692 * about the data.
693 *
694 * ARGUMENTS
695 * rdat = The float data in section 2. (Input)
696 * nrdat = Length of rdat. (Input)
697 * idat = The integer data in section 2. (Input)
698 * nidat = Length of idat. (Input)
699 * meta = The structure to fill. (Output)
700 *
701 * FILES/DATABASES: None
702 *
703 * RETURNS: int (could use errSprintf())
704 * 0 = OK
705 * -1 = nrdat or nidat is too small.
706 * -2 = unexpected values in rdat.
707 *
708 * HISTORY
709 * 2/2003 Arthur Taylor (MDL/RSIS): Created.
710 *
711 * NOTES
712 * In the extremely improbable case that there is both idat data and rdat
713 * data, we process the rdat data first.
714 *****************************************************************************
715 */
ParseSect2_Unknown(float * rdat,sInt4 nrdat,sInt4 * idat,sInt4 nidat,grib_MetaData * meta)716 static int ParseSect2_Unknown (float *rdat, sInt4 nrdat, sInt4 *idat,
717 sInt4 nidat, grib_MetaData *meta)
718 {
719 /* Used for easier access to answer. */
720 int loc; /* Where we currently are in idat. */
721 int ansLoc; /* Where we are in the answer data structure. */
722 sInt4 groupLen; /* Length of current group in idat. */
723 int j; /* Counter over the length of the current group. */
724
725 meta->pds2.sect2.unknown.dataLen = 0;
726 meta->pds2.sect2.unknown.data = nullptr;
727 ansLoc = 0;
728
729 /* Work with rdat data. */
730 loc = 0;
731 if (nrdat <= loc) {
732 errSprintf ("ERROR: Ran out of rdat data\n");
733 return -1;
734 }
735 groupLen = (sInt4) rdat[loc++];
736 loc++; /* Skip the decimal scale factor data. */
737 if (nrdat <= loc + groupLen) {
738 errSprintf ("ERROR: Ran out of rdat data\n");
739 return -1;
740 }
741 while (groupLen > 0) {
742 meta->pds2.sect2.unknown.dataLen += groupLen;
743 meta->pds2.sect2.unknown.data = (double *)
744 realloc ((void *) meta->pds2.sect2.unknown.data,
745 meta->pds2.sect2.unknown.dataLen * sizeof (double));
746 for (j = 0; j < groupLen; j++) {
747 meta->pds2.sect2.unknown.data[ansLoc++] = rdat[loc++];
748 }
749 if (nrdat <= loc) {
750 groupLen = 0;
751 } else {
752 groupLen = (sInt4) rdat[loc++];
753 if (groupLen != 0) {
754 loc++; /* Skip the decimal scale factor data. */
755 if (nrdat <= loc + groupLen) {
756 errSprintf ("ERROR: Ran out of rdat data\n");
757 return -1;
758 }
759 }
760 }
761 }
762
763 /* Work with idat data. */
764 loc = 0;
765 if (nidat <= loc) {
766 errSprintf ("ERROR: Ran out of idat data\n");
767 return -1;
768 }
769 groupLen = idat[loc++];
770 loc++; /* Skip the decimal scale factor data. */
771 if (nidat <= loc + groupLen) {
772 errSprintf ("ERROR: Ran out of idat data\n");
773 return -1;
774 }
775 while (groupLen > 0) {
776 meta->pds2.sect2.unknown.dataLen += groupLen;
777 meta->pds2.sect2.unknown.data = (double *)
778 realloc ((void *) meta->pds2.sect2.unknown.data,
779 meta->pds2.sect2.unknown.dataLen * sizeof (double));
780 for (j = 0; j < groupLen; j++) {
781 meta->pds2.sect2.unknown.data[ansLoc++] = idat[loc++];
782 }
783 if (nidat <= loc) {
784 groupLen = 0;
785 } else {
786 groupLen = idat[loc++];
787 if (groupLen != 0) {
788 loc++; /* Skip the decimal scale factor data. */
789 if (nidat <= loc + groupLen) {
790 errSprintf ("ERROR: Ran out of idat data\n");
791 return -1;
792 }
793 }
794 }
795 }
796 return 0;
797 }
798
799 /*****************************************************************************
800 * ParseSect3() --
801 *
802 * Arthur Taylor / MDL
803 *
804 * PURPOSE
805 * To verify and parse section 3 data.
806 *
807 * ARGUMENTS
808 * is3 = The unpacked section 3 array. (Input)
809 * ns3 = The size of section 3. (Input)
810 * meta = The structure to fill. (Output)
811 *
812 * FILES/DATABASES: None
813 *
814 * RETURNS: int (could use errSprintf())
815 * 0 = OK
816 * -1 = ns3 is too small.
817 * -2 = unexpected values in is3.
818 * -3 = un-supported map Projection.
819 *
820 * HISTORY
821 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
822 * 9/2003 AAT: Adjusted Radius Earth case 1,6 to be based on:
823 * Y * 10^D = R
824 * Where Y = original value, D is scale factor, R is scale value.
825 * 1/2004 AAT: Adjusted Radius Earth case 6 to always be 6371.229 km
826 *
827 * NOTES
828 * Need to add support for GS3_ORTHOGRAPHIC = 90,
829 * GS3_EQUATOR_EQUIDIST = 110, GS3_AZIMUTH_RANGE = 120
830 *****************************************************************************
831 */
ParseSect3(sInt4 * is3,sInt4 ns3,grib_MetaData * meta)832 static int ParseSect3 (sInt4 *is3, sInt4 ns3, grib_MetaData *meta)
833 {
834 double unit; /* Used to convert from stored value to degrees
835 * lat/lon. See GRIB2 Regulation 92.1.6 */
836 sInt4 angle; /* For Lat/Lon, 92.1.6 may not hold, in which case,
837 * angle != 0, and unit = angle/subdivision. */
838 sInt4 subdivision; /* see angle explanation. */
839
840 if (ns3 < 14) {
841 return -1;
842 }
843 if (is3[4] != 3) {
844 errSprintf ("ERROR IS3 not labeled correctly. %ld\n", is3[4]);
845 return -2;
846 }
847 if (is3[5] != 0) {
848 errSprintf ("Can not handle 'Source of Grid Definition' = %ld\n",
849 is3[5]);
850 errSprintf ("Can only handle grids defined in Code table 3.1\n");
851 // return -3;
852 }
853 meta->gds.numPts = is3[6];
854 if ((is3[10] != 0) || (is3[11] != 0)) {
855 errSprintf ("Un-supported Map Projection.\n All Supported "
856 "projections have 0 bytes following the template.\n");
857 // return -3;
858 }
859 meta->gds.projType = (uChar) is3[12];
860
861 // Do not refuse to convert the GRIB file if only the projection is unknown.
862
863 /*
864 if ((is3[12] != GS3_LATLON) && (is3[12] != GS3_MERCATOR) &&
865 (is3[12] != GS3_POLAR) && (is3[12] != GS3_LAMBERT)) {
866 errSprintf ("Un-supported Map Projection %ld\n", is3[12]);
867 return -3;
868 }
869 */
870
871 /*
872 * Handle variables common to the supported templates.
873 */
874 if (ns3 < 38) {
875 return -1;
876 }
877 /* Assert: is3[14] is the shape of the earth. */
878 meta->gds.hdatum = 0;
879 switch (is3[14]) {
880 case 0:
881 meta->gds.f_sphere = 1;
882 meta->gds.majEarth = 6367.47;
883 meta->gds.minEarth = 6367.47;
884 break;
885 case 6:
886 meta->gds.f_sphere = 1;
887 meta->gds.majEarth = 6371.229;
888 meta->gds.minEarth = 6371.229;
889 break;
890 case 1:
891 meta->gds.f_sphere = 1;
892 /* Following assumes scale factor and scale value refer to
893 * scientific notation. */
894 /* Incorrect Assumption (9/8/2003): scale factor / value are based
895 * on: Y * 10^D = R, where Y = original value, D = scale factor, ___
896 * R = scale value. */
897
898 if ((is3[16] != GRIB2MISSING_s4) && (is3[15] != GRIB2MISSING_s1)) {
899 /* Assumes data is given in m (not km). */
900 double denom = pow (10.0, is3[15]) * 1000.;
901 if( denom == 0 )
902 {
903 errSprintf ("Invalid radius.\n");
904 return -2;
905 }
906 meta->gds.majEarth = is3[16] / denom;
907 meta->gds.minEarth = meta->gds.majEarth;
908 } else {
909 errSprintf ("Missing info on radius of Earth.\n");
910 return -2;
911 }
912 /* Check if our m assumption was valid. If it was not, they give us
913 * 6371 km, which we convert to 6.371 < 6.4 */
914 if (meta->gds.majEarth < 6.4) {
915 meta->gds.majEarth = meta->gds.majEarth * 1000.;
916 meta->gds.minEarth = meta->gds.minEarth * 1000.;
917 }
918 break;
919 case 2:
920 meta->gds.f_sphere = 0;
921 meta->gds.majEarth = 6378.160;
922 meta->gds.minEarth = 6356.775;
923 break;
924 case 4: // GRS80
925 meta->gds.f_sphere = 0;
926 meta->gds.majEarth = 6378.137;
927 meta->gds.minEarth = meta->gds.majEarth * (1 - 1 / 298.257222101);
928 break;
929 case 5: // WGS84
930 meta->gds.f_sphere = 0;
931 meta->gds.majEarth = 6378.137;
932 meta->gds.minEarth = meta->gds.majEarth * (1 - 1 / 298.257223563);
933 break;
934 case 3:
935 meta->gds.f_sphere = 0;
936 /* Following assumes scale factor and scale value refer to
937 * scientific notation. */
938 /* Incorrect Assumption (9/8/2003): scale factor / value are based
939 * on: Y * 10^D = R, where Y = original value, D = scale factor, ___
940 * R = scale value. */
941 if ((is3[21] != GRIB2MISSING_s4) && (is3[20] != GRIB2MISSING_s1) &&
942 (is3[26] != GRIB2MISSING_s4) && (is3[25] != GRIB2MISSING_s1)) {
943 /* Assumes data is given in km (not m). */
944 double denomMaj = pow (10.0, is3[20]);
945 double denomMin = pow (10.0, is3[25]);
946 if( denomMaj == 0.0 || denomMin == 0.0 )
947 {
948 errSprintf ("Invalid major / minor axis.\n");
949 return -2;
950 }
951 meta->gds.majEarth = is3[21] / denomMaj;
952 meta->gds.minEarth = is3[26] / denomMin;
953 } else {
954 errSprintf ("Missing info on major / minor axis of Earth.\n");
955 return -2;
956 }
957 /* Check if our km assumption was valid. If it was not, they give us
958 * 6371000 m, which is > 6400. */
959 if (meta->gds.majEarth > 6400) {
960 meta->gds.majEarth = meta->gds.majEarth / 1000.;
961 }
962 if (meta->gds.minEarth > 6400) {
963 meta->gds.minEarth = meta->gds.minEarth / 1000.;
964 }
965 break;
966 case 7:
967 meta->gds.f_sphere = 0;
968 /* Following assumes scale factor and scale value refer to
969 * scientific notation. */
970 /* Incorrect Assumption (9/8/2003): scale factor / value are based
971 * on: Y * 10^D = R, where Y = original value, D = scale factor, ___
972 * R = scale value. */
973 if ((is3[21] != GRIB2MISSING_s4) && (is3[20] != GRIB2MISSING_s1) &&
974 (is3[26] != GRIB2MISSING_s4) && (is3[25] != GRIB2MISSING_s1)) {
975 /* Assumes data is given in m (not km). */
976 double denomMaj = pow (10.0, is3[20]) * 1000.;
977 double denomMin = pow (10.0, is3[25]) * 1000.;
978 if( denomMaj == 0.0 || denomMin == 0.0 )
979 {
980 errSprintf ("Invalid major / minor axis.\n");
981 return -2;
982 }
983 meta->gds.majEarth = is3[21] / denomMaj;
984 meta->gds.minEarth = is3[26] / denomMin;
985 } else {
986 errSprintf ("Missing info on major / minor axis of Earth.\n");
987 return -2;
988 }
989 /* Check if our m assumption was valid. If it was not, they give us
990 * 6371 km, which we convert to 6.371 < 6.4 */
991 if (meta->gds.majEarth < 6.4) {
992 meta->gds.majEarth = meta->gds.majEarth * 1000.;
993 }
994 if (meta->gds.minEarth < 6.4) {
995 meta->gds.minEarth = meta->gds.minEarth * 1000.;
996 }
997 break;
998 case 8:
999 meta->gds.f_sphere = 1;
1000 meta->gds.majEarth = 6371.2;
1001 meta->gds.minEarth = 6371.2;
1002 meta->gds.hdatum = 1;
1003 break;
1004 default:
1005 errSprintf ("Undefined shape of earth? %ld\n", is3[14]);
1006 return -2;
1007 }
1008 /* Validate the radEarth is reasonable. */
1009 if ((meta->gds.majEarth > 6400) || (meta->gds.majEarth < 6300) ||
1010 (meta->gds.minEarth > 6400) || (meta->gds.minEarth < 6300)) {
1011 errSprintf ("Bad shape of earth? %f %f\n", meta->gds.majEarth,
1012 meta->gds.minEarth);
1013 return -2;
1014 }
1015 meta->gds.Nx = is3[30];
1016 meta->gds.Ny = is3[34];
1017 if ((meta->gds.Nx != 0 && meta->gds.Ny > UINT_MAX / meta->gds.Nx) ||
1018 meta->gds.Nx * meta->gds.Ny != meta->gds.numPts) {
1019 errSprintf ("Nx * Ny != number of points?\n");
1020 return -2;
1021 }
1022
1023 /* Initialize variables prior to parsing the specific templates. */
1024 unit = 1e-6;
1025 meta->gds.center = 0;
1026 meta->gds.scaleLat1 = meta->gds.scaleLat2 = 0;
1027 meta->gds.southLat = meta->gds.southLon = 0;
1028 meta->gds.lat2 = meta->gds.lon2 = 0;
1029 switch (is3[12]) {
1030 case GS3_LATLON: /* 0: Regular lat/lon grid. */
1031 case GS3_ROTATED_LATLON: // 1: Rotated lat/lon grid
1032 case GS3_GAUSSIAN_LATLON: /* 40: Gaussian lat/lon grid. */
1033 if (ns3 < 72) {
1034 return -1;
1035 }
1036 angle = is3[38];
1037 subdivision = is3[42];
1038 if (angle != 0) {
1039 if (subdivision == 0) {
1040 errSprintf ("subdivision of 0? Could not determine unit"
1041 " for latlon grid\n");
1042 return -2;
1043 }
1044 unit = angle / (double) (subdivision);
1045 }
1046 if ((is3[46] == GRIB2MISSING_s4) || (is3[50] == GRIB2MISSING_s4) ||
1047 (is3[55] == GRIB2MISSING_s4) || (is3[59] == GRIB2MISSING_s4) ||
1048 (is3[63] == GRIB2MISSING_s4) || (is3[67] == GRIB2MISSING_s4)) {
1049 errSprintf ("Lat/Lon grid is not defined completely.\n");
1050 return -2;
1051 }
1052 meta->gds.lat1 = is3[46] * unit;
1053 meta->gds.lon1 = is3[50] * unit;
1054 meta->gds.resFlag = (uChar) is3[54];
1055 meta->gds.lat2 = is3[55] * unit;
1056 meta->gds.lon2 = is3[59] * unit;
1057 meta->gds.Dx = is3[63] * unit; /* degrees. */
1058 if (is3[12] == GS3_GAUSSIAN_LATLON) {
1059 int np = is3[67]; /* parallels between a pole and the equator */
1060 if( np == 0 )
1061 {
1062 errSprintf ("Gaussian Lat/Lon grid is not defined completely.\n");
1063 return -2;
1064 }
1065 meta->gds.Dy = 90.0 / np;
1066 } else
1067 meta->gds.Dy = is3[67] * unit; /* degrees. */
1068 meta->gds.scan = (uChar) is3[71];
1069 meta->gds.meshLat = 0;
1070 meta->gds.orientLon = 0;
1071 if( is3[12] == GS3_ROTATED_LATLON ) {
1072 if( ns3 < 84 ) {
1073 return -1;
1074 }
1075 meta->gds.f_typeLatLon = 3;
1076 meta->gds.southLat = is3[73-1] * unit;
1077 meta->gds.southLon = is3[77-1] * unit;
1078 meta->gds.angleRotate = is3[81-1] * unit;
1079 }
1080 /* Resolve resolution flag(bit 3,4). Copy Dx,Dy as appropriate. */
1081 if ((meta->gds.resFlag & GRIB2BIT_3) &&
1082 (!(meta->gds.resFlag & GRIB2BIT_4))) {
1083 meta->gds.Dy = meta->gds.Dx;
1084 } else if ((!(meta->gds.resFlag & GRIB2BIT_3)) &&
1085 (meta->gds.resFlag & GRIB2BIT_4)) {
1086 meta->gds.Dx = meta->gds.Dy;
1087 }
1088 break;
1089 case GS3_MERCATOR: /* 10: Mercator grid. */
1090 if (ns3 < 72) {
1091 return -1;
1092 }
1093 if ((is3[38] == GRIB2MISSING_s4) || (is3[42] == GRIB2MISSING_s4) ||
1094 (is3[47] == GRIB2MISSING_s4) || (is3[51] == GRIB2MISSING_s4) ||
1095 (is3[55] == GRIB2MISSING_s4) || (is3[60] == GRIB2MISSING_s4)) {
1096 errSprintf ("Mercator grid is not defined completely.\n");
1097 return -2;
1098 }
1099 meta->gds.lat1 = is3[38] * unit;
1100 meta->gds.lon1 = is3[42] * unit;
1101 meta->gds.resFlag = (uChar) is3[46];
1102 meta->gds.meshLat = is3[47] * unit;
1103 meta->gds.lat2 = is3[51] * unit;
1104 meta->gds.lon2 = is3[55] * unit;
1105 meta->gds.scan = (uChar) is3[59];
1106 meta->gds.orientLon = is3[60] * unit;
1107 meta->gds.Dx = is3[64] / 1000.; /* mm -> m */
1108 meta->gds.Dy = is3[68] / 1000.; /* mm -> m */
1109 /* Resolve resolution flag(bit 3,4). Copy Dx,Dy as appropriate. */
1110 if ((meta->gds.resFlag & GRIB2BIT_3) &&
1111 (!(meta->gds.resFlag & GRIB2BIT_4))) {
1112 if (is3[64] == GRIB2MISSING_s4) {
1113 errSprintf ("Mercator grid is not defined completely.\n");
1114 return -2;
1115 }
1116 meta->gds.Dy = meta->gds.Dx;
1117 } else if ((!(meta->gds.resFlag & GRIB2BIT_3)) &&
1118 (meta->gds.resFlag & GRIB2BIT_4)) {
1119 if (is3[68] == GRIB2MISSING_s4) {
1120 errSprintf ("Mercator grid is not defined completely.\n");
1121 return -2;
1122 }
1123 meta->gds.Dx = meta->gds.Dy;
1124 }
1125 break;
1126 case GS3_TRANSVERSE_MERCATOR: /* 12: Transverse mercator */
1127 if (ns3 < 84) {
1128 return -1;
1129 }
1130 meta->gds.latitude_of_origin = is3[38] * unit;
1131 meta->gds.central_meridian = is3[42] * unit;
1132 meta->gds.resFlag = (uChar) is3[46];
1133 {
1134 float fTemp;
1135 GUInt32 nTemp = is3[47] < 0 ? (-is3[47]) | 0x80000000U : is3[47];
1136 memcpy(&fTemp, &nTemp, 4);
1137 meta->gds.scaleLat1 = fTemp;
1138 }
1139 meta->gds.x0 = is3[51] / 100.0;
1140 meta->gds.y0 = is3[55] / 100.0;
1141 meta->gds.scan = (uChar) is3[59];
1142 meta->gds.Dx = is3[60] / 100.0;
1143 meta->gds.Dy = is3[64] / 100.0;
1144 meta->gds.x1 = is3[68] / 100.0;
1145 meta->gds.y1 = is3[72] / 100.0;
1146 meta->gds.x2 = is3[76] / 100.0;
1147 meta->gds.y2 = is3[80] / 100.0;
1148 break;
1149
1150 case GS3_POLAR: /* 20: Polar Stereographic grid. */
1151 if (ns3 < 65) {
1152 return -1;
1153 }
1154 if ((is3[38] == GRIB2MISSING_s4) || (is3[42] == GRIB2MISSING_s4) ||
1155 (is3[47] == GRIB2MISSING_s4) || (is3[51] == GRIB2MISSING_s4)) {
1156 errSprintf ("Polar Stereographic grid is not defined "
1157 "completely.\n");
1158 return -2;
1159 }
1160 meta->gds.lat1 = is3[38] * unit;
1161 meta->gds.lon1 = is3[42] * unit;
1162 meta->gds.resFlag = (uChar) is3[46];
1163 /* Note (1) resFlag (bit 3,4) not applicable. */
1164 meta->gds.meshLat = is3[47] * unit;
1165 meta->gds.orientLon = is3[51] * unit;
1166 meta->gds.Dx = is3[55] / 1000.; /* mm -> m */
1167 meta->gds.Dy = is3[59] / 1000.; /* mm -> m */
1168 meta->gds.center = (uChar) is3[63];
1169 if (meta->gds.center & GRIB2BIT_1) {
1170 /* South polar stereographic. */
1171 meta->gds.scaleLat1 = meta->gds.scaleLat2 = -90;
1172 } else {
1173 /* North polar stereographic. */
1174 meta->gds.scaleLat1 = meta->gds.scaleLat2 = 90;
1175 }
1176 if (meta->gds.center & GRIB2BIT_2) {
1177 errSprintf ("Note (4) specifies no 'bi-polar stereograhic"
1178 " projections'.\n");
1179 return -2;
1180 }
1181 meta->gds.scan = (uChar) is3[64];
1182 break;
1183 case GS3_LAMBERT: /* 30: Lambert Conformal grid. */
1184 case GS3_ALBERS_EQUAL_AREA: /* 31: Albers equal area */
1185 if (ns3 < 81) {
1186 return -1;
1187 }
1188 if ((is3[38] == GRIB2MISSING_s4) || (is3[42] == GRIB2MISSING_s4) ||
1189 (is3[47] == GRIB2MISSING_s4) || (is3[51] == GRIB2MISSING_s4) ||
1190 (is3[65] == GRIB2MISSING_s4) || (is3[69] == GRIB2MISSING_s4)) {
1191 if( is3[12] == GS3_LAMBERT )
1192 {
1193 errSprintf ("Lambert Conformal grid is not defined "
1194 "completely.\n");
1195 }
1196 else
1197 {
1198 errSprintf ("Albers Equal Area grid is not defined "
1199 "completely.\n");
1200 }
1201 return -2;
1202 }
1203 meta->gds.lat1 = is3[38] * unit;
1204 meta->gds.lon1 = is3[42] * unit;
1205 meta->gds.resFlag = (uChar) is3[46];
1206 /* Note (3) resFlag (bit 3,4) not applicable. */
1207 meta->gds.meshLat = is3[47] * unit;
1208 meta->gds.orientLon = is3[51] * unit;
1209 meta->gds.Dx = is3[55] / 1000.; /* mm -> m */
1210 meta->gds.Dy = is3[59] / 1000.; /* mm -> m */
1211 meta->gds.center = (uChar) is3[63];
1212 meta->gds.scan = (uChar) is3[64];
1213 meta->gds.scaleLat1 = is3[65] * unit;
1214 meta->gds.scaleLat2 = is3[69] * unit;
1215 if( (is3[73] == GRIB2MISSING_s4) || (is3[77] == GRIB2MISSING_s4) )
1216 {
1217 meta->gds.southLat = 0.0;
1218 meta->gds.southLon = 0.0;
1219 }
1220 else
1221 {
1222 meta->gds.southLat = is3[73] * unit;
1223 meta->gds.southLon = is3[77] * unit;
1224 }
1225 break;
1226 case GS3_ORTHOGRAPHIC: /* 90: Orthographic grid. */
1227 // Misusing gdsType elements (gdsType needs extension)
1228 meta->gds.lat1 = is3[38];
1229 meta->gds.lon1 = is3[42];
1230 meta->gds.resFlag = (uChar) is3[46];
1231 meta->gds.Dx = is3[47];
1232 meta->gds.Dy = is3[51];
1233
1234 meta->gds.lon2 = is3[55] / 1000.; /* xp - X-coordinateSub-satellite, mm -> m */
1235 meta->gds.lat2 = is3[59] / 1000.; /* yp - Y-coordinateSub-satellite, mm -> m */
1236 meta->gds.scan = (uChar) is3[63];
1237 meta->gds.orientLon = is3[64]; /* angle */
1238 meta->gds.stretchFactor = is3[68] * 1000000.; /* altitude */
1239
1240 meta->gds.southLon = is3[72]; /* x0 - X-coordinateOrigin */
1241 meta->gds.southLat = is3[76]; /* y0 - Y-coordinateOrigin */
1242 break;
1243 case GS3_LAMBERT_AZIMUTHAL: /* 140: Lambert Azimuthal Equal Area Projection */
1244 meta->gds.lat1 = is3[38] * unit;
1245 meta->gds.lon1 = is3[42] * unit;
1246 meta->gds.meshLat = is3[46] * unit;
1247 meta->gds.orientLon = is3[50] * unit;
1248 meta->gds.resFlag = (uChar) is3[54];
1249 meta->gds.Dx = is3[55] / 1000.; /* mm -> m */
1250 meta->gds.Dy = is3[59] / 1000.; /* mm -> m */
1251 meta->gds.scan = (uChar) is3[63];
1252 break;
1253 default:
1254 errSprintf ("Un-supported Map Projection. %ld\n", is3[12]);
1255 // Don't abandon the conversion only because of an unknown projection
1256 break;
1257 //return -3;
1258 }
1259 if (meta->gds.scan != GRIB2BIT_2) {
1260 #ifdef DEBUG
1261 printf ("Scan mode is expected to be 0100 (i.e. %d) not %u\n",
1262 GRIB2BIT_2, meta->gds.scan);
1263 printf ("The merged GRIB2 Library should return it in 0100\n");
1264 printf ("The merged library swaps both NCEP and MDL data to scan "
1265 "mode 0100\n");
1266 #endif
1267 /*
1268 errSprintf ("Scan mode is expected to be 0100 (i.e. %d) not %d",
1269 GRIB2BIT_2, meta->gds.scan);
1270 return -2;
1271 */
1272 }
1273 return 0;
1274 }
1275
1276 /*****************************************************************************
1277 * ParseSect4Time2secV1() --
1278 *
1279 * Arthur Taylor / MDL
1280 *
1281 * PURPOSE
1282 * Attempt to parse time data in units provided by GRIB1 table 4, to
1283 * seconds.
1284 *
1285 * ARGUMENTS
1286 * time = The delta time to convert. (Input)
1287 * unit = The unit to convert. (Input)
1288 * ans = The converted answer. (Output)
1289 *
1290 * FILES/DATABASES: None
1291 *
1292 * RETURNS: int
1293 * 0 = OK
1294 * -1 = could not determine.
1295 *
1296 * HISTORY
1297 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
1298 * 1/2005 AAT: Fixed unit2sec[] table to have element 10 be 10800 (3 hours)
1299 * instead of 0.
1300 *
1301 * NOTES
1302 * http://www.nco.ncep.noaa.gov/pmb/docs/on388/table4.html
1303 *****************************************************************************
1304 */
ParseSect4Time2secV1(sInt4 time,int unit,double * ans)1305 int ParseSect4Time2secV1 (sInt4 time, int unit, double *ans)
1306 {
1307 /* Following is a lookup table for unit conversion (see code table 4.4). */
1308 static const sInt4 unit2sec[] = {
1309 60, 3600, 86400L, 0, 0,
1310 0, 0, 0, 0, 0,
1311 10800, 21600L, 43200L
1312 };
1313 if ((unit >= 0) && (unit < 13)) {
1314 if (unit2sec[unit] != 0) {
1315 *ans = (double) (time) * unit2sec[unit];
1316 return 0;
1317 }
1318 } else if (unit == 254) {
1319 *ans = (double) (time);
1320 return 0;
1321 }
1322 *ans = 0;
1323 return -1;
1324 }
1325
1326 /*****************************************************************************
1327 * ParseSect4Time2sec() --
1328 *
1329 * Arthur Taylor / MDL
1330 *
1331 * PURPOSE
1332 * Attempt to parse time data in units provided by GRIB2 table 4.4, to
1333 * seconds.
1334 *
1335 * ARGUMENTS
1336 * refTime = To add "years / centuries / decades and normals", we need a
1337 * reference time.
1338 * delt = The delta time to convert. (Input)
1339 * unit = The unit to convert. (Input)
1340 * ans = The converted answer. (Output)
1341 *
1342 * FILES/DATABASES: None
1343 *
1344 * RETURNS: int
1345 * 0 = OK
1346 * -1 = could not determine.
1347 *
1348 * HISTORY
1349 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
1350 * 1/2005 AAT: Fixed unit2sec[] table to have element 10 be 10800 (3 hours)
1351 * instead of 0.
1352 *
1353 * NOTES
1354 *****************************************************************************
1355 */
ParseSect4Time2sec(double refTime,sInt4 delt,int unit,double * ans)1356 int ParseSect4Time2sec (double refTime, sInt4 delt, int unit, double *ans)
1357 {
1358 /* Following is a lookup table for unit conversion (see code table 4.4). */
1359 static const sInt4 unit2sec[] = {
1360 60, 3600, 86400L, 0, 0,
1361 0, 0, 0, 0, 0,
1362 10800, 21600L, 43200L, 1
1363 };
1364 if ((unit >= 0) && (unit < 14)) {
1365 if (unit2sec[unit] != 0) {
1366 *ans = (double) (delt) * unit2sec[unit];
1367 return 0;
1368 } else {
1369 /* The procedure returns number of seconds to adjust by, rather
1370 * than the new time, which is why we subtract refTime */
1371 switch (unit) {
1372 case 3: /* month */
1373 *ans = Clock_AddMonthYear (refTime, delt, 0) - refTime;
1374 return 0;
1375 case 4: /* year */
1376 *ans = Clock_AddMonthYear (refTime, 0, delt) - refTime;
1377 return 0;
1378 case 5: /* decade */
1379 if( delt < INT_MIN / 10 || delt > INT_MAX / 10 )
1380 return -1;
1381 *ans = Clock_AddMonthYear (refTime, 0, delt * 10) - refTime;
1382 return 0;
1383 case 6: /* normal (30 year) */
1384 if( delt < INT_MIN / 30 || delt > INT_MAX / 30 )
1385 return -1;
1386 *ans = Clock_AddMonthYear (refTime, 0, delt * 30) - refTime;
1387 return 0;
1388 case 7: /* century (100 year) */
1389 if( delt < INT_MIN / 100 || delt > INT_MAX / 100 )
1390 return -1;
1391 *ans = Clock_AddMonthYear (refTime, 0, delt * 100) - refTime;
1392 return 0;
1393 }
1394 }
1395 }
1396 *ans = 0;
1397 return -1;
1398 }
1399
1400 /*****************************************************************************
1401 * sbit_2Comp_fourByte() -- Arthur Taylor / MDL
1402 *
1403 * PURPOSE
1404 * The NCEP g2clib-1.0.2 library stored the lower limits and upper limits
1405 * of probabilities using unsigned ints, whereas version 1.0.4 used signed
1406 * ints. The reason for the change is because some thresholds were negative.
1407 * To encode a negative value using an unsigned int, 1.0.2 used "2's
1408 * complement + 1". To encode a negative value using signed an int, 1.0.4
1409 * used a "sign bit". Example -2 => FFFFFFFE (1.0.2) => 80000002 (1.0.4).
1410 * The problem (for backward compatibility sake) is to be able to read both
1411 * encodings and get -2. If one only read the new encoding method, then
1412 * archived data would not be handled.
1413 * The algorithm is: If the number is positive or missing, leave it alone.
1414 * If the number is negative, look at the 2's complement method, and the sign
1415 * bit method, and use the method which results in a smaller absolute value.
1416 *
1417 * ARGUMENTS
1418 * data = The number read by NCEP's library. (Input)
1419 *
1420 * RETURNS: sInt4
1421 * The value of treating the number as read by either method
1422 *
1423 * HISTORY
1424 * 10/2007 Arthur Taylor (MDL): Created.
1425 *
1426 * NOTES
1427 * 1) This algorithm will impact the possible range of values, by reducing it
1428 * from -2^31..(2^31-1) to -2^30..(2^31-1).
1429 * 2) The NCEP change also impacted large positive values. One originally
1430 * could encode 0..2^32-1. Some confusion could arrise if the value was
1431 * originally encoded by 1.0.2 was in the range of 2^31..2^32-1.
1432 ****************************************************************************/
sbit_2Comp_fourByte(sInt4 data)1433 sInt4 sbit_2Comp_fourByte(sInt4 data)
1434 {
1435 sInt4 x; /* The pos. 2's complement interpretation of data */
1436 sInt4 y; /* The pos. sign bit interpretation of data */
1437
1438 if ((data == GRIB2MISSING_s4) || (data >= 0)) {
1439 return data;
1440 }
1441 if( data == INT_MIN ) // doesn't make sense since it is negative 0 in sign bit logic
1442 return 0;
1443 x = ~data + 1;
1444 y = data & 0x7fffffff;
1445 if (x < y) {
1446 return -1 * x;
1447 } else {
1448 return -1 * y;
1449 }
1450 }
1451
1452 /*****************************************************************************
1453 * sbit_2Comp_oneByte() -- Arthur Taylor / MDL
1454 *
1455 * PURPOSE
1456 * The NCEP g2clib-1.0.2 library stored the lower limits and upper limits
1457 * of probabilities using unsigned ints, whereas version 1.0.4 used signed
1458 * ints. The reason for the change is because some thresholds were negative.
1459 * To encode a negative value using an unsigned int, 1.0.2 used "2's
1460 * complement + 1". To encode a negative value using signed an int, 1.0.4
1461 * used a "sign bit". Example -2 => 11111110 (1.0.2) => 10000010 (1.0.4).
1462 * The problem (for backward compatibility sake) is to be able to read both
1463 * encodings and get -2. If one only read the new encoding method, then
1464 * archived data would not be handled.
1465 * The algorithm is: If the number is positive or missing, leave it alone.
1466 * If the number is negative, look at the 2's complement method, and the sign
1467 * bit method, and use the method which results in a smaller absolute value.
1468 *
1469 * ARGUMENTS
1470 * data = The number read by NCEP's library. (Input)
1471 *
1472 * RETURNS: sChar
1473 * The value of treating the number as read by either method
1474 *
1475 * HISTORY
1476 * 10/2007 Arthur Taylor (MDL): Created.
1477 *
1478 * NOTES
1479 * 1) This algorithm will impact the possible range of values, by reducing it
1480 * from -128..127 to -64...127.
1481 * 2) The NCEP change also impacted large positive values. One originally
1482 * could encode 0..255. Some confusion could arrise if the value was
1483 * originally encoded by 1.0.2 was in the range of 128..255.
1484 ****************************************************************************/
sbit_2Comp_oneByte(sChar data)1485 sChar sbit_2Comp_oneByte(sChar data)
1486 {
1487 sChar x; /* The pos. 2's complement interpretation of data */
1488 sChar y; /* The pos. sign bit interpretation of data */
1489
1490 if ((data == GRIB2MISSING_s1) || (data >= 0)) {
1491 return data;
1492 }
1493 x = ~data + 1;
1494 y = data & 0x7f;
1495 if (x < y) {
1496 return -1 * x;
1497 } else {
1498 return -1 * y;
1499 }
1500 }
1501
1502 /*****************************************************************************
1503 * ParseSect4() --
1504 *
1505 * Arthur Taylor / MDL
1506 *
1507 * PURPOSE
1508 * To verify and parse section 4 data.
1509 *
1510 * ARGUMENTS
1511 * is4 = The unpacked section 4 array. (Input)
1512 * ns4 = The size of section 4. (Input)
1513 * meta = The structure to fill. (Output)
1514 *
1515 * FILES/DATABASES: None
1516 *
1517 * RETURNS: int (could use errSprintf())
1518 * 0 = OK
1519 * -1 = ns4 is too small.
1520 * -2 = unexpected values in is4.
1521 * -4 = un-supported Sect 4 template.
1522 * -5 = unsupported forecast time unit.
1523 * -6 = Ran out of memory.
1524 *
1525 * HISTORY
1526 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
1527 * 3/2003 AAT: Added support for GS4_SATELLITE.
1528 * 3/2003 AAT: Adjusted allocing of sect4.Interval (should be safer now).
1529 * 9/2003 AAT: Adjusted interpretation of scale factor / value.
1530 * 5/2004 AAT: Added some memory checks.
1531 * 3/2005 AAT: Added a cast to (uChar) when comparing to GRIB2MISSING_1
1532 * 3/2005 AAT: Added GS4_PROBABIL_PNT.
1533 *
1534 * NOTES
1535 * Need to add support for GS4_RADAR = 20
1536 *****************************************************************************
1537 */
ParseSect4(sInt4 * is4,sInt4 ns4,grib_MetaData * meta)1538 static int ParseSect4 (sInt4 *is4, sInt4 ns4, grib_MetaData *meta)
1539 {
1540 int i; /* Counter for time intervals in template 4.8, 4.9
1541 * (typically 1) or counter for satellite band in
1542 * template 4.30. */
1543 void *temp_ptr; /* A temporary pointer when reallocating memory. */
1544 char *msg; /* A pointer to the current error message. */
1545
1546 if (ns4 < 9) {
1547 return -1;
1548 }
1549 if (is4[4] != 4) {
1550 #ifdef DEBUG
1551 printf ("ERROR IS4 not labeled correctly. %d\n", is4[4]);
1552 #endif
1553 errSprintf ("ERROR IS4 not labeled correctly. %d\n", is4[4]);
1554 return -2;
1555 }
1556
1557 if ((is4[7] != GS4_ANALYSIS) && (is4[7] != GS4_ENSEMBLE) &&
1558 (is4[7] != GS4_DERIVED) && (is4[7] != GS4_PROBABIL_PNT) &&
1559 (is4[7] != GS4_PERCENT_PNT) && (is4[7] != GS4_ERROR) &&
1560 (is4[7] != GS4_STATISTIC) && (is4[7] != GS4_PROBABIL_TIME) &&
1561 (is4[7] != GS4_PERCENT_TIME) && (is4[7] != GS4_ENSEMBLE_STAT) &&
1562 (is4[7] != GS4_SATELLITE) && (is4[7] != GS4_SATELLITE_SYNTHETIC) &&
1563 (is4[7] != GS4_DERIVED_INTERVAL) && (is4[7] != GS4_STATISTIC_SPATIAL_AREA) &&
1564 (is4[7] != GS4_ANALYSIS_CHEMICAL) && (is4[7] != GS4_OPTICAL_PROPERTIES_AEROSOL)) {
1565 #ifdef DEBUG
1566 //printf ("Un-supported Template. %d\n", is4[7]);
1567 #endif
1568 errSprintf ("Un-supported Template. %d\n", is4[7]);
1569 return -4;
1570 }
1571 meta->pds2.sect4.templat = (unsigned short int) is4[7];
1572
1573 /*
1574 * Handle variables common to the supported templates.
1575 */
1576 if (ns4 < 34) {
1577 return -1;
1578 }
1579 meta->pds2.sect4.cat = (uChar) is4[9];
1580 meta->pds2.sect4.subcat = (uChar) is4[10];
1581 int nOffset = 0;
1582 if( is4[7] == GS4_ANALYSIS_CHEMICAL ) {
1583 nOffset = 16 - 14;
1584 }
1585 else if( is4[7] == GS4_OPTICAL_PROPERTIES_AEROSOL ) {
1586 nOffset = 38 - 14;
1587 }
1588 if (ns4 < 34 + nOffset) {
1589 return -1;
1590 }
1591 meta->pds2.sect4.genProcess = (uChar) is4[11 + nOffset];
1592
1593 /* Initialize variables prior to parsing the specific templates. */
1594 meta->pds2.sect4.typeEnsemble = 0;
1595 meta->pds2.sect4.perturbNum = 0;
1596 meta->pds2.sect4.numberFcsts = 0;
1597 meta->pds2.sect4.derivedFcst = (uChar)-1;
1598 meta->pds2.sect4.validTime = meta->pds2.refTime;
1599
1600 if (meta->pds2.sect4.templat == GS4_SATELLITE) {
1601 meta->pds2.sect4.genID = (uChar) is4[12];
1602 meta->pds2.sect4.numBands = (uChar) is4[13];
1603 meta->pds2.sect4.bands =
1604 (sect4_BandType *) realloc ((void *) meta->pds2.sect4.bands,
1605 meta->pds2.sect4.numBands *
1606 sizeof (sect4_BandType));
1607 for (i = 0; i < meta->pds2.sect4.numBands; i++) {
1608 if (ns4 < 20 + 10 * i + 1) {
1609 return -1;
1610 }
1611 meta->pds2.sect4.bands[i].series =
1612 (unsigned short int) is4[14 + 10 * i];
1613 meta->pds2.sect4.bands[i].numbers =
1614 (unsigned short int) is4[16 + 10 * i];
1615 meta->pds2.sect4.bands[i].instType = (uChar) is4[18 + 10 * i];
1616 meta->pds2.sect4.bands[i].centWaveNum.factor =
1617 (uChar) is4[19 + 10 * i];
1618 meta->pds2.sect4.bands[i].centWaveNum.value = is4[20 + 10 * i];
1619 }
1620
1621 meta->pds2.sect4.fstSurfType = GRIB2MISSING_u1;
1622 meta->pds2.sect4.fstSurfScale = GRIB2MISSING_s1;
1623 meta->pds2.sect4.fstSurfValue = 0;
1624 meta->pds2.sect4.sndSurfType = GRIB2MISSING_u1;
1625 meta->pds2.sect4.sndSurfScale = GRIB2MISSING_s1;
1626 meta->pds2.sect4.sndSurfValue = 0;
1627
1628 return 0;
1629 }
1630 if (meta->pds2.sect4.templat == GS4_SATELLITE_SYNTHETIC) {
1631 meta->pds2.sect4.genID = (uChar) is4[12];
1632 meta->pds2.sect4.numBands = (uChar) is4[22];
1633 meta->pds2.sect4.bands =
1634 (sect4_BandType *) realloc ((void *) meta->pds2.sect4.bands,
1635 meta->pds2.sect4.numBands *
1636 sizeof (sect4_BandType));
1637 for (i = 0; i < meta->pds2.sect4.numBands; i++) {
1638 if (ns4 < 30 + 11 * i + 1) {
1639 return -1;
1640 }
1641 meta->pds2.sect4.bands[i].series =
1642 (unsigned short int) is4[23 + 11 * i];
1643 meta->pds2.sect4.bands[i].numbers =
1644 (unsigned short int) is4[25 + 11 * i];
1645 meta->pds2.sect4.bands[i].instType = (uChar) is4[27 + 11 * i];
1646 meta->pds2.sect4.bands[i].centWaveNum.factor =
1647 (uChar) is4[29 + 11 * i];
1648 meta->pds2.sect4.bands[i].centWaveNum.value = is4[30 + 11 * i];
1649 }
1650
1651 meta->pds2.sect4.fstSurfType = GRIB2MISSING_u1;
1652 meta->pds2.sect4.fstSurfScale = GRIB2MISSING_s1;
1653 meta->pds2.sect4.fstSurfValue = 0;
1654 meta->pds2.sect4.sndSurfType = GRIB2MISSING_u1;
1655 meta->pds2.sect4.sndSurfScale = GRIB2MISSING_s1;
1656 meta->pds2.sect4.sndSurfValue = 0;
1657
1658 return 0;
1659 }
1660
1661 meta->pds2.sect4.bgGenID = (uChar) is4[12 + nOffset];
1662 meta->pds2.sect4.genID = (uChar) is4[13 + nOffset];
1663 if ((is4[14 + nOffset] == GRIB2MISSING_u2) || (is4[16 + nOffset] == GRIB2MISSING_u1)) {
1664 meta->pds2.sect4.f_validCutOff = 0;
1665 meta->pds2.sect4.cutOff = 0;
1666 } else {
1667 meta->pds2.sect4.f_validCutOff = 1;
1668 meta->pds2.sect4.cutOff = is4[14 + nOffset] * 3600 + is4[16 + nOffset] * 60;
1669 }
1670 if (is4[18 + nOffset] == GRIB2MISSING_s4) {
1671 errSprintf ("Missing 'forecast' time?\n");
1672 return -5;
1673 }
1674 meta->pds2.sect4.foreUnit = is4[17 + nOffset];
1675 if (ParseSect4Time2sec (meta->pds2.refTime, is4[18 + nOffset], is4[17 + nOffset],
1676 &(meta->pds2.sect4.foreSec)) != 0) {
1677 errSprintf ("Unable to convert this TimeUnit: %ld\n", is4[17 + nOffset]);
1678 return -5;
1679 }
1680
1681 meta->pds2.sect4.validTime = (time_t) (meta->pds2.refTime +
1682 meta->pds2.sect4.foreSec);
1683
1684 /*
1685 * Following is based on what was needed to get correct Radius of Earth in
1686 * section 3. (Hopefully they are consistent).
1687 */
1688 meta->pds2.sect4.fstSurfType = (uChar) is4[22 + nOffset];
1689 if ((is4[24 + nOffset] == GRIB2MISSING_s4) || (is4[23 + nOffset] == GRIB2MISSING_s1) ||
1690 (meta->pds2.sect4.fstSurfType == GRIB2MISSING_u1)) {
1691 meta->pds2.sect4.fstSurfScale = GRIB2MISSING_s1;
1692 meta->pds2.sect4.fstSurfValue = 0;
1693 } else {
1694 meta->pds2.sect4.fstSurfScale = is4[23 + nOffset];
1695 meta->pds2.sect4.fstSurfValue = is4[24 + nOffset] / pow (10.0, is4[23 + nOffset]);
1696 }
1697 meta->pds2.sect4.sndSurfType = (uChar) is4[28 + nOffset];
1698 if ((is4[30 + nOffset] == GRIB2MISSING_s4) || (is4[29 + nOffset] == GRIB2MISSING_s1) ||
1699 (meta->pds2.sect4.sndSurfType == GRIB2MISSING_u1)) {
1700 meta->pds2.sect4.sndSurfScale = GRIB2MISSING_s1;
1701 meta->pds2.sect4.sndSurfValue = 0;
1702 } else {
1703 meta->pds2.sect4.sndSurfScale = is4[29 + nOffset];
1704 meta->pds2.sect4.sndSurfValue = is4[30 + nOffset] / pow (10.0, is4[29 + nOffset]);
1705 }
1706 switch (meta->pds2.sect4.templat) {
1707 case GS4_ANALYSIS: /* 4.0 */
1708 case GS4_ERROR: /* 4.7 */
1709 break;
1710 case GS4_ENSEMBLE: /* 4.1 */
1711 if (ns4 < 37) {
1712 return -1;
1713 }
1714 meta->pds2.sect4.typeEnsemble = (uChar) is4[34];
1715 meta->pds2.sect4.perturbNum = (uChar) is4[35];
1716 meta->pds2.sect4.numberFcsts = (uChar) is4[36];
1717 break;
1718 case GS4_ENSEMBLE_STAT: /* 4.11 */
1719 if (ns4 < 46) {
1720 return -1;
1721 }
1722 meta->pds2.sect4.typeEnsemble = (uChar) is4[34];
1723 meta->pds2.sect4.perturbNum = (uChar) is4[35];
1724 meta->pds2.sect4.numberFcsts = (uChar) is4[36];
1725 if (ParseTime (&(meta->pds2.sect4.validTime), is4[37], is4[39],
1726 is4[40], is4[41], is4[42], is4[43]) != 0) {
1727 msg = errSprintf (nullptr);
1728 meta->pds2.sect4.numInterval = (uChar) is4[44];
1729 if (meta->pds2.sect4.numInterval != 1) {
1730 errSprintf ("ERROR: in call to ParseTime from ParseSect4\n%s",
1731 msg);
1732 errSprintf ("Most likely they didn't complete bytes 38-44 of "
1733 "Template 4.11\n");
1734 free (msg);
1735 meta->pds2.sect4.numInterval = 0;
1736 return -1;
1737 }
1738 printf ("Warning: in call to ParseTime from ParseSect4\n%s", msg);
1739 free (msg);
1740 meta->pds2.sect4.validTime = (time_t) (meta->pds2.refTime +
1741 meta->pds2.sect4.foreSec);
1742 printf ("Most likely they didn't complete bytes 38-44 of "
1743 "Template 4.11\n");
1744 } else {
1745 meta->pds2.sect4.numInterval = (uChar) is4[44];
1746 }
1747
1748 /* Added this check because some MOS grids didn't finish the
1749 * template. */
1750 if (meta->pds2.sect4.numInterval != 0) {
1751 temp_ptr = realloc ((void *) meta->pds2.sect4.Interval,
1752 meta->pds2.sect4.numInterval *
1753 sizeof (sect4_IntervalType));
1754 if (temp_ptr == nullptr) {
1755 printf ("Ran out of memory.\n");
1756 return -6;
1757 }
1758 meta->pds2.sect4.Interval = (sect4_IntervalType *) temp_ptr;
1759 meta->pds2.sect4.numMissing = is4[45];
1760 if (ns4 < 57 + (meta->pds2.sect4.numInterval-1)*12+1) {
1761 return -1;
1762 }
1763 for (i = 0; i < meta->pds2.sect4.numInterval; i++) {
1764 meta->pds2.sect4.Interval[i].processID =
1765 (uChar) is4[49 + i * 12];
1766 meta->pds2.sect4.Interval[i].incrType =
1767 (uChar) is4[50 + i * 12];
1768 meta->pds2.sect4.Interval[i].timeRangeUnit =
1769 (uChar) is4[51 + i * 12];
1770 meta->pds2.sect4.Interval[i].lenTime = is4[52 + i * 12];
1771 meta->pds2.sect4.Interval[i].incrUnit =
1772 (uChar) is4[56 + i * 12];
1773 meta->pds2.sect4.Interval[i].timeIncr =
1774 (uChar) is4[57 + i * 12];
1775 }
1776 } else {
1777 #ifdef DEBUG
1778 printf ("Caution: Template 4.11 had no Intervals.\n");
1779 #endif
1780 meta->pds2.sect4.numMissing = is4[45];
1781 }
1782 break;
1783 case GS4_DERIVED: /* 4.2 */
1784 if (ns4 < 36) {
1785 return -1;
1786 }
1787 meta->pds2.sect4.derivedFcst = (uChar) is4[34];
1788 meta->pds2.sect4.numberFcsts = (uChar) is4[35];
1789 break;
1790 case GS4_DERIVED_CLUSTER_RECTANGULAR_AREA: /* 4.3 */
1791 if (ns4 < 68) {
1792 return -1;
1793 }
1794 meta->pds2.sect4.derivedFcst = (uChar) is4[34];
1795 meta->pds2.sect4.numberFcsts = (uChar) is4[35];
1796 break;
1797 case GS4_DERIVED_CLUSTER_CIRCULAR_AREA: /* 4.4 */
1798 if (ns4 < 64) {
1799 return -1;
1800 }
1801 meta->pds2.sect4.derivedFcst = (uChar) is4[34];
1802 meta->pds2.sect4.numberFcsts = (uChar) is4[35];
1803 break;
1804 case GS4_DERIVED_INTERVAL: /* 4.12 */
1805 if (ns4 < 45) {
1806 return -1;
1807 }
1808 meta->pds2.sect4.derivedFcst = (uChar) is4[34];
1809 meta->pds2.sect4.numberFcsts = (uChar) is4[35];
1810
1811 if (ParseTime (&(meta->pds2.sect4.validTime), is4[36], is4[38],
1812 is4[39], is4[40], is4[41], is4[42]) != 0) {
1813 msg = errSprintf (nullptr);
1814 meta->pds2.sect4.numInterval = (uChar) is4[43];
1815 if (meta->pds2.sect4.numInterval != 1) {
1816 errSprintf ("ERROR: in call to ParseTime from ParseSect4\n%s",
1817 msg);
1818 errSprintf ("Most likely they didn't complete bytes 37-43 of "
1819 "Template 4.12\n");
1820 free (msg);
1821 meta->pds2.sect4.numInterval = 0;
1822 return -1;
1823 }
1824 printf ("Warning: in call to ParseTime from ParseSect4\n%s", msg);
1825 free (msg);
1826 meta->pds2.sect4.validTime = (time_t) (meta->pds2.refTime +
1827 meta->pds2.sect4.foreSec);
1828 printf ("Most likely they didn't complete bytes 37-43 of "
1829 "Template 4.12\n");
1830 } else {
1831 meta->pds2.sect4.numInterval = (uChar) is4[43];
1832 }
1833
1834 /* Added this check because some MOS grids didn't finish the
1835 * template. */
1836 if (meta->pds2.sect4.numInterval != 0) {
1837 temp_ptr = realloc ((void *) meta->pds2.sect4.Interval,
1838 meta->pds2.sect4.numInterval *
1839 sizeof (sect4_IntervalType));
1840 if (temp_ptr == nullptr) {
1841 printf ("Ran out of memory.\n");
1842 return -6;
1843 }
1844 meta->pds2.sect4.Interval = (sect4_IntervalType *) temp_ptr;
1845 meta->pds2.sect4.numMissing = is4[44];
1846 if (ns4 < 56 + (meta->pds2.sect4.numInterval-1)*12+1) {
1847 return -1;
1848 }
1849 for (i = 0; i < meta->pds2.sect4.numInterval; i++) {
1850 meta->pds2.sect4.Interval[i].processID =
1851 (uChar) is4[48 + i * 12];
1852 meta->pds2.sect4.Interval[i].incrType =
1853 (uChar) is4[49 + i * 12];
1854 meta->pds2.sect4.Interval[i].timeRangeUnit =
1855 (uChar) is4[50 + i * 12];
1856 meta->pds2.sect4.Interval[i].lenTime = is4[51 + i * 12];
1857 meta->pds2.sect4.Interval[i].incrUnit =
1858 (uChar) is4[55 + i * 12];
1859 meta->pds2.sect4.Interval[i].timeIncr =
1860 (uChar) is4[56 + i * 12];
1861 }
1862 } else {
1863 #ifdef DEBUG
1864 printf ("Caution: Template 4.12 had no Intervals.\n");
1865 #endif
1866 meta->pds2.sect4.numMissing = is4[44];
1867 }
1868 break;
1869 case GS4_DERIVED_INTERVAL_CLUSTER_RECTANGULAR_AREA: /* 4.13 */
1870 case GS4_DERIVED_INTERVAL_CLUSTER_CIRCULAR_AREA: /* 4.14 */
1871 if (ns4 < 36) {
1872 return -1;
1873 }
1874 meta->pds2.sect4.derivedFcst = (uChar) is4[34];
1875 meta->pds2.sect4.numberFcsts = (uChar) is4[35];
1876 break;
1877 case GS4_STATISTIC: /* 4.8 */
1878 if (ns4 < 43) {
1879 return -1;
1880 }
1881 if (ParseTime (&(meta->pds2.sect4.validTime), is4[34], is4[36],
1882 is4[37], is4[38], is4[39], is4[40]) != 0) {
1883 msg = errSprintf (nullptr);
1884 meta->pds2.sect4.numInterval = (uChar) is4[41];
1885 if (meta->pds2.sect4.numInterval != 1) {
1886 errSprintf ("ERROR: in call to ParseTime from ParseSect4\n%s",
1887 msg);
1888 errSprintf ("Most likely they didn't complete bytes 35-41 of "
1889 "Template 4.8\n");
1890 free (msg);
1891 meta->pds2.sect4.numInterval = 0;
1892 return -1;
1893 }
1894 printf ("Warning: in call to ParseTime from ParseSect4\n%s", msg);
1895 free (msg);
1896 meta->pds2.sect4.validTime = (time_t) (meta->pds2.refTime +
1897 meta->pds2.sect4.foreSec);
1898 printf ("Most likely they didn't complete bytes 35-41 of "
1899 "Template 4.8\n");
1900 } else {
1901 meta->pds2.sect4.numInterval = (uChar) is4[41];
1902 }
1903
1904 /* Added this check because some MOS grids didn't finish the
1905 * template. */
1906 if (meta->pds2.sect4.numInterval != 0) {
1907 temp_ptr = realloc ((void *) meta->pds2.sect4.Interval,
1908 meta->pds2.sect4.numInterval *
1909 sizeof (sect4_IntervalType));
1910 if (temp_ptr == nullptr) {
1911 printf ("Ran out of memory.\n");
1912 return -6;
1913 }
1914 meta->pds2.sect4.Interval = (sect4_IntervalType *) temp_ptr;
1915 meta->pds2.sect4.numMissing = is4[42];
1916 if (ns4 < 54 + (meta->pds2.sect4.numInterval-1)*12+1) {
1917 return -1;
1918 }
1919 for (i = 0; i < meta->pds2.sect4.numInterval; i++) {
1920 meta->pds2.sect4.Interval[i].processID =
1921 (uChar) is4[46 + i * 12];
1922 meta->pds2.sect4.Interval[i].incrType =
1923 (uChar) is4[47 + i * 12];
1924 meta->pds2.sect4.Interval[i].timeRangeUnit =
1925 (uChar) is4[48 + i * 12];
1926 meta->pds2.sect4.Interval[i].lenTime = is4[49 + i * 12];
1927 meta->pds2.sect4.Interval[i].incrUnit =
1928 (uChar) is4[53 + i * 12];
1929 meta->pds2.sect4.Interval[i].timeIncr =
1930 (uChar) is4[54 + i * 12];
1931 }
1932 } else {
1933 #ifdef DEBUG
1934 printf ("Caution: Template 4.8 had no Intervals.\n");
1935 #endif
1936 meta->pds2.sect4.numMissing = is4[42];
1937 }
1938 break;
1939 case GS4_PERCENT_PNT: /* 4.6 */
1940 if( ns4 < 35) {
1941 return -1;
1942 }
1943 meta->pds2.sect4.percentile = is4[34];
1944 break;
1945 case GS4_PERCENT_TIME: /* 4.10 */
1946 if (ns4 < 44) {
1947 return -1;
1948 }
1949 meta->pds2.sect4.percentile = is4[34];
1950 if (ParseTime (&(meta->pds2.sect4.validTime), is4[35], is4[37],
1951 is4[38], is4[39], is4[40], is4[41]) != 0) {
1952 msg = errSprintf (nullptr);
1953 meta->pds2.sect4.numInterval = (uChar) is4[42];
1954 if (meta->pds2.sect4.numInterval != 1) {
1955 errSprintf ("ERROR: in call to ParseTime from ParseSect4\n%s",
1956 msg);
1957 errSprintf ("Most likely they didn't complete bytes 35-41 of "
1958 "Template 4.10\n");
1959 free (msg);
1960 meta->pds2.sect4.numInterval = 0;
1961 return -1;
1962 }
1963 printf ("Warning: in call to ParseTime from ParseSect4\n%s", msg);
1964 free (msg);
1965 meta->pds2.sect4.validTime = (time_t) (meta->pds2.refTime +
1966 meta->pds2.sect4.foreSec);
1967 printf ("Most likely they didn't complete bytes 35-41 of "
1968 "Template 4.10\n");
1969 } else {
1970 meta->pds2.sect4.numInterval = (uChar) is4[42];
1971 }
1972
1973 /* Added this check because some MOS grids didn't finish the
1974 * template. */
1975 if (meta->pds2.sect4.numInterval != 0) {
1976 temp_ptr = realloc ((void *) meta->pds2.sect4.Interval,
1977 meta->pds2.sect4.numInterval *
1978 sizeof (sect4_IntervalType));
1979 if (temp_ptr == nullptr) {
1980 printf ("Ran out of memory.\n");
1981 return -6;
1982 }
1983 meta->pds2.sect4.Interval = (sect4_IntervalType *) temp_ptr;
1984 meta->pds2.sect4.numMissing = is4[43];
1985 if (ns4 < 55 + (meta->pds2.sect4.numInterval-1)*12+1) {
1986 return -1;
1987 }
1988 for (i = 0; i < meta->pds2.sect4.numInterval; i++) {
1989 meta->pds2.sect4.Interval[i].processID =
1990 (uChar) is4[47 + i * 12];
1991 meta->pds2.sect4.Interval[i].incrType =
1992 (uChar) is4[48 + i * 12];
1993 meta->pds2.sect4.Interval[i].timeRangeUnit =
1994 (uChar) is4[49 + i * 12];
1995 meta->pds2.sect4.Interval[i].lenTime = is4[50 + i * 12];
1996 meta->pds2.sect4.Interval[i].incrUnit =
1997 (uChar) is4[54 + i * 12];
1998 meta->pds2.sect4.Interval[i].timeIncr =
1999 (uChar) is4[55 + i * 12];
2000 }
2001 } else {
2002 #ifdef DEBUG
2003 printf ("Caution: Template 4.10 had no Intervals.\n");
2004 #endif
2005 meta->pds2.sect4.numMissing = is4[43];
2006 }
2007 break;
2008 case GS4_PROBABIL_PNT: /* 4.5 */
2009 if (ns4 < 44) {
2010 return -1;
2011 }
2012 meta->pds2.sect4.foreProbNum = (uChar) is4[34];
2013 meta->pds2.sect4.numForeProbs = (uChar) is4[35];
2014 meta->pds2.sect4.probType = (uChar) is4[36];
2015 meta->pds2.sect4.lowerLimit.factor =
2016 sbit_2Comp_oneByte((sChar) is4[37]);
2017 meta->pds2.sect4.lowerLimit.value = sbit_2Comp_fourByte(is4[38]);
2018 meta->pds2.sect4.upperLimit.factor =
2019 sbit_2Comp_oneByte((sChar) is4[42]);
2020 meta->pds2.sect4.upperLimit.value = sbit_2Comp_fourByte(is4[43]);
2021 break;
2022 case GS4_PROBABIL_TIME: /* 4.9 */
2023 if (ns4 < 56) {
2024 return -1;
2025 }
2026 meta->pds2.sect4.foreProbNum = (uChar) is4[34];
2027 meta->pds2.sect4.numForeProbs = (uChar) is4[35];
2028 meta->pds2.sect4.probType = (uChar) is4[36];
2029 meta->pds2.sect4.lowerLimit.factor =
2030 sbit_2Comp_oneByte((sChar) is4[37]);
2031 meta->pds2.sect4.lowerLimit.value = sbit_2Comp_fourByte(is4[38]);
2032 meta->pds2.sect4.upperLimit.factor =
2033 sbit_2Comp_oneByte((sChar) is4[42]);
2034 meta->pds2.sect4.upperLimit.value = sbit_2Comp_fourByte(is4[43]);
2035 if (ParseTime (&(meta->pds2.sect4.validTime), is4[47], is4[49],
2036 is4[50], is4[51], is4[52], is4[53]) != 0) {
2037 msg = errSprintf (nullptr);
2038 meta->pds2.sect4.numInterval = (uChar) is4[54];
2039 if (meta->pds2.sect4.numInterval != 1) {
2040 errSprintf ("ERROR: in call to ParseTime from ParseSect4\n%s",
2041 msg);
2042 errSprintf ("Most likely they didn't complete bytes 48-54 of "
2043 "Template 4.9\n");
2044 free (msg);
2045 meta->pds2.sect4.numInterval = 0;
2046 return -1;
2047 }
2048 printf ("Warning: in call to ParseTime from ParseSect4\n%s", msg);
2049 free (msg);
2050 meta->pds2.sect4.validTime = (time_t) (meta->pds2.refTime +
2051 meta->pds2.sect4.foreSec);
2052 printf ("Most likely they didn't complete bytes 48-54 of "
2053 "Template 4.9\n");
2054 } else {
2055 meta->pds2.sect4.numInterval = (uChar) is4[54];
2056 }
2057 temp_ptr = realloc ((void *) meta->pds2.sect4.Interval,
2058 meta->pds2.sect4.numInterval *
2059 sizeof (sect4_IntervalType));
2060 if (temp_ptr == nullptr) {
2061 printf ("Ran out of memory.\n");
2062 return -6;
2063 }
2064 meta->pds2.sect4.Interval = (sect4_IntervalType *) temp_ptr;
2065 meta->pds2.sect4.numMissing = is4[55];
2066 if (ns4 < 67 + (meta->pds2.sect4.numInterval-1)*12+1) {
2067 return -1;
2068 }
2069 for (i = 0; i < meta->pds2.sect4.numInterval; i++) {
2070 meta->pds2.sect4.Interval[i].processID = (uChar) is4[59 + i * 12];
2071 meta->pds2.sect4.Interval[i].incrType = (uChar) is4[60 + i * 12];
2072 meta->pds2.sect4.Interval[i].timeRangeUnit =
2073 (uChar) is4[61 + i * 12];
2074 meta->pds2.sect4.Interval[i].lenTime = is4[62 + i * 12];
2075 meta->pds2.sect4.Interval[i].incrUnit = (uChar) is4[66 + i * 12];
2076 meta->pds2.sect4.Interval[i].timeIncr = (uChar) is4[67 + i * 12];
2077 }
2078 break;
2079 case GS4_STATISTIC_SPATIAL_AREA: /* 4.15 */
2080 // TODO. Need to fetch
2081 // 35 Statistical process used within the spatial area defined by octet 36 (see Code Table 4.10)
2082 // 36 Type of spatial processing used to arrive at given data value from source data (see Code Table 4.15)
2083 // 37 Number of data points used in spatial processing defined in octet 36
2084 break;
2085 case GS4_ANALYSIS_CHEMICAL: /* 4.40 */
2086 // TODO
2087 break;
2088 case GS4_OPTICAL_PROPERTIES_AEROSOL: /* 4.48 */
2089 // TODO
2090 break;
2091 default:
2092 errSprintf ("Un-supported Template. %ld\n", is4[7]);
2093 return -4;
2094 }
2095
2096 /* Do only that check at the end so that other meta fields are properly set */
2097 /* otherwise we might do erroneous unit conversion as in */
2098 /* https://github.com/OSGeo/gdal/issues/3158 */
2099 if (is4[5] != 0) {
2100 #ifdef DEBUG
2101 printf ("Un-supported template.\n All Supported template "
2102 "have 0 coordinate vertical values after template.");
2103 #endif
2104 errSprintf ("Un-supported template.\n All Supported template "
2105 "have 0 coordinate vertical values after template.");
2106 return -4;
2107 }
2108
2109 return 0;
2110 }
2111
2112 /*****************************************************************************
2113 * ParseSect5() --
2114 *
2115 * Arthur Taylor / MDL
2116 *
2117 * PURPOSE
2118 * To verify and parse section 5 data.
2119 *
2120 * ARGUMENTS
2121 * is5 = The unpacked section 5 array. (Input)
2122 * ns5 = The size of section 5. (Input)
2123 * meta = The structure to fill. (Output)
2124 * xmissp = The primary missing value. (Input)
2125 * xmisss = The secondary missing value. (Input)
2126 *
2127 * FILES/DATABASES: None
2128 *
2129 * RETURNS: int (could use errSprintf())
2130 * 0 = OK
2131 * -1 = ns5 is too small.
2132 * -2 = unexpected values in is5.
2133 * -6 = unsupported packing.
2134 *
2135 * HISTORY
2136 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
2137 *
2138 * NOTES
2139 *****************************************************************************
2140 */
ParseSect5(sInt4 * is5,sInt4 ns5,grib_MetaData * meta,float xmissp,float xmisss)2141 static int ParseSect5 (sInt4 *is5, sInt4 ns5, grib_MetaData *meta,
2142 float xmissp, float xmisss)
2143 {
2144 if (ns5 < 22) {
2145 return -1;
2146 }
2147 if (is5[4] != 5) {
2148 errSprintf ("ERROR IS5 not labeled correctly. %ld\n", is5[5]);
2149 return -2;
2150 }
2151 if ((is5[9] != GS5_SIMPLE) && (is5[9] != GS5_CMPLX) &&
2152 (is5[9] != GS5_CMPLXSEC) && (is5[9] != GS5_IEEE) &&
2153 (is5[9] != GS5_SPECTRAL) &&
2154 (is5[9] != GS5_HARMONIC) && (is5[9] != GS5_JPEG2000) &&
2155 (is5[9] != GS5_PNG) && (is5[9] != GS5_JPEG2000_ORG) &&
2156 (is5[9] != GS5_PNG_ORG)) {
2157 errSprintf ("Un-supported Packing? %ld\n", is5[9]);
2158 return -6;
2159 }
2160 meta->gridAttrib.packType = (sInt4) is5[9];
2161 meta->gridAttrib.f_maxmin = 0;
2162 meta->gridAttrib.missPri = xmissp;
2163 meta->gridAttrib.missSec = xmisss;
2164 if ( (is5[9] == GS5_IEEE) || (is5[9] == GS5_SPECTRAL) || (is5[9] == GS5_HARMONIC)) {
2165 meta->gridAttrib.fieldType = 0;
2166 meta->gridAttrib.f_miss = 0;
2167 return 0;
2168 }
2169 if (is5[20] > 1) {
2170 errSprintf ("Invalid field type. %ld\n", is5[20]);
2171 return -2;
2172 }
2173 MEMCPY_BIG (&meta->gridAttrib.refVal, &(is5[11]), 4);
2174 meta->gridAttrib.ESF = is5[15];
2175 meta->gridAttrib.DSF = is5[17];
2176 meta->gridAttrib.fieldType = (uChar) is5[20];
2177 if ((is5[9] == GS5_SIMPLE) ||
2178 (is5[9] == GS5_JPEG2000) || (is5[9] == GS5_JPEG2000_ORG) ||
2179 (is5[9] == GS5_PNG) || (is5[9] == GS5_PNG_ORG)) {
2180 meta->gridAttrib.f_miss = 0;
2181 return 0;
2182 }
2183
2184 myAssert( (is5[9] == GS5_CMPLX) || (is5[9] == GS5_CMPLXSEC) );
2185
2186 if (ns5 < 23) {
2187 return -1;
2188 }
2189 if (is5[22] > 2) {
2190 errSprintf ("Invalid missing management type, f_miss = %ld\n",
2191 is5[22]);
2192 return -2;
2193 }
2194 meta->gridAttrib.f_miss = (uChar) is5[22];
2195
2196 return 0;
2197 }
2198
2199 /*****************************************************************************
2200 * MetaParse() --
2201 *
2202 * Arthur Taylor / MDL
2203 *
2204 * PURPOSE
2205 * To parse all the meta data from a grib2 message.
2206 *
2207 * ARGUMENTS
2208 * meta = The structure to fill. (Output)
2209 * is0 = The unpacked section 0 array. (Input)
2210 * ns0 = The size of section 0. (Input)
2211 * is1 = The unpacked section 1 array. (Input)
2212 * ns1 = The size of section 1. (Input)
2213 * is2 = The unpacked section 2 array. (Input)
2214 * ns2 = The size of section 2. (Input)
2215 * rdat = The float data in section 2. (Input)
2216 * nrdat = Length of rdat. (Input)
2217 * idat = The integer data in section 2. (Input)
2218 * nidat = Length of idat. (Input)
2219 * is3 = The unpacked section 3 array. (Input)
2220 * ns3 = The size of section 3. (Input)
2221 * is4 = The unpacked section 4 array. (Input)
2222 * ns4 = The size of section 4. (Input)
2223 * is5 = The unpacked section 5 array. (Input)
2224 * ns5 = The size of section 5. (Input)
2225 * grib_len = The length of the entire grib message. (Input)
2226 * xmissp = The primary missing value. (Input)
2227 * xmisss = The secondary missing value. (Input)
2228 * simpVer = The version of the simple weather code to use when parsing the
2229 * WxString (if applicable). (Input)
2230 *
2231 * FILES/DATABASES: None
2232 *
2233 * RETURNS: int (could use errSprintf())
2234 * 0 = OK
2235 * -1 = A dimension is too small.
2236 * -2 = unexpected values in a grib section.
2237 * -3 = un-supported map Projection.
2238 * -4 = un-supported Sect 4 template.
2239 * -5 = unsupported forecast time unit.
2240 * -6 = unsupported sect 5 packing.
2241 * -10 = Something the driver can't handle yet.
2242 * (prodType != 0, f_sphere != 1, etc)
2243 * -11 = Weather grid without a lookup table.
2244 *
2245 * HISTORY
2246 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
2247 *
2248 * NOTES
2249 *****************************************************************************
2250 */
MetaParse(grib_MetaData * meta,sInt4 * is0,sInt4 ns0,sInt4 * is1,sInt4 ns1,sInt4 * is2,sInt4 ns2,float * rdat,sInt4 nrdat,sInt4 * idat,sInt4 nidat,sInt4 * is3,sInt4 ns3,sInt4 * is4,sInt4 ns4,sInt4 * is5,sInt4 ns5,sInt4 grib_len,float xmissp,float xmisss,int simpVer,CPL_UNUSED int simpWWA)2251 int MetaParse (grib_MetaData *meta, sInt4 *is0, sInt4 ns0,
2252 sInt4 *is1, sInt4 ns1, sInt4 *is2, sInt4 ns2,
2253 float *rdat, sInt4 nrdat, sInt4 *idat, sInt4 nidat,
2254 sInt4 *is3, sInt4 ns3, sInt4 *is4, sInt4 ns4,
2255 sInt4 *is5, sInt4 ns5, sInt4 grib_len,
2256 float xmissp, float xmisss, int simpVer, CPL_UNUSED int simpWWA)
2257 {
2258 int ierr; /* The error code of a called routine */
2259 /* char *element; *//* Holds the name of the current variable. */
2260 /* char *comment; *//* Holds more comments about current variable. */
2261 /* char *unitName; *//* Holds the name of the unit [K] [%] .. etc */
2262 uChar probType; /* The probability type */
2263 double lowerProb; /* The lower limit on probability forecast if
2264 * template 4.5 or 4.9 */
2265 double upperProb; /* The upper limit on probability forecast if
2266 * template 4.5 or 4.9 */
2267 sInt4 lenTime; /* Length of time for element (see 4.8 and 4.9) */
2268 uChar timeRangeUnit = 1;
2269 uChar incrType;
2270 uChar statProcessID; /* Statistical process id or 255 for missing */
2271 uChar fstSurfType; /* Type of the first fixed surface. */
2272 sInt4 value; /* The scaled value from GRIB2 file. */
2273 sChar scale; /* Surface scale as opposed to probility factor. */
2274 double fstSurfValue; /* Value of first fixed surface. */
2275 sChar f_fstValue; /* flag if FstValue is valid. */
2276 uChar sndSurfType; /* Type of the second fixed surface. */
2277 double sndSurfValue; /* Value of second fixed surface. */
2278 sChar f_sndValue; /* flag if SndValue is valid. */
2279
2280 if ((ierr = ParseSect0 (is0, ns0, grib_len, meta)) != 0) {
2281 preErrSprintf ("Parse error Section 0\n");
2282 //return ierr;
2283 }
2284 if ((ierr = ParseSect1 (is1, ns1, meta)) != 0) {
2285 preErrSprintf ("Parse error Section 1\n");
2286 //return ierr;
2287 }
2288 if (ns2 < 7) {
2289 errSprintf ("ns2 was too small in MetaParse\n");
2290 //return -1;
2291 }
2292 meta->pds2.f_sect2 = (uChar) (is2[0] != 0);
2293 if (meta->pds2.f_sect2) {
2294 meta->pds2.sect2NumGroups = is2[7 - 1];
2295 } else {
2296 meta->pds2.sect2NumGroups = 0;
2297 }
2298 if ((ierr = ParseSect3 (is3, ns3, meta)) != 0) {
2299 preErrSprintf ("Parse error Section 3\n");
2300 //return ierr;
2301 }
2302 if (IsData_NDFD (meta->center, meta->subcenter)) {
2303 meta->gds.hdatum = 1;
2304 }
2305 if (meta->gds.f_sphere != 1) {
2306 errSprintf ("Driver Filter: Can only handle spheres.\n");
2307 //return -10;
2308 }
2309 if ((ierr = ParseSect4 (is4, ns4, meta)) != 0) {
2310 preErrSprintf ("Parse error Section 4\n");
2311 //return ierr;
2312 }
2313 if ((ierr = ParseSect5 (is5, ns5, meta, xmissp, xmisss)) != 0) {
2314 preErrSprintf ("Parse error Section 5\n");
2315 //return ierr;
2316 }
2317 /* Compute ElementName. */
2318 if (meta->element) {
2319 free (meta->element);
2320 meta->element = nullptr;
2321 }
2322 if (meta->unitName) {
2323 free (meta->unitName);
2324 meta->unitName = nullptr;
2325 }
2326 if (meta->comment) {
2327 free (meta->comment);
2328 meta->comment = nullptr;
2329 }
2330
2331 if ((meta->pds2.sect4.templat == GS4_PROBABIL_TIME) ||
2332 (meta->pds2.sect4.templat == GS4_PROBABIL_PNT)) {
2333 probType = meta->pds2.sect4.probType;
2334 lowerProb = meta->pds2.sect4.lowerLimit.value *
2335 pow (10.0, -1 * meta->pds2.sect4.lowerLimit.factor);
2336 upperProb = meta->pds2.sect4.upperLimit.value *
2337 pow (10.0, -1 * meta->pds2.sect4.upperLimit.factor);
2338 } else {
2339 probType = 0;
2340 lowerProb = 0;
2341 upperProb = 0;
2342 }
2343 if (meta->pds2.sect4.numInterval > 0) {
2344 /* Try to convert lenTime to hourly. */
2345 timeRangeUnit = meta->pds2.sect4.Interval[0].timeRangeUnit;
2346 if (meta->pds2.sect4.Interval[0].timeRangeUnit == 255) {
2347 lenTime = (sInt4) ((meta->pds2.sect4.validTime -
2348 meta->pds2.sect4.foreSec -
2349 meta->pds2.refTime) / 3600);
2350 } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 0) {
2351 lenTime = (sInt4) (meta->pds2.sect4.Interval[0].lenTime / 60.);
2352 timeRangeUnit = 1;
2353 } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 1) {
2354 lenTime = meta->pds2.sect4.Interval[0].lenTime;
2355 timeRangeUnit = 1;
2356 } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 2) {
2357 lenTime = meta->pds2.sect4.Interval[0].lenTime * 24;
2358 timeRangeUnit = 1;
2359 } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 10) {
2360 lenTime = meta->pds2.sect4.Interval[0].lenTime * 3;
2361 timeRangeUnit = 1;
2362 } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 11) {
2363 lenTime = meta->pds2.sect4.Interval[0].lenTime * 6;
2364 timeRangeUnit = 1;
2365 } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 12) {
2366 lenTime = meta->pds2.sect4.Interval[0].lenTime * 12;
2367 timeRangeUnit = 1;
2368 } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 13) {
2369 lenTime = (sInt4) (meta->pds2.sect4.Interval[0].lenTime / 3600.);
2370 timeRangeUnit = 1;
2371 } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 3) { /* month */
2372 lenTime = meta->pds2.sect4.Interval[0].lenTime;
2373 timeRangeUnit = 3;
2374 } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 4) { /* year */
2375 lenTime = meta->pds2.sect4.Interval[0].lenTime;
2376 timeRangeUnit = 4;
2377 } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 5) { /* decade */
2378 lenTime = meta->pds2.sect4.Interval[0].lenTime * 10;
2379 timeRangeUnit = 4;
2380 } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 6) { /* normal */
2381 lenTime = meta->pds2.sect4.Interval[0].lenTime * 30;
2382 timeRangeUnit = 4;
2383 } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 7) { /* century */
2384 lenTime = meta->pds2.sect4.Interval[0].lenTime * 100;
2385 timeRangeUnit = 4;
2386 } else {
2387 lenTime = 0;
2388 printf ("Can't handle this timeRangeUnit\n");
2389 myAssert (meta->pds2.sect4.Interval[0].timeRangeUnit == 1);
2390 }
2391 /*
2392 } else {
2393 lenTime = 255;
2394 }
2395 if (lenTime == 255) {
2396 lenTime = (meta->pds2.sect4.validTime - meta->pds2.sect4.foreSec -
2397 meta->pds2.refTime) / 3600;
2398 }
2399 */
2400 if (lenTime == GRIB2MISSING_s4) {
2401 lenTime = 0;
2402 }
2403 incrType = meta->pds2.sect4.Interval[0].incrType;
2404 statProcessID = meta->pds2.sect4.Interval[0].processID;
2405 } else {
2406 lenTime = 0;
2407 timeRangeUnit = 1;
2408 incrType = 255;
2409 statProcessID = 255;
2410 }
2411
2412 if ((meta->pds2.sect4.templat == GS4_RADAR) || (meta->pds2.sect4.templat == GS4_SATELLITE)
2413 || (meta->pds2.sect4.templat == 254) || (meta->pds2.sect4.templat == 1000) || (meta->pds2.sect4.templat == 1001)
2414 || (meta->pds2.sect4.templat == 1002)) {
2415 fstSurfValue = 0;
2416 f_fstValue = 0;
2417 fstSurfType = 0;
2418 sndSurfValue = 0;
2419 f_sndValue = 0;
2420 } else {
2421 fstSurfType = meta->pds2.sect4.fstSurfType;
2422 scale = meta->pds2.sect4.fstSurfScale;
2423 if (meta->pds2.sect4.fstSurfValue >= std::numeric_limits<sInt4>::max() ||
2424 meta->pds2.sect4.fstSurfValue <= std::numeric_limits<sInt4>::min()) {
2425 // Out of range, so just call it missing.
2426 preErrSprintf ("fstSurfValue out of range\n");
2427 value = GRIB2MISSING_s4;
2428 } else {
2429 value = static_cast<sInt4>(meta->pds2.sect4.fstSurfValue);
2430 }
2431 if ((value == GRIB2MISSING_s4) || (scale == GRIB2MISSING_s1) ||
2432 (fstSurfType == GRIB2MISSING_u1)) {
2433 fstSurfValue = 0;
2434 f_fstValue = 1;
2435 } else {
2436 fstSurfValue = value * pow (10.0, -1 * scale);
2437 f_fstValue = 1;
2438 }
2439 sndSurfType = meta->pds2.sect4.sndSurfType;
2440 scale = meta->pds2.sect4.sndSurfScale;
2441 if (meta->pds2.sect4.sndSurfValue < std::numeric_limits<int>::max() &&
2442 meta->pds2.sect4.sndSurfValue > std::numeric_limits<int>::min()) {
2443 value = static_cast<int>(meta->pds2.sect4.sndSurfValue);
2444 } else {
2445 // sndSurfValue is out of range, so just call it missing.
2446 // TODO(schwehr): Consider using a tmp double if the scale will
2447 // make the resulting sndSurfValue be within range.
2448 preErrSprintf ("sndSurfValue out of range\n");
2449 value = GRIB2MISSING_s4;
2450 }
2451 if ((value == GRIB2MISSING_s4) || (scale == GRIB2MISSING_s1) ||
2452 (sndSurfType == GRIB2MISSING_u1)) {
2453 sndSurfValue = 0;
2454 f_sndValue = 0;
2455 } else {
2456 sndSurfValue = value * pow (10.0, -1 * scale);
2457 f_sndValue = 1;
2458 }
2459 }
2460
2461 ParseElemName (meta->pds2.mstrVersion, meta->center, meta->subcenter, meta->pds2.prodType,
2462 meta->pds2.sect4.templat, meta->pds2.sect4.cat,
2463 meta->pds2.sect4.subcat, lenTime, timeRangeUnit, statProcessID,
2464 incrType, meta->pds2.sect4.genID, probType, lowerProb, upperProb,
2465 meta->pds2.sect4.derivedFcst,
2466 &(meta->element), &(meta->comment), &(meta->unitName),
2467 &(meta->convert), meta->pds2.sect4.percentile,
2468 meta->pds2.sect4.genProcess,
2469 f_fstValue, fstSurfValue, f_sndValue, sndSurfValue);
2470 #ifdef DEBUG
2471 /*
2472 printf ("Element: %s\nunitName: %s\ncomment: %s\n", meta->element,
2473 meta->comment, meta->unitName);
2474 */
2475 #endif
2476
2477 if (! f_fstValue) {
2478 reallocSprintf (&(meta->shortFstLevel), "0 undefined");
2479 reallocSprintf (&(meta->longFstLevel), "0.000[-] undefined ()");
2480 } else {
2481 ParseLevelName (meta->center, meta->subcenter, fstSurfType,
2482 fstSurfValue, f_sndValue, sndSurfValue,
2483 &(meta->shortFstLevel), &(meta->longFstLevel));
2484 }
2485
2486 /* Continue parsing section 2 data. */
2487 if (meta->pds2.f_sect2) {
2488 MetaSect2Free (meta);
2489 if (strcmp (meta->element, "Wx") == 0) {
2490 meta->pds2.sect2.ptrType = GS2_WXTYPE;
2491 if ((ierr = ParseSect2_Wx (rdat, nrdat, idat, nidat,
2492 &(meta->pds2.sect2.wx), simpVer)) != 0) {
2493 preErrSprintf ("Parse error Section 2 : Weather Data\n");
2494 return ierr;
2495 }
2496 } else if (strcmp (meta->element, "WWA") == 0) {
2497 meta->pds2.sect2.ptrType = GS2_HAZARD;
2498 if ((ierr = ParseSect2_Hazard (rdat, nrdat, idat, nidat,
2499 &(meta->pds2.sect2.hazard), simpWWA)) != 0) {
2500 preErrSprintf ("Parse error Section 2 : Hazard Data\n");
2501 return ierr;
2502 }
2503 } else {
2504 meta->pds2.sect2.ptrType = GS2_UNKNOWN;
2505 if ((ierr = ParseSect2_Unknown (rdat, nrdat, idat, nidat, meta))
2506 != 0) {
2507 preErrSprintf ("Parse error Section 2 : Unknown Data type\n");
2508 //return ierr;
2509 }
2510 }
2511 } else {
2512 if (strcmp (meta->element, "Wx") == 0) {
2513 errSprintf ("Weather grid does not have look up table?");
2514 //return -11;
2515 }
2516 if (strcmp (meta->element, "WWA") == 0) {
2517 errSprintf ("Hazard grid does not have look up table?");
2518 //return -11;
2519 }
2520 }
2521 return 0;
2522 }
2523
2524 /*****************************************************************************
2525 * ParseGridNoMiss() --
2526 *
2527 * Arthur Taylor / MDL
2528 *
2529 * PURPOSE
2530 * A helper function for ParseGrid. In this particular case it is dealing
2531 * with a field that has NO missing value type.
2532 * Walks through either a float or an integer grid, computing the min/max
2533 * values in the grid, and converts the units. It uses gridAttrib info for the
2534 * missing values and it updates gridAttrib with the observed min/max values.
2535 *
2536 * ARGUMENTS
2537 * attrib = Grid Attribute structure already filled in (Input/Output)
2538 * grib_Data = The place to store the grid data. (Output)
2539 * Nx, Ny = The dimensions of the grid (Input)
2540 * iain = Place to find data if it is an Integer (or float). (Input)
2541 * unitM = M in unit conversion equation y(new) = m x(orig) + b (Input)
2542 * unitB = B in unit conversion equation y(new) = m x(orig) + b (Input)
2543 * f_txtType = true if we have a valid wx/hazard type. (Input)
2544 * txt_dataLen = Length of text table
2545 * txt_f_valid = whether that entry is used/valid. (Input)
2546 * startX = The start of the X values. (Input)
2547 * startY = The start of the Y values. (Input)
2548 * subNx = The Nx dimension of the subgrid (Input)
2549 * subNy = The Ny dimension of the subgrid (Input)
2550 *
2551 * FILES/DATABASES: None
2552 *
2553 * RETURNS: void
2554 *
2555 * HISTORY
2556 * 12/2002 Arthur Taylor (MDL/RSIS): Created to optimize part of ParseGrid.
2557 * 5/2003 AAT: Added ability to see if wxType occurs. If so sets table
2558 * valid to 2, otherwise leaves it at 1. If table valid is 0 then
2559 * sets value to missing value (if applicable).
2560 * 2/2004 AAT: Added the subgrid capability.
2561 *
2562 * NOTES
2563 * 1) Don't have to check if value became missing value, because we can check
2564 * if missing falls in the range of the min/max converted units. If
2565 * missing does fall in that range we need to move missing.
2566 * (See f_readjust in ParseGrid)
2567 *****************************************************************************
2568 */
ParseGridNoMiss(gridAttribType * attrib,double * grib_Data,sInt4 Nx,sInt4 Ny,sInt4 * iain,double unitM,double unitB,uChar f_txtType,uInt4 txt_dataLen,uChar * txt_f_valid,int startX,int startY,int subNx,int subNy)2569 static void ParseGridNoMiss (gridAttribType *attrib, double *grib_Data,
2570 sInt4 Nx, sInt4 Ny, sInt4 *iain,
2571 double unitM, double unitB,
2572 uChar f_txtType, uInt4 txt_dataLen,
2573 uChar *txt_f_valid, int startX, int startY,
2574 int subNx, int subNy)
2575 {
2576 sInt4 x, y; /* Where we are in the grid. */
2577 double value; /* The data in the new units. */
2578 uChar f_maxmin = 0; /* Flag if max/min is valid yet. */
2579 uInt4 index; /* Current index into Wx table. */
2580 sInt4 *itemp = nullptr;
2581 float *ftemp = nullptr;
2582
2583 /* Resolve possibility that the data is an integer or a float and find
2584 * max/min values. (see note 1) */
2585 for (y = 0; y < subNy; y++) {
2586 if (((startY + y - 1) < 0) || ((startY + y - 1) >= Ny)) {
2587 for (x = 0; x < subNx; x++) {
2588 *grib_Data++ = 9999;
2589 }
2590 } else {
2591 if (attrib->fieldType) {
2592 itemp = iain + (startY + y - 1) * Nx + (startX - 1);
2593 } else {
2594 ftemp = ((float *) iain) + (startY + y - 1) * Nx + (startX - 1);
2595 }
2596 for (x = 0; x < subNx; x++) {
2597 if (((startX + x - 1) < 0) || ((startX + x - 1) >= Nx)) {
2598 *grib_Data++ = 9999;
2599 } else {
2600 /* Convert the units. */
2601 if (attrib->fieldType) {
2602 if (unitM == -10) {
2603 value = pow (10.0, (*itemp++));
2604 } else {
2605 value = unitM * (*itemp++) + unitB;
2606 }
2607 } else {
2608 if (unitM == -10) {
2609 value = pow (10.0, (double) (*ftemp++));
2610 } else {
2611 value = unitM * (*ftemp++) + unitB;
2612 }
2613 }
2614 if (f_txtType) {
2615 index = (uInt4) value;
2616 if (index < txt_dataLen) {
2617 if (txt_f_valid[index] == 1) {
2618 txt_f_valid[index] = 2;
2619 } else if (txt_f_valid[index] == 0) {
2620 /* Table is not valid here so set value to missing? */
2621 /* No missing value, so use index = WxType->dataLen? */
2622 /* No... set f_valid to 3 so we know we used this
2623 * invalid element, then handle it in degrib2.c ::
2624 * ReadGrib2Record() where we set it back to 0. */
2625 txt_f_valid[index] = 3;
2626 }
2627 }
2628 }
2629 if (f_maxmin) {
2630 if (value < attrib->min) {
2631 attrib->min = value;
2632 } else if (value > attrib->max) {
2633 attrib->max = value;
2634 }
2635 } else {
2636 attrib->min = attrib->max = value;
2637 f_maxmin = 1;
2638 }
2639 *grib_Data++ = value;
2640 }
2641 }
2642 }
2643 }
2644 attrib->f_maxmin = f_maxmin;
2645 }
2646
2647 /*****************************************************************************
2648 * ParseGridPrimMiss() --
2649 *
2650 * Arthur Taylor / MDL
2651 *
2652 * PURPOSE
2653 * A helper function for ParseGrid. In this particular case it is dealing
2654 * with a field that has primary missing value type.
2655 * Walks through either a float or an integer grid, computing the min/max
2656 * values in the grid, and converts the units. It uses gridAttrib info for the
2657 * missing values and it updates gridAttrib with the observed min/max values.
2658 *
2659 * ARGUMENTS
2660 * attrib = sect 5 structure already filled in by ParseSect5 (In/Output)
2661 * grib_Data = The place to store the grid data. (Output)
2662 * Nx, Ny = The dimensions of the grid (Input)
2663 * iain = Place to find data if it is an Integer (or float). (Input)
2664 * unitM = M in unit conversion equation y(new) = m x(orig) + b (Input)
2665 * unitB = B in unit conversion equation y(new) = m x(orig) + b (Input)
2666 * f_txtType = true if we have a valid wx/hazard type. (Input)
2667 * txt_dataLen = Length of text table
2668 * txt_f_valid = whether that entry is used/valid. (Input)
2669 * startX = The start of the X values. (Input)
2670 * startY = The start of the Y values. (Input)
2671 * subNx = The Nx dimension of the subgrid (Input)
2672 * subNy = The Ny dimension of the subgrid (Input)
2673 *
2674 * FILES/DATABASES: None
2675 *
2676 * RETURNS: void
2677 *
2678 * HISTORY
2679 * 12/2002 Arthur Taylor (MDL/RSIS): Created to optimize part of ParseGrid.
2680 * 5/2003 AAT: Added ability to see if wxType occurs. If so sets table
2681 * valid to 2, otherwise leaves it at 1. If table valid is 0 then
2682 * sets value to missing value (if applicable).
2683 * 2/2004 AAT: Added the subgrid capability.
2684 *
2685 * NOTES
2686 * 1) Don't have to check if value became missing value, because we can check
2687 * if missing falls in the range of the min/max converted units. If
2688 * missing does fall in that range we need to move missing.
2689 * (See f_readjust in ParseGrid)
2690 *****************************************************************************
2691 */
ParseGridPrimMiss(gridAttribType * attrib,double * grib_Data,sInt4 Nx,sInt4 Ny,sInt4 * iain,double unitM,double unitB,sInt4 * missCnt,uChar f_txtType,uInt4 txt_dataLen,uChar * txt_f_valid,int startX,int startY,int subNx,int subNy)2692 static void ParseGridPrimMiss (gridAttribType *attrib, double *grib_Data,
2693 sInt4 Nx, sInt4 Ny, sInt4 *iain,
2694 double unitM, double unitB, sInt4 *missCnt,
2695 uChar f_txtType, uInt4 txt_dataLen,
2696 uChar *txt_f_valid,
2697 int startX, int startY, int subNx, int subNy)
2698 {
2699 sInt4 x, y; /* Where we are in the grid. */
2700 double value; /* The data in the new units. */
2701 uChar f_maxmin = 0; /* Flag if max/min is valid yet. */
2702 uInt4 index; /* Current index into Wx table. */
2703 sInt4 *itemp = nullptr;
2704 float *ftemp = nullptr;
2705 /* float *ain = (float *) iain;*/
2706
2707 /* Resolve possibility that the data is an integer or a float and find
2708 * max/min values. (see note 1) */
2709 for (y = 0; y < subNy; y++) {
2710 if (((startY + y - 1) < 0) || ((startY + y - 1) >= Ny)) {
2711 for (x = 0; x < subNx; x++) {
2712 *grib_Data++ = attrib->missPri;
2713 (*missCnt)++;
2714 }
2715 } else {
2716 if (attrib->fieldType) {
2717 itemp = iain + (startY + y - 1) * Nx + (startX - 1);
2718 } else {
2719 ftemp = ((float *) iain) + (startY + y - 1) * Nx + (startX - 1);
2720 }
2721 for (x = 0; x < subNx; x++) {
2722 if (((startX + x - 1) < 0) || ((startX + x - 1) >= Nx)) {
2723 *grib_Data++ = attrib->missPri;
2724 (*missCnt)++;
2725 } else {
2726 if (attrib->fieldType) {
2727 value = (*itemp++);
2728 } else {
2729 value = (*ftemp++);
2730 }
2731
2732 /* Make sure value is not a missing value when converting
2733 * units, and while computing max/min. */
2734 if (value == attrib->missPri) {
2735 (*missCnt)++;
2736 } else {
2737 /* Convert the units. */
2738 if (unitM == -10) {
2739 value = pow (10.0, value);
2740 } else {
2741 value = unitM * value + unitB;
2742 }
2743 if (f_txtType) {
2744 index = (uInt4) value;
2745 if (index < txt_dataLen) {
2746 if (txt_f_valid[index]) {
2747 txt_f_valid[index] = 2;
2748 } else {
2749 /* Table is not valid here so set value to missPri
2750 */
2751 value = attrib->missPri;
2752 (*missCnt)++;
2753 }
2754 }
2755 }
2756 if ((!f_txtType) || (value != attrib->missPri)) {
2757 if (f_maxmin) {
2758 if (value < attrib->min) {
2759 attrib->min = value;
2760 } else if (value > attrib->max) {
2761 attrib->max = value;
2762 }
2763 } else {
2764 attrib->min = attrib->max = value;
2765 f_maxmin = 1;
2766 }
2767 }
2768 }
2769 *grib_Data++ = value;
2770 }
2771 }
2772 }
2773 }
2774 attrib->f_maxmin = f_maxmin;
2775 }
2776
2777 /*****************************************************************************
2778 * ParseGridSecMiss() --
2779 *
2780 * Arthur Taylor / MDL
2781 *
2782 * PURPOSE
2783 * A helper function for ParseGrid. In this particular case it is dealing
2784 * with a field that has NO missing value type.
2785 * Walks through either a float or an integer grid, computing the min/max
2786 * values in the grid, and converts the units. It uses gridAttrib info for the
2787 * missing values and it updates gridAttrib with the observed min/max values.
2788 *
2789 * ARGUMENTS
2790 * attrib = sect 5 structure already filled in by ParseSect5 (In/Output)
2791 * grib_Data = The place to store the grid data. (Output)
2792 * Nx, Ny = The dimensions of the grid (Input)
2793 * iain = Place to find data if it is an Integer (or float). (Input)
2794 * unitM = M in unit conversion equation y(new) = m x(orig) + b (Input)
2795 * unitB = B in unit conversion equation y(new) = m x(orig) + b (Input)
2796 * f_txtType = true if we have a valid wx/hazard type. (Input)
2797 * txt_dataLen = Length of text table
2798 * txt_f_valid = whether that entry is used/valid. (Input)
2799 * startX = The start of the X values. (Input)
2800 * startY = The start of the Y values. (Input)
2801 * subNx = The Nx dimension of the subgrid (Input)
2802 * subNy = The Ny dimension of the subgrid (Input)
2803 *
2804 * FILES/DATABASES: None
2805 *
2806 * RETURNS: void
2807 *
2808 * HISTORY
2809 * 12/2002 Arthur Taylor (MDL/RSIS): Created to optimize part of ParseGrid.
2810 * 5/2003 AAT: Added ability to see if wxType occurs. If so sets table
2811 * valid to 2, otherwise leaves it at 1. If table valid is 0 then
2812 * sets value to missing value (if applicable).
2813 * 2/2004 AAT: Added the subgrid capability.
2814 *
2815 * NOTES
2816 * 1) Don't have to check if value became missing value, because we can check
2817 * if missing falls in the range of the min/max converted units. If
2818 * missing does fall in that range we need to move missing.
2819 * (See f_readjust in ParseGrid)
2820 *****************************************************************************
2821 */
ParseGridSecMiss(gridAttribType * attrib,double * grib_Data,sInt4 Nx,sInt4 Ny,sInt4 * iain,double unitM,double unitB,sInt4 * missCnt,uChar f_txtType,uInt4 txt_dataLen,uChar * txt_f_valid,int startX,int startY,int subNx,int subNy)2822 static void ParseGridSecMiss (gridAttribType *attrib, double *grib_Data,
2823 sInt4 Nx, sInt4 Ny, sInt4 *iain,
2824 double unitM, double unitB, sInt4 *missCnt,
2825 uChar f_txtType, uInt4 txt_dataLen,
2826 uChar *txt_f_valid,
2827 int startX, int startY, int subNx, int subNy)
2828 {
2829 sInt4 x, y; /* Where we are in the grid. */
2830 double value; /* The data in the new units. */
2831 uChar f_maxmin = 0; /* Flag if max/min is valid yet. */
2832 uInt4 index; /* Current index into Wx table. */
2833 sInt4 *itemp = nullptr;
2834 float *ftemp = nullptr;
2835 /* float *ain = (float *) iain;*/
2836
2837 /* Resolve possibility that the data is an integer or a float and find
2838 * max/min values. (see note 1) */
2839 for (y = 0; y < subNy; y++) {
2840 if (((startY + y - 1) < 0) || ((startY + y - 1) >= Ny)) {
2841 for (x = 0; x < subNx; x++) {
2842 *grib_Data++ = attrib->missPri;
2843 (*missCnt)++;
2844 }
2845 } else {
2846 if (attrib->fieldType) {
2847 itemp = iain + (startY + y - 1) * Nx + (startX - 1);
2848 } else {
2849 ftemp = ((float *) iain) + (startY + y - 1) * Nx + (startX - 1);
2850 }
2851 for (x = 0; x < subNx; x++) {
2852 if (((startX + x - 1) < 0) || ((startX + x - 1) >= Nx)) {
2853 *grib_Data++ = attrib->missPri;
2854 (*missCnt)++;
2855 } else {
2856 if (attrib->fieldType) {
2857 value = (*itemp++);
2858 } else {
2859 value = (*ftemp++);
2860 }
2861
2862 /* Make sure value is not a missing value when converting
2863 * units, and while computing max/min. */
2864 if ((value == attrib->missPri) || (value == attrib->missSec)) {
2865 (*missCnt)++;
2866 } else {
2867 /* Convert the units. */
2868 if (unitM == -10) {
2869 value = pow (10.0, value);
2870 } else {
2871 value = unitM * value + unitB;
2872 }
2873 if (f_txtType) {
2874 index = (uInt4) value;
2875 if (index < txt_dataLen) {
2876 if (txt_f_valid[index]) {
2877 txt_f_valid[index] = 2;
2878 } else {
2879 /* Table is not valid here so set value to missPri
2880 */
2881 value = attrib->missPri;
2882 (*missCnt)++;
2883 }
2884 }
2885 }
2886 if ((!f_txtType) || (value != attrib->missPri)) {
2887 if (f_maxmin) {
2888 if (value < attrib->min) {
2889 attrib->min = value;
2890 } else if (value > attrib->max) {
2891 attrib->max = value;
2892 }
2893 } else {
2894 attrib->min = attrib->max = value;
2895 f_maxmin = 1;
2896 }
2897 }
2898 }
2899 *grib_Data++ = value;
2900 }
2901 }
2902 }
2903 }
2904 attrib->f_maxmin = f_maxmin;
2905 }
2906
2907 /*****************************************************************************
2908 * ParseGrid() --
2909 *
2910 * Arthur Taylor / MDL
2911 *
2912 * PURPOSE
2913 * To walk through the 2 possible grids (and possible bitmap) created by
2914 * UNPK_GRIB2, and combine the info into 1 grid, at the same time computing
2915 * the min/max values in the grid. It uses gridAttrib info for the missing values
2916 * and it then updates the gridAttrib structure for the min/max values that it
2917 * found.
2918 * It also uses scan, and ScanIndex2XY, to parse the data and organize the
2919 * Grib_Data so that 0,0 is the lower left part of the grid, it then traverses
2920 * the row and then moved up to the next row starting on the left.
2921 *
2922 * ARGUMENTS
2923 * attrib = sect 5 structure already filled in by ParseSect5 (In/Output)
2924 * Grib_Data = The place to store the grid data. (Output)
2925 * grib_DataLen = The current size of Grib_Data (can increase) (Input/Output)
2926 * Nx, Ny = The dimensions of the grid (Input)
2927 * scan = How to walk through the original grid. (Input)
2928 * iain = Place to find data if it is an Integer (or float). (Input)
2929 * ibitmap = Flag stating the data has a bitmap for missing values (In)
2930 * ib = Where to find the bitmap if we have one (Input)
2931 * unitM = M in unit conversion equation y(new) = m x(orig) + b (Input)
2932 * unitB = B in unit conversion equation y(new) = m x(orig) + b (Input)
2933 * f_txtType = true if we have a valid wx/hazard type. (Input)
2934 * txt_dataLen = Length of text table
2935 * txt_f_valid = whether that entry is used/valid. (Input)
2936 * f_subGrid = True if we have a subgrid, false if not. (Input)
2937 * startX stopX = The bounds of the subgrid in X. (0,-1) means full grid (In)
2938 * startY stopY = The bounds of the subgrid in Y. (0,-1) means full grid (In)
2939 *
2940 * FILES/DATABASES: None
2941 *
2942 * RETURNS: void
2943 *
2944 * HISTORY
2945 * 9/2002 Arthur Taylor (MDL/RSIS): Created.
2946 * 11/2002 AAT: Added unit conversion to metaparse.c
2947 * 12/2002 AAT: Optimized first loop to make it assume scan 0100 (64)
2948 * (valid 99.9%), but still have slow loop for generic case.
2949 * 5/2003 AAT: Added ability to see if wxType occurs. If so sets table
2950 * valid to 2, otherwise leaves it at 1. If table valid is 0 then
2951 * sets value to missing value (if applicable).
2952 * 7/2003 AAT: added check if f_maxmin before checking if missing was in
2953 * range of max, min for "readjust" check.
2954 * 2/2004 AAT: Added startX / startY / stopX / stopY
2955 * 5/2004 AAT: Found out that I used the opposite definition for bitmap
2956 * 0 = missing, 1 = valid.
2957 *
2958 * NOTES
2959 *****************************************************************************
2960 */
ParseGrid(VSILFILE * fp,gridAttribType * attrib,double ** Grib_Data,uInt4 * grib_DataLen,uInt4 Nx,uInt4 Ny,int scan,sInt4 nd2x3,sInt4 * iain,sInt4 ibitmap,sInt4 * ib,double unitM,double unitB,uChar f_txtType,uInt4 txt_dataLen,uChar * txt_f_valid,CPL_UNUSED uChar f_subGrid,int startX,int startY,int stopX,int stopY)2961 void ParseGrid (VSILFILE *fp, gridAttribType *attrib, double **Grib_Data,
2962 uInt4 *grib_DataLen, uInt4 Nx, uInt4 Ny, int scan,
2963 sInt4 nd2x3, sInt4 *iain, sInt4 ibitmap, sInt4 *ib, double unitM,
2964 double unitB, uChar f_txtType, uInt4 txt_dataLen,
2965 uChar *txt_f_valid,
2966 CPL_UNUSED uChar f_subGrid,
2967 int startX, int startY, int stopX, int stopY)
2968 {
2969 double xmissp; /* computed missing value needed for ibitmap = 1,
2970 * Also used if unit conversion causes confusion
2971 * over_ missing values. */
2972 double xmisss; /* Used if unit conversion causes confusion over
2973 * missing values. */
2974 uChar f_readjust; /* True if unit conversion caused confusion over
2975 * missing values. */
2976 uInt4 scanIndex; /* Where we are in the original grid. */
2977 sInt4 x, y; /* Where we are in a grid of scan value 0100 */
2978 uInt4 newIndex; /* x,y in a 1 dimensional array. */
2979 double value; /* The data in the new units. */
2980 /* A pointer to Grib_Data for ease of manipulation. */
2981 double *grib_Data = nullptr;
2982 sInt4 missCnt = 0; /* Number of detected missing values. */
2983 uInt4 index; /* Current index into Wx table. */
2984 float *ain = (float *) iain;
2985 uInt4 subNx; /* The Nx dimension of the subgrid. */
2986 uInt4 subNy; /* The Ny dimension of the subgrid. */
2987
2988 subNx = stopX - startX + 1;
2989 subNy = stopY - startY + 1;
2990
2991 myAssert (((!f_subGrid) && (subNx == Nx)) || (f_subGrid));
2992 myAssert (((!f_subGrid) && (subNy == Ny)) || (f_subGrid));
2993
2994 if( subNy == 0 || subNx > UINT_MAX / subNy )
2995 {
2996 errSprintf ("Too large raster");
2997 *grib_DataLen = 0;
2998 *Grib_Data = nullptr;
2999 return;
3000 }
3001
3002 const uInt4 subNxNy = subNx * subNy;
3003 if (subNxNy > *grib_DataLen) {
3004
3005 if( subNxNy > 100 * 1024 * 1024 )
3006 {
3007 vsi_l_offset curPos = VSIFTellL(fp);
3008 VSIFSeekL(fp, 0, SEEK_END);
3009 vsi_l_offset fileSize = VSIFTellL(fp);
3010 VSIFSeekL(fp, curPos, SEEK_SET);
3011 // allow a compression ratio of 1:1000
3012 if( subNxNy / 1000 > fileSize )
3013 {
3014 errSprintf ("ERROR: File too short\n");
3015 *grib_DataLen = 0;
3016 *Grib_Data = nullptr;
3017 return;
3018 }
3019 }
3020
3021 double* newData = nullptr;
3022 const size_t nBufferSize = subNxNy * sizeof (double);
3023 if( nBufferSize / sizeof(double) == subNxNy )
3024 {
3025 newData = (double *) realloc ((void *) (*Grib_Data), nBufferSize);
3026 }
3027 if( newData == nullptr )
3028 {
3029 errSprintf ("Memory allocation failed");
3030 free(*Grib_Data);
3031 *Grib_Data = nullptr;
3032 *grib_DataLen = 0;
3033 return;
3034 }
3035 *grib_DataLen = subNxNy;
3036 *Grib_Data = newData;
3037 }
3038 grib_Data = *Grib_Data;
3039
3040 /* Resolve possibility that the data is an integer or a float, find
3041 * max/min values, and do unit conversion. (see note 1) */
3042 if (scan == 64) {
3043 if (attrib->f_miss == 0) {
3044 ParseGridNoMiss (attrib, grib_Data, Nx, Ny, iain, unitM, unitB,
3045 f_txtType, txt_dataLen, txt_f_valid, startX, startY, subNx, subNy);
3046 } else if (attrib->f_miss == 1) {
3047 ParseGridPrimMiss (attrib, grib_Data, Nx, Ny, iain, unitM, unitB,
3048 &missCnt, f_txtType, txt_dataLen, txt_f_valid, startX, startY,
3049 subNx, subNy);
3050 } else if (attrib->f_miss == 2) {
3051 ParseGridSecMiss (attrib, grib_Data, Nx, Ny, iain, unitM, unitB,
3052 &missCnt, f_txtType, txt_dataLen, txt_f_valid, startX, startY, subNx,
3053 subNy);
3054 }
3055 } else {
3056 /* Internally we use scan = 0100. Scan is usually 0100 from the
3057 * unpacker library, but if scan is not, the following code converts
3058 * it. We optimized the previous (scan 0100) case by calling a
3059 * dedicated procedure. Here we don't since for scan != 0100, we
3060 * would_ need a different unpacker library, which is extremely
3061 * unlikely. */
3062 for (scanIndex = 0; scanIndex < (uInt4)nd2x3 && scanIndex < Nx * Ny; scanIndex++) {
3063 if (attrib->fieldType) {
3064 value = iain[scanIndex];
3065 } else {
3066 value = ain[scanIndex];
3067 }
3068 /* Make sure value is not a missing value when converting units, and
3069 * while computing max/min. */
3070 if ((attrib->f_miss == 0) ||
3071 ((attrib->f_miss == 1) && (value != attrib->missPri)) ||
3072 ((attrib->f_miss == 2) && (value != attrib->missPri) &&
3073 (value != attrib->missSec))) {
3074 /* Convert the units. */
3075 if (unitM == -10) {
3076 value = pow (10.0, value);
3077 } else {
3078 value = unitM * value + unitB;
3079 }
3080 /* Don't have to check if value became missing value, because we
3081 * can check if missing falls in the range of min/max. If
3082 * missing does fall in that range we need to move missing. See
3083 * f_readjust */
3084 if (f_txtType) {
3085 index = (uInt4) value;
3086 if (index < txt_dataLen) {
3087 if (txt_f_valid[index] == 1) {
3088 txt_f_valid[index] = 2;
3089 } else if (txt_f_valid[index] == 0) {
3090 /* Table is not valid here so set value to missPri */
3091 if (attrib->f_miss != 0) {
3092 value = attrib->missPri;
3093 missCnt++;
3094 } else {
3095 /* No missing value, so use index = WxType->dataLen */
3096 /* No... set f_valid to 3 so we know we used this
3097 * invalid element, then handle it in degrib2.c ::
3098 * ReadGrib2Record() where we set it back to 0. */
3099 txt_f_valid[index] = 3;
3100 }
3101 }
3102 }
3103 }
3104 if ((!f_txtType) ||
3105 ((attrib->f_miss == 0) || (value != attrib->missPri))) {
3106 if (attrib->f_maxmin) {
3107 if (value < attrib->min) {
3108 attrib->min = value;
3109 } else if (value > attrib->max) {
3110 attrib->max = value;
3111 }
3112 } else {
3113 attrib->min = attrib->max = value;
3114 attrib->f_maxmin = 1;
3115 }
3116 }
3117 } else {
3118 missCnt++;
3119 }
3120 ScanIndex2XY (scanIndex, &x, &y, scan, Nx, Ny);
3121 /* ScanIndex returns value as if scan was 0100 */
3122 newIndex = (uInt4)(x - 1) + (uInt4)(y - 1) * Nx;
3123 grib_Data[newIndex] = value;
3124 }
3125 }
3126
3127 /* Deal with possibility that unit conversion ended up with valid numbers
3128 * being interpreted as missing. */
3129 f_readjust = 0;
3130 xmissp = attrib->missPri;
3131 xmisss = attrib->missSec;
3132 if (attrib->f_maxmin) {
3133 if ((attrib->f_miss == 1) || (attrib->f_miss == 2)) {
3134 if ((attrib->missPri >= attrib->min) &&
3135 (attrib->missPri <= attrib->max)) {
3136 xmissp = attrib->max + 1;
3137 f_readjust = 1;
3138 }
3139 if (attrib->f_miss == 2) {
3140 if ((attrib->missSec >= attrib->min) &&
3141 (attrib->missSec <= attrib->max)) {
3142 xmisss = attrib->max + 2;
3143 f_readjust = 1;
3144 }
3145 }
3146 }
3147 }
3148
3149 /* Walk through the grid, resetting the missing values, as determined by
3150 * the original grid. */
3151 if (f_readjust) {
3152 for (scanIndex = 0; scanIndex < (uInt4)nd2x3 && scanIndex < Nx * Ny; scanIndex++) {
3153 ScanIndex2XY (scanIndex, &x, &y, scan, Nx, Ny);
3154 /* ScanIndex returns value as if scan was 0100 */
3155 newIndex = (uInt4)(x - 1) + (uInt4)(y - 1) * Nx;
3156 if (attrib->fieldType) {
3157 value = iain[scanIndex];
3158 } else {
3159 value = ain[scanIndex];
3160 }
3161 if (value == attrib->missPri) {
3162 grib_Data[newIndex] = xmissp;
3163 } else if ((attrib->f_miss == 2) && (value == attrib->missSec)) {
3164 grib_Data[newIndex] = xmisss;
3165 }
3166 }
3167 attrib->missPri = xmissp;
3168 if (attrib->f_miss == 2) {
3169 attrib->missSec = xmisss;
3170 }
3171 }
3172
3173 /* Resolve bitmap (if there is one) in the data. */
3174 if (ibitmap) {
3175 attrib->f_maxmin = 0;
3176 if ((attrib->f_miss != 1) && (attrib->f_miss != 2)) {
3177 missCnt = 0;
3178 /* Figure out a missing value. */
3179 xmissp = 9999;
3180 #ifdef unreachable
3181 if (attrib->f_maxmin) {
3182 if ((xmissp <= attrib->max) && (xmissp >= attrib->min)) {
3183 xmissp = attrib->max + 1;
3184 }
3185 }
3186 #endif
3187 /* embed the missing value. */
3188 for (scanIndex = 0; scanIndex < (uInt4)nd2x3 && scanIndex < Nx * Ny; scanIndex++) {
3189 ScanIndex2XY (scanIndex, &x, &y, scan, Nx, Ny);
3190 /* ScanIndex returns value as if scan was 0100 */
3191 newIndex = (uInt4)(x - 1) + (uInt4)(y - 1) * Nx;
3192 /* Corrected this on 5/10/2004 */
3193 if (ib[scanIndex] != 1) {
3194 grib_Data[newIndex] = xmissp;
3195 missCnt++;
3196 } else {
3197 if (!attrib->f_maxmin) {
3198 attrib->f_maxmin = 1;
3199 attrib->max = attrib->min = grib_Data[newIndex];
3200 } else {
3201 if (attrib->max < grib_Data[newIndex])
3202 attrib->max = grib_Data[newIndex];
3203 if (attrib->min > grib_Data[newIndex])
3204 attrib->min = grib_Data[newIndex];
3205 }
3206 }
3207 }
3208 attrib->f_miss = 1;
3209 attrib->missPri = xmissp;
3210 }
3211 if (!attrib->f_maxmin) {
3212 attrib->f_maxmin = 1;
3213 attrib->max = attrib->min = xmissp;
3214 }
3215 }
3216 attrib->numMiss = missCnt;
3217 }
3218
3219 #ifdef unused_by_GDAL
3220 typedef struct {
3221 double value;
3222 int cnt;
3223 } freqType;
3224
freqCompare(const void * A,const void * B)3225 static int freqCompare (const void *A, const void *B)
3226 {
3227 const freqType *a = (freqType *) A;
3228 const freqType *b = (freqType *) B;
3229
3230 if (a->value < b->value)
3231 return -1;
3232 if (a->value > b->value)
3233 return 1;
3234 return 0;
3235 }
3236
FreqPrint(char ** ans,double * Data,sInt4 DataLen,sInt4 Nx,sInt4 Ny,sChar decimal,char * comment)3237 void FreqPrint (char **ans, double *Data, sInt4 DataLen, sInt4 Nx,
3238 sInt4 Ny, sChar decimal, char *comment)
3239 {
3240 int x, y, i;
3241 double *ptr = NULL;
3242 double value;
3243 freqType *freq = NULL;
3244 int numFreq = 0;
3245 char format[20];
3246
3247 myAssert (*ans == NULL);
3248
3249 if ((Nx < 0) || (Ny < 0) || (Nx * Ny > DataLen)) {
3250 return;
3251 }
3252
3253 ptr = Data;
3254 for (y = 0; y < Ny; y++) {
3255 for (x = 0; x < Nx; x++) {
3256 /* 2/28/2006 Introduced value to round before putting the data in
3257 * the Freq table. */
3258 value = myRound (*ptr, decimal);
3259 for (i = 0; i < numFreq; i++) {
3260 if (value == freq[i].value) {
3261 freq[i].cnt++;
3262 break;
3263 }
3264 }
3265 if (i == numFreq) {
3266 numFreq++;
3267 freq = (freqType *) realloc (freq, numFreq * sizeof (freqType));
3268 freq[i].value = value;
3269 freq[i].cnt = 1;
3270 }
3271 ptr++;
3272 }
3273 }
3274
3275 if( freq )
3276 qsort (freq, numFreq, sizeof (freq[0]), freqCompare);
3277
3278 mallocSprintf (ans, "%s | count\n", comment);
3279 snprintf (format, sizeof(format), "%%.%df | %%d\n", decimal);
3280 for (i = 0; i < numFreq; i++) {
3281 reallocSprintf (ans, format, myRound (freq[i].value, decimal),
3282 freq[i].cnt);
3283 }
3284 free (freq);
3285 }
3286 #endif
3287