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