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