1 : #include <stdio.h>
2 : #include <stdlib.h>
3 : #include "grib2.h"
4 :
5 : g2int g2_unpack1(unsigned char *,g2int *,g2int **,g2int *);
6 : g2int g2_unpack2(unsigned char *,g2int *,g2int *,unsigned char **);
7 : g2int g2_unpack3(unsigned char *,g2int *,g2int **,g2int **,
8 : g2int *,g2int **,g2int *);
9 : g2int g2_unpack4(unsigned char *,g2int *,g2int *,g2int **,
10 : g2int *,g2float **,g2int *);
11 : g2int g2_unpack5(unsigned char *,g2int *,g2int *,g2int *, g2int **,g2int *);
12 : g2int g2_unpack6(unsigned char *,g2int *,g2int ,g2int *, g2int **);
13 : g2int g2_unpack7(unsigned char *,g2int *,g2int ,g2int *,
14 : g2int ,g2int *,g2int ,g2float **);
15 :
16 2 : g2int g2_getfld(unsigned char *cgrib,g2int ifldnum,g2int unpack,g2int expand,
17 : gribfield **gfld)
18 : //$$$ SUBPROGRAM DOCUMENTATION BLOCK
19 : // . . . .
20 : // SUBPROGRAM: g2_getfld
21 : // PRGMMR: Gilbert ORG: W/NP11 DATE: 2002-10-28
22 : //
23 : // ABSTRACT: This subroutine returns all the metadata, template values,
24 : // Bit-map ( if applicable ), and the unpacked data for a given data
25 : // field. All of the information returned is stored in a gribfield
26 : // structure, which is defined in file grib2.h.
27 : // Users of this routine will need to include "grib2.h" in their source
28 : // code that calls this routine. Each component of the gribfield
29 : // struct is also described in the OUTPUT ARGUMENTS section below.
30 : //
31 : // Since there can be multiple data fields packed into a GRIB2
32 : // message, the calling routine indicates which field is being requested
33 : // with the ifldnum argument.
34 : //
35 : // PROGRAM HISTORY LOG:
36 : // 2002-10-28 Gilbert
37 : //
38 : // USAGE: #include "grib2.h"
39 : // int g2_getfld(unsigned char *cgrib,g2int ifldnum,g2int unpack,
40 : // g2int expand,gribfield **gfld)
41 : // INPUT ARGUMENTS:
42 : // cgrib - Character pointer to the GRIB2 message
43 : // ifldnum - Specifies which field in the GRIB2 message to return.
44 : // unpack - Boolean value indicating whether to unpack bitmap/data field
45 : // 1 = unpack bitmap (if present) and data values
46 : // 0 = do not unpack bitmap and data values
47 : // expand - Boolean value indicating whether the data points should be
48 : // expanded to the correspond grid, if a bit-map is present.
49 : // 1 = if possible, expand data field to grid, inserting zero
50 : // values at gridpoints that are bitmapped out.
51 : // (SEE REMARKS2)
52 : // 0 = do not expand data field, leaving it an array of
53 : // consecutive data points for each "1" in the bitmap.
54 : // This argument is ignored if unpack == 0 OR if the
55 : // returned field does not contain a bit-map.
56 : //
57 : // OUTPUT ARGUMENT:
58 : // gribfield gfld; - pointer to structure gribfield containing
59 : // all decoded data for the data field.
60 : //
61 : // gfld->version = GRIB edition number ( currently 2 )
62 : // gfld->discipline = Message Discipline ( see Code Table 0.0 )
63 : // gfld->idsect = Contains the entries in the Identification
64 : // Section ( Section 1 )
65 : // This element is a pointer to an array
66 : // that holds the data.
67 : // gfld->idsect[0] = Identification of originating Centre
68 : // ( see Common Code Table C-1 )
69 : // 7 - US National Weather Service
70 : // gfld->idsect[1] = Identification of originating Sub-centre
71 : // gfld->idsect[2] = GRIB Master Tables Version Number
72 : // ( see Code Table 1.0 )
73 : // 0 - Experimental
74 : // 1 - Initial operational version number
75 : // gfld->idsect[3] = GRIB Local Tables Version Number
76 : // ( see Code Table 1.1 )
77 : // 0 - Local tables not used
78 : // 1-254 - Number of local tables version used
79 : // gfld->idsect[4] = Significance of Reference Time (Code Table 1.2)
80 : // 0 - Analysis
81 : // 1 - Start of forecast
82 : // 2 - Verifying time of forecast
83 : // 3 - Observation time
84 : // gfld->idsect[5] = Year ( 4 digits )
85 : // gfld->idsect[6] = Month
86 : // gfld->idsect[7) = Day
87 : // gfld->idsect[8] = Hour
88 : // gfld->idsect[9] = Minute
89 : // gfld->idsect[10] = Second
90 : // gfld->idsect[11] = Production status of processed data
91 : // ( see Code Table 1.3 )
92 : // 0 - Operational products
93 : // 1 - Operational test products
94 : // 2 - Research products
95 : // 3 - Re-analysis products
96 : // gfld->idsect[12] = Type of processed data ( see Code Table 1.4 )
97 : // 0 - Analysis products
98 : // 1 - Forecast products
99 : // 2 - Analysis and forecast products
100 : // 3 - Control forecast products
101 : // 4 - Perturbed forecast products
102 : // 5 - Control and perturbed forecast products
103 : // 6 - Processed satellite observations
104 : // 7 - Processed radar observations
105 : // gfld->idsectlen = Number of elements in gfld->idsect[].
106 : // gfld->local = Pointer to character array containing contents
107 : // of Local Section 2, if included
108 : // gfld->locallen = length of array gfld->local[]
109 : // gfld->ifldnum = field number within GRIB message
110 : // gfld->griddef = Source of grid definition (see Code Table 3.0)
111 : // 0 - Specified in Code table 3.1
112 : // 1 - Predetermined grid Defined by originating centre
113 : // gfld->ngrdpts = Number of grid points in the defined grid.
114 : // gfld->numoct_opt = Number of octets needed for each
115 : // additional grid points definition.
116 : // Used to define number of
117 : // points in each row ( or column ) for
118 : // non-regular grids.
119 : // = 0, if using regular grid.
120 : // gfld->interp_opt = Interpretation of list for optional points
121 : // definition. (Code Table 3.11)
122 : // gfld->igdtnum = Grid Definition Template Number (Code Table 3.1)
123 : // gfld->igdtmpl = Contains the data values for the specified Grid
124 : // Definition Template ( NN=gfld->igdtnum ). Each
125 : // element of this integer array contains an entry (in
126 : // the order specified) of Grid Defintion Template 3.NN
127 : // This element is a pointer to an array
128 : // that holds the data.
129 : // gfld->igdtlen = Number of elements in gfld->igdtmpl[]. i.e. number of
130 : // entries in Grid Defintion Template 3.NN
131 : // ( NN=gfld->igdtnum ).
132 : // gfld->list_opt = (Used if gfld->numoct_opt .ne. 0) This array
133 : // contains the number of grid points contained in
134 : // each row ( or column ). (part of Section 3)
135 : // This element is a pointer to an array
136 : // that holds the data. This pointer is nullified
137 : // if gfld->numoct_opt=0.
138 : // gfld->num_opt = (Used if gfld->numoct_opt .ne. 0)
139 : // The number of entries
140 : // in array ideflist. i.e. number of rows ( or columns )
141 : // for which optional grid points are defined. This value
142 : // is set to zero, if gfld->numoct_opt=0.
143 : // gfdl->ipdtnum = Product Definition Template Number(see Code Table 4.0)
144 : // gfld->ipdtmpl = Contains the data values for the specified Product
145 : // Definition Template ( N=gfdl->ipdtnum ). Each element
146 : // of this integer array contains an entry (in the
147 : // order specified) of Product Defintion Template 4.N.
148 : // This element is a pointer to an array
149 : // that holds the data.
150 : // gfld->ipdtlen = Number of elements in gfld->ipdtmpl[]. i.e. number of
151 : // entries in Product Defintion Template 4.N
152 : // ( N=gfdl->ipdtnum ).
153 : // gfld->coord_list = Real array containing floating point values
154 : // intended to document the vertical discretisation
155 : // associated to model data on hybrid coordinate
156 : // vertical levels. (part of Section 4)
157 : // This element is a pointer to an array
158 : // that holds the data.
159 : // gfld->num_coord = number of values in array gfld->coord_list[].
160 : // gfld->ndpts = Number of data points unpacked and returned.
161 : // gfld->idrtnum = Data Representation Template Number
162 : // ( see Code Table 5.0)
163 : // gfld->idrtmpl = Contains the data values for the specified Data
164 : // Representation Template ( N=gfld->idrtnum ). Each
165 : // element of this integer array contains an entry
166 : // (in the order specified) of Product Defintion
167 : // Template 5.N.
168 : // This element is a pointer to an array
169 : // that holds the data.
170 : // gfld->idrtlen = Number of elements in gfld->idrtmpl[]. i.e. number
171 : // of entries in Data Representation Template 5.N
172 : // ( N=gfld->idrtnum ).
173 : // gfld->unpacked = logical value indicating whether the bitmap and
174 : // data values were unpacked. If false,
175 : // gfld->bmap and gfld->fld pointers are nullified.
176 : // gfld->expanded = Logical value indicating whether the data field
177 : // was expanded to the grid in the case where a
178 : // bit-map is present. If true, the data points in
179 : // gfld->fld match the grid points and zeros were
180 : // inserted at grid points where data was bit-mapped
181 : // out. If false, the data values in gfld->fld were
182 : // not expanded to the grid and are just a consecutive
183 : // array of data points corresponding to each value of
184 : // "1" in gfld->bmap.
185 : // gfld->ibmap = Bitmap indicator ( see Code Table 6.0 )
186 : // 0 = bitmap applies and is included in Section 6.
187 : // 1-253 = Predefined bitmap applies
188 : // 254 = Previously defined bitmap applies to this field
189 : // 255 = Bit map does not apply to this product.
190 : // gfld->bmap = integer array containing decoded bitmap,
191 : // if gfld->ibmap=0 or gfld->ibap=254. Otherwise nullified
192 : // This element is a pointer to an array
193 : // that holds the data.
194 : // gfld->fld = Array of gfld->ndpts unpacked data points.
195 : // This element is a pointer to an array
196 : // that holds the data.
197 : //
198 : //
199 : // RETURN VALUES:
200 : // ierr - Error return code.
201 : // 0 = no error
202 : // 1 = Beginning characters "GRIB" not found.
203 : // 2 = GRIB message is not Edition 2.
204 : // 3 = The data field request number was not positive.
205 : // 4 = End string "7777" found, but not where expected.
206 : // 6 = GRIB message did not contain the requested number of
207 : // data fields.
208 : // 7 = End string "7777" not found at end of message.
209 : // 8 = Unrecognized Section encountered.
210 : // 9 = Data Representation Template 5.NN not yet implemented.
211 : // 15 = Error unpacking Section 1.
212 : // 16 = Error unpacking Section 2.
213 : // 10 = Error unpacking Section 3.
214 : // 11 = Error unpacking Section 4.
215 : // 12 = Error unpacking Section 5.
216 : // 13 = Error unpacking Section 6.
217 : // 14 = Error unpacking Section 7.
218 : // 17 = Previous bitmap specified, yet none exists.
219 : //
220 : // REMARKS: Note that struct gribfield is allocated by this routine and it
221 : // also contains pointers to many arrays of data that were allocated
222 : // during decoding. Users are encouraged to free up this memory,
223 : // when it is no longer needed, by an explicit call to routine g2_free.
224 : // EXAMPLE:
225 : // #include "grib2.h"
226 : // gribfield *gfld;
227 : // ret=g2_getfld(cgrib,1,1,1,&gfld);
228 : // ...
229 : // g2_free(gfld);
230 : //
231 : // Routine g2_info can be used to first determine
232 : // how many data fields exist in a given GRIB message.
233 : //
234 : // REMARKS2: It may not always be possible to expand a bit-mapped data field.
235 : // If a pre-defined bit-map is used and not included in the GRIB2
236 : // message itself, this routine would not have the necessary
237 : // information to expand the data. In this case, gfld->expanded would
238 : // would be set to 0 (false), regardless of the value of input
239 : // argument expand.
240 : //
241 : // ATTRIBUTES:
242 : // LANGUAGE: C
243 : // MACHINE:
244 : //
245 : //$$$
246 : {
247 :
248 : g2int have3,have4,have5,have6,have7,ierr,jerr;
249 : g2int numfld,j,n,istart,iofst,ipos;
250 : g2int disc,ver,lensec0,lengrib,lensec,isecnum;
251 : g2int *igds;
252 : g2int *bmpsave;
253 : g2float *newfld;
254 : gribfield *lgfld;
255 :
256 2 : have3=0;
257 2 : have4=0;
258 2 : have5=0;
259 2 : have6=0;
260 2 : have7=0;
261 2 : ierr=0;
262 2 : numfld=0;
263 :
264 2 : lgfld=(gribfield *)malloc(sizeof(gribfield));
265 2 : *gfld=lgfld;
266 :
267 2 : lgfld->locallen=0;
268 2 : lgfld->idsect=0;
269 2 : lgfld->local=0;
270 2 : lgfld->list_opt=0;
271 2 : lgfld->igdtmpl=0;
272 2 : lgfld->ipdtmpl=0;
273 2 : lgfld->idrtmpl=0;
274 2 : lgfld->coord_list=0;
275 2 : lgfld->bmap=0;
276 2 : lgfld->fld=0;
277 : //
278 : // Check for valid request number
279 : //
280 2 : if (ifldnum <= 0) {
281 0 : printf("g2_getfld: Request for field number must be positive.\n");
282 0 : ierr=3;
283 0 : return(ierr);
284 : }
285 : //
286 : // Check for beginning of GRIB message in the first 100 bytes
287 : //
288 2 : istart=-1;
289 2 : for (j=0;j<100;j++) {
290 4 : if (cgrib[j]=='G' && cgrib[j+1]=='R' &&cgrib[j+2]=='I' &&
291 2 : cgrib[j+3]=='B') {
292 2 : istart=j;
293 2 : break;
294 : }
295 : }
296 2 : if (istart == -1) {
297 0 : printf("g2_getfld: Beginning characters GRIB not found.\n");
298 0 : ierr=1;
299 0 : return(ierr);
300 : }
301 : //
302 : // Unpack Section 0 - Indicator Section
303 : //
304 2 : iofst=8*(istart+6);
305 2 : gbit(cgrib,&disc,iofst,8); // Discipline
306 2 : iofst=iofst+8;
307 2 : gbit(cgrib,&ver,iofst,8); // GRIB edition number
308 2 : iofst=iofst+8;
309 2 : iofst=iofst+32;
310 2 : gbit(cgrib,&lengrib,iofst,32); // Length of GRIB message
311 2 : iofst=iofst+32;
312 2 : lensec0=16;
313 2 : ipos=istart+lensec0;
314 : //
315 : // Currently handles only GRIB Edition 2.
316 : //
317 2 : if (ver != 2) {
318 0 : printf("g2_getfld: can only decode GRIB edition 2.\n");
319 0 : ierr=2;
320 0 : return(ierr);
321 : }
322 : //
323 : // Loop through the remaining sections keeping track of the
324 : // length of each. Also keep the latest Grid Definition Section info.
325 : // Unpack the requested field number.
326 : //
327 : for (;;) {
328 : // Check to see if we are at end of GRIB message
329 12 : if (cgrib[ipos]=='7' && cgrib[ipos+1]=='7' && cgrib[ipos+2]=='7' &&
330 0 : cgrib[ipos+3]=='7') {
331 0 : ipos=ipos+4;
332 : // If end of GRIB message not where expected, issue error
333 0 : if (ipos != (istart+lengrib)) {
334 0 : printf("g2_getfld: '7777' found, but not where expected.\n");
335 0 : ierr=4;
336 0 : return(ierr);
337 : }
338 : break;
339 : }
340 : // Get length of Section and Section number
341 12 : iofst=(ipos-1)*8;
342 12 : iofst=ipos*8;
343 12 : gbit(cgrib,&lensec,iofst,32); // Get Length of Section
344 12 : iofst=iofst+32;
345 12 : gbit(cgrib,&isecnum,iofst,8); // Get Section number
346 12 : iofst=iofst+8;
347 : //printf(" lensec= %ld secnum= %ld \n",lensec,isecnum);
348 : //
349 : // Check to see if section number is valid
350 : //
351 12 : if ( isecnum<1 || isecnum>7 ) {
352 0 : printf("g2_getfld: Unrecognized Section Encountered=%d\n",isecnum);
353 0 : ierr=8;
354 0 : return(ierr);
355 : }
356 : //
357 : // If found Section 1, decode elements in Identification Section
358 : //
359 12 : if (isecnum == 1) {
360 2 : iofst=iofst-40; // reset offset to beginning of section
361 2 : jerr=g2_unpack1(cgrib,&iofst,&lgfld->idsect,&lgfld->idsectlen);
362 2 : if (jerr !=0 ) {
363 0 : ierr=15;
364 0 : return(ierr);
365 : }
366 : }
367 : //
368 : // If found Section 2, Grab local section
369 : // Save in case this is the latest one before the requested field.
370 : //
371 12 : if (isecnum == 2) {
372 0 : iofst=iofst-40; // reset offset to beginning of section
373 0 : if (lgfld->local!=0) free(lgfld->local);
374 0 : jerr=g2_unpack2(cgrib,&iofst,&lgfld->locallen,&lgfld->local);
375 0 : if (jerr != 0) {
376 0 : ierr=16;
377 0 : return(ierr);
378 : }
379 : }
380 : //
381 : // If found Section 3, unpack the GDS info using the
382 : // appropriate template. Save in case this is the latest
383 : // grid before the requested field.
384 : //
385 12 : if (isecnum == 3) {
386 2 : iofst=iofst-40; // reset offset to beginning of section
387 2 : if (lgfld->igdtmpl!=0) free(lgfld->igdtmpl);
388 2 : if (lgfld->list_opt!=0) free(lgfld->list_opt);
389 2 : jerr=g2_unpack3(cgrib,&iofst,&igds,&lgfld->igdtmpl,
390 : &lgfld->igdtlen,&lgfld->list_opt,&lgfld->num_opt);
391 2 : if (jerr == 0) {
392 2 : have3=1;
393 2 : lgfld->griddef=igds[0];
394 2 : lgfld->ngrdpts=igds[1];
395 2 : lgfld->numoct_opt=igds[2];
396 2 : lgfld->interp_opt=igds[3];
397 2 : lgfld->igdtnum=igds[4];
398 2 : free( igds );
399 : }
400 : else {
401 0 : ierr=10;
402 0 : return(ierr);
403 : }
404 : }
405 : //
406 : // If found Section 4, check to see if this field is the
407 : // one requested.
408 : //
409 12 : if (isecnum == 4) {
410 2 : numfld=numfld+1;
411 2 : if (numfld == ifldnum) {
412 2 : lgfld->discipline=disc;
413 2 : lgfld->version=ver;
414 2 : lgfld->ifldnum=ifldnum;
415 2 : lgfld->unpacked=unpack;
416 2 : lgfld->expanded=0;
417 2 : iofst=iofst-40; // reset offset to beginning of section
418 2 : jerr=g2_unpack4(cgrib,&iofst,&lgfld->ipdtnum,
419 : &lgfld->ipdtmpl,&lgfld->ipdtlen,&lgfld->coord_list,
420 : &lgfld->num_coord);
421 2 : if (jerr == 0)
422 2 : have4=1;
423 : else {
424 0 : ierr=11;
425 0 : return(ierr);
426 : }
427 : }
428 : }
429 : //
430 : // If found Section 5, check to see if this field is the
431 : // one requested.
432 : //
433 12 : if (isecnum == 5 && numfld == ifldnum) {
434 2 : iofst=iofst-40; // reset offset to beginning of section
435 2 : jerr=g2_unpack5(cgrib,&iofst,&lgfld->ndpts,&lgfld->idrtnum,
436 : &lgfld->idrtmpl,&lgfld->idrtlen);
437 2 : if (jerr == 0)
438 2 : have5=1;
439 : else {
440 0 : ierr=12;
441 0 : return(ierr);
442 : }
443 : }
444 : //
445 : // If found Section 6, Unpack bitmap.
446 : // Save in case this is the latest
447 : // bitmap before the requested field.
448 : //
449 12 : if (isecnum == 6) {
450 2 : if (unpack) { // unpack bitmap
451 2 : iofst=iofst-40; // reset offset to beginning of section
452 2 : bmpsave=lgfld->bmap; // save pointer to previous bitmap
453 2 : jerr=g2_unpack6(cgrib,&iofst,lgfld->ngrdpts,&lgfld->ibmap,
454 : &lgfld->bmap);
455 2 : if (jerr == 0) {
456 2 : have6=1;
457 2 : if (lgfld->ibmap == 254) // use previously specified bitmap
458 0 : if( bmpsave!=0 )
459 0 : lgfld->bmap=bmpsave;
460 : else {
461 0 : printf("g2_getfld: Prev bit-map specified, but none exist.\n");
462 0 : ierr=17;
463 0 : return(ierr);
464 : }
465 : else // get rid of it
466 2 : if( bmpsave!=0 ) free(bmpsave);
467 : }
468 : else {
469 0 : ierr=13;
470 0 : return(ierr);
471 : }
472 : }
473 : else { // do not unpack bitmap
474 0 : gbit(cgrib,&lgfld->ibmap,iofst,8); // Get BitMap Indicator
475 0 : have6=1;
476 : }
477 : }
478 : //
479 : // If found Section 7, check to see if this field is the
480 : // one requested.
481 : //
482 12 : if (isecnum==7 && numfld==ifldnum && unpack) {
483 2 : iofst=iofst-40; // reset offset to beginning of section
484 2 : jerr=g2_unpack7(cgrib,&iofst,lgfld->igdtnum,lgfld->igdtmpl,
485 : lgfld->idrtnum,lgfld->idrtmpl,lgfld->ndpts,
486 : &lgfld->fld);
487 2 : if (jerr == 0) {
488 2 : have7=1;
489 : // If bitmap is used with this field, expand data field
490 : // to grid, if possible.
491 2 : if ( lgfld->ibmap != 255 && lgfld->bmap != 0 ) {
492 0 : if ( expand == 1 ) {
493 0 : n=0;
494 0 : newfld=(g2float *)calloc(lgfld->ngrdpts,sizeof(g2float));
495 0 : for (j=0;j<lgfld->ngrdpts;j++) {
496 0 : if (lgfld->bmap[j]==1) newfld[j]=lgfld->fld[n++];
497 : }
498 0 : free(lgfld->fld);
499 0 : lgfld->fld=newfld;
500 0 : lgfld->expanded=1;
501 : }
502 : else {
503 0 : lgfld->expanded=0;
504 : }
505 : }
506 : else {
507 2 : lgfld->expanded=1;
508 : }
509 : }
510 : else {
511 0 : printf("g2_getfld: return from g2_unpack7 = %d \n",(int)jerr);
512 0 : ierr=14;
513 0 : return(ierr);
514 : }
515 : }
516 : //
517 : // Check to see if we read pass the end of the GRIB
518 : // message and missed the terminator string '7777'.
519 : //
520 12 : ipos=ipos+lensec; // Update beginning of section pointer
521 12 : if (ipos > (istart+lengrib)) {
522 0 : printf("g2_getfld: '7777' not found at end of GRIB message.\n");
523 0 : ierr=7;
524 0 : return(ierr);
525 : }
526 : //
527 : // If unpacking requested, return when all sections have been
528 : // processed
529 : //
530 12 : if (unpack && have3 && have4 && have5 && have6 && have7)
531 2 : return(ierr);
532 : //
533 : // If unpacking is not requested, return when sections
534 : // 3 through 6 have been processed
535 : //
536 10 : if ((! unpack) && have3 && have4 && have5 && have6)
537 0 : return(ierr);
538 :
539 10 : }
540 :
541 : //
542 : // If exited from above loop, the end of the GRIB message was reached
543 : // before the requested field was found.
544 : //
545 0 : printf("g2_getfld: GRIB message contained %d different fields.\n",numfld);
546 0 : printf("g2_getfld: The request was for field %d.\n",ifldnum);
547 0 : ierr=6;
548 :
549 0 : return(ierr);
550 :
551 : }
|
