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