Line data Source code
1 : /*****************************************************************************
2 : * grib2api.c
3 : *
4 : * DESCRIPTION
5 : * This file contains the API to the GRIB2 libraries which is as close as
6 : * possible to the "official" NWS GRIB2 Library's API. The reason for this
7 : * is so we can use NCEP's unofficial GRIB2 library while having minimal
8 : * impact on the already written drivers.
9 : *
10 : * HISTORY
11 : * 12/2003 Arthur Taylor (MDL / RSIS): Created.
12 : *
13 : * NOTES
14 : *****************************************************************************
15 : */
16 : #include <limits.h>
17 : #include <math.h>
18 : #include <stdlib.h>
19 : #include <string.h>
20 : #include "grib2api.h"
21 : #include "grib2.h"
22 : #include "scan.h"
23 : #include "tendian.h"
24 : #include "myassert.h"
25 : #include "gridtemplates.h"
26 : #include "pdstemplates.h"
27 : #include "drstemplates.h"
28 :
29 : #include "cpl_port.h"
30 :
31 : /* Commented out by GDAL: we can actually include gribtemplates.h */
32 : #if 0
33 :
34 : /* Would prefer to include "gribtemplates.h" but can't since NCEP decided to
35 : * put a const struct in it. The result is if I include it, I get 2 copies
36 : * of the data.
37 : */
38 : extern sInt4 getgridindex(sInt4 number);
39 :
40 : /* have now put this in grib2api.h... may bring it back. */
41 : #define MAXGRIDTEMP 23 /* maximum number of templates */
42 : #define MAXGRIDMAPLEN 200 /* maximum template map length */
43 : struct gridtemplate {
44 : sInt4 template_num;
45 : sInt4 mapgridlen;
46 : sInt4 needext;
47 : sInt4 mapgrid[MAXGRIDMAPLEN];
48 : };
49 : extern const struct gridtemplate templatesgrid[MAXGRIDTEMP];
50 :
51 : /* Would prefer to include "pdstemplates.h" but can't since NCEP decided to
52 : * put a const struct in it. The result is if I include it, I get 2 copies
53 : * of the data.
54 : */
55 : extern sInt4 getpdsindex(sInt4 number);
56 :
57 : /* have now put this in grib2api.h... may bring it back. */
58 : #define MAXPDSTEMP 23 /* maximum number of templates */
59 : #define MAXPDSMAPLEN 200 /* maximum template map length */
60 : struct pdstemplate {
61 : sInt4 template_num;
62 : sInt4 mappdslen;
63 : sInt4 needext;
64 : sInt4 mappds[MAXPDSMAPLEN];
65 : };
66 : extern const struct pdstemplate templatespds[MAXPDSTEMP];
67 :
68 : /* Would prefer to include "drstemplates.h" but can't since NCEP decided to
69 : * put a const struct in it. The result is if I include it, I get 2 copies
70 : * of the data.
71 : */
72 : extern sInt4 getdrsindex(sInt4 number);
73 :
74 : #define MAXDRSTEMP 8 /* maximum number of templates */
75 : #define MAXDRSMAPLEN 200 /* maximum template map length */
76 : struct drstemplate {
77 : sInt4 template_num;
78 : sInt4 mapdrslen;
79 : sInt4 needext;
80 : sInt4 mapdrs[MAXDRSMAPLEN];
81 : };
82 : extern const struct drstemplate templatesdrs[MAXDRSTEMP];
83 : #endif
84 :
85 546952 : static sInt4 FloatToSInt4Clamp(float val) {
86 546952 : if ((double)val >= (double)INT_MAX) return INT_MAX;
87 546952 : if ((double)val <= (double)INT_MIN) return INT_MIN;
88 546943 : if (CPLIsNan(val)) return 0;
89 546914 : return (sInt4)val;
90 : }
91 :
92 : /*****************************************************************************
93 : * mdl_LocalUnpack() --
94 : *
95 : * Arthur Taylor / MDL
96 : *
97 : * PURPOSE
98 : * Unpack the local use data assuming that it was packed using the MDL
99 : * encoder. This assumes "local" starts at octet 6 (i.e. skipping over the
100 : * length and section ID octets)
101 : *
102 : * In Section 2, GRIB2 provides for local use data. The MDL encoder packs
103 : * that data, and signifies that it has done so by setting octet 6 to 1.
104 : * This may have been inadvisable, since the GRIB2 specs did not state that
105 : * octet 6 was special.
106 : *
107 : * ARGUMENTS
108 : * local = The section 2 data prior to being unpacked. (Input)
109 : * locallen = The length of "local" (Input)
110 : * idat = Unpacked MDL local data (if it was an integer) (Output)
111 : * nidat = length of idat. (Input)
112 : * rdat = Unpacked MDL local data (if it was a float) (Output)
113 : * nrdat = length of rdat. (Input)
114 : *
115 : * FILES/DATABASES: None
116 : *
117 : * RETURNS: int
118 : * 0 = Ok.
119 : * 1 = Mixed local data types?
120 : * 2 = nrdat is not large enough.
121 : * 3 = nidat is not large enough.
122 : * 4 = Too many bits to unpack, should be a max of 32 bits.
123 : * 5 = Locallen is too small
124 : *
125 : * HISTORY
126 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
127 : *
128 : * NOTES
129 : *****************************************************************************
130 : */
131 : #define GRIB_UNSIGN_INT2(a,b) ((a<<8)+b)
132 0 : static int mdl_LocalUnpack(unsigned char *local, sInt4 locallen,
133 : sInt4 *idat, sInt4 *nidat, float *rdat,
134 : sInt4 *nrdat)
135 : {
136 0 : int BytesUsed = 0; /* How many bytes have been used. */
137 : unsigned int numGroup; /* Number of groups */
138 : int i; /* Counter over the groups. */
139 : sInt4 numVal; /* Number of values in a given group. */
140 : float refVal; /* The reference value in a given group. */
141 : unsigned int scale; /* The power of 10 scale value in the group. */
142 : sInt4 recScale10; /* 1 / 10**scale.. For faster computations. */
143 : unsigned char numBits; /* # of bits for a single element in the group. */
144 : char f_dataType; /* If the local data is a float or integer. */
145 0 : char f_firstType = 0; /* The type of the first group of local data. */
146 0 : int curIndex = 0; /* Where to store the current data. */
147 : int j; /* Counter over the number of values in a group. */
148 : size_t numUsed; /* Number of bytes used in call to memBitRead. */
149 : uChar bufLoc; /* Where to read for more bits in the data stream. */
150 : uInt4 uli_temp; /* Temporary storage to hold unpacked data. */
151 :
152 0 : if (locallen < BytesUsed + 3) {
153 : #ifdef DEBUG
154 0 : printf("Locallen is too small.\n");
155 : #endif
156 0 : return 5;
157 : }
158 : /* The calling routine should check octet 6, which is local[0], to be 1,
159 : * so we just assert it is 1. */
160 0 : myAssert(local[0] == 1);
161 0 : numGroup = GRIB_UNSIGN_INT2(local[1], local[2]);
162 0 : local += 3;
163 0 : BytesUsed += 3;
164 0 : myAssert(*nrdat > 1);
165 0 : myAssert(*nidat > 1);
166 0 : idat[0] = 0;
167 0 : rdat[0] = 0;
168 :
169 0 : for (i = 0; (unsigned int)i < numGroup; i++) {
170 0 : if (locallen < BytesUsed + 12) {
171 : #ifdef DEBUG
172 0 : printf("Locallen is too small.\n");
173 : #endif
174 0 : return 5;
175 : }
176 0 : MEMCPY_BIG(&numVal, local, sizeof (sInt4));
177 0 : MEMCPY_BIG(&refVal, local + 4, sizeof (float));
178 0 : scale = GRIB_UNSIGN_INT2(local[8], local[9]);
179 0 : recScale10 = (sInt4)(1 / pow (10.0, scale));
180 0 : numBits = local[10];
181 0 : if (numBits >= 32) {
182 : #ifdef DEBUG
183 0 : printf("Too many bits too unpack.\n");
184 : #endif
185 0 : return 4;
186 : }
187 0 : f_dataType = local[11];
188 0 : local += 12;
189 0 : BytesUsed += 12;
190 0 : if (locallen < BytesUsed + ((numBits * numVal) + 7) / 8) {
191 : #ifdef DEBUG
192 0 : printf("Locallen is too small.\n");
193 : #endif
194 0 : return 5;
195 : }
196 0 : if (i == 0) {
197 0 : f_firstType = f_dataType;
198 0 : } else if (f_firstType != f_dataType) {
199 : #ifdef DEBUG
200 0 : printf("Local use data has mixed float/integer type?\n");
201 : #endif
202 0 : return 1;
203 : }
204 0 : bufLoc = 8;
205 0 : if (f_dataType == 0) {
206 : /* Floating point data. */
207 0 : if (*nrdat < (curIndex + numVal + 3)) {
208 : #ifdef DEBUG
209 0 : printf("nrdat is not large enough.\n");
210 : #endif
211 0 : return 2;
212 : }
213 0 : rdat[curIndex] = numVal;
214 0 : curIndex++;
215 0 : rdat[curIndex] = scale;
216 0 : curIndex++;
217 0 : for (j = 0; j < numVal; j++) {
218 0 : memBitRead(&uli_temp, sizeof(sInt4), local, numBits,
219 : &bufLoc, &numUsed);
220 0 : local += numUsed;
221 0 : BytesUsed += (int) numUsed;
222 0 : rdat[curIndex] = (refVal + uli_temp) * recScale10;
223 0 : curIndex++;
224 : }
225 0 : rdat[curIndex] = 0;
226 : } else {
227 : /* Integer point data. */
228 0 : if (*nidat < (curIndex + numVal + 3)) {
229 : #ifdef DEBUG
230 0 : printf("nidat is not large enough.\n");
231 : #endif
232 0 : return 3;
233 : }
234 0 : idat[curIndex] = numVal;
235 0 : curIndex++;
236 0 : idat[curIndex] = scale;
237 0 : curIndex++;
238 0 : for (j = 0; j < numVal; j++) {
239 0 : memBitRead(&uli_temp, sizeof(sInt4), local, numBits,
240 : &bufLoc, &numUsed);
241 0 : local += numUsed;
242 0 : BytesUsed += (int) numUsed;
243 0 : idat[curIndex] = (sInt4)((refVal + uli_temp) * recScale10);
244 0 : curIndex++;
245 : }
246 0 : idat[curIndex] = 0;
247 : }
248 : }
249 0 : return 0;
250 : }
251 :
252 : /*****************************************************************************
253 : * fillOutSectLen() --
254 : *
255 : * Arthur Taylor / MDL
256 : *
257 : * PURPOSE
258 : * The GRIB2 API returns the lengths of each GRIB2 section along with the
259 : * meta data. NCEP's routines did not, so this routine fills in that part.
260 : * This routine has to consider which grid is being unpacked.
261 : *
262 : * c_ipack is passed in after section 1 (since section 1 is not repeated)
263 : *
264 : * ARGUMENTS
265 : * c_ipack = Complete GRIB2 message to look for the section lengths in. (In)
266 : * lenCpack = Length of c_ipack (Input)
267 : * subgNum = Which sub rid to find the section lengths of. (Input)
268 : * is2 = Section 2 data (Output)
269 : * is3 = Section 3 data (Output)
270 : * is4 = Section 4 data (Output)
271 : * is5 = Section 5 data (Output)
272 : * is6 = Section 6 data (Output)
273 : * is7 = Section 7 data (Output)
274 : *
275 : * FILES/DATABASES: None
276 : *
277 : * RETURNS: int
278 : * 0 = Ok.
279 : * 1 = c_ipack is not large enough
280 : * 2 = invalid section number
281 : *
282 : * HISTORY
283 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
284 : *
285 : * NOTES
286 : *****************************************************************************
287 : */
288 460 : static int fillOutSectLen(unsigned char *c_ipack, int lenCpack,
289 : int subgNum, sInt4 *is2, sInt4 *is3,
290 : sInt4 *is4, sInt4 *is5, sInt4 *is6, sInt4 *is7)
291 : {
292 460 : sInt4 offset = 0; /* How far in c_ipack we have read. */
293 : sInt4 sectLen; /* The length of the current section. */
294 : int sectId; /* The current section number. */
295 460 : unsigned int gNum = 0; /* Which sub group we are currently working with. */
296 :
297 460 : if (lenCpack < 5) {
298 : #ifdef DEBUG
299 0 : printf("Cpack is not large enough.\n");
300 : #endif
301 0 : return 1;
302 : }
303 : /* assert that we start with data in either section 2 or 3. */
304 460 : myAssert((c_ipack[4] == 2) || (c_ipack[4] == 3));
305 3196 : while (gNum <= (unsigned int)subgNum) {
306 2736 : if (lenCpack < offset + 5) {
307 : #ifdef DEBUG
308 0 : printf("Cpack is not large enough.\n");
309 : #endif
310 0 : return 1;
311 : }
312 2736 : MEMCPY_BIG(§Len, c_ipack + offset, sizeof(sInt4));
313 : /* Check if we just read section 8. If so, then it is "7777" =
314 : * 926365495 regardless of endian'ness. */
315 2736 : if (sectLen == 926365495L) {
316 : #ifdef DEBUG
317 0 : printf("Shouldn't see sect 8. Should stop after correct sect 7\n");
318 : #endif
319 0 : return 2;
320 : }
321 2736 : sectId = c_ipack[offset + 4];
322 2736 : switch (sectId) {
323 424 : case 2:
324 424 : is2[0] = sectLen;
325 424 : break;
326 460 : case 3:
327 460 : is3[0] = sectLen;
328 460 : break;
329 463 : case 4:
330 463 : is4[0] = sectLen;
331 463 : break;
332 463 : case 5:
333 463 : is5[0] = sectLen;
334 463 : break;
335 463 : case 6:
336 463 : is6[0] = sectLen;
337 463 : break;
338 463 : case 7:
339 463 : is7[0] = sectLen;
340 463 : gNum++;
341 463 : break;
342 0 : default:
343 : #ifdef DEBUG
344 0 : printf("Invalid section id %d.\n", sectId);
345 : #endif
346 0 : return 2;
347 : }
348 2736 : offset += sectLen;
349 : }
350 460 : return 0;
351 : }
352 :
353 : /*****************************************************************************
354 : * TransferInt() --
355 : *
356 : * Arthur Taylor / MDL
357 : *
358 : * PURPOSE
359 : * To transfer the data from "fld" and "bmap" to "iain" and "ib". The API
360 : * attempts to rearrange the order, so that anything returned from it has
361 : * scan mode 0100????
362 : *
363 : * ARGUMENTS
364 : * fld = The expanded grid from NCEPs routines (Input)
365 : * ngrdpts = Length of fld (Input)
366 : * ibitmap = flag if we have a bitmap or not. (Input)
367 : * bmap = bitmap from NCEPs routines. (Input)
368 : * f_ignoreScan = Flag to ignore the attempt at changing the scan (Input)
369 : * scan = The scan orientation of fld/bmap/iain/ib (Input/Output)
370 : * nx, ny = The dimensions of the grid. (Input)
371 : * iclean = 1 means the user wants the unpacked data returned without
372 : * missing values in it. 0 means embed the missing values. (In)
373 : * xmissp = The primary missing value to use if iclean = 0. (Input).
374 : * iain = The grid to return. (Output)
375 : * nd2x3 = The length of iain. (Input)
376 : * ib = Bitmap if it was packed (Output)
377 : *
378 : * FILES/DATABASES: None
379 : *
380 : * RETURNS: int
381 : * 0 = Ok.
382 : * 1 = nd2x3 is too small.
383 : * 2 = nx*ny != ngrdpts
384 : *
385 : * HISTORY
386 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
387 : *
388 : * NOTES
389 : * May want to disable the scan adjustment in the future.
390 : *****************************************************************************
391 : */
392 67 : static int TransferInt(float *fld, sInt4 ngrdpts, sInt4 ibitmap,
393 : sInt4 *bmap, char f_ignoreScan, sInt4 *scan,
394 : sInt4 nx, sInt4 ny, sInt4 iclean, float xmissp,
395 : sInt4 *iain, sInt4 nd2x3, sInt4 *ib)
396 : {
397 : int i; /* loop counter over all grid points. */
398 : sInt4 x, y; /* Where we are in a grid of scan value 0100???? */
399 : int curIndex; /* Where in iain to store the current data. */
400 :
401 67 : if (nd2x3 < ngrdpts) {
402 : #ifdef DEBUG
403 0 : printf("nd2x3(%d) is < ngrdpts(%d)\n", nd2x3, ngrdpts);
404 : #endif
405 0 : return 1;
406 : }
407 67 : if (f_ignoreScan || ((*scan & 0xf0) == 64)) {
408 67 : if (ibitmap) {
409 0 : for (i = 0; i < ngrdpts; i++) {
410 0 : ib[i] = bmap[i];
411 : /* Check if we are supposed to insert xmissp into the field */
412 0 : if ((iclean != 0) && (ib[i] == 0)) {
413 0 : iain[i] = FloatToSInt4Clamp(xmissp);
414 : } else {
415 0 : iain[i] = FloatToSInt4Clamp(fld[i]);
416 : }
417 : }
418 : } else {
419 546457 : for (i = 0; i < ngrdpts; i++) {
420 546390 : iain[i] = FloatToSInt4Clamp(fld[i]);
421 : }
422 : }
423 : } else {
424 0 : if( nx <= 0 || ny <= 0 )
425 : {
426 0 : return 1;
427 : }
428 0 : if (ny != ngrdpts / nx) {
429 : #ifdef DEBUG
430 0 : printf("nx(%d) * ny(%d) != ngrdpts(%d)\n", nx, ny, ngrdpts);
431 : #endif
432 0 : return 2;
433 : }
434 0 : if (ibitmap) {
435 0 : for (i = 0; i < ngrdpts; i++) {
436 0 : ScanIndex2XY(i, &x, &y, *scan, nx, ny);
437 : /* ScanIndex returns value as if scan was 0100(0000) */
438 0 : curIndex = (x - 1) + (y - 1) * nx;
439 0 : if( curIndex < 0 || curIndex >= nd2x3 )
440 0 : return 1;
441 0 : ib[curIndex] = bmap[i];
442 : /* Check if we are supposed to insert xmissp into the field */
443 0 : if ((iclean != 0) && (ib[curIndex] == 0)) {
444 0 : iain[i] = FloatToSInt4Clamp(xmissp);
445 : } else {
446 0 : iain[curIndex] = FloatToSInt4Clamp(fld[i]);
447 : }
448 : }
449 : } else {
450 0 : for (i = 0; i < ngrdpts; i++) {
451 0 : ScanIndex2XY (i, &x, &y, *scan, nx, ny);
452 : /* ScanIndex returns value as if scan was 0100(0000) */
453 0 : curIndex = (x - 1) + (y - 1) * nx;
454 0 : if( curIndex < 0 || curIndex >= nd2x3 )
455 0 : return 1;
456 0 : iain[curIndex] = FloatToSInt4Clamp(fld[i]);
457 : }
458 : }
459 0 : *scan = 64 + (*scan & 0x0f);
460 : }
461 67 : return 0;
462 : }
463 :
464 : /*****************************************************************************
465 : * TransferFloat() --
466 : *
467 : * Arthur Taylor / MDL
468 : *
469 : * PURPOSE
470 : * To transfer the data from "fld" and "bmap" to "ain" and "ib". The API
471 : * attempts to rearrange the order, so that anything returned from it has
472 : * scan mode 0100????
473 : *
474 : * ARGUMENTS
475 : * fld = The expanded grid from NCEPs routines (Input)
476 : * ngrdpts = Length of fld (Input)
477 : * ibitmap = flag if we have a bitmap or not. (Input)
478 : * bmap = bitmap from NCEPs routines. (Input)
479 : * scan = The scan orientation of fld/bmap/iain/ib (Input/Output)
480 : * f_ignoreScan = Flag to ignore the attempt at changing the scan (Input)
481 : * nx, ny = The dimensions of the grid. (Input)
482 : * iclean = 1 means the user wants the unpacked data returned without
483 : * missing values in it. 0 means embed the missing values. (In)
484 : * xmissp = The primary missing value to use if iclean = 0. (Input).
485 : * ain = The grid to return. (Output)
486 : * nd2x3 = The length of iain. (Input)
487 : * ib = Bitmap if it was packed (Output)
488 : *
489 : * FILES/DATABASES: None
490 : *
491 : * RETURNS: int
492 : * 0 = Ok.
493 : * 1 = nd2x3 is too small.
494 : * 2 = nx*ny != ngrdpts
495 : *
496 : * HISTORY
497 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
498 : *
499 : * NOTES
500 : * May want to disable the scan adjustment in the future.
501 : *****************************************************************************
502 : */
503 90 : static int TransferFloat(float *fld, sInt4 ngrdpts, sInt4 ibitmap,
504 : sInt4 *bmap, char f_ignoreScan, sInt4 *scan,
505 : sInt4 nx, sInt4 ny, sInt4 iclean, float xmissp,
506 : float *ain, sInt4 nd2x3, sInt4 *ib)
507 : {
508 : int i; /* loop counter over all grid points. */
509 : sInt4 x, y; /* Where we are in a grid of scan value 0100???? */
510 : uInt4 curIndex; /* Where in ain to store the current data. */
511 :
512 90 : if (nd2x3 < ngrdpts) {
513 : #ifdef DEBUG
514 0 : printf("nd2x3(%d) is < ngrdpts(%d)\n", nd2x3, ngrdpts);
515 : #endif
516 0 : return 1;
517 : }
518 90 : if (f_ignoreScan || ((*scan & 0xf0) == 64)) {
519 76 : if( ain ) {
520 76 : if (ibitmap) {
521 802 : for (i = 0; i < ngrdpts; i++) {
522 800 : ib[i] = bmap[i];
523 : /* Check if we are supposed to insert xmissp into the field */
524 800 : if ((iclean != 0) && (ib[i] == 0)) {
525 0 : ain[i] = xmissp;
526 : } else {
527 800 : ain[i] = fld[i];
528 : }
529 : }
530 : } else {
531 2216200 : for (i = 0; i < ngrdpts; i++) {
532 2216130 : ain[i] = fld[i];
533 : }
534 : }
535 : }
536 : } else {
537 14 : if( nx <= 0 || ny <= 0 )
538 : {
539 0 : return 1;
540 : }
541 14 : if (ny != ngrdpts / nx) {
542 : #ifdef DEBUG
543 0 : printf("nx(%d) * ny(%d) != ngrdpts(%d)\n", nx, ny, ngrdpts);
544 : #endif
545 0 : return 2;
546 : }
547 14 : if (ibitmap) {
548 0 : for (i = 0; i < ngrdpts; i++) {
549 0 : ScanIndex2XY(i, &x, &y, *scan, nx, ny);
550 0 : if (x <= 0 || y <= 0 )
551 0 : return 1;
552 : /* ScanIndex returns value as if scan was 0100(0000) */
553 0 : curIndex = (uInt4)(x - 1) + (uInt4)(y - 1) * (uInt4)nx;
554 0 : if (curIndex >= (uInt4)nd2x3)
555 0 : return 1;
556 0 : ib[curIndex] = bmap[i];
557 0 : if( ain )
558 : {
559 : /* Check if we are supposed to insert xmissp into the field */
560 0 : if ((iclean != 0) && (ib[curIndex] == 0)) {
561 0 : ain[i] = xmissp;
562 : } else {
563 0 : ain[curIndex] = fld[i];
564 : }
565 : }
566 : }
567 : } else {
568 5408080 : for (i = 0; i < ngrdpts; i++) {
569 5408060 : ScanIndex2XY(i, &x, &y, *scan, nx, ny);
570 5408060 : if (x <= 0 || y <= 0 )
571 0 : return 1;
572 : /* ScanIndex returns value as if scan was 0100(0000) */
573 5408060 : curIndex = (uInt4)(x - 1) + (uInt4)(y - 1) * (uInt4)nx;
574 5408060 : if( curIndex >= (uInt4)nd2x3 )
575 0 : return 1;
576 5408060 : if( ain )
577 5408060 : ain[curIndex] = fld[i];
578 : }
579 : }
580 14 : *scan = 64 + (*scan & 0x0f);
581 : }
582 : /*
583 : if (1==1) {
584 : int i;
585 : for (i=0; i < ngrdpts; i++) {
586 : if (ain[i] != 9999) {
587 : fprintf (stderr, "grib2api.c: ain[%d] %f, fld[%d] %f\n", i, ain[i],
588 : i, fld[i]);
589 : }
590 : }
591 : }
592 : */
593 90 : return 0;
594 : }
595 :
596 : #ifdef PKNCEP
597 : int pk_g2ncep(sInt4 *kfildo, float *ain, sInt4 *iain, sInt4 *nx,
598 : sInt4 *ny, sInt4 *idat, sInt4 *nidat, float *rdat,
599 : sInt4 *nrdat, sInt4 *is0, sInt4 *ns0, sInt4 *is1,
600 : sInt4 *ns1, sInt4 *is3, sInt4 *ns3, sInt4 *is4,
601 : sInt4 *ns4, sInt4 *is5, sInt4 *ns5, sInt4 *is6,
602 : sInt4 *ns6, sInt4 *is7, sInt4 *ns7, sInt4 *ib,
603 : sInt4 *ibitmap, unsigned char *cgrib, sInt4 *nd5,
604 : sInt4 *missp, float *xmissp, sInt4 *misss,
605 : float *xmisss, sInt4 *inew, sInt4 *minpk, sInt4 *iclean,
606 : sInt4 *l3264b, sInt4 *jer, sInt4 *ndjer, sInt4 *kjer)
607 : {
608 : g2int listsec0[2];
609 : g2int listsec1[13];
610 : g2int igds[5];
611 : g2int igdstmpl[];
612 : int ierr; /* Holds the error code from a called routine. */
613 : int i;
614 :
615 : listsec0[0] = is0[6];
616 : listsec0[1] = is0[7];
617 :
618 : listsec1[0] = is1[5];
619 : listsec1[1] = is1[7];
620 : listsec1[2] = is1[9];
621 : listsec1[3] = is1[10];
622 : listsec1[4] = is1[11];
623 : listsec1[5] = is1[12];
624 : listsec1[6] = is1[14];
625 : listsec1[7] = is1[15];
626 : listsec1[8] = is1[16];
627 : listsec1[9] = is1[17];
628 : listsec1[10] = is1[18];
629 : listsec1[11] = is1[19];
630 : listsec1[12] = is1[20];
631 :
632 : ierr = g2_create(cgrib, listsec0, listsec1);
633 : printf("Length = %d\n", ierr);
634 :
635 : if ((idat[0] != 0) || (rdat[0] != 0)) {
636 : printf("Don't handle this yet.\n");
637 : /*
638 : ierr = g2_addlocal (cgrib, unsigned char *csec2, g2int lcsec2);
639 : */
640 : }
641 :
642 : igds[0] = is3[5];
643 : igds[1] = is3[6];
644 : igds[2] = is3[10];
645 : igds[3] = is3[11];
646 : igds[4] = is3[12];
647 :
648 : IS3(15) - IS3(nn) = Grid Definition Template, stored in bytes 15 - nn(*)
649 :
650 : ierr = g2_addgrid(cgrib, igds, g2int * igdstmpl, g2int * ideflist,
651 : g2int idefnum)
652 :
653 :
654 :
655 :
656 : return 0;
657 : /*
658 :
659 : To start a new GRIB2 message, call function g2_create. G2_create
660 : encodes Sections 0 and 1 at the beginning of the message. This routine
661 : must be used to create each message.
662 :
663 : Routine g2_addlocal can be used to add a Local Use Section ( Section 2 ).
664 : Note that this section is optional and need not appear in a GRIB2 message.
665 :
666 : Function g2_addgrid is used to encode a grid definition into Section 3.
667 : This grid definition defines the geometry of the data values in the
668 : fields that follow it. g2_addgrid can be called again to change the grid
669 : definition describing subsequent data fields.
670 :
671 : Each data field is added to the GRIB2 message using routine g2_addfield,
672 : which adds Sections 4, 5, 6, and 7 to the message.
673 :
674 : After all desired data fields have been added to the GRIB2 message, a
675 : call to function g2_gribend is needed to add the final section 8 to the
676 : message and to update the length of the message. A call to g2_gribend
677 : is required for each GRIB2 message.
678 : */
679 :
680 : }
681 : #endif
682 :
683 : /*****************************************************************************
684 : * unpk_g2ncep() --
685 : *
686 : * Arthur Taylor / MDL
687 : *
688 : * PURPOSE
689 : * This procedure is a wrapper around the NCEP GRIB2 routines to interface
690 : * their routines with the "official NWS" GRIB2 API. The reason for this is
691 : * so drivers that have been written to use the "official NWS" GRIB2 API, can
692 : * use the NCEP library with minimal disruption.
693 : *
694 : * ARGUMENTS
695 : * kfildo = Unit number for output diagnostics (C ignores this). (Input)
696 : * ain = contains the data if the original data was float data.
697 : * (size = nd2x3) (Output)
698 : * iain = contains the data if the original data was integer data.
699 : * (size = nd2x3) (Output)
700 : * nd2x3 = length of ain, iain, ib (Input)
701 : * (at least the size of num grid points)
702 : * idat = local use data if any that were unpacked from Section 2. (Output)
703 : * nidat = length of idat. (Input)
704 : * rdat = floating point local use data (Output)
705 : * nrdat = length of rdat. (Input)
706 : * is0 = Section 0 data (Output)
707 : * ns0 = length of is0 (16 is fine) (Input)
708 : * is1 = Section 1 data (Output)
709 : * ns1 = length of is1 (21 is fine) (Input)
710 : * is2 = Section 2 data (Output)
711 : * ns2 = length of is2 () (Input)
712 : * is3 = Section 3 data (Output)
713 : * ns3 = length of is3 (96 or 1600) (Input)
714 : * is4 = Section 4 data (Output)
715 : * ns4 = length of is4 (60) (Input)
716 : * is5 = Section 5 data (Output)
717 : * ns5 = length of is5 (49 is fine) (Input)
718 : * is6 = Section 6 data (Output)
719 : * ns6 = length of is6 (6 is fine) (Input)
720 : * is7 = Section 7 data (Output)
721 : * ns7 = length of is7 (8 is fine) (Input)
722 : * ib = Bitmap if user requested it, and it was packed (Output)
723 : * ibitmap = 0 means ib is invalid, 1 means ib is valid. (Output)
724 : * c_ipack = The message to unpack (Input)
725 : * nd5 = 1/4 the size of c_ipack (Input)
726 : * xmissp = The floating point representation for the primary missing value.
727 : * (Input/Output)
728 : * xmisss = The floating point representation for the secondary missing
729 : * value (Input/Output)
730 : * new = 1 if this is the first grid to be unpacked, else 0. (Input)
731 : * iclean = 1 means the user wants the unpacked data returned without
732 : * missing values in it. 0 means embed the missing values. (Input)
733 : * l3264b = Integer word length in bits (32 or 64) (Input)
734 : * iendpk = A 1 means no more grids in this message, a 0 means more grids.
735 : * (Output)
736 : * jer(ndjer,2) = error codes along with severity. (Output)
737 : * ndjer = 1/2 length of jer. (>= 15) (Input)
738 : * kjer = number of error messages stored in jer.
739 : *
740 : * FILES/DATABASES: None
741 : *
742 : * RETURNS: void
743 : * "jer" is used to store error messages, and kjer is used to denote how
744 : * many there were. Jer is set up as a 2 column array, with the first
745 : * column in the first half of the array, and denoting the GRIB2 section
746 : * an error occurred in. The second column denoting the severity.
747 : *
748 : * The possibilities from unpk_g2ncep() are as follows:
749 : * ker=1 jer[0,0]=0 jer[0,1]=2: Message is not formed correctly
750 : * Request for an invalid subgrid
751 : * Problems unpacking the data.
752 : * problems expanding the data.
753 : * Calling dimensions were too small.
754 : * ker=2 jer[1,0]=100 jer[1,1]=2: Error unpacking section 1.
755 : * ker=3 jer[2,0]=200 jer[2,1]=2: Error unpacking section 2.
756 : * ker=4 jer[3,0]=300 jer[3,1]=2: Error unpacking section 3.
757 : * ker=5 jer[4,0]=400 jer[4,1]=2: Error unpacking section 4.
758 : * ker=6 jer[5,0]=500 jer[5,1]=2: Error unpacking section 5.
759 : * Data Template not implemented.
760 : * During Transfer, nx * ny != ngrdpts.
761 : * ker=7 jer[6,0]=600 jer[6,1]=2: Error unpacking section 6.
762 : * ker=8 jer[7,0]=700 jer[7,1]=2: Error unpacking section 7.
763 : * ker=9 jer[8,0]=2001 jer[8,1]=2: nd2x3 is not large enough.
764 : * ker=9 jer[8,0]=2003 jer[8,1]=2: undefined sect 3 template
765 : * (see gridtemplates.h).
766 : * ker=9 jer[8,0]=2004 jer[8,1]=2: undefined sect 4 template
767 : * (see pdstemplates.h).
768 : * ker=9 jer[8,0]=2005 jer[8,1]=2: undefined sect 5 template
769 : * (see drstemplates.h).
770 : * ker=9 jer[8,0]=9999 jer[8,1]=2: NCEP returns an unrecognized error.
771 : *
772 : * HISTORY
773 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
774 : *
775 : * NOTES
776 : * MDL handles is5[12], is5[23], and is5[27] in an "interesting" manner.
777 : * MDL attempts to always return grids in scan mode 0100????
778 : *
779 : * ToDo: Check length of parameters better.
780 : * ToDo: Probably shouldn't abort if they have problems expanding the data.
781 : * ToDo: Better handling of:
782 : * gfld->list_opt = (Used if gfld->numoct_opt .ne. 0) This array contains
783 : * the number of grid points contained in each row (or column). (part of
784 : * Section 3) This element is a pointer to an array that holds the data.
785 : * This pointer is nullified if gfld->numoct_opt=0.
786 : * gfld->num_opt = (Used if gfld->numoct_opt .ne. 0) The number of entries in
787 : * array ideflist. i.e. number of rows (or columns) for which optional
788 : * grid points are defined. This value is set to zero, if
789 : * gfld->numoct_opt=0.
790 : * gfld->coord_list = Real array containing floating point values intended to
791 : * document the vertical discretisation associated to model data on
792 : * hybrid coordinate vertical levels. (part of Section 4) This element
793 : * is a pointer to an array that holds the data.
794 : * gfld->num_coord = number of values in array gfld->coord_list[].
795 : *****************************************************************************
796 : */
797 460 : void unpk_g2ncep(CPL_UNUSED sInt4 *kfildo, float *ain, sInt4 *iain, sInt4 *nd2x3,
798 : sInt4 *idat, sInt4 *nidat, float *rdat, sInt4 *nrdat,
799 : sInt4 *is0, CPL_UNUSED sInt4 *ns0, sInt4 *is1, CPL_UNUSED sInt4 *ns1,
800 : sInt4 *is2, sInt4 *ns2, sInt4 *is3, sInt4 *ns3,
801 : sInt4 *is4, sInt4 *ns4, sInt4 *is5, CPL_UNUSED sInt4 *ns5,
802 : sInt4 *is6, CPL_UNUSED sInt4 *ns6, sInt4 *is7, CPL_UNUSED sInt4 *ns7,
803 : sInt4 *ib, sInt4 *ibitmap, unsigned char *c_ipack,
804 : sInt4 *nd5, float *xmissp, float *xmisss,
805 : sInt4 *inew, sInt4 *iclean, CPL_UNUSED sInt4 *l3264b,
806 : sInt4 *iendpk, sInt4 *jer, sInt4 *ndjer, sInt4 *kjer)
807 : {
808 : int i; /* A counter used for a number of purposes. */
809 : static unsigned int subgNum = 0; /* The sub grid we read most recently.
810 : * This is primarily to help with the
811 : * inew option. */
812 : int ierr; /* Holds the error code from a called routine. */
813 : sInt4 listsec0[3];
814 : sInt4 listsec1[13];
815 : static sInt4 numfields = 1; /* Number of sub Grids in this message */
816 : sInt4 numlocal; /* Number of local sections in this message. */
817 : int unpack; /* Tell g2_getfld to unpack the message. */
818 : int expand; /* Tell g2_getflt to attempt to expand the bitmap. */
819 : gribfield *gfld; /* Holds the data after g2_getfld unpacks it. */
820 : sInt4 gridIndex; /* index in templatesgrid[] for this sect 3 templat */
821 : sInt4 pdsIndex; /* index in templatespds[] for this sect 4 template */
822 : sInt4 drsIndex; /* index in templatesdrs[] for this sect 5 template */
823 : int curIndex; /* Where in is3, is4, or is5 to store meta data */
824 : int scanIndex; /* Where in is3 to find the scan mode. */
825 : int nxIndex; /* Where in is3 to find the number of x values. */
826 : int nyIndex; /* Where in is3 to find the number of y values. */
827 : float f_temp; /* Assist with handling weird MDL behavior in is5[] */
828 : char f_ignoreScan; /* Flag to ignore the attempt at changing the scan */
829 : sInt4 dummyScan; /* Dummy place holder for call to Transfer routines
830 : * if ignoring scan. */
831 460 : const struct gridtemplate *templatesgrid = get_templatesgrid();
832 :
833 460 : myAssert(*ndjer >= 8);
834 : /* Init the error handling array. */
835 460 : memset((void *)jer, 0, 2 * *ndjer * sizeof(sInt4));
836 4140 : for (i = 0; i < 8; i++) {
837 3680 : jer[i] = i * 100;
838 : }
839 460 : *kjer = 8;
840 :
841 : /* The first time in, figure out how many grids there are, and store it in
842 : * numfields for subsequent calls with inew != 1. */
843 460 : if (*inew == 1) {
844 457 : subgNum = 0;
845 457 : ierr = g2_info(c_ipack, listsec0, listsec1, &numfields, &numlocal);
846 457 : if (ierr != 0) {
847 0 : switch (ierr) {
848 0 : case 1: /* Beginning characters "GRIB" not found. */
849 : case 2: /* GRIB message is not Edition 2. */
850 : case 3: /* Could not find Section 1, where expected. */
851 : case 4: /* End string "7777" found, but not where expected. */
852 : case 5: /* End string "7777" not found at end of message. */
853 : case 6: /* Invalid section number found. */
854 0 : jer[0 + *ndjer] = 2;
855 0 : *kjer = 1;
856 0 : break;
857 0 : default:
858 0 : jer[8 + *ndjer] = 2;
859 0 : jer[8] = 9999; /* Unknown error message. */
860 0 : *kjer = 9;
861 : }
862 303 : return;
863 : }
864 : } else {
865 3 : if (subgNum + 1 >= (unsigned int)numfields) {
866 : /* Field request error. */
867 0 : jer[0 + *ndjer] = 2;
868 0 : *kjer = 1;
869 0 : return;
870 : }
871 3 : subgNum++;
872 : }
873 :
874 : /* Expand the desired subgrid. */
875 460 : unpack = ain != NULL || iain != NULL;
876 460 : expand = 1;
877 : /* The size of c_ipack is *nd5 * sizeof(sInt4) */
878 460 : ierr = g2_getfld(c_ipack, *nd5 * sizeof(sInt4), subgNum + 1, unpack, expand, &gfld);
879 460 : if (ierr != 0) {
880 0 : switch (ierr) {
881 0 : case 1: /* Beginning characters "GRIB" not found. */
882 : case 2: /* GRIB message is not Edition 2. */
883 : case 3: /* The data field request number was not positive. */
884 : case 4: /* End string "7777" found, but not where expected. */
885 : case 6: /* message did not contain requested # of fields */
886 : case 7: /* End string "7777" not found at end of message. */
887 : case 8: /* Unrecognized Section encountered. */
888 0 : jer[0 + *ndjer] = 2;
889 0 : *kjer = 1;
890 0 : break;
891 0 : case 15: /* Error unpacking Section 1. */
892 0 : jer[1 + *ndjer] = 2;
893 0 : *kjer = 2;
894 0 : break;
895 0 : case 16: /* Error unpacking Section 2. */
896 0 : jer[2 + *ndjer] = 2;
897 0 : *kjer = 3;
898 0 : break;
899 0 : case 10: /* Error unpacking Section 3. */
900 0 : jer[3 + *ndjer] = 2;
901 0 : *kjer = 4;
902 0 : break;
903 0 : case 11: /* Error unpacking Section 4. */
904 0 : jer[4 + *ndjer] = 2;
905 0 : *kjer = 5;
906 0 : break;
907 0 : case 9: /* Data Template 5.NN not implemented. */
908 : case 12: /* Error unpacking Section 5. */
909 0 : jer[5 + *ndjer] = 2;
910 0 : *kjer = 6;
911 0 : break;
912 0 : case 13: /* Error unpacking Section 6. */
913 0 : jer[6 + *ndjer] = 2;
914 0 : *kjer = 7;
915 0 : break;
916 0 : case 14: /* Error unpacking Section 7. */
917 0 : jer[7 + *ndjer] = 2;
918 0 : *kjer = 8;
919 0 : break;
920 0 : default:
921 0 : jer[8 + *ndjer] = 2;
922 0 : jer[8] = 9999; /* Unknown error message. */
923 0 : *kjer = 9;
924 0 : break;
925 : }
926 0 : g2_free(gfld);
927 0 : return;
928 : }
929 : /* Check if data was not unpacked. */
930 460 : if (unpack && !gfld->unpacked) {
931 0 : jer[0 + *ndjer] = 2;
932 0 : *kjer = 1;
933 0 : g2_free(gfld);
934 0 : return;
935 : }
936 :
937 : /* Start going through the gfld structure and converting it to the needed
938 : * data output formats. */
939 460 : myAssert(*ns0 >= 16);
940 460 : MEMCPY_BIG(&(is0[0]), c_ipack, sizeof(sInt4));
941 460 : is0[6] = gfld->discipline;
942 460 : is0[7] = gfld->version;
943 460 : MEMCPY_BIG(&(is0[8]), c_ipack + 8, sizeof(sInt4));
944 : /* The following assert fails only if the GRIB message is more that 4
945 : * giga-bytes large, which I think would break the fortran library. */
946 460 : myAssert(is0[8] == 0);
947 460 : MEMCPY_BIG(&(is0[8]), c_ipack + 12, sizeof(sInt4));
948 :
949 460 : myAssert(*ns1 >= 21);
950 460 : myAssert(gfld->idsectlen >= 13);
951 460 : MEMCPY_BIG(&(is1[0]), c_ipack + 16, sizeof(sInt4));
952 460 : is1[4] = c_ipack[20];
953 460 : is1[5] = gfld->idsect[0];
954 460 : is1[7] = gfld->idsect[1];
955 460 : is1[9] = gfld->idsect[2];
956 460 : is1[10] = gfld->idsect[3];
957 460 : is1[11] = gfld->idsect[4];
958 460 : is1[12] = gfld->idsect[5]; /* Year */
959 460 : is1[14] = gfld->idsect[6]; /* Month */
960 460 : is1[15] = gfld->idsect[7]; /* Day */
961 460 : is1[16] = gfld->idsect[8]; /* Hour */
962 460 : is1[17] = gfld->idsect[9]; /* Min */
963 460 : is1[18] = gfld->idsect[10]; /* Sec */
964 460 : is1[19] = gfld->idsect[11];
965 460 : is1[20] = gfld->idsect[12];
966 :
967 : /* Fill out section lengths (separate procedure because of possibility of
968 : * having multiple grids. Should combine fillOutSectLen g2_info, and
969 : * g2_getfld into one procedure to optimize it. */
970 460 : fillOutSectLen(c_ipack + 16 + is1[0], 4 * *nd5 - 15 - is1[0], subgNum,
971 : is2, is3, is4, is5, is6, is7);
972 :
973 : /* Check if there is section 2 data. */
974 460 : if (gfld->locallen > 0) {
975 : /* The + 1 is so we don't overwrite the section length */
976 2 : memset((void *)(is2 + 1), 0, (*ns2 - 1) * sizeof(sInt4));
977 2 : is2[4] = 2;
978 2 : is2[5] = gfld->local[0];
979 : /* check if MDL Local use simple packed data */
980 2 : if (is2[5] == 1) {
981 0 : mdl_LocalUnpack(gfld->local, gfld->locallen, idat, nidat, rdat,
982 : nrdat);
983 : } else {
984 : /* local use section was not MDL packed, return it in is2. */
985 26 : for (i = 0; i < gfld->locallen; i++) {
986 24 : is2[i + 5] = gfld->local[i];
987 : }
988 : }
989 : } else {
990 : /* API specified that is2[0] = 0; idat[0] = 0, and rdat[0] = 0 */
991 458 : is2[0] = 0;
992 458 : idat[0] = 0;
993 458 : rdat[0] = 0;
994 : }
995 :
996 460 : is3[4] = 3;
997 460 : is3[5] = gfld->griddef;
998 460 : is3[6] = gfld->ngrdpts;
999 460 : if (*nd2x3 < gfld->ngrdpts) {
1000 0 : jer[8 + *ndjer] = 2;
1001 0 : jer[8] = 2001; /* nd2x3 is not large enough */
1002 0 : *kjer = 9;
1003 0 : g2_free(gfld);
1004 0 : return;
1005 : }
1006 460 : is3[10] = gfld->numoct_opt;
1007 460 : is3[11] = gfld->interp_opt;
1008 460 : is3[12] = gfld->igdtnum;
1009 460 : gridIndex = getgridindex(gfld->igdtnum);
1010 460 : if (gridIndex == -1) {
1011 0 : jer[8 + *ndjer] = 2;
1012 0 : jer[8] = 2003; /* undefined sect 3 template */
1013 0 : *kjer = 9;
1014 0 : g2_free(gfld);
1015 0 : return;
1016 : }
1017 460 : if( gfld->igdtlen > templatesgrid[gridIndex].mapgridlen )
1018 : {
1019 0 : jer[8 + *ndjer] = 2;
1020 0 : jer[8] = 2003; /* undefined sect 3 template: FIXME: wrong code probably */
1021 0 : *kjer = 9;
1022 0 : g2_free (gfld);
1023 0 : return;
1024 : }
1025 460 : curIndex = 14;
1026 9817 : for (i = 0; i < gfld->igdtlen; i++) {
1027 9357 : if( curIndex < 0 || curIndex >= *ns3 )
1028 : {
1029 0 : jer[8 + *ndjer] = 2;
1030 0 : jer[8] = 2003; /* undefined sect 3 template: FIXME: wrong code probably */
1031 0 : *kjer = 9;
1032 0 : g2_free (gfld);
1033 0 : return;
1034 : }
1035 9357 : is3[curIndex] = gfld->igdtmpl[i];
1036 9357 : curIndex += abs(templatesgrid[gridIndex].mapgrid[i]);
1037 : }
1038 : /* API attempts to return grid in scan mode 0100????. Find the necessary
1039 : * indexes into the is3 array for the attempt. */
1040 460 : switch (gfld->igdtnum) {
1041 186 : case 0:
1042 : case 1:
1043 : case 2:
1044 : case 3:
1045 : case 40:
1046 : case 41:
1047 : case 42:
1048 : case 43:
1049 186 : scanIndex = 72 - 1;
1050 186 : nxIndex = 31 - 1;
1051 186 : nyIndex = 35 - 1;
1052 186 : break;
1053 241 : case 10:
1054 : case 12: // Transverse Mercator
1055 241 : scanIndex = 60 - 1;
1056 241 : nxIndex = 31 - 1;
1057 241 : nyIndex = 35 - 1;
1058 241 : break;
1059 26 : case 20:
1060 : case 30:
1061 : case 31:
1062 26 : scanIndex = 65 - 1;
1063 26 : nxIndex = 31 - 1;
1064 26 : nyIndex = 35 - 1;
1065 26 : break;
1066 0 : case 90:
1067 0 : scanIndex = 64 - 1;
1068 0 : nxIndex = 31 - 1;
1069 0 : nyIndex = 35 - 1;
1070 0 : break;
1071 0 : case 110:
1072 0 : scanIndex = 57 - 1;
1073 0 : nxIndex = 31 - 1;
1074 0 : nyIndex = 35 - 1;
1075 0 : break;
1076 7 : case 50:
1077 : case 51:
1078 : case 52:
1079 : case 53:
1080 : case 100:
1081 : case 120:
1082 : case 1000:
1083 : case 1200:
1084 : default:
1085 7 : scanIndex = -1;
1086 7 : nxIndex = -1;
1087 7 : nyIndex = -1;
1088 : }
1089 :
1090 460 : is4[4] = 4;
1091 460 : is4[5] = gfld->num_coord;
1092 460 : is4[7] = gfld->ipdtnum;
1093 460 : pdsIndex = getpdsindex(gfld->ipdtnum);
1094 : #if 0
1095 : if (pdsIndex == -1) {
1096 : jer[8 + *ndjer] = 2;
1097 : jer[8] = 2004; /* undefined sect 4 template */
1098 : *kjer = 9;
1099 : g2_free(gfld);
1100 : return;
1101 : }
1102 : #endif
1103 460 : if( pdsIndex >= 0 )
1104 : {
1105 457 : curIndex = 9;
1106 8056 : for (i = 0; i < gfld->ipdtlen; i++) {
1107 7599 : const struct pdstemplate *templatespds = get_templatespds();
1108 7599 : if( curIndex >= *ns4 )
1109 : {
1110 : /* Should we error out instead ? */
1111 0 : break;
1112 : }
1113 7599 : is4[curIndex] = gfld->ipdtmpl[i];
1114 7599 : if( i == MAXPDSMAPLEN )
1115 : {
1116 0 : jer[8 + *ndjer] = 2;
1117 0 : jer[8] = 2004; /* undefined sect 4 template */
1118 0 : *kjer = 9;
1119 0 : g2_free (gfld);
1120 0 : return;
1121 : }
1122 : /* 2020-08-04 Taylor: In 2015-10-07, Vuong added the ability for the PDS
1123 : * template to have negative forecast times, so data encoded with the
1124 : * update has -1 * abs(value). Unfortunately negative values encoded via
1125 : * older versions used 2s complement. So the updated decoder has to deal
1126 : * with both correctly and incorrectly decoded neg values.
1127 : *
1128 : * To resolve this, negative numbers that are < -1 * 2^30-1 are assumed
1129 : * to have been a 2s complement number that was incorrectly decoded.
1130 : * */
1131 7599 : if (curIndex == 18) { /* forecast time is stored in octet 18-22 */
1132 454 : if (gfld->ipdtmpl[i] < -1 * (0x3fffffff)) {
1133 : /* Undoing the incorrect decoding of the negative number. */
1134 0 : is4[curIndex] = -1 * (int)(((unsigned)is4[curIndex])^(0x80000000));
1135 : }
1136 : }
1137 7599 : curIndex += abs(templatespds[pdsIndex].mappds[i]);
1138 : }
1139 : }
1140 :
1141 460 : is5[4] = 5;
1142 460 : is5[5] = gfld->ndpts;
1143 460 : is5[9] = gfld->idrtnum;
1144 460 : drsIndex = getdrsindex(gfld->idrtnum);
1145 460 : if (drsIndex == -1) {
1146 0 : jer[8 + *ndjer] = 2;
1147 0 : jer[8] = 2005; /* undefined sect 5 template */
1148 0 : *kjer = 9;
1149 0 : g2_free(gfld);
1150 0 : return;
1151 : }
1152 460 : curIndex = 11;
1153 3397 : for (i = 0; i < gfld->idrtlen; i++) {
1154 2937 : const struct drstemplate *templatesdrs = get_templatesdrs();
1155 2937 : is5[curIndex] = gfld->idrtmpl[i];
1156 2937 : curIndex += abs(templatesdrs[drsIndex].mapdrs[i]);
1157 : }
1158 : /* Mimic MDL's handling of reference value (IS5(12)) */
1159 460 : memcpy(&f_temp, &(is5[11]), sizeof (float));
1160 460 : is5[11] = FloatToSInt4Clamp(f_temp);
1161 460 : if ((is5[9] == 2) || (is5[9] == 3)) {
1162 87 : if (is5[20] == 0) {
1163 51 : memcpy(&(f_temp), &(is5[23]), sizeof (float));
1164 51 : *xmissp = f_temp;
1165 51 : is5[23] = FloatToSInt4Clamp(f_temp);
1166 51 : memcpy(&(f_temp), &(is5[27]), sizeof (float));
1167 51 : *xmisss = f_temp;
1168 51 : is5[27] = FloatToSInt4Clamp(f_temp);
1169 : } else {
1170 36 : *xmissp = is5[23];
1171 36 : *xmisss = is5[27];
1172 : }
1173 : }
1174 :
1175 460 : is6[4] = 6;
1176 460 : is6[5] = gfld->ibmap;
1177 460 : is7[4] = 7;
1178 :
1179 460 : if (subgNum + 1 == (unsigned int)numfields) {
1180 451 : *iendpk = 1;
1181 : } else {
1182 9 : *iendpk = 0;
1183 : }
1184 :
1185 460 : if ((gfld->ibmap == 0) || (gfld->ibmap == 254)) {
1186 9 : *ibitmap = 1;
1187 : } else {
1188 451 : *ibitmap = 0;
1189 : }
1190 :
1191 : /* Check type of original field, before transferring the memory. */
1192 460 : myAssert(*ns5 > 20);
1193 :
1194 460 : if( ain == NULL && iain == NULL )
1195 : {
1196 : /* Simulate effect of TransferInt / TransferFloat on the scan flag */
1197 303 : if( scanIndex >= 0 && (is3[scanIndex] & 0xf0) != GRIB2BIT_2 )
1198 : {
1199 16 : is3[scanIndex] = GRIB2BIT_2 + (is3[scanIndex] & 0x0f);
1200 : }
1201 :
1202 303 : g2_free(gfld);
1203 303 : return;
1204 : }
1205 :
1206 : /* Check if NCEP had problems expanding the data. If so we currently
1207 : * abort. May need to revisit this behavior. */
1208 157 : if (!gfld->expanded) {
1209 0 : jer[0 + *ndjer] = 2;
1210 0 : *kjer = 1;
1211 0 : g2_free(gfld);
1212 0 : return;
1213 : }
1214 157 : f_ignoreScan = 0;
1215 : /* Check if integer type... is5[20] == 1 implies integer (code table 5.1),
1216 : * but only for certain templates. (0,1,2,3,40,41,40000,40010). (not 4,50,51)
1217 : */
1218 157 : if ((is5[20] == 1) && ((is5[9] != 4) && (is5[9] != 50) && (is5[9] != 51))) {
1219 : /* integer data, use iain */
1220 67 : if ((scanIndex < 0) || (nxIndex < 0) || (nyIndex < 0)) {
1221 1 : ierr = TransferInt(gfld->fld, gfld->ngrdpts, *ibitmap, gfld->bmap,
1222 : 1, &(dummyScan), 0, 0, *iclean, *xmissp,
1223 : iain, *nd2x3, ib);
1224 : } else {
1225 66 : ierr = TransferInt(gfld->fld, gfld->ngrdpts, *ibitmap, gfld->bmap,
1226 66 : f_ignoreScan, &(is3[scanIndex]), is3[nxIndex],
1227 66 : is3[nyIndex], *iclean, *xmissp, iain, *nd2x3, ib);
1228 : }
1229 : } else {
1230 : /* float data, use ain */
1231 90 : if ((scanIndex < 0) || (nxIndex < 0) || (nyIndex < 0)) {
1232 0 : ierr = TransferFloat(gfld->fld, gfld->ngrdpts, *ibitmap,
1233 0 : gfld->bmap, 1, &(dummyScan), 0, 0, *iclean,
1234 : *xmissp, ain, *nd2x3, ib);
1235 : } else {
1236 90 : ierr = TransferFloat(gfld->fld, gfld->ngrdpts, *ibitmap,
1237 90 : gfld->bmap, f_ignoreScan, &(is3[scanIndex]),
1238 90 : is3[nxIndex], is3[nyIndex], *iclean, *xmissp,
1239 : ain, *nd2x3, ib);
1240 : }
1241 : }
1242 157 : if (ierr != 0) {
1243 0 : switch (ierr) {
1244 0 : case 1: /* nd2x3 is too small */
1245 0 : jer[0 + *ndjer] = 2;
1246 0 : *kjer = 1;
1247 0 : break;
1248 0 : case 2: /* nx * ny != ngrdpts */
1249 0 : jer[5 + *ndjer] = 2;
1250 0 : *kjer = 6;
1251 0 : break;
1252 0 : default:
1253 0 : jer[8 + *ndjer] = 2;
1254 0 : jer[8] = 9999; /* Unknown error message. */
1255 0 : *kjer = 9;
1256 : }
1257 0 : g2_free(gfld);
1258 0 : return;
1259 : }
1260 157 : g2_free(gfld);
1261 : }
1262 :
1263 : /*****************************************************************************
1264 : * BigByteCpy() --
1265 : *
1266 : * Arthur Taylor / MDL
1267 : *
1268 : * PURPOSE
1269 : * This is so we can copy up to 4 bytes from a big endian 4 byte int data
1270 : * stream.
1271 : *
1272 : * The reason this is needed is because the GRIB2 API required the GRIB2
1273 : * message to be passed in as a big endian 4 byte int data stream instead of
1274 : * something more reasonable (such as a stream of 1 byte char's) By having
1275 : * this routine we reduce the number of memswaps of the message on a little
1276 : * endian system.
1277 : *
1278 : * ARGUMENTS
1279 : * dst = Where to copy the data to. (Output)
1280 : * ipack = The 4 byte int data stream. (Input)
1281 : * nd5 = length of ipack (Input)
1282 : * startInt = Which integer to start reading in ipack. (Input)
1283 : * startByte = Which byte in the startInt to start reading from. (Input)
1284 : * numByte = How many bytes to read. (Input)
1285 : *
1286 : * FILES/DATABASES: None
1287 : *
1288 : * RETURNS: NULL
1289 : *
1290 : * HISTORY
1291 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
1292 : *
1293 : * NOTES
1294 : *****************************************************************************
1295 : */
1296 : #if 0
1297 : static void BigByteCpy(sInt4 *dst, sInt4 *ipack, sInt4 nd5,
1298 : unsigned int startInt, unsigned int startByte,
1299 : int numByte)
1300 : {
1301 : static char Lshift[] = { 0, 8, 16, 24 }; /* Amounts to shift left by. */
1302 : unsigned int intIndex; /* Where in ipack to read from. */
1303 : unsigned int byteIndex; /* Where in intIndex to read from. */
1304 : uInt4 curInt; /* An unsigned version of the current int. */
1305 : unsigned int curByte; /* The current byte we have read. */
1306 : int i; /* Loop counter over number of bytes to read. */
1307 :
1308 : myAssert(numByte <= 4);
1309 : myAssert(startByte < 4);
1310 : *dst = 0;
1311 : intIndex = startInt;
1312 : byteIndex = startByte;
1313 : for (i = 0; i < numByte; i++) {
1314 : if( intIndex >= (unsigned)nd5 )
1315 : {
1316 : /* TODO should error out */
1317 : return;
1318 : }
1319 : curInt = (uInt4)ipack[intIndex];
1320 : curByte = (curInt << Lshift[byteIndex]) >> 24;
1321 : *dst = (sInt4)((unsigned)*dst << 8) + curByte;
1322 : byteIndex++;
1323 : if (byteIndex == 4) {
1324 : byteIndex = 0;
1325 : intIndex++;
1326 : }
1327 : }
1328 : }
1329 : #endif
1330 :
1331 : /*****************************************************************************
1332 : * FindTemplateIDs() --
1333 : *
1334 : * Arthur Taylor / MDL
1335 : *
1336 : * PURPOSE
1337 : * This is so we can find out which templates are in this GRIB2 message.
1338 : * Which allows us to determine if we can call the MDL routines, or if we
1339 : * should call the NCEP routines instead.
1340 : *
1341 : * ARGUMENTS
1342 : * ipack = The 4 byte int data stream. (Input)
1343 : * nd5 = length of ipack (Input)
1344 : * subgNum = Which subgrid to look at. (Input)
1345 : * gdsTmpl = The gds template number for this subgrid. (Output)
1346 : * pdsTmpl = The pds template number for this subgrid. (Output)
1347 : * drsTmpl = The drs template number for this subgrid. (Output)
1348 : * f_noBitmap = 0 has a bitmap, else doesn't have a bitmap. (Output)
1349 : * orderDiff = The order of the differencing in template 5.3 (1st, 2nd) (out)
1350 : *
1351 : * FILES/DATABASES: None
1352 : *
1353 : * RETURNS: int
1354 : * 0 = Ok.
1355 : * 2 = invalid section number
1356 : *
1357 : * HISTORY
1358 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
1359 : *
1360 : * NOTES
1361 : *****************************************************************************
1362 : */
1363 : #if 0
1364 : static int FindTemplateIDs(sInt4 *ipack, sInt4 nd5, int subgNum,
1365 : sInt4 *gdsTmpl, sInt4 *pdsTmpl,
1366 : sInt4 *drsTmpl, sInt4 *numGrps,
1367 : uChar *f_noBitmap, sInt4 *orderDiff)
1368 : {
1369 : unsigned int gNum = 0; /* Which sub group we are currently working with. */
1370 : sInt4 offset; /* Where in the bytes stream we currently are. */
1371 : sInt4 sectLen; /* The length of the current section. */
1372 : sInt4 sectId; /* The current section number. */
1373 : sInt4 li_temp; /* A temporary holder for the bitmap flag. */
1374 :
1375 : /* Jump over section 0. */
1376 : offset = 16;
1377 : while (gNum <= (unsigned int)subgNum) {
1378 : BigByteCpy(§Len, ipack, nd5, (offset / 4), (offset % 4), 4);
1379 : /* Check if we just read section 8. If so, then it is "7777" =
1380 : * 926365495 regardless of endian'ness. */
1381 : if (sectLen == 926365495L) {
1382 : #ifdef DEBUG
1383 : printf("Shouldn't see sect 8. Should stop after correct sect 5\n");
1384 : #endif
1385 : return 2;
1386 : }
1387 : BigByteCpy(§Id, ipack, nd5, ((offset + 4) / 4),
1388 : ((offset + 4) % 4), 1);
1389 : switch (sectId) {
1390 : case 1:
1391 : case 2:
1392 : break;
1393 : case 3:
1394 : BigByteCpy(gdsTmpl, ipack, nd5, ((offset + 12) / 4),
1395 : ((offset + 12) % 4), 2);
1396 : break;
1397 : case 4:
1398 : BigByteCpy(pdsTmpl, ipack, nd5, ((offset + 7) / 4),
1399 : ((offset + 7) % 4), 2);
1400 : break;
1401 : case 5:
1402 : BigByteCpy(drsTmpl, ipack, nd5, ((offset + 9) / 4),
1403 : ((offset + 9) % 4), 2);
1404 : if ((*drsTmpl == 2) || (*drsTmpl == 3)) {
1405 : BigByteCpy(numGrps, ipack, nd5, ((offset + 31) / 4),
1406 : ((offset + 31) % 4), 4);
1407 : } else {
1408 : *numGrps = 0;
1409 : }
1410 : if (*drsTmpl == 3) {
1411 : BigByteCpy(&li_temp, ipack, nd5, ((offset + 44) / 4),
1412 : ((offset + 44) % 4), 1);
1413 : *orderDiff = li_temp;
1414 : } else {
1415 : *orderDiff = 0;
1416 : }
1417 : break;
1418 : case 6:
1419 : BigByteCpy(&li_temp, ipack, nd5, ((offset + 5) / 4),
1420 : ((offset + 5) % 4), 1);
1421 : if (li_temp == 255) {
1422 : *f_noBitmap = 1;
1423 : }
1424 : gNum++;
1425 : break;
1426 : case 7:
1427 : break;
1428 : default:
1429 : #ifdef DEBUG
1430 : printf("Invalid section id %d.\n", sectId);
1431 : #endif
1432 : return 2;
1433 : }
1434 : offset += sectLen;
1435 : }
1436 : return 0;
1437 : }
1438 : #endif
1439 :
1440 : /*****************************************************************************
1441 : * unpk_grib2() --
1442 : *
1443 : * Arthur Taylor / MDL
1444 : *
1445 : * PURPOSE
1446 : * This procedure is the main API for decoding GRIB2 messages. It is
1447 : * intended to be a branching routine to call either the MDL GRIB2 library,
1448 : * or the NCEP GRIB2 library, depending on which template it sees in the
1449 : * message.
1450 : *
1451 : * ARGUMENTS
1452 : * kfildo = Unit number for output diagnostics (C ignores this). (Input)
1453 : * ain = contains the data if the original data was float data.
1454 : * (size = nd2x3) (Output)
1455 : * iain = contains the data if the original data was integer data.
1456 : * (size = nd2x3) (Output)
1457 : * nd2x3 = length of ain, iain, ib (Input)
1458 : * (at least the size of num grid points)
1459 : * idat = local use data if any that were unpacked from Section 2. (Output)
1460 : * nidat = length of idat. (Input)
1461 : * rdat = floating point local use data (Output)
1462 : * nrdat = length of rdat. (Input)
1463 : * is0 = Section 0 data (Output)
1464 : * ns0 = length of is0 (16 is fine) (Input)
1465 : * is1 = Section 1 data (Output)
1466 : * ns1 = length of is1 (21 is fine) (Input)
1467 : * is2 = Section 2 data (Output)
1468 : * ns2 = length of is2 () (Input)
1469 : * is3 = Section 3 data (Output)
1470 : * ns3 = length of is3 (96 or 1600) (Input)
1471 : * is4 = Section 4 data (Output)
1472 : * ns4 = length of is4 (60) (Input)
1473 : * is5 = Section 5 data (Output)
1474 : * ns5 = length of is5 (49 is fine) (Input)
1475 : * is6 = Section 6 data (Output)
1476 : * ns6 = length of is6 (6 is fine) (Input)
1477 : * is7 = Section 7 data (Output)
1478 : * ns7 = length of is7 (8 is fine) (Input)
1479 : * ib = Bitmap if user requested it, and it was packed (Output)
1480 : * ibitmap = 0 means ib is invalid, 1 means ib is valid. (Output)
1481 : * ipack = The message to unpack (This is assumed to be Big endian) (Input)
1482 : * nd5 = The size of ipack (Input)
1483 : * xmissp = The floating point representation for the primary missing value.
1484 : * (Input/Output)
1485 : * xmisss = The floating point representation for the secondary missing
1486 : * value (Input/Output)
1487 : * new = 1 if this is the first grid to be unpacked, else 0. (Input)
1488 : * iclean = 1 means the user wants the unpacked data returned without
1489 : * missing values in it. 0 means embed the missing values. (Input)
1490 : * l3264b = Integer word length in bits (32 or 64) (Input)
1491 : * iendpk = A 1 means no more grids in this message, a 0 means more grids.
1492 : * (Output)
1493 : * jer(ndjer,2) = error codes along with severity. (Output)
1494 : * ndjer = 1/2 length of jer. (>= 15) (Input)
1495 : * kjer = number of error messages stored in jer. (Output)
1496 : *
1497 : * FILES/DATABASES: None
1498 : *
1499 : * RETURNS: void
1500 : *
1501 : * HISTORY
1502 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
1503 : *
1504 : * NOTES
1505 : * Inefficiencies: have to memswap ipack multiple times.
1506 : * MDL handles is5[12], is5[23], and is5[27] in an "interesting" manner.
1507 : * MDL attempts to always return grids in scan mode 0100????
1508 : * ToDo: Check length of parameters better.
1509 : *
1510 : * According to MDL's unpk_grib2.f, it currently supports (so for others we
1511 : * call the NCEP routines):
1512 : * TEMPLATE 3.0 EQUIDISTANT CYLINDRICAL LATITUDE/LONGITUDE
1513 : * TEMPLATE 3.10 MERCATOR
1514 : * TEMPLATE 3.20 POLAR STEREOGRAPHIC
1515 : * TEMPLATE 3.30 LAMBERT
1516 : * TEMPLATE 3.90 ORTHOGRAPHIC SPACE VIEW
1517 : * TEMPLATE 3.110 EQUATORIAL AZIMUTHAL EQUIDISTANT
1518 : * TEMPLATE 3.120 AZIMUTH-RANGE (RADAR)
1519 : *
1520 : * TEMPLATE 4.0 ANALYSIS OR FORECAST AT A LEVEL AND POINT
1521 : * TEMPLATE 4.1 INDIVIDUAL ENSEMBLE
1522 : * TEMPLATE 4.2 DERIVED FORECAST BASED ON ENSEMBLES
1523 : * TEMPLATE 4.8 AVERAGE, ACCUMULATION, EXTREMES
1524 : * TEMPLATE 4.20 RADAR
1525 : * TEMPLATE 4.30 SATELLITE
1526 : *
1527 : * TEMPLATE 5.0 SIMPLE PACKING
1528 : * TEMPLATE 5.2 COMPLEX PACKING
1529 : * TEMPLATE 5.3 COMPLEX PACKING AND SPATIAL DIFFERENCING
1530 : *
1531 : * Correction to "unpk_grib2.f" : It also supports:
1532 : * TEMPLATE 4.9 Probability forecast in a time interval
1533 : *
1534 : *****************************************************************************
1535 : */
1536 : #ifdef unused_by_GDAL
1537 : void unpk_grib2(sInt4 *kfildo, float *ain, sInt4 *iain, sInt4 *nd2x3,
1538 : sInt4 *idat, sInt4 *nidat, float *rdat, sInt4 *nrdat,
1539 : sInt4 *is0, sInt4 *ns0, sInt4 *is1, sInt4 *ns1,
1540 : sInt4 *is2, sInt4 *ns2, sInt4 *is3, sInt4 *ns3,
1541 : sInt4 *is4, sInt4 *ns4, sInt4 *is5, sInt4 *ns5,
1542 : sInt4 *is6, sInt4 *ns6, sInt4 *is7, sInt4 *ns7,
1543 : sInt4 *ib, sInt4 *ibitmap, sInt4 *ipack, sInt4 *nd5,
1544 : float *xmissp, float *xmisss, sInt4 *inew,
1545 : sInt4 *iclean, sInt4 *l3264b, sInt4 *iendpk, sInt4 *jer,
1546 : sInt4 *ndjer, sInt4 *kjer)
1547 : {
1548 : unsigned char *c_ipack; /* The compressed data as char instead of sInt4 so
1549 : * it is easier to work with. */
1550 : #if 0
1551 : char f_useMDL = 0; /* Instructed 3/8/2005 10:30 to not use MDL. */
1552 : #endif
1553 :
1554 : /* Since NCEP's code works with byte level stuff, we need to un-do the
1555 : * byte swap of the (int *) data, then cast it to an (unsigned char *) */
1556 : #ifndef WORDS_BIGENDIAN
1557 : /* Can't make this dependent on inew, since they could have a sequence of
1558 : * get first message... do stuff, get second message, which unfortunately
1559 : * means they would have to get the first message again, causing 2 swaps,
1560 : * and breaking their second request for the first message (as well as
1561 : * their second message). */
1562 : memswp(ipack, sizeof(sInt4), *nd5);
1563 : #endif
1564 : c_ipack = (unsigned char *)ipack;
1565 : unpk_g2ncep(kfildo, ain, iain, nd2x3, idat, nidat, rdat, nrdat, is0,
1566 : ns0, is1, ns1, is2, ns2, is3, ns3, is4, ns4, is5, ns5,
1567 : is6, ns6, is7, ns7, ib, ibitmap, c_ipack, nd5, xmissp,
1568 : xmisss, inew, iclean, l3264b, iendpk, jer, ndjer, kjer);
1569 :
1570 : #ifndef WORDS_BIGENDIAN
1571 : /* Swap back because we could be called again for the subgrid data. */
1572 : memswp(ipack, sizeof(sInt4), *nd5);
1573 : #endif
1574 : }
1575 : #endif
1576 :
1577 : /*****************************************************************************
1578 : * C_pkGrib2() --
1579 : *
1580 : * Arthur Taylor / MDL
1581 : *
1582 : * PURPOSE
1583 : * This procedure is the main API for encoding GRIB2 messages. It is
1584 : * intended to call NCEP's GRIB2 library.
1585 : *
1586 : * RETURNS: void
1587 : *
1588 : * HISTORY
1589 : * 1/2004 Arthur Taylor (MDL/RSIS): Created.
1590 : *
1591 : * NOTES
1592 : *****************************************************************************
1593 : */
1594 :
1595 : #ifdef unused_by_GDAL
1596 : int C_pkGrib2(unsigned char *cgrib, sInt4 *sec0, sInt4 *sec1,
1597 : unsigned char *csec2, sInt4 lcsec2,
1598 : sInt4 *igds, sInt4 *igdstmpl, sInt4 *ideflist,
1599 : sInt4 idefnum, sInt4 ipdsnum, sInt4 *ipdstmpl,
1600 : float *coordlist, sInt4 numcoord, sInt4 idrsnum,
1601 : sInt4 *idrstmpl, float *fld, sInt4 ngrdpts,
1602 : sInt4 ibmap, sInt4 *bmap)
1603 : {
1604 : int ierr; /* error value from grib2 library. */
1605 :
1606 : if ((ierr = g2_create(cgrib, sec0, sec1)) == -1) {
1607 : /* Tried to use for version other than GRIB Ed 2 */
1608 : return -1;
1609 : }
1610 :
1611 : if ((ierr = g2_addlocal(cgrib, csec2, lcsec2)) < 0) {
1612 : /* Some how got a bad section 2. Should be impossible unless an assert
1613 : * was broken. */
1614 : return -2;
1615 : }
1616 :
1617 : if ((ierr = g2_addgrid(cgrib, igds, igdstmpl, ideflist, idefnum)) < 0) {
1618 : /* Some how got a bad section 3. Only way would be should be with an
1619 : * unsupported template number unless an assert was broken. */
1620 : return -3;
1621 : }
1622 :
1623 : if ((ierr = g2_addfield(cgrib, ipdsnum, ipdstmpl, coordlist, numcoord,
1624 : idrsnum, idrstmpl, fld, ngrdpts, ibmap,
1625 : bmap)) < 0) {
1626 : return -4;
1627 : }
1628 :
1629 : if ((ierr = g2_gribend(cgrib)) < 0) {
1630 : return -5;
1631 : }
1632 :
1633 : return ierr;
1634 : }
1635 : #endif /* unused_by_GDAL */
|
