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