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 6 : void MetaInit (grib_MetaData *meta)
53 : {
54 6 : meta->element = NULL;
55 6 : meta->comment = NULL;
56 6 : meta->unitName = NULL;
57 6 : meta->convert = 0;
58 6 : meta->shortFstLevel = NULL;
59 6 : meta->longFstLevel = NULL;
60 6 : meta->pds2.sect2.ptrType = GS2_NONE;
61 :
62 6 : meta->pds2.sect2.wx.data = NULL;
63 6 : meta->pds2.sect2.wx.dataLen = 0;
64 6 : meta->pds2.sect2.wx.maxLen = 0;
65 6 : meta->pds2.sect2.wx.ugly = NULL;
66 6 : meta->pds2.sect2.unknown.data = NULL;
67 6 : meta->pds2.sect2.unknown.dataLen = 0;
68 :
69 6 : meta->pds2.sect4.numInterval = 0;
70 6 : meta->pds2.sect4.Interval = NULL;
71 6 : meta->pds2.sect4.numBands = 0;
72 6 : 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 6 : void MetaSect2Free (grib_MetaData *meta)
99 : {
100 : size_t i; /* Counter for use when freeing Wx data. */
101 :
102 6 : 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 6 : free (meta->pds2.sect2.wx.ugly);
107 6 : meta->pds2.sect2.wx.ugly = NULL;
108 6 : free (meta->pds2.sect2.wx.data);
109 6 : meta->pds2.sect2.wx.data = NULL;
110 6 : meta->pds2.sect2.wx.dataLen = 0;
111 6 : meta->pds2.sect2.wx.maxLen = 0;
112 6 : meta->pds2.sect2.ptrType = GS2_NONE;
113 :
114 6 : free (meta->pds2.sect2.wx.data);
115 6 : meta->pds2.sect2.unknown.data = NULL;
116 6 : meta->pds2.sect2.unknown.dataLen = 0;
117 6 : }
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 6 : void MetaFree (grib_MetaData *meta)
141 : {
142 6 : free (meta->pds2.sect4.bands);
143 6 : meta->pds2.sect4.bands = NULL;
144 6 : meta->pds2.sect4.numBands = 0;
145 6 : free (meta->pds2.sect4.Interval);
146 6 : meta->pds2.sect4.Interval = NULL;
147 6 : meta->pds2.sect4.numInterval = 0;
148 6 : MetaSect2Free (meta);
149 6 : free (meta->unitName);
150 6 : meta->unitName = NULL;
151 6 : meta->convert = 0;
152 6 : free (meta->comment);
153 6 : meta->comment = NULL;
154 6 : free (meta->element);
155 6 : meta->element = NULL;
156 6 : free (meta->shortFstLevel);
157 6 : meta->shortFstLevel = NULL;
158 6 : free (meta->longFstLevel);
159 6 : meta->longFstLevel = NULL;
160 6 : }
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 24 : 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 24 : 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 24 : 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 24 : Clock_ScanDate (AnsTime, year, mon, day);
218 24 : *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 24 : 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 2 : static int ParseSect0 (sInt4 *is0, sInt4 ns0, sInt4 grib_len,
273 : grib_MetaData *meta)
274 : {
275 2 : if (ns0 < 9) {
276 0 : return -1;
277 : }
278 2 : 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 2 : meta->pds2.prodType = (uChar) is0[6];
285 2 : 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 2 : static int ParseSect1 (sInt4 *is1, sInt4 ns1, grib_MetaData *meta)
315 : {
316 2 : if (ns1 < 21) {
317 0 : return -1;
318 : }
319 2 : if (is1[4] != 1) {
320 0 : errSprintf ("ERROR IS1 not labeled correctly. %ld\n", is1[4]);
321 0 : return -2;
322 : }
323 2 : meta->center = (unsigned short int) is1[5];
324 2 : meta->subcenter = (unsigned short int) is1[7];
325 2 : meta->pds2.mstrVersion = (uChar) is1[9];
326 2 : meta->pds2.lclVersion = (uChar) is1[10];
327 2 : 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 2 : meta->pds2.sigTime = (uChar) is1[11];
341 4 : if (ParseTime (&(meta->pds2.refTime), is1[12], is1[14], is1[15], is1[16],
342 4 : is1[17], is1[18]) != 0) {
343 0 : preErrSprintf ("Error in call to ParseTime from ParseSect1 (GRIB2)");
344 0 : return -2;
345 : }
346 2 : meta->pds2.operStatus = (uChar) is1[19];
347 2 : meta->pds2.dataType = (uChar) is1[20];
348 2 : 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 2 : 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 2 : if (ns3 < 14) {
665 0 : return -1;
666 : }
667 2 : if (is3[4] != 3) {
668 0 : errSprintf ("ERROR IS3 not labeled correctly. %ld\n", is3[4]);
669 0 : return -2;
670 : }
671 2 : 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 2 : meta->gds.numPts = is3[6];
678 2 : 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 2 : 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 2 : if (ns3 < 38) {
698 0 : return -1;
699 : }
700 : /* Assert: is3[14] is the shape of the earth. */
701 2 : 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 2 : 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 2 : if ((is3[16] != GRIB2MISSING_s4) && (is3[15] != GRIB2MISSING_s1)) {
721 : /* Assumes data is given in m (not km). */
722 2 : meta->gds.majEarth = is3[16] / (pow (10.0, is3[15]) * 1000.);
723 2 : 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 2 : 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 2 : 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 2 : 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 2 : meta->gds.Nx = is3[30];
812 2 : meta->gds.Ny = is3[34];
813 2 : 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 2 : unit = 1e-6;
820 2 : meta->gds.center = 0;
821 2 : meta->gds.scaleLat1 = meta->gds.scaleLat2 = 0;
822 2 : meta->gds.southLat = meta->gds.southLon = 0;
823 2 : meta->gds.lat2 = meta->gds.lon2 = 0;
824 2 : 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 2 : if (ns3 < 72) {
871 0 : return -1;
872 : }
873 10 : if ((is3[38] == GRIB2MISSING_s4) || (is3[42] == GRIB2MISSING_s4) ||
874 4 : (is3[47] == GRIB2MISSING_s4) || (is3[51] == GRIB2MISSING_s4) ||
875 4 : (is3[55] == GRIB2MISSING_s4) || (is3[60] == GRIB2MISSING_s4)) {
876 0 : errSprintf ("Mercator grid is not defined completely.\n");
877 0 : return -2;
878 : }
879 2 : meta->gds.lat1 = is3[38] * unit;
880 2 : meta->gds.lon1 = is3[42] * unit;
881 2 : meta->gds.resFlag = (uChar) is3[46];
882 2 : meta->gds.meshLat = is3[47] * unit;
883 2 : meta->gds.lat2 = is3[51] * unit;
884 2 : meta->gds.lon2 = is3[55] * unit;
885 2 : meta->gds.scan = (uChar) is3[59];
886 2 : meta->gds.orientLon = is3[60] * unit;
887 2 : meta->gds.Dx = is3[64] / 1000.; /* mm -> m */
888 2 : meta->gds.Dy = is3[68] / 1000.; /* mm -> m */
889 : /* Resolve resolution flag(bit 3,4). Copy Dx,Dy as appropriate. */
890 2 : 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 2 : } 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 2 : 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 2 : if (meta->gds.scan != GRIB2BIT_2) {
990 : #ifdef DEBUG
991 : printf ("Scan mode is expected to be 0100 (ie %d) not %d\n",
992 : GRIB2BIT_2, meta->gds.scan);
993 : 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 : "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 2 : 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 44 : 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 44 : if ((unit >= 0) && (unit < 13)) {
1043 44 : if (unit2sec[unit] != 0) {
1044 44 : *ans = (double) (time * unit2sec[unit]);
1045 44 : 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 4 : 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 4 : if ((unit >= 0) && (unit < 14)) {
1092 4 : if (unit2sec[unit] != 0) {
1093 4 : *ans = (double) (time * unit2sec[unit]);
1094 4 : 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 2 : 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 2 : if (ns4 < 9) {
1146 0 : return -1;
1147 : }
1148 2 : if (is4[4] != 4) {
1149 : #ifdef DEBUG
1150 : 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 2 : if (is4[5] != 0) {
1156 : #ifdef DEBUG
1157 : printf ("Un-supported template.\n All Supported template "
1158 : "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 8 : if ((is4[7] != GS4_ANALYSIS) && (is4[7] != GS4_ENSEMBLE) &&
1165 4 : (is4[7] != GS4_DERIVED) && (is4[7] != GS4_PROBABIL_PNT) &&
1166 2 : (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 : 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 2 : meta->pds2.sect4.templat = (unsigned short int) is4[7];
1176 :
1177 : /*
1178 : * Handle variables common to the supported templates.
1179 : */
1180 2 : if (ns4 < 34) {
1181 0 : return -1;
1182 : }
1183 2 : meta->pds2.sect4.cat = (uChar) is4[9];
1184 2 : meta->pds2.sect4.subcat = (uChar) is4[10];
1185 2 : meta->pds2.sect4.genProcess = (uChar) is4[11];
1186 :
1187 : /* Initialize variables prior to parsing the specific templates. */
1188 2 : meta->pds2.sect4.typeEnsemble = 0;
1189 2 : meta->pds2.sect4.perturbNum = 0;
1190 2 : meta->pds2.sect4.numberFcsts = 0;
1191 2 : meta->pds2.sect4.derivedFcst = 0;
1192 2 : meta->pds2.sect4.validTime = meta->pds2.refTime;
1193 :
1194 2 : 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 2 : meta->pds2.sect4.bgGenID = (uChar) is4[12];
1222 2 : meta->pds2.sect4.genID = (uChar) is4[13];
1223 4 : if ((is4[14] == GRIB2MISSING_u2) || (is4[16] == GRIB2MISSING_u1)) {
1224 2 : meta->pds2.sect4.f_validCutOff = 0;
1225 2 : 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 2 : if (is4[18] == GRIB2MISSING_s4) {
1231 0 : errSprintf ("Missing 'forecast' time?\n");
1232 0 : return -5;
1233 : }
1234 2 : 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 2 : 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 2 : meta->pds2.sect4.fstSurfType = (uChar) is4[22];
1248 2 : 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 2 : meta->pds2.sect4.fstSurfScale = is4[23];
1254 2 : meta->pds2.sect4.fstSurfValue = is4[24] / pow (10.0, is4[23]);
1255 : }
1256 2 : meta->pds2.sect4.sndSurfType = (uChar) is4[28];
1257 4 : if ((is4[30] == GRIB2MISSING_s4) || (is4[29] == GRIB2MISSING_s1) ||
1258 : (meta->pds2.sect4.sndSurfType == GRIB2MISSING_u1)) {
1259 2 : meta->pds2.sect4.sndSurfScale = GRIB2MISSING_s1;
1260 2 : 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 2 : 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 : meta->pds2.sect4.numInterval = (uChar) is4[44];
1281 0 : if (meta->pds2.sect4.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 : printf ("Warning: in call to ParseTime from ParseSect4\n%s", msg);
1290 0 : free (msg);
1291 : meta->pds2.sect4.validTime = (time_t) (meta->pds2.refTime +
1292 0 : meta->pds2.sect4.foreSec);
1293 : printf ("Most likely they didn't complete bytes 38-44 of "
1294 0 : "Template 4.11\n");
1295 : } else {
1296 0 : meta->pds2.sect4.numInterval = (uChar) is4[44];
1297 : }
1298 :
1299 : /* Added this check because some MOS grids didn't finish the
1300 : * template. */
1301 0 : if (meta->pds2.sect4.numInterval != 0) {
1302 : temp_ptr = realloc ((void *) meta->pds2.sect4.Interval,
1303 : meta->pds2.sect4.numInterval *
1304 0 : sizeof (sect4_IntervalType));
1305 0 : if (temp_ptr == NULL) {
1306 0 : printf ("Ran out of memory.\n");
1307 0 : return -6;
1308 : }
1309 0 : meta->pds2.sect4.Interval = (sect4_IntervalType *) temp_ptr;
1310 0 : meta->pds2.sect4.numMissing = is4[45];
1311 0 : for (i = 0; i < meta->pds2.sect4.numInterval; i++) {
1312 0 : meta->pds2.sect4.Interval[i].processID =
1313 0 : (uChar) is4[49 + i * 12];
1314 0 : meta->pds2.sect4.Interval[i].incrType =
1315 0 : (uChar) is4[50 + i * 12];
1316 0 : meta->pds2.sect4.Interval[i].timeRangeUnit =
1317 0 : (uChar) is4[51 + i * 12];
1318 0 : meta->pds2.sect4.Interval[i].lenTime = is4[52 + i * 12];
1319 0 : meta->pds2.sect4.Interval[i].incrUnit =
1320 0 : (uChar) is4[56 + i * 12];
1321 0 : meta->pds2.sect4.Interval[i].timeIncr =
1322 0 : (uChar) is4[57 + i * 12];
1323 : }
1324 : } else {
1325 : #ifdef DEBUG
1326 : printf ("Caution: Template 4.11 had no Intervals.\n");
1327 : #endif
1328 0 : meta->pds2.sect4.numMissing = is4[45];
1329 : }
1330 0 : break;
1331 : case GS4_DERIVED: /* 4.2 */
1332 0 : meta->pds2.sect4.derivedFcst = (uChar) is4[34];
1333 0 : meta->pds2.sect4.numberFcsts = (uChar) is4[35];
1334 0 : break;
1335 : case GS4_DERIVED_INTERVAL: /* 4.12 */
1336 0 : meta->pds2.sect4.derivedFcst = (uChar) is4[34];
1337 0 : meta->pds2.sect4.numberFcsts = (uChar) is4[35];
1338 :
1339 0 : if (ParseTime (&(meta->pds2.sect4.validTime), is4[36], is4[38],
1340 0 : is4[39], is4[40], is4[41], is4[42]) != 0) {
1341 0 : msg = errSprintf (NULL);
1342 0 : meta->pds2.sect4.numInterval = (uChar) is4[43];
1343 0 : if (meta->pds2.sect4.numInterval != 1) {
1344 : errSprintf ("ERROR: in call to ParseTime from ParseSect4\n%s",
1345 0 : msg);
1346 : errSprintf ("Most likely they didn't complete bytes 37-43 of "
1347 0 : "Template 4.12\n");
1348 0 : free (msg);
1349 0 : return -1;
1350 : }
1351 0 : printf ("Warning: in call to ParseTime from ParseSect4\n%s", msg);
1352 0 : free (msg);
1353 : meta->pds2.sect4.validTime = (time_t) (meta->pds2.refTime +
1354 0 : meta->pds2.sect4.foreSec);
1355 : printf ("Most likely they didn't complete bytes 37-43 of "
1356 0 : "Template 4.12\n");
1357 : } else {
1358 0 : meta->pds2.sect4.numInterval = (uChar) is4[43];
1359 : }
1360 :
1361 : /* Added this check because some MOS grids didn't finish the
1362 : * template. */
1363 0 : if (meta->pds2.sect4.numInterval != 0) {
1364 : temp_ptr = realloc ((void *) meta->pds2.sect4.Interval,
1365 : meta->pds2.sect4.numInterval *
1366 0 : sizeof (sect4_IntervalType));
1367 0 : if (temp_ptr == NULL) {
1368 0 : printf ("Ran out of memory.\n");
1369 0 : return -6;
1370 : }
1371 0 : meta->pds2.sect4.Interval = (sect4_IntervalType *) temp_ptr;
1372 0 : meta->pds2.sect4.numMissing = is4[44];
1373 0 : for (i = 0; i < meta->pds2.sect4.numInterval; i++) {
1374 0 : meta->pds2.sect4.Interval[i].processID =
1375 0 : (uChar) is4[48 + i * 12];
1376 0 : meta->pds2.sect4.Interval[i].incrType =
1377 0 : (uChar) is4[49 + i * 12];
1378 0 : meta->pds2.sect4.Interval[i].timeRangeUnit =
1379 0 : (uChar) is4[50 + i * 12];
1380 0 : meta->pds2.sect4.Interval[i].lenTime = is4[51 + i * 12];
1381 0 : meta->pds2.sect4.Interval[i].incrUnit =
1382 0 : (uChar) is4[55 + i * 12];
1383 0 : meta->pds2.sect4.Interval[i].timeIncr =
1384 0 : (uChar) is4[56 + i * 12];
1385 : }
1386 : } else {
1387 : #ifdef DEBUG
1388 : printf ("Caution: Template 4.12 had no Intervals.\n");
1389 : #endif
1390 0 : meta->pds2.sect4.numMissing = is4[44];
1391 : }
1392 0 : break;
1393 : case GS4_STATISTIC: /* 4.8 */
1394 8 : if (ParseTime (&(meta->pds2.sect4.validTime), is4[34], is4[36],
1395 8 : is4[37], is4[38], is4[39], is4[40]) != 0) {
1396 0 : msg = errSprintf (NULL);
1397 0 : meta->pds2.sect4.numInterval = (uChar) is4[41];
1398 0 : if (meta->pds2.sect4.numInterval != 1) {
1399 : errSprintf ("ERROR: in call to ParseTime from ParseSect4\n%s",
1400 0 : msg);
1401 : errSprintf ("Most likely they didn't complete bytes 35-41 of "
1402 0 : "Template 4.8\n");
1403 0 : free (msg);
1404 0 : return -1;
1405 : }
1406 0 : printf ("Warning: in call to ParseTime from ParseSect4\n%s", msg);
1407 0 : free (msg);
1408 : meta->pds2.sect4.validTime = (time_t) (meta->pds2.refTime +
1409 0 : meta->pds2.sect4.foreSec);
1410 : printf ("Most likely they didn't complete bytes 35-41 of "
1411 0 : "Template 4.8\n");
1412 : } else {
1413 2 : meta->pds2.sect4.numInterval = (uChar) is4[41];
1414 : }
1415 :
1416 : /* Added this check because some MOS grids didn't finish the
1417 : * template. */
1418 2 : if (meta->pds2.sect4.numInterval != 0) {
1419 : temp_ptr = realloc ((void *) meta->pds2.sect4.Interval,
1420 : meta->pds2.sect4.numInterval *
1421 2 : sizeof (sect4_IntervalType));
1422 2 : if (temp_ptr == NULL) {
1423 0 : printf ("Ran out of memory.\n");
1424 0 : return -6;
1425 : }
1426 2 : meta->pds2.sect4.Interval = (sect4_IntervalType *) temp_ptr;
1427 2 : meta->pds2.sect4.numMissing = is4[42];
1428 4 : for (i = 0; i < meta->pds2.sect4.numInterval; i++) {
1429 2 : meta->pds2.sect4.Interval[i].processID =
1430 2 : (uChar) is4[46 + i * 12];
1431 2 : meta->pds2.sect4.Interval[i].incrType =
1432 2 : (uChar) is4[47 + i * 12];
1433 2 : meta->pds2.sect4.Interval[i].timeRangeUnit =
1434 2 : (uChar) is4[48 + i * 12];
1435 2 : meta->pds2.sect4.Interval[i].lenTime = is4[49 + i * 12];
1436 2 : meta->pds2.sect4.Interval[i].incrUnit =
1437 2 : (uChar) is4[53 + i * 12];
1438 2 : meta->pds2.sect4.Interval[i].timeIncr =
1439 2 : (uChar) is4[54 + i * 12];
1440 : }
1441 : } else {
1442 : #ifdef DEBUG
1443 : printf ("Caution: Template 4.8 had no Intervals.\n");
1444 : #endif
1445 0 : meta->pds2.sect4.numMissing = is4[42];
1446 : }
1447 2 : break;
1448 : case GS4_PERCENTILE: /* 4.10 */
1449 0 : meta->pds2.sect4.percentile = is4[34];
1450 0 : if (ParseTime (&(meta->pds2.sect4.validTime), is4[35], is4[37],
1451 0 : is4[38], is4[39], is4[40], is4[41]) != 0) {
1452 0 : msg = errSprintf (NULL);
1453 0 : meta->pds2.sect4.numInterval = (uChar) is4[42];
1454 0 : if (meta->pds2.sect4.numInterval != 1) {
1455 : errSprintf ("ERROR: in call to ParseTime from ParseSect4\n%s",
1456 0 : msg);
1457 : errSprintf ("Most likely they didn't complete bytes 35-41 of "
1458 0 : "Template 4.8\n");
1459 0 : free (msg);
1460 0 : return -1;
1461 : }
1462 0 : printf ("Warning: in call to ParseTime from ParseSect4\n%s", msg);
1463 0 : free (msg);
1464 : meta->pds2.sect4.validTime = (time_t) (meta->pds2.refTime +
1465 0 : meta->pds2.sect4.foreSec);
1466 : printf ("Most likely they didn't complete bytes 35-41 of "
1467 0 : "Template 4.8\n");
1468 : } else {
1469 0 : meta->pds2.sect4.numInterval = (uChar) is4[42];
1470 : }
1471 :
1472 : /* Added this check because some MOS grids didn't finish the
1473 : * template. */
1474 0 : if (meta->pds2.sect4.numInterval != 0) {
1475 : temp_ptr = realloc ((void *) meta->pds2.sect4.Interval,
1476 : meta->pds2.sect4.numInterval *
1477 0 : sizeof (sect4_IntervalType));
1478 0 : if (temp_ptr == NULL) {
1479 0 : printf ("Ran out of memory.\n");
1480 0 : return -6;
1481 : }
1482 0 : meta->pds2.sect4.Interval = (sect4_IntervalType *) temp_ptr;
1483 0 : meta->pds2.sect4.numMissing = is4[43];
1484 0 : for (i = 0; i < meta->pds2.sect4.numInterval; i++) {
1485 0 : meta->pds2.sect4.Interval[i].processID =
1486 0 : (uChar) is4[47 + i * 12];
1487 0 : meta->pds2.sect4.Interval[i].incrType =
1488 0 : (uChar) is4[48 + i * 12];
1489 0 : meta->pds2.sect4.Interval[i].timeRangeUnit =
1490 0 : (uChar) is4[49 + i * 12];
1491 0 : meta->pds2.sect4.Interval[i].lenTime = is4[50 + i * 12];
1492 0 : meta->pds2.sect4.Interval[i].incrUnit =
1493 0 : (uChar) is4[54 + i * 12];
1494 0 : meta->pds2.sect4.Interval[i].timeIncr =
1495 0 : (uChar) is4[55 + i * 12];
1496 : }
1497 : } else {
1498 : #ifdef DEBUG
1499 : printf ("Caution: Template 4.10 had no Intervals.\n");
1500 : #endif
1501 0 : meta->pds2.sect4.numMissing = is4[43];
1502 : }
1503 0 : break;
1504 : case GS4_PROBABIL_PNT: /* 4.5 */
1505 0 : meta->pds2.sect4.foreProbNum = (uChar) is4[34];
1506 0 : meta->pds2.sect4.numForeProbs = (uChar) is4[35];
1507 0 : meta->pds2.sect4.probType = (uChar) is4[36];
1508 0 : meta->pds2.sect4.lowerLimit.factor = (sChar) is4[37];
1509 0 : meta->pds2.sect4.lowerLimit.value = is4[38];
1510 0 : meta->pds2.sect4.upperLimit.factor = (sChar) is4[42];
1511 0 : meta->pds2.sect4.upperLimit.value = is4[43];
1512 0 : break;
1513 : case GS4_PROBABIL_TIME: /* 4.9 */
1514 0 : meta->pds2.sect4.foreProbNum = (uChar) is4[34];
1515 0 : meta->pds2.sect4.numForeProbs = (uChar) is4[35];
1516 0 : meta->pds2.sect4.probType = (uChar) is4[36];
1517 0 : meta->pds2.sect4.lowerLimit.factor = (sChar) is4[37];
1518 0 : meta->pds2.sect4.lowerLimit.value = is4[38];
1519 0 : meta->pds2.sect4.upperLimit.factor = (sChar) is4[42];
1520 0 : meta->pds2.sect4.upperLimit.value = is4[43];
1521 0 : if (ParseTime (&(meta->pds2.sect4.validTime), is4[47], is4[49],
1522 0 : is4[50], is4[51], is4[52], is4[53]) != 0) {
1523 0 : msg = errSprintf (NULL);
1524 0 : meta->pds2.sect4.numInterval = (uChar) is4[54];
1525 0 : if (meta->pds2.sect4.numInterval != 1) {
1526 : errSprintf ("ERROR: in call to ParseTime from ParseSect4\n%s",
1527 0 : msg);
1528 : errSprintf ("Most likely they didn't complete bytes 48-54 of "
1529 0 : "Template 4.9\n");
1530 0 : free (msg);
1531 0 : return -1;
1532 : }
1533 0 : printf ("Warning: in call to ParseTime from ParseSect4\n%s", msg);
1534 0 : free (msg);
1535 : meta->pds2.sect4.validTime = (time_t) (meta->pds2.refTime +
1536 0 : meta->pds2.sect4.foreSec);
1537 : printf ("Most likely they didn't complete bytes 48-54 of "
1538 0 : "Template 4.9\n");
1539 : } else {
1540 0 : meta->pds2.sect4.numInterval = (uChar) is4[54];
1541 : }
1542 : temp_ptr = realloc ((void *) meta->pds2.sect4.Interval,
1543 : meta->pds2.sect4.numInterval *
1544 0 : sizeof (sect4_IntervalType));
1545 0 : if (temp_ptr == NULL) {
1546 0 : printf ("Ran out of memory.\n");
1547 0 : return -6;
1548 : }
1549 0 : meta->pds2.sect4.Interval = (sect4_IntervalType *) temp_ptr;
1550 0 : meta->pds2.sect4.numMissing = is4[55];
1551 0 : for (i = 0; i < meta->pds2.sect4.numInterval; i++) {
1552 0 : meta->pds2.sect4.Interval[i].processID = (uChar) is4[59 + i * 12];
1553 0 : meta->pds2.sect4.Interval[i].incrType = (uChar) is4[60 + i * 12];
1554 0 : meta->pds2.sect4.Interval[i].timeRangeUnit =
1555 0 : (uChar) is4[61 + i * 12];
1556 0 : meta->pds2.sect4.Interval[i].lenTime = is4[62 + i * 12];
1557 0 : meta->pds2.sect4.Interval[i].incrUnit = (uChar) is4[66 + i * 12];
1558 0 : meta->pds2.sect4.Interval[i].timeIncr = (uChar) is4[67 + i * 12];
1559 : }
1560 0 : break;
1561 : default:
1562 0 : errSprintf ("Un-supported Template. %ld\n", is4[7]);
1563 0 : return -4;
1564 : }
1565 2 : return 0;
1566 : }
1567 :
1568 : /*****************************************************************************
1569 : * ParseSect5() --
1570 : *
1571 : * Arthur Taylor / MDL
1572 : *
1573 : * PURPOSE
1574 : * To verify and parse section 5 data.
1575 : *
1576 : * ARGUMENTS
1577 : * is5 = The unpacked section 5 array. (Input)
1578 : * ns5 = The size of section 5. (Input)
1579 : * meta = The structure to fill. (Output)
1580 : * xmissp = The primary missing value. (Input)
1581 : * xmisss = The secondary missing value. (Input)
1582 : *
1583 : * FILES/DATABASES: None
1584 : *
1585 : * RETURNS: int (could use errSprintf())
1586 : * 0 = OK
1587 : * -1 = ns5 is too small.
1588 : * -2 = unexpected values in is5.
1589 : * -6 = unsupported packing.
1590 : *
1591 : * HISTORY
1592 : * 9/2002 Arthur Taylor (MDL/RSIS): Created.
1593 : *
1594 : * NOTES
1595 : *****************************************************************************
1596 : */
1597 2 : static int ParseSect5 (sInt4 *is5, sInt4 ns5, grib_MetaData *meta,
1598 : float xmissp, float xmisss)
1599 : {
1600 2 : if (ns5 < 22) {
1601 0 : return -1;
1602 : }
1603 2 : if (is5[4] != 5) {
1604 0 : errSprintf ("ERROR IS5 not labeled correctly. %ld\n", is5[5]);
1605 0 : return -2;
1606 : }
1607 4 : if ((is5[9] != GS5_SIMPLE) && (is5[9] != GS5_CMPLX) &&
1608 2 : (is5[9] != GS5_CMPLXSEC) && (is5[9] != GS5_SPECTRAL) &&
1609 0 : (is5[9] != GS5_HARMONIC) && (is5[9] != GS5_JPEG2000) &&
1610 0 : (is5[9] != GS5_PNG) && (is5[9] != GS5_JPEG2000_ORG) &&
1611 0 : (is5[9] != GS5_PNG_ORG)) {
1612 0 : errSprintf ("Un-supported Packing? %ld\n", is5[9]);
1613 0 : return -6;
1614 : }
1615 2 : meta->gridAttrib.packType = (sInt4) is5[9];
1616 2 : meta->gridAttrib.f_maxmin = 0;
1617 2 : meta->gridAttrib.missPri = xmissp;
1618 2 : meta->gridAttrib.missSec = xmisss;
1619 2 : if ((is5[9] == GS5_SPECTRAL) || (is5[9] == GS5_HARMONIC)) {
1620 0 : meta->gridAttrib.fieldType = 0;
1621 0 : meta->gridAttrib.f_miss = 0;
1622 0 : return 0;
1623 : }
1624 2 : if (is5[20] > 1) {
1625 0 : errSprintf ("Invalid field type. %ld\n", is5[20]);
1626 0 : return -2;
1627 : }
1628 2 : MEMCPY_BIG (&meta->gridAttrib.refVal, &(is5[11]), 4);
1629 2 : meta->gridAttrib.ESF = is5[15];
1630 2 : meta->gridAttrib.DSF = is5[17];
1631 2 : meta->gridAttrib.fieldType = (uChar) is5[20];
1632 6 : if ((is5[9] == GS5_JPEG2000) || (is5[9] == GS5_JPEG2000_ORG) ||
1633 4 : (is5[9] == GS5_PNG) || (is5[9] == GS5_PNG_ORG)) {
1634 0 : meta->gridAttrib.f_miss = 0;
1635 0 : return 0;
1636 : }
1637 2 : if (meta->gridAttrib.packType == 0) {
1638 0 : meta->gridAttrib.f_miss = 0;
1639 : } else {
1640 2 : if (ns5 < 23) {
1641 0 : return -1;
1642 : }
1643 2 : if (is5[22] > 2) {
1644 : errSprintf ("Invalid missing management type, f_miss = %ld\n",
1645 0 : is5[22]);
1646 0 : return -2;
1647 : }
1648 2 : meta->gridAttrib.f_miss = (uChar) is5[22];
1649 : }
1650 2 : return 0;
1651 : }
1652 :
1653 : /*****************************************************************************
1654 : * MetaParse() --
1655 : *
1656 : * Arthur Taylor / MDL
1657 : *
1658 : * PURPOSE
1659 : * To parse all the meta data from a grib2 message.
1660 : *
1661 : * ARGUMENTS
1662 : * meta = The structure to fill. (Output)
1663 : * is0 = The unpacked section 0 array. (Input)
1664 : * ns0 = The size of section 0. (Input)
1665 : * is1 = The unpacked section 1 array. (Input)
1666 : * ns1 = The size of section 1. (Input)
1667 : * is2 = The unpacked section 2 array. (Input)
1668 : * ns2 = The size of section 2. (Input)
1669 : * rdat = The float data in section 2. (Input)
1670 : * nrdat = Length of rdat. (Input)
1671 : * idat = The integer data in section 2. (Input)
1672 : * nidat = Length of idat. (Input)
1673 : * is3 = The unpacked section 3 array. (Input)
1674 : * ns3 = The size of section 3. (Input)
1675 : * is4 = The unpacked section 4 array. (Input)
1676 : * ns4 = The size of section 4. (Input)
1677 : * is5 = The unpacked section 5 array. (Input)
1678 : * ns5 = The size of section 5. (Input)
1679 : * grib_len = The length of the entire grib message. (Input)
1680 : * xmissp = The primary missing value. (Input)
1681 : * xmisss = The secondary missing value. (Input)
1682 : * simpVer = The version of the simple weather code to use when parsing the
1683 : * WxString (if applicable). (Input)
1684 : *
1685 : * FILES/DATABASES: None
1686 : *
1687 : * RETURNS: int (could use errSprintf())
1688 : * 0 = OK
1689 : * -1 = A dimension is too small.
1690 : * -2 = unexpected values in a grib section.
1691 : * -3 = un-supported map Projection.
1692 : * -4 = un-supported Sect 4 template.
1693 : * -5 = unsupported forecast time unit.
1694 : * -6 = unsupported sect 5 packing.
1695 : * -10 = Something the driver can't handle yet.
1696 : * (prodType != 0, f_sphere != 1, etc)
1697 : * -11 = Weather grid without a lookup table.
1698 : *
1699 : * HISTORY
1700 : * 9/2002 Arthur Taylor (MDL/RSIS): Created.
1701 : *
1702 : * NOTES
1703 : *****************************************************************************
1704 : */
1705 2 : int MetaParse (grib_MetaData *meta, sInt4 *is0, sInt4 ns0,
1706 : sInt4 *is1, sInt4 ns1, sInt4 *is2, sInt4 ns2,
1707 : float *rdat, sInt4 nrdat, sInt4 *idat, sInt4 nidat,
1708 : sInt4 *is3, sInt4 ns3, sInt4 *is4, sInt4 ns4,
1709 : sInt4 *is5, sInt4 ns5, sInt4 grib_len,
1710 : float xmissp, float xmisss, int simpVer)
1711 : {
1712 : int ierr; /* The error code of a called routine */
1713 : /* char *element; *//* Holds the name of the current variable. */
1714 : /* char *comment; *//* Holds more comments about current variable. */
1715 : /* char *unitName; *//* Holds the name of the unit [K] [%] .. etc */
1716 : uChar probType; /* The probability type */
1717 : double lowerProb; /* The lower limit on probability forecast if
1718 : * template 4.5 or 4.9 */
1719 : double upperProb; /* The upper limit on probability forecast if
1720 : * template 4.5 or 4.9 */
1721 : sInt4 lenTime; /* Length of time for element (see 4.8 and 4.9) */
1722 :
1723 2 : if ((ierr = ParseSect0 (is0, ns0, grib_len, meta)) != 0) {
1724 0 : preErrSprintf ("Parse error Section 0\n");
1725 : //return ierr;
1726 : }
1727 2 : if ((ierr = ParseSect1 (is1, ns1, meta)) != 0) {
1728 0 : preErrSprintf ("Parse error Section 1\n");
1729 : //return ierr;
1730 : }
1731 2 : if (ns2 < 7) {
1732 0 : errSprintf ("ns2 was too small in MetaParse\n");
1733 : //return -1;
1734 : }
1735 2 : meta->pds2.f_sect2 = (uChar) (is2[0] != 0);
1736 2 : if (meta->pds2.f_sect2) {
1737 0 : meta->pds2.sect2NumGroups = is2[7 - 1];
1738 : } else {
1739 2 : meta->pds2.sect2NumGroups = 0;
1740 : }
1741 2 : if ((ierr = ParseSect3 (is3, ns3, meta)) != 0) {
1742 0 : preErrSprintf ("Parse error Section 3\n");
1743 : //return ierr;
1744 : }
1745 2 : if (meta->gds.f_sphere != 1) {
1746 0 : errSprintf ("Driver Filter: Can only handle spheres.\n");
1747 : //return -10;
1748 : }
1749 2 : if ((ierr = ParseSect4 (is4, ns4, meta)) != 0) {
1750 0 : preErrSprintf ("Parse error Section 4\n");
1751 : //return ierr;
1752 : }
1753 2 : if ((ierr = ParseSect5 (is5, ns5, meta, xmissp, xmisss)) != 0) {
1754 0 : preErrSprintf ("Parse error Section 5\n");
1755 : //return ierr;
1756 : }
1757 : /* Compute ElementName. */
1758 2 : if (meta->element) {
1759 0 : free (meta->element);
1760 0 : meta->element = NULL;
1761 : }
1762 2 : if (meta->unitName) {
1763 0 : free (meta->unitName);
1764 0 : meta->unitName = NULL;
1765 : }
1766 2 : if (meta->comment) {
1767 0 : free (meta->comment);
1768 0 : meta->comment = NULL;
1769 : }
1770 :
1771 2 : if ((meta->pds2.sect4.templat == GS4_PROBABIL_TIME) ||
1772 : (meta->pds2.sect4.templat == GS4_PROBABIL_PNT)) {
1773 0 : probType = meta->pds2.sect4.probType;
1774 : lowerProb = meta->pds2.sect4.lowerLimit.value *
1775 0 : pow (10.0, -1 * meta->pds2.sect4.lowerLimit.factor);
1776 : upperProb = meta->pds2.sect4.upperLimit.value *
1777 0 : pow (10.0, -1 * meta->pds2.sect4.upperLimit.factor);
1778 : } else {
1779 2 : probType = 0;
1780 2 : lowerProb = 0;
1781 2 : upperProb = 0;
1782 : }
1783 2 : if (meta->pds2.sect4.numInterval > 0) {
1784 : /* Try to convert lenTime to hourly. */
1785 2 : if (meta->pds2.sect4.Interval[0].timeRangeUnit == 255) {
1786 : lenTime = (sInt4) ((meta->pds2.sect4.validTime -
1787 : meta->pds2.sect4.foreSec -
1788 0 : meta->pds2.refTime) / 3600);
1789 2 : } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 0) {
1790 0 : lenTime = (sInt4) (meta->pds2.sect4.Interval[0].lenTime / 60.);
1791 2 : } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 1) {
1792 2 : lenTime = meta->pds2.sect4.Interval[0].lenTime;
1793 0 : } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 2) {
1794 0 : lenTime = meta->pds2.sect4.Interval[0].lenTime * 24;
1795 0 : } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 10) {
1796 0 : lenTime = meta->pds2.sect4.Interval[0].lenTime * 3;
1797 0 : } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 11) {
1798 0 : lenTime = meta->pds2.sect4.Interval[0].lenTime * 6;
1799 0 : } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 12) {
1800 0 : lenTime = meta->pds2.sect4.Interval[0].lenTime * 12;
1801 0 : } else if (meta->pds2.sect4.Interval[0].timeRangeUnit == 13) {
1802 0 : lenTime = (sInt4) (meta->pds2.sect4.Interval[0].lenTime / 3600.);
1803 : } else {
1804 0 : lenTime = 0;
1805 0 : printf ("Can't handle this timeRangeUnit\n");
1806 : myAssert (meta->pds2.sect4.Interval[0].timeRangeUnit == 1);
1807 : }
1808 : /*
1809 : } else {
1810 : lenTime = 255;
1811 : }
1812 : if (lenTime == 255) {
1813 : lenTime = (meta->pds2.sect4.validTime - meta->pds2.sect4.foreSec -
1814 : meta->pds2.refTime) / 3600;
1815 : }
1816 : */
1817 2 : if (lenTime == GRIB2MISSING_s4) {
1818 0 : lenTime = 0;
1819 : }
1820 : ParseElemName (meta->center, meta->subcenter,
1821 : meta->pds2.prodType, meta->pds2.sect4.templat,
1822 : meta->pds2.sect4.cat, meta->pds2.sect4.subcat,
1823 2 : lenTime, meta->pds2.sect4.Interval[0].incrType,
1824 : meta->pds2.sect4.genID, probType, lowerProb,
1825 : upperProb, &(meta->element), &(meta->comment),
1826 : &(meta->unitName), &(meta->convert),
1827 4 : meta->pds2.sect4.percentile);
1828 : } else {
1829 : ParseElemName (meta->center, meta->subcenter,
1830 : meta->pds2.prodType, meta->pds2.sect4.templat,
1831 : meta->pds2.sect4.cat, meta->pds2.sect4.subcat, 0, 255,
1832 : meta->pds2.sect4.genID, probType, lowerProb, upperProb,
1833 : &(meta->element), &(meta->comment), &(meta->unitName),
1834 0 : &(meta->convert), meta->pds2.sect4.percentile);
1835 : }
1836 : #ifdef DEBUG
1837 : /*
1838 : printf ("Element: %s\nunitName: %s\ncomment: %s\n", meta->element,
1839 : meta->comment, meta->unitName);
1840 : */
1841 : #endif
1842 :
1843 : /*
1844 : if (strcmp (element, "") == 0) {
1845 : meta->element = (char *) realloc ((void *) (meta->element),
1846 : (1 + strlen ("unknown")) *
1847 : sizeof (char));
1848 : strcpy (meta->element, "unknown");
1849 : } else {
1850 : if (IsData_MOS (meta->pds2.center, meta->pds2.subcenter)) {
1851 : * See : http://www.nco.ncep.noaa.gov/pmb/docs/on388/tablea.html *
1852 : if (meta->pds2.sect4.genID == 96) {
1853 : meta->element = (char *) realloc ((void *) (meta->element),
1854 : (1 + 7 + strlen (element)) *
1855 : sizeof (char));
1856 : sprintf (meta->element, "MOSGFS-%s", element);
1857 : } else {
1858 : meta->element = (char *) realloc ((void *) (meta->element),
1859 : (1 + 4 + strlen (element)) *
1860 : sizeof (char));
1861 : sprintf (meta->element, "MOS-%s", element);
1862 : }
1863 : } else {
1864 : meta->element = (char *) realloc ((void *) (meta->element),
1865 : (1 + strlen (element)) *
1866 : sizeof (char));
1867 : strcpy (meta->element, element);
1868 : }
1869 : }
1870 : meta->unitName = (char *) realloc ((void *) (meta->unitName),
1871 : (1 + 2 + strlen (unitName)) *
1872 : sizeof (char));
1873 : sprintf (meta->unitName, "[%s]", unitName);
1874 : meta->comment = (char *) realloc ((void *) (meta->comment),
1875 : (1 + strlen (comment) +
1876 : strlen (unitName)
1877 : + 2 + 1) * sizeof (char));
1878 : sprintf (meta->comment, "%s [%s]", comment, unitName);
1879 : */
1880 4 : if ((meta->pds2.sect4.sndSurfScale == GRIB2MISSING_s1) ||
1881 : (meta->pds2.sect4.sndSurfType == GRIB2MISSING_u1)) {
1882 : /*
1883 : if ((meta->pds2.sect4.fstSurfScale == GRIB2MISSING_s1) ||
1884 : (meta->pds2.sect4.fstSurfType == GRIB2MISSING_u1)) {
1885 : ParseLevelName (meta->center, meta->subcenter,
1886 : meta->pds2.sect4.fstSurfType, 0, 0, 0,
1887 : &(meta->shortFstLevel), &(meta->longFstLevel));
1888 : } else {
1889 : */
1890 : ParseLevelName (meta->center, meta->subcenter,
1891 : meta->pds2.sect4.fstSurfType,
1892 : meta->pds2.sect4.fstSurfValue, 0, 0,
1893 2 : &(meta->shortFstLevel), &(meta->longFstLevel));
1894 : /*
1895 : }
1896 : */
1897 : } else {
1898 : ParseLevelName (meta->center, meta->subcenter,
1899 : meta->pds2.sect4.fstSurfType,
1900 : meta->pds2.sect4.fstSurfValue, 1,
1901 : meta->pds2.sect4.sndSurfValue, &(meta->shortFstLevel),
1902 0 : &(meta->longFstLevel));
1903 : }
1904 :
1905 : /* Continue parsing section 2 data. */
1906 2 : if (meta->pds2.f_sect2) {
1907 0 : MetaSect2Free (meta);
1908 0 : if (strcmp (meta->element, "Wx") == 0) {
1909 0 : meta->pds2.sect2.ptrType = GS2_WXTYPE;
1910 0 : if ((ierr = ParseSect2_Wx (rdat, nrdat, idat, nidat,
1911 : &(meta->pds2.sect2.wx), simpVer)) != 0) {
1912 0 : preErrSprintf ("Parse error Section 2 : Weather Data\n");
1913 : //return ierr;
1914 : }
1915 : } else {
1916 0 : meta->pds2.sect2.ptrType = GS2_UNKNOWN;
1917 0 : if ((ierr = ParseSect2_Unknown (rdat, nrdat, idat, nidat, meta))
1918 : != 0) {
1919 0 : preErrSprintf ("Parse error Section 2 : Unknown Data type\n");
1920 : //return ierr;
1921 : }
1922 : }
1923 : } else {
1924 2 : if (strcmp (meta->element, "Wx") == 0) {
1925 0 : errSprintf ("Weather grid does not have look up table?");
1926 : //return -11;
1927 : }
1928 : }
1929 2 : return 0;
1930 : }
1931 :
1932 : /*****************************************************************************
1933 : * ParseGridNoMiss() --
1934 : *
1935 : * Arthur Taylor / MDL
1936 : *
1937 : * PURPOSE
1938 : * A helper function for ParseGrid. In this particular case it is dealing
1939 : * with a field that has NO missing value type.
1940 : * Walks through either a float or an integer grid, computing the min/max
1941 : * values in the grid, and converts the units. It uses gridAttrib info for the
1942 : * missing values and it updates gridAttrib with the observed min/max values.
1943 : *
1944 : * ARGUMENTS
1945 : * attrib = Grid Attribute structure already filled in (Input/Output)
1946 : * grib_Data = The place to store the grid data. (Output)
1947 : * Nx, Ny = The dimensions of the grid (Input)
1948 : * iain = Place to find data if it is an Integer (or float). (Input)
1949 : * unitM = M in unit conversion equation y(new) = m x(orig) + b (Input)
1950 : * unitB = B in unit conversion equation y(new) = m x(orig) + b (Input)
1951 : * f_wxType = true if we have a valid wx type. (Input)
1952 : * WxType = table to look up values in. (Input)
1953 : * startX = The start of the X values. (Input)
1954 : * startY = The start of the Y values. (Input)
1955 : * subNx = The Nx dimmension of the subgrid (Input)
1956 : * subNy = The Ny dimmension of the subgrid (Input)
1957 : *
1958 : * FILES/DATABASES: None
1959 : *
1960 : * RETURNS: void
1961 : *
1962 : * HISTORY
1963 : * 12/2002 Arthur Taylor (MDL/RSIS): Created to optimize part of ParseGrid.
1964 : * 5/2003 AAT: Added ability to see if wxType occurs. If so sets table
1965 : * valid to 2, otherwise leaves it at 1. If table valid is 0 then
1966 : * sets value to missing value (if applicable).
1967 : * 2/2004 AAT: Added the subgrid capability.
1968 : *
1969 : * NOTES
1970 : * 1) Don't have to check if value became missing value, because we can check
1971 : * if missing falls in the range of the min/max converted units. If
1972 : * missing does fall in that range we need to move missing.
1973 : * (See f_readjust in ParseGrid)
1974 : *****************************************************************************
1975 : */
1976 0 : static void ParseGridNoMiss (gridAttribType *attrib, double *grib_Data,
1977 : sInt4 Nx, sInt4 Ny, sInt4 *iain,
1978 : double unitM, double unitB, uChar f_wxType,
1979 : sect2_WxType *WxType, int startX, int startY,
1980 : int subNx, int subNy)
1981 : {
1982 : sInt4 x, y; /* Where we are in the grid. */
1983 : double value; /* The data in the new units. */
1984 0 : uChar f_maxmin = 0; /* Flag if max/min is valid yet. */
1985 : uInt4 index; /* Current index into Wx table. */
1986 0 : sInt4 *itemp = NULL;
1987 0 : float *ftemp = NULL;
1988 :
1989 : /* Resolve possibility that the data is an integer or a float and find
1990 : * max/min values. (see note 1) */
1991 0 : for (y = 0; y < subNy; y++) {
1992 0 : if (((startY + y - 1) < 0) || ((startY + y - 1) >= Ny)) {
1993 0 : for (x = 0; x < subNx; x++) {
1994 0 : *grib_Data++ = 9999;
1995 : }
1996 : } else {
1997 0 : if (attrib->fieldType) {
1998 0 : itemp = iain + (startY + y - 1) * Nx + (startX - 1);
1999 : } else {
2000 0 : ftemp = ((float *) iain) + (startY + y - 1) * Nx + (startX - 1);
2001 : }
2002 0 : for (x = 0; x < subNx; x++) {
2003 0 : if (((startX + x - 1) < 0) || ((startX + x - 1) >= Nx)) {
2004 0 : *grib_Data++ = 9999;
2005 : } else {
2006 : /* Convert the units. */
2007 0 : if (attrib->fieldType) {
2008 0 : if (unitM == -10) {
2009 0 : value = pow (10.0, (*itemp++));
2010 : } else {
2011 0 : value = unitM * (*itemp++) + unitB;
2012 : }
2013 : } else {
2014 0 : if (unitM == -10) {
2015 0 : value = pow (10.0, (double) (*ftemp++));
2016 : } else {
2017 0 : value = unitM * (*ftemp++) + unitB;
2018 : }
2019 : }
2020 0 : if (f_wxType) {
2021 0 : index = (uInt4) value;
2022 0 : if (index < WxType->dataLen) {
2023 0 : if (WxType->ugly[index].f_valid == 1) {
2024 0 : WxType->ugly[index].f_valid = 2;
2025 0 : } else if (WxType->ugly[index].f_valid == 0) {
2026 : /* Table is not valid here so set value to missing? */
2027 : /* No missing value, so use index = WxType->dataLen? */
2028 : /* No... set f_valid to 3 so we know we used this
2029 : * invalid element, then handle it in degrib2.c ::
2030 : * ReadGrib2Record() where we set it back to 0. */
2031 0 : WxType->ugly[index].f_valid = 3;
2032 : }
2033 : }
2034 : }
2035 0 : if (f_maxmin) {
2036 0 : if (value < attrib->min) {
2037 0 : attrib->min = value;
2038 0 : } else if (value > attrib->max) {
2039 0 : attrib->max = value;
2040 : }
2041 : } else {
2042 0 : attrib->min = attrib->max = value;
2043 0 : f_maxmin = 1;
2044 : }
2045 0 : *grib_Data++ = value;
2046 : }
2047 : }
2048 : }
2049 : }
2050 0 : attrib->f_maxmin = f_maxmin;
2051 0 : }
2052 :
2053 : /*****************************************************************************
2054 : * ParseGridPrimMiss() --
2055 : *
2056 : * Arthur Taylor / MDL
2057 : *
2058 : * PURPOSE
2059 : * A helper function for ParseGrid. In this particular case it is dealing
2060 : * with a field that has primary missing value type.
2061 : * Walks through either a float or an integer grid, computing the min/max
2062 : * values in the grid, and converts the units. It uses gridAttrib info for the
2063 : * missing values and it updates gridAttrib with the observed min/max values.
2064 : *
2065 : * ARGUMENTS
2066 : * attrib = sect 5 structure already filled in by ParseSect5 (In/Output)
2067 : * grib_Data = The place to store the grid data. (Output)
2068 : * Nx, Ny = The dimensions of the grid (Input)
2069 : * iain = Place to find data if it is an Integer (or float). (Input)
2070 : * unitM = M in unit conversion equation y(new) = m x(orig) + b (Input)
2071 : * unitB = B in unit conversion equation y(new) = m x(orig) + b (Input)
2072 : * f_wxType = true if we have a valid wx type. (Input)
2073 : * WxType = table to look up values in. (Input)
2074 : * startX = The start of the X values. (Input)
2075 : * startY = The start of the Y values. (Input)
2076 : * subNx = The Nx dimmension of the subgrid (Input)
2077 : * subNy = The Ny dimmension of the subgrid (Input)
2078 : *
2079 : * FILES/DATABASES: None
2080 : *
2081 : * RETURNS: void
2082 : *
2083 : * HISTORY
2084 : * 12/2002 Arthur Taylor (MDL/RSIS): Created to optimize part of ParseGrid.
2085 : * 5/2003 AAT: Added ability to see if wxType occurs. If so sets table
2086 : * valid to 2, otherwise leaves it at 1. If table valid is 0 then
2087 : * sets value to missing value (if applicable).
2088 : * 2/2004 AAT: Added the subgrid capability.
2089 : *
2090 : * NOTES
2091 : * 1) Don't have to check if value became missing value, because we can check
2092 : * if missing falls in the range of the min/max converted units. If
2093 : * missing does fall in that range we need to move missing.
2094 : * (See f_readjust in ParseGrid)
2095 : *****************************************************************************
2096 : */
2097 2 : static void ParseGridPrimMiss (gridAttribType *attrib, double *grib_Data,
2098 : sInt4 Nx, sInt4 Ny, sInt4 *iain,
2099 : double unitM, double unitB, sInt4 *missCnt,
2100 : uChar f_wxType, sect2_WxType *WxType,
2101 : int startX, int startY, int subNx, int subNy)
2102 : {
2103 : sInt4 x, y; /* Where we are in the grid. */
2104 : double value; /* The data in the new units. */
2105 2 : uChar f_maxmin = 0; /* Flag if max/min is valid yet. */
2106 : uInt4 index; /* Current index into Wx table. */
2107 2 : sInt4 *itemp = NULL;
2108 2 : float *ftemp = NULL;
2109 : /* float *ain = (float *) iain;*/
2110 :
2111 : /* Resolve possibility that the data is an integer or a float and find
2112 : * max/min values. (see note 1) */
2113 260 : for (y = 0; y < subNy; y++) {
2114 258 : if (((startY + y - 1) < 0) || ((startY + y - 1) >= Ny)) {
2115 0 : for (x = 0; x < subNx; x++) {
2116 0 : *grib_Data++ = attrib->missPri;
2117 0 : (*missCnt)++;
2118 : }
2119 : } else {
2120 258 : if (attrib->fieldType) {
2121 0 : itemp = iain + (startY + y - 1) * Nx + (startX - 1);
2122 : } else {
2123 258 : ftemp = ((float *) iain) + (startY + y - 1) * Nx + (startX - 1);
2124 : }
2125 45924 : for (x = 0; x < subNx; x++) {
2126 45666 : if (((startX + x - 1) < 0) || ((startX + x - 1) >= Nx)) {
2127 0 : *grib_Data++ = attrib->missPri;
2128 0 : (*missCnt)++;
2129 : } else {
2130 45666 : if (attrib->fieldType) {
2131 0 : value = (*itemp++);
2132 : } else {
2133 45666 : value = (*ftemp++);
2134 : }
2135 :
2136 : /* Make sure value is not a missing value when converting
2137 : * units, and while computing max/min. */
2138 45666 : if (value == attrib->missPri) {
2139 7512 : (*missCnt)++;
2140 : } else {
2141 : /* Convert the units. */
2142 38154 : if (unitM == -10) {
2143 0 : value = pow (10.0, value);
2144 : } else {
2145 38154 : value = unitM * value + unitB;
2146 : }
2147 38154 : if (f_wxType) {
2148 0 : index = (uInt4) value;
2149 0 : if (index < WxType->dataLen) {
2150 0 : if (WxType->ugly[index].f_valid) {
2151 0 : WxType->ugly[index].f_valid = 2;
2152 : } else {
2153 : /* Table is not valid here so set value to missPri
2154 : */
2155 0 : value = attrib->missPri;
2156 0 : (*missCnt)++;
2157 : }
2158 : }
2159 : }
2160 38154 : if ((!f_wxType) || (value != attrib->missPri)) {
2161 38154 : if (f_maxmin) {
2162 38152 : if (value < attrib->min) {
2163 30 : attrib->min = value;
2164 38122 : } else if (value > attrib->max) {
2165 0 : attrib->max = value;
2166 : }
2167 : } else {
2168 2 : attrib->min = attrib->max = value;
2169 2 : f_maxmin = 1;
2170 : }
2171 : }
2172 : }
2173 45666 : *grib_Data++ = value;
2174 : }
2175 : }
2176 : }
2177 : }
2178 2 : attrib->f_maxmin = f_maxmin;
2179 2 : }
2180 :
2181 : /*****************************************************************************
2182 : * ParseGridSecMiss() --
2183 : *
2184 : * Arthur Taylor / MDL
2185 : *
2186 : * PURPOSE
2187 : * A helper function for ParseGrid. In this particular case it is dealing
2188 : * with a field that has NO missing value type.
2189 : * Walks through either a float or an integer grid, computing the min/max
2190 : * values in the grid, and converts the units. It uses gridAttrib info for the
2191 : * missing values and it updates gridAttrib with the observed min/max values.
2192 : *
2193 : * ARGUMENTS
2194 : * attrib = sect 5 structure already filled in by ParseSect5 (In/Output)
2195 : * grib_Data = The place to store the grid data. (Output)
2196 : * Nx, Ny = The dimensions of the grid (Input)
2197 : * iain = Place to find data if it is an Integer (or float). (Input)
2198 : * unitM = M in unit conversion equation y(new) = m x(orig) + b (Input)
2199 : * unitB = B in unit conversion equation y(new) = m x(orig) + b (Input)
2200 : * f_wxType = true if we have a valid wx type. (Input)
2201 : * WxType = table to look up values in. (Input)
2202 : * startX = The start of the X values. (Input)
2203 : * startY = The start of the Y values. (Input)
2204 : * subNx = The Nx dimmension of the subgrid (Input)
2205 : * subNy = The Ny dimmension of the subgrid (Input)
2206 : *
2207 : * FILES/DATABASES: None
2208 : *
2209 : * RETURNS: void
2210 : *
2211 : * HISTORY
2212 : * 12/2002 Arthur Taylor (MDL/RSIS): Created to optimize part of ParseGrid.
2213 : * 5/2003 AAT: Added ability to see if wxType occurs. If so sets table
2214 : * valid to 2, otherwise leaves it at 1. If table valid is 0 then
2215 : * sets value to missing value (if applicable).
2216 : * 2/2004 AAT: Added the subgrid capability.
2217 : *
2218 : * NOTES
2219 : * 1) Don't have to check if value became missing value, because we can check
2220 : * if missing falls in the range of the min/max converted units. If
2221 : * missing does fall in that range we need to move missing.
2222 : * (See f_readjust in ParseGrid)
2223 : *****************************************************************************
2224 : */
2225 0 : static void ParseGridSecMiss (gridAttribType *attrib, double *grib_Data,
2226 : sInt4 Nx, sInt4 Ny, sInt4 *iain,
2227 : double unitM, double unitB, sInt4 *missCnt,
2228 : uChar f_wxType, sect2_WxType *WxType,
2229 : int startX, int startY, int subNx, int subNy)
2230 : {
2231 : sInt4 x, y; /* Where we are in the grid. */
2232 : double value; /* The data in the new units. */
2233 0 : uChar f_maxmin = 0; /* Flag if max/min is valid yet. */
2234 : uInt4 index; /* Current index into Wx table. */
2235 0 : sInt4 *itemp = NULL;
2236 0 : float *ftemp = NULL;
2237 : /* float *ain = (float *) iain;*/
2238 :
2239 : /* Resolve possibility that the data is an integer or a float and find
2240 : * max/min values. (see note 1) */
2241 0 : for (y = 0; y < subNy; y++) {
2242 0 : if (((startY + y - 1) < 0) || ((startY + y - 1) >= Ny)) {
2243 0 : for (x = 0; x < subNx; x++) {
2244 0 : *grib_Data++ = attrib->missPri;
2245 0 : (*missCnt)++;
2246 : }
2247 : } else {
2248 0 : if (attrib->fieldType) {
2249 0 : itemp = iain + (startY + y - 1) * Nx + (startX - 1);
2250 : } else {
2251 0 : ftemp = ((float *) iain) + (startY + y - 1) * Nx + (startX - 1);
2252 : }
2253 0 : for (x = 0; x < subNx; x++) {
2254 0 : if (((startX + x - 1) < 0) || ((startX + x - 1) >= Nx)) {
2255 0 : *grib_Data++ = attrib->missPri;
2256 0 : (*missCnt)++;
2257 : } else {
2258 0 : if (attrib->fieldType) {
2259 0 : value = (*itemp++);
2260 : } else {
2261 0 : value = (*ftemp++);
2262 : }
2263 :
2264 : /* Make sure value is not a missing value when converting
2265 : * units, and while computing max/min. */
2266 0 : if ((value == attrib->missPri) || (value == attrib->missSec)) {
2267 0 : (*missCnt)++;
2268 : } else {
2269 : /* Convert the units. */
2270 0 : if (unitM == -10) {
2271 0 : value = pow (10.0, value);
2272 : } else {
2273 0 : value = unitM * value + unitB;
2274 : }
2275 0 : if (f_wxType) {
2276 0 : index = (uInt4) value;
2277 0 : if (index < WxType->dataLen) {
2278 0 : if (WxType->ugly[index].f_valid) {
2279 0 : WxType->ugly[index].f_valid = 2;
2280 : } else {
2281 : /* Table is not valid here so set value to missPri
2282 : */
2283 0 : value = attrib->missPri;
2284 0 : (*missCnt)++;
2285 : }
2286 : }
2287 : }
2288 0 : if ((!f_wxType) || (value != attrib->missPri)) {
2289 0 : if (f_maxmin) {
2290 0 : if (value < attrib->min) {
2291 0 : attrib->min = value;
2292 0 : } else if (value > attrib->max) {
2293 0 : attrib->max = value;
2294 : }
2295 : } else {
2296 0 : attrib->min = attrib->max = value;
2297 0 : f_maxmin = 1;
2298 : }
2299 : }
2300 : }
2301 0 : *grib_Data++ = value;
2302 : }
2303 : }
2304 : }
2305 : }
2306 0 : attrib->f_maxmin = f_maxmin;
2307 0 : }
2308 :
2309 : /*****************************************************************************
2310 : * ParseGrid() --
2311 : *
2312 : * Arthur Taylor / MDL
2313 : *
2314 : * PURPOSE
2315 : * To walk through the 2 possible grids (and possible bitmap) created by
2316 : * UNPK_GRIB2, and combine the info into 1 grid, at the same time computing
2317 : * the min/max values in the grid. It uses gridAttrib info for the missing values
2318 : * and it then updates the gridAttrib structure for the min/max values that it
2319 : * found.
2320 : * It also uses scan, and ScanIndex2XY, to parse the data and organize the
2321 : * Grib_Data so that 0,0 is the lower left part of the grid, it then traverses
2322 : * the row and then moved up to the next row starting on the left.
2323 : *
2324 : * ARGUMENTS
2325 : * attrib = sect 5 structure already filled in by ParseSect5 (In/Output)
2326 : * Grib_Data = The place to store the grid data. (Output)
2327 : * grib_DataLen = The current size of Grib_Data (can increase) (Input/Output)
2328 : * Nx, Ny = The dimensions of the grid (Input)
2329 : * scan = How to walk through the original grid. (Input)
2330 : * iain = Place to find data if it is an Integer (or float). (Input)
2331 : * ibitmap = Flag stating the data has a bitmap for missing values (In)
2332 : * ib = Where to find the bitmap if we have one (Input)
2333 : * unitM = M in unit conversion equation y(new) = m x(orig) + b (Input)
2334 : * unitB = B in unit conversion equation y(new) = m x(orig) + b (Input)
2335 : * f_wxType = true if we have a valid wx type. (Input)
2336 : * WxType = table to look up values in. (Input)
2337 : * f_subGrid = True if we have a subgrid, false if not. (Input)
2338 : * startX stopX = The bounds of the subgrid in X. (0,-1) means full grid (In)
2339 : * startY stopY = The bounds of the subgrid in Y. (0,-1) means full grid (In)
2340 : *
2341 : * FILES/DATABASES: None
2342 : *
2343 : * RETURNS: void
2344 : *
2345 : * HISTORY
2346 : * 9/2002 Arthur Taylor (MDL/RSIS): Created.
2347 : * 11/2002 AAT: Added unit conversion to metaparse.c
2348 : * 12/2002 AAT: Optimized first loop to make it assume scan 0100 (64)
2349 : * (valid 99.9%), but still have slow loop for generic case.
2350 : * 5/2003 AAT: Added ability to see if wxType occurs. If so sets table
2351 : * valid to 2, otherwise leaves it at 1. If table valid is 0 then
2352 : * sets value to missing value (if applicable).
2353 : * 7/2003 AAT: added check if f_maxmin before checking if missing was in
2354 : * range of max, min for "readjust" check.
2355 : * 2/2004 AAT: Added startX / startY / stopX / stopY
2356 : * 5/2004 AAT: Found out that I used the opposite definition for bitmap
2357 : * 0 = missing, 1 = valid.
2358 : *
2359 : * NOTES
2360 : *****************************************************************************
2361 : */
2362 2 : void ParseGrid (gridAttribType *attrib, double **Grib_Data,
2363 : uInt4 *grib_DataLen, uInt4 Nx, uInt4 Ny, int scan,
2364 : sInt4 *iain, sInt4 ibitmap, sInt4 *ib, double unitM,
2365 : double unitB, uChar f_wxType, sect2_WxType *WxType,
2366 : uChar f_subGrid, int startX, int startY, int stopX, int stopY)
2367 : {
2368 : double xmissp; /* computed missing value needed for ibitmap = 1,
2369 : * Also used if unit conversion causes confusion
2370 : * over_ missing values. */
2371 : double xmisss; /* Used if unit conversion causes confusion over
2372 : * missing values. */
2373 : uChar f_readjust; /* True if unit conversion caused confusion over
2374 : * missing values. */
2375 : uInt4 scanIndex; /* Where we are in the original grid. */
2376 : sInt4 x, y; /* Where we are in a grid of scan value 0100 */
2377 : sInt4 newIndex; /* x,y in a 1 dimensional array. */
2378 : double value; /* The data in the new units. */
2379 : double *grib_Data; /* A pointer to Grib_Data for ease of manipulation. */
2380 2 : sInt4 missCnt = 0; /* Number of detected missing values. */
2381 : uInt4 index; /* Current index into Wx table. */
2382 2 : float *ain = (float *) iain;
2383 : uInt4 subNx; /* The Nx dimmension of the subgrid. */
2384 : uInt4 subNy; /* The Ny dimmension of the subgrid. */
2385 :
2386 2 : subNx = stopX - startX + 1;
2387 2 : subNy = stopY - startY + 1;
2388 :
2389 : myAssert (((!f_subGrid) && (subNx == Nx)) || (f_subGrid));
2390 : myAssert (((!f_subGrid) && (subNy == Ny)) || (f_subGrid));
2391 :
2392 2 : if (subNx * subNy > *grib_DataLen) {
2393 2 : *grib_DataLen = subNx * subNy;
2394 : *Grib_Data = (double *) realloc ((void *) (*Grib_Data),
2395 2 : (*grib_DataLen) * sizeof (double));
2396 : }
2397 2 : grib_Data = *Grib_Data;
2398 :
2399 : /* Resolve possibility that the data is an integer or a float, find
2400 : * max/min values, and do unit conversion. (see note 1) */
2401 2 : if (scan == 64) {
2402 2 : if (attrib->f_miss == 0) {
2403 : ParseGridNoMiss (attrib, grib_Data, Nx, Ny, iain, unitM, unitB,
2404 0 : f_wxType, WxType, startX, startY, subNx, subNy);
2405 2 : } else if (attrib->f_miss == 1) {
2406 : ParseGridPrimMiss (attrib, grib_Data, Nx, Ny, iain, unitM, unitB,
2407 : &missCnt, f_wxType, WxType, startX, startY,
2408 2 : subNx, subNy);
2409 0 : } else if (attrib->f_miss == 2) {
2410 : ParseGridSecMiss (attrib, grib_Data, Nx, Ny, iain, unitM, unitB,
2411 : &missCnt, f_wxType, WxType, startX, startY, subNx,
2412 0 : subNy);
2413 : }
2414 : } else {
2415 : /* Internally we use scan = 0100. Scan is usually 0100 from the
2416 : * unpacker library, but if scan is not, the following code converts
2417 : * it. We optimized the previous (scan 0100) case by calling a
2418 : * dedicated procedure. Here we don't since for scan != 0100, we
2419 : * would_ need a different unpacker library, which is extremely
2420 : * unlikely. */
2421 0 : for (scanIndex = 0; scanIndex < Nx * Ny; scanIndex++) {
2422 0 : if (attrib->fieldType) {
2423 0 : value = iain[scanIndex];
2424 : } else {
2425 0 : value = ain[scanIndex];
2426 : }
2427 : /* Make sure value is not a missing value when converting units, and
2428 : * while computing max/min. */
2429 0 : if ((attrib->f_miss == 0) ||
2430 : ((attrib->f_miss == 1) && (value != attrib->missPri)) ||
2431 : ((attrib->f_miss == 2) && (value != attrib->missPri) &&
2432 : (value != attrib->missSec))) {
2433 : /* Convert the units. */
2434 0 : if (unitM == -10) {
2435 0 : value = pow (10.0, value);
2436 : } else {
2437 0 : value = unitM * value + unitB;
2438 : }
2439 : /* Don't have to check if value became missing value, because we
2440 : * can check if missing falls in the range of min/max. If
2441 : * missing does fall in that range we need to move missing. See
2442 : * f_readjust */
2443 0 : if (f_wxType) {
2444 0 : index = (uInt4) value;
2445 0 : if (index < WxType->dataLen) {
2446 0 : if (WxType->ugly[index].f_valid == 1) {
2447 0 : WxType->ugly[index].f_valid = 2;
2448 0 : } else if (WxType->ugly[index].f_valid == 0) {
2449 : /* Table is not valid here so set value to missPri */
2450 0 : if (attrib->f_miss != 0) {
2451 0 : value = attrib->missPri;
2452 0 : missCnt++;
2453 : } else {
2454 : /* No missing value, so use index = WxType->dataLen */
2455 : /* No... set f_valid to 3 so we know we used this
2456 : * invalid element, then handle it in degrib2.c ::
2457 : * ReadGrib2Record() where we set it back to 0. */
2458 0 : WxType->ugly[index].f_valid = 3;
2459 : }
2460 : }
2461 : }
2462 : }
2463 0 : if ((!f_wxType) ||
2464 : ((attrib->f_miss == 0) || (value != attrib->missPri))) {
2465 0 : if (attrib->f_maxmin) {
2466 0 : if (value < attrib->min) {
2467 0 : attrib->min = value;
2468 0 : } else if (value > attrib->max) {
2469 0 : attrib->max = value;
2470 : }
2471 : } else {
2472 0 : attrib->min = attrib->max = value;
2473 0 : attrib->f_maxmin = 1;
2474 : }
2475 : }
2476 : } else {
2477 0 : missCnt++;
2478 : }
2479 0 : ScanIndex2XY (scanIndex, &x, &y, scan, Nx, Ny);
2480 : /* ScanIndex returns value as if scan was 0100 */
2481 0 : newIndex = (x - 1) + (y - 1) * Nx;
2482 0 : grib_Data[newIndex] = value;
2483 : }
2484 : }
2485 :
2486 : /* Deal with possibility that unit conversion ended up with valid numbers
2487 : * being interpreted as missing. */
2488 2 : f_readjust = 0;
2489 2 : xmissp = attrib->missPri;
2490 2 : xmisss = attrib->missSec;
2491 2 : if (attrib->f_maxmin) {
2492 2 : if ((attrib->f_miss == 1) || (attrib->f_miss == 2)) {
2493 2 : if ((attrib->missPri >= attrib->min) &&
2494 : (attrib->missPri <= attrib->max)) {
2495 0 : xmissp = attrib->max + 1;
2496 0 : f_readjust = 1;
2497 : }
2498 2 : if (attrib->f_miss == 2) {
2499 0 : if ((attrib->missSec >= attrib->min) &&
2500 : (attrib->missSec <= attrib->max)) {
2501 0 : xmisss = attrib->max + 2;
2502 0 : f_readjust = 1;
2503 : }
2504 : }
2505 : }
2506 : }
2507 :
2508 : /* Walk through the grid, resetting the missing values, as determined by
2509 : * the original grid. */
2510 2 : if (f_readjust) {
2511 0 : for (scanIndex = 0; scanIndex < Nx * Ny; scanIndex++) {
2512 0 : ScanIndex2XY (scanIndex, &x, &y, scan, Nx, Ny);
2513 : /* ScanIndex returns value as if scan was 0100 */
2514 0 : newIndex = (x - 1) + (y - 1) * Nx;
2515 0 : if (attrib->fieldType) {
2516 0 : value = iain[scanIndex];
2517 : } else {
2518 0 : value = ain[scanIndex];
2519 : }
2520 0 : if (value == attrib->missPri) {
2521 0 : grib_Data[newIndex] = xmissp;
2522 0 : } else if ((attrib->f_miss == 2) && (value == attrib->missSec)) {
2523 0 : grib_Data[newIndex] = xmisss;
2524 : }
2525 : }
2526 0 : attrib->missPri = xmissp;
2527 0 : if (attrib->f_miss == 2) {
2528 0 : attrib->missSec = xmisss;
2529 : }
2530 : }
2531 :
2532 : /* Resolve bitmap (if there is one) in the data. */
2533 2 : if (ibitmap) {
2534 0 : attrib->f_maxmin = 0;
2535 0 : if ((attrib->f_miss != 1) && (attrib->f_miss != 2)) {
2536 0 : missCnt = 0;
2537 : /* Figure out a missing value. */
2538 0 : xmissp = 9999;
2539 0 : if (attrib->f_maxmin) {
2540 0 : if ((xmissp <= attrib->max) && (xmissp >= attrib->min)) {
2541 0 : xmissp = attrib->max + 1;
2542 : }
2543 : }
2544 : /* embed the missing value. */
2545 0 : for (scanIndex = 0; scanIndex < Nx * Ny; scanIndex++) {
2546 0 : ScanIndex2XY (scanIndex, &x, &y, scan, Nx, Ny);
2547 : /* ScanIndex returns value as if scan was 0100 */
2548 0 : newIndex = (x - 1) + (y - 1) * Nx;
2549 : /* Corrected this on 5/10/2004 */
2550 0 : if (ib[scanIndex] != 1) {
2551 0 : grib_Data[newIndex] = xmissp;
2552 0 : missCnt++;
2553 : } else {
2554 0 : if (!attrib->f_maxmin) {
2555 0 : attrib->f_maxmin = 1;
2556 0 : attrib->max = attrib->min = grib_Data[newIndex];
2557 : } else {
2558 0 : if (attrib->max < grib_Data[newIndex])
2559 0 : attrib->max = grib_Data[newIndex];
2560 0 : if (attrib->min > grib_Data[newIndex])
2561 0 : attrib->min = grib_Data[newIndex];
2562 : }
2563 : }
2564 : }
2565 0 : attrib->f_miss = 1;
2566 0 : attrib->missPri = xmissp;
2567 : }
2568 0 : if (!attrib->f_maxmin) {
2569 0 : attrib->f_maxmin = 1;
2570 0 : attrib->max = attrib->min = xmissp;
2571 : }
2572 : }
2573 2 : attrib->numMiss = missCnt;
2574 2 : }
2575 :
2576 : typedef struct {
2577 : double value;
2578 : int cnt;
2579 : } freqType;
2580 :
2581 0 : int freqCompare (const void *A, const void *B)
2582 : {
2583 0 : const freqType *a = (freqType *) A;
2584 0 : const freqType *b = (freqType *) B;
2585 :
2586 0 : if (a->value < b->value)
2587 0 : return -1;
2588 0 : if (a->value > b->value)
2589 0 : return 1;
2590 0 : return 0;
2591 : }
2592 :
2593 0 : void FreqPrint (char **ans, double *Data, sInt4 DataLen, sInt4 Nx,
2594 : sInt4 Ny, sChar decimal, char *comment)
2595 : {
2596 : int x, y, i;
2597 : double *ptr;
2598 : double value;
2599 0 : freqType *freq = NULL;
2600 0 : int numFreq = 0;
2601 : char format[20];
2602 :
2603 : myAssert (*ans == NULL);
2604 :
2605 0 : if ((Nx < 0) || (Ny < 0) || (Nx * Ny > DataLen)) {
2606 0 : return;
2607 : }
2608 :
2609 0 : ptr = Data;
2610 0 : for (y = 0; y < Ny; y++) {
2611 0 : for (x = 0; x < Nx; x++) {
2612 : /* 2/28/2006 Introduced value to round before putting the data in
2613 : * the Freq table. */
2614 0 : value = myRound (*ptr, decimal);
2615 0 : for (i = 0; i < numFreq; i++) {
2616 0 : if (value == freq[i].value) {
2617 0 : freq[i].cnt++;
2618 0 : break;
2619 : }
2620 : }
2621 0 : if (i == numFreq) {
2622 0 : numFreq++;
2623 0 : freq = (freqType *) realloc (freq, numFreq * sizeof (freqType));
2624 0 : freq[i].value = value;
2625 0 : freq[i].cnt = 1;
2626 : }
2627 0 : ptr++;
2628 : }
2629 : }
2630 :
2631 0 : qsort (freq, numFreq, sizeof (freq[0]), freqCompare);
2632 :
2633 0 : mallocSprintf (ans, "%s | count\n", comment);
2634 0 : sprintf (format, "%%.%df | %%d\n", decimal);
2635 0 : for (i = 0; i < numFreq; i++) {
2636 0 : reallocSprintf (ans, format, myRound (freq[i].value, decimal),
2637 0 : freq[i].cnt);
2638 : }
2639 0 : free (freq);
2640 : }
|
