1 /*****************************************************************************
2 * engrib.c
3 *
4 * DESCRIPTION
5 * This file contains simple tools to fill out meta data section prior to
6 * calling NCEP GRIB2 encoding routines.
7 *
8 * HISTORY
9 * 4/2006 Arthur Taylor (MDL): Created.
10 *
11 * NOTES
12 *****************************************************************************
13 */
14 #include <stdio.h>
15 #include <stdlib.h>
16 #include <string.h>
17 #include <math.h>
18
19 #include "grib2api.h"
20 #include "engribapi.h"
21 #include "myassert.h"
22 #include "gridtemplates.h"
23 #include "pdstemplates.h"
24 #include "drstemplates.h"
25
26 #ifdef MEMWATCH
27 #include "memwatch.h"
28 #endif
29
30 #define GRIB2MISSING_u1 (uChar) (0xff)
31 #define GRIB2MISSING_s1 (sChar) -1 * (0x7f)
32 #define GRIB2MISSING_u2 (uShort2) (0xffff)
33 #define GRIB2MISSING_s2 (sShort2) -1 * (0x7fff)
34 #define GRIB2MISSING_u4 (uInt4) (0xffffffff)
35 /* following is -1 * 2&31 because of the way signed integers are stored in
36 GRIB2. */
37 #define GRIB2MISSING_s4 (sInt4) -2147483647
38 /*
39 #define GRIB2MISSING_1 (int) (0xff)
40 #define GRIB2MISSING_2 (int) (0xffff)
41 #define GRIB2MISSING_4 (sInt4) (0xffffffff)
42 */
43
44 /*****************************************************************************
45 * NearestInt() -- Arthur Taylor / MDL
46 *
47 * PURPOSE
48 * Find the nearest integer to the given double.
49 *
50 * ARGUMENTS
51 * a = The given float. (Input)
52 *
53 * RETURNS: sInt4 (The nearest integer)
54 *
55 * 4/2006 Arthur Taylor (MDL): Commented (copied from pack.c).
56 *
57 * NOTES:
58 *****************************************************************************
59 */
NearestInt(double a)60 static sInt4 NearestInt(double a)
61 {
62 return (sInt4)floor(a + .5);
63 }
64
65 /*****************************************************************************
66 * AdjustLon() -- Arthur Taylor / MDL
67 *
68 * PURPOSE
69 * Adjust the longitude so that it is in the range of 0..360.
70 *
71 * ARGUMENTS
72 * lon = The given longitude. (Input)
73 *
74 * RETURNS: double (a longitude in the correct range)
75 *
76 * 4/2006 Arthur Taylor (MDL): Created
77 *
78 * NOTES:
79 *****************************************************************************
80 */
AdjustLon(double lon)81 static double AdjustLon (double lon)
82 {
83 while (lon < 0)
84 lon += 360;
85 while (lon > 360)
86 lon -= 360;
87 return lon;
88 }
89
90 /*****************************************************************************
91 * initEnGribMeta() -- Arthur Taylor / MDL
92 *
93 * PURPOSE
94 * Initialize the dynamic memory in the enGribMeta data structure.
95 *
96 * ARGUMENTS
97 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
98 *
99 * RETURNS: void
100 *
101 * 4/2006 Arthur Taylor (MDL): Created.
102 *
103 * NOTES:
104 *****************************************************************************
105 */
initEnGribMeta(enGribMeta * en)106 void initEnGribMeta(enGribMeta * en)
107 {
108 en->sec2 = NULL;
109 en->lenSec2 = 0;
110 en->gdsTmpl = NULL;
111 en->lenGdsTmpl = 0;
112 en->idefList = NULL;
113 en->idefnum = 0;
114 en->pdsTmpl = NULL;
115 en->lenPdsTmpl = 0;
116 en->coordlist = NULL;
117 en->numcoord = 0;
118 en->drsTmpl = NULL;
119 en->lenDrsTmpl = 0;
120 en->fld = NULL;
121 en->ngrdpts = 0;
122 en->bmap = NULL;
123 en->ibmap = GRIB2MISSING_u1;
124 }
125
126 /*****************************************************************************
127 * freeEnGribMeta() -- Arthur Taylor / MDL
128 *
129 * PURPOSE
130 * Free the dynamic memory in the enGribMeta data structure.
131 *
132 * ARGUMENTS
133 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
134 *
135 * RETURNS: void
136 *
137 * 4/2006 Arthur Taylor (MDL): Created.
138 *
139 * NOTES:
140 *****************************************************************************
141 */
freeEnGribMeta(enGribMeta * en)142 void freeEnGribMeta(enGribMeta * en)
143 {
144 if (en->sec2 != NULL) {
145 free(en->sec2);
146 en->sec2 = NULL;
147 }
148 en->lenSec2 = 0;
149 if (en->gdsTmpl != NULL) {
150 free(en->gdsTmpl);
151 en->gdsTmpl = NULL;
152 }
153 en->lenGdsTmpl = 0;
154 if (en->idefList != NULL) {
155 free(en->idefList);
156 en->idefList = NULL;
157 }
158 en->idefnum = 0;
159 if (en->pdsTmpl != NULL) {
160 free(en->pdsTmpl);
161 en->pdsTmpl = NULL;
162 }
163 en->lenPdsTmpl = 0;
164 if (en->coordlist != NULL) {
165 free(en->coordlist);
166 en->coordlist = NULL;
167 }
168 en->numcoord = 0;
169 if (en->drsTmpl != NULL) {
170 free(en->drsTmpl);
171 en->drsTmpl = NULL;
172 }
173 en->lenDrsTmpl = 0;
174 if (en->fld != NULL) {
175 free(en->fld);
176 en->fld = NULL;
177 }
178 en->ngrdpts = 0;
179 if (en->bmap != NULL) {
180 free(en->bmap);
181 en->bmap = NULL;
182 }
183 en->ibmap = GRIB2MISSING_u1;
184 }
185
186 /*****************************************************************************
187 * fillSect0() -- Arthur Taylor / MDL
188 *
189 * PURPOSE
190 * Complete section 0 data.
191 *
192 * ARGUMENTS
193 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
194 * prodType = Discipline-GRIB Master Table [Code:0.0] (Input)
195 *
196 * RETURNS: void
197 *
198 * 4/2006 Arthur Taylor (MDL): Created.
199 *
200 * NOTES:
201 *****************************************************************************
202 */
fillSect0(enGribMeta * en,uChar prodType)203 void fillSect0(enGribMeta * en, uChar prodType)
204 {
205 en->sec0[0] = prodType;
206 en->sec0[1] = 2;
207 }
208
209 /*****************************************************************************
210 * fillSect1() -- Arthur Taylor / MDL
211 *
212 * PURPOSE
213 * Complete section 1 data.
214 *
215 * ARGUMENTS
216 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
217 * center = center of origin of data. (Input)
218 * subCenter = subCenter of origin of data. (Input)
219 * mstrVer = GRIB2 master table Version (currently 3) (Input)
220 * lclVer = GRIB2 local table version (typically 0) (Input)
221 * refCode = Significance of refTime [Code:1.2] (Input)
222 * refYear = The year of the reference time. (Input)
223 * refMonth = The month of the reference time. (Input)
224 * refDay = The day of the reference time. (Input)
225 * refHour = The hour of the reference time. (Input)
226 * refMin = The min of the reference time. (Input)
227 * refSec = The sec of the reference time. (Input)
228 * prodStat = Production Status of data [Code:1.3] (Input)
229 * typeData = Type of Data [Code:1.4] (Input)
230 *
231 * RETURNS: void
232 *
233 * 4/2006 Arthur Taylor (MDL): Created.
234 *
235 * NOTES:
236 *****************************************************************************
237 */
fillSect1(enGribMeta * en,uShort2 center,uShort2 subCenter,uChar mstrVer,uChar lclVer,uChar refCode,sInt4 refYear,int refMonth,int refDay,int refHour,int refMin,int refSec,uChar prodStat,uChar typeData)238 void fillSect1(enGribMeta * en, uShort2 center, uShort2 subCenter,
239 uChar mstrVer, uChar lclVer, uChar refCode, sInt4 refYear,
240 int refMonth, int refDay, int refHour, int refMin, int refSec,
241 uChar prodStat, uChar typeData)
242 {
243 en->sec1[0] = center;
244 en->sec1[1] = subCenter;
245 en->sec1[2] = mstrVer;
246 en->sec1[3] = lclVer;
247 en->sec1[4] = refCode;
248 en->sec1[5] = refYear;
249 en->sec1[6] = refMonth;
250 en->sec1[7] = refDay;
251 en->sec1[8] = refHour;
252 en->sec1[9] = refMin;
253 en->sec1[10] = refSec;
254 en->sec1[11] = prodStat;
255 en->sec1[12] = typeData;
256 }
257
258 /*****************************************************************************
259 * fillSect2() -- Arthur Taylor / MDL
260 *
261 * PURPOSE
262 * Complete section 2 data.
263 *
264 * ARGUMENTS
265 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
266 * sec2 = Character array with info to be added in Sect 2. (Input)
267 * lenSec2 = Num of bytes of sec2 to be added to Sect 2. (Input)
268 *
269 * RETURNS: void
270 *
271 * 4/2006 Arthur Taylor (MDL): Created.
272 *
273 * NOTES:
274 *****************************************************************************
275 */
fillSect2(enGribMeta * en,uChar * sec2,sInt4 lenSec2)276 void fillSect2(enGribMeta * en, uChar *sec2, sInt4 lenSec2)
277 {
278 if (lenSec2 == 0) {
279 if (en->sec2 != NULL) {
280 free(en->sec2);
281 en->sec2 = NULL;
282 }
283 en->lenSec2 = 0;
284 return;
285 }
286 if (lenSec2 > en->lenSec2) {
287 if (en->sec2 != NULL) {
288 free(en->sec2);
289 }
290 en->sec2 = (uChar *)malloc(lenSec2 * sizeof(char));
291 }
292 en->lenSec2 = lenSec2;
293 memcpy(en->sec2, sec2, lenSec2);
294 }
295
296 /*****************************************************************************
297 * getShpEarth() -- Arthur Taylor / MDL
298 *
299 * PURPOSE
300 * Given a major Earth axis and a minor Earth axis, determine how to store
301 * it in GRIB2.
302 *
303 * ARGUMENTS
304 * majEarth = major axis of earth in km (Input)
305 * minEarth = minor axis of earth in km (Input)
306 * shapeEarth = [Code:3.2] shape of the Earth defined by GRIB2 (Output)
307 * factRad = Scale factor of radius of spherical Earth. (Output)
308 * valRad = Value of radius of spherical Earth (Output)
309 * factMaj = Scale factor of major axis of elliptical Earth. (Output)
310 * valMaj = Value of major axis of elliptical Earth (Output)
311 * factMin = Scale factor of minor axis of elliptical Earth. (Output)
312 * valMin = Value of minor axis of elliptical Earth (Output)
313 *
314 * RETURNS: void
315 *
316 * 4/2006 Arthur Taylor (MDL): Created.
317 *
318 * NOTES:
319 *****************************************************************************
320 */
getShpEarth(double majEarth,double minEarth,sInt4 * shapeEarth,sInt4 * factRad,sInt4 * valRad,sInt4 * factMaj,sInt4 * valMaj,sInt4 * factMin,sInt4 * valMin)321 static void getShpEarth(double majEarth, double minEarth, sInt4 *shapeEarth,
322 sInt4 *factRad, sInt4 *valRad, sInt4 *factMaj,
323 sInt4 *valMaj, sInt4 *factMin, sInt4 *valMin)
324 {
325 *factRad = 0;
326 *factMaj = 0;
327 *factMin = 0;
328 *valRad = 0;
329 *valMaj = 0;
330 *valMin = 0;
331 if (majEarth == minEarth) {
332 if (majEarth == 6367.47) {
333 *shapeEarth = 0;
334 *valRad = 6367470;
335 } else if (majEarth == 6371.229) {
336 *shapeEarth = 6;
337 *valRad = 6371229;
338 } else {
339 *shapeEarth = 1;
340 *valRad = NearestInt(majEarth * 1000);
341 }
342 } else {
343 if ((majEarth == 6378.16) && (minEarth == 6356.775)) {
344 *shapeEarth = 2;
345 *valMaj = 6378160;
346 *valMin = 6356775;
347 } else if ((majEarth == 6378.137) && (minEarth == 6356.752314)) {
348 *shapeEarth = 4;
349 *valMaj = 6378137;
350 /* Should this be 3 or -3? */
351 *factMin = 2;
352 /* 6 356 752 314 > Max unsigned Long Int */
353 *valMin = 635675231;
354 } else {
355 *shapeEarth = 7;
356 *valMaj = NearestInt(majEarth * 1000);
357 *valMin = NearestInt(majEarth * 1000);
358 }
359 }
360 }
361
362 /*****************************************************************************
363 * fillSect3() -- Arthur Taylor / MDL
364 *
365 * PURPOSE
366 * Complete section 3 data.
367 *
368 * ARGUMENTS
369 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
370 * tmplNum = [Code Table 3.1] Grid Definition Template Number. (Input)
371 * majEarth = major axis of earth in km (Input)
372 * minEarth = minor axis of earth in km (Input)
373 * Nx = Number of X coordinates (Input)
374 * Ny = Number of Y coordinates (Input)
375 * lat1 = latitude of first grid point (Input)
376 * lon1 = longitude of first grid point (Input)
377 * [tmplNum=0,10]
378 * lat2 = latitude of last grid point (Input)
379 * lon2 = longitude of last grid point (Input)
380 * [tmplNum=*]
381 * Dx = Grid Length in X direction (Input)
382 * Dy = Grid length in Y direction (Input)
383 * resFlag = [Code Table 3.3] Resolution and components flag (Input)
384 * (bit 3 = flag of i direction given)
385 * (bit 4 = flag of j direction given)
386 * (bit 5 = flag of u/v relative to grid) (typically 0)
387 * scanFlag = [Code Table 3.4] Scanning mode flag (Input)
388 * (bit 1 = flag of flow in -i direction)
389 * (bit 2 = flag of flow in +j direction)
390 * (bit 3 = column oriented (varies by column faster than row)
391 * (bit 4 = flag of boustophotonic) (typically 64)
392 * [tmplNum=20,30]
393 * centerFlag = [Code Table 3.5] Center flag (Input)
394 * (bit 1 = flag of south pole on plane)
395 * (bit 2 = flag of bi-polar projection) (typically 0)
396 * [tmplNum=0]
397 * angle = rule 92.1.6 may not hold, in which case, angle != 0, and
398 * unit = angle/subdivision. Typically 0 (Input)
399 * subDivis = 0 or see angle explanation. (Input)
400 * [tmplNum=10,20,30]
401 * meshLat = Latitude where Dx and Dy are specified (Input)
402 * orientLon = Orientation of the grid (Input)
403 * [tmplNum=30]
404 * scaleLat1 = The tangent latitude. If differs from scaleLat2, then they
405 * the two are the latitudes where the scale should be equal.
406 * One can compute a tangent latitude so that the scale at
407 * scaleLat1 is the same as the scale at scaleLat2. (Input)
408 * scaleLat2 = see scaleLat1. (Input)
409 * southLat = latitude of the south pole of the projection (Input)
410 * southLon = longitude of the south pole of the projection (Input)
411 *
412 * RETURNS: int
413 * > 0 (length of section 3).
414 * -1 if numExOctet != 0 (can't handle idefList yet)
415 * -2 not in list of templates supported by NCEP
416 * -3 can't determine the unit. (see angle / subDivis)
417 * -4 haven't finished mapping this projection to the template.
418 *
419 * 4/2006 Arthur Taylor (MDL): Created.
420 *
421 * NOTES:
422 *****************************************************************************
423 */
fillSect3(enGribMeta * en,uShort2 tmplNum,double majEarth,double minEarth,sInt4 Nx,sInt4 Ny,double lat1,double lon1,double lat2,double lon2,double Dx,double Dy,uChar resFlag,uChar scanFlag,uChar centerFlag,sInt4 angle,sInt4 subDivis,double meshLat,double orientLon,double scaleLat1,double scaleLat2,double southLat,double southLon)424 int fillSect3(enGribMeta * en, uShort2 tmplNum, double majEarth,
425 double minEarth, sInt4 Nx, sInt4 Ny, double lat1, double lon1,
426 double lat2, double lon2, double Dx, double Dy, uChar resFlag,
427 uChar scanFlag, uChar centerFlag, sInt4 angle, sInt4 subDivis,
428 double meshLat, double orientLon, double scaleLat1,
429 double scaleLat2, double southLat, double southLon)
430 {
431 const struct gridtemplate *templatesgrid = get_templatesgrid();
432 int i; /* loop counter over number of GDS templates. */
433 double unit; /* Used to convert from stored value to degrees
434 * lat/lon. See GRIB2 Regulation 92.1.6 */
435
436 if (tmplNum == 65535) {
437 /* can't handle lack of a grid definition template */
438 return -1;
439 }
440 /* srcGridDef = [Code:3.0] 0 => Use a grid template * 1 => predetermined
441 * grid (may not have grid template) * 255 => means no grid applies (no
442 * grid def applies) * for 1,255 templateNum = 65535 means no grid
443 * template. */
444 en->gds[0] = 0;
445 en->gds[1] = Nx * Ny;
446 /* numExOctet = Number of octets needed for each additional grid points
447 * definition. Used to define number of points in each row (or column)
448 * for non-regular grids. 0, if using regular grid. */
449 en->gds[2] = 0;
450 /* interpList = [Code Table 3.11] Interpretation of list for optional
451 * points definition. 0 if no appended list. */
452 en->gds[3] = 0;
453 en->gds[4] = tmplNum;
454
455 /* Find NCEP's template match */
456 for (i = 0; i < MAXGRIDTEMP; i++) {
457 if (templatesgrid[i].template_num == tmplNum) {
458 break;
459 }
460 }
461 if (i == MAXGRIDTEMP) {
462 /* not in list of templates supported by NCEP */
463 return -2;
464 }
465 if (templatesgrid[i].needext) {
466 /* can't handle idefList yet. */
467 return -1;
468 }
469
470 if (en->lenGdsTmpl < templatesgrid[i].mapgridlen) {
471 if (en->gdsTmpl != NULL) {
472 free(en->gdsTmpl);
473 }
474 en->gdsTmpl = (sInt4 *)malloc(templatesgrid[i].mapgridlen *
475 sizeof(sInt4));
476 }
477 en->lenGdsTmpl = templatesgrid[i].mapgridlen;
478
479 /* using 1 / 10^-6 to reduce division later */
480 unit = 1e6;
481 /* lat/lon grid */
482 if (tmplNum == 0) {
483 getShpEarth(majEarth, minEarth, &(en->gdsTmpl[0]), &(en->gdsTmpl[1]),
484 &(en->gdsTmpl[2]), &(en->gdsTmpl[3]), &(en->gdsTmpl[4]),
485 &(en->gdsTmpl[5]), &(en->gdsTmpl[6]));
486 en->gdsTmpl[7] = Nx;
487 en->gdsTmpl[8] = Ny;
488 en->gdsTmpl[9] = angle;
489 en->gdsTmpl[10] = subDivis;
490 if (angle != 0) {
491 if (subDivis == 0) {
492 /* can't determine the unit. */
493 return -3;
494 }
495 /* using 1 / (angle / subdivis) to reduce division later */
496 unit = subDivis / (double)angle;
497 }
498 en->gdsTmpl[11] = NearestInt(lat1 * unit);
499 en->gdsTmpl[12] = NearestInt(AdjustLon(lon1) * unit);
500 en->gdsTmpl[13] = resFlag;
501 en->gdsTmpl[14] = NearestInt(lat2 * unit);
502 en->gdsTmpl[15] = NearestInt(AdjustLon(lon2) * unit);
503 en->gdsTmpl[16] = NearestInt(Dx * unit);
504 en->gdsTmpl[17] = NearestInt(Dy * unit);
505 en->gdsTmpl[18] = scanFlag;
506 return 72;
507 /* mercator grid */
508 } else if (tmplNum == 10) {
509 getShpEarth(majEarth, minEarth, &(en->gdsTmpl[0]), &(en->gdsTmpl[1]),
510 &(en->gdsTmpl[2]), &(en->gdsTmpl[3]), &(en->gdsTmpl[4]),
511 &(en->gdsTmpl[5]), &(en->gdsTmpl[6]));
512 en->gdsTmpl[7] = Nx;
513 en->gdsTmpl[8] = Ny;
514 en->gdsTmpl[9] = NearestInt(lat1 * unit);
515 en->gdsTmpl[10] = NearestInt(AdjustLon(lon1) * unit);
516 en->gdsTmpl[11] = resFlag;
517 en->gdsTmpl[12] = NearestInt(meshLat * unit);
518 en->gdsTmpl[13] = NearestInt(lat2 * unit);
519 en->gdsTmpl[14] = NearestInt(AdjustLon(lon2) * unit);
520 en->gdsTmpl[15] = scanFlag;
521 en->gdsTmpl[16] = NearestInt(AdjustLon(orientLon) * unit);
522 en->gdsTmpl[17] = NearestInt(Dx * 1000.);
523 en->gdsTmpl[18] = NearestInt(Dy * 1000.);
524 return 72;
525 /* polar grid */
526 } else if (tmplNum == 20) {
527 getShpEarth(majEarth, minEarth, &(en->gdsTmpl[0]), &(en->gdsTmpl[1]),
528 &(en->gdsTmpl[2]), &(en->gdsTmpl[3]), &(en->gdsTmpl[4]),
529 &(en->gdsTmpl[5]), &(en->gdsTmpl[6]));
530 en->gdsTmpl[7] = Nx;
531 en->gdsTmpl[8] = Ny;
532 en->gdsTmpl[9] = NearestInt(lat1 * unit);
533 en->gdsTmpl[10] = NearestInt(AdjustLon(lon1) * unit);
534 en->gdsTmpl[11] = resFlag;
535 en->gdsTmpl[12] = NearestInt(meshLat * unit);
536 en->gdsTmpl[13] = NearestInt(AdjustLon(orientLon) * unit);
537 en->gdsTmpl[14] = NearestInt(Dx * 1000.);
538 en->gdsTmpl[15] = NearestInt(Dy * 1000.);
539 en->gdsTmpl[16] = centerFlag;
540 en->gdsTmpl[17] = scanFlag;
541 return 65;
542 /* lambert grid */
543 } else if (tmplNum == 30) {
544 getShpEarth(majEarth, minEarth, &(en->gdsTmpl[0]), &(en->gdsTmpl[1]),
545 &(en->gdsTmpl[2]), &(en->gdsTmpl[3]), &(en->gdsTmpl[4]),
546 &(en->gdsTmpl[5]), &(en->gdsTmpl[6]));
547 en->gdsTmpl[7] = Nx;
548 en->gdsTmpl[8] = Ny;
549 en->gdsTmpl[9] = NearestInt(lat1 * unit);
550 en->gdsTmpl[10] = NearestInt(AdjustLon(lon1) * unit);
551 en->gdsTmpl[11] = resFlag;
552 en->gdsTmpl[12] = NearestInt(meshLat * unit);
553 en->gdsTmpl[13] = NearestInt(AdjustLon(orientLon) * unit);
554 en->gdsTmpl[14] = NearestInt(Dx * 1000.);
555 en->gdsTmpl[15] = NearestInt(Dy * 1000.);
556 en->gdsTmpl[16] = centerFlag;
557 en->gdsTmpl[17] = scanFlag;
558 en->gdsTmpl[18] = NearestInt(scaleLat1 * unit);
559 en->gdsTmpl[19] = NearestInt(scaleLat2 * unit);
560 en->gdsTmpl[20] = NearestInt(southLat * unit);
561 en->gdsTmpl[21] = NearestInt(AdjustLon(southLon) * unit);
562 return 81;
563 }
564 /* Haven't finished mapping this projection to the template. */
565 return -4;
566 }
567
568 /*****************************************************************************
569 * getCodedTime() -- Arthur Taylor / MDL
570 *
571 * PURPOSE
572 * Change from seconds to the time units provided in timeCode.
573 *
574 * ARGUMENTS
575 * timeCode = The time units to convert into (see code table 4.4). (Input)
576 * time = The time in seconds to convert. (Input)
577 * ans = The converted answer. (Output)
578 *
579 * RETURNS: int
580 * 0 = OK
581 * -1 = could not determine.
582 *
583 * 4/2006 Arthur Taylor (MDL): Created.
584 *
585 * NOTES
586 *****************************************************************************
587 */
getCodedTime(uChar timeCode,double time,sInt4 * ans)588 static int getCodedTime(uChar timeCode, double time, sInt4 *ans)
589 {
590 /* Following is a lookup table for unit conversion (see code table 4.4). */
591 static const sInt4 unit2sec[] = {
592 60, 3600, 86400L, 0, 0,
593 0, 0, 0, 0, 0,
594 10800, 21600L, 43200L, 1
595 };
596
597 if (timeCode < 14) {
598 if (unit2sec[timeCode] != 0) {
599 *ans = NearestInt(time / unit2sec[timeCode]);
600 return 0;
601 }
602 }
603 *ans = 0;
604 return -1;
605 }
606
607 /*****************************************************************************
608 * fillSect4_0() -- Arthur Taylor / MDL
609 *
610 * PURPOSE
611 * Complete section 4 (using template 0) data.
612 *
613 * ARGUMENTS
614 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
615 * tmplNum = [Code Table 4.0] Product Definition Template Number (Input)
616 * cat = [Code Table 4.1] General category of Meteo Product. (Input)
617 * subCat = [Code Table 4.2] Specific subcategory of Meteo Product. (In)
618 * genProcess = [Code Table 4.3] type of generating process (Analysis,
619 * Forecast, Probability Forecast, etc) (Input)
620 * bgGenID = Background generating process id. (Input)
621 * genID = Analysis/Forecast generating process id. (Input)
622 * f_valCutOff = Flag if we have a valid cutoff time (Input)
623 * cutOff = Cut off time for forecast (Input)
624 * timeCode = [Code Table 4.4] Unit of time to store in (Input)
625 * foreSec = Forecast time in seconds (Input)
626 * surfType1 = [Code Table 4.5] Type of the first surface (Input)
627 * surfScale1 = scale amount for the first surface (Input)
628 * dSurfVal1 = value of the first surface (before scaling) (Input)
629 * surfType2 = [Code Table 4.5] Type of the second surface (Input)
630 * surfScale2 = scale amount for the second surface (Input)
631 * dSurfVal2 = value of the second surface (before scaling) (Input)
632 *
633 * RETURNS: int
634 * > 0 (length of section 4).
635 * -1 This is specifically for template 4.0 (1,2,5,8,8,12)
636 * -2 not in list of templates supported by NCEP
637 * -3 can't handle the timeCode.
638 *
639 * 4/2006 Arthur Taylor (MDL): Created.
640 *
641 * NOTES:
642 *****************************************************************************
643 */
fillSect4_0(enGribMeta * en,uShort2 tmplNum,uChar cat,uChar subCat,uChar genProcess,uChar bgGenID,uChar genID,uChar f_valCutOff,sInt4 cutOff,uChar timeCode,double foreSec,uChar surfType1,sChar surfScale1,double dSurfVal1,uChar surfType2,sChar surfScale2,double dSurfVal2)644 int fillSect4_0(enGribMeta * en, uShort2 tmplNum, uChar cat, uChar subCat,
645 uChar genProcess, uChar bgGenID, uChar genID,
646 uChar f_valCutOff, sInt4 cutOff, uChar timeCode,
647 double foreSec, uChar surfType1, sChar surfScale1,
648 double dSurfVal1, uChar surfType2, sChar surfScale2,
649 double dSurfVal2)
650 {
651 int i; /* loop counter over number of PDS templates. */
652 const struct pdstemplate *templatespds = get_templatespds();
653
654 /* analysis template (0) */
655 /* In addition templates (1, 2, 5, 8, 9, 12) begin with 4.0 info. */
656 if ((tmplNum != 0) && (tmplNum != 1) && (tmplNum != 2) && (tmplNum != 5) &&
657 (tmplNum != 8) && (tmplNum != 9) && (tmplNum != 10) &&
658 (tmplNum != 12)) {
659 /* This is specifically for template 4.0 (1,2,5,8,9,10,12) */
660 return -1;
661 }
662 en->ipdsnum = tmplNum;
663
664 /* Find NCEP's template match */
665 for (i = 0; i < MAXPDSTEMP; i++) {
666 if (templatespds[i].template_num == tmplNum) {
667 break;
668 }
669 }
670 if (i == MAXPDSTEMP) {
671 /* not in list of templates supported by NCEP */
672 return -2;
673 }
674 /* Allocate memory for it. */
675 if (en->lenPdsTmpl < templatespds[i].mappdslen) {
676 if (en->pdsTmpl != NULL) {
677 free(en->pdsTmpl);
678 }
679 en->pdsTmpl = (sInt4 *)malloc(templatespds[i].mappdslen *
680 sizeof(sInt4));
681 }
682 en->lenPdsTmpl = templatespds[i].mappdslen;
683
684 en->pdsTmpl[0] = cat;
685 en->pdsTmpl[1] = subCat;
686 en->pdsTmpl[2] = genProcess;
687 en->pdsTmpl[3] = bgGenID;
688 en->pdsTmpl[4] = genID;
689 if (f_valCutOff) {
690 en->pdsTmpl[5] = cutOff / 3600;
691 en->pdsTmpl[6] = (cutOff % 3600) / 60;
692 } else {
693 en->pdsTmpl[5] = GRIB2MISSING_u2;
694 en->pdsTmpl[6] = GRIB2MISSING_u1;
695 }
696 en->pdsTmpl[7] = timeCode;
697 if (getCodedTime(timeCode, foreSec, &(en->pdsTmpl[8])) != 0) {
698 /* can't handle this time code yet. */
699 return -3;
700 }
701 en->pdsTmpl[9] = surfType1;
702 if (surfType1 == GRIB2MISSING_u1) {
703 en->pdsTmpl[10] = GRIB2MISSING_s1;
704 en->pdsTmpl[11] = GRIB2MISSING_s4;
705 } else {
706 en->pdsTmpl[10] = surfScale1;
707 en->pdsTmpl[11] = NearestInt (dSurfVal1 * pow (10.0, surfScale1));
708 }
709 en->pdsTmpl[12] = surfType2;
710 if (surfType2 == GRIB2MISSING_u1) {
711 en->pdsTmpl[13] = GRIB2MISSING_s1;
712 en->pdsTmpl[14] = GRIB2MISSING_s4;
713 } else {
714 en->pdsTmpl[13] = surfScale2;
715 en->pdsTmpl[14] = NearestInt (dSurfVal2 * pow (10.0, surfScale2));
716 }
717 return 34;
718 }
719
720 /*****************************************************************************
721 * fillSect4_1() -- Arthur Taylor / MDL
722 *
723 * PURPOSE
724 * Complete section 4 (using template 1) data. Call fillSect4_0 first.
725 *
726 * ARGUMENTS
727 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
728 * tmplNum = [Code Table 4.0] Product Definition Template Number (Input)
729 * typeEnsemble = [Code Table 4.6] Type of ensemble (Input)
730 * perturbNum = Perturbation number (Input)
731 * numFcsts = number of forecasts (Input)
732 *
733 * RETURNS: int
734 * > 0 (length of section 4).
735 * -1 if not template 4.1, or fillSect4_0 was not already called.
736 *
737 * 4/2006 Arthur Taylor (MDL): Created.
738 *
739 * NOTES:
740 *****************************************************************************
741 */
fillSect4_1(enGribMeta * en,uShort2 tmplNum,uChar typeEnsemble,uChar perturbNum,uChar numFcsts)742 int fillSect4_1(enGribMeta * en, uShort2 tmplNum, uChar typeEnsemble,
743 uChar perturbNum, uChar numFcsts)
744 {
745 /* Ensemble template (1) */
746 if (tmplNum != 1) {
747 /* This is specifically for template 4.1 */
748 return -1;
749 }
750 if (en->ipdsnum != tmplNum) {
751 /* Didn't call fillSect4_0 first */
752 return -1;
753 }
754 en->pdsTmpl[15] = typeEnsemble;
755 en->pdsTmpl[16] = perturbNum;
756 en->pdsTmpl[17] = numFcsts;
757 return 37;
758 }
759
760 /*****************************************************************************
761 * fillSect4_2() -- Arthur Taylor / MDL
762 *
763 * PURPOSE
764 * Complete section 4 (using template 2) data. Call fillSect4_0 first.
765 *
766 * ARGUMENTS
767 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
768 * tmplNum = [Code Table 4.0] Product Definition Template Number (Input)
769 * numFcsts = number of forecasts (Input)
770 * derivedFcst = [Code Table 4.7] Derived forecast type (Input)
771 *
772 * RETURNS: int
773 * > 0 (length of section 4).
774 * -1 if not template 4.2, or fillSect4_0 was not already called.
775 *
776 * 4/2006 Arthur Taylor (MDL): Created.
777 *
778 * NOTES:
779 *****************************************************************************
780 */
fillSect4_2(enGribMeta * en,uShort2 tmplNum,uChar numFcsts,uChar derivedFcst)781 int fillSect4_2(enGribMeta * en, uShort2 tmplNum, uChar numFcsts,
782 uChar derivedFcst)
783 {
784 /* derived template (2) */
785 if (tmplNum != 2) {
786 /* This is specifically for template 4.2 */
787 return -1;
788 }
789 if (en->ipdsnum != tmplNum) {
790 /* Didn't call fillSect4_0 first */
791 return -1;
792 }
793 en->pdsTmpl[15] = derivedFcst;
794 en->pdsTmpl[16] = numFcsts;
795 return 36;
796 }
797
798 /*****************************************************************************
799 * fillSect4_5() -- Arthur Taylor / MDL
800 *
801 * PURPOSE
802 * Complete section 4 (using template 5) data. Call fillSect4_0 first.
803 *
804 * ARGUMENTS
805 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
806 * tmplNum = [Code Table 4.0] Product Definition Template Number (Input)
807 * numFcsts = number of forecasts (Input)
808 * foreProbNum = Forecast Probability number (Input)
809 * probType = [Code Table 4.9] (Input)
810 * 0 probability of event below lower limit
811 * 1 probability of event above upper limit
812 * 2 probability of event between lower (inclusive) and upper
813 * 3 probability of event above lower limit
814 * 4 probability of event below upper limit
815 * lowScale = scale amount for the lower limit (Input)
816 * dlowVal = value of the lower limit (before scaling) (Input)
817 * upScale = scale amount for the upper limit (Input)
818 * dupVal = value of the upper limit (before scaling) (Input)
819 *
820 * RETURNS: int
821 * > 0 (length of section 4).
822 * -1 if not template 4.5, or fillSect4_0 was not already called.
823 *
824 * 4/2006 Arthur Taylor (MDL): Created.
825 *
826 * NOTES:
827 *****************************************************************************
828 */
fillSect4_5(enGribMeta * en,uShort2 tmplNum,uChar numFcsts,uChar foreProbNum,uChar probType,sChar lowScale,double dlowVal,sChar upScale,double dupVal)829 int fillSect4_5(enGribMeta * en, uShort2 tmplNum, uChar numFcsts,
830 uChar foreProbNum, uChar probType, sChar lowScale,
831 double dlowVal, sChar upScale, double dupVal)
832 {
833 /* Point Probability template */
834 if (tmplNum != 5) {
835 /* This is specifically for template 4.5 */
836 return -1;
837 }
838 if (en->ipdsnum != tmplNum) {
839 /* Didn't call fillSect4_0 first */
840 return -1;
841 }
842 en->pdsTmpl[15] = foreProbNum;
843 en->pdsTmpl[16] = numFcsts;
844 en->pdsTmpl[17] = probType;
845 if (lowScale == GRIB2MISSING_s1) {
846 en->pdsTmpl[18] = GRIB2MISSING_s1;
847 en->pdsTmpl[19] = GRIB2MISSING_s4;
848 } else {
849 en->pdsTmpl[18] = lowScale;
850 en->pdsTmpl[19] = NearestInt (dlowVal * pow (10.0, lowScale));
851 }
852 if (upScale == GRIB2MISSING_s1) {
853 en->pdsTmpl[20] = GRIB2MISSING_s1;
854 en->pdsTmpl[21] = GRIB2MISSING_s4;
855 } else {
856 en->pdsTmpl[20] = upScale;
857 en->pdsTmpl[21] = NearestInt (dupVal * pow (10.0, upScale));
858 }
859 return 47;
860 }
861
862 /*****************************************************************************
863 * fillSect4_8() -- Arthur Taylor / MDL
864 *
865 * PURPOSE
866 * Complete section 4 (using template 8) data. Call fillSect4_0 first.
867 *
868 * ARGUMENTS
869 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
870 * tmplNum = [Code Table 4.0] Product Definition Template Number (Input)
871 * endYear = The year of the end time (valid Time). (Input)
872 * endMonth = The month of the end time (valid Time). (Input)
873 * endDay = The day of the end time (valid Time). (Input)
874 * endHour = The hour of the end time (valid Time). (Input)
875 * endMin = The min of the end time (valid Time). (Input)
876 * endSec = The sec of the end time (valid Time). (Input)
877 * numInterval = num of time range specifications (Has to = 1) (Input)
878 * numMissing = total num of missing values in statistical process (Input)
879 * interval = time range intervals.
880 *
881 * RETURNS: int
882 * > 0 (length of section 4).
883 * -1 if not template 4.8, or fillSect4_0 was not already called.
884 * -4 can only handle 1 and only 1 time interval
885 *
886 * 4/2006 Arthur Taylor (MDL): Created.
887 *
888 * NOTES:
889 *****************************************************************************
890 */
fillSect4_8(enGribMeta * en,uShort2 tmplNum,sInt4 endYear,int endMonth,int endDay,int endHour,int endMin,int endSec,uChar numInterval,sInt4 numMissing,sect4IntervalType * interval)891 int fillSect4_8 (enGribMeta *en, uShort2 tmplNum, sInt4 endYear, int endMonth,
892 int endDay, int endHour, int endMin, int endSec,
893 uChar numInterval, sInt4 numMissing,
894 sect4IntervalType * interval)
895 {
896 int j; /* loop counter over number of intervals. */
897
898 /* statistic template (8) */
899 if (tmplNum != 8) {
900 /* This is specifically for template 4.8 */
901 return -1;
902 }
903 if (en->ipdsnum != tmplNum) {
904 /* Didn't call fillSect4_0 first */
905 return -1;
906 }
907 en->pdsTmpl[15] = endYear;
908 en->pdsTmpl[16] = endMonth;
909 en->pdsTmpl[17] = endDay;
910 en->pdsTmpl[18] = endHour;
911 en->pdsTmpl[19] = endMin;
912 en->pdsTmpl[20] = endSec;
913 en->pdsTmpl[21] = numInterval;
914 if (numInterval != 1) {
915 /* can only handle 1 and only 1 time interval */
916 return -4;
917 }
918 en->pdsTmpl[22] = numMissing;
919 for (j = 0; j < numInterval; j++) {
920 en->pdsTmpl[23] = interval[j].processID;
921 en->pdsTmpl[24] = interval[j].incrType;
922 en->pdsTmpl[25] = interval[j].timeRangeUnit;
923 en->pdsTmpl[26] = interval[j].lenTime;
924 en->pdsTmpl[27] = interval[j].incrUnit;
925 en->pdsTmpl[28] = interval[j].timeIncr;
926 }
927 return 58;
928 }
929
930 /*****************************************************************************
931 * fillSect4_9() -- Arthur Taylor / MDL
932 *
933 * PURPOSE
934 * Complete section 4 (using template 9) data. Call fillSect4_0 first.
935 *
936 * ARGUMENTS
937 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
938 * tmplNum = [Code Table 4.0] Product Definition Template Number (Input)
939 * numFcsts = number of forecasts (Input)
940 * foreProbNum = Forecast Probability number (Input)
941 * probType = [Code Table 4.9] (Input)
942 * 0 probability of event below lower limit
943 * 1 probability of event above upper limit
944 * 2 probability of event between lower (inclusive) and upper
945 * 3 probability of event above lower limit
946 * 4 probability of event below upper limit
947 * lowScale = scale amount for the lower limit (Input)
948 * dlowVal = value of the lower limit (before scaling) (Input)
949 * upScale = scale amount for the upper limit (Input)
950 * dupVal = value of the upper limit (before scaling) (Input)
951 * endYear = The year of the end time (valid Time). (Input)
952 * endMonth = The month of the end time (valid Time). (Input)
953 * endDay = The day of the end time (valid Time). (Input)
954 * endHour = The hour of the end time (valid Time). (Input)
955 * endMin = The min of the end time (valid Time). (Input)
956 * endSec = The sec of the end time (valid Time). (Input)
957 * numInterval = num of time range specifications (Has to = 1) (Input)
958 * numMissing = total num of missing values in statistical process (Input)
959 * interval = time range intervals.
960 *
961 * RETURNS: int
962 * > 0 (length of section 4).
963 * -1 if not template 4.9, or fillSect4_0 was not already called.
964 * -4 can only handle 1 and only 1 time interval
965 *
966 * 4/2006 Arthur Taylor (MDL): Created.
967 *
968 * NOTES:
969 *****************************************************************************
970 */
fillSect4_9(enGribMeta * en,uShort2 tmplNum,uChar numFcsts,uChar foreProbNum,uChar probType,sChar lowScale,double dlowVal,sChar upScale,double dupVal,sInt4 endYear,int endMonth,int endDay,int endHour,int endMin,int endSec,uChar numInterval,sInt4 numMissing,sect4IntervalType * interval)971 int fillSect4_9(enGribMeta * en, uShort2 tmplNum, uChar numFcsts,
972 uChar foreProbNum, uChar probType, sChar lowScale,
973 double dlowVal, sChar upScale, double dupVal, sInt4 endYear,
974 int endMonth, int endDay, int endHour, int endMin,
975 int endSec, uChar numInterval, sInt4 numMissing,
976 sect4IntervalType * interval)
977 {
978 int j; /* loop counter over number of intervals. */
979
980 /* probability time template (9) */
981 if (tmplNum != 9) {
982 /* This is specifically for template 4.9 */
983 return -1;
984 }
985 if (en->ipdsnum != tmplNum) {
986 /* Didn't call fillSect4_0 first */
987 return -1;
988 }
989 en->pdsTmpl[15] = foreProbNum;
990 en->pdsTmpl[16] = numFcsts;
991 en->pdsTmpl[17] = probType;
992 if (lowScale == GRIB2MISSING_s1) {
993 en->pdsTmpl[18] = GRIB2MISSING_s1;
994 en->pdsTmpl[19] = GRIB2MISSING_s4;
995 } else {
996 en->pdsTmpl[18] = lowScale;
997 en->pdsTmpl[19] = NearestInt (dlowVal * pow (10.0, lowScale));
998 }
999 if (upScale == GRIB2MISSING_s1) {
1000 en->pdsTmpl[20] = GRIB2MISSING_s1;
1001 en->pdsTmpl[21] = GRIB2MISSING_s4;
1002 } else {
1003 en->pdsTmpl[20] = upScale;
1004 en->pdsTmpl[21] = NearestInt (dupVal * pow (10.0, upScale));
1005 }
1006 en->pdsTmpl[22] = endYear;
1007 en->pdsTmpl[23] = endMonth;
1008 en->pdsTmpl[24] = endDay;
1009 en->pdsTmpl[25] = endHour;
1010 en->pdsTmpl[26] = endMin;
1011 en->pdsTmpl[27] = endSec;
1012 en->pdsTmpl[28] = numInterval;
1013 if (numInterval != 1) {
1014 /* can only handle 1 and only 1 time interval */
1015 return -4;
1016 }
1017 en->pdsTmpl[29] = numMissing;
1018 for (j = 0; j < numInterval; j++) {
1019 en->pdsTmpl[30] = interval[j].processID;
1020 en->pdsTmpl[31] = interval[j].incrType;
1021 en->pdsTmpl[32] = interval[j].timeRangeUnit;
1022 en->pdsTmpl[33] = interval[j].lenTime;
1023 en->pdsTmpl[34] = interval[j].incrUnit;
1024 en->pdsTmpl[35] = interval[j].timeIncr;
1025 }
1026 return 71;
1027 }
1028
1029 /*****************************************************************************
1030 * fillSect4_10() -- Arthur Taylor / MDL
1031 *
1032 * PURPOSE
1033 * Complete section 4 (using template 10) data. Call fillSect4_0 first.
1034 *
1035 * ARGUMENTS
1036 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
1037 * tmplNum = [Code Table 4.0] Product Definition Template Number (Input)
1038 * percentile = Percentile value. (Input)
1039 * endYear = The year of the end time (valid Time). (Input)
1040 * endMonth = The month of the end time (valid Time). (Input)
1041 * endDay = The day of the end time (valid Time). (Input)
1042 * endHour = The hour of the end time (valid Time). (Input)
1043 * endMin = The min of the end time (valid Time). (Input)
1044 * endSec = The sec of the end time (valid Time). (Input)
1045 * numInterval = num of time range specifications (Has to = 1) (Input)
1046 * numMissing = total num of missing values in statistical process (Input)
1047 * interval = time range intervals.
1048 *
1049 * RETURNS: int
1050 * > 0 (length of section 4).
1051 * -1 if not template 4.9, or fillSect4_0 was not already called.
1052 * -4 can only handle 1 and only 1 time interval
1053 *
1054 * 5/2006 Arthur Taylor (MDL): Created.
1055 *
1056 * NOTES:
1057 *****************************************************************************
1058 */
fillSect4_10(enGribMeta * en,uShort2 tmplNum,int percentile,sInt4 endYear,int endMonth,int endDay,int endHour,int endMin,int endSec,uChar numInterval,sInt4 numMissing,sect4IntervalType * interval)1059 int fillSect4_10(enGribMeta * en, uShort2 tmplNum, int percentile,
1060 sInt4 endYear, int endMonth, int endDay, int endHour,
1061 int endMin, int endSec, uChar numInterval, sInt4 numMissing,
1062 sect4IntervalType * interval)
1063 {
1064 int j; /* loop counter over number of intervals. */
1065
1066 /* percentile template (10) */
1067 if (tmplNum != 10) {
1068 /* This is specifically for template 4.10 */
1069 return -1;
1070 }
1071 if (en->ipdsnum != tmplNum) {
1072 /* Didn't call fillSect4_0 first */
1073 return -1;
1074 }
1075 en->pdsTmpl[15] = percentile;
1076 en->pdsTmpl[16] = endYear;
1077 en->pdsTmpl[17] = endMonth;
1078 en->pdsTmpl[18] = endDay;
1079 en->pdsTmpl[19] = endHour;
1080 en->pdsTmpl[20] = endMin;
1081 en->pdsTmpl[21] = endSec;
1082 en->pdsTmpl[22] = numInterval;
1083 if (numInterval != 1) {
1084 /* can only handle 1 and only 1 time interval */
1085 return -4;
1086 }
1087 en->pdsTmpl[23] = numMissing;
1088 for (j = 0; j < numInterval; j++) {
1089 en->pdsTmpl[24] = interval[j].processID;
1090 en->pdsTmpl[25] = interval[j].incrType;
1091 en->pdsTmpl[26] = interval[j].timeRangeUnit;
1092 en->pdsTmpl[27] = interval[j].lenTime;
1093 en->pdsTmpl[28] = interval[j].incrUnit;
1094 en->pdsTmpl[29] = interval[j].timeIncr;
1095 }
1096 return 59;
1097 }
1098
1099 /*****************************************************************************
1100 * fillSect4_12() -- Arthur Taylor / MDL
1101 *
1102 * PURPOSE
1103 * Complete section 4 (using template 12) data. Call fillSect4_0 first.
1104 *
1105 * ARGUMENTS
1106 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
1107 * tmplNum = [Code Table 4.0] Product Definition Template Number (Input)
1108 * numFcsts = number of forecasts (Input)
1109 * derivedFcst = [Code Table 4.7] Derived forecast type (Input)
1110 * endYear = The year of the end time (valid Time). (Input)
1111 * endMonth = The month of the end time (valid Time). (Input)
1112 * endDay = The day of the end time (valid Time). (Input)
1113 * endHour = The hour of the end time (valid Time). (Input)
1114 * endMin = The min of the end time (valid Time). (Input)
1115 * endSec = The sec of the end time (valid Time). (Input)
1116 * numInterval = num of time range specifications (Has to = 1) (Input)
1117 * numMissing = total num of missing values in statistical process (Input)
1118 * interval = time range intervals.
1119 *
1120 * RETURNS: int
1121 * > 0 (length of section 4).
1122 * -1 if not template 4.12, or fillSect4_0 was not already called.
1123 * -4 can only handle 1 and only 1 time interval
1124 *
1125 * 4/2006 Arthur Taylor (MDL): Created.
1126 *
1127 * NOTES:
1128 *****************************************************************************
1129 */
fillSect4_12(enGribMeta * en,uShort2 tmplNum,uChar numFcsts,uChar derivedFcst,sInt4 endYear,int endMonth,int endDay,int endHour,int endMin,int endSec,uChar numInterval,sInt4 numMissing,sect4IntervalType * interval)1130 int fillSect4_12(enGribMeta * en, uShort2 tmplNum, uChar numFcsts,
1131 uChar derivedFcst, sInt4 endYear, int endMonth, int endDay,
1132 int endHour, int endMin, int endSec, uChar numInterval,
1133 sInt4 numMissing, sect4IntervalType * interval)
1134 {
1135 int j; /* loop counter over number of intervals. */
1136
1137 /* derived interval template (12) */
1138 if (tmplNum != 12) {
1139 /* This is specifically for template 4.12 */
1140 return -1;
1141 }
1142 if (en->ipdsnum != tmplNum) {
1143 /* Didn't call fillSect4_0 first */
1144 return -1;
1145 }
1146 en->pdsTmpl[15] = derivedFcst;
1147 en->pdsTmpl[16] = numFcsts;
1148 en->pdsTmpl[17] = endYear;
1149 en->pdsTmpl[18] = endMonth;
1150 en->pdsTmpl[19] = endDay;
1151 en->pdsTmpl[20] = endHour;
1152 en->pdsTmpl[21] = endMin;
1153 en->pdsTmpl[22] = endSec;
1154 en->pdsTmpl[23] = numInterval;
1155 if (numInterval != 1) {
1156 /* can only handle 1 and only 1 time interval */
1157 return -4;
1158 }
1159 en->pdsTmpl[24] = numMissing;
1160 for (j = 0; j < numInterval; j++) {
1161 en->pdsTmpl[25] = interval[j].processID;
1162 en->pdsTmpl[26] = interval[j].incrType;
1163 en->pdsTmpl[27] = interval[j].timeRangeUnit;
1164 en->pdsTmpl[28] = interval[j].lenTime;
1165 en->pdsTmpl[29] = interval[j].incrUnit;
1166 en->pdsTmpl[30] = interval[j].timeIncr;
1167 }
1168 return 60;
1169 }
1170
1171 /*****************************************************************************
1172 * fillSect5() -- Arthur Taylor / MDL
1173 *
1174 * PURPOSE
1175 * Complete section 5 data.
1176 *
1177 * ARGUMENTS
1178 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
1179 * tmplNum = [Code Table 5.0] Product Definition Template Number (Input)
1180 * BSF = Binary scale factor (Input)
1181 * DSF = Decimal scale factor (Input)
1182 * [tmplNum=0,2,3,40,41]
1183 * fieldType = [Code Table 5.1] type of original field values. (Input)
1184 * [tmplNum=2,3]
1185 * f_miss = [Code Table 5.5] missing value management used. (Input)
1186 * missPri = Primary missing value (Input)
1187 * missSec = Secondary missing value (Input)
1188 * [tmplNum=3]
1189 * orderOfDiff = Order of differencing (1 or 2) (Input)
1190 *
1191 * RETURNS: int
1192 * > 0 (length of section 5).
1193 * -1 can't handle extended lists yet
1194 * -2 not in list of templates supported by NCEP
1195 * -3 can't handle this order of differencing.
1196 * -4 haven't finished mapping this projection to the template.
1197 *
1198 * 4/2006 Arthur Taylor (MDL): Created.
1199 *
1200 * NOTES:
1201 *****************************************************************************
1202 */
fillSect5(enGribMeta * en,uShort2 tmplNum,sShort2 BSF,sShort2 DSF,uChar fieldType,uChar f_miss,float missPri,float missSec,uChar orderOfDiff)1203 int fillSect5(enGribMeta * en, uShort2 tmplNum, sShort2 BSF, sShort2 DSF,
1204 uChar fieldType, uChar f_miss, float missPri, float missSec,
1205 uChar orderOfDiff)
1206 {
1207 int i; /* loop counter over number of DRS templates. */
1208 const struct drstemplate *templatesdrs = get_templatesdrs();
1209
1210 /* Find NCEP's template match */
1211 for (i = 0; i < MAXDRSTEMP; i++) {
1212 if (templatesdrs[i].template_num == tmplNum) {
1213 break;
1214 }
1215 }
1216 if (i == MAXDRSTEMP) {
1217 /* not in list of templates supported by NCEP */
1218 return -2;
1219 }
1220 if (templatesdrs[i].needext) {
1221 /* can't handle extended data yet. */
1222 return -1;
1223 }
1224
1225 if (en->lenDrsTmpl < templatesdrs[i].mapdrslen) {
1226 if (en->drsTmpl != NULL) {
1227 free(en->drsTmpl);
1228 }
1229 en->drsTmpl = (sInt4 *)malloc(templatesdrs[i].mapdrslen *
1230 sizeof(sInt4));
1231 }
1232 en->lenDrsTmpl = templatesdrs[i].mapdrslen;
1233
1234 en->idrsnum = tmplNum;
1235 /* simple packing */
1236 if (tmplNum == 0) {
1237 en->drsTmpl[0] = 9999; /* missing for Ref value (set later) */
1238 en->drsTmpl[1] = BSF;
1239 en->drsTmpl[2] = DSF;
1240 en->drsTmpl[3] = 9999; /* missing for numBits used (set later) */
1241 en->drsTmpl[4] = fieldType; /* code table 5.1 */
1242 return 21;
1243 /* complex packing */
1244 } else if (tmplNum == 2) {
1245 en->drsTmpl[0] = 9999; /* missing for Ref value (set later) */
1246 en->drsTmpl[1] = BSF;
1247 en->drsTmpl[2] = DSF;
1248 en->drsTmpl[3] = 9999; /* missing for numBits used (set later) */
1249 en->drsTmpl[4] = fieldType; /* code table 5.1 */
1250 en->drsTmpl[5] = 9999; /* missing for group splitting method used */
1251 en->drsTmpl[6] = f_miss;
1252 if (fieldType == 1) {
1253 en->drsTmpl[7] = (sInt4)missPri;
1254 en->drsTmpl[8] = (sInt4)missSec;
1255 } else {
1256 memcpy(&(en->drsTmpl[7]), &missPri, sizeof(float));
1257 memcpy(&(en->drsTmpl[8]), &missSec, sizeof(float));
1258 }
1259 en->drsTmpl[9] = 9999; /* number of groups */
1260 en->drsTmpl[10] = 9999; /* group widths */
1261 en->drsTmpl[11] = 9999; /* numBits for group widths */
1262 en->drsTmpl[12] = 9999; /* ref for group len */
1263 en->drsTmpl[13] = 9999; /* len increment for group lengths */
1264 en->drsTmpl[14] = 9999; /* true len of last group */
1265 en->drsTmpl[15] = 9999; /* numBits used for scaled group lens */
1266 return 47;
1267 /* complex spatial packing */
1268 } else if (tmplNum == 3) {
1269 en->drsTmpl[0] = 9999; /* missing for Ref value (set later) */
1270 en->drsTmpl[1] = BSF;
1271 en->drsTmpl[2] = DSF;
1272 en->drsTmpl[3] = 9999; /* missing for numBits used (set later) */
1273 en->drsTmpl[4] = fieldType; /* code table 5.1 */
1274 en->drsTmpl[5] = 9999; /* missing for group splitting method used */
1275 en->drsTmpl[6] = f_miss;
1276 if (fieldType == 1) {
1277 en->drsTmpl[7] = (sInt4)missPri;
1278 en->drsTmpl[8] = (sInt4)missSec;
1279 } else {
1280 memcpy(&(en->drsTmpl[7]), &missPri, sizeof(float));
1281 memcpy(&(en->drsTmpl[8]), &missSec, sizeof(float));
1282 }
1283 en->drsTmpl[9] = 9999; /* number of groups */
1284 en->drsTmpl[10] = 9999; /* group widths */
1285 en->drsTmpl[11] = 9999; /* numBits for group widths */
1286 en->drsTmpl[12] = 9999; /* ref for group len */
1287 en->drsTmpl[13] = 9999; /* len increment for group lengths */
1288 en->drsTmpl[14] = 9999; /* true len of last group */
1289 en->drsTmpl[15] = 9999; /* numBits used for scaled group lens */
1290 if (orderOfDiff > 2) {
1291 /* NCEP can not handle order of differencing > 2 */
1292 return -3;
1293 }
1294 en->drsTmpl[16] = orderOfDiff;
1295 en->drsTmpl[17] = 9999; /* num extra octets need for spatial differ */
1296 return 49;
1297 /* jpeg2000 packing */
1298 } else if ((tmplNum == 40) || (tmplNum == 40000)) {
1299 en->drsTmpl[0] = 9999; /* missing for Ref value (set later) */
1300 en->drsTmpl[1] = BSF;
1301 en->drsTmpl[2] = DSF;
1302 en->drsTmpl[3] = 9999; /* depth of grayscale image (set later) */
1303 en->drsTmpl[4] = fieldType; /* code table 5.1 */
1304 en->drsTmpl[5] = 9999; /* type of compression used (0 is lossless)
1305 * (code table 5.40) */
1306 en->drsTmpl[6] = 9999; /* compression ratio */
1307 return 23;
1308 /* png packing */
1309 } else if ((tmplNum == 41) || (tmplNum == 40010)) {
1310 en->drsTmpl[0] = 9999; /* missing for Ref value (set later) */
1311 en->drsTmpl[1] = BSF;
1312 en->drsTmpl[2] = DSF;
1313 en->drsTmpl[3] = 9999; /* depth of grayscale image (set later) */
1314 en->drsTmpl[4] = fieldType; /* code table 5.1 */
1315 return 21;
1316 /* spectral packing */
1317 } else if (tmplNum == 50) {
1318 en->drsTmpl[0] = 9999; /* missing for Ref value (set later) */
1319 en->drsTmpl[1] = BSF;
1320 en->drsTmpl[2] = DSF;
1321 en->drsTmpl[3] = 9999; /* num bits used for each packed value */
1322 en->drsTmpl[4] = 9999; /* real part of (0,0) coefficient */
1323 return 24;
1324 /* harmonic packing */
1325 } else if (tmplNum == 51) {
1326 en->drsTmpl[0] = 9999; /* missing for Ref value (set later) */
1327 en->drsTmpl[1] = BSF;
1328 en->drsTmpl[2] = DSF;
1329 en->drsTmpl[3] = 9999; /* num bits used for each packed value */
1330 en->drsTmpl[4] = 9999; /* P - Laplacian scaling factor */
1331 en->drsTmpl[5] = 9999; /* Js - pentagonal resolution parameter */
1332 en->drsTmpl[6] = 9999; /* Ks - pentagonal resolution parameter */
1333 en->drsTmpl[7] = 9999; /* Ms - pentagonal resolution parameter */
1334 en->drsTmpl[8] = 9999; /* Ts - total num values in subset */
1335 en->drsTmpl[9] = 9999; /* Precision of unpacked subset */
1336 return 35;
1337 }
1338 /* Haven't finished mapping this drs to a template. */
1339 return -4;
1340 }
1341
1342 /*****************************************************************************
1343 * fillGrid() -- Arthur Taylor / MDL
1344 *
1345 * PURPOSE
1346 * Completes the data portion. If f_boustify, then it walks through the
1347 * data winding back and forth. Note it does this in a row oriented fashion
1348 * If you need a column oriented fashion because your grid is defined the
1349 * other way, then swap your Nx and Ny in your call.
1350 *
1351 * ARGUMENTS
1352 * en = A pointer to meta data to pass to GRIB2 encoder. (Output)
1353 * data = Data array to add. (Input)
1354 * lenData = Length of Data array. (Input)
1355 * Nx = Number of X coordinates (Input)
1356 * Ny = Number of Y coordinates (Input)
1357 * ibmap = [Code 6.0] Bitmap indicator (Input)
1358 * 0 = bitmap applies and is included in Section 6.
1359 * 1-253 = Predefined bitmap applies
1360 * 254 = Previously defined bitmap applies to this field
1361 * 255 = Bit map does not apply to this product.
1362 * f_boustify = true if we should re-Wrap the grid. (Input)
1363 * f_miss = 1 if missPri valid, 2 if missSec valid. (Input)
1364 * missPri = Primary missing value (Input)
1365 * missSec = Secondary missing value (Input)
1366 *
1367 * RETURNS: int
1368 * > 0 (max length of sect 6 and sect 7).
1369 * -1 Can't handle this kind of bitmap (pre-defined).
1370 * -2 No missing value when trying to create the bmap.
1371 * -3 Can't handle Nx * Ny != lenData.
1372 *
1373 * 4/2006 Arthur Taylor (MDL): Created.
1374 *
1375 * NOTES:
1376 *****************************************************************************
1377 */
fillGrid(enGribMeta * en,double * data,sInt4 lenData,sInt4 Nx,sInt4 Ny,sInt4 ibmap,sChar f_boustify,uChar f_miss,float missPri,float missSec)1378 int fillGrid(enGribMeta * en, double *data, sInt4 lenData, sInt4 Nx, sInt4 Ny,
1379 sInt4 ibmap, sChar f_boustify, uChar f_miss, float missPri,
1380 float missSec)
1381 {
1382 uChar f_flip; /* Used to help keep track of the direction when
1383 * "boustifying" the data. */
1384 sInt4 x; /* loop counter over Nx. */
1385 sInt4 y; /* loop counter over Ny. */
1386 sInt4 ind1; /* index to copy to. */
1387 sInt4 ind2; /* index to copy from. */
1388
1389 if ((ibmap != 0) && (ibmap != 255)) {
1390 /* Can't handle this kind of bitmap (pre-defined). */
1391 return -1;
1392 }
1393 if ((ibmap == 0) && (f_miss != 1) && (f_miss != 2)) {
1394 /* No missing value when trying to create the bmap. */
1395 return -2;
1396 }
1397 if (Nx * Ny != lenData) {
1398 /* Can't handle Nx * Ny != lenData. */
1399 return -3;
1400 }
1401
1402 if (en->ngrdpts < lenData) {
1403 if (en->fld != NULL) {
1404 free(en->fld);
1405 }
1406 en->fld = (float *)malloc(lenData * sizeof(float));
1407 if (ibmap == 0) {
1408 if (en->bmap != NULL) {
1409 free(en->bmap);
1410 }
1411 en->bmap = (sInt4 *)malloc(lenData * sizeof(sInt4));
1412 }
1413 }
1414 en->ngrdpts = lenData;
1415 en->ibmap = ibmap;
1416
1417 /* Now need to walk over data and boustify it and create bmap. */
1418
1419 if (ibmap == 0) {
1420 /* boustify uses row oriented boustification, however for column
1421 * oriented, swap the Ny and Nx in the call to the procedure. */
1422 if (f_boustify) {
1423 f_flip = 0;
1424 for (y = 0; y < Ny; y++) {
1425 for (x = 0; x < Nx; x++) {
1426 ind1 = x + y * Nx;
1427 if (!f_flip) {
1428 ind2 = ind1;
1429 } else {
1430 ind2 = (Nx - x - 1) + y * Nx;
1431 }
1432 en->fld[ind1] = (float)data[ind2];
1433 if ((data[ind2] == missPri) ||
1434 ((f_miss == 2) && (data[ind2] == missSec))) {
1435 en->bmap[ind1] = 0;
1436 } else {
1437 en->bmap[ind1] = 1;
1438 }
1439 }
1440 f_flip = (!f_flip);
1441 }
1442 } else {
1443 for (ind1 = 0; ind1 < lenData; ind1++) {
1444 en->fld[ind1] = (float)data[ind1];
1445 if ((data[ind1] == missPri) ||
1446 ((f_miss == 2) && (data[ind1] == missSec))) {
1447 en->bmap[ind1] = 0;
1448 } else {
1449 en->bmap[ind1] = 1;
1450 }
1451 }
1452 }
1453 /* len(sect6) < 6 + (lenData/8 + 1), len(sect7) < 5 + lenData * 4 */
1454 return (6 + lenData / 8 + 1) + (5 + lenData * 4);
1455 } else {
1456 /* boustify uses row oriented boustification, however for column
1457 * oriented, swap the Ny and Nx in the call to the procedure. */
1458 if (f_boustify) {
1459 f_flip = 0;
1460 for (y = 0; y < Ny; y++) {
1461 for (x = 0; x < Nx; x++) {
1462 ind1 = x + y * Nx;
1463 if (!f_flip) {
1464 ind2 = ind1;
1465 } else {
1466 ind2 = (Nx - x - 1) + y * Nx;
1467 }
1468 en->fld[ind1] = (float)data[ind2];
1469 }
1470 f_flip = (!f_flip);
1471 }
1472 } else {
1473 for (ind1 = 0; ind1 < lenData; ind1++) {
1474 en->fld[ind1] = (float)data[ind1];
1475 }
1476 }
1477 /* len(sect6) = 6, len(sect7) < 5 + lenData * 4 */
1478 return 6 + (5 + lenData * 4);
1479 }
1480 }
1481
1482