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