1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <math.h>
5 #include "tdlpack.h"
6 #include "myerror.h"
7 #include "meta.h"
8 #include "memendian.h"
9 #include "fileendian.h"
10 #include "myassert.h"
11 #include "myutil.h"
12 #include "clock.h"
13 #ifdef MEMWATCH
14 #include "memwatch.h"
15 #endif
16
17 #define GRIB_UNSIGN_INT3(a,b,c) ((a<<16)+(b<<8)+c)
18 #define GRIB_UNSIGN_INT2(a,b) ((a<<8)+b)
19 #define GRIB_SIGN_INT3(a,b,c) ((1-(int) ((unsigned) (a & 0x80) >> 6)) * (int) (((a & 127) << 16)+(b<<8)+c))
20 #define GRIB_SIGN_INT2(a,b) ((1-(int) ((unsigned) (a & 0x80) >> 6)) * (int) (((a & 0x7f) << 8) + b))
21
22 /* *INDENT-OFF* */
23 TDLP_TableType TDLP_B_Table[5] = {
24 /* 0 */ {0, "Continuous field"},
25 /* 1 */ {1, "Point Binary - cumulative from above"},
26 /* 2 */ {2, "Point Binary - cumulative from below"},
27 /* 3 */ {3, "Point Binary - discrete"},
28 /* 4 */ {5, "Grid Binary - "},
29 };
30
31 TDLP_TableType TDLP_DD_Table[9] = {
32 /* 0 */ {0, "Independent of model"},
33 /* 1 */ {6, "NGM model"},
34 /* 2 */ {7, "Eta model"},
35 /* 3 */ {8, "AVN model"},
36 /* 4 */ {9, "MRF model"},
37 /* 5 */ {79, "NDFD forecast"},
38 /* 6 */ {80, "MOS AEV forecast"},
39 /* 7 */ {81, "Local AEV Firecasts"},
40 /* 8 */ {82, "Obs matching AEV Forecasts"},
41 };
42
43 TDLP_TableType TDLP_V_Table[4] = {
44 /* 0 */ {0, "No vertical processing"},
45 /* 1 */ {1, "Difference Levels (UUUU - LLLL)"},
46 /* 2 */ {2, "Sum Levels (UUUU + LLLL)"},
47 /* 3 */ {3, "Mean Levels (UUUU + LLLL) / 2."},
48 };
49
50 TDLP_TableType TDLP_T_Table[3] = {
51 /* 0 */ {0, "No nolinear tranform"},
52 /* 1 */ {1, "Square transform"},
53 /* 2 */ {2, "Square root tranform"},
54 };
55
56 TDLP_TableType TDLP_Oper_Table[9] = {
57 /* 0 */ {0, "No time operator"},
58 /* 1 */ {1, "Mean (Var 1, Var 2)"},
59 /* 2 */ {2, "Difference (Var 1 - Var 2)"},
60 /* 3 */ {3, "Maximum (Var 1, Var 2)"},
61 /* 4 */ {4, "Minimum (Var 1, Var 2)"},
62 /* 5 */ {5, "Mean (Var 1, Var 2)"},
63 /* 6 */ {6, "Difference (Var 2 - Var 1)"},
64 /* 7 */ {7, "Maximum (Var 1, Var 2)"},
65 /* 8 */ {8, "Minimum (Var 1, Var 2)"},
66 };
67
68 TDLP_TableType TDLP_I_Table[4] = {
69 /* 0 */ {0, "No interpolation"},
70 /* 1 */ {1, "Bi-quadratic interpolation"},
71 /* 2 */ {2, "Bi-linear interpolation"},
72 /* 3 */ {3, "Special interpolation for QPF"},
73 };
74
75 TDLP_TableType TDLP_S_Table[6] = {
76 /* 0 */ {0, "No smoothing"},
77 /* 1 */ {1, "5-point smoothing"},
78 /* 2 */ {2, "9-point smoothing"},
79 /* 3 */ {3, "25-point smoothing"},
80 /* 4 */ {4, "81-point smoothing"},
81 /* 5 */ {5, "168-point smoothing"},
82 };
83 /* *INDENT-ON* */
84
85 typedef struct {
86 sInt4 min; /* Minimum value in a group. Can cause problems if
87 * this is unsigned. */
88 uChar bit; /* # of bits in a group. 2^31 is the largest # of
89 * bits that can be used to hold # of bits in a
90 * group. However, the # of bits for a number can't
91 * be > 64, and is probably < 32, so bit is < 6 and
92 * probably < 5. */
93 uInt4 num; /* number of values in the group. May need this to be
94 * signed. */
95 sInt4 max; /* Max value for the group */
96 uInt4 start; /* index in Data where group starts. */
97 uChar f_trySplit; /* Flag, if we have tried to split this group. */
98 uChar f_tryShift; /* Flag, if we have tried to shift this group. */
99 } TDLGroupType;
100
101 /*****************************************************************************
102 * ReadTDLPSect1() --
103 *
104 * Arthur Taylor / MDL
105 *
106 * PURPOSE
107 * Parses the TDLP "Product Definition Section" or section 1, filling out
108 * the pdsMeta data structure.
109 *
110 * ARGUMENTS
111 * pds = The compressed part of the message dealing with "PDS". (Input)
112 * gribLen = The total length of the TDLP message. (Input)
113 * curLoc = Current location in the TDLP message. (Output)
114 * pdsMeta = The filled out pdsMeta data structure. (Output)
115 * f_gds = boolean if there is a Grid Definition Section. (Output)
116 * f_bms = boolean if there is a Bitmap Section. (Output)
117 * DSF = Decimal Scale Factor for unpacking the data. (Output)
118 * BSF = Binary Scale Factor for unpacking the data. (Output)
119 *
120 * FILES/DATABASES: None
121 *
122 * RETURNS: int (could use errSprintf())
123 * 0 = OK
124 * -1 = gribLen is too small.
125 *
126 * HISTORY
127 * 10/2004 Arthur Taylor (MDL): Created.
128 *
129 * NOTES
130 *****************************************************************************
131 */
ReadTDLPSect1(uChar * pds,sInt4 tdlpLen,sInt4 * curLoc,pdsTDLPType * pdsMeta,char * f_gds,char * f_bms,short int * DSF,short int * BSF)132 static int ReadTDLPSect1 (uChar *pds, sInt4 tdlpLen, sInt4 *curLoc,
133 pdsTDLPType * pdsMeta, char *f_gds, char *f_bms,
134 short int *DSF, short int *BSF)
135 {
136 char sectLen; /* Length of section. */
137 sInt4 li_temp; /* Temporary variable. */
138 int W, XXXX, YY; /* Helps with Threshold calculation. */
139 int year, t_year; /* The reference year, and a consistency test */
140 uChar month, t_month; /* The reference month, and a consistency test */
141 uChar day, t_day; /* The reference day, and a consistency test */
142 uChar hour, t_hour; /* The reference hour, and a consistency test */
143 uChar min; /* The reference minute */
144 uShort2 project_hr; /* The projection in hours. */
145 sInt4 tau; /* Used to cross check project_hr */
146 int lenPL; /* Length of the Plain Language descriptor. */
147
148 sectLen = *(pds++);
149 *curLoc += sectLen;
150 if (*curLoc > tdlpLen) {
151 errSprintf ("Ran out of data in PDS (TDLP Section 1)\n");
152 return -1;
153 }
154 myAssert (sectLen <= 71);
155 if (sectLen < 39) {
156 errSprintf ("TDLP Section 1 is too small.\n");
157 return -1;
158 }
159 *f_bms = (GRIB2BIT_7 & *pds) ? 1 : 0;
160 *f_gds = (GRIB2BIT_8 & *pds) ? 1 : 0;
161 pds++;
162 year = GRIB_UNSIGN_INT2 (*pds, pds[1]);
163 pds += 2;
164 month = *(pds++);
165 day = *(pds++);
166 hour = *(pds++);
167 min = *(pds++);
168 MEMCPY_BIG (&li_temp, pds, sizeof (sInt4));
169 pds += 4;
170 t_year = li_temp / 1000000L;
171 li_temp -= t_year * 1000000L;
172 t_month = li_temp / 10000L;
173 li_temp -= t_month * 10000L;
174 t_day = li_temp / 100;
175 t_hour = li_temp - t_day * 100;
176 if ((t_year != year) || (t_month != month) || (t_day != day) ||
177 (t_hour != hour)) {
178 errSprintf ("Error Inconsistent Times in ReadTDLPSect1.\n");
179 return -1;
180 }
181 if (ParseTime (&(pdsMeta->refTime), year, month, day, hour, min, 0) != 0) {
182 preErrSprintf ("Error In call to ParseTime in ReadTDLPSect1.\n");
183 return -1;
184 }
185 MEMCPY_BIG (&(li_temp), pds, sizeof (sInt4));
186 pds += 4;
187 pdsMeta->ID1 = li_temp;
188 pdsMeta->CCC = li_temp / 1000000L;
189 li_temp -= pdsMeta->CCC * 1000000L;
190 pdsMeta->FFF = li_temp / 1000;
191 li_temp -= pdsMeta->FFF * 1000;
192 pdsMeta->B = li_temp / 100;
193 pdsMeta->DD = li_temp - pdsMeta->B * 100;
194 MEMCPY_BIG (&(li_temp), pds, sizeof (sInt4));
195 pds += 4;
196 pdsMeta->ID2 = li_temp;
197 pdsMeta->V = li_temp / 100000000L;
198 li_temp -= pdsMeta->V * 100000000L;
199 pdsMeta->LLLL = li_temp / 10000;
200 pdsMeta->UUUU = li_temp - pdsMeta->LLLL * 10000;
201 MEMCPY_BIG (&(li_temp), pds, sizeof (sInt4));
202 pds += 4;
203 pdsMeta->ID3 = li_temp;
204 pdsMeta->T = li_temp / 100000000L;
205 li_temp -= pdsMeta->T * 100000000L;
206 pdsMeta->RR = li_temp / 1000000L;
207 li_temp -= pdsMeta->RR * 1000000L;
208 pdsMeta->Oper = li_temp / 100000L;
209 li_temp -= pdsMeta->Oper * 100000L;
210 pdsMeta->HH = li_temp / 1000;
211 pdsMeta->ttt = li_temp - pdsMeta->HH * 1000;
212 MEMCPY_BIG (&(li_temp), pds, sizeof (sInt4));
213 pds += 4;
214 pdsMeta->ID4 = li_temp;
215 W = li_temp / 1000000000L;
216 li_temp -= W * 1000000000L;
217 XXXX = li_temp / 100000L;
218 li_temp -= XXXX * 100000L;
219 if (W) {
220 XXXX = -1 * XXXX;
221 }
222 YY = li_temp / 1000L;
223 li_temp -= YY * 1000L;
224 if (YY >= 50) {
225 YY = -1 * (YY - 50);
226 }
227 pdsMeta->thresh = (XXXX / 10000.) * pow (10.0, YY);
228 pdsMeta->I = li_temp / 100;
229 li_temp -= pdsMeta->I * 100L;
230 pdsMeta->S = li_temp / 10;
231 pdsMeta->G = li_temp - pdsMeta->S * 10;
232 project_hr = GRIB_UNSIGN_INT2 (*pds, pds[1]);
233 tau = pdsMeta->ID3 - ((pdsMeta->ID3 / 1000) * 1000);
234 if (tau != project_hr) {
235 printf ("Warning: Inconsistent Projections in hours in "
236 "ReadTDLPSect1 (%d vs %d)\n", tau, project_hr);
237 /*
238 errSprintf ("Warning: Inconsistent Projections in hours in "
239 "ReadTDLPSect1 (%ld vs %d)\n", tau, project_hr);
240 */
241 project_hr = tau;
242 /* return -1; */
243 }
244 pds += 2;
245 pdsMeta->project = project_hr * 3600 + (*(pds++)) * 60;
246 pdsMeta->procNum = (*(pds++));
247 pdsMeta->seqNum = (*(pds++));
248 *DSF = (*pds > 128) ? 128 - (*(pds++)) : (*(pds++));
249 *BSF = (*pds > 128) ? 128 - (*(pds++)) : (*(pds++));
250 if ((*pds != 0) || (pds[1] != 0) || (pds[2] != 0)) {
251 errSprintf ("Error Reserved was not set to 0 in ReadTDLPSect1.\n");
252 return -1;
253 }
254 pds += 3;
255 lenPL = (*(pds++));
256 if (sectLen - lenPL != 39) {
257 errSprintf ("Error sectLen(%d) - lenPL(%d) != 39 in ReadTDLPSect1.\n",
258 sectLen, lenPL);
259 return -1;
260 }
261 if (lenPL > 32) {
262 lenPL = 32;
263 }
264 strncpy (pdsMeta->Descriptor, (char *) pds, lenPL);
265 pdsMeta->Descriptor[lenPL] = '\0';
266 strTrim (pdsMeta->Descriptor);
267 return 0;
268 }
269
270 /*****************************************************************************
271 * TDLP_TableLookUp() --
272 *
273 * Arthur Taylor / MDL
274 *
275 * PURPOSE
276 * To look up TDL Ids information in a given table.
277 *
278 * ARGUMENTS
279 * table = The Table to look up the Id in. (Input)
280 * tableLen = The length of the Table. (Input)
281 * index = The index into the Table. (Input)
282 *
283 * FILES/DATABASES: None
284 *
285 * RETURNS: char *
286 * Returns the meta data that is associated with index in the table.
287 *
288 * HISTORY
289 * 10/2004 Arthur Taylor (MDL): Created.
290 *
291 * NOTES
292 *****************************************************************************
293 */
TDLP_TableLookUp(TDLP_TableType * table,int tableLen,int index)294 static const char *TDLP_TableLookUp(TDLP_TableType * table, int tableLen,
295 int index)
296 {
297 int i; /* Loop counter. */
298
299 for (i = 0; i < tableLen; i++) {
300 if (table[i].index == index)
301 return (table[i].data);
302 }
303 return "Unknown";
304 }
305
306 /*****************************************************************************
307 * PrintPDS_TDLP() --
308 *
309 * Arthur Taylor / MDL
310 *
311 * PURPOSE
312 * To generate the message for the Product Definition Sections of the TDLP
313 * Message.
314 *
315 * ARGUMENTS
316 * pds = The TDLP Product Definition Section to print. (Input)
317 *
318 * FILES/DATABASES: None
319 *
320 * RETURNS: void
321 *
322 * HISTORY
323 * 10/2004 Arthur Taylor (MDL): Created.
324 *
325 * NOTES
326 *****************************************************************************
327 */
PrintPDS_TDLP(pdsTDLPType * pds)328 void PrintPDS_TDLP (pdsTDLPType * pds)
329 {
330 char buffer[25]; /* Stores format of pds1->refTime. */
331
332 /*
333 strftime (buffer, 25, "%m/%d/%Y %H:%M:%S UTC", gmtime (&(pds->refTime)));
334 */
335 Clock_Print (buffer, 25, pds->refTime, "%m/%d/%Y %H:%M:%S UTC", 0);
336
337 Print ("PDS-TDLP", "Reference Time", Prt_S, buffer);
338 Print ("PDS-TDLP", "Plain Language", Prt_S, pds->Descriptor);
339 sprintf (buffer, "%09d", pds->ID1);
340 Print ("PDS-TDLP", "ID 1", Prt_S, buffer);
341 sprintf (buffer, "%09d", pds->ID2);
342 Print ("PDS-TDLP", "ID 2", Prt_S, buffer);
343 sprintf (buffer, "%09d", pds->ID3);
344 Print ("PDS-TDLP", "ID 3", Prt_S, buffer);
345 Print ("PDS-TDLP", "ID 4", Prt_D, pds->ID4);
346 Print ("PDS-TDLP", "Model or Process Number", Prt_D, pds->procNum);
347 Print ("PDS-TDLP", "Sequence Number", Prt_D, pds->seqNum);
348
349 sprintf (buffer, "%03d", pds->CCC);
350 Print ("PDS-TDLP", "ID1-CCC", Prt_S, buffer);
351 sprintf (buffer, "%03d", pds->FFF);
352 Print ("PDS-TDLP", "ID1-FFF", Prt_S, buffer);
353 Print ("PDS-TDLP", "ID1-B", Prt_DS, pds->B,
354 TDLP_TableLookUp (TDLP_B_Table, sizeof (TDLP_B_Table), pds->B));
355 sprintf (buffer, "%02d", pds->DD);
356 Print ("PDS-TDLP", "ID1-DD", Prt_SS, buffer,
357 TDLP_TableLookUp (TDLP_DD_Table, sizeof (TDLP_DD_Table), pds->DD));
358
359 Print ("PDS-TDLP", "ID2-V", Prt_DS, pds->V,
360 TDLP_TableLookUp (TDLP_V_Table, sizeof (TDLP_V_Table), pds->V));
361 sprintf (buffer, "%04d", pds->LLLL);
362 Print ("PDS-TDLP", "ID2-LLLL", Prt_S, buffer);
363 sprintf (buffer, "%04d", pds->UUUU);
364 Print ("PDS-TDLP", "ID2-UUUU", Prt_S, buffer);
365
366 if (pds->Oper != 0) {
367 Print ("PDS-TDLP", "ID3-T", Prt_DS, pds->T,
368 TDLP_TableLookUp (TDLP_T_Table, sizeof (TDLP_T_Table), pds->T));
369 sprintf (buffer, "%02d", pds->RR);
370 Print ("PDS-TDLP", "ID3-RR", Prt_SS, buffer,
371 "Run time offset in hours");
372 Print ("PDS-TDLP", "ID3-Oper", Prt_DS, pds->Oper,
373 TDLP_TableLookUp (TDLP_Oper_Table, sizeof (TDLP_Oper_Table),
374 pds->Oper));
375 sprintf (buffer, "%02d", pds->HH);
376 Print ("PDS-TDLP", "ID3-HH", Prt_SS, buffer,
377 "Number of hours between variables");
378 } else {
379 Print ("PDS-TDLP", "ID3-Oper", Prt_DS, pds->Oper,
380 TDLP_TableLookUp (TDLP_Oper_Table, sizeof (TDLP_Oper_Table),
381 pds->Oper));
382 }
383 sprintf (buffer, "%03d", pds->ttt);
384 Print ("PDS-TDLP", "ID3-ttt", Prt_SS, buffer, "Forecast Projection");
385
386 Print ("PDS-TDLP", "ID4-thresh", Prt_F, pds->thresh);
387 Print ("PDS-TDLP", "ID4-I", Prt_DS, pds->I,
388 TDLP_TableLookUp (TDLP_I_Table, sizeof (TDLP_I_Table), pds->I));
389 Print ("PDS-TDLP", "ID4-S", Prt_DS, pds->S,
390 TDLP_TableLookUp (TDLP_S_Table, sizeof (TDLP_S_Table), pds->S));
391 Print ("PDS-TDLP", "ID4-G", Prt_D, pds->G);
392 }
393
394 /*****************************************************************************
395 * TDLP_ElemSurfUnit() --
396 *
397 * Arthur Taylor / MDL
398 *
399 * PURPOSE
400 * Deal with element name, unit, and comment for TDLP items. Also deals
401 * with the surface information stored in ID2.
402 *
403 * ARGUMENTS
404 * pds = The TDLP Product Definition Section to parse. (Input)
405 * element = The resulting element name. (Output)
406 * unitName = The resulting unit name. (Output)
407 * comment = The resulting comment. (Output)
408 * shortFstLevel = The resulting short forecast level. (Output)
409 * longFstLevel = The resulting long forecast level. (Output)
410 *
411 * FILES/DATABASES: None
412 *
413 * RETURNS: void
414 *
415 * HISTORY
416 * 10/2004 Arthur Taylor (MDL): Created.
417 *
418 * NOTES
419 *****************************************************************************
420 */
TDLP_ElemSurfUnit(pdsTDLPType * pds,char ** element,char ** unitName,char ** comment,char ** shortFstLevel,char ** longFstLevel)421 static void TDLP_ElemSurfUnit (pdsTDLPType * pds, char **element,
422 char **unitName, char **comment,
423 char **shortFstLevel, char **longFstLevel)
424 {
425 char *ptr; /* Help guess unitName, and clean the elemName. */
426 char *ptr2; /* Help guess unitName, and clean the elemName. */
427
428 myAssert (*element == NULL);
429 myAssert (*unitName == NULL);
430 myAssert (*comment == NULL);
431 myAssert (*shortFstLevel == NULL);
432 myAssert (*longFstLevel == NULL);
433
434 *element = (char *) malloc (1 + strlen (pds->Descriptor) * sizeof (char));
435 strcpy (*element, pds->Descriptor);
436 (*element)[strlen (pds->Descriptor)] = '\0';
437
438 ptr = strchr (*element, '(');
439 if (ptr != NULL) {
440 ptr2 = strchr (ptr, ')');
441 *ptr2 = '\0';
442 if (strcmp (ptr + 1, "unofficial id") == 0) {
443 *unitName = (char *) malloc ((1 + 3) * sizeof (char));
444 strcpy (*unitName, "[-]");
445 } else {
446 reallocSprintf (unitName, "[%s]", ptr + 1);
447 }
448 /* Trim any parens from element. */
449 *ptr = '\0';
450 strTrimRight (*element, ' ');
451 } else {
452 *unitName = (char *) malloc ((1 + 3) * sizeof (char));
453 strcpy (*unitName, "[-]");
454 }
455 ptr = *element;
456 while (*ptr != '\0') {
457 if (*ptr == ' ') {
458 *ptr = '-';
459 }
460 ptr++;
461 }
462 strCompact (*element, '-');
463
464 reallocSprintf (comment, "%09ld-%09ld-%09ld-%ld %s", pds->ID1,
465 pds->ID2, pds->ID3, pds->ID4, *unitName);
466 reallocSprintf (shortFstLevel, "%09ld", pds->ID2);
467 reallocSprintf (longFstLevel, "%09ld", pds->ID2);
468 }
469
470 /*****************************************************************************
471 * TDLP_Inventory() --
472 *
473 * Arthur Taylor / MDL
474 *
475 * PURPOSE
476 * Parses the TDLP "Product Definition Section" for enough information to
477 * fill out the inventory data structure so we can do a simple inventory on
478 * the file in a similar way to how we did it for GRIB1 and GRIB2.
479 *
480 * ARGUMENTS
481 * fp = An opened TDLP file already at the correct message. (Input)
482 * tdlpLen = The total length of the TDLP message. (Input)
483 * inv = The inventory data structure that we need to fill. (Output)
484 *
485 * FILES/DATABASES: None
486 *
487 * RETURNS: int (could use errSprintf())
488 * 0 = OK
489 * -1 = tdlpLen is too small.
490 *
491 * HISTORY
492 * 10/2004 Arthur Taylor (MDL): Created.
493 *
494 * NOTES
495 * Speed improvements...
496 * 1) pds doen't need to be allocated each time.
497 * 2) Not all data is needed, do something like TDLP_RefTime
498 * 3) TDLP_ElemSurfUnit may be slow?
499 *****************************************************************************
500 */
TDLP_Inventory(DataSource & fp,sInt4 tdlpLen,inventoryType * inv)501 int TDLP_Inventory (DataSource &fp, sInt4 tdlpLen, inventoryType *inv)
502 {
503 sInt4 curLoc; /* Where we are in the current TDLP message. */
504 int i_temp;
505 uChar sectLen; /* Length of section. */
506 uChar *pds; /* The part of the message dealing with the PDS. */
507 pdsTDLPType pdsMeta; /* The pds parsed into a usable data structure. */
508 char f_gds; /* flag if there is a GDS section. */
509 char f_bms; /* flag if there is a BMS section. */
510 short int DSF; /* Decimal Scale Factor for unpacking the data. */
511 short int BSF; /* Binary Scale Factor for unpacking the data. */
512
513 curLoc = 8;
514 if ((i_temp = fp.DataSourceFgetc()) == EOF) {
515 errSprintf ("Ran out of file in PDS (TDLP_Inventory).\n");
516 return -1;
517 }
518 sectLen = (uChar) i_temp;
519 curLoc += sectLen;
520 if (curLoc > tdlpLen) {
521 errSprintf ("Ran out of data in PDS (TDLP_Inventory)\n");
522 return -1;
523 }
524 pds = (uChar *) malloc (sectLen * sizeof (uChar));
525 *pds = sectLen;
526 if (fp.DataSourceFread (pds + 1, sizeof (char), sectLen - 1) + 1 != sectLen) {
527 errSprintf ("Ran out of file.\n");
528 free (pds);
529 return -1;
530 }
531
532 if (ReadTDLPSect1 (pds, tdlpLen, &curLoc, &pdsMeta, &f_gds, &f_bms,
533 &DSF, &BSF) != 0) {
534 preErrSprintf ("Inside TDLP_Inventory\n");
535 free (pds);
536 return -1;
537 }
538 free (pds);
539
540 inv->element = NULL;
541 inv->unitName = NULL;
542 inv->comment = NULL;
543 free (inv->shortFstLevel);
544 inv->shortFstLevel = NULL;
545 free (inv->longFstLevel);
546 inv->longFstLevel = NULL;
547 TDLP_ElemSurfUnit (&pdsMeta, &(inv->element), &(inv->unitName),
548 &(inv->comment), &(inv->shortFstLevel),
549 &(inv->longFstLevel));
550
551 inv->refTime = pdsMeta.refTime;
552 inv->validTime = pdsMeta.refTime + pdsMeta.project;
553 inv->foreSec = pdsMeta.project;
554 return 0;
555 }
556
557 /*****************************************************************************
558 * TDLP_RefTime() --
559 *
560 * Arthur Taylor / MDL
561 *
562 * PURPOSE
563 * Return just the reference time of a TDLP message.
564 *
565 * ARGUMENTS
566 * fp = An opened TDLP file already at the correct message. (Input)
567 * tdlpLen = The total length of the TDLP message. (Input)
568 * refTime = The reference time of interest. (Output)
569 *
570 * FILES/DATABASES: None
571 *
572 * RETURNS: int (could use errSprintf())
573 * 0 = OK
574 * -1 = tdlpLen is too small.
575 *
576 * HISTORY
577 * 10/2004 Arthur Taylor (MDL): Created.
578 *
579 * NOTES
580 *****************************************************************************
581 */
TDLP_RefTime(DataSource & fp,sInt4 tdlpLen,double * refTime)582 int TDLP_RefTime (DataSource &fp, sInt4 tdlpLen, double *refTime)
583 {
584 int sectLen; /* Length of section. */
585 sInt4 curLoc; /* Where we are in the current TDLP message. */
586 int c_temp; /* Temporary variable for use with fgetc */
587 short int si_temp; /* Temporary variable. */
588 int year, t_year; /* The reference year, and a consistency test */
589 uChar month, t_month; /* The reference month, and a consistency test */
590 uChar day, t_day; /* The reference day, and a consistency test */
591 uChar hour, t_hour; /* The reference hour, and a consistency test */
592 uChar min; /* The reference minute */
593 sInt4 li_temp; /* Temporary variable. */
594
595 if ((sectLen = fp.DataSourceFgetc ()) == EOF)
596 goto error;
597 curLoc = 8 + sectLen;
598 if (curLoc > tdlpLen) {
599 errSprintf ("Ran out of data in PDS (TDLP_RefTime)\n");
600 return -1;
601 }
602 myAssert (sectLen <= 71);
603 if (sectLen < 39) {
604 errSprintf ("TDLP Section 1 is too small.\n");
605 return -1;
606 }
607
608 if ((c_temp = fp.DataSourceFgetc()) == EOF)
609 goto error;
610 if (FREAD_BIG (&si_temp, sizeof (short int), 1, fp) != 1)
611 goto error;
612 year = si_temp;
613 if ((c_temp = fp.DataSourceFgetc()) == EOF)
614 goto error;
615 month = c_temp;
616 if ((c_temp = fp.DataSourceFgetc()) == EOF)
617 goto error;
618 day = c_temp;
619 if ((c_temp = fp.DataSourceFgetc()) == EOF)
620 goto error;
621 hour = c_temp;
622 if ((c_temp = fp.DataSourceFgetc()) == EOF)
623 goto error;
624 min = c_temp;
625
626 if (FREAD_BIG (&li_temp, sizeof (sInt4), 1, fp) != 1)
627 goto error;
628 t_year = li_temp / 1000000L;
629 li_temp -= t_year * 1000000L;
630 t_month = li_temp / 10000L;
631 li_temp -= t_month * 10000L;
632 t_day = li_temp / 100;
633 t_hour = li_temp - t_day * 100;
634
635 if ((t_year != year) || (t_month != month) || (t_day != day) ||
636 (t_hour != hour)) {
637 errSprintf ("Error Inconsistent Times in TDLP_RefTime.\n");
638 return -1;
639 }
640 if (ParseTime (refTime, year, month, day, hour, min, 0) != 0) {
641 preErrSprintf ("Error In call to ParseTime in TDLP_RefTime.\n");
642 return -1;
643 }
644
645 /* Get to the end of the TDLP message. */
646 /* (inventory.c : GRIB2Inventory), is responsible for this. */
647 /* fseek (fp, gribLen - sectLen, SEEK_CUR); */
648 return 0;
649 error:
650 errSprintf ("Ran out of file in PDS (TDLP_RefTime).\n");
651 return -1;
652 }
653
654 /*****************************************************************************
655 * ReadTDLPSect2() --
656 *
657 * Arthur Taylor / MDL
658 *
659 * PURPOSE
660 * Parses the TDLP "Grid Definition Section" or section 2, filling out
661 * the gdsMeta data structure.
662 *
663 * ARGUMENTS
664 * gds = The compressed part of the message dealing with "GDS". (Input)
665 * tdlpLen = The total length of the TDLP message. (Input)
666 * curLoc = Current location in the TDLP message. (Output)
667 * gdsMeta = The filled out gdsMeta data structure. (Output)
668 *
669 * FILES/DATABASES: None
670 *
671 * RETURNS: int (could use errSprintf())
672 * 0 = OK
673 * -1 = tdlpLen is too small.
674 * -2 = unexpected values in gds.
675 *
676 * HISTORY
677 * 10/2004 Arthur Taylor (MDL): Created.
678 *
679 * NOTES
680 *****************************************************************************
681 */
ReadTDLPSect2(uChar * gds,sInt4 tdlpLen,sInt4 * curLoc,gdsType * gdsMeta)682 static int ReadTDLPSect2 (uChar *gds, sInt4 tdlpLen, sInt4 *curLoc,
683 gdsType *gdsMeta)
684 {
685 char sectLen; /* Length of section. */
686 int gridType; /* Which type of grid. (see enumerated types). */
687 sInt4 li_temp; /* Temporary variable. */
688
689 sectLen = *(gds++);
690 *curLoc += sectLen;
691 if (*curLoc > tdlpLen) {
692 errSprintf ("Ran out of data in GDS (TDLP Section 2)\n");
693 return -1;
694 }
695 myAssert (sectLen == 28);
696
697 gridType = *(gds++);
698 gdsMeta->Nx = GRIB_UNSIGN_INT2 (*gds, gds[1]);
699 gds += 2;
700 gdsMeta->Ny = GRIB_UNSIGN_INT2 (*gds, gds[1]);
701 gds += 2;
702 gdsMeta->lat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) / 10000.0;
703 gds += 3;
704 gdsMeta->lon1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) / 10000.0;
705 gdsMeta->lon1 = 360 - gdsMeta->lon1;
706 if (gdsMeta->lon1 < 0) {
707 gdsMeta->lon1 += 360;
708 }
709 if (gdsMeta->lon1 > 360) {
710 gdsMeta->lon1 -= 360;
711 }
712 gds += 3;
713 gdsMeta->orientLon = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) / 10000.0;
714 gdsMeta->orientLon = 360 - gdsMeta->orientLon;
715 if (gdsMeta->orientLon < 0) {
716 gdsMeta->orientLon += 360;
717 }
718 if (gdsMeta->orientLon > 360) {
719 gdsMeta->orientLon -= 360;
720 }
721 gds += 3;
722 MEMCPY_BIG (&li_temp, gds, sizeof (sInt4));
723 gdsMeta->Dx = li_temp / 1000.0;
724 gds += 4;
725 gdsMeta->meshLat = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) / 10000.0;
726 gds += 3;
727 if ((*gds != 0) || (gds[1] != 0) || (gds[2] != 0) || (gds[3] != 0) ||
728 (gds[4] != 0) || (gds[5] != 0)) {
729 errSprintf ("Error Reserved was not set to 0 in ReadTDLPSect2.\n");
730 return -1;
731 }
732
733 gdsMeta->numPts = gdsMeta->Nx * gdsMeta->Ny;
734 gdsMeta->f_sphere = 1;
735 gdsMeta->majEarth = 6371.2;
736 gdsMeta->minEarth = 6371.2;
737 gdsMeta->Dy = gdsMeta->Dx;
738 gdsMeta->resFlag = 0;
739 gdsMeta->center = 0;
740 gdsMeta->scan = 64;
741 gdsMeta->lat2 = 0;
742 gdsMeta->lon2 = 0;
743 gdsMeta->scaleLat1 = gdsMeta->meshLat;
744 gdsMeta->scaleLat2 = gdsMeta->meshLat;
745 gdsMeta->southLat = 0;
746 gdsMeta->southLon = 0;
747 switch (gridType) {
748 case TDLP_POLAR:
749 gdsMeta->projType = GS3_POLAR;
750 /* 4/24/2006 Added the following. */
751 gdsMeta->scaleLat1 = 90;
752 gdsMeta->scaleLat2 = 90;
753 break;
754 case TDLP_LAMBERT:
755 gdsMeta->projType = GS3_LAMBERT;
756 break;
757 case TDLP_MERCATOR:
758 gdsMeta->projType = GS3_MERCATOR;
759 /* 4/24/2006 Added the following. */
760 gdsMeta->scaleLat1 = 0;
761 gdsMeta->scaleLat2 = 0;
762 break;
763 default:
764 errSprintf ("Grid projection number is %d\n", gridType);
765 return -2;
766 }
767 return 0;
768 }
769
770 /*****************************************************************************
771 * ReadTDLPSect3() --
772 *
773 * Arthur Taylor / MDL
774 *
775 * PURPOSE
776 * Parses the TDLP "Bit Map Section" or section 3, filling out the bitmap
777 * as needed.
778 *
779 * ARGUMENTS
780 * bms = The compressed part of the message dealing with "BMS". (Input)
781 * tdlpLen = The total length of the TDLP message. (Input)
782 * curLoc = Current location in the TDLP message. (Output)
783 * bitmap = The extracted bitmap. (Output)
784 * NxNy = The total size of the grid. (Input)
785 *
786 * FILES/DATABASES: None
787 *
788 * RETURNS: int (could use errSprintf())
789 * 0 = OK
790 * -1 = tdlpLen is too small.
791 * -2 = unexpected values in bms.
792 *
793 * HISTORY
794 * 10/2004 Arthur Taylor (MDL): Created.
795 *
796 * NOTES
797 *****************************************************************************
798 */
ReadTDLPSect3(CPL_UNUSED uChar * bms,CPL_UNUSED sInt4 tdlpLen,CPL_UNUSED sInt4 * curLoc,CPL_UNUSED uChar * bitmap,CPL_UNUSED sInt4 NxNy)799 static int ReadTDLPSect3 (CPL_UNUSED uChar *bms,
800 CPL_UNUSED sInt4 tdlpLen,
801 CPL_UNUSED sInt4 *curLoc,
802 CPL_UNUSED uChar *bitmap,
803 CPL_UNUSED sInt4 NxNy)
804 {
805 errSprintf ("Bitmap data is Not Supported\n");
806 return -1;
807 }
808
809 /*****************************************************************************
810 * ReadTDLPSect4() --
811 *
812 * Arthur Taylor / MDL
813 *
814 * PURPOSE
815 * Unpacks the "Binary Data Section" or section 4.
816 *
817 * ARGUMENTS
818 * bds = The compressed part of the message dealing with "BDS". (Input)
819 * tdlpLen = The total length of the TDLP message. (Input)
820 * curLoc = Current location in the TDLP message. (Output)
821 * DSF = Decimal Scale Factor for unpacking the data. (Input)
822 * BSF = Binary Scale Factor for unpacking the data. (Input)
823 * data = The extracted grid. (Output)
824 * meta = The meta data associated with the grid (Input/Output)
825 * unitM = The M unit conversion value in equation y = Mx + B. (Input)
826 * unitB = The B unit conversion value in equation y = Mx + B. (Input)
827 *
828 * FILES/DATABASES: None
829 *
830 * RETURNS: int (could use errSprintf())
831 * 0 = OK
832 * -1 = gribLen is too small.
833 * -2 = unexpected values in bds.
834 *
835 * HISTORY
836 * 4/2003 Arthur Taylor (MDL/RSIS): Created
837 * 3/2004 AAT: Switched {# Pts * (# Bits in a Group) +
838 * # of unused bits != # of available bits} to a warning from an
839 * error.
840 * 2/2005 AAT: Found bug: memBitRead grp[i].bit was sizeof sInt4 instead
841 * of uChar.
842 * 2/2005 AAT: Second order diff, no miss value bug (lastData - 1) should
843 * be lastData.
844 * 2/2005 AAT: Added test to see if the number of bits needed matches the
845 * section length.
846 *
847 * NOTES
848 * 1) See metaparse.c : ParseGrid()
849 *****************************************************************************
850 */
ReadTDLPSect4(uChar * bds,sInt4 tdlpLen,sInt4 * curLoc,short int DSF,short int BSF,double * data,grib_MetaData * meta,CPL_UNUSED double unitM,CPL_UNUSED double unitB)851 static int ReadTDLPSect4 (uChar *bds, sInt4 tdlpLen, sInt4 *curLoc,
852 short int DSF, short int BSF, double *data,
853 grib_MetaData *meta,
854 CPL_UNUSED double unitM,
855 CPL_UNUSED double unitB)
856 {
857 uInt4 sectLen; /* Length in bytes of the current section. */
858 uChar f_notGridPnt; /* Not Grid point data? */
859 uChar f_complexPack; /* Complex packing? */
860 uChar f_sndOrder; /* Second order differencing? */
861 uChar f_primMiss; /* Primary missing value? */
862 uChar f_secMiss; /* Secondary missing value? */
863 uInt4 numPack; /* Number of points packed. */
864 sInt4 li_temp; /* Temporary variable. */
865 uInt4 uli_temp; /* Temporary variable. */
866 uChar bufLoc; /* Keeps track of where to start getting more data
867 * out of the packed data stream. */
868 uChar f_negative; /* used to help with signs of numbers. */
869 size_t numUsed; /* How many bytes were used in a given call to
870 * memBitRead. */
871 sInt4 origVal = 0; /* Original value. */
872 uChar mbit; /* # of bits for abs (first first order difference) */
873 sInt4 fstDiff = 0; /* First first order difference. */
874 sInt4 diff = 0; /* general first order difference. */
875 uChar nbit; /* # of bits for abs (overall min value) */
876 sInt4 minVal; /* Minimum value. */
877 size_t LX; /* Number of groups. */
878 uChar ibit; /* # of bits for group min values. */
879 uChar jbit; /* # of bits for # of bits for group. */
880 uChar kbit; /* # of bits for # values in a group. */
881 TDLGroupType *grp; /* Holds the info about each group. */
882 size_t i, j; /* Loop counters. */
883 uInt4 t_numPack; /* Used to total number of values in a group to check
884 * the numPack value. */
885 uInt4 t_numBits; /* Used to total number of bits used in the groups. */
886 uInt4 t_numBytes; /* Used to total number of bytes used to compare to
887 * sectLen. */
888 sInt4 maxVal; /* The max value in a group. */
889 uInt4 dataCnt; /* How many values (miss or othewise) we have read. */
890 uInt4 lastData; /* Index to last actual data. */
891 uInt4 numVal; /* # of actual (non-missing values) we have. */
892 double scale; /* Amount to scale values by. */
893 uInt4 dataInd; /* Index into data for this value (used to switch
894 * from a11..a1n,a2n..a21 to normal grid of
895 * a11..a1n,a21..a2n. */
896 uChar f_missing; /* Used to help with primary missing values, and the
897 * 0 bit possibility. */
898 #ifdef DEBUG
899 sInt4 t_UK1 = 0; /* Used to test theories about un defined values. */
900 sInt4 t_UK2 = 0; /* Used to test theories about un defined values. */
901 #endif
902
903 sectLen = GRIB_UNSIGN_INT3 (*bds, bds[1], bds[2]);
904 *curLoc += sectLen;
905 if (*curLoc > tdlpLen) {
906 errSprintf ("Ran out of data in BDS (TDLP Section 4)\n");
907 return -1;
908 }
909 bds += 3;
910 t_numBytes = 3;
911 f_notGridPnt = (GRIB2BIT_4 & *bds) ? 1 : 0;
912 f_complexPack = (GRIB2BIT_5 & *bds) ? 1 : 0;
913 f_sndOrder = (GRIB2BIT_6 & *bds) ? 1 : 0;
914 f_primMiss = (GRIB2BIT_7 & *bds) ? 1 : 0;
915 f_secMiss = (GRIB2BIT_8 & *bds) ? 1 : 0;
916
917 if (f_secMiss && (!f_primMiss)) {
918 errSprintf ("Secondary missing value without a primary!\n");
919 return -1;
920 }
921 if (f_complexPack) {
922 if (!f_sndOrder) {
923 meta->gridAttrib.packType = GS5_CMPLX;
924 } else {
925 meta->gridAttrib.packType = GS5_CMPLXSEC;
926 }
927 } else {
928 errSprintf ("Simple pack is not supported at this time.\n");
929 return -1;
930 }
931 bds++;
932 MEMCPY_BIG (&numPack, bds, sizeof (sInt4));
933 bds += 4;
934 t_numBytes += 5;
935 if (!f_notGridPnt) {
936 if (numPack != meta->gds.numPts) {
937 errSprintf ("Number packed %d != number of points %d\n", numPack,
938 meta->gds.numPts);
939 return -1;
940 }
941 }
942 meta->gridAttrib.DSF = DSF;
943 meta->gridAttrib.ESF = BSF;
944 meta->gridAttrib.fieldType = 0;
945 if (!f_primMiss) {
946 meta->gridAttrib.f_miss = 0;
947 meta->gridAttrib.missPri = 0;
948 meta->gridAttrib.missSec = 0;
949 } else {
950 MEMCPY_BIG (&li_temp, bds, sizeof (sInt4));
951 bds += 4;
952 t_numBytes += 4;
953 meta->gridAttrib.missPri = li_temp / 10000.0;
954 if (!f_secMiss) {
955 meta->gridAttrib.f_miss = 1;
956 meta->gridAttrib.missSec = 0;
957 } else {
958 MEMCPY_BIG (&li_temp, bds, sizeof (sInt4));
959 bds += 4;
960 t_numBytes += 4;
961 meta->gridAttrib.missSec = li_temp / 10000.0;
962 meta->gridAttrib.f_miss = 2;
963 }
964 }
965 /* Init the buffer location. */
966 bufLoc = 8;
967 /* The origValue and fstDiff are only present if sndOrder packed. */
968 if (f_sndOrder) {
969 memBitRead (&f_negative, sizeof (f_negative), bds, 1, &bufLoc,
970 &numUsed);
971 memBitRead (&uli_temp, sizeof (sInt4), bds, 31, &bufLoc, &numUsed);
972 myAssert (numUsed == 4);
973 bds += numUsed;
974 t_numBytes += numUsed;
975 origVal = (f_negative) ? -1 * uli_temp : uli_temp;
976 memBitRead (&mbit, sizeof (mbit), bds, 5, &bufLoc, &numUsed);
977 memBitRead (&f_negative, sizeof (f_negative), bds, 1, &bufLoc,
978 &numUsed);
979 myAssert (numUsed == 0);
980 myAssert ((mbit > 0) && (mbit < 32));
981 memBitRead (&uli_temp, sizeof (sInt4), bds, mbit, &bufLoc, &numUsed);
982 bds += numUsed;
983 t_numBytes += numUsed;
984 fstDiff = (f_negative) ? -1 * uli_temp : uli_temp;
985 }
986 memBitRead (&nbit, sizeof (nbit), bds, 5, &bufLoc, &numUsed);
987 bds += numUsed;
988 t_numBytes += numUsed;
989 memBitRead (&f_negative, sizeof (f_negative), bds, 1, &bufLoc, &numUsed);
990 bds += numUsed;
991 t_numBytes += numUsed;
992 myAssert ((nbit > 0) && (nbit < 32));
993 memBitRead (&uli_temp, sizeof (sInt4), bds, nbit, &bufLoc, &numUsed);
994 bds += numUsed;
995 t_numBytes += numUsed;
996 minVal = (f_negative) ? -1 * uli_temp : uli_temp;
997 memBitRead (&LX, sizeof (LX), bds, 16, &bufLoc, &numUsed);
998 bds += numUsed;
999 t_numBytes += numUsed;
1000 grp = (TDLGroupType *) malloc (LX * sizeof (TDLGroupType));
1001 memBitRead (&ibit, sizeof (ibit), bds, 5, &bufLoc, &numUsed);
1002 bds += numUsed;
1003 t_numBytes += numUsed;
1004 memBitRead (&jbit, sizeof (jbit), bds, 5, &bufLoc, &numUsed);
1005 /* Following assert is because it is the # of bits of # of bits. Which
1006 * means that # of bits of value that has a max of 64. */
1007 myAssert (jbit < 6);
1008 bds += numUsed;
1009 t_numBytes += numUsed;
1010 memBitRead (&kbit, sizeof (kbit), bds, 5, &bufLoc, &numUsed);
1011 bds += numUsed;
1012 t_numBytes += numUsed;
1013 myAssert (ibit < 33);
1014 for (i = 0; i < LX; i++) {
1015 if (ibit == 0) {
1016 grp[i].min = 0;
1017 } else {
1018 memBitRead (&(grp[i].min), sizeof (sInt4), bds, ibit, &bufLoc,
1019 &numUsed);
1020 bds += numUsed;
1021 t_numBytes += numUsed;
1022 }
1023 }
1024 myAssert (jbit < 8);
1025 for (i = 0; i < LX; i++) {
1026 if (jbit == 0) {
1027 grp[i].bit = 0;
1028 } else {
1029 myAssert (jbit <= sizeof (uChar) * 8);
1030 memBitRead (&(grp[i].bit), sizeof (uChar), bds, jbit, &bufLoc,
1031 &numUsed);
1032 bds += numUsed;
1033 t_numBytes += numUsed;
1034 }
1035 myAssert (grp[i].bit < 32);
1036 }
1037 myAssert (kbit < 33);
1038 t_numPack = 0;
1039 t_numBits = 0;
1040 for (i = 0; i < LX; i++) {
1041 if (kbit == 0) {
1042 grp[i].num = 0;
1043 } else {
1044 memBitRead (&(grp[i].num), sizeof (sInt4), bds, kbit, &bufLoc,
1045 &numUsed);
1046 bds += numUsed;
1047 t_numBytes += numUsed;
1048 }
1049 t_numPack += grp[i].num;
1050 t_numBits += grp[i].num * grp[i].bit;
1051 }
1052 if (t_numPack != numPack) {
1053 errSprintf ("Number packed %d != number of values in groups %d\n",
1054 numPack, t_numPack);
1055 free (grp);
1056 return -1;
1057 }
1058 if ((t_numBytes + ceil (t_numBits / 8.)) > sectLen) {
1059 errSprintf ("# bytes in groups %ld (%ld + %ld / 8) > sectLen %ld\n",
1060 (sInt4) (t_numBytes + ceil (t_numBits / 8.)),
1061 t_numBytes, t_numBits, sectLen);
1062 free (grp);
1063 return -1;
1064 }
1065 dataCnt = 0;
1066 dataInd = 0;
1067
1068 #ifdef DEBUG
1069 printf ("nbit %d, ibit %d, jbit %d, kbit %d\n", nbit, ibit, jbit, kbit);
1070 if ((t_numBytes + ceil (t_numBits / 8.)) != sectLen) {
1071 printf ("Caution: # bytes in groups %d (%d + %d / 8) != "
1072 "sectLen %d\n", (sInt4) (t_numBytes + ceil (t_numBits / 8.)),
1073 t_numBytes, t_numBits, sectLen);
1074 }
1075 #endif
1076
1077 /* Binary scale factor in TDLP has reverse sign from GRIB definition. */
1078 scale = pow (10.0, -1 * DSF) * pow (2.0, -1 * BSF);
1079
1080 meta->gridAttrib.f_maxmin = 0;
1081 /* Work with Second order complex packed data. */
1082 if (f_sndOrder) {
1083 /* *INDENT-OFF* */
1084 /* The algorithm appears to be:
1085 * Data: a1 a2 a3 a4 a5 ...
1086 * 1st diff: 0 b2 b3 b4 b5 ...
1087 * 2nd diff: UK1 UK2 c3 c4 c5 ...
1088 * We already know a1 and b2, and unpack a stream of UK1 UK2 c3 c4
1089 * The problem is that UK1, UK2 is undefined. Originally I thought
1090 * this was 0, or c3, but it appears that if b2 != 0, then
1091 * UK2 = c3 + 2 b2, and UK1 = c3 + 1 * b2, otherwise it appears that
1092 * UK1 == UK2, and typically UK1 == c3 (but not always). */
1093 /* *INDENT-ON* */
1094 myAssert (numPack >= 2);
1095 if (f_secMiss) {
1096 numVal = 0;
1097 lastData = 0;
1098 for (i = 0; i < LX; i++) {
1099 maxVal = (1 << grp[i].bit) - 1;
1100 for (j = 0; j < grp[i].num; j++) {
1101 /* signed int. */
1102 memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1103 &bufLoc, &numUsed);
1104 bds += numUsed;
1105 if (li_temp == maxVal) {
1106 data[dataInd] = meta->gridAttrib.missPri;
1107 } else if (li_temp == (maxVal - 1)) {
1108 data[dataInd] = meta->gridAttrib.missSec;
1109 } else {
1110 if (numVal > 1) {
1111 #ifdef DEBUG
1112 if (numVal == 2) {
1113 if (fstDiff != 0) {
1114 /*
1115 myAssert (t_UK1 == li_temp + fstDiff);
1116 myAssert (t_UK2 == li_temp + 2 * fstDiff);
1117 */
1118 } else {
1119 myAssert (t_UK1 == t_UK2);
1120 }
1121 }
1122 #endif
1123 diff += (li_temp + grp[i].min + minVal);
1124 data[dataInd] = data[lastData] + diff * scale;
1125 lastData = dataInd;
1126 } else if (numVal == 1) {
1127 data[dataInd] = (origVal + fstDiff) * scale;
1128 lastData = dataInd;
1129 diff = fstDiff;
1130 #ifdef DEBUG
1131 t_UK2 = li_temp;
1132 #endif
1133 } else {
1134 data[dataInd] = origVal * scale;
1135 #ifdef DEBUG
1136 t_UK1 = li_temp;
1137 #endif
1138 }
1139 numVal++;
1140 if (!meta->gridAttrib.f_maxmin) {
1141 meta->gridAttrib.min = data[dataInd];
1142 meta->gridAttrib.max = data[dataInd];
1143 meta->gridAttrib.f_maxmin = 1;
1144 } else {
1145 if (data[dataInd] < meta->gridAttrib.min) {
1146 meta->gridAttrib.min = data[dataInd];
1147 }
1148 if (data[dataInd] > meta->gridAttrib.max) {
1149 meta->gridAttrib.max = data[dataInd];
1150 }
1151 }
1152 }
1153 dataCnt++;
1154 dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1155 (2 * (dataCnt / meta->gds.Nx) + 1) *
1156 meta->gds.Nx - dataCnt - 1 : dataCnt;
1157 myAssert ((dataInd < numPack) || (dataCnt == numPack));
1158 }
1159 }
1160 } else if (f_primMiss) {
1161 numVal = 0;
1162 lastData = 0;
1163 for (i = 0; i < LX; i++) {
1164 maxVal = (1 << grp[i].bit) - 1;
1165 for (j = 0; j < grp[i].num; j++) {
1166 /* signed int. */
1167 memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1168 &bufLoc, &numUsed);
1169 bds += numUsed;
1170 f_missing = 0;
1171 if (li_temp == maxVal) {
1172 data[dataInd] = meta->gridAttrib.missPri;
1173 /* In the case of grp[i].bit == 0, if grp[i].min == 0, then
1174 * it is the missing value, otherwise regular value. Only
1175 * need to be concerned for primary missing values. */
1176 f_missing = 1;
1177 if ((grp[i].bit == 0) && (grp[i].min != 0)) {
1178 #ifdef DEBUG
1179 printf ("This doesn't happen often.\n");
1180 printf ("%d %d %d\n", (int) i, grp[i].bit, grp[i].min);
1181 #endif
1182 myAssert (1 == 2);
1183 f_missing = 0;
1184 }
1185 }
1186 if (!f_missing) {
1187 if (numVal > 1) {
1188 #ifdef DEBUG
1189 if (numVal == 2) {
1190 if (fstDiff != 0) {
1191 /*
1192 myAssert (t_UK1 == li_temp + fstDiff);
1193 myAssert (t_UK2 == li_temp + 2 * fstDiff);
1194 */
1195 } else {
1196 myAssert (t_UK1 == t_UK2);
1197 }
1198 }
1199 #endif
1200 diff += (li_temp + grp[i].min + minVal);
1201 data[dataInd] = data[lastData] + diff * scale;
1202 lastData = dataInd;
1203 } else if (numVal == 1) {
1204 data[dataInd] = (origVal + fstDiff) * scale;
1205 lastData = dataInd;
1206 diff = fstDiff;
1207 #ifdef DEBUG
1208 t_UK2 = li_temp;
1209 #endif
1210 } else {
1211 data[dataInd] = origVal * scale;
1212 #ifdef DEBUG
1213 t_UK1 = li_temp;
1214 #endif
1215 }
1216 numVal++;
1217 if (!meta->gridAttrib.f_maxmin) {
1218 meta->gridAttrib.min = data[dataInd];
1219 meta->gridAttrib.max = data[dataInd];
1220 meta->gridAttrib.f_maxmin = 1;
1221 } else {
1222 if (data[dataInd] < meta->gridAttrib.min) {
1223 meta->gridAttrib.min = data[dataInd];
1224 }
1225 if (data[dataInd] > meta->gridAttrib.max) {
1226 meta->gridAttrib.max = data[dataInd];
1227 }
1228 }
1229 }
1230 dataCnt++;
1231 dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1232 (2 * (dataCnt / meta->gds.Nx) + 1) *
1233 meta->gds.Nx - dataCnt - 1 : dataCnt;
1234 myAssert ((dataInd < numPack) || (dataCnt == numPack));
1235 }
1236 }
1237 } else {
1238 lastData = 0;
1239 for (i = 0; i < LX; i++) {
1240 for (j = 0; j < grp[i].num; j++) {
1241 /* signed int. */
1242 memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1243 &bufLoc, &numUsed);
1244 bds += numUsed;
1245 if (dataCnt > 1) {
1246 #ifdef DEBUG
1247 if (dataCnt == 2) {
1248 if (fstDiff != 0) {
1249 /*
1250 myAssert (t_UK1 == li_temp + fstDiff);
1251 myAssert (t_UK2 == li_temp + 2 * fstDiff);
1252 */
1253 } else {
1254 myAssert (t_UK1 == t_UK2);
1255 }
1256 }
1257 #endif
1258 diff += (li_temp + grp[i].min + minVal);
1259 data[dataInd] = data[lastData] + diff * scale;
1260 lastData = dataInd;
1261 } else if (dataCnt == 1) {
1262 data[dataInd] = (origVal + fstDiff) * scale;
1263 lastData = dataInd;
1264 diff = fstDiff;
1265 #ifdef DEBUG
1266 t_UK2 = li_temp;
1267 #endif
1268 } else {
1269 data[dataInd] = origVal * scale;
1270 #ifdef DEBUG
1271 t_UK1 = li_temp;
1272 #endif
1273 }
1274 #ifdef DEBUG
1275 /*
1276 if (i >= 4153) {
1277 */
1278 /*
1279 if ((data[dataInd] > 100) || (data[dataInd] < -100)) {
1280 */
1281 /*
1282 if ((diff > 50) || (diff < -50)) {
1283 printf ("li_temp :: %ld, diff = %ld\n", li_temp, diff);
1284 printf ("data[dataInd] :: %f\n", data[dataInd]);
1285 printf ("Group # %d element %d, grp[i].min %ld, "
1286 "grp[i].bit %d, minVal %ld\n", i, j, grp[i].min,
1287 grp[i].bit, minVal);
1288 }
1289 */
1290 #endif
1291 if (!meta->gridAttrib.f_maxmin) {
1292 meta->gridAttrib.min = data[dataInd];
1293 meta->gridAttrib.max = data[dataInd];
1294 meta->gridAttrib.f_maxmin = 1;
1295 } else {
1296 if (data[dataInd] < meta->gridAttrib.min) {
1297 meta->gridAttrib.min = data[dataInd];
1298 }
1299 if (data[dataInd] > meta->gridAttrib.max) {
1300 meta->gridAttrib.max = data[dataInd];
1301 }
1302 }
1303 dataCnt++;
1304 dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1305 (2 * (dataCnt / meta->gds.Nx) + 1) *
1306 meta->gds.Nx - dataCnt - 1 : dataCnt;
1307 myAssert ((dataInd < numPack) || (dataCnt == numPack));
1308 }
1309 }
1310 numVal = dataCnt;
1311 }
1312
1313 /* Work with regular complex packed data. */
1314 } else {
1315 #ifdef DEBUG
1316 /*
1317 printf ("Work with regular complex packed data\n");
1318 */
1319 #endif
1320 if (f_secMiss) {
1321 numVal = 0;
1322 for (i = 0; i < LX; i++) {
1323 maxVal = (1 << grp[i].bit) - 1;
1324 for (j = 0; j < grp[i].num; j++) {
1325 /* signed int. */
1326 memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1327 &bufLoc, &numUsed);
1328 bds += numUsed;
1329 if (li_temp == maxVal) {
1330 data[dataInd] = meta->gridAttrib.missPri;
1331 } else if (li_temp == (maxVal - 1)) {
1332 data[dataInd] = meta->gridAttrib.missSec;
1333 } else {
1334 data[dataInd] = (li_temp + grp[i].min + minVal) * scale;
1335 numVal++;
1336 if (!meta->gridAttrib.f_maxmin) {
1337 meta->gridAttrib.min = data[dataInd];
1338 meta->gridAttrib.max = data[dataInd];
1339 meta->gridAttrib.f_maxmin = 1;
1340 } else {
1341 if (data[dataInd] < meta->gridAttrib.min) {
1342 meta->gridAttrib.min = data[dataInd];
1343 }
1344 if (data[dataInd] > meta->gridAttrib.max) {
1345 meta->gridAttrib.max = data[dataInd];
1346 }
1347 }
1348 }
1349 dataCnt++;
1350 dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1351 (2 * (dataCnt / meta->gds.Nx) + 1) *
1352 meta->gds.Nx - dataCnt - 1 : dataCnt;
1353 myAssert ((dataInd < numPack) || (dataCnt == numPack));
1354 }
1355 }
1356 } else if (f_primMiss) {
1357 #ifdef DEBUG
1358 /*
1359 printf ("Work with primary missing data\n");
1360 */
1361 #endif
1362 numVal = 0;
1363 for (i = 0; i < LX; i++) {
1364 maxVal = (1 << grp[i].bit) - 1;
1365 for (j = 0; j < grp[i].num; j++) {
1366 memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1367 &bufLoc, &numUsed);
1368 bds += numUsed;
1369 f_missing = 0;
1370 if (li_temp == maxVal) {
1371 data[dataInd] = meta->gridAttrib.missPri;
1372 /* In the case of grp[i].bit == 0, if grp[i].min == 0, then
1373 * it is the missing value, otherwise regular value. Only
1374 * need to be concerned for primary missing values. */
1375 f_missing = 1;
1376 if ((grp[i].bit == 0) && (grp[i].min != 0)) {
1377 #ifdef DEBUG
1378 printf ("This doesn't happen often.\n");
1379 printf ("%d %d %d\n", (int) i, grp[i].bit, grp[i].min);
1380 myAssert (1 == 2);
1381 #endif
1382 f_missing = 0;
1383 }
1384 }
1385 if (!f_missing) {
1386 data[dataInd] = (li_temp + grp[i].min + minVal) * scale;
1387 numVal++;
1388 if (!meta->gridAttrib.f_maxmin) {
1389 meta->gridAttrib.min = data[dataInd];
1390 meta->gridAttrib.max = data[dataInd];
1391 meta->gridAttrib.f_maxmin = 1;
1392 } else {
1393 if (data[dataInd] < meta->gridAttrib.min) {
1394 meta->gridAttrib.min = data[dataInd];
1395 }
1396 if (data[dataInd] > meta->gridAttrib.max) {
1397 meta->gridAttrib.max = data[dataInd];
1398 }
1399 }
1400 }
1401 dataCnt++;
1402 dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1403 (2 * (dataCnt / meta->gds.Nx) + 1) *
1404 meta->gds.Nx - dataCnt - 1 : dataCnt;
1405 myAssert ((dataInd < numPack) || (dataCnt == numPack));
1406 }
1407 }
1408 } else {
1409 for (i = 0; i < LX; i++) {
1410 for (j = 0; j < grp[i].num; j++) {
1411 memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1412 &bufLoc, &numUsed);
1413 bds += numUsed;
1414 data[dataInd] = (li_temp + grp[i].min + minVal) * scale;
1415 if (!meta->gridAttrib.f_maxmin) {
1416 meta->gridAttrib.min = data[dataInd];
1417 meta->gridAttrib.max = data[dataInd];
1418 meta->gridAttrib.f_maxmin = 1;
1419 } else {
1420 if (data[dataInd] < meta->gridAttrib.min) {
1421 meta->gridAttrib.min = data[dataInd];
1422 }
1423 if (data[dataInd] > meta->gridAttrib.max) {
1424 meta->gridAttrib.max = data[dataInd];
1425 }
1426 }
1427 dataCnt++;
1428 dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1429 (2 * (dataCnt / meta->gds.Nx) + 1) *
1430 meta->gds.Nx - dataCnt - 1 : dataCnt;
1431 myAssert ((dataInd < numPack) || (dataCnt == numPack));
1432 }
1433 }
1434 numVal = dataCnt;
1435 }
1436 }
1437 meta->gridAttrib.numMiss = dataCnt - numVal;
1438 meta->gridAttrib.refVal = minVal * scale;
1439
1440 free (grp);
1441 return 0;
1442 }
1443
1444 /*****************************************************************************
1445 * ReadTDLPRecord() --
1446 *
1447 * Arthur Taylor / MDL
1448 *
1449 * PURPOSE
1450 * Reads in a TDLP message, and parses the data into various data
1451 * structures, for use with other code.
1452 *
1453 * ARGUMENTS
1454 * fp = An opened TDLP file already at the correct message. (Input)
1455 * TDLP_Data = The read in TDLP data. (Output)
1456 * tdlp_DataLen = Size of TDLP_Data. (Output)
1457 * meta = A filled in meta structure (Output)
1458 * IS = The structure containing all the arrays that the
1459 * unpacker uses (Output)
1460 * sect0 = Already read in section 0 data. (Input)
1461 * tdlpLen = Length of the TDLP message. (Input)
1462 * majEarth = Used to override the TDLP major axis of earth. (Input)
1463 * minEarth = Used to override the TDLP minor axis of earth. (Input)
1464 *
1465 * FILES/DATABASES:
1466 * An already opened file pointing to the desired TDLP message.
1467 *
1468 * RETURNS: int (could use errSprintf())
1469 * 0 = OK
1470 * -1 = Problems reading in the PDS.
1471 * -2 = Problems reading in the GDS.
1472 * -3 = Problems reading in the BMS.
1473 * -4 = Problems reading in the BDS.
1474 * -5 = Problems reading the closing section.
1475 *
1476 * HISTORY
1477 * 10/2004 Arthur Taylor (MDL): Created
1478 *
1479 * NOTES
1480 *****************************************************************************
1481 */
ReadTDLPRecord(DataSource & fp,double ** TDLP_Data,uInt4 * tdlp_DataLen,grib_MetaData * meta,IS_dataType * IS,sInt4 sect0[SECT0LEN_WORD],uInt4 tdlpLen,double majEarth,double minEarth)1482 int ReadTDLPRecord (DataSource &fp, double **TDLP_Data, uInt4 *tdlp_DataLen,
1483 grib_MetaData *meta, IS_dataType *IS,
1484 sInt4 sect0[SECT0LEN_WORD], uInt4 tdlpLen,
1485 double majEarth, double minEarth)
1486 {
1487 sInt4 nd5; /* Size of TDLP message rounded up to the nearest *
1488 * sInt4. */
1489 uChar *c_ipack; /* A char ptr to the message stored in IS->ipack */
1490 sInt4 curLoc; /* Current location in the GRIB message. */
1491 char f_gds; /* flag if there is a GDS section. */
1492 char f_bms; /* flag if there is a BMS section. */
1493 short int DSF; /* Decimal Scale Factor for unpacking the data. */
1494 short int BSF; /* Binary Scale Factor for unpacking the data. */
1495 double *tdlp_Data; /* A pointer to TDLP_Data for ease of manipulation. */
1496 double unitM = 1; /* M in y = Mx + B, for unit conversion. */
1497 double unitB = 0; /* B in y = Mx + B, for unit conversion. */
1498 sInt4 li_temp; /* Used to make sure section 5 is 7777. */
1499 size_t pad; /* Number of bytes to pad the message to get to the
1500 * correct byte boundary. */
1501 char buffer[24]; /* Read the trailing bytes in the TDLPack record. */
1502 uChar *bitmap; /* Would contain bitmap data if it was supported. */
1503
1504 /* Make room for entire message, and read it in. */
1505 /* nd5 needs to be tdlpLen in (sInt4) units rounded up. */
1506 nd5 = (tdlpLen + 3) / 4;
1507 if (nd5 > IS->ipackLen) {
1508 IS->ipackLen = nd5;
1509 IS->ipack = (sInt4 *) realloc ((void *) (IS->ipack),
1510 (IS->ipackLen) * sizeof (sInt4));
1511 }
1512 c_ipack = (uChar *) IS->ipack;
1513 /* Init last sInt4 to 0, to make sure that the padded bytes are 0. */
1514 IS->ipack[nd5 - 1] = 0;
1515 /* Init first 2 sInt4 to sect0. */
1516 memcpy (c_ipack, sect0, SECT0LEN_WORD * 2);
1517 /* Read in the rest of the message. */
1518 if (fp.DataSourceFread(c_ipack + SECT0LEN_WORD * 2, sizeof (char),
1519 (tdlpLen - SECT0LEN_WORD * 2)) + SECT0LEN_WORD * 2 != tdlpLen) {
1520 errSprintf ("Ran out of file\n");
1521 return -1;
1522 }
1523
1524 /* Preceeding was in degrib2, next part is specific to TDLP. */
1525 curLoc = 8;
1526 if (ReadTDLPSect1 (c_ipack + curLoc, tdlpLen, &curLoc, &(meta->pdsTdlp),
1527 &f_gds, &f_bms, &DSF, &BSF) != 0) {
1528 preErrSprintf ("Inside ReadGrib1Record\n");
1529 return -1;
1530 }
1531
1532 /* Figure out some basic stuff about the grid. */
1533 free (meta->element);
1534 meta->element = NULL;
1535 free (meta->unitName);
1536 meta->unitName = NULL;
1537 free (meta->comment);
1538 meta->comment = NULL;
1539 free (meta->shortFstLevel);
1540 meta->shortFstLevel = NULL;
1541 free (meta->longFstLevel);
1542 meta->longFstLevel = NULL;
1543 TDLP_ElemSurfUnit (&(meta->pdsTdlp), &(meta->element), &(meta->unitName),
1544 &(meta->comment), &(meta->shortFstLevel),
1545 &(meta->longFstLevel));
1546 meta->center = 7; /* US NWS, NCEP */
1547 meta->subcenter = 14; /* NWS Meteorological Development Laboratory */
1548
1549 /* strftime (meta->refTime, 20, "%Y%m%d%H%M",
1550 gmtime (&(meta->pdsTdlp.refTime)));
1551 */
1552 Clock_Print (meta->refTime, 20, meta->pdsTdlp.refTime, "%Y%m%d%H%M", 0);
1553
1554 /*
1555 validTime = meta->pdsTdlp.refTime + meta->pdsTdlp.project;
1556 strftime (meta->validTime, 20, "%Y%m%d%H%M", gmtime (&(validTime)));
1557 */
1558 Clock_Print (meta->validTime, 20, meta->pdsTdlp.refTime +
1559 meta->pdsTdlp.project, "%Y%m%d%H%M", 0);
1560
1561 meta->deltTime = meta->pdsTdlp.project;
1562
1563 /* Get the Grid Definition Section. */
1564 if (f_gds) {
1565 if (ReadTDLPSect2 (c_ipack + curLoc, tdlpLen, &curLoc,
1566 &(meta->gds)) != 0) {
1567 preErrSprintf ("Inside ReadGrib1Record\n");
1568 return -2;
1569 }
1570 } else {
1571 errSprintf ("Don't know how to handle vector data yet.\n");
1572 return -2;
1573 }
1574
1575 /* Allow over ride of the earth radii. */
1576 if ((majEarth > 6300) && (majEarth < 6400)) {
1577 if ((minEarth > 6300) && (minEarth < 6400)) {
1578 meta->gds.f_sphere = 0;
1579 meta->gds.majEarth = majEarth;
1580 meta->gds.minEarth = minEarth;
1581 if (majEarth == minEarth) {
1582 meta->gds.f_sphere = 1;
1583 }
1584 } else {
1585 meta->gds.f_sphere = 1;
1586 meta->gds.majEarth = majEarth;
1587 meta->gds.minEarth = majEarth;
1588 }
1589 }
1590
1591 /* Allocate memory for the grid. */
1592 if (meta->gds.numPts > *tdlp_DataLen) {
1593 *tdlp_DataLen = meta->gds.numPts;
1594 *TDLP_Data = (double *) realloc ((void *) (*TDLP_Data),
1595 (*tdlp_DataLen) * sizeof (double));
1596 }
1597 tdlp_Data = *TDLP_Data;
1598
1599 /* Get the Bit Map Section. */
1600 if (f_bms) {
1601 /* errSprintf ("Bitmap data is Not Supported\n");*/
1602 /* Need to allocate bitmap when this is implemented. */
1603 bitmap = NULL;
1604 ReadTDLPSect3 (c_ipack + curLoc, tdlpLen, &curLoc, bitmap,
1605 meta->gds.numPts);
1606 return -1;
1607 }
1608
1609 /* Read the GRID. */
1610 if (ReadTDLPSect4 (c_ipack + curLoc, tdlpLen, &curLoc, DSF, BSF,
1611 tdlp_Data, meta, unitM, unitB) != 0) {
1612 preErrSprintf ("Inside ReadTDLPRecord\n");
1613 return -4;
1614 }
1615
1616 /* Read section 5. If it is "7777" == 926365495 we are done. */
1617 memcpy (&li_temp, c_ipack + curLoc, 4);
1618 if (li_temp != 926365495L) {
1619 errSprintf ("Did not find the end of the message.\n");
1620 return -5;
1621 }
1622 curLoc += 4;
1623 /* Read the trailing part of the message. */
1624 /* TDLPack uses 4 bytes for FORTRAN record size, then another 8 bytes for
1625 * the size of the record (so FORTRAN can see it), then the data rounded
1626 * up to an 8 byte boundary, then a trailing 4 bytes for a final FORTRAN
1627 * record size. However it only stores in the gribLen the non-rounded
1628 * amount, so we need to take care of the rounding, and the trailing 4
1629 * bytes here. */
1630 pad = ((sInt4) (ceil (tdlpLen / 8.0))) * 8 - tdlpLen + 4;
1631 if (fp.DataSourceFread(buffer, sizeof (char), pad) != pad) {
1632 errSprintf ("Ran out of file\n");
1633 return -1;
1634 }
1635 return 0;
1636 }
1637
1638 /*****************************************************************************
1639 * TDL_ScaleData() --
1640 *
1641 * Arthur Taylor / MDL
1642 *
1643 * PURPOSE
1644 * Deal with scaling while excluding the primary and secondary missing
1645 * values. After this, dst should contain scaled data + primary or secondary
1646 * missing values
1647 * "tdlpack library"::pack2d.f line 257 or search for:
1648 "the above statement"
1649 *
1650 * ARGUMENTS
1651 * Src = The original data. (Input)
1652 * Dst = The scaled data. (Output)
1653 * numData = The number of elements in data. (Input)
1654 * DSF = Decimal Scale Factor for scaling the data. (Input)
1655 * BSF = Binary Scale Factor for scaling the data. (Input)
1656 * f_primMiss = Flag saying if we have a primary missing value (In/Out)
1657 * primMiss = primary missing value. (In/Out)
1658 * f_secMiss = Flag saying if we have a secondary missing value (In/Out)
1659 * secMiss = secondary missing value. (In/Out)
1660 * f_min = Flag saying if we have the minimum value. (Output)
1661 * min = minimum scaled value in the grid. (Output)
1662 *
1663 * FILES/DATABASES: None
1664 *
1665 * RETURNS: void
1666 *
1667 * HISTORY
1668 * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
1669 *
1670 * NOTES
1671 *****************************************************************************
1672 */
1673 #define SCALE_MISSING 10000L
TDL_ScaleData(double * Src,sInt4 * Dst,sInt4 numData,int DSF,int BSF,char * f_primMiss,double * primMiss,char * f_secMiss,double * secMiss,char * f_min,sInt4 * min)1674 static void TDL_ScaleData (double *Src, sInt4 *Dst, sInt4 numData,
1675 int DSF, int BSF, char *f_primMiss,
1676 double *primMiss, char *f_secMiss,
1677 double *secMiss, char *f_min, sInt4 *min)
1678 {
1679 sInt4 cnt;
1680 double *src = Src;
1681 sInt4 *dst = Dst;
1682 double scale = pow (10.0, -1 * DSF) * pow (2.0, -1 * BSF);
1683 char f_actualPrim = 0;
1684 char f_actualSec = 0;
1685 sInt4 li_primMiss = (sInt4) (*primMiss * SCALE_MISSING + .5);
1686 sInt4 li_secMiss = (sInt4) (*secMiss * SCALE_MISSING + .5);
1687
1688 *f_min = 0;
1689 for (cnt = 0; cnt < numData; cnt++) {
1690 if (((*f_primMiss) || (*f_secMiss)) && (*src == *primMiss)) {
1691 *(dst++) = li_primMiss;
1692 src++;
1693 f_actualPrim = 1;
1694 } else if ((*f_secMiss) && (*src == *secMiss)) {
1695 *(dst++) = li_secMiss;
1696 src++;
1697 f_actualSec = 1;
1698 } else {
1699 *(dst) = (long int) (floor ((*(src++) / scale) + .5));
1700 /* Check if scaled value == primary missing value. */
1701 if (((*f_primMiss) || (*f_secMiss)) && (*dst == li_primMiss)) {
1702 *dst = *dst - 1;
1703 }
1704 /* Check if scaled value == secondary missing value. */
1705 if ((*f_secMiss) && (*dst == li_secMiss)) {
1706 *dst = *dst - 1;
1707 /* Check if adjustment caused scaled value == primary missing. */
1708 if (*dst == li_primMiss) {
1709 *dst = *dst - 1;
1710 }
1711 }
1712 if (!(*f_min)) {
1713 *min = *dst;
1714 *f_min = 1;
1715 } else if (*min > *dst) {
1716 *min = *dst;
1717 }
1718 dst++;
1719 }
1720 }
1721 if ((*f_secMiss) && (!f_actualSec)) {
1722 *f_secMiss = 0;
1723 }
1724 if (((*f_secMiss) || (*f_primMiss)) && (!f_actualPrim)) {
1725 *f_primMiss = 0;
1726 /* Check consistency. */
1727 if (*f_secMiss) {
1728 *f_secMiss = 0;
1729 *f_primMiss = 1;
1730 *primMiss = *secMiss;
1731 }
1732 }
1733 }
1734
1735 /*****************************************************************************
1736 * TDL_ReorderGrid() --
1737 *
1738 * Arthur Taylor / MDL
1739 *
1740 * PURPOSE
1741 * Loop through the data, so that
1742 * data is: "a1,1 ... a1,n a2,n ... a2,1 ..."
1743 * instead of: "a1,1 ... a1,n a2,1 ... a2,n ..."
1744 *
1745 * ARGUMENTS
1746 * Src = The data. (Input/Output)
1747 * NX = The number of X values. (Input)
1748 * NY = The number of Y values. (Input)
1749 *
1750 * FILES/DATABASES: None
1751 *
1752 * RETURNS: void
1753 *
1754 * HISTORY
1755 * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
1756 *
1757 * NOTES
1758 *****************************************************************************
1759 */
TDL_ReorderGrid(sInt4 * Src,short int NX,short int NY)1760 static void TDL_ReorderGrid (sInt4 *Src, short int NX, short int NY)
1761 {
1762 int i, j;
1763 sInt4 *src1, *src2;
1764 sInt4 li_temp;
1765
1766 for (j = 1; j < NY; j += 2) {
1767 src1 = Src + j * NX;
1768 src2 = Src + (j + 1) * NX - 1;
1769 for (i = 0; i < (NX / 2); i++) {
1770 li_temp = *src1;
1771 *(src1++) = *src2;
1772 *(src2--) = li_temp;
1773 }
1774 }
1775 }
1776
1777 /*****************************************************************************
1778 * TDL_GetSecDiff() --
1779 *
1780 * Arthur Taylor / MDL
1781 *
1782 * PURPOSE
1783 * Get the second order difference where we have special values for missing,
1784 * and for actual data we have the following scheme.
1785 * Data: a1 a2 a3 a4 a5 ...
1786 * 1st diff: 0 b2 b3 b4 b5 ...
1787 * 2nd diff: UK1 UK2 c3 c4 c5 ...
1788 * where UK1 = c3 + b2, and UK2 = c3 + 2 * b2. Note: The choice of UK1, and
1789 * UK2 doesn't matter because of the following FORTRAN unpacking code:
1790 * IWORK(1)=IFIRST
1791 * IWORK(2)=IWORK(1)+IFOD
1792 * ISUM=IFOD
1793 * DO 385 K=3,IS4(3)
1794 * ISUM=IWORK(K)+ISUM
1795 * IWORK(K)=IWORK(K-1)+ISUM
1796 * 385 CONTINUE
1797 * So ISUM is a function of IWORK(3), not IWORK(1).
1798 *
1799 * ARGUMENTS
1800 * Data = The data. (Input)
1801 * numData = The number of elements in data. (Input)
1802 * SecDiff = The secondary differences of the data. (Output)
1803 * f_primMiss = Flag saying if we have a primary missing value (Input)
1804 * li_primMiss = Scaled primary missing value. (Input)
1805 * a1 = First non-missing value in the field. (Output)
1806 * b2 = First non-missing value in the 1st order delta field (Out)
1807 * min = The minimum value (Input).
1808 *
1809 * FILES/DATABASES: None
1810 *
1811 * RETURNS: int
1812 * 0 = Success
1813 * 1 = Couldn't find second differences (don't use).
1814 *
1815 * HISTORY
1816 * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
1817 *
1818 * NOTES
1819 *****************************************************************************
1820 */
TDL_GetSecDiff(sInt4 * Data,int numData,sInt4 * SecDiff,char f_primMiss,sInt4 li_primMiss,sInt4 * a1,sInt4 * b2,sInt4 * min)1821 static int TDL_GetSecDiff (sInt4 *Data, int numData, sInt4 *SecDiff,
1822 char f_primMiss, sInt4 li_primMiss,
1823 sInt4 *a1, sInt4 *b2, sInt4 *min)
1824 {
1825 int i;
1826 char f_min = 0;
1827 sInt4 last = 0, before_last = 0;
1828 int a1Index = -1;
1829 int a2Index = -1;
1830
1831 if (numData < 3) {
1832 return 1;
1833 }
1834 if (f_primMiss) {
1835 for (i = 0; i < numData; i++) {
1836 if (Data[i] == li_primMiss) {
1837 SecDiff[i] = li_primMiss;
1838 } else if (a1Index == -1) {
1839 a1Index = i;
1840 *a1 = Data[a1Index];
1841 } else if (a2Index == -1) {
1842 a2Index = i;
1843 *b2 = Data[a2Index] - Data[a1Index];
1844 before_last = Data[a1Index];
1845 last = Data[a2Index];
1846 } else {
1847 SecDiff[i] = Data[i] - 2 * last + before_last;
1848 before_last = last;
1849 last = Data[i];
1850 if (!f_min) {
1851 /* Set the UK1, UK2 values. */
1852 *min = SecDiff[i];
1853 f_min = 1;
1854 SecDiff[a1Index] = SecDiff[i] + *b2;
1855 SecDiff[a2Index] = SecDiff[i] + 2 * (*b2);
1856 } else if (*min > SecDiff[i]) {
1857 *min = SecDiff[i];
1858 }
1859 }
1860 }
1861 if (!f_min) {
1862 return 1;
1863 }
1864 } else {
1865 *a1 = Data[0];
1866 *b2 = Data[1] - Data[0];
1867 for (i = 3; i < numData; i++) {
1868 SecDiff[i] = Data[i] - 2 * Data[i - 1] - Data[i - 2];
1869 if (i == 3) {
1870 *min = SecDiff[i];
1871 /* Set the UK1, UK2 values. */
1872 SecDiff[0] = SecDiff[i] + *b2;
1873 SecDiff[1] = SecDiff[i] + 2 * (*b2);
1874 } else if (*min > SecDiff[i]) {
1875 *min = SecDiff[i];
1876 }
1877 }
1878 }
1879 return 0;
1880 }
1881
1882 /*****************************************************************************
1883 * TDL_UseSecDiff_Prim() --
1884 *
1885 * Arthur Taylor / MDL
1886 *
1887 * PURPOSE
1888 * Checks if the average range of 2nd order differences < average range of
1889 * 0 order differnces, to determine if we should use second order differences
1890 * This deals with the case when we have primary missing values.
1891 *
1892 * ARGUMENTS
1893 * Data = The data. (Input)
1894 * numData = The number of elements in data. (Input)
1895 * SecDiff = The secondary differences of the data. (Input)
1896 * li_primMiss = Scaled primary missing value. (Input)
1897 * minGroup = The minimum group size. (Input)
1898 *
1899 * FILES/DATABASES: None
1900 *
1901 * RETURNS: int
1902 * 0 = Don't use 2nd order differences.
1903 * 1 = Use 2nd order differences.
1904 *
1905 * HISTORY
1906 * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
1907 *
1908 * NOTES
1909 *****************************************************************************
1910 */
TDL_UseSecDiff_Prim(sInt4 * Data,sInt4 numData,sInt4 * SecDiff,sInt4 li_primMiss,int minGroup)1911 static int TDL_UseSecDiff_Prim (sInt4 *Data, sInt4 numData,
1912 sInt4 *SecDiff, sInt4 li_primMiss,
1913 int minGroup)
1914 {
1915 int i, locCnt;
1916 int range0, range2;
1917 int tot0, tot2;
1918 char f_min;
1919 sInt4 min = 0, max = 0;
1920
1921 locCnt = 0;
1922 range0 = 0;
1923 tot0 = 0;
1924 f_min = 0;
1925 /* Compute scores for no differences */
1926 for (i = 0; i < numData; i++) {
1927 if (Data[i] != li_primMiss) {
1928 if (!f_min) {
1929 min = Data[i];
1930 max = Data[i];
1931 f_min = 1;
1932 } else {
1933 if (min > Data[i])
1934 min = Data[i];
1935 if (max < Data[i])
1936 max = Data[i];
1937 }
1938 }
1939 locCnt++;
1940 /* Fake a "group" by using the minimum group size. */
1941 if (locCnt == minGroup) {
1942 if (f_min) {
1943 range0 += (max - min);
1944 tot0++;
1945 f_min = 0;
1946 }
1947 locCnt = 0;
1948 }
1949 }
1950 if (locCnt != 0) {
1951 range0 += (max - min);
1952 tot0++;
1953 }
1954
1955 locCnt = 0;
1956 range2 = 0;
1957 tot2 = 0;
1958 f_min = 0;
1959 /* Compute scores for second order differences */
1960 for (i = 0; i < numData; i++) {
1961 if (SecDiff[i] != li_primMiss) {
1962 if (!f_min) {
1963 min = SecDiff[i];
1964 max = SecDiff[i];
1965 f_min = 1;
1966 } else {
1967 if (min > SecDiff[i])
1968 min = SecDiff[i];
1969 if (max < SecDiff[i])
1970 max = SecDiff[i];
1971 }
1972 }
1973 locCnt++;
1974 /* Fake a "group" by using the minimum group size. */
1975 if (locCnt == minGroup) {
1976 if (f_min) {
1977 range2 += (max - min);
1978 tot2++;
1979 f_min = 0;
1980 }
1981 locCnt = 0;
1982 }
1983 }
1984 if (locCnt != 0) {
1985 range2 += (max - min);
1986 tot2++;
1987 }
1988
1989 /* Compare average group size of no differencing to second order. */
1990 if ((range0 / (tot0 + 0.0)) <= (range2 / (tot2 + 0.0))) {
1991 return 0;
1992 } else {
1993 return 1;
1994 }
1995 }
1996
1997 /*****************************************************************************
1998 * TDL_UseSecDiff() --
1999 *
2000 * Arthur Taylor / MDL
2001 *
2002 * PURPOSE
2003 * Checks if the average range of 2nd order differences < average range of
2004 * 0 order differnces, to determine if we should use second order differences
2005 *
2006 * ARGUMENTS
2007 * Data = The data. (Input)
2008 * numData = The number of elements in data. (Input)
2009 * SecDiff = The secondary differences of the data. (Input)
2010 * minGroup = The minimum group size. (Input)
2011 *
2012 * FILES/DATABASES: None
2013 *
2014 * RETURNS: int
2015 * 0 = Don't use 2nd order differences.
2016 * 1 = Use 2nd order differences.
2017 *
2018 * HISTORY
2019 * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2020 *
2021 * NOTES
2022 *****************************************************************************
2023 */
TDL_UseSecDiff(sInt4 * Data,sInt4 numData,sInt4 * SecDiff,int minGroup)2024 static int TDL_UseSecDiff (sInt4 *Data, sInt4 numData,
2025 sInt4 *SecDiff, int minGroup)
2026 {
2027 int i, locCnt;
2028 int range0, range2;
2029 int tot0, tot2;
2030 sInt4 min = 0, max = 0;
2031
2032 locCnt = 0;
2033 range0 = 0;
2034 tot0 = 0;
2035 /* Compute scores for no differences */
2036 for (i = 0; i < numData; i++) {
2037 if (locCnt == 0) {
2038 min = Data[i];
2039 max = Data[i];
2040 } else {
2041 if (min > Data[i])
2042 min = Data[i];
2043 if (max < Data[i])
2044 max = Data[i];
2045 }
2046 locCnt++;
2047 /* Fake a "group" by using the minimum group size. */
2048 if (locCnt == minGroup) {
2049 range0 += (max - min);
2050 tot0++;
2051 locCnt = 0;
2052 }
2053 }
2054 if (locCnt != 0) {
2055 range0 += (max - min);
2056 tot0++;
2057 }
2058
2059 locCnt = 0;
2060 range2 = 0;
2061 tot2 = 0;
2062 /* Compute scores for second order differences */
2063 for (i = 0; i < numData; i++) {
2064 if (locCnt == 0) {
2065 min = SecDiff[i];
2066 max = SecDiff[i];
2067 } else {
2068 if (min > SecDiff[i])
2069 min = SecDiff[i];
2070 if (max < SecDiff[i])
2071 max = SecDiff[i];
2072 }
2073 locCnt++;
2074 /* Fake a "group" by using the minimum group size. */
2075 if (locCnt == minGroup) {
2076 range2 += (max - min);
2077 tot2++;
2078 locCnt = 0;
2079 }
2080 }
2081 if (locCnt != 0) {
2082 range2 += (max - min);
2083 tot2++;
2084 }
2085
2086 /* Compare average group size of no differencing to second order. */
2087 if ((range0 / (tot0 + 0.0)) <= (range2 / (tot2 + 0.0))) {
2088 return 0;
2089 } else {
2090 return 1;
2091 }
2092 }
2093
2094 /*****************************************************************************
2095 * power() --
2096 *
2097 * Arthur Taylor / MDL
2098 *
2099 * PURPOSE
2100 * Calculate the number of bits required to store a given positive number.
2101 *
2102 * ARGUMENTS
2103 * val = The number to store (Input)
2104 * extra = number of slots to allocate for prim/sec missing values (Input)
2105 *
2106 * FILES/DATABASES: None
2107 *
2108 * RETURNS: int
2109 * The number of bits needed to store this number.
2110 *
2111 * HISTORY
2112 * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2113 *
2114 * NOTES
2115 *****************************************************************************
2116 */
power(uInt4 val,int extra)2117 static int power (uInt4 val, int extra)
2118 {
2119 int i;
2120
2121 val += extra;
2122 if (val == 0) {
2123 return 1;
2124 }
2125 for (i = 0; val != 0; i++) {
2126 val = val >> 1;
2127 }
2128 return i;
2129 }
2130
2131 /*****************************************************************************
2132 * findMaxMin2() --
2133 *
2134 * Arthur Taylor / MDL
2135 *
2136 * PURPOSE
2137 * Find the min/max value between start/stop index values in the data
2138 * Assuming primary and secondary missing values.
2139 *
2140 * ARGUMENTS
2141 * Data = The data. (Input)
2142 * start = The starting index in data (Input)
2143 * stop = The stoping index in data (Input)
2144 * li_primMiss = scaled primary missing value (Input)
2145 * li_secMiss = scaled secondary missing value (Input)
2146 * min = The min value found (Output)
2147 * max = The max value found (Output)
2148 *
2149 * FILES/DATABASES: None
2150 *
2151 * RETURNS: char
2152 * Flag if min/max are valid.
2153 *
2154 * HISTORY
2155 * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2156 *
2157 * NOTES
2158 *****************************************************************************
2159 */
findMaxMin2(sInt4 * Data,int start,int stop,sInt4 li_primMiss,sInt4 li_secMiss,sInt4 * min,sInt4 * max)2160 static char findMaxMin2 (sInt4 *Data, int start, int stop,
2161 sInt4 li_primMiss, sInt4 li_secMiss,
2162 sInt4 *min, sInt4 *max)
2163 {
2164 char f_min = 0; /* Flag if we found the max/min values */
2165 int i; /* Loop counter. */
2166
2167 *max = *min = Data[start];
2168 for (i = start; i < stop; i++) {
2169 if ((Data[i] != li_secMiss) && (Data[i] != li_primMiss)) {
2170 if (!f_min) {
2171 *max = Data[i];
2172 *min = Data[i];
2173 f_min = 1;
2174 } else {
2175 if (*max < Data[i]) {
2176 *max = Data[i];
2177 } else if (*min > Data[i]) {
2178 *min = Data[i];
2179 }
2180 }
2181 }
2182 }
2183 return f_min;
2184 }
2185
2186 /*****************************************************************************
2187 * findMaxMin1() --
2188 *
2189 * Arthur Taylor / MDL
2190 *
2191 * PURPOSE
2192 * Find the min/max value between start/stop index values in the data
2193 * Assuming primary missing values.
2194 *
2195 * ARGUMENTS
2196 * Data = The data. (Input)
2197 * start = The starting index in data (Input)
2198 * stop = The stoping index in data (Input)
2199 * li_primMiss = scaled primary missing value (Input)
2200 * min = The min value found (Output)
2201 * max = The max value found (Output)
2202 *
2203 * FILES/DATABASES: None
2204 *
2205 * RETURNS: char
2206 * Flag if min/max are valid.
2207 *
2208 * HISTORY
2209 * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2210 *
2211 * NOTES
2212 *****************************************************************************
2213 */
findMaxMin1(sInt4 * Data,int start,int stop,sInt4 li_primMiss,sInt4 * min,sInt4 * max)2214 static char findMaxMin1 (sInt4 *Data, int start, int stop,
2215 sInt4 li_primMiss, sInt4 *min, sInt4 *max)
2216 {
2217 char f_min = 0; /* Flag if we found the max/min values */
2218 int i; /* Loop counter. */
2219
2220 *max = *min = Data[start];
2221 for (i = start; i < stop; i++) {
2222 if (Data[i] != li_primMiss) {
2223 if (!f_min) {
2224 *max = Data[i];
2225 *min = Data[i];
2226 f_min = 1;
2227 } else {
2228 if (*max < Data[i]) {
2229 *max = Data[i];
2230 } else if (*min > Data[i]) {
2231 *min = Data[i];
2232 }
2233 }
2234 }
2235 }
2236 return f_min;
2237 }
2238
2239 /*****************************************************************************
2240 * findMaxMin0() --
2241 *
2242 * Arthur Taylor / MDL
2243 *
2244 * PURPOSE
2245 * Find the min/max value between start/stop index values in the data
2246 * Assuming no missing values.
2247 *
2248 * ARGUMENTS
2249 * Data = The data. (Input)
2250 * start = The starting index in data (Input)
2251 * stop = The stoping index in data (Input)
2252 * min = The min value found (Output)
2253 * max = The max value found (Output)
2254 *
2255 * FILES/DATABASES: None
2256 *
2257 * RETURNS: void
2258 *
2259 * HISTORY
2260 * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2261 *
2262 * NOTES
2263 *****************************************************************************
2264 */
findMaxMin0(sInt4 * Data,int start,int stop,sInt4 * min,sInt4 * max)2265 static void findMaxMin0 (sInt4 *Data, int start, int stop, sInt4 *min,
2266 sInt4 *max)
2267 {
2268 int i; /* Loop counter. */
2269
2270 *max = *min = Data[start];
2271 for (i = start + 1; i < stop; i++) {
2272 if (*max < Data[i]) {
2273 *max = Data[i];
2274 } else if (*min > Data[i]) {
2275 *min = Data[i];
2276 }
2277 }
2278 }
2279
2280 /*****************************************************************************
2281 * findGroup2() --
2282 *
2283 * Arthur Taylor / MDL
2284 *
2285 * PURPOSE
2286 * Find "split" so that the numbers between start and split are within
2287 * "range" of each other... stops if it reaches "stop".
2288 * Assumes primary and secondary missing values.
2289 *
2290 * ARGUMENTS
2291 * Data = The data. (Input)
2292 * start = The starting index in data (Input)
2293 * stop = The stoping index in data (Input)
2294 * li_primMiss = scaled primary missing value (Input)
2295 * li_secMiss = scaled secondary missing value (Input)
2296 * range = The range to use (Input)
2297 * split = The first index that is out of the range (Output)
2298 * min = The min value for the group. (Output)
2299 * max = The max value for the group. (Output)
2300 *
2301 * FILES/DATABASES: None
2302 *
2303 * RETURNS: void
2304 *
2305 * HISTORY
2306 * 1/2005 Arthur Taylor (MDL): Created
2307 *
2308 * NOTES
2309 *****************************************************************************
2310 */
findGroup2(sInt4 * Data,int start,int stop,sInt4 li_primMiss,sInt4 li_secMiss,sInt4 range,int * split,sInt4 * min,sInt4 * max)2311 static void findGroup2 (sInt4 *Data, int start, int stop,
2312 sInt4 li_primMiss, sInt4 li_secMiss,
2313 sInt4 range, int *split, sInt4 *min, sInt4 *max)
2314 {
2315 char f_min = 0; /* Flag if we found the max/min values */
2316 int i; /* Loop counter. */
2317
2318 *min = *max = 0;
2319 for (i = start; i < stop; i++) {
2320 if ((Data[i] != li_secMiss) && (Data[i] != li_primMiss)) {
2321 if (!f_min) {
2322 *max = *min = Data[i];
2323 f_min = 1;
2324 } else {
2325 if (*max < Data[i]) {
2326 if ((Data[i] - *min) > range) {
2327 *split = i;
2328 return;
2329 }
2330 *max = Data[i];
2331 } else if (*min > Data[i]) {
2332 if ((*max - Data[i]) > range) {
2333 *split = i;
2334 return;
2335 }
2336 *min = Data[i];
2337 }
2338 }
2339 }
2340 }
2341 *split = stop;
2342 }
2343
2344 /*****************************************************************************
2345 * findGroup1() --
2346 *
2347 * Arthur Taylor / MDL
2348 *
2349 * PURPOSE
2350 * Find "split" so that the numbers between start and split are within
2351 * "range" of each other... stops if it reaches "stop".
2352 * Assumes primary missing values.
2353 *
2354 * ARGUMENTS
2355 * Data = The data. (Input)
2356 * start = The starting index in data (Input)
2357 * stop = The stoping index in data (Input)
2358 * li_primMiss = scaled primary missing value (Input)
2359 * range = The range to use (Input)
2360 * split = The first index that is out of the range (Output)
2361 * min = The min value for the group. (Output)
2362 * max = The max value for the group. (Output)
2363 *
2364 * FILES/DATABASES: None
2365 *
2366 * RETURNS: void
2367 *
2368 * HISTORY
2369 * 1/2005 Arthur Taylor (MDL): Created
2370 *
2371 * NOTES
2372 *****************************************************************************
2373 */
findGroup1(sInt4 * Data,int start,int stop,sInt4 li_primMiss,sInt4 range,int * split,sInt4 * min,sInt4 * max)2374 static void findGroup1 (sInt4 *Data, int start, int stop,
2375 sInt4 li_primMiss, sInt4 range, int *split,
2376 sInt4 *min, sInt4 *max)
2377 {
2378 char f_min = 0; /* Flag if we found the max/min values */
2379 int i; /* Loop counter. */
2380
2381 *min = *max = 0;
2382 for (i = start; i < stop; i++) {
2383 if (Data[i] != li_primMiss) {
2384 if (!f_min) {
2385 *max = *min = Data[i];
2386 f_min = 1;
2387 } else {
2388 if (*max < Data[i]) {
2389 if ((Data[i] - *min) > range) {
2390 *split = i;
2391 return;
2392 }
2393 *max = Data[i];
2394 } else if (*min > Data[i]) {
2395 if ((*max - Data[i]) > range) {
2396 *split = i;
2397 return;
2398 }
2399 *min = Data[i];
2400 }
2401 }
2402 }
2403 }
2404 *split = stop;
2405 }
2406
2407 /*****************************************************************************
2408 * findGroup0() --
2409 *
2410 * Arthur Taylor / MDL
2411 *
2412 * PURPOSE
2413 * Find "split" so that the numbers between start and split are within
2414 * "range" of each other... stops if it reaches "stop".
2415 * Assumes no missing values.
2416 *
2417 * ARGUMENTS
2418 * Data = The data. (Input)
2419 * start = The starting index in data (Input)
2420 * stop = The stoping index in data (Input)
2421 * range = The range to use (Input)
2422 * split = The first index that is out of the range (Output)
2423 * min = The min value for the group. (Output)
2424 * max = The max value for the group. (Output)
2425 *
2426 * FILES/DATABASES: None
2427 *
2428 * RETURNS: void
2429 *
2430 * HISTORY
2431 * 1/2005 Arthur Taylor (MDL): Created
2432 *
2433 * NOTES
2434 *****************************************************************************
2435 */
findGroup0(sInt4 * Data,int start,int stop,sInt4 range,int * split,sInt4 * min,sInt4 * max)2436 static void findGroup0 (sInt4 *Data, int start, int stop,
2437 sInt4 range, int *split, sInt4 *min, sInt4 *max)
2438 {
2439 int i; /* Loop counter. */
2440
2441 *max = *min = Data[0];
2442 for (i = start + 1; i < stop; i++) {
2443 if (*max < Data[i]) {
2444 if ((Data[i] - *min) > range) {
2445 *split = i;
2446 return;
2447 }
2448 *max = Data[i];
2449 } else if (*min > Data[i]) {
2450 if ((*max - Data[i]) > range) {
2451 *split = i;
2452 return;
2453 }
2454 *min = Data[i];
2455 }
2456 }
2457 *split = stop;
2458 }
2459
2460 /*****************************************************************************
2461 * findGroupRev2() --
2462 *
2463 * Arthur Taylor / MDL
2464 *
2465 * PURPOSE
2466 * Find "split" so that the numbers between split and stop are within
2467 * "range" of each other... stops if it reaches "start".
2468 * Assumes primary and secondary missing values.
2469 *
2470 * ARGUMENTS
2471 * Data = The data. (Input)
2472 * start = The starting index in data (Input)
2473 * stop = The stoping index in data (Input)
2474 * li_primMiss = scaled primary missing value (Input)
2475 * li_secMiss = scaled secondary missing value (Input)
2476 * range = The range to use (Input)
2477 * split = The first index that is still in the range (Output)
2478 * min = The min value for the group. (Output)
2479 * max = The max value for the group. (Output)
2480 *
2481 * FILES/DATABASES: None
2482 *
2483 * RETURNS: void
2484 *
2485 * HISTORY
2486 * 1/2005 Arthur Taylor (MDL): Created
2487 *
2488 * NOTES
2489 *****************************************************************************
2490 */
findGroupRev2(sInt4 * Data,int start,int stop,sInt4 li_primMiss,sInt4 li_secMiss,sInt4 range,int * split,sInt4 * min,sInt4 * max)2491 static void findGroupRev2 (sInt4 *Data, int start, int stop,
2492 sInt4 li_primMiss, sInt4 li_secMiss,
2493 sInt4 range, int *split, sInt4 *min, sInt4 *max)
2494 {
2495 char f_min = 0; /* Flag if we found the max/min values */
2496 int i; /* Loop counter. */
2497
2498 *min = *max = 0;
2499 for (i = stop - 1; i >= start; i--) {
2500 if ((Data[i] != li_secMiss) && (Data[i] != li_primMiss)) {
2501 if (!f_min) {
2502 *max = *min = Data[i];
2503 f_min = 1;
2504 } else {
2505 if (*max < Data[i]) {
2506 if ((Data[i] - *min) > range) {
2507 *split = i + 1;
2508 return;
2509 }
2510 *max = Data[i];
2511 } else if (*min > Data[i]) {
2512 if ((*max - Data[i]) > range) {
2513 *split = i + 1;
2514 return;
2515 }
2516 *min = Data[i];
2517 }
2518 }
2519 }
2520 }
2521 *split = start;
2522 }
2523
2524 /*****************************************************************************
2525 * findGroupRev1() --
2526 *
2527 * Arthur Taylor / MDL
2528 *
2529 * PURPOSE
2530 * Find "split" so that the numbers between split and stop are within
2531 * "range" of each other... stops if it reaches "start".
2532 * Assumes primary missing values.
2533 *
2534 * ARGUMENTS
2535 * Data = The data. (Input)
2536 * start = The starting index in data (Input)
2537 * stop = The stoping index in data (Input)
2538 * li_primMiss = scaled primary missing value (Input)
2539 * range = The range to use (Input)
2540 * split = The first index that is still in the range (Output)
2541 * min = The min value for the group. (Output)
2542 * max = The max value for the group. (Output)
2543 *
2544 * FILES/DATABASES: None
2545 *
2546 * RETURNS: void
2547 *
2548 * HISTORY
2549 * 1/2005 Arthur Taylor (MDL): Created
2550 *
2551 * NOTES
2552 *****************************************************************************
2553 */
findGroupRev1(sInt4 * Data,int start,int stop,sInt4 li_primMiss,sInt4 range,int * split,sInt4 * min,sInt4 * max)2554 static void findGroupRev1 (sInt4 *Data, int start, int stop,
2555 sInt4 li_primMiss, sInt4 range, int *split,
2556 sInt4 *min, sInt4 *max)
2557 {
2558 char f_min = 0; /* Flag if we found the max/min values */
2559 int i; /* Loop counter. */
2560
2561 *min = *max = 0;
2562 for (i = stop - 1; i >= start; i--) {
2563 if (Data[i] != li_primMiss) {
2564 if (!f_min) {
2565 *max = *min = Data[i];
2566 f_min = 1;
2567 } else {
2568 if (*max < Data[i]) {
2569 if ((Data[i] - *min) > range) {
2570 *split = i + 1;
2571 return;
2572 }
2573 *max = Data[i];
2574 } else if (*min > Data[i]) {
2575 if ((*max - Data[i]) > range) {
2576 *split = i + 1;
2577 return;
2578 }
2579 *min = Data[i];
2580 }
2581 }
2582 }
2583 }
2584 *split = start;
2585 }
2586
2587 /*****************************************************************************
2588 * findGroupRev0() --
2589 *
2590 * Arthur Taylor / MDL
2591 *
2592 * PURPOSE
2593 * Find "split" so that the numbers between split and stop are within
2594 * "range" of each other... stops if it reaches "start".
2595 * Assumes no missing values.
2596 *
2597 * ARGUMENTS
2598 * Data = The data. (Input)
2599 * start = The starting index in data (Input)
2600 * stop = The stoping index in data (Input)
2601 * li_primMiss = scaled primary missing value (Input)
2602 * range = The range to use (Input)
2603 * split = The first index that is still in the range (Output)
2604 * min = The min value for the group. (Output)
2605 * max = The max value for the group. (Output)
2606 *
2607 * FILES/DATABASES: None
2608 *
2609 * RETURNS: void
2610 *
2611 * HISTORY
2612 * 1/2005 Arthur Taylor (MDL): Created
2613 *
2614 * NOTES
2615 *****************************************************************************
2616 */
findGroupRev0(sInt4 * Data,int start,int stop,sInt4 range,int * split,sInt4 * min,sInt4 * max)2617 static void findGroupRev0 (sInt4 *Data, int start, int stop,
2618 sInt4 range, int *split, sInt4 *min, sInt4 *max)
2619 {
2620 int i; /* Loop counter. */
2621
2622 *max = *min = Data[stop - 1];
2623 for (i = stop - 2; i >= start; i--) {
2624 if (*max < Data[i]) {
2625 if ((Data[i] - *min) > range) {
2626 *split = i + 1;
2627 return;
2628 }
2629 *max = Data[i];
2630 } else if (*min > Data[i]) {
2631 if ((*max - Data[i]) > range) {
2632 *split = i + 1;
2633 return;
2634 }
2635 *min = Data[i];
2636 }
2637 }
2638 *split = start;
2639 }
2640
2641 /*****************************************************************************
2642 * shiftGroup2() --
2643 *
2644 * Arthur Taylor / MDL
2645 *
2646 * PURPOSE
2647 * Find "split" so that the numbers between split and start1 are all inside
2648 * the range defined by max, min and bit. It allows max and min to change,
2649 * as long as it doesn't exceed the range defined by bit.
2650 * This is very similar to findGroupRev?() but here we already know
2651 * information about the min/max values, and are just expanding the group a
2652 * little, while the other one knew nothing about the group, and just wanted
2653 * a group of a given range.
2654 * Assumes primary and secondary missing values.
2655 *
2656 * ARGUMENTS
2657 * Data = The data. (Input)
2658 * start1 = The starting index in data (Input)
2659 * start2 = The starting index of the earlier group (ie don't go to any
2660 * earlier indicies than this. (Input)
2661 * li_primMiss = scaled primary missing value (Input)
2662 * li_secMiss = scaled secondary missing value (Input)
2663 * bit = The range we are allowed to store this in. (Input)
2664 * min = The min value for the group. (Input/Output)
2665 * max = The max value for the group. (Input/Output)
2666 * split = The first index that is still in the range (Output)
2667 *
2668 * FILES/DATABASES: None
2669 *
2670 * RETURNS: void
2671 *
2672 * HISTORY
2673 * 1/2005 Arthur Taylor (MDL): Created
2674 *
2675 * NOTES
2676 *****************************************************************************
2677 */
shiftGroup2(sInt4 * Data,int start1,int start2,sInt4 li_primMiss,sInt4 li_secMiss,int bit,sInt4 * min,sInt4 * max,size_t * split)2678 static void shiftGroup2 (sInt4 *Data, int start1, int start2,
2679 sInt4 li_primMiss, sInt4 li_secMiss, int bit,
2680 sInt4 *min, sInt4 *max, size_t *split)
2681 {
2682 int i; /* Loop counter. */
2683 int range; /* The range defined by bit. */
2684
2685 range = (int) (pow (2.0, bit) - 1) - 1;
2686 myAssert (start2 <= start1);
2687 for (i = start1; i >= start2; i--) {
2688 if ((Data[i] != li_primMiss) && (Data[i] != li_secMiss)) {
2689 if (Data[i] > *max) {
2690 if ((Data[i] - *min) <= range) {
2691 *max = Data[i];
2692 } else {
2693 *split = i + 1;
2694 return;
2695 }
2696 } else if (Data[i] < *min) {
2697 if ((*max - Data[i]) <= range) {
2698 *min = Data[i];
2699 } else {
2700 *split = i + 1;
2701 return;
2702 }
2703 }
2704 }
2705 }
2706 *split = start2;
2707 }
2708
2709 /*****************************************************************************
2710 * shiftGroup1() --
2711 *
2712 * Arthur Taylor / MDL
2713 *
2714 * PURPOSE
2715 * Find "split" so that the numbers between split and start1 are all inside
2716 * the range defined by max, min and bit. It allows max and min to change,
2717 * as long as it doesn't exceed the range defined by bit.
2718 * This is very similar to findGroupRev?() but here we already know
2719 * information about the min/max values, and are just expanding the group a
2720 * little, while the other one knew nothing about the group, and just wanted
2721 * a group of a given range.
2722 * Assumes primary missing values.
2723 *
2724 * ARGUMENTS
2725 * Data = The data. (Input)
2726 * start1 = The starting index in data (Input)
2727 * start2 = The starting index of the earlier group (ie don't go to any
2728 * earlier indicies than this. (Input)
2729 * li_primMiss = scaled primary missing value (Input)
2730 * bit = The range we are allowed to store this in. (Input)
2731 * min = The min value for the group. (Input/Output)
2732 * max = The max value for the group. (Input/Output)
2733 * split = The first index that is still in the range (Output)
2734 *
2735 * FILES/DATABASES: None
2736 *
2737 * RETURNS: void
2738 *
2739 * HISTORY
2740 * 1/2005 Arthur Taylor (MDL): Created
2741 *
2742 * NOTES
2743 *****************************************************************************
2744 */
shiftGroup1(sInt4 * Data,int start1,int start2,sInt4 li_primMiss,int bit,sInt4 * min,sInt4 * max,size_t * split)2745 static void shiftGroup1 (sInt4 *Data, int start1, int start2,
2746 sInt4 li_primMiss, int bit,
2747 sInt4 *min, sInt4 *max, size_t *split)
2748 {
2749 int i; /* Loop counter. */
2750 int range; /* The range defined by bit. */
2751
2752 range = (int) (pow (2.0, bit) - 1) - 1;
2753 myAssert (start2 <= start1);
2754 for (i = start1; i >= start2; i--) {
2755 if (Data[i] != li_primMiss) {
2756 if (Data[i] > *max) {
2757 if ((Data[i] - *min) <= range) {
2758 *max = Data[i];
2759 } else {
2760 *split = i + 1;
2761 return;
2762 }
2763 } else if (Data[i] < *min) {
2764 if ((*max - Data[i]) <= range) {
2765 *min = Data[i];
2766 } else {
2767 *split = i + 1;
2768 return;
2769 }
2770 }
2771 }
2772 }
2773 *split = start2;
2774 }
2775
2776 /*****************************************************************************
2777 * shiftGroup0() --
2778 *
2779 * Arthur Taylor / MDL
2780 *
2781 * PURPOSE
2782 * Find "split" so that the numbers between split and start1 are all inside
2783 * the range defined by max, min and bit. It allows max and min to change,
2784 * as long as it doesn't exceed the range defined by bit.
2785 * This is very similar to findGroupRev?() but here we already know
2786 * information about the min/max values, and are just expanding the group a
2787 * little, while the other one knew nothing about the group, and just wanted
2788 * a group of a given range.
2789 * Assumes no missing values.
2790 *
2791 * ARGUMENTS
2792 * Data = The data. (Input)
2793 * start1 = The starting index in data (Input)
2794 * start2 = The starting index of the earlier group (ie don't go to any
2795 * earlier indicies than this. (Input)
2796 * bit = The range we are allowed to store this in. (Input)
2797 * min = The min value for the group. (Input/Output)
2798 * max = The max value for the group. (Input/Output)
2799 * split = The first index that is still in the range (Output)
2800 *
2801 * FILES/DATABASES: None
2802 *
2803 * RETURNS: void
2804 *
2805 * HISTORY
2806 * 1/2005 Arthur Taylor (MDL): Created
2807 *
2808 * NOTES
2809 *****************************************************************************
2810 */
shiftGroup0(sInt4 * Data,int start1,int start2,int bit,sInt4 * min,sInt4 * max,size_t * split)2811 static void shiftGroup0 (sInt4 *Data, int start1, int start2, int bit,
2812 sInt4 *min, sInt4 *max, size_t *split)
2813 {
2814 int i; /* Loop counter. */
2815 int range; /* The range defined by bit. */
2816
2817 range = (int) (pow (2.0, bit) - 1) - 0;
2818 myAssert (start2 <= start1);
2819 for (i = start1; i >= start2; i--) {
2820 if (Data[i] > *max) {
2821 if ((Data[i] - *min) <= range) {
2822 *max = Data[i];
2823 } else {
2824 *split = i + 1;
2825 return;
2826 }
2827 } else if (Data[i] < *min) {
2828 if ((*max - Data[i]) <= range) {
2829 *min = Data[i];
2830 } else {
2831 *split = i + 1;
2832 return;
2833 }
2834 }
2835 }
2836 *split = start2;
2837 }
2838
2839 /*****************************************************************************
2840 * doSplit() --
2841 *
2842 * Arthur Taylor / MDL
2843 *
2844 * PURPOSE
2845 * Reduce the "bit range", and create groups that grab as much as they can
2846 * to the right. Then reduce those groups if they improve the score.
2847 *
2848 * ARGUMENTS
2849 * Data = The data. (Input)
2850 * numData = The number of elements in data. (Input)
2851 * G = The group to split. (Input)
2852 * lclGroup = The resulting groups (Output)
2853 * numLclGroup = The number of resulting groups. (Output)
2854 * f_primMiss = Flag if we have a primary missing value (Input)
2855 * li_primMiss = scaled primary missing value (Input)
2856 * f_secMiss = Flag if we have a secondary missing value (Input)
2857 * li_secMiss = scaled secondary missing value (Input)
2858 * xFactor = Estimate of cost (in bits) of a group. (Input)
2859 *
2860 * FILES/DATABASES: None
2861 *
2862 * RETURNS: void
2863 *
2864 * HISTORY
2865 * 1/2005 Arthur Taylor (MDL): Created.
2866 *
2867 * NOTES
2868 *****************************************************************************
2869 */
doSplit(sInt4 * Data,CPL_UNUSED int numData,TDLGroupType * G,TDLGroupType ** lclGroup,int * numLclGroup,char f_primMiss,sInt4 li_primMiss,char f_secMiss,sInt4 li_secMiss,int xFactor)2870 static void doSplit (sInt4 *Data,
2871 CPL_UNUSED int numData,
2872 TDLGroupType * G,
2873 TDLGroupType ** lclGroup, int *numLclGroup,
2874 char f_primMiss, sInt4 li_primMiss,
2875 char f_secMiss, sInt4 li_secMiss, int xFactor)
2876 {
2877 int start; /* Where to start the current group. */
2878 int range; /* The range to make the groups. */
2879 int final; /* One more than the last index in the group G. */
2880 int split; /* Where to split the group. */
2881 TDLGroupType G1; /* The current group to add. */
2882 TDLGroupType G2; /* The group if we evaporated the previous group. */
2883 int evaporate; /* How many groups we have "evaporated". */
2884 int i; /* Loop counter to help "evaporate" groups. */
2885 sInt4 scoreA; /* The original score for 2 groups */
2886 sInt4 scoreB; /* The new score (having evaporated a group) */
2887 int GroupLen; /* Actual alloc'ed group len. */
2888
2889 /* *INDENT-OFF* */
2890 /* The (pow (2, ..) -1) is because 2^n - 1 is max range of a group.
2891 * Example n = 1, 2^1 -1 = 1 range of (1,0) is 1.
2892 * Example n = 2, 2^2 -1 = 3 range of (3,2,1,0) is 3.
2893 * The G.bit - 1 is because we are trying to reduce the range. */
2894 /* *INDENT-ON* */
2895 range = (int) (pow (2.0, G->bit - 1) - 1) - (f_secMiss + f_primMiss);
2896 split = G->start;
2897 start = G->start;
2898 final = G->start + G->num;
2899 myAssert (final <= numData);
2900 *numLclGroup = 0;
2901 GroupLen = 1;
2902 *lclGroup = (TDLGroupType *) malloc (GroupLen * sizeof (TDLGroupType));
2903 while (split < final) {
2904 if (f_secMiss) {
2905 findGroup2 (Data, start, final, li_primMiss, li_secMiss, range,
2906 &split, &(G1.min), &(G1.max));
2907 } else if (f_primMiss) {
2908 findGroup1 (Data, start, final, li_primMiss, range, &split,
2909 &(G1.min), &(G1.max));
2910 } else {
2911 findGroup0 (Data, start, final, range, &split, &(G1.min), &(G1.max));
2912 }
2913 G1.bit = (char) power ((uInt4) (G1.max - G1.min),
2914 f_secMiss + f_primMiss);
2915 G1.num = split - start;
2916 G1.start = start;
2917 G1.f_trySplit = 1;
2918 G1.f_tryShift = 1;
2919 /* Test if we should add to previous group, or create a new group. */
2920 if (*numLclGroup == 0) {
2921 *numLclGroup = 1;
2922 (*lclGroup)[0] = G1;
2923 } else {
2924 G2.start = (*lclGroup)[*numLclGroup - 1].start;
2925 G2.num = (*lclGroup)[*numLclGroup - 1].num + G1.num;
2926 G2.min = ((*lclGroup)[*numLclGroup - 1].min < G1.min) ?
2927 (*lclGroup)[*numLclGroup - 1].min : G1.min;
2928 G2.max = ((*lclGroup)[*numLclGroup - 1].max > G1.max) ?
2929 (*lclGroup)[*numLclGroup - 1].max : G1.max;
2930 G2.bit = (char) power ((uInt4) (G2.max - G2.min),
2931 f_secMiss + f_primMiss);
2932 G2.f_trySplit = 1;
2933 G2.f_tryShift = 1;
2934 scoreA = ((*lclGroup)[*numLclGroup - 1].bit *
2935 (*lclGroup)[*numLclGroup - 1].num) + xFactor;
2936 scoreA += G1.bit * G1.num + xFactor;
2937 scoreB = G2.bit * G2.num + xFactor;
2938 if (scoreB < scoreA) {
2939 (*lclGroup)[*numLclGroup - 1] = G2;
2940 /* See if we can evaporate any of the old groups */
2941 evaporate = 0;
2942 for (i = *numLclGroup - 1; i > 0; i--) {
2943 G1.start = (*lclGroup)[i - 1].start;
2944 G1.num = (*lclGroup)[i].num + (*lclGroup)[i - 1].num;
2945 G1.min = ((*lclGroup)[i].min < (*lclGroup)[i - 1].min) ?
2946 (*lclGroup)[i].min : (*lclGroup)[i - 1].min;
2947 G1.max = ((*lclGroup)[i].max > (*lclGroup)[i - 1].max) ?
2948 (*lclGroup)[i].max : (*lclGroup)[i - 1].max;
2949 G1.bit = (char) power ((uInt4) (G1.max - G1.min),
2950 f_secMiss + f_primMiss);
2951 G1.f_trySplit = 1;
2952 G1.f_tryShift = 1;
2953 scoreA = (*lclGroup)[i].bit * (*lclGroup)[i].num + xFactor;
2954 scoreA += ((*lclGroup)[i - 1].bit * (*lclGroup)[i - 1].num +
2955 xFactor);
2956 scoreB = G1.bit * G1.num + xFactor;
2957 if (scoreB < scoreA) {
2958 evaporate++;
2959 (*lclGroup)[i - 1] = G1;
2960 } else {
2961 break;
2962 }
2963 }
2964 if (evaporate != 0) {
2965 *numLclGroup = *numLclGroup - evaporate;
2966 }
2967 } else {
2968 *numLclGroup = *numLclGroup + 1;
2969 if (*numLclGroup > GroupLen) {
2970 GroupLen = *numLclGroup;
2971 *lclGroup = (TDLGroupType *) realloc ((void *) *lclGroup,
2972 GroupLen *
2973 sizeof (TDLGroupType));
2974 }
2975 (*lclGroup)[*numLclGroup - 1] = G1;
2976 }
2977 }
2978 start = split;
2979 }
2980 }
2981
2982 /*****************************************************************************
2983 * doSplitRight() --
2984 *
2985 * Arthur Taylor / MDL
2986 *
2987 * PURPOSE
2988 * Break into two groups right has range n - 1, left has range n.
2989 *
2990 * ARGUMENTS
2991 * Data = The data. (Input)
2992 * numData = The number of elements in data. (Input)
2993 * G = The group to split. (Input)
2994 * G1 = The right most group (Output)
2995 * G2 = The remainder. (Output)
2996 * f_primMiss = Flag if we have a primary missing value (Input)
2997 * li_primMiss = scaled primary missing value (Input)
2998 * f_secMiss = Flag if we have a secondary missing value (Input)
2999 * li_secMiss = scaled secondary missing value (Input)
3000 *
3001 * FILES/DATABASES: None
3002 *
3003 * RETURNS: void
3004 *
3005 * HISTORY
3006 * 1/2005 Arthur Taylor (MDL): Created.
3007 *
3008 * NOTES
3009 *****************************************************************************
3010 */
doSplitRight(sInt4 * Data,CPL_UNUSED int numData,TDLGroupType * G,TDLGroupType * G1,TDLGroupType * G2,char f_primMiss,sInt4 li_primMiss,char f_secMiss,sInt4 li_secMiss)3011 static void doSplitRight (sInt4 *Data,
3012 CPL_UNUSED int numData,
3013 TDLGroupType * G,
3014 TDLGroupType * G1, TDLGroupType * G2,
3015 char f_primMiss, sInt4 li_primMiss,
3016 char f_secMiss, sInt4 li_secMiss)
3017 {
3018 int range; /* The range to make the right most group. */
3019 int final; /* One more than the last index in the group. */
3020 int split; /* Where to split the group. */
3021
3022 /* *INDENT-OFF* */
3023 /* The (pow (2, ..) -1) is because 2^n - 1 is max range of a group.
3024 * Example n = 1, 2^1 -1 = 1 range of (1,0) is 1.
3025 * Example n = 2, 2^2 -1 = 3 range of (3,2,1,0) is 3.
3026 * The G.bit - 1 is because we are trying to reduce the range. */
3027 /* *INDENT-ON* */
3028 range = (int) (pow (2.0, G->bit - 1) - 1) - (f_secMiss + f_primMiss);
3029 final = G->start + G->num;
3030 split = final;
3031 myAssert (final <= numData);
3032
3033 if (f_secMiss) {
3034 findGroupRev2 (Data, G->start, final, li_primMiss, li_secMiss, range,
3035 &split, &(G1->min), &(G1->max));
3036 findMaxMin2 (Data, G->start, split, li_primMiss, li_secMiss,
3037 &(G2->min), &(G2->max));
3038 } else if (f_primMiss) {
3039 findGroupRev1 (Data, G->start, final, li_primMiss, range, &split,
3040 &(G1->min), &(G1->max));
3041 findMaxMin1 (Data, G->start, split, li_primMiss, &(G2->min),
3042 &(G2->max));
3043 } else {
3044 findGroupRev0 (Data, G->start, final, range, &split, &(G1->min),
3045 &(G1->max));
3046 findMaxMin0 (Data, G->start, split, &(G2->min), &(G2->max));
3047 }
3048
3049 G1->bit = (char) power ((uInt4) (G1->max - G1->min),
3050 f_secMiss + f_primMiss);
3051 G2->bit = (char) power ((uInt4) (G2->max - G2->min),
3052 f_secMiss + f_primMiss);
3053 G1->start = split;
3054 G2->start = G->start;
3055 G1->num = final - split;
3056 G2->num = split - G->start;
3057 G1->f_trySplit = 1;
3058 G1->f_tryShift = 1;
3059 G2->f_trySplit = 1;
3060 G2->f_tryShift = 1;
3061 }
3062
3063 /*****************************************************************************
3064 * ComputeGroupSize() --
3065 *
3066 * Arthur Taylor / MDL
3067 *
3068 * PURPOSE
3069 * Compute the number of bits needed for the various elements of the groups
3070 * as well as the total number of bits needed.
3071 *
3072 * ARGUMENTS
3073 * group = Groups (Input)
3074 * numGroup = Number of groups (Input)
3075 * ibit = Number of bits needed for the minimum values of each group.
3076 * Find max absolute value of group mins. Note: all group mins
3077 * are positive (Output)
3078 * jbit = Number of bits needed for the number of bits for each group.
3079 * Find max absolute value of number of bits. (Output)
3080 * kbit = Number of bits needed for the number of values for each group.
3081 * Find max absolute value of number of values. (Output)
3082 *
3083 * FILES/DATABASES: None
3084 *
3085 * RETURNS: sInt4
3086 * number of bits needed by the groups
3087 *
3088 * HISTORY
3089 * 12/2004 Arthur Taylor (MDL): Updated from "tdlpack.c" in "C" tdlpack code
3090 *
3091 * NOTES
3092 *****************************************************************************
3093 */
ComputeGroupSize(TDLGroupType * group,int numGroup,size_t * ibit,size_t * jbit,size_t * kbit)3094 static sInt4 ComputeGroupSize (TDLGroupType * group, int numGroup,
3095 size_t *ibit, size_t *jbit, size_t *kbit)
3096 {
3097 int i; /* loop counter. */
3098 sInt4 ans = 0; /* The number of bits needed. */
3099 sInt4 maxMin = 0; /* The largest min value in the groups */
3100 uChar maxBit = 0; /* The largest needed bits in the groups */
3101 uInt4 maxNum = 0; /* The largest number of values in the groups. */
3102
3103 for (i = 0; i < numGroup; i++) {
3104 ans += group[i].bit * group[i].num;
3105 if (group[i].min > maxMin) {
3106 maxMin = group[i].min;
3107 }
3108 if (group[i].bit > maxBit) {
3109 maxBit = group[i].bit;
3110 }
3111 if (group[i].num > maxNum) {
3112 maxNum = group[i].num;
3113 }
3114 }
3115 /* This only works for pos numbers... */
3116 for (i = 0; (maxMin != 0); i++) {
3117 maxMin = maxMin >> 1;
3118 }
3119 /* Allow 0 bits for min. Assumes that decoder allows 0 bits */
3120 *ibit = i;
3121 /* This only works for pos numbers... */
3122 for (i = 0; (maxBit != 0); i++) {
3123 maxBit = maxBit >> 1;
3124 }
3125 /* Allow 0 bits for min. Assumes that decoder allows 0 bits */
3126 *jbit = i;
3127 /* This only works for pos numbers... */
3128 for (i = 0; (maxNum != 0); i++) {
3129 maxNum = maxNum >> 1;
3130 }
3131 /* Allow 0 bits for min. Assumes that decoder allows 0 bits */
3132 *kbit = i;
3133 ans += ((*ibit) + (*jbit) + (*kbit)) * numGroup;
3134 return ans;
3135 }
3136
3137 /*****************************************************************************
3138 * splitGroup() --
3139 *
3140 * Arthur Taylor / MDL
3141 *
3142 * PURPOSE
3143 * Tries to reduce (split) each group by 1 bit. It does this by:
3144 * A) reduce the "bit range", and create groups that grab as much as they
3145 * can to the right. Then reduce those groups if they improve the score.
3146 * B) reduce the bit range and grab the left most group only, leaving the
3147 * rest unchanged.
3148 * C) reduce the bit range and grab the right most group only, leaving the
3149 * rest unchanged.
3150 *
3151 * ARGUMENTS
3152 * Data = The data. (Input)
3153 * numData = The number of elements in data. (Input)
3154 * group = The groups using the "best minGroup. (Input)
3155 * numGroup = Number of groups (Input)
3156 * lclGroup = The local copy of the groups (Output)
3157 * numLclGroup = Number of local groups (Output)
3158 * f_primMiss = Flag if we have a primary missing value (Input)
3159 * li_primMiss = scaled primary missing value (Input)
3160 * f_secMiss = Flag if we have a secondary missing value (Input)
3161 * li_secMiss = scaled secondary missing value (Input)
3162 * xFactor = Estimate of cost (in bits) of a group. (Input)
3163 *
3164 * FILES/DATABASES: None
3165 *
3166 * RETURNS: void
3167 *
3168 * HISTORY
3169 * 1/2005 Arthur Taylor (MDL): Created.
3170 *
3171 * NOTES
3172 *****************************************************************************
3173 */
splitGroup(sInt4 * Data,int numData,TDLGroupType * group,int numGroup,TDLGroupType ** lclGroup,int * numLclGroup,char f_primMiss,sInt4 li_primMiss,char f_secMiss,sInt4 li_secMiss,size_t xFactor)3174 static int splitGroup (sInt4 *Data, int numData, TDLGroupType * group,
3175 int numGroup, TDLGroupType ** lclGroup,
3176 int *numLclGroup, char f_primMiss,
3177 sInt4 li_primMiss, char f_secMiss,
3178 sInt4 li_secMiss, size_t xFactor)
3179 {
3180 uInt4 minBit; /* The fewest # of bits, with no subdivision. */
3181 TDLGroupType *subGroup; /* The subgroups that we tried splitting the
3182 * primary group into. */
3183 int numSubGroup; /* The number of groups in subGroup. */
3184 sInt4 A_max; /* Max value of a given group. */
3185 sInt4 A_min; /* Min value of a given group. */
3186 sInt4 scoreA; /* The original score for a given group */
3187 sInt4 scoreB; /* The new score */
3188 int f_adjust = 0; /* Flag if group has changed. */
3189 int f_keep; /* Flag to keep the subgroup instead of original */
3190 int i; /* Loop counters */
3191 int sub; /* loop counter over the sub group. */
3192 int lclIndex; /* Used to help copy data from subGroup to answer. */
3193 int GroupLen; /* Actual alloc'ed group len. */
3194 int extra; /* Used to reduce number of allocs. */
3195
3196 /* Figure out how few bits a group can have without being able to further
3197 * divide it. */
3198 if (f_secMiss) {
3199 /* 11 = primMiss 10 = secMiss 01, 00 = data. */
3200 minBit = 2;
3201 } else if (f_primMiss) {
3202 /* 1 = primMiss 0 = data. */
3203 /* might try minBit = 1 here. */
3204 minBit = 1;
3205 } else {
3206 /* 1, 0 = data. */
3207 minBit = 1;
3208 }
3209
3210 *numLclGroup = 0;
3211 *lclGroup = (TDLGroupType *) malloc (numGroup * sizeof (TDLGroupType));
3212 GroupLen = numGroup;
3213 extra = 0;
3214 for (i = 0; i < numGroup; i++) {
3215 /* Check if we have already tried to split this group, or it has too
3216 * few members, or it doesn't have enough bits to split. If so, skip
3217 * this group. */
3218 if ((group[i].f_trySplit) && (group[i].num > xFactor) &&
3219 (group[i].bit > minBit)) {
3220 f_keep = 0;
3221 doSplit (Data, numData, &(group[i]), &subGroup, &numSubGroup,
3222 f_primMiss, li_primMiss, f_secMiss, li_secMiss, xFactor);
3223 if (numSubGroup != 1) {
3224 scoreA = group[i].bit * group[i].num + xFactor;
3225 scoreB = 0;
3226 for (sub = 0; sub < numSubGroup; sub++) {
3227 scoreB += subGroup[sub].bit * subGroup[sub].num + xFactor;
3228 }
3229 if (scoreB < scoreA) {
3230 f_keep = 1;
3231 } else if (numSubGroup > 2) {
3232 /* We can do "doSplitLeft" (which is breaking it into 2 groups
3233 * the first having range n - 1, the second having range n,
3234 * using what we know from doSplit. */
3235 subGroup[1].num = group[i].num - subGroup[0].num;
3236 if (f_secMiss) {
3237 findMaxMin2 (Data, subGroup[1].start,
3238 subGroup[1].start + subGroup[1].num,
3239 li_primMiss, li_secMiss, &A_min, &A_max);
3240 } else if (f_primMiss) {
3241 findMaxMin1 (Data, subGroup[1].start,
3242 subGroup[1].start + subGroup[1].num,
3243 li_primMiss, &A_min, &A_max);
3244 } else {
3245 findMaxMin0 (Data, subGroup[1].start,
3246 subGroup[1].start + subGroup[1].num,
3247 &A_min, &A_max);
3248 }
3249 subGroup[1].min = A_min;
3250 subGroup[1].max = A_max;
3251 subGroup[1].bit = power ((uInt4) (A_max - A_min),
3252 f_secMiss + f_primMiss);
3253 subGroup[1].f_trySplit = 1;
3254 subGroup[1].f_tryShift = 1;
3255 numSubGroup = 2;
3256 scoreB = subGroup[0].bit * subGroup[0].num + xFactor;
3257 scoreB += subGroup[1].bit * subGroup[1].num + xFactor;
3258 if (scoreB < scoreA) {
3259 f_keep = 1;
3260 }
3261 }
3262 }
3263 if (!f_keep) {
3264 if (numSubGroup == 1) {
3265 subGroup = (TDLGroupType *) realloc (subGroup,
3266 2 *
3267 sizeof (TDLGroupType));
3268 }
3269 numSubGroup = 2;
3270 doSplitRight (Data, numData, &(group[i]), &(subGroup[1]),
3271 &(subGroup[0]), f_primMiss, li_primMiss, f_secMiss,
3272 li_secMiss);
3273 scoreA = group[i].bit * group[i].num + xFactor;
3274 scoreB = subGroup[0].bit * subGroup[0].num + xFactor;
3275 scoreB += subGroup[1].bit * subGroup[1].num + xFactor;
3276 if (scoreB < scoreA) {
3277 f_keep = 1;
3278 }
3279 }
3280 if (f_keep) {
3281 lclIndex = *numLclGroup;
3282 *numLclGroup = *numLclGroup + numSubGroup;
3283 if (*numLclGroup > GroupLen) {
3284 GroupLen += extra;
3285 extra = 0;
3286 if (*numLclGroup > GroupLen) {
3287 GroupLen = *numLclGroup;
3288 }
3289 *lclGroup = (TDLGroupType *) realloc ((void *) *lclGroup,
3290 GroupLen *
3291 sizeof (TDLGroupType));
3292 } else {
3293 extra += numSubGroup - 1;
3294 }
3295 memcpy ((*lclGroup) + lclIndex, subGroup,
3296 numSubGroup * sizeof (TDLGroupType));
3297 f_adjust = 1;
3298 } else {
3299 *numLclGroup = *numLclGroup + 1;
3300 if (*numLclGroup > GroupLen) {
3301 GroupLen += extra;
3302 extra = 0;
3303 if (*numLclGroup > GroupLen) {
3304 GroupLen = *numLclGroup;
3305 }
3306 *lclGroup = (TDLGroupType *) realloc ((void *) *lclGroup,
3307 GroupLen *
3308 sizeof (TDLGroupType));
3309 }
3310 (*lclGroup)[*numLclGroup - 1] = group[i];
3311 (*lclGroup)[*numLclGroup - 1].f_trySplit = 0;
3312 }
3313 free (subGroup);
3314 subGroup = NULL;
3315 } else {
3316 *numLclGroup = *numLclGroup + 1;
3317 if (*numLclGroup > GroupLen) {
3318 GroupLen += extra;
3319 extra = 0;
3320 if (*numLclGroup > GroupLen) {
3321 GroupLen = *numLclGroup;
3322 }
3323 *lclGroup = (TDLGroupType *) realloc ((void *) *lclGroup,
3324 GroupLen *
3325 sizeof (TDLGroupType));
3326 }
3327 (*lclGroup)[*numLclGroup - 1] = group[i];
3328 (*lclGroup)[*numLclGroup - 1].f_trySplit = 0;
3329 }
3330 }
3331 myAssert (GroupLen == *numLclGroup);
3332 return f_adjust;
3333 }
3334
3335 /*****************************************************************************
3336 * shiftGroup() --
3337 *
3338 * Arthur Taylor / MDL
3339 *
3340 * PURPOSE
3341 * Tries to shift / join the groups together. It does this by first
3342 * calculating if a group should still exist. If it should, then it has
3343 * each group grab as much as it can to the left without increasing its "bit
3344 * range".
3345 *
3346 * ARGUMENTS
3347 * Data = The data. (Input)
3348 * numData = The number of elements in data. (Input)
3349 * group = The resulting groups. (Output)
3350 * numGroup = Number of groups (Output)
3351 * f_primMiss = Flag if we have a primary missing value (Input)
3352 * li_primMiss = scaled primary missing value (Input)
3353 * f_secMiss = Flag if we have a secondary missing value (Input)
3354 * li_secMiss = scaled secondary missing value (Input)
3355 * xFactor = Estimate of cost (in bits) of a group. (Input)
3356 *
3357 * FILES/DATABASES: None
3358 *
3359 * RETURNS: void
3360 *
3361 * HISTORY
3362 * 1/2005 Arthur Taylor (MDL): Created.
3363 *
3364 * NOTES
3365 *****************************************************************************
3366 */
shiftGroup(sInt4 * Data,CPL_UNUSED int numData,TDLGroupType ** Group,size_t * NumGroup,char f_primMiss,sInt4 li_primMiss,char f_secMiss,sInt4 li_secMiss,int xFactor)3367 static void shiftGroup (sInt4 *Data,
3368 CPL_UNUSED int numData,
3369 TDLGroupType ** Group,
3370 size_t *NumGroup, char f_primMiss, sInt4 li_primMiss,
3371 char f_secMiss, sInt4 li_secMiss, int xFactor)
3372 {
3373 TDLGroupType *group = (*Group); /* Local pointer to Group. */
3374 int numGroup = (*NumGroup); /* # elements in group. */
3375 int i, j; /* loop counters. */
3376 sInt4 A_max; /* Max value of a given group. */
3377 sInt4 A_min; /* Min value of a given group. */
3378 size_t begin; /* New start to the group. */
3379 sInt4 scoreA; /* The original score for group i and i - 1 */
3380 sInt4 scoreB; /* The new score for group i and i - 1 */
3381 TDLGroupType G1; /* The "new" group[i - 1]. */
3382 TDLGroupType G2; /* The "new" group[i]. */
3383 int evaporate = 0; /* number of groups that "evaporated" */
3384 int index; /* index while getting rid of "evaporated groups". */
3385
3386 for (i = numGroup - 1; i > 0; i--) {
3387 myAssert (group[i].num > 0);
3388 /* See if we can evaporate the group n - 1 */
3389 G1.start = group[i - 1].start;
3390 G1.num = group[i].num + group[i - 1].num;
3391 G1.min = (group[i].min < group[i - 1].min) ?
3392 group[i].min : group[i - 1].min;
3393 G1.max = (group[i].max > group[i - 1].max) ?
3394 group[i].max : group[i - 1].max;
3395 G1.bit = (char) power ((uInt4) (G1.max - G1.min),
3396 f_secMiss + f_primMiss);
3397 G1.f_trySplit = 1;
3398 G1.f_tryShift = 1;
3399 scoreA = group[i].bit * group[i].num + xFactor;
3400 scoreA += group[i - 1].bit * group[i - 1].num + xFactor;
3401 scoreB = G1.bit * G1.num + xFactor;
3402 if (scoreB < scoreA) {
3403 /* One of the groups evaporated. */
3404 evaporate++;
3405 group[i - 1] = G1;
3406 group[i].num = 0;
3407 /* See if that affects any of the previous groups. */
3408 for (j = i + 1; j < numGroup; j++) {
3409 if (group[j].num != 0) {
3410 G1.start = G1.start;
3411 G1.num = group[i - 1].num + group[j].num;
3412 G1.min = (group[i - 1].min < group[j].min) ?
3413 group[i - 1].min : group[j].min;
3414 G1.max = (group[i - 1].max > group[j].max) ?
3415 group[i - 1].max : group[j].max;
3416 G1.bit = (char) power ((uInt4) (G1.max - G1.min),
3417 f_secMiss + f_primMiss);
3418 G1.f_trySplit = 1;
3419 G1.f_tryShift = 1;
3420 scoreA = group[i - 1].bit * group[i - 1].num + xFactor;
3421 scoreA += group[j].bit * group[j].num + xFactor;
3422 scoreB = G1.bit * G1.num + xFactor;
3423 if (scoreB < scoreA) {
3424 evaporate++;
3425 group[i - 1] = G1;
3426 group[j].num = 0;
3427 } else {
3428 break;
3429 }
3430 }
3431 }
3432 } else if (group[i].f_tryShift) {
3433 /* Group did not evaporate, so do the "grabby" algorithm. */
3434 if ((group[i].bit != 0) && (group[i - 1].bit >= group[i].bit)) {
3435 if (f_secMiss) {
3436 A_max = group[i].max;
3437 A_min = group[i].min;
3438 shiftGroup2 (Data, group[i].start - 1, group[i - 1].start,
3439 li_primMiss, li_secMiss, group[i].bit, &A_min,
3440 &A_max, &begin);
3441 } else if (f_primMiss) {
3442 A_max = group[i].max;
3443 A_min = group[i].min;
3444 shiftGroup1 (Data, group[i].start - 1, group[i - 1].start,
3445 li_primMiss, group[i].bit, &A_min, &A_max,
3446 &begin);
3447 } else {
3448 A_max = group[i].max;
3449 A_min = group[i].min;
3450 shiftGroup0 (Data, group[i].start - 1, group[i - 1].start,
3451 group[i].bit, &A_min, &A_max, &begin);
3452 }
3453 if (begin != group[i].start) {
3454 /* Re-Calculate min/max of group[i - 1], since it could have
3455 * moved to group[i] */
3456 G1 = group[i - 1];
3457 G2 = group[i];
3458 G2.min = A_min;
3459 G2.max = A_max;
3460 G1.num -= (group[i].start - begin);
3461 if (f_secMiss) {
3462 findMaxMin2 (Data, G1.start, G1.start + G1.num,
3463 li_primMiss, li_secMiss, &A_min, &A_max);
3464 } else if (f_primMiss) {
3465 findMaxMin1 (Data, G1.start, G1.start + G1.num,
3466 li_primMiss, &A_min, &A_max);
3467 } else {
3468 findMaxMin0 (Data, G1.start, G1.start + G1.num,
3469 &A_min, &A_max);
3470 }
3471 if ((A_min != G1.min) || (A_max != G1.max)) {
3472 G1.min = A_min;
3473 G1.max = A_max;
3474 G1.bit = (char) power ((uInt4) (A_max - A_min),
3475 f_secMiss + f_primMiss);
3476 G1.f_trySplit = 1;
3477 G1.f_tryShift = 1;
3478 }
3479 G2.num += (group[i].start - begin);
3480 G2.start = begin;
3481 G2.f_trySplit = 1;
3482 G2.f_tryShift = 1;
3483 scoreA = group[i].bit * group[i].num + xFactor;
3484 scoreA += group[i - 1].bit * group[i - 1].num + xFactor;
3485 scoreB = G2.bit * G2.num + xFactor;
3486 if (G1.num != 0) {
3487 scoreB += G1.bit * G1.num + xFactor;
3488 }
3489 #ifdef DEBUG
3490 if (scoreB > scoreA) {
3491 printf ("Made score worse!\n");
3492 }
3493 #endif
3494 if (begin == group[i - 1].start) {
3495 /* Grabby algorithm evaporated a group. */
3496 evaporate++;
3497 myAssert (G1.num == 0);
3498 /* Switch the evaporating group to other side so we have
3499 * potential to continue evaporating the next group down. */
3500 group[i - 1] = G2;
3501 group[i] = G1;
3502 } else {
3503 group[i - 1] = G1;
3504 group[i] = G2;
3505 }
3506 } else {
3507 group[i].f_tryShift = 0;
3508 }
3509 }
3510 }
3511 }
3512 /* Loop through the grid removing the evaporated groups. */
3513 if (evaporate != 0) {
3514 index = 0;
3515 for (i = 0; i < numGroup; i++) {
3516 if (group[i].num != 0) {
3517 group[index] = group[i];
3518 index++;
3519 }
3520 }
3521 *NumGroup = numGroup - evaporate;
3522 *Group = (TDLGroupType *) realloc ((void *) (*Group),
3523 *NumGroup * sizeof (TDLGroupType));
3524 }
3525 }
3526
3527 /*****************************************************************************
3528 * GroupIt() --
3529 *
3530 * Arthur Taylor / MDL
3531 *
3532 * PURPOSE
3533 * Attempts to find groups for packing the data. It starts by preparing
3534 * the data, by removing the overall min value.
3535 *
3536 * Next it Creates any 0 bit groups (primary missing: missings are all 0
3537 * bit groups, const values are 0 bit groups. No missing: const values are
3538 * 0 bit groups.)
3539 *
3540 * Next it tries to reduce (split) each group by 1 bit. It does this by:
3541 * A) reduce the "bit range", and create groups that grab as much as they
3542 * can to the right. Then reduce those groups if they improve the score.
3543 * B) reduce the bit range and grab the left most group only, leaving the
3544 * rest unchanged.
3545 * C) reduce the bit range and grab the right most group only, leaving the
3546 * rest unchanged.
3547 *
3548 * Next it tries to shift / join those groups together. It does this by
3549 * first calculating if a group should still exist. If it should, then it
3550 * has each group grab as much as it can to the left without increasing its
3551 * "bit range".
3552 *
3553 * ARGUMENTS
3554 * OverallMin = The overall min value in the data. (Input)
3555 * Data = The data. (Input)
3556 * numData = The number of elements in data. (Input)
3557 * group = The resulting groups. (Output)
3558 * numGroup = Number of groups (Output)
3559 * f_primMiss = Flag if we have a primary missing value (Input)
3560 * li_primMiss = scaled primary missing value (Input)
3561 * f_secMiss = Flag if we have a secondary missing value (Input)
3562 * li_secMiss = scaled secondary missing value (Input)
3563 * groupSize = How many bytes the groups and data will take. (Output)
3564 * ibit = # of bits for largest minimum value in groups (Output)
3565 * jbit = # of bits for largest # bits in groups (Output)
3566 * kbit = # of bits for largest # values in groups (Output)
3567 *
3568 * FILES/DATABASES: None
3569 *
3570 * RETURNS: void
3571 *
3572 * HISTORY
3573 * 1/2005 Arthur Taylor (MDL): Created.
3574 *
3575 * NOTES
3576 * 1) Have not implemented const 0 bit groups for prim miss or no miss.
3577 * 2) MAX_GROUP_LEN (found experimentally) is useful to keep cost of groups
3578 * down.
3579 *****************************************************************************
3580 */
3581 #define MAX_GROUP_LEN 255
GroupIt(sInt4 OverallMin,sInt4 * Data,size_t numData,TDLGroupType ** group,size_t * numGroup,char f_primMiss,sInt4 li_primMiss,char f_secMiss,sInt4 li_secMiss,sInt4 * groupSize,size_t * ibit,size_t * jbit,size_t * kbit)3582 static void GroupIt (sInt4 OverallMin, sInt4 *Data, size_t numData,
3583 TDLGroupType ** group, size_t *numGroup, char f_primMiss,
3584 sInt4 li_primMiss, char f_secMiss,
3585 sInt4 li_secMiss, sInt4 *groupSize, size_t *ibit,
3586 size_t *jbit, size_t *kbit)
3587 {
3588 sInt4 A_max; /* Max value of a given group. */
3589 sInt4 A_min; /* Min value of a given group. */
3590 TDLGroupType G; /* Used to init the groups. */
3591 int f_adjust; /* Flag if we have changed the groups. */
3592 TDLGroupType *lclGroup; /* A temporary copy of the groups. */
3593 int numLclGroup; /* # of groups in lclGroup. */
3594 size_t xFactor; /* Estimate of cost (in bits) of a group. */
3595 size_t i; /* loop counter. */
3596
3597 /* Subtract the Overall Min Value. */
3598 if (OverallMin != 0) {
3599 if (f_secMiss) {
3600 for (i = 0; i < numData; i++) {
3601 if ((Data[i] != li_secMiss) && (Data[i] != li_primMiss)) {
3602 Data[i] -= OverallMin;
3603 /* Check if we accidently adjusted to prim or sec, if so add
3604 * 1. */
3605 if ((Data[i] == li_secMiss) || (Data[i] == li_primMiss)) {
3606 myAssert (1 == 2);
3607 Data[i]++;
3608 if ((Data[i] == li_secMiss) || (Data[i] == li_primMiss)) {
3609 myAssert (1 == 2);
3610 Data[i]++;
3611 }
3612 }
3613 }
3614 }
3615 } else if (f_primMiss) {
3616 for (i = 0; i < numData; i++) {
3617 if (Data[i] != li_primMiss) {
3618 Data[i] -= OverallMin;
3619 /* Check if we accidently adjusted to prim or sec, if so add
3620 * 1. */
3621 if (Data[i] == li_primMiss) {
3622 myAssert (1 == 2);
3623 Data[i]++;
3624 }
3625 }
3626 }
3627 } else {
3628 for (i = 0; i < numData; i++) {
3629 Data[i] -= OverallMin;
3630 }
3631 }
3632 }
3633
3634 myAssert ((f_secMiss == 0) || (f_secMiss == 1));
3635 myAssert ((f_primMiss == 0) || (f_primMiss == 1));
3636
3637 /* Create zero groups. */
3638 *numGroup = 0;
3639 *group = NULL;
3640 if (f_primMiss) {
3641 G.min = Data[0];
3642 G.max = Data[0];
3643 G.num = 1;
3644 G.start = 0;
3645 for (i = 1; i < numData; i++) {
3646 if (G.min == li_primMiss) {
3647 if (Data[i] == li_primMiss) {
3648 G.num++;
3649 if (G.num == (MAX_GROUP_LEN + 1)) {
3650 /* Close a missing group */
3651 G.f_trySplit = 0;
3652 G.f_tryShift = 1;
3653 G.bit = 0;
3654 G.min = 0;
3655 G.max = 0;
3656 G.num = MAX_GROUP_LEN;
3657 (*numGroup)++;
3658 *group = (TDLGroupType *) realloc ((void *) *group,
3659 *numGroup *
3660 sizeof (TDLGroupType));
3661 (*group)[(*numGroup) - 1] = G;
3662 /* Init a missing group. */
3663 G.min = Data[i];
3664 G.max = Data[i];
3665 G.num = 1;
3666 G.start = i;
3667 }
3668 } else {
3669 /* Close a missing group */
3670 G.f_trySplit = 0;
3671 G.f_tryShift = 1;
3672 G.bit = 0;
3673 G.min = 0;
3674 G.max = 0;
3675 (*numGroup)++;
3676 *group = (TDLGroupType *) realloc ((void *) *group,
3677 *numGroup *
3678 sizeof (TDLGroupType));
3679 (*group)[(*numGroup) - 1] = G;
3680 /* Init a non-missing group. */
3681 G.min = Data[i];
3682 G.max = Data[i];
3683 G.num = 1;
3684 G.start = i;
3685 }
3686 } else {
3687 if (Data[i] == li_primMiss) {
3688 /* Close a non-missing group */
3689 G.f_trySplit = 1;
3690 G.f_tryShift = 1;
3691 G.bit = (char) power ((uInt4) (G.max - G.min),
3692 f_secMiss + f_primMiss);
3693 myAssert (G.bit != 0);
3694 if ((G.min == 0) && (G.bit == 0) && (f_primMiss == 1)) {
3695 printf ("Warning: potential confusion between const value "
3696 "and prim-missing.\n");
3697 G.bit = 1;
3698 }
3699 (*numGroup)++;
3700 *group = (TDLGroupType *) realloc ((void *) *group,
3701 *numGroup *
3702 sizeof (TDLGroupType));
3703 (*group)[(*numGroup) - 1] = G;
3704 /* Init a missing group. */
3705 G.min = Data[i];
3706 G.max = Data[i];
3707 G.num = 1;
3708 G.start = i;
3709 } else {
3710 if (G.min > Data[i]) {
3711 G.min = Data[i];
3712 } else if (G.max < Data[i]) {
3713 G.max = Data[i];
3714 }
3715 G.num++;
3716 }
3717 }
3718 }
3719 if (G.min == li_primMiss) {
3720 /* Close a missing group */
3721 G.f_trySplit = 0;
3722 G.f_tryShift = 1;
3723 G.bit = 0;
3724 G.min = 0;
3725 G.max = 0;
3726 (*numGroup)++;
3727 *group = (TDLGroupType *) realloc ((void *) *group,
3728 *numGroup *
3729 sizeof (TDLGroupType));
3730 (*group)[(*numGroup) - 1] = G;
3731 } else {
3732 /* Close a non-missing group */
3733 G.f_trySplit = 1;
3734 G.f_tryShift = 1;
3735 G.bit = (char) power ((uInt4) (G.max - G.min),
3736 f_secMiss + f_primMiss);
3737 myAssert (G.bit != 0);
3738 if ((G.min == 0) && (G.bit == 0) && (f_primMiss == 1)) {
3739 printf ("Warning: potential confusion between const value and "
3740 "prim-missing.\n");
3741 G.bit = 1;
3742 }
3743 (*numGroup)++;
3744 *group = (TDLGroupType *) realloc ((void *) *group,
3745 *numGroup *
3746 sizeof (TDLGroupType));
3747 (*group)[(*numGroup) - 1] = G;
3748 }
3749 } else {
3750 /* Already handled the f_primMiss case */
3751 if (f_secMiss) {
3752 findMaxMin2 (Data, 0, numData, li_primMiss, li_secMiss, &A_min,
3753 &A_max);
3754 } else {
3755 findMaxMin0 (Data, 0, numData, &A_min, &A_max);
3756 }
3757 G.start = 0;
3758 G.num = numData;
3759 G.min = A_min;
3760 G.max = A_max;
3761 G.bit = (char) power ((uInt4) (A_max - A_min), f_secMiss + f_primMiss);
3762 G.f_trySplit = 1;
3763 G.f_tryShift = 1;
3764 *numGroup = 1;
3765 *group = (TDLGroupType *) malloc (sizeof (TDLGroupType));
3766 (*group)[0] = G;
3767 }
3768
3769 lclGroup = NULL;
3770 numLclGroup = 0;
3771 *groupSize = ComputeGroupSize (*group, *numGroup, ibit, jbit, kbit);
3772 xFactor = *ibit + *jbit + *kbit;
3773 #ifdef DEBUG
3774 /*
3775 printf ("NumGroup = %d: Bytes = %ld: XFactor %d\n", *numGroup,
3776 (*groupSize / 8) + 1, xFactor);
3777 */
3778 #endif
3779
3780 f_adjust = 1;
3781 while (f_adjust) {
3782 f_adjust = splitGroup (Data, numData, *group, *numGroup, &lclGroup,
3783 &numLclGroup, f_primMiss, li_primMiss,
3784 f_secMiss, li_secMiss, xFactor);
3785 free (*group);
3786 *group = lclGroup;
3787 *numGroup = numLclGroup;
3788
3789 if (f_adjust) {
3790 shiftGroup (Data, numData, group, numGroup, f_primMiss,
3791 li_primMiss, f_secMiss, li_secMiss, xFactor);
3792 *groupSize = ComputeGroupSize (*group, *numGroup, ibit, jbit, kbit);
3793 if (xFactor != *ibit + *jbit + *kbit) {
3794 for (i = 0; i < *numGroup; i++) {
3795 if (((*group)[i].num > *ibit + *jbit + *kbit) &&
3796 ((*group)[i].num <= xFactor)) {
3797 (*group)[i].f_trySplit = 1;
3798 }
3799 }
3800 }
3801 xFactor = *ibit + *jbit + *kbit;
3802 #ifdef DEBUG
3803 /*
3804 printf ("NumGroup = %d: Bytes = %ld: XFactor %d\n", *numGroup,
3805 (*groupSize / 8) + 1, xFactor);
3806 fflush (stdout);
3807 */
3808 #endif
3809 }
3810 }
3811 }
3812
3813 /*****************************************************************************
3814 * GroupPack() --
3815 *
3816 * Arthur Taylor / MDL
3817 *
3818 * PURPOSE
3819 * To compute groups for packing the data using complex or second order
3820 * complex packing.
3821 *
3822 * ARGUMENTS
3823 * Src = The original data. (Input)
3824 * Dst = The scaled data. (Output)
3825 * numData = The number of elements in data. (Input)
3826 * DSF = Decimal Scale Factor for scaling the data. (Input)
3827 * BSF = Binary Scale Factor for scaling the data. (Input)
3828 * f_primMiss = Flag saying if we have a primary missing value (In/Out)
3829 * primMiss = primary missing value. (In/Out)
3830 * f_secMiss = Flag saying if we have a secondary missing value (In/Out)
3831 * secMiss = secondary missing value. (In/Out)
3832 * f_grid = Flag if this is grid data (or vector) (Input)
3833 * NX = The number of X values. (Input)
3834 * NY = The number of Y values. (Input)
3835 * f_sndOrder = Flag if we should do second order diffencing (Output)
3836 * group = Resulting groups. (Output)
3837 * numGroup = Number of groups. (Output)
3838 * Min = Overall minimum. (Output)
3839 * a1 = if f_sndOrder, the first first order difference (Output)
3840 * b2 = if f_sndOrder, the first second order difference (Output)
3841 * groupSize = How many bytes the groups and data will take. (Output)
3842 * ibit = # of bits for largest minimum value in groups (Output)
3843 * jbit = # of bits for largest # bits in groups (Output)
3844 * kbit = # of bits for largest # values in groups (Output)
3845 *
3846 * FILES/DATABASES: None
3847 *
3848 * RETURNS: int (could use errSprintf())
3849 * 0 = OK
3850 * -1 = Primary or Secondary missing value == 0.
3851 *
3852 * HISTORY
3853 * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
3854 * 1/2005 AAT: Cleaned up.
3855 *
3856 * NOTES
3857 *****************************************************************************
3858 */
GroupPack(double * Src,sInt4 ** Dst,sInt4 numData,int DSF,int BSF,char * f_primMiss,double * primMiss,char * f_secMiss,double * secMiss,char f_grid,short int NX,short int NY,char * f_sndOrder,TDLGroupType ** group,size_t * numGroup,sInt4 * Min,sInt4 * a1,sInt4 * b2,sInt4 * groupSize,size_t * ibit,size_t * jbit,size_t * kbit)3859 static int GroupPack (double *Src, sInt4 **Dst, sInt4 numData,
3860 int DSF, int BSF, char *f_primMiss, double *primMiss,
3861 char *f_secMiss, double *secMiss, char f_grid,
3862 short int NX, short int NY, char *f_sndOrder,
3863 TDLGroupType ** group, size_t *numGroup,
3864 sInt4 *Min, sInt4 *a1, sInt4 *b2,
3865 sInt4 *groupSize, size_t *ibit, size_t *jbit,
3866 size_t *kbit)
3867 {
3868 sInt4 *SecDiff = NULL; /* Consists of the 2nd order differences if *
3869 * requested. */
3870 sInt4 *Data; /* The scaled data. */
3871 char f_min; /* Flag saying overallMin is valid. */
3872 sInt4 overallMin = 0; /* The overall min of the scaled data. */
3873 sInt4 secMin; /* The overall min of the 2nd order differences */
3874 sInt4 li_primMiss = 0; /* The scaled primary missing value */
3875 sInt4 li_secMiss = 0; /* The scaled secondary missing value */
3876 int minGroup = 20; /* The minimum group size. Equivalent to xFactor?
3877 * Chose 20 because that was a good estimate of
3878 * XFactor. */
3879
3880 /* Check consistency of f_primMiss and f_secMiss. */
3881 if (*primMiss == *secMiss) {
3882 *f_secMiss = 0;
3883 }
3884 if ((*f_secMiss) && (!(*f_primMiss))) {
3885 *f_primMiss = *f_secMiss;
3886 *primMiss = *secMiss;
3887 *f_secMiss = 0;
3888 }
3889 if (*f_secMiss && (*secMiss == 0)) {
3890 errSprintf ("Error: Secondary missing value not allowed to = 0.\n");
3891 return -1;
3892 }
3893 if (*f_primMiss && (*primMiss == 0)) {
3894 errSprintf ("Error: Primary missing value not allowed to = 0.\n");
3895 return -1;
3896 }
3897
3898 /* Check minGroup size. */
3899 if (minGroup > numData) {
3900 minGroup = numData;
3901 }
3902
3903 /* Scale the data and check if we can change f_prim or f_sec. */
3904 /* Note: if we use sec_diff, we have a different overall min. */
3905 f_min = 0;
3906 Data = (sInt4 *) malloc (numData * sizeof (sInt4));
3907 TDL_ScaleData (Src, Data, numData, DSF, BSF, f_primMiss, primMiss,
3908 f_secMiss, secMiss, &f_min, &overallMin);
3909 /* Note: ScaleData also scales missing values. */
3910 if (*f_primMiss) {
3911 li_primMiss = (sInt4) (*primMiss * SCALE_MISSING + .5);
3912 }
3913 if (*f_secMiss) {
3914 li_secMiss = (sInt4) (*secMiss * SCALE_MISSING + .5);
3915 }
3916
3917 /* Reason this is after TDL_ScaleData is we don't want to reoder the
3918 * caller's copy of the data. */
3919 if (f_grid) {
3920 TDL_ReorderGrid (Data, NX, NY);
3921 } else {
3922 /* TDLPack requires the following (see pack2d.f and pack1d.f) */
3923 *f_sndOrder = 0;
3924 }
3925 /* TDLPack requires the following (see pack2d.f) */
3926 if (*f_secMiss) {
3927 *f_sndOrder = 0;
3928 }
3929 /* If overallMin is "invalid" then they are all prim or sec values. */
3930 /* See pack2d.f line 336 */
3931 /* IF ALL VALUES ARE MISSING, JUST PACK; DON'T CONSIDER 2ND ORDER
3932 * DIFFERENCES. */
3933 if (!f_min) {
3934 *f_sndOrder = 0;
3935 }
3936
3937 /* This has to be after TDL_ReorderGrid */
3938 if (*f_sndOrder) {
3939 SecDiff = (sInt4 *) malloc (numData * sizeof (sInt4));
3940 if (TDL_GetSecDiff (Data, numData, SecDiff, *f_primMiss, li_primMiss,
3941 a1, b2, &secMin)) {
3942 /* Problem finding SecDiff, so we don't bother with it. */
3943 *f_sndOrder = 0;
3944 } else {
3945 /* Check if it is worth doing second order differences. */
3946 if (*f_primMiss) {
3947 *f_sndOrder = TDL_UseSecDiff_Prim (Data, numData, SecDiff,
3948 li_primMiss, minGroup);
3949 } else {
3950 *f_sndOrder = TDL_UseSecDiff (Data, numData, SecDiff, minGroup);
3951 }
3952 }
3953 }
3954
3955 /* Side affect of GroupIt2: it subtracts OverallMin from Data. */
3956 if (!(*f_sndOrder)) {
3957 GroupIt (overallMin, Data, numData, group, numGroup, *f_primMiss,
3958 li_primMiss, *f_secMiss, li_secMiss, groupSize, ibit, jbit,
3959 kbit);
3960 *Min = overallMin;
3961 *a1 = 0;
3962 *b2 = 0;
3963 *Dst = Data;
3964 free (SecDiff);
3965 } else {
3966 GroupIt (secMin, SecDiff, numData, group, numGroup, *f_primMiss,
3967 li_primMiss, *f_secMiss, li_secMiss, groupSize, ibit, jbit,
3968 kbit);
3969 *Min = secMin;
3970 *Dst = SecDiff;
3971 free (Data);
3972 }
3973 return 0;
3974 }
3975
3976 /*****************************************************************************
3977 * WriteTDLPRecord() --
3978 *
3979 * Arthur Taylor / MDL
3980 *
3981 * PURPOSE
3982 * Writes a TDLP message to file.
3983 *
3984 * ARGUMENTS
3985 * fp = An opened TDLP file already at the correct location. (Input)
3986 * Data = The data to write. (Input)
3987 * DataLen = Length of Data. (Input)
3988 * DSF = Decimal scale factor to apply to the data (Input)
3989 * BSF = Binary scale factor to apply to the data (Input)
3990 * f_primMiss = Flag saying if we have a primary missing value (Input)
3991 * primMiss = primary missing value. (Input)
3992 * f_secMiss = Flag saying if we have a secondary missing value (Input)
3993 * secMiss = secondary missing value. (Input)
3994 * gds = The grid definition section (Input)
3995 * comment = Describes the kind of data (max 32 bytes). (Input)
3996 * refTime = The reference (creation) time of this message. (Input)
3997 * ID1 = TDLPack ID1 (Input)
3998 * ID2 = TDLPack ID2 (Input)
3999 * ID3 = TDLPack ID3 (Input)
4000 * ID4 = TDLPack ID4 (Input)
4001 * projSec = The projection in seconds (Input)
4002 * processNum = The process number that created it (Input)
4003 * seqNum = The sequence number that created it (Input)
4004 *
4005 * FILES/DATABASES:
4006 * An already opened file pointing to the desired TDLP message.
4007 *
4008 * RETURNS: int (could use errSprintf())
4009 * 0 = OK
4010 * -1 = comment is too long.
4011 * -2 = projHr is inconsistent with ID3.
4012 * -3 = Type of map projection that TDLP can't handle.
4013 * -4 = Primary or Secondary missing value == 0.
4014 *
4015 * HISTORY
4016 * 12/2004 Arthur Taylor (MDL): Created
4017 * 1/2005 AAT: Cleaned up.
4018 *
4019 * NOTES
4020 *****************************************************************************
4021 */
WriteTDLPRecord(FILE * fp,double * Data,sInt4 DataLen,int DSF,int BSF,char f_primMiss,double primMiss,char f_secMiss,double secMiss,gdsType * gds,char * comment,double refTime,sInt4 ID1,sInt4 ID2,sInt4 ID3,sInt4 ID4,sInt4 projSec,sInt4 processNum,sInt4 seqNum)4022 int WriteTDLPRecord (FILE * fp, double *Data, sInt4 DataLen, int DSF,
4023 int BSF, char f_primMiss, double primMiss,
4024 char f_secMiss, double secMiss, gdsType *gds,
4025 char *comment, double refTime, sInt4 ID1,
4026 sInt4 ID2, sInt4 ID3, sInt4 ID4,
4027 sInt4 projSec, sInt4 processNum, sInt4 seqNum)
4028 {
4029 sInt4 *Scaled; /* The scaled data. */
4030 TDLGroupType *group; /* The groups used to pack the data. */
4031 size_t numGroup; /* Number of groups. */
4032 char f_grid = 1; /* Flag if this is gridded data. In theory can handle
4033 * vector data, but haven't tested it. */
4034 char f_sndOrder; /* Flag if we should try second order packing. */
4035
4036 /* TODO: Trace overallMin to figure out how it could be used uninitialized */
4037 sInt4 overallMin; /* Overall min value of the scaled data. */
4038
4039 sInt4 a1; /* 2nd order difference: 1st value. */
4040 sInt4 b2; /* 2nd order difference: 1st 1st order difference */
4041 sInt4 li_primMiss; /* Scaled primary missing value. */
4042 sInt4 li_secMiss; /* Scaled secondary missing value. */
4043 int mbit; /* # of bits for b2. */
4044 int nbit; /* # of bits for overallMin. */
4045 size_t ibit; /* # of bits for largest minimum value in groups */
4046 size_t jbit; /* # of bits for largest # bits in groups */
4047 size_t kbit; /* # of bits for largest # values in groups */
4048 int sec1Len; /* Length of section 1. */
4049 size_t pad; /* Number of bytes to pad the message to get to the
4050 * correct byte boundary. */
4051 sInt4 groupSize; /* How many bytes the groups and data will take. */
4052 sInt4 sec4Len; /* Length of section 4. */
4053 sInt4 tdlRecSize; /* Size of the TDLP message. */
4054 sInt4 recSize; /* Actual record size (including FORTRAN bytes). */
4055 int commentLen; /* Length of comment */
4056 short int projHr; /* The hours part of the forecast projection. */
4057 char projMin; /* The minutes part of the forecast projection. */
4058 sInt4 year; /* The reference year. */
4059 int month, day; /* The reference month day. */
4060 int hour, min; /* The reference hour minute. */
4061 double sec; /* The reference second. */
4062 char f_bitmap = 0; /* Bitmap flag: not implemented in specs. */
4063 char f_simple = 0; /* Simple Pack flag: not implemented in specs. */
4064 int gridType; /* Which type of grid. (Polar, Mercator, Lambert). */
4065 int dataCnt; /* Keeps track of which element we are writing. */
4066 sInt4 max0; /* The max value in a group. Represents primary or *
4067 * secondary missing value depending on scheme. */
4068 sInt4 max1; /* The next to max value in a group. Represents *
4069 * secondary missing value. */
4070 size_t i, j; /* loop counters */
4071 sInt4 li_temp; /* Temporary variable (sInt4). */
4072 short int si_temp; /* Temporary variable (short int). */
4073 double d_temp; /* Temporary variable (double). */
4074 char buffer[6]; /* Used to write reserved values */
4075 uChar pbuf; /* A buffer of bits that weren't written to disk */
4076 sChar pbufLoc; /* Where in pbuf to add more bits. */
4077
4078 commentLen = strlen (comment);
4079 if (commentLen > 32) {
4080 errSprintf ("Error: '%s' is > 32 bytes long\n", comment);
4081 return -1;
4082 }
4083 projHr = projSec / 3600;
4084 projMin = (projSec % 3600) / 60;
4085 if (projHr != (ID3 - ((ID3 / 1000) * 1000))) {
4086 errSprintf ("Error: projHr = %d is inconsistent with ID3 = %ld\n",
4087 projHr, ID3);
4088 return -2;
4089 }
4090 if (f_grid) {
4091 switch (gds->projType) {
4092 case GS3_POLAR:
4093 gridType = TDLP_POLAR;
4094 break;
4095 case GS3_LAMBERT:
4096 gridType = TDLP_LAMBERT;
4097 break;
4098 case GS3_MERCATOR:
4099 gridType = TDLP_MERCATOR;
4100 break;
4101 default:
4102 errSprintf ("TDLPack can't handle GRIB projection type %d\n",
4103 gds->projType);
4104 return -3;
4105 }
4106 }
4107
4108 if (GroupPack (Data, &Scaled, DataLen, DSF, BSF, &f_primMiss, &primMiss,
4109 &f_secMiss, &secMiss, f_grid, gds->Nx, gds->Ny,
4110 &f_sndOrder, &group, &numGroup, &overallMin, &a1, &b2,
4111 &groupSize, &ibit, &jbit, &kbit) != 0) {
4112 return -4;
4113 }
4114
4115 /* Make sure missing data is properly scaled. */
4116 if (f_primMiss) {
4117 li_primMiss = (sInt4) (primMiss * SCALE_MISSING + .5);
4118 }
4119 if (f_secMiss) {
4120 li_secMiss = (sInt4) (secMiss * SCALE_MISSING + .5);
4121 }
4122
4123 /* Compute TDL record size. */
4124 /* *INDENT-OFF* */
4125 /* TDL Record size
4126 * 8 (section 0),
4127 * 39 + strlen(comment) (section 1),
4128 * 0 or 28 (depending on if you have a section 2),
4129 * 0 (section 3),
4130 * 16 + 5 + 1 + nbit(min val) + 16 + 5 + 5 + 5 + GroupSize() (section 4)
4131 * 4 (section 5)
4132 * pad (to even 8 bytes...)
4133 * Group size uses:
4134 * ibit * num_groups + jbit * num_groups +
4135 * kbit * num_groups + group.nbit * group.nvalues */
4136 /* *INDENT-ON* */
4137 sec1Len = 39 + commentLen;
4138 if (overallMin < 0) {
4139 nbit = power (-1 * overallMin, 0);
4140 } else {
4141 nbit = power (overallMin, 0);
4142 }
4143 if (!f_sndOrder) {
4144 if (f_secMiss) {
4145 sec4Len = 16 + (sInt4) (ceil ((5 + 1 + nbit + 16 + 5 + 5 + 5
4146 + groupSize) / 8.));
4147 } else if (f_primMiss) {
4148 sec4Len = 12 + (sInt4) (ceil ((5 + 1 + nbit + 16 + 5 + 5 + 5
4149 + groupSize) / 8.));
4150 } else {
4151 sec4Len = 8 + (sInt4) (ceil ((5 + 1 + nbit + 16 + 5 + 5 + 5
4152 + groupSize) / 8.));
4153 }
4154 } else {
4155 if (b2 < 0) {
4156 mbit = power (-1 * b2, 0);
4157 } else {
4158 mbit = power (b2, 0);
4159 }
4160 if (f_secMiss) {
4161 sec4Len =
4162 16 +
4163 (sInt4) (ceil
4164 ((32 + 5 + 1 + mbit + 5 + 1 + nbit + 16 + 5 + 5 + 5 +
4165 groupSize) / 8.));
4166 } else if (f_primMiss) {
4167 sec4Len =
4168 12 +
4169 (sInt4) (ceil
4170 ((32 + 5 + 1 + mbit + 5 + 1 + nbit + 16 + 5 + 5 + 5 +
4171 groupSize) / 8.));
4172 } else {
4173 sec4Len =
4174 8 +
4175 (sInt4) (ceil
4176 ((32 + 5 + 1 + mbit + 5 + 1 + nbit + 16 + 5 + 5 + 5 +
4177 groupSize) / 8.));
4178 }
4179 }
4180 if (f_grid) {
4181 tdlRecSize = 8 + sec1Len + 28 + 0 + sec4Len + 4;
4182 } else {
4183 tdlRecSize = 8 + sec1Len + 0 + 0 + sec4Len + 4;
4184 }
4185 /* Actual recSize is 8 + round(tdlRecSize) to nearest 8 byte boundary. */
4186 recSize = 8 + (sInt4) (ceil (tdlRecSize / 8.0)) * 8;
4187 pad = (int) ((recSize - 8) - tdlRecSize);
4188
4189 /* --- Start Writing record. --- */
4190 /* First write FORTRAN record information */
4191 const size_t read_size = fread(&(recSize), sizeof (sInt4), 1, fp);
4192 if (read_size != 1) {
4193 fprintf(stderr, "WARNING: tdlpack.cpp read of recSize failed.\n");
4194 }
4195 li_temp = 0;
4196 FWRITE_BIG (&(li_temp), sizeof (sInt4), 1, fp);
4197 li_temp = recSize - 8; /* FORTRAN rec length. */
4198 FWRITE_BIG (&(li_temp), sizeof (sInt4), 1, fp);
4199
4200 /* Now write section 0. */
4201 fwrite ("TDLP", sizeof (char), 4, fp);
4202 FWRITE_ODDINT_BIG (&(tdlRecSize), 3, fp);
4203 /* version number */
4204 fputc (0, fp);
4205
4206 /* Now write section 1. */
4207 fputc (sec1Len, fp);
4208 /* output type specification... */
4209 i = 0;
4210 if (f_grid) {
4211 i |= 1;
4212 }
4213 if (f_bitmap) {
4214 i |= 2;
4215 }
4216 fputc (i, fp);
4217 /* tempTime = gmtime (&(refTime));*/
4218 Clock_PrintDate (refTime, &year, &month, &day, &hour, &min, &sec);
4219 /* year = tempTime->tm_year + 1900;
4220 month = tempTime->tm_mon + 1;
4221 day = tempTime->tm_mday;
4222 hour = tempTime->tm_hour;
4223 min = tempTime->tm_min;
4224 */
4225 si_temp = year;
4226 FWRITE_BIG (&si_temp, sizeof (short int), 1, fp);
4227 fputc (month, fp);
4228 fputc (day, fp);
4229 fputc (hour, fp);
4230 fputc (min, fp);
4231 li_temp = (year * 1000000L + month * 10000L + day * 100 + hour);
4232 FWRITE_BIG (&li_temp, sizeof (sInt4), 1, fp);
4233 FWRITE_BIG (&ID1, sizeof (sInt4), 1, fp);
4234 FWRITE_BIG (&ID2, sizeof (sInt4), 1, fp);
4235 FWRITE_BIG (&ID3, sizeof (sInt4), 1, fp);
4236 FWRITE_BIG (&ID4, sizeof (sInt4), 1, fp);
4237 FWRITE_BIG (&projHr, sizeof (short int), 1, fp);
4238 fputc (projMin, fp);
4239 fputc (processNum, fp);
4240 fputc (seqNum, fp);
4241 i = (DSF < 0) ? 128 - DSF : DSF;
4242 fputc (i, fp);
4243 i = (BSF < 0) ? 128 - BSF : BSF;
4244 fputc (i, fp);
4245 /* Reserved: 3 bytes of 0. */
4246 li_temp = 0;
4247 fwrite (&li_temp, sizeof (char), 3, fp);
4248 fputc (commentLen, fp);
4249 fwrite (comment, sizeof (char), commentLen, fp);
4250
4251 /* Now write section 2. */
4252 if (f_grid) {
4253 fputc (28, fp);
4254 fputc (gridType, fp);
4255 si_temp = gds->Nx;
4256 FWRITE_BIG (&si_temp, sizeof (short int), 1, fp);
4257 si_temp = gds->Ny;
4258 FWRITE_BIG (&si_temp, sizeof (short int), 1, fp);
4259 li_temp = (sInt4) (gds->lat1 * 10000. + .5);
4260 if (li_temp < 0) {
4261 pbuf = 128;
4262 li_temp = -1 * li_temp;
4263 } else {
4264 pbuf = 0;
4265 }
4266 pbufLoc = 7;
4267 fileBitWrite (&(li_temp), sizeof (li_temp), 23, fp, &pbuf, &pbufLoc);
4268 myAssert (pbufLoc == 8);
4269 myAssert (pbuf == 0);
4270 d_temp = 360 - gds->lon1;
4271 if (d_temp < 0)
4272 d_temp += 360;
4273 if (d_temp > 360)
4274 d_temp -= 360;
4275 li_temp = (sInt4) (d_temp * 10000. + .5);
4276 if (li_temp < 0) {
4277 pbuf = 128;
4278 li_temp = -1 * li_temp;
4279 }
4280 pbufLoc = 7;
4281 fileBitWrite (&(li_temp), sizeof (li_temp), 23, fp, &pbuf, &pbufLoc);
4282 myAssert (pbufLoc == 8);
4283 myAssert (pbuf == 0);
4284 d_temp = 360 - gds->orientLon;
4285 if (d_temp < 0)
4286 d_temp += 360;
4287 if (d_temp > 360)
4288 d_temp -= 360;
4289 li_temp = (sInt4) (d_temp * 10000. + .5);
4290 if (li_temp < 0) {
4291 pbuf = 128;
4292 li_temp = -1 * li_temp;
4293 }
4294 pbufLoc = 7;
4295 fileBitWrite (&(li_temp), sizeof (li_temp), 23, fp, &pbuf, &pbufLoc);
4296 myAssert (pbufLoc == 8);
4297 myAssert (pbuf == 0);
4298 li_temp = (sInt4) (gds->Dx * 1000. + .5);
4299 FWRITE_BIG (&li_temp, sizeof (sInt4), 1, fp);
4300 li_temp = (sInt4) (gds->meshLat * 10000. + .5);
4301 if (li_temp < 0) {
4302 pbuf = 128;
4303 li_temp = -1 * li_temp;
4304 }
4305 pbufLoc = 7;
4306 fileBitWrite (&(li_temp), sizeof (li_temp), 23, fp, &pbuf, &pbufLoc);
4307 myAssert (pbufLoc == 8);
4308 myAssert (pbuf == 0);
4309 memset (buffer, 0, 6);
4310 fwrite (buffer, sizeof (char), 6, fp);
4311 }
4312
4313 /* Now write section 3. */
4314 /* Bitmap is not supported, skipping. */
4315 myAssert (!f_bitmap);
4316
4317 /* Now write section 4. */
4318 FWRITE_ODDINT_BIG (&(sec4Len), 3, fp);
4319 i = 0;
4320 if (f_secMiss)
4321 i |= 1;
4322 if (f_primMiss)
4323 i |= 2;
4324 if (f_sndOrder)
4325 i |= 4;
4326 if (!f_simple)
4327 i |= 8;
4328 if (!f_grid)
4329 i |= 16;
4330 fputc (i, fp);
4331 li_temp = DataLen;
4332 FWRITE_BIG (&li_temp, sizeof (sInt4), 1, fp);
4333 if (f_primMiss) {
4334 FWRITE_BIG (&(li_primMiss), sizeof (sInt4), 1, fp);
4335 if (f_secMiss) {
4336 FWRITE_BIG (&(li_secMiss), sizeof (sInt4), 1, fp);
4337 }
4338 }
4339 if (f_sndOrder) {
4340 if (a1 < 0) {
4341 pbuf = 128;
4342 li_temp = -1 * a1;
4343 } else {
4344 pbuf = 0;
4345 li_temp = a1;
4346 }
4347 pbufLoc = 7;
4348 fileBitWrite (&(li_temp), sizeof (li_temp), 31, fp, &pbuf, &pbufLoc);
4349 myAssert (pbufLoc == 8);
4350 myAssert (pbuf == 0);
4351 fileBitWrite (&mbit, sizeof (mbit), 5, fp, &pbuf, &pbufLoc);
4352 if (b2 < 0) {
4353 i = 1;
4354 li_temp = -1 * b2;
4355 } else {
4356 i = 0;
4357 li_temp = b2;
4358 }
4359 fileBitWrite (&i, sizeof (i), 1, fp, &pbuf, &pbufLoc);
4360 myAssert (pbufLoc == 2);
4361 fileBitWrite (&li_temp, sizeof (li_temp), (unsigned short int) mbit,
4362 fp, &pbuf, &pbufLoc);
4363 }
4364 fileBitWrite (&nbit, sizeof (nbit), 5, fp, &pbuf, &pbufLoc);
4365 if (overallMin < 0) {
4366 i = 1;
4367 li_temp = -1 * overallMin;
4368 } else {
4369 i = 0;
4370 li_temp = overallMin;
4371 }
4372 fileBitWrite (&i, sizeof (i), 1, fp, &pbuf, &pbufLoc);
4373 fileBitWrite (&li_temp, sizeof (li_temp), (unsigned short int) nbit, fp,
4374 &pbuf, &pbufLoc);
4375 fileBitWrite (&numGroup, sizeof (numGroup), 16, fp, &pbuf, &pbufLoc);
4376 fileBitWrite (&ibit, sizeof (ibit), 5, fp, &pbuf, &pbufLoc);
4377 fileBitWrite (&jbit, sizeof (jbit), 5, fp, &pbuf, &pbufLoc);
4378 fileBitWrite (&kbit, sizeof (kbit), 5, fp, &pbuf, &pbufLoc);
4379 for (i = 0; i < numGroup; i++) {
4380 fileBitWrite (&(group[i].min), sizeof (sInt4),
4381 (unsigned short int) ibit, fp, &pbuf, &pbufLoc);
4382 }
4383 for (i = 0; i < numGroup; i++) {
4384 fileBitWrite (&(group[i].bit), sizeof (char),
4385 (unsigned short int) jbit, fp, &pbuf, &pbufLoc);
4386 }
4387 #ifdef DEBUG
4388 li_temp = 0;
4389 #endif
4390 for (i = 0; i < numGroup; i++) {
4391 fileBitWrite (&(group[i].num), sizeof (sInt4),
4392 (unsigned short int) kbit, fp, &pbuf, &pbufLoc);
4393 #ifdef DEBUG
4394 li_temp += group[i].num;
4395 #endif
4396 }
4397 #ifdef DEBUG
4398 /* Sanity check ! */
4399 if (li_temp != DataLen) {
4400 printf ("Total packed in groups %d != DataLen %d\n", li_temp,
4401 DataLen);
4402 }
4403 myAssert (li_temp == DataLen);
4404 #endif
4405
4406 dataCnt = 0;
4407 /* Start going through the data. Grid data has already been reordered.
4408 * Already taken care of 2nd order differences. Data at this point should
4409 * contain a stream of either 2nd order differences or 0 order
4410 * differences. We have also removed overallMin from all non-missing
4411 * elements in the stream */
4412 if (f_secMiss) {
4413 for (i = 0; i < numGroup; i++) {
4414 max0 = (1 << group[i].bit) - 1;
4415 max1 = (1 << group[i].bit) - 2;
4416 for (j = 0; j < group[i].num; j++) {
4417 if (Scaled[dataCnt] == li_primMiss) {
4418 li_temp = max0;
4419 } else if (Scaled[dataCnt] == li_secMiss) {
4420 li_temp = max1;
4421 } else {
4422 li_temp = Scaled[dataCnt] - group[i].min;
4423 }
4424 fileBitWrite (&(li_temp), sizeof (sInt4), group[i].bit, fp,
4425 &pbuf, &pbufLoc);
4426 dataCnt++;
4427 }
4428 }
4429 } else if (f_primMiss) {
4430 for (i = 0; i < numGroup; i++) {
4431 /* see what happens when bit == 0. */
4432 max0 = (1 << group[i].bit) - 1;
4433 for (j = 0; j < group[i].num; j++) {
4434 if (group[i].bit != 0) {
4435 if (Scaled[dataCnt] == li_primMiss) {
4436 li_temp = max0;
4437 } else {
4438 li_temp = Scaled[dataCnt] - group[i].min;
4439 }
4440 fileBitWrite (&(li_temp), sizeof (sInt4), group[i].bit, fp,
4441 &pbuf, &pbufLoc);
4442 }
4443 dataCnt++;
4444 }
4445 }
4446 } else {
4447 for (i = 0; i < numGroup; i++) {
4448 for (j = 0; j < group[i].num; j++) {
4449 li_temp = Scaled[dataCnt] - group[i].min;
4450 if (group[i].bit != 0) {
4451 fileBitWrite (&(li_temp), sizeof (sInt4), group[i].bit, fp,
4452 &pbuf, &pbufLoc);
4453 }
4454 dataCnt++;
4455 }
4456 }
4457 }
4458 myAssert (dataCnt == DataLen);
4459 /* flush the PutBit buffer... */
4460 if (pbufLoc != 8) {
4461 fputc ((int) pbuf, fp);
4462 }
4463
4464 /* Now write section 5. */
4465 fwrite ("7777", sizeof (char), 4, fp);
4466
4467 /* Deal with padding at end of record... */
4468 for (i = 0; i < pad; i++) {
4469 fputc (0, fp);
4470 }
4471 /* Now write FORTRAN record information */
4472 FWRITE_BIG (&(recSize), sizeof (sInt4), 1, fp);
4473
4474 free (Scaled);
4475 free (group);
4476 return 0;
4477 }
4478