1 /**********************************************************************
2 zyGrib: meteorological GRIB file viewer
3 Copyright (C) 2008 - Jacques Zaninetti - http://www.zygrib.org
4 
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 GNU General Public License for more details.
14 
15 You should have received a copy of the GNU General Public License
16 along with this program.  If not, see <http://www.gnu.org/licenses/>.
17 ***********************************************************************/
18 /*
19 ** File: unpackgrib2.c
20 **
21 ** Author:  Bob Dattore
22 **          NCAR/DSS
23 **          dattore@ucar.edu
24 **          (303) 497-1825
25 **	latest
26 **          14 Aug 2015:
27 **             DRS Template 5.3 (complex packing and spatial differencing)
28 
29     copyright ?
30 */
31 
32 #define __STDC_LIMIT_MACROS
33 
34 #include "wx/wxprec.h"
35 
36 #ifndef  WX_PRECOMP
37 #include "wx/wx.h"
38 #endif //precompiled headers
39 
40 #include <stdlib.h>
41 
42 #include "GribV2Record.h"
43 
44 #ifdef JASPER
45 #include <jasper/jasper.h>
46 #endif
47 
48 const double GRIB_MISSING_VALUE = GRIB_NOTDEF;
49 
50 class GRIBStatproc {
51 public:
52     int proc_code;
53     int incr_type;
54     int time_unit;
55     int time_length;
56     int incr_unit;
57     int incr_length;
58 };
59 
60 class GRIBMetadata {
61 public:
GRIBMetadata()62   GRIBMetadata() : bitmap(0), bms(0) {
63     stat_proc.t = 0;
64     lvl1_type = 0;
65     lvl2_type = 0;
66     lvl1 = 0.;
67     lvl2 = 0.;
68   };
69 
~GRIBMetadata()70   ~GRIBMetadata() {
71     delete [] stat_proc.t;
72     delete [] bitmap;
73     delete [] bms;
74   };
75 
76   int gds_templ_num;
77 
78   int earth_shape;
79   unsigned char earth_sphere_scale_factor; // Scale factor of radius of spherical Earth
80   int earth_sphere_scale_value;            // Scale value of radius of spherical Earth
81 
82   unsigned char earth_major_scale_factor; // Scale factor of major axis of oblate spheroid Earth
83   int earth_major_scale_value;            // Scaled value of major axis of oblate spheroid Earth
84 
85   unsigned char earth_minor_scale_factor; // Scale factor of minor axis of oblate spheroid Earth
86   int earth_minor_scale_value;            // Scaled value of minor axis of oblate spheroid Earth
87 
88   int nx,ny;
89   double slat,slon,latin1,latin2,splat,splon;
90   double latD;
91   union {
92     double elat;
93     double lad;
94   } lats;
95   union {
96     double elon;
97     double lov;
98   } lons;
99   union {
100     double loinc;
101     double dxinc;
102   } xinc;
103   union {
104     double lainc;
105     double dyinc;
106   } yinc;
107   int rescomp,scan_mode,proj_flag;
108   int pds_templ_num;
109   int param_cat,param_num,gen_proc,time_unit,fcst_time;
110   int ens_type,perturb_num,derived_fcst_code,nfcst_in_ensemble;
111   int lvl1_type,lvl2_type;
112   double lvl1,lvl2;
113   struct {
114     int eyr,emo,edy,etime;
115     int num_ranges, nmiss;
116     GRIBStatproc *t;
117   } stat_proc;
118   struct {
119     int stat_proc,type,num_points;
120   } spatial_proc;
121   struct {
122     int split_method,miss_val_mgmt;
123     unsigned int num_groups;
124     float primary_miss_sub,secondary_miss_sub;
125     struct {
126 	int ref,pack_width;
127     } width;
128     struct {
129 	unsigned int ref,incr,last,pack_width;
130     } length;
131     struct {
132 	unsigned int order,order_vals_width;
133     } spatial_diff;
134   } complex_pack;
135   int drs_templ_num;
136   int precision;
137   float R;
138   int E,D,num_packed,pack_width,orig_val_type;
139   int bms_ind;
140   unsigned char *bitmap;
141   ///
142   zuchar	*bms;
143   int		bmssize;
144 } ;
145 
146 class  GRIB2Grid {
147 public:
GRIB2Grid()148   GRIB2Grid() : gridpoints(0) { };
~GRIB2Grid()149   ~GRIB2Grid() { delete [] gridpoints; };
150 
151   double *gridpoints;
152 };
153 
154 class  GRIBMessage {
155 public:
GRIBMessage()156   GRIBMessage() : buffer(0) {};
~GRIBMessage()157   ~GRIBMessage() {
158       delete [] buffer;
159   };
160   unsigned char *buffer;
161   int offset;  /* offset in bytes to next GRIB2 section */
162   int total_len,disc,ed_num;
163   int center_id,sub_center_id,table_ver,local_table_ver,ref_time_type;
164   int yr,mo,dy,time;
165   int prod_status,data_type;
166   GRIBMetadata md;
167   size_t num_grids;
168   GRIB2Grid grids;
169 };
170 
171 
172 #ifdef JASPER
dec_jpeg2000(char * injpc,int bufsize,int * outfld)173 static int dec_jpeg2000(char *injpc,int bufsize,int *outfld)
174 /*$$$  SUBPROGRAM DOCUMENTATION BLOCK
175 *                .      .    .                                       .
176 * SUBPROGRAM:    dec_jpeg2000      Decodes JPEG2000 code stream
177 *   PRGMMR: Gilbert          ORG: W/NP11     DATE: 2002-12-02
178 *
179 * ABSTRACT: This Function decodes a JPEG2000 code stream specified in the
180 *   JPEG2000 Part-1 standard (i.e., ISO/IEC 15444-1) using JasPer
181 *   Software version 1.500.4 (or 1.700.2) written by the University of British
182 *   Columbia and Image Power Inc, and others.
183 *   JasPer is available at http://www.ece.uvic.ca/~mdadams/jasper/.
184 *
185 * PROGRAM HISTORY LOG:
186 * 2002-12-02  Gilbert
187 *
188 * USAGE:     int dec_jpeg2000(char *injpc,int bufsize,int *outfld)
189 *
190 *   INPUT ARGUMENTS:
191 *      injpc - Input JPEG2000 code stream.
192 *    bufsize - Length (in bytes) of the input JPEG2000 code stream.
193 *
194 *   OUTPUT ARGUMENTS:
195 *     outfld - Output matrix of grayscale image values.
196 *
197 *   RETURN VALUES :
198 *          0 = Successful decode
199 *         -2 = no memory Error.
200 *         -3 = Error decode jpeg2000 code stream.
201 *         -5 = decoded image had multiple color components.
202 *              Only grayscale is expected.
203 *
204 * REMARKS:
205 *
206 *      Requires JasPer Software version 1.500.4 or 1.700.2
207 *
208 * ATTRIBUTES:
209 *   LANGUAGE: C
210 *   MACHINE:  IBM SP
211 *
212 *$$$*/
213 
214 {
215     int ier;
216     int i,j,k;
217     jas_image_t *image=0;
218     jas_stream_t *jpcstream;
219     jas_image_cmpt_t *pcmpt;
220     char *opts=0;
221     jas_matrix_t *data;
222 
223 //    jas_init();
224 
225     ier=0;
226 //
227 //     Create jas_stream_t containing input JPEG200 codestream in memory.
228 //
229 
230     jpcstream=jas_stream_memopen(injpc,bufsize);
231     if (jpcstream == nullptr) {
232         printf(" dec_jpeg2000: no memory\n");
233         return -2;
234     }
235 //
236 //     Decode JPEG200 codestream into jas_image_t structure.
237 //
238     image=jpc_decode(jpcstream,opts);
239     if ( image == 0 ) {
240        printf(" jpc_decode return = %d \n",ier);
241        return -3;
242     }
243 
244     pcmpt=image->cmpts_[0];
245 
246 //   Expecting jpeg2000 image to be grayscale only.
247 //   No color components.
248 //
249     if (image->numcmpts_ != 1 ) {
250        printf("dec_jpeg2000: Found color image.  Grayscale expected.\n");
251        return (-5);
252     }
253 
254 //
255 //    Create a data matrix of grayscale image values decoded from
256 //    the jpeg2000 codestream.
257 //
258     data=jas_matrix_create(jas_image_height(image), jas_image_width(image));
259     jas_image_readcmpt(image,0,0,0,jas_image_width(image),
260                        jas_image_height(image),data);
261 //
262 //    Copy data matrix to output integer array.
263 //
264     k=0;
265     for (i=0;i<pcmpt->height_;i++)
266       for (j=0;j<pcmpt->width_;j++)
267         outfld[k++]=data->rows_[i][j];
268 //
269 //     Clean up JasPer work structures.
270 //
271     jas_matrix_destroy(data);
272     ier=jas_stream_close(jpcstream);
273     jas_image_destroy(image);
274 
275     return 0;
276 
277 }
278 #endif
279 
uint2(unsigned char const * p)280 static unsigned int uint2(unsigned char const *p) {
281     return (p[0] << 8) + p[1];
282 }
283 
uint4(unsigned const char * p)284 static unsigned int uint4(unsigned const char *p) {
285     return ((p[0] << 24) + (p[1] << 16) + (p[2] << 8) + p[3]);
286 }
287 
int2(unsigned const char * p)288 static int int2(unsigned const char *p) {
289     int i;
290     if ((p[0] & 0x80)) {
291         i = -(((p[0] & 0x7f) << 8) + p[1]);
292     }
293     else {
294         i = (p[0] << 8) + p[1];
295     }
296     return i;
297 }
298 
int4(unsigned const char * p)299 static int int4(unsigned const char *p) {
300     int i;
301     if ((p[0] & 0x80)) {
302         i = -(((p[0] & 0x7f) << 24) + (p[1] << 16) + (p[2] << 8) + p[3]);
303     }
304     else {
305         i = (p[0] << 24) + (p[1] << 16) + (p[2] << 8) + p[3];
306     }
307     return i;
308 }
309 
ieee2flt(unsigned const char * ieee)310 static float ieee2flt(unsigned const char *ieee) {
311     double fmant;
312     int exp;
313 
314     if ((ieee[0] & 127) == 0 && ieee[1] == 0 && ieee[2] == 0 && ieee[3] == 0)
315         return (float) 0.0;
316 
317     exp = ((ieee[0] & 127) << 1) + (ieee[1] >> 7);
318     fmant = (double) ((int) ieee[3] + (int) (ieee[2] << 8) + (int) ((ieee[1] | 128) << 16));
319     if (ieee[0] & 128)
320         fmant = -fmant;
321 
322     return (float) (ldexp(fmant, (int) (exp - 128 - 22)));
323 }
324 
getBits(unsigned const char * buf,int * loc,size_t first,size_t nbBits)325 static inline void getBits(unsigned const char *buf, int *loc, size_t first, size_t nbBits)
326 {
327     if (nbBits == 0) {
328         // x >> 32 is undefined behavior, on x86 it returns x
329         *loc = 0;
330         return;
331     }
332 
333     zuint oct = first / 8;
334     zuint bit = first % 8;
335 
336     zuint val = (buf[oct]<<24) + (buf[oct+1]<<16) + (buf[oct+2]<<8) + (buf[oct+3]);
337     val = val << bit;
338     val = val >> (32-nbBits);
339     *loc = val;
340 }
341 
342 //-------------------------------------------------------------------------------
343 // Lecture depuis un fichier
344 //-------------------------------------------------------------------------------
unpackIDS(GRIBMessage * grib_msg)345 static void unpackIDS(GRIBMessage *grib_msg)
346 {
347   int length;
348   int hh,mm,ss;
349   size_t ofs = grib_msg->offset/8;
350   unsigned char *b = grib_msg->buffer +ofs;
351 
352   length = uint4(b);   /* length of the IDS */
353 
354   grib_msg->center_id       = uint2(b +5); /* center ID */
355   grib_msg->sub_center_id   = uint2(b +7); /* sub-center ID */
356   grib_msg->table_ver       = b[9];        /* table version */
357   grib_msg->local_table_ver = b[10];       /* local table version */
358   grib_msg->ref_time_type   = b[11];       /* significance of reference time */
359   grib_msg->yr              = uint2(b+12); /* year */
360   grib_msg->mo              = b[14];       /* month */
361   grib_msg->dy              = b[15];       /* day */
362   hh                        = b[16];       /* hours */
363   mm                        = b[17];       /* minutes */
364   ss                        = b[18];       /* seconds */
365   grib_msg->time=hh*10000+mm*100+ss;
366   grib_msg->prod_status = b[19];  /* production status */
367   grib_msg->data_type   = b[20]; /* type of data */
368   grib_msg->offset += length*8;
369 }
370 
unpackLUS(GRIBMessage * grib_msg)371 static bool unpackLUS(GRIBMessage *grib_msg)
372 {
373     return true;
374 }
375 
parse_earth(GRIBMessage * grib_msg)376 static void parse_earth(GRIBMessage *grib_msg)
377 {
378   size_t ofs = grib_msg->offset/8;
379   unsigned char *b = grib_msg->buffer +ofs;
380 
381   grib_msg->md.earth_shape = b[14];         // shape of the earth
382 
383   grib_msg->md.earth_sphere_scale_factor = b[15];        // Scale factor of radius of spherical Earth
384   grib_msg->md.earth_sphere_scale_value  = uint4(b +16); // Scale value of radius of spherical Earth
385 
386   grib_msg->md.earth_major_scale_factor = b[20];         // Scale factor of major axis of oblate spheroid Earth
387   grib_msg->md.earth_major_scale_value  = uint4(b +21);  // Scaled value of major axis of oblate spheroid Earth
388 
389   grib_msg->md.earth_minor_scale_factor = b[25];         // Scale factor of minor axis of oblate spheroid Earth
390   grib_msg->md.earth_minor_scale_value = uint4(b +26);   // Scaled value of minor axis of oblate spheroid Earth
391 
392 }
393 
unpackGDS(GRIBMessage * grib_msg)394 static bool unpackGDS(GRIBMessage *grib_msg)
395 {
396   int src,num_in_list;
397   size_t ofs = grib_msg->offset/8;
398   unsigned char *b = grib_msg->buffer +ofs;
399 
400   src = b[5]; /* source of grid definition */
401   if (src != 0) {
402     fprintf(stderr,"Don't recognize predetermined grid definitions");
403     return false;
404   }
405 
406   num_in_list = b[10];  /* quasi-regular grid indication */
407   if (num_in_list > 0) {
408     fprintf(stderr,"Unable to unpack quasi-regular grids");
409     return false;
410   }
411 
412   /* grid definition template number Table 3.1 */
413   grib_msg->md.gds_templ_num = uint2(b +12);
414   switch (grib_msg->md.gds_templ_num) {
415     case 0:   /* Latitude/Longitude Also called Equidistant Cylindrical or Plate Caree */
416     case 40:  /* Gaussian Latitude/Longitude  */
417         parse_earth(grib_msg);
418 
419         grib_msg->md.nx          = uint4(b +30);  /* number of latitudes */
420         grib_msg->md.ny          = uint4(b +34);  /* number of longitudes */
421 
422 	grib_msg->md.slat = int4(b +46)/1000000.; /* latitude of first gridpoint */
423 	grib_msg->md.slon = int4(b +50)/1000000.; /* longitude of first gridpoint */
424 
425 	grib_msg->md.rescomp = b[54];             /* resolution and component flags */
426 
427 	grib_msg->md.lats.elat = int4(b +55)/1000000.; /* latitude of last gridpoint */
428 	grib_msg->md.lons.elon = int4(b +59)/1000000.; /* longitude of last gridpoint */
429 
430 	grib_msg->md.xinc.loinc = uint4(b +63)/1000000.; /* longitude increment */
431 
432 	if (grib_msg->md.gds_templ_num == 0)
433 	  grib_msg->md.yinc.lainc = uint4(b +67)/1000000.; /* latitude increment */
434 
435 	grib_msg->md.scan_mode = b[71]; /* scanning mode flag */
436 	break;
437     case 10: /* Mercator */
438         parse_earth(grib_msg);
439 
440 	grib_msg->md.nx = uint4(b +30); 	/* number of points along a parallel */
441 	grib_msg->md.ny = uint4(b +34); 	/* number of points along a meridian */
442 
443 	grib_msg->md.slat = int4(b + 38)/1000000.; /* latitude of first gridpoint */
444 	grib_msg->md.slon = int4(b + 42)/1000000.; /* longitude of first gridpoint */
445 
446 	grib_msg->md.rescomp = b[46]; /* resolution and component flags */
447 	grib_msg->md.latD    = int4(b +47)/1000000.; /* latitude at which the Mercator projection intersects the Earth
448                                                         (Latitude where Di and Dj are specified)   */
449 
450 	grib_msg->md.lats.elat = int4(b +51)/1000000.; /* latitude of last gridpoint */
451 	grib_msg->md.lons.elon = int4(b +55)/1000000.; /* longitude of last gridpoint */
452 
453         grib_msg->md.scan_mode = b[59]; /* scanning mode flag */
454 
455 	grib_msg->md.xinc.loinc = uint4(b +64)/1000.; /* longitude increment */
456 	grib_msg->md.yinc.lainc = uint4(b +68)/1000.; /* latitude increment */
457         break;
458     case 30: /* Lambert conformal grid */
459         parse_earth(grib_msg);
460 
461 	grib_msg->md.nx = uint4(b +30); /* number of points along a parallel */
462 	grib_msg->md.ny = uint4(b +34); /* number of points along a meridian */
463 
464 	grib_msg->md.slat = int4(b + 38)/1000000.; /* latitude of first gridpoint */
465 	grib_msg->md.slon = int4(b + 42)/1000000.; /* longitude of first gridpoint */
466 
467 	grib_msg->md.rescomp = b[46]; /* resolution and component flags */
468 
469 	grib_msg->md.lats.lad = int4(b +47)/1000000.; /* LaD */
470 	grib_msg->md.lons.lov = int4(b +51)/1000000.; /* LoV */
471 
472 	grib_msg->md.xinc.dxinc = int4(b +55)/1000.; /* x-direction increment */
473 	grib_msg->md.yinc.dyinc = int4(b +59)/1000.; /* y-direction increment */
474 
475 	grib_msg->md.proj_flag = b[63]; /* projection center flag */
476         grib_msg->md.scan_mode = b[64]; /* scanning mode flag */
477 
478 	grib_msg->md.latin1 = int4(b +65)/1000000.; /* latin1 */
479 	grib_msg->md.latin2 = int4(b +69)/1000000.; /* latin2 */
480 
481 	grib_msg->md.splat = int4(b +73)/1000000.; /* latitude of southern pole of projection */
482 	grib_msg->md.splon = int4(b +77)/1000000.; /* longitude of southern pole of projection */
483 	break;
484     default:
485 	fprintf(stderr,"Grid template %d is not understood\n",grib_msg->md.gds_templ_num);
486 	return false;
487   }
488   return true;
489 }
490 
unpack_stat_proc(GRIBMessage * grib_msg,unsigned const char * b)491 static void unpack_stat_proc(GRIBMessage *grib_msg, unsigned const char *b)
492 {
493   int hh,mm,ss;
494   size_t n, off;
495 
496   grib_msg->md.stat_proc.eyr = uint2(b);
497   grib_msg->md.stat_proc.emo = b[2];
498   grib_msg->md.stat_proc.edy = b[3];
499   hh = b[4];
500   mm = b[5];
501   ss = b[6];
502   grib_msg->md.stat_proc.etime = hh*10000+mm*100+ss;
503 
504   grib_msg->md.stat_proc.num_ranges = b[7];        /* number of time range specifications */
505   grib_msg->md.stat_proc.nmiss      = uint4(b +8); /* number of values missing from process */
506 
507   if (grib_msg->md.stat_proc.t != 0) {
508       delete [] grib_msg->md.stat_proc.t;
509   }
510   grib_msg->md.stat_proc.t= new GRIBStatproc[grib_msg->md.stat_proc.num_ranges];
511   off=12;
512   for (n=0; n < (size_t)grib_msg->md.stat_proc.num_ranges; n++) {
513       grib_msg->md.stat_proc.t[n].proc_code   = b[off];
514       grib_msg->md.stat_proc.t[n].incr_type   = b[off +1];
515       grib_msg->md.stat_proc.t[n].time_unit   = b[off +2];
516       grib_msg->md.stat_proc.t[n].time_length = uint4(b + off +3);
517       grib_msg->md.stat_proc.t[n].incr_unit   = b[off +7];
518       grib_msg->md.stat_proc.t[n].incr_length = uint4(b +off +8);
519       off += 12;
520   }
521 }
522 
523 // Section 4: Product Definition Section
unpackPDS(GRIBMessage * grib_msg)524 static bool unpackPDS(GRIBMessage *grib_msg)
525 {
526   int num_coords,factor;
527   size_t ofs = grib_msg->offset/8;
528   unsigned char *b = grib_msg->buffer +ofs;
529 
530   num_coords = uint2(b +5);/* indication of hybrid coordinate system */
531   if (num_coords > 0) {
532     fprintf(stderr,"Unable to decode hybrid coordinates");
533     return false;
534   }
535 
536   grib_msg->md.pds_templ_num = uint2(b +7); /* product definition template number */
537   grib_msg->md.stat_proc.num_ranges=0;
538   switch (grib_msg->md.pds_templ_num) {
539     case 0:
540     case 1:
541     case 2:
542     case 8:   // Average, accumulation, extreme values
543     case 11:
544     case 12:
545     case 15:
546 	grib_msg->md.ens_type=-1;
547 	grib_msg->md.derived_fcst_code=-1;
548 	grib_msg->md.spatial_proc.type=-1;
549 	grib_msg->md.param_cat = b[9];         /* parameter category */
550 	grib_msg->md.param_num = b[10];        /* parameter number */
551 	grib_msg->md.gen_proc  = b[11];        /* generating process */
552 
553 	grib_msg->md.time_unit = b[17];        /* time range indicator*/
554 	grib_msg->md.fcst_time = uint4(b +18); /* forecast time */
555 
556 	grib_msg->md.lvl1_type = b[22];        /* type of first level */
557 	factor = b[23];			       /* value of first level */
558 	grib_msg->md.lvl1 = int4(b +24)/pow(10.,(double)factor);
559 
560 	grib_msg->md.lvl2_type = b[28];        /* type of second level */
561 	factor = b[29];                        /* value of second level */
562 	grib_msg->md.lvl2 = int4(b +30)/pow(10.,(double)factor);
563 
564 	switch (grib_msg->md.pds_templ_num) {
565 	  case 1:
566 	  case 11:
567 	    grib_msg->md.ens_type    = b[34];
568 	    grib_msg->md.perturb_num = b[35];
569 	    grib_msg->md.nfcst_in_ensemble = b[36];
570 
571 	    switch (grib_msg->md.pds_templ_num) {
572 		case 11:
573 		  unpack_stat_proc(grib_msg, b +37);
574 		  break;
575 	    }
576 	    break;
577 	  case 2:
578 	  case 12:
579 	    grib_msg->md.derived_fcst_code = b[34];
580 	    grib_msg->md.nfcst_in_ensemble = b[35];
581 
582 	    switch (grib_msg->md.pds_templ_num) {
583 		case 12:
584 		  unpack_stat_proc(grib_msg, b +36);
585 		  break;
586 	    }
587 	    break;
588 	  case 8:
589             unpack_stat_proc(grib_msg, b +34);
590 	    break;
591 	  case 15:
592 	    grib_msg->md.spatial_proc.stat_proc  = b[34];
593 	    grib_msg->md.spatial_proc.type       = b[35];
594 	    grib_msg->md.spatial_proc.num_points = b[36];
595 	    break;
596 	}
597 	break;
598     default:
599 	fprintf(stderr,"Product Definition Template %d is not understood\n",grib_msg->md.pds_templ_num);
600 	return false;
601   }
602   return true;
603 }
604 
605 //  Section 5: Data Representation Section
unpackDRS(GRIBMessage * grib_msg)606 static bool unpackDRS(GRIBMessage *grib_msg)
607 {
608   size_t ofs = grib_msg->offset/8;
609   unsigned char *b = grib_msg->buffer +ofs;
610 
611   grib_msg->md.num_packed = uint4(b +5);    /* number of packed values */
612   grib_msg->md.drs_templ_num = uint2(b +9); /* data representation template number */
613 
614   switch (grib_msg->md.drs_templ_num) { // Table 5.0
615     case 4:        // Grid Point Data - Simple Packing
616         grib_msg->md.precision = b[11];
617         break;
618     case 0:        // Grid Point Data - Simple Packing
619     case 2:        // Grid Point Data - Complex Packing
620     case 3:        // Grid Point Data - Complex Packing and Spatial Differencing
621 #ifdef JASPER
622     case 40:       // Grid Point Data - JPEG2000 Compression
623     case 40000:
624 #endif
625    /* cf http://www.wmo.int/pages/prog/www/WMOCodes/Guides/GRIB/GRIB2_062006.pdf p. 36*/
626 	grib_msg->md.R = ieee2flt(b+ 11);
627 	grib_msg->md.E = int2(b +15);
628 	grib_msg->md.D = int2(b +17);
629 	grib_msg->md.R /= pow(10.,grib_msg->md.D);
630 
631 	grib_msg->md.pack_width = b[19];
632 	grib_msg->md.orig_val_type = b[20];
633 
634 	if (grib_msg->md.drs_templ_num == 3 || grib_msg->md.drs_templ_num == 2) {
635 	  grib_msg->md.complex_pack.split_method = b[21];
636 	  grib_msg->md.complex_pack.miss_val_mgmt = b[22];
637 	  if (grib_msg->md.orig_val_type == 0) { // Table 5.1
638 	    grib_msg->md.complex_pack.primary_miss_sub   = ieee2flt(b+ 23);
639 	    grib_msg->md.complex_pack.secondary_miss_sub = ieee2flt(b+ 27);
640 	  }
641 	  else if (grib_msg->md.orig_val_type == 1) {
642 	    grib_msg->md.complex_pack.primary_miss_sub   = uint4(b +23);
643 	    grib_msg->md.complex_pack.secondary_miss_sub = uint4(b +27);
644 	  }
645 	  else {
646 	    fprintf(stderr,"Unable to decode missing value substitutes for original value type %d\n",grib_msg->md.orig_val_type);
647 	    return false;
648 	  }
649 	  grib_msg->md.complex_pack.num_groups = uint4(b +31);
650 
651 	  grib_msg->md.complex_pack.width.ref        = b[35];
652 	  grib_msg->md.complex_pack.width.pack_width = b[36];
653 
654 	  grib_msg->md.complex_pack.length.ref        = uint4(b +37);
655 	  grib_msg->md.complex_pack.length.incr       = b[41];
656 	  grib_msg->md.complex_pack.length.last       = uint4(b +42);
657 	  grib_msg->md.complex_pack.length.pack_width = b[46];
658 	}
659 	if (grib_msg->md.drs_templ_num == 3) {
660 	  grib_msg->md.complex_pack.spatial_diff.order = b[47];
661 	  grib_msg->md.complex_pack.spatial_diff.order_vals_width = b[48];
662 	}
663 	else {
664 	  grib_msg->md.complex_pack.spatial_diff.order = 0;
665 	  grib_msg->md.complex_pack.spatial_diff.order_vals_width = 0;
666 	}
667 	break;
668     default:
669 	fprintf(stderr,"Data template %d is not understood\n",grib_msg->md.drs_templ_num);
670 	return false;
671   }
672   return true;
673 }
674 
675 //  Section 6: Bit-Map Section
unpackBMS(GRIBMessage * grib_msg)676 static bool unpackBMS(GRIBMessage *grib_msg)
677 {
678   int ind,len,n,bit;
679   size_t ofs = grib_msg->offset/8;
680   unsigned char *b = grib_msg->buffer +ofs;
681 
682   ind = b[5]; /* bit map indicator */
683   switch (ind) {
684     case 0: // A bit map applies to this product and is specified in this section.
685 	len = uint4(b);
686 	if (len < 7)
687 	    return false;
688         len -=6;
689 	grib_msg->md.bmssize = len;
690 	len *= 8;
691         delete [] grib_msg->md.bitmap;
692         delete [] grib_msg->md.bms;
693 	grib_msg->md.bitmap = new unsigned char[len];
694 	grib_msg->md.bms = new zuchar[grib_msg->md.bmssize];
695 	memcpy (grib_msg->md.bms, b + 6, grib_msg->md.bmssize);
696 	for (n=0; n < len; n++) {
697 	  getBits(grib_msg->buffer, &bit, grib_msg->offset+48+n, 1);
698 	  grib_msg->md.bitmap[n]=bit;
699 	}
700 	break;
701     case 254: // A bit map previously defined in the same GRIB2 message applies to this product.
702 	break;
703     case 255:  // A bit map does not apply to this product.
704         delete [] grib_msg->md.bitmap;
705 	grib_msg->md.bitmap=NULL;
706         delete [] grib_msg->md.bms;
707 	grib_msg->md.bms=NULL;
708 	grib_msg->md.bmssize = 0;
709 	break;
710     default:
711 	fprintf(stderr,"This code is not currently set up to deal with predefined bit-maps\n");
712 	return false;
713   }
714   return true;
715 }
716 
717 // Section 7: Data Section
unpackDS(GRIBMessage * grib_msg)718 static bool unpackDS(GRIBMessage *grib_msg)
719 {
720   int off,pval, l;
721   unsigned int n, m;
722 
723   struct {
724     int *ref_vals,*widths;
725     int *lengths;
726     int *first_vals = 0,sign,omin;
727     long long miss_val,group_miss_val;
728     int max_length;
729   } groups;
730   float lastgp,D=pow(10.,grib_msg->md.D),E=pow(2.,grib_msg->md.E);
731 
732   groups.omin = 0;
733   groups.first_vals = nullptr;
734 
735   off= grib_msg->offset+40;
736   int npoints = grib_msg->md.ny *grib_msg->md.nx;
737   switch (grib_msg->md.drs_templ_num) {
738     case 0:
739 	grib_msg->grids.gridpoints = new double[npoints];
740 	for (l=0; l < npoints; l++) {
741 	  if (grib_msg->md.bitmap == NULL || grib_msg->md.bitmap[l] == 1) {
742 	    getBits(grib_msg->buffer,&pval,off,grib_msg->md.pack_width);
743 	    grib_msg->grids.gridpoints[l]=grib_msg->md.R+pval*E/D;
744 	    off+=grib_msg->md.pack_width;
745 	  }
746 	  else
747 	    grib_msg->grids.gridpoints[l]=GRIB_MISSING_VALUE;
748 	}
749 	break;
750     case 3:
751 	if (grib_msg->md.complex_pack.num_groups > 0) {
752           if (grib_msg->md.complex_pack.spatial_diff.order) {
753 	      groups.first_vals= new int[grib_msg->md.complex_pack.spatial_diff.order];
754 	      for (n=0; n < grib_msg->md.complex_pack.spatial_diff.order; ++n) {
755 	          getBits(grib_msg->buffer,&groups.first_vals[n],off,grib_msg->md.complex_pack.spatial_diff.order_vals_width*8);
756 	          off+=grib_msg->md.complex_pack.spatial_diff.order_vals_width*8;
757               }
758 	  }
759 	  getBits(grib_msg->buffer,&groups.sign,off,1);
760 	  getBits(grib_msg->buffer,&groups.omin,off+1,grib_msg->md.complex_pack.spatial_diff.order_vals_width*8-1);
761 	  if (groups.sign == 1) {
762 	    groups.omin=-groups.omin;
763 	  }
764 	  off+=grib_msg->md.complex_pack.spatial_diff.order_vals_width*8;
765 	}
766 	// fall through
767     case 2:
768 	grib_msg->grids.gridpoints=new double[npoints];
769 	if (grib_msg->md.complex_pack.num_groups == 0) {
770 	  for (l = 0; l < npoints; ++l) {
771  	    grib_msg->grids.gridpoints[l] = GRIB_MISSING_VALUE;
772 	  }
773 	  break;
774 	}
775         if (grib_msg->md.complex_pack.miss_val_mgmt > 0) {
776 	  groups.miss_val=pow(2.,grib_msg->md.pack_width)-1;
777 	}
778 	else {
779 	  groups.miss_val=GRIB_MISSING_VALUE;
780 	}
781 
782 	groups.ref_vals = new int[grib_msg->md.complex_pack.num_groups];
783 	groups.widths   = new int[grib_msg->md.complex_pack.num_groups];
784 	groups.lengths  = new int[grib_msg->md.complex_pack.num_groups];
785 
786 	for (n=0; n < grib_msg->md.complex_pack.num_groups; ++n) {
787 	  getBits(grib_msg->buffer,&groups.ref_vals[n],off,grib_msg->md.pack_width);
788 	  off+=grib_msg->md.pack_width;
789 	}
790 	off = (off + 7) & ~7; // byte boundary padding
791 
792 	for (n=0; n < grib_msg->md.complex_pack.num_groups; ++n) {
793 	  getBits(grib_msg->buffer,&groups.widths[n],off,grib_msg->md.complex_pack.width.pack_width);
794           groups.widths[n] += grib_msg->md.complex_pack.width.ref;
795 	  off+=grib_msg->md.complex_pack.width.pack_width;
796 	}
797 	off = (off + 7) & ~7;
798 
799 	for (n=0; n < grib_msg->md.complex_pack.num_groups; ++n) {
800 	  getBits(grib_msg->buffer,&groups.lengths[n],off,grib_msg->md.complex_pack.length.pack_width);
801 	  off+=grib_msg->md.complex_pack.length.pack_width;
802 	}
803 	off = (off + 7) & ~7;
804 
805 	groups.max_length=0;
806 	for (n=0; n < grib_msg->md.complex_pack.num_groups-1; ++n) {
807 	  groups.lengths[n]=grib_msg->md.complex_pack.length.ref+groups.lengths[n]*grib_msg->md.complex_pack.length.incr;
808 	  if (groups.lengths[n] > groups.max_length) {
809 		groups.max_length=groups.lengths[n];
810 	  }
811 	}
812 	groups.lengths[n]=grib_msg->md.complex_pack.length.last;
813 	if (groups.lengths[n] > groups.max_length) {
814 	  groups.max_length=groups.lengths[n];
815 	}
816 	// unpack the field of differences
817 	for (n=0,l=0; n < grib_msg->md.complex_pack.num_groups; ++n) {
818 	  if (groups.widths[n] > 0) {
819 		if (grib_msg->md.complex_pack.miss_val_mgmt > 0) {
820 		  groups.group_miss_val=pow(2.,groups.widths[n])-1;
821 		}
822 		else {
823 		  groups.group_miss_val=GRIB_MISSING_VALUE;
824 		}
825 		for (int i = 0; i < groups.lengths[n]; ) {
826 		  if (grib_msg->md.bitmap != NULL && grib_msg->md.bitmap[l] == 0) {
827                     grib_msg->grids.gridpoints[l]=GRIB_MISSING_VALUE;
828                   }
829                   else {
830                     getBits(grib_msg->buffer,&pval,off,groups.widths[n]);
831                     off+=groups.widths[n];
832 		      if (pval == groups.group_miss_val) {
833                          grib_msg->grids.gridpoints[l]=GRIB_MISSING_VALUE;
834 		      }
835 		      else {
836 		          grib_msg->grids.gridpoints[l]=pval+groups.ref_vals[n]+groups.omin;
837                     }
838                     ++i;
839 		  }
840 		  ++l;
841 		}
842 	  }
843 	  else { // constant group XXX bitmap?
844             for (int i=0; i < groups.lengths[n]; ) {
845 		  if (grib_msg->md.bitmap != NULL && grib_msg->md.bitmap[l] == 0) {
846 		    grib_msg->grids.gridpoints[l]=GRIB_MISSING_VALUE;
847 		  }
848 		  else {
849 		    if (groups.ref_vals[n] == groups.miss_val) {
850                       grib_msg->grids.gridpoints[l]=GRIB_MISSING_VALUE;
851                     }
852                     else {
853 		        grib_msg->grids.gridpoints[l]=groups.ref_vals[n]+groups.omin;
854                     }
855 		    ++i;
856 		  }
857 		  ++l;
858             }
859 	  }
860 	}
861 
862 	for (; l < npoints; ++l) {
863 	  grib_msg->grids.gridpoints[l]=GRIB_MISSING_VALUE;
864 	}
865 
866 	if (grib_msg->md.drs_templ_num == 3) {
867       	   if (groups.first_vals != nullptr) {
868       	      for (n=grib_msg->md.complex_pack.spatial_diff.order-1; n > 0; --n) {
869   	         lastgp=groups.first_vals[n]-groups.first_vals[n-1];
870   	         for (l=0,m=0; l < grib_msg->md.nx*grib_msg->md.ny; ++l) {
871   	            if (grib_msg->grids.gridpoints[l] != GRIB_MISSING_VALUE) {
872                        if (m >= grib_msg->md.complex_pack.spatial_diff.order) {
873                           grib_msg->grids.gridpoints[l]+=lastgp;
874                           lastgp=grib_msg->grids.gridpoints[l];
875                        }
876                        ++m;
877                     }
878   	         }
879               }
880   	   }
881   	   for (l=0,m=0,lastgp=0; l < npoints; ++l) {
882   	     if (grib_msg->grids.gridpoints[l] != GRIB_MISSING_VALUE) {
883   	   	if (m < grib_msg->md.complex_pack.spatial_diff.order) {
884   	   	  grib_msg->grids.gridpoints[l]=grib_msg->md.R+groups.first_vals[m]*E/D;
885   	   	  lastgp=grib_msg->md.R*D/E+groups.first_vals[m];
886   	   	}
887   	   	else {
888   	   	  lastgp+=grib_msg->grids.gridpoints[l];
889   	   	  grib_msg->grids.gridpoints[l]=lastgp*E/D;
890   	   	}
891   	   	++m;
892   	     }
893   	   }
894   	   delete [] groups.first_vals;
895 	}
896 	else for (l=0; l < npoints; ++l) {
897   	   if (grib_msg->grids.gridpoints[l] != GRIB_MISSING_VALUE) {
898                 grib_msg->grids.gridpoints[l]= grib_msg->md.R+grib_msg->grids.gridpoints[l]*E/D;
899            }
900 	}
901 	delete [] groups.ref_vals;
902 	delete [] groups.widths;
903 	delete [] groups.lengths;
904  	break;
905     case 4:
906         {
907 
908  	    // Grid point data - IEEE Floating Point Data
909  	    if (grib_msg->md.precision == 1) { // IEEE754 single precision
910      	        grib_msg->grids.gridpoints = new double[npoints];
911      	        for (int l=0; l < npoints; l++) {
912  	            if (grib_msg->md.bitmap == nullptr || grib_msg->md.bitmap[l] == 1) {
913  	                grib_msg->grids.gridpoints[l]= ieee2flt(grib_msg->buffer +off/8);
914  	                off += 32;
915                     }
916                     else
917                         grib_msg->grids.gridpoints[l]=GRIB_MISSING_VALUE;
918                 }
919             }
920  	    else if (grib_msg->md.precision == 2) { // IEEE754 single precision
921                 static const int one = 1;
922                 bool const is_lsb = *((char*)&one) == 1;
923      	        grib_msg->grids.gridpoints = new double[npoints];
924      	        for (l=0; l < npoints; l++) {
925  	            if (grib_msg->md.bitmap == nullptr || grib_msg->md.bitmap[l] == 1) {
926  	                double d;
927  	                if( is_lsb ) {
928                             unsigned char temp[8];
929  	                    for(int j = 0; j < 8; j++ ) {
930  	                        temp[j] = grib_msg->buffer[off/8 + 7 - j];
931                             }
932                             memcpy(&d, temp, 8);
933                         }
934                         else {
935                             memcpy(&d, grib_msg->buffer +off/8, 8);
936                         }
937  	                grib_msg->grids.gridpoints[l]= d;
938  	                off += 64;
939                     }
940                     else
941                         grib_msg->grids.gridpoints[l]=GRIB_MISSING_VALUE;
942                 }
943             }
944             else {
945                 fprintf(stderr,"g2_unpack7: Invalid precision=%d for Data Section 5.4.\n", grib_msg->md.precision);
946                 return false;
947             }
948         }
949  	break;
950 #ifdef JASPER
951     case 40:
952     case 40000:
953         int len, *jvals, cnt;
954 	getBits(grib_msg->buffer,&len,grib_msg->offset,32);
955 	if (len < 5)
956 	    return false;
957 	len=len-5;
958 	jvals= new int[npoints];
959 	grib_msg->grids.gridpoints= new double[npoints];
960 	if (len > 0)
961 	  dec_jpeg2000((char *)&grib_msg->buffer[grib_msg->offset/8+5],len,jvals);
962 	cnt=0;
963 	for (l=0; l < npoints; l++) {
964 	  if (grib_msg->md.bitmap == NULL || grib_msg->md.bitmap[l] == 1) {
965 	    if (len == 0)
966 		jvals[cnt]=0;
967 	    grib_msg->grids.gridpoints[l]=grib_msg->md.R+jvals[cnt++]*E/D;
968 	  }
969 	  else
970 	    grib_msg->grids.gridpoints[l]=GRIB_MISSING_VALUE;
971 	}
972 	delete [] jvals;
973 	break;
974 #endif
975     default:
976         erreur("Unknown packing %d", grib_msg->md.drs_templ_num);
977         break;
978   }
979   return true;
980 }
981 
GRBV2_TO_DATA(int productDiscipline,int dataCat,int dataNum)982 static zuchar GRBV2_TO_DATA(int productDiscipline, int dataCat, int dataNum)
983 {
984     zuchar ret = 255;
985     // printf("search %d %d %d\n", productDiscipline, dataCat,  dataNum);
986     switch (productDiscipline) { // TABLE 4.2
987     case 0:      // Meteorological products
988         switch (dataCat) {
989         case 0:  // Temperature
990             switch (dataNum) {
991             case 0: ret = GRB_TEMP; break; // DATA_TO_GRBV2[DATA_TEMP] = grb2DataType(0,0,0);
992             case 2: ret= GRB_TPOT; break;  // DATA_TO_GRBV2[DATA_TEMP_POT] = grb2DataType(0,0,2);
993             case 4: ret = GRB_TMAX; break; // DATA_TO_GRBV2[DATA_TMAX] = grb2DataType(0,0,4);
994             case 5: ret = GRB_TMIN; break; // DATA_TO_GRBV2[DATA_TMIN] = grb2DataType(0,0,5);
995             case 6: ret = GRB_DEWPOINT; break; //DATA_TO_GRBV2[DATA_DEWPOINT] = grb2DataType(0,0,6);
996             }
997             break;
998         case 1: // dataCat Moisture
999             switch (dataNum) {
1000             case 0: ret = GRB_HUMID_SPEC; break; //DATA_TO_GRBV2[DATA_HUMID_SPEC] = grb2DataType(0,1,0);
1001             case 1: ret = GRB_HUMID_REL; break; // DATA_TO_GRBV2[DATA_HUMID_REL] = grb2DataType(0,1,1);
1002             case 7: ret= GRB_PRECIP_RATE; break; // DATA_TO_GRBV2[DATA_PRECIP_RATE] = grb2DataType(0,1,7);
1003             case 49:							 // Total Water Precipitation (Meteo France Arome 0.01
1004             case 52:                             // Total precipitation rate kg m–2 s–1
1005             case 8: ret = GRB_PRECIP_TOT; break; // DATA_TO_GRBV2[DATA_PRECIP_TOT] = grb2DataType(0,1,8);
1006             case 11: ret = GRB_SNOW_DEPTH; break; // DATA_TO_GRBV2[DATA_SNOW_DEPTH] = grb2DataType(0,1,11);
1007             case 193: ret = GRB_FRZRAIN_CATEG; break; // DATA_TO_GRBV2[DATA_FRZRAIN_CATEG] = grb2DataType(0,1,193);
1008             case 195: ret = GRB_SNOW_CATEG; break; //DATA_TO_GRBV2[DATA_SNOW_CATEG] = grb2DataType(0,1,195);
1009             }
1010             break;
1011         case 2: // dataCat  Momentum
1012             switch (dataNum) {
1013             case 0: ret = GRB_WIND_DIR; break;
1014             case 1: ret = GRB_WIND_SPEED; break;
1015             case 2: ret = GRB_WIND_VX; break; // DATA_TO_GRBV2[DATA_WIND_VX] = grb2DataType(0,2,2);
1016             case 3: ret = GRB_WIND_VY; break; // DATA_TO_GRBV2[DATA_WIND_VY] = grb2DataType(0,2,3);
1017             case 22: ret = GRB_WIND_GUST; break; //
1018             }
1019             break;
1020         case 3: // dataCat mass
1021             switch (dataNum) {
1022             case 0: ret = GRB_PRESSURE; break; //DATA_TO_GRBV2[DATA_PRESSURE] = grb2DataType(0,3,0);
1023             case 1: ret = GRB_PRESSURE; break; // PRSMSL //DATA_TO_GRBV2[DATA_PRESSURE] = grb2DataType(0,3,0);
1024             case 5: ret = GRB_GEOPOT_HGT; break; // DATA_TO_GRBV2[DATA_GEOPOT_HGT]= grb2DataType(0,3,5);
1025 
1026             case 192: ret = GRB_PRESSURE; break; //DATA_TO_GRBV2[DATA_MSLET] = grb2DataType(0,3,192);
1027             }
1028             break;
1029         case 6: // dataCat
1030             switch (dataNum) {
1031             case 1: ret = GRB_CLOUD_TOT; break; // DATA_TO_GRBV2[DATA_CLOUD_TOT] = grb2DataType(0,6,1);
1032             }
1033             break;
1034         case 7:// dataCat
1035             switch (dataNum) {
1036             // case 7: ret =  GRB_?; break // DATA_TO_GRBV2[DATA_CIN] = grb2DataType(0,7,7);
1037             case 6: ret = GRB_CAPE; break; // DATA_TO_GRBV2[DATA_CAPE] = grb2DataType(0,7,6);
1038             }
1039             break;
1040         case 16: // Meteorological products, Forecast Radar Imagery category
1041             switch (dataNum) {
1042             case 196: ret =  GRB_COMP_REFL; break; // = grb2DataType(0,16, 196);
1043             }
1044             break;
1045         }
1046         break;
1047     case 10: // productDiscipline Oceanographic products
1048         switch (dataCat) {
1049         case 0:         // waves
1050 #if 0
1051             switch (dataNum) {
1052             case 3: ret= GRB_WVHGT; break; //DATA_TO_GRBV2[DATA_WAVES_SIG_HGT_COMB] = grb2DataType(10,0,3);
1053             DATA_TO_GRBV2[DATA_WAVES_WND_DIR] = grb2DataType(10,0,4);
1054             DATA_TO_GRBV2[DATA_WAVES_WND_HGT] = grb2DataType(10,0,5);
1055             DATA_TO_GRBV2[DATA_WAVES_WND_PERIOD] = grb2DataType(10,0,6);
1056             DATA_TO_GRBV2[DATA_WAVES_SWL_DIR] = grb2DataType(10,0,7);
1057             DATA_TO_GRBV2[DATA_WAVES_SWL_HGT] = grb2DataType(10,0,8);
1058             DATA_TO_GRBV2[DATA_WAVES_SWL_PERIOD] = grb2DataType(10,0,9);
1059             DATA_TO_GRBV2[DATA_WAVES_PRIM_DIR] = grb2DataType(10,0,10);
1060             DATA_TO_GRBV2[DATA_WAVES_PRIM_PERIOD] = grb2DataType(10,0,11);
1061             DATA_TO_GRBV2[DATA_WAVES_SEC_DIR] = grb2DataType(10,0,12);
1062             DATA_TO_GRBV2[DATA_WAVES_SEC_PERIOD] = grb2DataType(10,0,13);
1063             }
1064 #endif
1065 
1066             switch (dataNum) {
1067                 case 3: ret= GRB_HTSGW; break; // Significant Height of Combined Wind Waves and Swell
1068                 case 4: ret= GRB_WVDIR; break; // Direction of Wind Waves
1069                 case 5: ret= GRB_WVHGT; break; // Significant Height of Wind Waves
1070                 case 6: ret= GRB_WVPER; break; // Mean Period of Wind Waves
1071                 case 14: ret= GRB_DIR; break;  // Direction of Combined Wind Waves and Swell
1072                 case 15: ret= GRB_PER; break;  // Mean Period of Combined Wind Waves and Swell
1073             }
1074             break;
1075 
1076         case 1: // Currents
1077             switch (dataNum) {
1078                 case 0: ret = GRB_CUR_DIR; break;
1079                 case 1: ret = GRB_CUR_SPEED; break;
1080                 case 2: ret = GRB_UOGRD; break; // DATA_TO_GRBV2[DATA_CURRENT_VX] = grb2DataType(10,1,2);
1081                 case 3: ret = GRB_VOGRD; break; // DATA_TO_GRBV2[DATA_CURRENT_VY] = grb2DataType(10,1,3);
1082             }
1083             break;
1084          case 3: // Surface Properties
1085             switch (dataNum) {
1086                 case 0: ret = GRB_WTMP; break; // DATA_TO_GRBV2[DATA_CURRENT_VX] = grb2DataType(10,1,2);
1087             }
1088             break;
1089 
1090         }
1091         break;
1092     }
1093 #if 1
1094     if (ret == 255) {
1095         erreur("unknown Discipline %d dataCat %d dataNum %d", productDiscipline,  dataCat, dataNum);
1096     }
1097 #endif
1098     return ret;
1099 }
1100 
1101 /** Return UINT_MAX on errors. */
mapStatisticalEndTime(GRIBMessage * grid)1102 static int mapStatisticalEndTime(GRIBMessage *grid)
1103 {
1104    // lovely md.fcst_time is in grid->md.time_unit but md.stat_proc.t[0].time_length is in grid->md.stat_proc.t[0].time_unit
1105    // not always the same.
1106   if (grid->md.time_unit == grid->md.stat_proc.t[0].time_unit) switch (grid->md.time_unit) { // table 4.4
1107     case 0:  // minute
1108 	// return (grid->md.stat_proc.etime/100 % 100)-(grid->time/100 % 100);
1109     case 1:  // hour
1110          return grid->md.fcst_time +grid->md.stat_proc.t[0].time_length;
1111 	 // return (grid->md.stat_proc.etime/10000- grid->time/10000);
1112     case 2:  // Day
1113 	return (grid->md.stat_proc.edy -grid->dy);
1114     case 3:
1115 	return (grid->md.stat_proc.emo -grid->mo);
1116     case 4:
1117 	return (grid->md.stat_proc.eyr -grid->yr);
1118     default:
1119 	fprintf(stderr,"Unable to map end time with units %d to GRIB1\n",grid->md.time_unit);
1120 	return UINT_MAX;
1121   }
1122 
1123   if (grid->md.time_unit == 0 && grid->md.stat_proc.t[0].time_unit == 1) {
1124          // in minute + hourly increment
1125          return grid->md.fcst_time +grid->md.stat_proc.t[0].time_length *60;
1126   }
1127 
1128   if (grid->md.time_unit == 1 && grid->md.stat_proc.t[0].time_unit == 0 && (grid->md.stat_proc.t[0].time_unit  % 60) != 0 ) {
1129           // convert in hour
1130          return grid->md.fcst_time +grid->md.stat_proc.t[0].time_length /60;
1131   }
1132 
1133   fprintf(stderr, "Unable to map end time %d %d %d %d \n", grid->md.time_unit, grid->md.stat_proc.t[0].time_unit, grid->md.fcst_time,
1134             grid->md.stat_proc.t[0].time_length);
1135   return UINT_MAX;
1136 }
1137 
1138 // map GRIB2 msg time to GRIB1 P1 and P2 in sec
mapTimeRange(GRIBMessage * grid,zuint * p1,zuint * p2,zuchar * t_range,int * n_avg,int * n_missing,int center)1139 static bool mapTimeRange(GRIBMessage *grid, zuint *p1, zuint *p2, zuchar *t_range,int *n_avg,int *n_missing, int center)
1140 {
1141   switch (grid->md.pds_templ_num) {
1142     case 0:
1143     case 1:
1144     case 2:
1145     case 15:
1146 	*t_range=0;
1147 	*p1=grid->md.fcst_time;
1148 	*p2=0;
1149 	*n_avg=*n_missing=0;
1150 	break;
1151     case 8:
1152     case 11:
1153     case 12:
1154 	if (grid->md.stat_proc.num_ranges > 1) {
1155 	  if (center == 7 && grid->md.stat_proc.num_ranges == 2) {
1156 /* NCEP CFSR monthly grids */
1157 	    *p2=grid->md.stat_proc.t[0].incr_length;
1158 	    *p1=*p2 -grid->md.stat_proc.t[1].time_length;
1159 	    *n_avg=grid->md.stat_proc.t[0].time_length;
1160 	    switch (grid->md.stat_proc.t[0].proc_code) {
1161 		case 193:
1162 		  *t_range=113;
1163 		  break;
1164 		case 194:
1165 		  *t_range=123;
1166 		  break;
1167 		case 195:
1168 		  *t_range=128;
1169 		  break;
1170 		case 196:
1171 		  *t_range=129;
1172 		  break;
1173 		case 197:
1174 		  *t_range=130;
1175 		  break;
1176 		case 198:
1177 		  *t_range=131;
1178 		  break;
1179 		case 199:
1180 		  *t_range=132;
1181 		  break;
1182 		case 200:
1183 		  *t_range=133;
1184 		  break;
1185 		case 201:
1186 		  *t_range=134;
1187 		  break;
1188 		case 202:
1189 		  *t_range=135;
1190 		  break;
1191 		case 203:
1192 		  *t_range=136;
1193 		  break;
1194 		case 204:
1195 		  *t_range=137;
1196 		  break;
1197 		case 205:
1198 		  *t_range=138;
1199 		  break;
1200 		case 206:
1201 		  *t_range=139;
1202 		  break;
1203 		case 207:
1204 		  *t_range=140;
1205 		  break;
1206 		default:
1207 		  fprintf(stderr,"Unable to map NCEP statistical process code %d to GRIB1\n",grid->md.stat_proc.t[0].proc_code);
1208 		  return false;
1209 	    }
1210 	  }
1211 	  else {
1212 	    fprintf(stderr,"Unable to map multiple statistical processes to GRIB1\n");
1213 	    return false;
1214 	  }
1215 	}
1216 	else {
1217 	  switch (grid->md.stat_proc.t[0].proc_code) {
1218 	    case 0:
1219 	    case 1:
1220 	    case 4:
1221 		switch (grid->md.stat_proc.t[0].proc_code) {
1222 		  case 0: /* average */
1223 		    *t_range=3;
1224 		    break;
1225 		  case 1: /* accumulation */
1226 		    *t_range=4;
1227 		    break;
1228 		  case 4: /* difference */
1229 		    *t_range=5;
1230 		    break;
1231 		}
1232 		*p1=grid->md.fcst_time;
1233 		*p2=mapStatisticalEndTime(grid);
1234                 if (*p2 == UINT_MAX) {
1235                     return false;
1236                 }
1237 		if (grid->md.stat_proc.t[0].incr_length == 0)
1238 		  *n_avg=0;
1239 		else {
1240 		  fprintf(stderr,"Unable to map discrete processing to GRIB1\n");
1241 		  return false;
1242 		}
1243 		break;
1244 
1245 	    case 2: // maximum
1246 	    case 3: // minimum
1247 		*t_range=2;
1248 		*p1=grid->md.fcst_time;
1249 		*p2=mapStatisticalEndTime(grid);
1250                 if (*p2 == UINT_MAX) {
1251                     return false;
1252                 }
1253 		if (grid->md.stat_proc.t[0].incr_length == 0)
1254 		  *n_avg=0;
1255 		else {
1256 		  fprintf(stderr,"Unable to map discrete processing to GRIB1\n");
1257 		  return false;
1258 		}
1259 		break;
1260 	    default:
1261 // patch for NCEP grids
1262 		if (grid->md.stat_proc.t[0].proc_code == 255 && center == 7) {
1263  		  if (grid->disc == 0) {
1264 		    if (grid->md.param_cat == 0) {
1265 			switch (grid->md.param_num) {
1266 			  case 4:
1267 			  case 5:
1268 			    *t_range=2;
1269 			    *p1=grid->md.fcst_time;
1270 			    *p2=mapStatisticalEndTime(grid);
1271                             if (*p2 == UINT_MAX) {
1272                                 return false;
1273                             }
1274 			    if (grid->md.stat_proc.t[0].incr_length == 0)
1275 				*n_avg=0;
1276 			    else {
1277 				fprintf(stderr,"Unable to map discrete processing to GRIB1\n");
1278 				return false;
1279 			    }
1280 			    break;
1281 			}
1282 		    }
1283 		  }
1284 		}
1285 		else {
1286 		  fprintf(stderr,"Unable to map statistical process %d to GRIB1\n",grid->md.stat_proc.t[0].proc_code);
1287 		  return false;
1288 		}
1289 	  }
1290 	}
1291 	*n_missing=grid->md.stat_proc.nmiss;
1292 	break;
1293     default:
1294 	fprintf(stderr,"Unable to map time range for Product Definition Template %d into GRIB1\n",grid->md.pds_templ_num);
1295 	return false;
1296   }
1297   return true;
1298 }
1299 
1300 //-------------------------------------------------------------------------------
1301 // Adjust data type from different mete center
1302 //-------------------------------------------------------------------------------
translateDataType()1303 void  GribV2Record::translateDataType()
1304 {
1305     this->knownData = true;
1306     dataCenterModel = OTHER_DATA_CENTER;
1307     //------------------------
1308     // NOAA GFS
1309     //------------------------
1310     if (dataType == GRB_PRECIP_RATE) {	// mm/s -> mm/h
1311         multiplyAllData( 3600.0 );
1312     }
1313     if ( idCenter==7 && idModel==2 )		// NOAA
1314     {
1315         dataCenterModel = NOAA_GFS;
1316         // altitude level (entire atmosphere vs entire atmosphere considered as 1 level)
1317         if (levelType == LV_ATMOS_ENT) {
1318             levelType = LV_ATMOS_ALL;
1319         }
1320         if (dataType == GRB_TEMP          //gfs Water surface Temperature
1321             && levelType == LV_GND_SURF
1322             && levelValue == 0) dataType = GRB_WTMP;
1323     }
1324     //------------------------
1325 	//DNMI-NEurope.grb
1326     //------------------------
1327     else if ( idCenter==7 && idModel==88 && idGrid==255 ) {  // saildocs
1328         dataCenterModel = NOAA_NCEP_WW3;
1329     }
1330     //----------------------------
1331     //NOAA RTOFS
1332     //--------------------------------
1333     else if(idCenter==7 && idModel==45 && idGrid==255) {
1334         dataCenterModel = NOAA_RTOFS;
1335     }
1336     //----------------------------------------------
1337     // NCEP sea surface temperature
1338     //----------------------------------------------
1339     else if ((idCenter==7 && idModel==44 && idGrid==173)
1340         || (idCenter==7 && idModel==44 && idGrid==235))
1341     {
1342         dataCenterModel = NOAA_NCEP_SST;
1343     }
1344     //----------------------------------------------
1345     // FNMOC WW3 mediterranean sea
1346     //----------------------------------------------
1347     else if (idCenter==58 && idModel==111 && idGrid==179)
1348     {
1349         dataCenterModel = FNMOC_WW3_MED;
1350     }
1351     //----------------------------------------------
1352     // FNMOC WW3
1353     //----------------------------------------------
1354     else if (idCenter==58 && idModel==110 && idGrid==240)
1355     {
1356         dataCenterModel = FNMOC_WW3_GLB;
1357     }
1358 	//------------------------
1359 	// Meteorem (Scannav)
1360 	//------------------------
1361     else if (idCenter==59 && idModel==78 && idGrid==255)
1362     {
1363         //dataCenterModel = ??
1364 		if ( (getDataType()==GRB_WIND_VX || getDataType()==GRB_WIND_VY)
1365 			&& getLevelType()==LV_MSL
1366 			&& getLevelValue()==0)
1367 		{
1368 			levelType  = LV_ABOV_GND;
1369 			levelValue = 10;
1370 		}
1371 		if ( getDataType()==GRB_PRECIP_TOT
1372 			&& getLevelType()==LV_MSL
1373 			&& getLevelValue()==0)
1374 		{
1375 			levelType  = LV_GND_SURF;
1376 			levelValue = 0;
1377 		}
1378 	}
1379     else if (idCenter==84 && idModel <= 5 && idGrid==0)
1380     {
1381     }
1382 
1383 	//------------------------
1384 	// Unknown center
1385 	//------------------------
1386 	else
1387 	{
1388         dataCenterModel = OTHER_DATA_CENTER;
1389 //		printf("Uncorrected GribRecord: ");
1390 //		this->print();
1391 //		this->knownData = false;
1392 
1393 	}
1394     //translate significant wave height and dir
1395     if (this->knownData) {
1396         switch (levelType) {
1397             case 100: // LV_ISOBARIC
1398                 /* GRIB1 is in hectoPascal
1399                    GRIB2 in Pascal, convert to GRIB1
1400                 */
1401                 levelValue = levelValue /100;
1402                 break;
1403             case 103: levelType = LV_ABOV_GND;break;
1404             case 101: levelType = LV_MSL;break;
1405         }
1406         switch (getDataType()) {
1407             case GRB_WIND_GUST:
1408                 levelType  = LV_GND_SURF;
1409                 levelValue = 0;
1410                 break;
1411             case GRB_UOGRD:
1412             case GRB_VOGRD:
1413                 break;
1414             case GRB_WVHGT:
1415             case GRB_HTSGW:
1416             case GRB_WVDIR:
1417             case GRB_WVPER:
1418             case GRB_DIR:
1419             case GRB_PER:
1420                 levelType  = LV_GND_SURF;
1421                 levelValue = 0;
1422                 break;
1423         }
1424     }
1425     // this->print();
1426 }
1427 
1428 // -------------------------------------
readDataSet(ZUFILE * file)1429 void GribV2Record::readDataSet(ZUFILE* file)
1430 {
1431     bool skip = false;
1432     bool DS = false;
1433     int len, sec_num;
1434 
1435     data    = NULL;
1436     BMSbits = NULL;
1437     hasBMS = false;
1438     knownData = false;
1439     IsDuplicated = false;
1440 
1441     while (strncmp(&((char *)grib_msg->buffer)[grib_msg->offset/8],"7777",4) != 0) {
1442         DS = false;
1443         getBits(grib_msg->buffer, &len, grib_msg->offset, 32);
1444         getBits(grib_msg->buffer, &sec_num, grib_msg->offset +4*8, 8);
1445         switch (sec_num) {
1446 	case 2: //  Section 2: Local Use Section
1447 	     if (skip == true)  break;
1448              ok = unpackLUS(grib_msg);
1449              break;
1450 	case 3: //  Section 3: Grid Definition Section
1451 	     if (skip == true)  break;
1452 	     ok = unpackGDS(grib_msg);
1453 	     if (ok) {
1454 	         Ni = grib_msg->md.nx;
1455 	         Nj = grib_msg->md.ny;
1456 	         La1 = grib_msg->md.slat;
1457 	         Lo1 = grib_msg->md.slon;
1458 	         La2 = grib_msg->md.lats.elat;
1459 	         Lo2 = grib_msg->md.lons.elon;
1460 	         Di = grib_msg->md.xinc.loinc;
1461 	         Dj = grib_msg->md.yinc.lainc;
1462                  scanFlags = grib_msg->md.scan_mode;
1463                  isScanIpositive = (scanFlags&0x80) ==0;
1464                  isScanJpositive = (scanFlags&0x40) !=0;
1465                  isAdjacentI     = (scanFlags&0x20) ==0;
1466                  if (Lo1>=0 && Lo1<=180 && Lo2<0)
1467                      Lo2 += 360.0;    // cross the 180 deg meridien,beetwen alaska and russia
1468 
1469 	         if (isScanIpositive) while ( Lo1> Lo2 ) {   // horizontal size > 360 °
1470 	             Lo1 -= 360.0;
1471                  }
1472                  if (Lo2 > Lo1) {
1473                      lonMin = Lo1;
1474                      lonMax = Lo2;
1475                   }
1476                   else {
1477                      lonMin = Lo2;
1478                      lonMax = Lo1;
1479                   }
1480                   if (La2 > La1) {
1481                      latMin = La1;
1482                      latMax = La2;
1483                   }
1484                   else {
1485                      latMin = La2;
1486                      latMax = La1;
1487                   }
1488                   if (Ni<=1 || Nj<=1) {
1489                       erreur("Record %d: Ni=%d Nj=%d",id,Ni,Nj);
1490                       ok = false;
1491                   }
1492                   else {
1493                       Di = (Lo2-Lo1) / (Ni-1);
1494                       Dj = (La2-La1) / (Nj-1);
1495                   }
1496 	     }
1497 	     break;
1498 	case 4: //  Section 4: Product Definition Section
1499 	     if (skip == true)  break;
1500 	     ok = unpackPDS(grib_msg);
1501 	     if (ok) {
1502 	         // printf("template %d 0 meteo data cat %d data num %d\n", grib_msg->md.pds_templ_num, grib_msg->md.param_cat, grib_msg->md.param_num);
1503 	         productTemplate = grib_msg->md.pds_templ_num;
1504 	         dataCat = grib_msg->md.param_cat;
1505                  dataNum = grib_msg->md.param_num;
1506                  dataType= GRBV2_TO_DATA(productDiscipline,dataCat,dataNum);
1507                  if (dataType == 255) {
1508                      //printf("unused data type, skip\n");
1509                      skip = true;
1510                      break;
1511                  }
1512 
1513 	         levelType = grib_msg->md.lvl1_type;
1514 	         levelValue = grib_msg->md.lvl1;
1515 	         if (grib_msg->md.lvl2_type == 8 && grib_msg->md.lvl1_type == 1) {
1516 	             // cf table 4.5:  8 Nominal top of the atmosphere
1517 	             levelType = LV_ATMOS_ALL;
1518 	             levelValue = 0.;
1519 	         }
1520 	         int n_avg, n_missing;
1521 
1522 	         if (!mapTimeRange(grib_msg, &periodP1 , &periodP2, &timeRange , &n_avg, &n_missing, idCenter)) {
1523 	             skip = true;
1524 	             break;
1525                  }
1526 	         periodsec = periodSeconds(grib_msg->md.time_unit, periodP1, periodP2, timeRange);
1527 	         setRecordCurrentDate(makeDate(refyear,refmonth,refday,refhour,refminute,periodsec));
1528 	         //printf("%d %d %d %d %d %d \n", refyear,refmonth,refday,refhour,refminute,periodsec);
1529 	         //printf("%d Periode %d P1=%d p2=%d %s\n", grib_msg->md.time_unit, periodsec, periodP1,periodP2, strCurDate);
1530 
1531 	     }
1532 	     break;
1533 	case 5: //  Section 5: Data Representation Section
1534 	     if (skip == true)  break;
1535 	     ok = unpackDRS(grib_msg);
1536 	     break;
1537 	case 6: //  Section 6: Bit-Map Section
1538 	     if (skip == true)  break;
1539 	     ok = unpackBMS(grib_msg);
1540 	     if (ok) {
1541 	        if (grib_msg->md.bmssize != 0) {
1542 	             hasBMS = true;
1543 	             BMSsize = grib_msg->md.bmssize;
1544 	             BMSbits = new zuchar[grib_msg->md.bmssize];
1545 	             memcpy (BMSbits, grib_msg->md.bms, grib_msg->md.bmssize);
1546                 }
1547 	     }
1548 	     break;
1549 	case 7:  // Section 7: Data Section
1550 	     if (skip == false) {
1551     	         ok = unpackDS(grib_msg);
1552     	         if (ok) {
1553 	             data = grib_msg->grids.gridpoints;
1554 	             grib_msg->grids.gridpoints = 0;
1555                  }
1556 	     }
1557 	     if (grib_msg->num_grids != 1)
1558     	         DS = true;
1559 	     break;
1560         }
1561         grib_msg->offset += len*8;
1562         if (ok == false || DS == true )
1563             break;
1564     }
1565 
1566     //ok = false;
1567 if (false) {
1568 //if (true) {
1569 printf("==== GV2 %d\n", ok);
1570 printf("Lo1=%f Lo2=%f    La1=%f La2=%f\n", Lo1,Lo2,La1,La2);
1571 printf("Lo1=%f Lo2=%f    La1=%f La2=%f\n", grib_msg->md.slon, grib_msg->md.lons.elon, grib_msg->md.slat, grib_msg->md.lats.elat);
1572 printf("Ni=%d Nj=%d\n", Ni,Nj);
1573 printf("hasDiDj=%d Di,Dj=(%f %f)\n", hasDiDj, Di,Dj);
1574 printf("isScanIpositive=%d isScanJpositive=%d isAdjacentI=%d\n",isScanIpositive,isScanJpositive,isAdjacentI );
1575 printf("hasBMS=%d\n", hasBMS);
1576 }
1577     if (ok) {
1578         if (!skip)
1579         {
1580 		translateDataType();
1581 		setDataType(dataType);
1582         }
1583     }
1584     if (!ok || !DS || strncmp(&((char *)grib_msg->buffer)[grib_msg->offset/8],"7777",4) == 0) {
1585         delete grib_msg;
1586         grib_msg = 0;
1587     }
1588 }
1589 
1590 // -----------------
GribV2Record(ZUFILE * file,int id_)1591 GribV2Record::GribV2Record(ZUFILE* file, int id_)
1592 {
1593     id = id_;
1594     seekStart = zu_tell(file);           // moved to section 0 read
1595     data    = NULL;
1596     BMSsize = 0;
1597     BMSbits = NULL;
1598     hasBMS = false;
1599     eof     = false;
1600     knownData = false;
1601     IsDuplicated = false;
1602     long start = seekStart;
1603 
1604     grib_msg = new GRIBMessage();
1605 
1606     //      Pre read 4 bytes to check for length adder needed for some GRIBS (like WRAMS and NAM)
1607     char strgrib[5];
1608     if (zu_read(file, strgrib, 4) != 4)
1609     {
1610             ok = false;
1611             eof = true;
1612             return;
1613     }
1614 
1615     bool b_haveReadGRIB = false;         // already read the "GRIB" of section 0 ??
1616 
1617     if (strncmp(strgrib, "GRIB", 4) != 0)
1618             b_len_add_8 = true;
1619     else
1620     {
1621             b_len_add_8 = false;
1622             b_haveReadGRIB = true;
1623     }
1624 
1625     // Another special case, where zero padding is used between records.
1626     if((strgrib[0] == 0) &&
1627         (strgrib[1] == 0) &&
1628         (strgrib[2] == 0) &&
1629         (strgrib[3] == 0))
1630     {
1631           b_len_add_8 = false;
1632           b_haveReadGRIB = false;
1633     }
1634     ok = readGribSection0_IS(file, b_haveReadGRIB ); // Section 0: Indicator Section
1635 
1636     int len, sec_num;
1637     if (ok) {
1638         unpackIDS(grib_msg);  // Section 1: Identification Section
1639         int off;
1640         /* find out how many grids are in this message */
1641         off = grib_msg->offset /8;
1642         while (strncmp(&((char *)grib_msg->buffer)[off], "7777", 4) != 0) {
1643             len = uint4(grib_msg->buffer +off);
1644             sec_num = grib_msg->buffer[off+4];
1645             if (sec_num == 7)
1646                 grib_msg->num_grids++;
1647             off += len;
1648         }
1649     }
1650     else {
1651         // seek back if V1
1652         (void)zu_seek(file, start, SEEK_SET);
1653         return;
1654     }
1655     refyear  = grib_msg->yr;
1656     refmonth = grib_msg->mo;
1657     refday   = grib_msg->dy;
1658     refhour  = grib_msg->time /10000;
1659     refminute = (grib_msg->time/100) % 100;
1660     refDate = makeDate(refyear,refmonth,refday,refhour,refminute,0);
1661     sprintf(strRefDate, "%04d-%02d-%02d %02d:%02d", refyear,refmonth,refday,refhour,refminute);
1662     idCenter = grib_msg->center_id;
1663     idModel  = grib_msg->table_ver;
1664     idGrid   = 0; // FIXME data1[6];
1665     productDiscipline = grib_msg->disc;
1666     readDataSet(file);
1667 }
1668 
1669 // ---------------------------------------
hasMoreDataSet() const1670 bool GribV2Record::hasMoreDataSet() const
1671 {
1672     return grib_msg && grib_msg->num_grids != 1?true:false;
1673 }
1674 
1675 // ---------------------------------------
GribV2NextDataSet(ZUFILE * file,int id_)1676 GribV2Record *GribV2Record::GribV2NextDataSet(ZUFILE* file, int id_)
1677 {
1678     GribV2Record *rec1 = new GribV2Record(*this);
1679     // XXX should have a shallow copy constructor
1680     delete [] rec1->data;
1681     delete [] rec1->BMSbits;
1682     // new records take ownership
1683     this->grib_msg = 0;
1684     rec1->id = id_;
1685     rec1->readDataSet(file);
1686     return rec1;
1687 }
1688 
1689 //-------------------------------------------------------------------------------
1690 // Constructeur de recopie
1691 //-------------------------------------------------------------------------------
1692 #pragma warning(disable: 4717)
GribV2Record(const GribRecord & rec)1693 GribV2Record::GribV2Record(const GribRecord &rec) : GribRecord(rec)
1694 {
1695     *this = rec;
1696     #pragma warning(default: 4717)
1697 }
1698 
~GribV2Record()1699 GribV2Record::~GribV2Record()
1700 {
1701     delete grib_msg;
1702 }
1703 
1704 //==============================================================
1705 // Lecture des données
1706 //==============================================================
1707 //----------------------------------------------
1708 // SECTION 0: THE INDICATOR SECTION (IS)
1709 //----------------------------------------------
unpackIS(ZUFILE * fp,GRIBMessage * grib_msg)1710 static bool unpackIS(ZUFILE* fp, GRIBMessage *grib_msg)
1711 {
1712   unsigned char temp[16];
1713   int status;
1714   size_t num;
1715 
1716   if (grib_msg->buffer != NULL) {
1717     delete [] grib_msg->buffer;
1718     grib_msg->buffer=NULL;
1719   }
1720   grib_msg->num_grids = 0;
1721 
1722   if ( (status = zu_read(fp, &temp[4], 12)) != 12)
1723   {
1724     return false;
1725   }
1726   grib_msg->disc = temp[6];
1727   grib_msg->ed_num = temp[7];
1728 
1729   //  Bail out early if this is not GRIB2
1730   if(grib_msg->ed_num != 2)
1731       return false;
1732 
1733   getBits(temp,&grib_msg->total_len,96,32);
1734   // too small or overflow
1735   if ( grib_msg->total_len < 16 || grib_msg->total_len > (INT_MAX - 4))
1736       return false;
1737 
1738   grib_msg->md.nx = grib_msg->md.ny = 0;
1739   grib_msg->buffer = new unsigned char[grib_msg->total_len+4];
1740   memcpy(grib_msg->buffer,temp,16);
1741   num = grib_msg->total_len -16;
1742 
1743   status = zu_read(fp, &grib_msg->buffer[16], num);
1744   if (status != (int)num)
1745     return false;
1746 
1747   if (strncmp(&((char *)grib_msg->buffer)[grib_msg->total_len-4],"7777",4) != 0)
1748         fprintf(stderr,"Warning: no end section found\n");
1749 
1750   grib_msg->offset=128;
1751   return true;
1752 }
1753 
readGribSection0_IS(ZUFILE * file,bool b_skip_initial_GRIB)1754 bool GribV2Record::readGribSection0_IS(ZUFILE* file, bool b_skip_initial_GRIB) {
1755     char strgrib[4];
1756     fileOffset0 = zu_tell(file);
1757 
1758     if(!b_skip_initial_GRIB)
1759     {
1760             // Cherche le 1er 'G'
1761       while ( (zu_read(file, strgrib, 1) == 1)
1762                               &&  (strgrib[0] != 'G') )
1763             { }
1764 
1765       if (strgrib[0] != 'G') {
1766             ok = false;
1767             eof = true;
1768             return false;
1769       }
1770       if (zu_read(file, strgrib+1, 3) != 3) {
1771             ok = false;
1772             eof = true;
1773             return false;
1774       }
1775       if (strncmp(strgrib, "GRIB", 4) != 0)  {
1776             printf("readGribSection0_IS(): Unknown file header : %c%c%c%c\n", strgrib[0],strgrib[1],strgrib[2],strgrib[3]);
1777             ok = false;
1778             eof = true;
1779             return false;
1780       }
1781     }
1782 
1783     seekStart = zu_tell(file) - 4;
1784     // totalSize = readInt3(file);
1785     if (unpackIS(file , grib_msg) == false) {
1786         ok = false;
1787         eof = true;
1788         return false;
1789     }
1790 
1791     editionNumber = grib_msg->ed_num;
1792     if (editionNumber != 2)  {
1793         ok = false;
1794         eof = true;
1795         return false;
1796     }
1797 
1798     return true;
1799 }
1800 
1801 //==============================================================
1802 // Fonctions utiles
1803 //==============================================================
periodSeconds(zuchar unit,zuint P1,zuint P2,zuchar range)1804 zuint GribV2Record::periodSeconds(zuchar unit,zuint P1,zuint P2,zuchar range) {
1805     zuint res, dur;
1806 
1807     switch (unit) {
1808         case 0: //      Minute
1809             res = 60; break;
1810         case 1: //      Hour
1811             res = 3600; break;
1812         case 2: //      Day
1813             res = 86400; break;
1814         case 10: //     3 hours
1815             res = 10800; break;
1816         case 11: //     6 hours
1817             res = 21600; break;
1818         case 12: //     12 hours
1819             res = 43200; break;
1820         case 13: // Second
1821             res = 1; break;
1822         case 3: //      Month
1823         case 4: //      Year
1824         case 5: //      Decade (10 years)
1825         case 6: //      Normal (30 years)
1826         case 7: //      Century (100 years)
1827         default:
1828             erreur("id=%d: unknown time unit in PDS b18=%d",id,unit);
1829             res = 0;
1830             ok = false;
1831     }
1832     grib_debug("id=%d: PDS unit %d (time range) b21=%d %d P1=%d P2=%d\n",id,unit, range,res,P1,P2);
1833     dur = 0;
1834 
1835     switch (range) {
1836         case 0:
1837             dur = (zuint)P1; break;
1838         case 1:
1839             dur = 0; break;
1840 
1841         case 2:
1842         case 3: // Average  (reference time + P1 to reference time + P2)
1843             // dur = ((zuint)P1+(zuint)P2)/2; break;     // TODO
1844             dur = (zuint)P2; break;
1845 
1846          case 4: // Accumulation  (reference time + P1 to reference time + P2)
1847             dur = (zuint)P2; break;
1848 
1849         case 10:
1850             dur = ((zuint)P1<<8) + (zuint)P2; break;
1851         default:
1852             erreur("id=%d: unknown time range in PDS b21=%d",id,range);
1853             dur = 0;
1854             ok = false;
1855     }
1856     return res*dur;
1857 }
1858 
1859 
1860 //===============================================================================================
1861 
1862