1 #ifdef HAVE_CONFIG_H
2 #include "config.h"
3 #endif
4
5 #ifdef HAVE_LIBGRIB_API
6
7 #include "dmemory.h"
8 #include "cdi.h"
9 #include "cdi_int.h"
10 #include "file.h"
11 #include "gribapi_utilities.h"
12 #include "stream_scan.h"
13 #include "stream_grb.h"
14 #include "stream_gribapi.h"
15 #include "varscan.h"
16 #include "datetime.h"
17 #include "vlist.h"
18 #include "calendar.h"
19 #include "subtype.h"
20
21 #include "cgribex.h" /* gribGetSize, gribRead, gribGetZip, GRIB1_LTYPE_99 */
22 #include "gribapi.h"
23
24 #include <grib_api.h>
25
26 extern int CDI_Inventory_Mode;
27
28 static const var_tile_t dummy_tiles = { 0, -1, -1, -1, -1, -1 };
29
30
31 typedef struct {
32 int param;
33 int level1;
34 int level2;
35 int ltype;
36 int tsteptype;
37 size_t gridsize;
38 char name[32];
39 VarScanKeys scanKeys;
40 var_tile_t tiles;
41 } compvar2_t;
42
43
44 static
gribapiGetZaxisType(long editionNumber,int grib_ltype)45 int gribapiGetZaxisType(long editionNumber, int grib_ltype)
46 {
47 return (editionNumber <= 1) ? grib1ltypeToZaxisType(grib_ltype) : grib2ltypeToZaxisType(grib_ltype);
48 }
49
50 static
getTimeunits(const long unitsOfTime)51 int getTimeunits(const long unitsOfTime)
52 {
53 switch (unitsOfTime)
54 {
55 case 13: return TUNIT_SECOND;
56 case 0: return TUNIT_MINUTE;
57 case 1: return TUNIT_HOUR;
58 case 10: return TUNIT_3HOURS;
59 case 11: return TUNIT_6HOURS;
60 case 12: return TUNIT_12HOURS;
61 case 2: return TUNIT_DAY;
62 }
63
64 return TUNIT_HOUR;
65 }
66
67 static
timeunit_factor(const int tu1,const int tu2)68 double timeunit_factor(const int tu1, const int tu2)
69 {
70 if (tu2 == TUNIT_HOUR)
71 {
72 switch (tu1)
73 {
74 case TUNIT_SECOND: return 3600;
75 case TUNIT_MINUTE: return 60;
76 case TUNIT_HOUR: return 1;
77 case TUNIT_3HOURS: return 1./3;
78 case TUNIT_6HOURS: return 1./6;
79 case TUNIT_12HOURS: return 1./12;
80 case TUNIT_DAY: return 1./24;
81 }
82 }
83
84 return 1;
85 }
86
87 static
gribapiGetTimeUnits(grib_handle * gh)88 int gribapiGetTimeUnits(grib_handle *gh)
89 {
90 long unitsOfTime = -1;
91 grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
92
93 GRIB_CHECK(my_grib_set_long(gh, "stepUnits", unitsOfTime), 0);
94
95 return getTimeunits(unitsOfTime);
96 }
97
98 static
gribapiGetSteps(grib_handle * gh,int timeunits,int * startStep,int * endStep)99 void gribapiGetSteps(grib_handle *gh, int timeunits, int *startStep, int *endStep)
100 {
101 long unitsOfTime;
102 int status = grib_get_long(gh, "stepUnits", &unitsOfTime);
103 const int timeunits2 = (status == 0) ? getTimeunits(unitsOfTime) : timeunits;
104 //timeunits2 = gribapiGetTimeUnits(gh);
105
106 long lpar;
107 status = grib_get_long(gh, "forecastTime", &lpar);
108 if ( status == 0 ) *startStep = (int) lpar;
109 else
110 {
111 status = grib_get_long(gh, "startStep", &lpar);
112 if ( status == 0 )
113 *startStep = (int) (((double)lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
114 }
115
116 *endStep = *startStep;
117 status = grib_get_long(gh, "endStep", &lpar);
118 if ( status == 0 )
119 *endStep = (int) (((double)lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
120 // printf("%d %d %d %d %d %g\n", *startStep, *endStep, lpar, timeunits, timeunits2, timeunit_factor(timeunits, timeunits2));
121 }
122
123 static
gribapiGetDataDateTime(grib_handle * gh,int64_t * datadate,int * datatime)124 void gribapiGetDataDateTime(grib_handle *gh, int64_t *datadate, int *datatime)
125 {
126 long date;
127 GRIB_CHECK(grib_get_long(gh, "dataDate", &date), 0);
128 *datadate = (int64_t) date;
129 long hour, minute, second;
130 GRIB_CHECK(grib_get_long(gh, "hour", &hour), 0);
131 GRIB_CHECK(grib_get_long(gh, "minute", &minute), 0);
132 GRIB_CHECK(grib_get_long(gh, "second", &second), 0);
133 *datatime = cdiEncodeTime((int)hour, (int)minute, (int)second);
134 }
135
136 static
gribapiSetDataDateTime(grib_handle * gh,int64_t datadate,int datatime)137 void gribapiSetDataDateTime(grib_handle *gh, int64_t datadate, int datatime)
138 {
139 GRIB_CHECK(my_grib_set_long(gh, "dataDate", (long)datadate), 0);
140 int hour, minute, second;
141 cdiDecodeTime(datatime, &hour, &minute, &second);
142 GRIB_CHECK(my_grib_set_long(gh, "hour", hour), 0);
143 GRIB_CHECK(my_grib_set_long(gh, "minute", minute), 0);
144 GRIB_CHECK(my_grib_set_long(gh, "second", second), 0);
145 }
146
147 static
gribapiGetTimeUnitFactor(const int timeUnits)148 int gribapiGetTimeUnitFactor(const int timeUnits)
149 {
150 static bool lprint = true;
151 switch (timeUnits)
152 {
153 case TUNIT_SECOND: return 1;
154 case TUNIT_MINUTE: return 60;
155 case TUNIT_HOUR: return 3600;
156 case TUNIT_3HOURS: return 10800;
157 case TUNIT_6HOURS: return 21600;
158 case TUNIT_12HOURS: return 43200;
159 case TUNIT_DAY: return 86400;
160 default:
161 if ( lprint )
162 {
163 Warning("Time unit %d unsupported", timeUnits);
164 lprint = false;
165 }
166 }
167
168 return 0;
169 }
170
171 static
gribapiGetValidityDateTime(grib_handle * gh,int64_t * vdate,int * vtime,int64_t * sDate,int * sTime)172 void gribapiGetValidityDateTime(grib_handle *gh, int64_t *vdate, int *vtime, int64_t *sDate, int *sTime)
173 {
174 *sDate = 0;
175 *sTime = 0;
176
177 long sigofrtime = 3;
178 if ( gribEditionNumber(gh) > 1 )
179 GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
180 else
181 GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &sigofrtime), 0);
182
183 if ( sigofrtime == 3 ) //XXX: This looks like a bug to me, because timeRangeIndicator == 3 does not seem to have the same meaning as significanceOfReferenceTime == 3. I would recommend replacing this condition with `if(!gribapiTimeIsFC())`.
184 {
185 gribapiGetDataDateTime(gh, vdate, vtime);
186 }
187 else
188 {
189 int64_t rdate;
190 int rtime;
191 gribapiGetDataDateTime(gh, &rdate, &rtime);
192
193 const int timeUnits = gribapiGetTimeUnits(gh);
194 int startStep = 0, endStep = 0;
195 gribapiGetSteps(gh, timeUnits, &startStep, &endStep);
196
197 {
198 int ryear, rmonth, rday, rhour, rminute, rsecond;
199 cdiDecodeDate(rdate, &ryear, &rmonth, &rday);
200 cdiDecodeTime(rtime, &rhour, &rminute, &rsecond);
201
202 if ( rday > 0 )
203 {
204 extern int CGRIBEX_grib_calendar;
205 int64_t julday;
206 int secofday;
207 encode_caldaysec(CGRIBEX_grib_calendar, ryear, rmonth, rday, rhour, rminute, rsecond, &julday, &secofday);
208
209 int64_t time_unit_factor = gribapiGetTimeUnitFactor(timeUnits);
210
211 //if (startStep > 0)
212 {
213 int64_t julday_x = julday;
214 int secofday_x = secofday;
215 julday_add_seconds(time_unit_factor * startStep, &julday_x, &secofday_x);
216 decode_caldaysec(CGRIBEX_grib_calendar, julday_x, secofday_x, &ryear, &rmonth, &rday, &rhour, &rminute, &rsecond);
217 *sDate = cdiEncodeDate(ryear, rmonth, rday);
218 *sTime = cdiEncodeTime(rhour, rminute, 0);
219 }
220
221 julday_add_seconds(time_unit_factor * endStep, &julday, &secofday);
222 decode_caldaysec(CGRIBEX_grib_calendar, julday, secofday, &ryear, &rmonth, &rday, &rhour, &rminute, &rsecond);
223 }
224
225 *vdate = cdiEncodeDate(ryear, rmonth, rday);
226 *vtime = cdiEncodeTime(rhour, rminute, rsecond);
227 }
228 }
229 }
230
231 static
grib1GetLevel(grib_handle * gh,int * leveltype,int * lbounds,int * level1,int * level2)232 void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
233 {
234 *leveltype = 0;
235 *lbounds = 0;
236 *level1 = 0;
237 *level2 = 0;
238
239 long lpar;
240 if ( !grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar) ) //1 byte
241 {
242 *leveltype = (int) lpar;
243
244 switch (*leveltype)
245 {
246 case GRIB1_LTYPE_SIGMA_LAYER:
247 case GRIB1_LTYPE_HYBRID_LAYER:
248 case GRIB1_LTYPE_LANDDEPTH_LAYER:
249 { *lbounds = 1; break; }
250 }
251
252 if ( *lbounds )
253 {
254 GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0); //1 byte
255 if (lpar == GRIB_MISSING_LONG) lpar = 0;
256 *level1 = (int)lpar;
257 GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0); //1 byte
258 if (lpar == GRIB_MISSING_LONG) lpar = 0;
259 *level2 = (int)lpar;
260 }
261 else
262 {
263 double dlevel;
264 GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0); //2 byte
265 if ( *leveltype == GRIB1_LTYPE_ISOBARIC ) dlevel *= 100;
266 if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
267 if ( *leveltype == GRIB1_LTYPE_99 || *leveltype == GRIB1_LTYPE_ISOBARIC_PA ) *leveltype = GRIB1_LTYPE_ISOBARIC;
268
269 *level1 = (int) dlevel;
270 *level2 = 0;
271 }
272 }
273 }
274
275 static
grib2ScaleFactor(const long factor)276 double grib2ScaleFactor(const long factor)
277 {
278 switch(factor)
279 {
280 case GRIB_MISSING_LONG: return 1;
281 case -9: return 1000000000;
282 case -8: return 100000000;
283 case -7: return 10000000;
284 case -6: return 1000000;
285 case -5: return 100000;
286 case -4: return 10000;
287 case -3: return 1000;
288 case -2: return 100;
289 case -1: return 10;
290 case 0: return 1;
291 case 1: return 0.1;
292 case 2: return 0.01;
293 case 3: return 0.001;
294 case 4: return 0.0001;
295 case 5: return 0.00001;
296 case 6: return 0.000001;
297 case 7: return 0.0000001;
298 case 8: return 0.00000001;
299 case 9: return 0.000000001;
300 default: return 0;
301 }
302 }
303
304 static
calcLevel(int level_sf,long factor,long level)305 int calcLevel(int level_sf, long factor, long level)
306 {
307 double result = 0;
308 if (level != GRIB_MISSING_LONG) result = (double)level*grib2ScaleFactor(factor);
309 if (level_sf) result *= level_sf;
310 return (int)result;
311 }
312
313 static
grib2GetLevel(grib_handle * gh,int * leveltype1,int * leveltype2,int * lbounds,int * level1,int * level2,int * level_sf,int * level_unit)314 void grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1,
315 int *level2, int *level_sf, int *level_unit)
316 {
317 *leveltype1 = 0;
318 *leveltype2 = -1;
319 *lbounds = 0;
320 *level1 = 0;
321 *level2 = 0;
322 *level_sf = 0;
323 *level_unit = 0;
324
325 long lpar;
326 int status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar); //1 byte
327 if (status == 0)
328 {
329 *leveltype1 = (int) lpar;
330
331 status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar); //1 byte
332 /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
333 if (status == 0) *leveltype2 = (int)lpar;
334
335 if ( *leveltype1 != 255 && *leveltype2 != 255 && *leveltype2 > 0 ) *lbounds = 1;
336 switch (*leveltype1)
337 {
338 case GRIB2_LTYPE_REFERENCE:
339 if(*leveltype2 == 1) *lbounds = 0;
340 break;
341 case GRIB2_LTYPE_LANDDEPTH:
342 *level_sf = 1000;
343 *level_unit = CDI_UNIT_M;
344 break;
345 case GRIB2_LTYPE_ISOBARIC:
346 *level_sf = 1000;
347 *level_unit = CDI_UNIT_PA;
348 break;
349 case GRIB2_LTYPE_SIGMA:
350 *level_sf = 1000;
351 *level_unit = 0;
352 break;
353 }
354
355 long factor, llevel;
356 GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0); //1 byte
357 GRIB_CHECK(grib_get_long(gh, "scaledValueOfFirstFixedSurface", &llevel), 0); //4 byte
358 *level1 = calcLevel(*level_sf, factor, llevel);
359
360 if ( *lbounds )
361 {
362 GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0); //1 byte
363 GRIB_CHECK(grib_get_long(gh, "scaledValueOfSecondFixedSurface", &llevel), 0); //4 byte
364 *level2 = calcLevel(*level_sf, factor, llevel);
365 }
366 }
367 }
368
369 static
gribGetLevel(grib_handle * gh,int * leveltype1,int * leveltype2,int * lbounds,int * level1,int * level2,int * level_sf,int * level_unit,var_tile_t * tiles)370 void gribGetLevel(grib_handle *gh, int* leveltype1, int* leveltype2, int* lbounds, int* level1, int* level2, int* level_sf, int* level_unit, var_tile_t* tiles)
371 {
372 if ( gribEditionNumber(gh) <= 1 )
373 {
374 grib1GetLevel(gh, leveltype1, lbounds, level1, level2);
375 *leveltype2 = -1;
376 *level_sf = 0;
377 *level_unit = 0;
378 }
379 else
380 {
381 grib2GetLevel(gh, leveltype1, leveltype2, lbounds, level1, level2, level_sf, level_unit);
382
383 // read in tiles attributes (if there are any)
384 tiles->tileindex = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TILEINDEX], 0);
385 tiles->totalno_of_tileattr_pairs = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TOTALNO_OF_TILEATTR_PAIRS], -1);
386 tiles->tileClassification = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TILE_CLASSIFICATION], -1);
387 tiles->numberOfTiles = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_NUMBER_OF_TILES], -1);
388 tiles->numberOfAttributes = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_NUMBER_OF_ATTR], -1);
389 tiles->attribute = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TILEATTRIBUTE], -1);
390 }
391 }
392
393 static
gribapiGetString(grib_handle * gh,const char * key,char * string,size_t length)394 void gribapiGetString(grib_handle *gh, const char *key, char *string, size_t length)
395 {
396 string[0] = 0;
397
398 int ret = grib_get_string(gh, key, string, &length);
399 if (ret != 0)
400 {
401 fprintf(stderr, "grib_get_string(gh, \"%s\", ...) failed!\n", key);
402 GRIB_CHECK(ret, 0);
403 }
404 if ( length == 8 && memcmp(string, "unknown", length) == 0 ) string[0] = 0;
405 else if ( length == 2 && memcmp(string, "~", length) == 0 ) string[0] = 0;
406 }
407
408 static
gribapiGetEnsembleInfo(grib_handle * gh,long * typeOfEnsembleForecast,long * numberOfForecastsInEnsemble,long * perturbationNumber)409 int gribapiGetEnsembleInfo(grib_handle *gh, long *typeOfEnsembleForecast, long *numberOfForecastsInEnsemble, long *perturbationNumber)
410 {
411 int status = 0;
412 if (grib_get_long(gh, "typeOfEnsembleForecast", typeOfEnsembleForecast) == 0)
413 {
414 GRIB_CHECK(grib_get_long(gh, "numberOfForecastsInEnsemble", numberOfForecastsInEnsemble), 0);
415 GRIB_CHECK(grib_get_long(gh, "perturbationNumber", perturbationNumber), 0);
416 if (*perturbationNumber > 0) status = 1;
417 }
418
419 if (status == 0)
420 {
421 *typeOfEnsembleForecast = 0;
422 *perturbationNumber = 0;
423 *numberOfForecastsInEnsemble = 0;
424 }
425
426 return status;
427 }
428
429 static
gribapiGetScanKeys(grib_handle * gh)430 VarScanKeys gribapiGetScanKeys(grib_handle *gh)
431 {
432 VarScanKeys scanKeys;
433 varScanKeysInit(&scanKeys);
434
435 long typeOfEnsembleForecast = 0, numberOfForecastsInEnsemble = 0, perturbationNumber = 0;
436 gribapiGetEnsembleInfo(gh, &typeOfEnsembleForecast, &numberOfForecastsInEnsemble, &perturbationNumber);
437 scanKeys.perturbationNumber = (short)perturbationNumber;
438
439 long typeOfGeneratingProcess = 0;
440 if (grib_get_long(gh, "typeOfGeneratingProcess", &typeOfGeneratingProcess) == 0)
441 scanKeys.typeOfGeneratingProcess = (short)typeOfGeneratingProcess;
442
443 return scanKeys;
444 }
445
446 static
gribapiGetNameKeys(grib_handle * gh,int varID)447 void gribapiGetNameKeys(grib_handle *gh, int varID)
448 {
449 char string[CDI_MAX_NAME];
450
451 size_t vlen = CDI_MAX_NAME;
452 gribapiGetString(gh, "name", string, vlen); // longname
453 if (string[0]) varDefKeyString(varID, CDI_KEY_LONGNAME, string);
454
455 gribapiGetString(gh, "units", string, vlen);
456 if (string[0]) varDefKeyString(varID, CDI_KEY_UNITS, string);
457
458 string[0] = 0;
459 const int status = grib_get_string(gh, "cfName", string, &vlen);
460 if ( status != 0 || vlen <= 1 || strncmp(string, "unknown", 7) == 0 ) string[0] = 0;
461 if (string[0]) varDefKeyString(varID, CDI_KEY_STDNAME, string);
462 }
463
464 static
gribapiGetKeys(grib_handle * gh,int varID)465 void gribapiGetKeys(grib_handle *gh, int varID)
466 {
467 long tablesVersion = 0;
468 if ( grib_get_long(gh, "tablesVersion", &tablesVersion) == 0 )
469 varDefKeyInt(varID, CDI_KEY_TABLESVERSION, (int) tablesVersion);
470
471 long localTablesVersion = 0;
472 if ( grib_get_long(gh, "localTablesVersion", &localTablesVersion) == 0 )
473 varDefKeyInt(varID, CDI_KEY_LOCALTABLESVERSION, (int) localTablesVersion);
474
475 long typeOfGeneratingProcess = 0;
476 if ( grib_get_long(gh, "typeOfGeneratingProcess", &typeOfGeneratingProcess) == 0 )
477 varDefKeyInt(varID, CDI_KEY_TYPEOFGENERATINGPROCESS, (int) typeOfGeneratingProcess);
478
479 long productDefinitionTemplate = 0;
480 if ( grib_get_long(gh, "productDefinitionTemplateNumber", &productDefinitionTemplate) == 0 )
481 varDefKeyInt(varID, CDI_KEY_PRODUCTDEFINITIONTEMPLATE, (int) productDefinitionTemplate);
482
483 long typeOfProcessedData = 0;
484 if ( grib_get_long(gh, "typeOfProcessedData", &typeOfProcessedData) == 0 )
485 varDefKeyInt(varID, CDI_KEY_TYPEOFPROCESSEDDATA, (int) typeOfProcessedData);
486
487 long shapeOfTheEarth = 0;
488 if ( grib_get_long(gh, "shapeOfTheEarth", &shapeOfTheEarth) == 0 )
489 varDefKeyInt(varID, CDI_KEY_SHAPEOFTHEEARTH, (int) shapeOfTheEarth);
490
491 long backgroundProcess = 0;
492 if ( grib_get_long(gh, "backgroundProcess", &backgroundProcess) == 0 )
493 varDefKeyInt(varID, CDI_KEY_BACKGROUNDPROCESS, (int) backgroundProcess);
494
495 long typeOfTimeIncrement = 0;
496 if ( grib_get_long(gh, "typeOfTimeIncrement", &typeOfTimeIncrement) == 0 )
497 varDefKeyInt(varID, CDI_KEY_TYPEOFTIMEINCREMENT, (int) typeOfTimeIncrement);
498 /*
499 long constituentType = 0;
500 if ( grib_get_long(gh, "constituentType", &constituentType) == 0 )
501 varDefKeyInt(varID, CDI_KEY_CONSTITUENTTYPE, (int) constituentType);
502 */
503
504 /*
505 Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure.
506 Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars"
507 */
508 long typeOfEnsembleForecast = 0, numberOfForecastsInEnsemble = 0, perturbationNumber = 0;
509 gribapiGetEnsembleInfo(gh, &typeOfEnsembleForecast, &numberOfForecastsInEnsemble, &perturbationNumber);
510 if ( perturbationNumber > 0 )
511 {
512 varDefKeyInt(varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, (int) typeOfEnsembleForecast);
513 varDefKeyInt(varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, (int) numberOfForecastsInEnsemble);
514 varDefKeyInt(varID, CDI_KEY_PERTURBATIONNUMBER, (int) perturbationNumber);
515 }
516
517 long section2Length = 0;
518 int status = grib_get_long(gh, "section2Length", §ion2Length);
519 if (status == 0 && section2Length > 0)
520 {
521 long grib2LocalSectionNumber;
522 long mpimType, mpimClass, mpimUser;
523 status = grib_get_long(gh, "grib2LocalSectionNumber", &grib2LocalSectionNumber);
524 if (status == 0)
525 {
526 size_t section2PaddingLength = 0;
527 status = grib_get_size(gh, "section2Padding", §ion2PaddingLength);
528 if (status == 0 && section2PaddingLength > 0)
529 {
530 varDefKeyInt(varID, CDI_KEY_GRIB2LOCALSECTIONNUMBER, (int) grib2LocalSectionNumber);
531 varDefKeyInt(varID, CDI_KEY_SECTION2PADDINGLENGTH, (int) section2PaddingLength);
532 unsigned char *section2Padding = (unsigned char*) Malloc(section2PaddingLength);
533 grib_get_bytes(gh, "section2Padding", section2Padding, §ion2PaddingLength);
534 varDefKeyBytes(varID, CDI_KEY_SECTION2PADDING, section2Padding, (int)section2PaddingLength);
535 Free(section2Padding);
536 }
537 else if ( grib_get_long(gh, "mpimType", &mpimType) == 0 &&
538 grib_get_long(gh, "mpimClass", &mpimClass) == 0 &&
539 grib_get_long(gh, "mpimUser", &mpimUser) == 0)
540 {
541 varDefKeyInt(varID, CDI_KEY_MPIMTYPE, mpimType);
542 varDefKeyInt(varID, CDI_KEY_MPIMCLASS, mpimClass);
543 varDefKeyInt(varID, CDI_KEY_MPIMUSER, mpimUser);
544
545 size_t revNumLen = 20;
546 unsigned char revNumber[revNumLen];
547 if ( grib_get_bytes(gh, "revNumber", revNumber, &revNumLen) == 0 )
548 varDefKeyBytes(varID, CDI_KEY_REVNUMBER, revNumber, (int)revNumLen);
549
550 long revStatus;
551 grib_get_long(gh, "revStatus", &revStatus);
552 varDefKeyInt(varID, CDI_KEY_REVSTATUS, revStatus);
553 }
554 }
555 }
556 }
557
558 static
gribapiDefProjRLL(grib_handle * gh,int gridID)559 void gribapiDefProjRLL(grib_handle *gh, int gridID)
560 {
561 double xpole = 0, ypole = 0, angle = 0;
562 grib_get_double(gh, "latitudeOfSouthernPoleInDegrees", &ypole);
563 grib_get_double(gh, "longitudeOfSouthernPoleInDegrees", &xpole);
564 grib_get_double(gh, "angleOfRotation", &angle);
565 xpole -= 180;
566 if ( fabs(ypole) > 0 ) ypole = -ypole; // change from south to north pole
567 if ( fabs(angle) > 0 ) angle = -angle;
568
569 gridDefParamRLL(gridID, xpole, ypole, angle);
570 }
571
572 static
shapeOfTheEarthToRadius(long shapeOfTheEarth,double * a,double * b,double * rf)573 void shapeOfTheEarthToRadius(long shapeOfTheEarth, double *a, double *b, double *rf)
574 {
575 // clang-format off
576 switch (shapeOfTheEarth)
577 {
578 case 2: *a = 6378160.0; *b = 6356775.0; *rf = 297.0; break;
579 case 4: *a = 6378137.0; *b = 6356752.314; *rf = 298.257222101; break;
580 case 6: *a = 6371229.0; break;
581 case 8: *a = 6371200.0; break;
582 case 0: *a = 6367470.0; break;
583 default: *a = 6367470.0;
584 }
585 // clang-format on
586 }
587
588 static
gribapiDefProjLCC(grib_handle * gh,int gridID)589 void gribapiDefProjLCC(grib_handle *gh, int gridID)
590 {
591 struct CDI_GridProjParams gpp;
592 gridProjParamsInit(&gpp);
593
594 long shapeOfTheEarth;
595 grib_get_long(gh, "shapeOfTheEarth", &shapeOfTheEarth);
596 shapeOfTheEarthToRadius(shapeOfTheEarth, &gpp.a, &gpp.b, &gpp.rf);
597
598 long projflag = 0;
599 grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &gpp.xval_0);
600 grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &gpp.yval_0);
601 grib_get_double(gh, "LoVInDegrees", &gpp.lon_0);
602 grib_get_double(gh, "Latin1InDegrees", &gpp.lat_1);
603 grib_get_double(gh, "Latin2InDegrees", &gpp.lat_2);
604 grib_get_long(gh, "projectionCentreFlag", &projflag);
605 bool lsouth = gribbyte_get_bit((int)projflag, 1);
606 if (lsouth) { gpp.lat_1 = -gpp.lat_1; gpp.lat_2 = -gpp.lat_2; }
607
608 gpp.lat_0 = gpp.lat_2;
609
610 if (proj_lonlat_to_lcc_func)
611 {
612 double x_0 = gpp.xval_0, y_0 = gpp.yval_0;
613 proj_lonlat_to_lcc_func(gpp, (size_t)1, &x_0, &y_0);
614 if (IS_NOT_EQUAL(x_0, gpp.mv) && IS_NOT_EQUAL(y_0, gpp.mv))
615 { gpp.x_0 = -x_0; gpp.y_0 = -y_0; }
616 }
617
618 gridDefParamsLCC(gridID, gpp);
619 }
620
621
622 static
gribapiDefProjSTERE(grib_handle * gh,int gridID)623 void gribapiDefProjSTERE(grib_handle *gh, int gridID)
624 {
625 struct CDI_GridProjParams gpp;
626 gridProjParamsInit(&gpp);
627
628 long shapeOfTheEarth;
629 grib_get_long(gh, "shapeOfTheEarth", &shapeOfTheEarth);
630 shapeOfTheEarthToRadius(shapeOfTheEarth, &gpp.a, &gpp.b, &gpp.rf);
631
632 grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &gpp.xval_0);
633 grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &gpp.yval_0);
634 grib_get_double(gh, "LaDInDegrees", &gpp.lat_1);
635 grib_get_double(gh, "orientationOfTheGridInDegrees", &gpp.lon_0);
636
637 long southPoleOnProjectionPlane;
638 grib_get_long(gh, "southPoleOnProjectionPlane", &southPoleOnProjectionPlane);
639 gpp.lat_0 = southPoleOnProjectionPlane ? -90.0 : 90.0;
640
641 if (proj_lonlat_to_stere_func)
642 {
643 double x_0 = gpp.xval_0, y_0 = gpp.yval_0;
644 proj_lonlat_to_stere_func(gpp, (size_t)1, &x_0, &y_0);
645 if (IS_NOT_EQUAL(x_0, gpp.mv) && IS_NOT_EQUAL(y_0, gpp.mv))
646 { gpp.x_0 = -x_0; gpp.y_0 = -y_0; }
647 }
648
649 gridDefParamsSTERE(gridID, gpp);
650 }
651
652 static
gribapiAddRecord(stream_t * streamptr,int param,grib_handle * gh,size_t recsize,off_t position,int datatype,int comptype,const char * varname,int leveltype1,int leveltype2,int lbounds,int level1,int level2,int level_sf,int level_unit,VarScanKeys * scanKeys,const var_tile_t * tiles,int lread_additional_keys,void * fdbItem)653 void gribapiAddRecord(stream_t *streamptr, int param, grib_handle *gh,
654 size_t recsize, off_t position, int datatype, int comptype, const char *varname,
655 int leveltype1, int leveltype2, int lbounds, int level1, int level2, int level_sf, int level_unit,
656 VarScanKeys *scanKeys, const var_tile_t *tiles, int lread_additional_keys, void *fdbItem)
657 {
658 const int vlistID = streamptr->vlistID;
659 const int tsID = streamptr->curTsID;
660 const int recID = recordNewEntry(streamptr, tsID);
661 record_t *record = &streamptr->tsteps[tsID].records[recID];
662
663 const int tsteptype = gribapiGetTsteptype(gh);
664 int numavg = 0;
665
666 // fprintf(stderr, "param %d %d %d %d\n", param, level1, level2, leveltype1);
667
668 record->size = recsize;
669 record->position = position;
670 record->param = param;
671 record->ilevel = level1;
672 record->ilevel2 = level2;
673 record->ltype = leveltype1;
674 record->tsteptype = (short)tsteptype;
675 record->gridsize = gribapiGetGridsize(gh);
676 record->scanKeys = *scanKeys;
677 record->tiles = tiles ? *tiles : dummy_tiles;
678 #ifdef HAVE_LIBFDB5
679 record->fdbItem = fdbItem;
680 #endif
681
682 strncpy(record->varname, varname, sizeof(record->varname)-1);
683 record->varname[sizeof(record->varname) - 1] = 0;
684
685 grid_t *grid = (grid_t *)Malloc(sizeof(*grid));
686 const bool uvRelativeToGrid = gribapiGetGrid(gh, grid);
687
688 struct addIfNewRes gridAdded = cdiVlistAddGridIfNew(vlistID, grid, 0);
689 const int gridID = gridAdded.Id;
690 if ( !gridAdded.isNew ) Free(grid);
691 else if ( grid->projtype == CDI_PROJ_RLL ) gribapiDefProjRLL(gh, gridID);
692 else if ( grid->projtype == CDI_PROJ_LCC ) gribapiDefProjLCC(gh, gridID);
693 else if ( grid->projtype == CDI_PROJ_STERE ) gribapiDefProjSTERE(gh, gridID);
694
695 const int zaxistype = gribapiGetZaxisType(gribEditionNumber(gh), leveltype1);
696
697 switch (zaxistype)
698 {
699 case ZAXIS_HYBRID:
700 case ZAXIS_HYBRID_HALF:
701 {
702 long lpar;
703 GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
704 /* FIXME: assert(lpar >= 0) */
705 const size_t vctsize = (size_t)lpar;
706 if ( vctsize > 0 )
707 {
708 double *vctptr = (double *) Malloc(vctsize*sizeof(double));
709 size_t dummy = vctsize;
710 GRIB_CHECK(grib_get_double_array(gh, "pv", vctptr, &dummy), 0);
711 varDefVCT(vctsize, vctptr);
712 Free(vctptr);
713 }
714 break;
715 }
716 case ZAXIS_REFERENCE:
717 {
718 unsigned char uuid[CDI_UUID_SIZE];
719 long lpar;
720 GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
721 if ( lpar != 6 ) fprintf(stderr, "Warning ...\n");
722 GRIB_CHECK(grib_get_long(gh, "nlev", &lpar), 0);
723 int nhlev = (int)lpar;
724 GRIB_CHECK(grib_get_long(gh, "numberOfVGridUsed", &lpar), 0);
725 int nvgrid = (int)lpar;
726 size_t len = (size_t)CDI_UUID_SIZE;
727 memset(uuid, 0, CDI_UUID_SIZE);
728 GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", uuid, &len), 0);
729 varDefZAxisReference(nhlev, nvgrid, uuid);
730 break;
731 }
732 }
733
734 // if ( datatype > 32 ) datatype = CDI_DATATYPE_PACK32;
735 if ( datatype < 0 ) datatype = CDI_DATATYPE_PACK;
736
737 // add the previously read record data to the (intermediate) list of records
738 int tile_index = 0;
739 int varID = 0, levelID = 0;
740 varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, level_sf, level_unit,
741 datatype, &varID, &levelID, tsteptype, leveltype1, leveltype2,
742 varname, scanKeys, tiles, &tile_index);
743
744 record->varID = (short)varID;
745 record->levelID = levelID;
746
747 varDefCompType(varID, comptype);
748
749 if ( uvRelativeToGrid ) varDefKeyInt(varID, CDI_KEY_UVRELATIVETOGRID, 1);
750
751 if (varname[0]) gribapiGetNameKeys(gh, varID);
752 gribapiGetKeys(gh, varID);
753
754 if (lread_additional_keys)
755 {
756 long lval;
757 double dval;
758 for ( int i = 0; i < cdiNAdditionalGRIBKeys; i++ )
759 {
760 // note: if the key is not defined, we do not throw an error!
761 if ( grib_get_long(gh, cdiAdditionalGRIBKeys[i], &lval) == 0 )
762 varDefOptGribInt(varID, tile_index, lval, cdiAdditionalGRIBKeys[i]);
763 if ( grib_get_double(gh, cdiAdditionalGRIBKeys[i], &dval) == 0 )
764 varDefOptGribDbl(varID, tile_index, dval, cdiAdditionalGRIBKeys[i]);
765 }
766 }
767
768 if ( varInqInst(varID) == CDI_UNDEFID )
769 {
770 long center, subcenter;
771 GRIB_CHECK(grib_get_long(gh, "centre", ¢er), 0);
772 GRIB_CHECK(grib_get_long(gh, "subCentre", &subcenter), 0);
773 int instID = institutInq((int)center, (int)subcenter, NULL, NULL);
774 if ( instID == CDI_UNDEFID )
775 instID = institutDef((int)center, (int)subcenter, NULL, NULL);
776 varDefInst(varID, instID);
777 }
778
779 if ( varInqModel(varID) == CDI_UNDEFID )
780 {
781 long processID;
782 if ( grib_get_long(gh, "generatingProcessIdentifier", &processID) == 0 )
783 {
784 /* FIXME: assert(processID >= INT_MIN && processID <= INT_MAX) */
785 int modelID = modelInq(varInqInst(varID), (int)processID, NULL);
786 if ( modelID == CDI_UNDEFID )
787 modelID = modelDef(varInqInst(varID), (int)processID, NULL);
788 varDefModel(varID, modelID);
789 }
790 }
791
792 if ( varInqTable(varID) == CDI_UNDEFID )
793 {
794 int pdis, pcat, pnum;
795 cdiDecodeParam(param, &pnum, &pcat, &pdis);
796
797 if ( pdis == 255 )
798 {
799 int tabnum = pcat;
800 int tableID = tableInq(varInqModel(varID), tabnum, NULL);
801 if ( tableID == CDI_UNDEFID )
802 tableID = tableDef(varInqModel(varID), tabnum, NULL);
803 varDefTable(varID, tableID);
804 }
805 }
806
807 streamptr->tsteps[tsID].nallrecs++;
808 streamptr->nrecs++;
809
810 if ( CDI_Debug )
811 Message("varID = %d param = %d zaxistype = %d gridID = %d levelID = %d",
812 varID, param, zaxistype, gridID, levelID);
813 }
814
815 static
gribapiVarSet(int param,int level1,int level2,int leveltype,int tsteptype,size_t gridsize,char * name,VarScanKeys scanKeys,var_tile_t tiles_data)816 compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, int tsteptype,
817 size_t gridsize, char *name, VarScanKeys scanKeys, var_tile_t tiles_data)
818 {
819 compvar2_t compVar;
820 memset(&compVar, 0, sizeof(compvar2_t));
821 const size_t maxlen = sizeof(compVar.name);
822 size_t len = strlen(name);
823 if ( len > maxlen ) len = maxlen;
824
825 compVar.param = param;
826 compVar.level1 = level1;
827 compVar.level2 = level2;
828 compVar.ltype = leveltype;
829 compVar.tsteptype = tsteptype;
830 compVar.gridsize = gridsize;
831 //memset(compVar.name, 0, maxlen);
832 memcpy(compVar.name, name, len);
833 compVar.scanKeys = scanKeys;
834 compVar.tiles = tiles_data;
835
836 return compVar;
837 }
838
839 static
gribapiVarCompare(compvar2_t compVar,record_t record,int flag)840 int gribapiVarCompare(compvar2_t compVar, record_t record, int flag)
841 {
842 compvar2_t compVar0;
843 memset(&compVar0, 0, sizeof(compvar2_t));
844 compVar0.param = record.param;
845 compVar0.level1 = record.ilevel;
846 compVar0.level2 = record.ilevel2;
847 compVar0.ltype = record.ltype;
848 compVar0.tsteptype = record.tsteptype;
849 compVar0.gridsize = record.gridsize;
850 memcpy(compVar0.name, record.varname, sizeof(compVar.name));
851
852 if ( flag == 0 )
853 {
854 if ( compVar0.tsteptype == TSTEP_INSTANT && compVar.tsteptype == TSTEP_INSTANT3 ) compVar0.tsteptype = TSTEP_INSTANT3;
855 if ( compVar0.tsteptype == TSTEP_INSTANT3 && compVar.tsteptype == TSTEP_INSTANT ) compVar0.tsteptype = TSTEP_INSTANT;
856 }
857
858 compVar0.scanKeys = record.scanKeys;
859 compVar0.tiles = record.tiles;
860
861 return memcmp(&compVar0, &compVar, sizeof(compvar2_t));
862 }
863
864 static
gribapiGetDiskRepresentation(size_t recsize,size_t * buffersize,void ** gribbuffer,int * outDatatype,int * outCompressionType)865 grib_handle *gribapiGetDiskRepresentation(size_t recsize, size_t *buffersize, void **gribbuffer, int *outDatatype, int *outCompressionType)
866 {
867 const int gribversion = (int)((char*)*gribbuffer)[7];
868
869 if ( gribversion <= 1 ) *outCompressionType = grbDecompress(recsize, buffersize, gribbuffer);
870
871 grib_handle *gh = grib_handle_new_from_message(NULL, *gribbuffer, recsize);
872
873 bool lieee = false;
874
875 if ( gribversion > 1 )
876 {
877 size_t len = 256;
878 char typeOfPacking[256];
879 if ( grib_get_string(gh, "packingType", typeOfPacking, &len) == 0 )
880 {
881 // fprintf(stderr, "packingType %zu %s\n", len, typeOfPacking);
882 if ( strncmp(typeOfPacking, "grid_jpeg", len) == 0 ) *outCompressionType = CDI_COMPRESS_JPEG;
883 else if ( strncmp(typeOfPacking, "grid_ccsds", len) == 0 ) *outCompressionType = CDI_COMPRESS_AEC;
884 else if ( strncmp(typeOfPacking, "grid_ieee", len) == 0 ) lieee = true;
885 }
886 }
887
888 if ( lieee )
889 {
890 long precision;
891 const int status = grib_get_long(gh, "precision", &precision);
892 *outDatatype = (status == 0 && precision == 1) ? CDI_DATATYPE_FLT32 : CDI_DATATYPE_FLT64;
893 }
894 else
895 {
896 *outDatatype = CDI_DATATYPE_PACK;
897 long bitsPerValue;
898 if ( grib_get_long(gh, "bitsPerValue", &bitsPerValue) == 0 )
899 {
900 if ( bitsPerValue > 0 && bitsPerValue <= 32 ) *outDatatype = (int)bitsPerValue;
901 }
902 }
903
904 return gh;
905 }
906
907 typedef enum { CHECKTIME_OK, CHECKTIME_SKIP, CHECKTIME_STOP, CHECKTIME_INCONSISTENT } checkTimeResult;
908
checkTime(stream_t * streamptr,compvar2_t compVar,const DateTime * verificationTime,const DateTime * expectedVTime)909 static checkTimeResult checkTime(stream_t* streamptr, compvar2_t compVar, const DateTime* verificationTime, const DateTime* expectedVTime) {
910 // First determine whether the current record exists already.
911 int recID = 0;
912 for ( ; recID < streamptr->nrecs; recID++ )
913 {
914 if ( gribapiVarCompare(compVar, streamptr->tsteps[0].records[recID], 1) == 0 ) break;
915 }
916 const bool recordExists = recID < streamptr->nrecs;
917
918 // Then we need to know whether the verification time is consistent.
919 const bool consistentTime = !memcmp(verificationTime, expectedVTime, sizeof(*verificationTime));
920
921 // Finally, we make a decision.
922 if ( CDI_Inventory_Mode == 1 )
923 {
924 if ( recordExists ) return CHECKTIME_STOP;
925 if ( !consistentTime ) return CHECKTIME_INCONSISTENT;
926 }
927 else
928 {
929 if ( !consistentTime ) return CHECKTIME_STOP;
930 if ( recordExists ) return CHECKTIME_SKIP;
931 }
932
933 return CHECKTIME_OK;
934 }
935
936 #define gribWarning(text, nrecs, timestep, varname, param, level1, level2) do \
937 { \
938 char paramstr[32]; \
939 cdiParamToString(param, paramstr, sizeof(paramstr)); \
940 Warning("Record %2d (name=%s id=%s lev1=%d lev2=%d) timestep %d: %s", nrecs, varname, paramstr, level1, level2, timestep, text); \
941 } \
942 while(0)
943
944
945 #ifdef HAVE_LIBFDB5
946 #include "cdi_fdb.h"
947 #endif
948
fdbScanTimesteps(stream_t * streamptr)949 int fdbScanTimesteps(stream_t *streamptr)
950 {
951 #ifdef HAVE_LIBFDB5
952 void *gribbuffer = NULL;
953 size_t buffersize = 0;
954 grib_handle *gh = NULL;
955
956 fdb_handle_t *fdbHandle = streamptr->fdbHandle;
957
958 char **fdbItemList = NULL;
959 fdb_request_t *request = create_fdb_request(streamptr->filename);
960 int numItems = fdb_fill_itemlist(fdbHandle, request, &fdbItemList);
961 fdb_delete_request(request);
962 if (numItems == 0) Error("FDB request does not find any database entries!");
963 if (CDI_Debug) for (int i = 0; i < numItems; ++i) printf("item[%d] = %s\n", i + 1, fdbItemList[i]);
964
965 KeyValueEntry *keyValueList = (KeyValueEntry *) malloc(numItems * sizeof(KeyValueEntry));
966 for (int i = 0; i < numItems; ++i) keyValueList[i].item = NULL;
967
968 RecordInfoEntry *recordInfoList = (RecordInfoEntry *) malloc(numItems * sizeof(RecordInfoEntry));
969 for (int i = 0; i < numItems; ++i) record_info_entry_init(&recordInfoList[i]);
970
971 for (int i = 0; i < numItems; ++i) decode_fdbitem(fdbItemList[i], &keyValueList[i]);
972 if (check_keyvalueList(numItems, keyValueList) != 0) return 1;
973 for (int i = 0; i < numItems; ++i) decode_keyvalue(&keyValueList[i], &recordInfoList[i]);
974
975 int numRecords = get_num_records(numItems, recordInfoList);
976 int numTimesteps = numItems / numRecords;
977 if (CDI_Debug) Message("numRecords=%d numTimesteps=%d", numRecords, numTimesteps);
978
979 for (int i = 0; i < numItems; ++i) if (keyValueList[i].item) free(keyValueList[i].item);
980 free(keyValueList);
981
982 if (numRecords == 0) return CDI_EUFSTRUCT;
983
984 int *timestepRecordOffset = (int *) malloc(numTimesteps * sizeof(int));
985 for (int i = 0; i < numTimesteps; i++) timestepRecordOffset[i] = i * numRecords;
986 numTimesteps = remove_duplicate_timesteps(recordInfoList, numRecords, numTimesteps, timestepRecordOffset);
987 if (CDI_Debug) Message("numRecords=%d numTimesteps=%d", numRecords, numTimesteps);
988
989 DateTime datetime0 = { .date = 10101, .time = 0 };
990 int fcast = 0;
991
992 streamptr->curTsID = 0;
993
994 const int tsIDnew = tstepsNewEntry(streamptr);
995 if (tsIDnew != 0) Error("Internal problem! tstepsNewEntry returns %d", tsIDnew);
996
997 taxis_t *taxis = &streamptr->tsteps[tsIDnew].taxis;
998
999 for (int recID = 0; recID < numRecords; recID++)
1000 {
1001 const long recsize = fdb_read_record(fdbHandle, fdbItemList[recID], &buffersize, &gribbuffer);
1002
1003 int datatype, comptype = 0;
1004 gh = gribapiGetDiskRepresentation(recsize, &buffersize, &gribbuffer, &datatype, &comptype);
1005
1006 GRIB_CHECK(my_grib_set_double(gh, "missingValue", CDI_Default_Missval), 0);
1007
1008 const int param = gribapiGetParam(gh);
1009 int leveltype1 = -1, leveltype2 = -1, lbounds, level_sf, level_unit;
1010 var_tile_t tiles = dummy_tiles;
1011 int level1 = 0, level2 = 0;
1012 gribGetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles);
1013
1014 char varname[256];
1015 gribapiGetString(gh, "shortName", varname, sizeof(varname));
1016
1017 int64_t vdate, sdate;
1018 int vtime, stime;
1019 gribapiGetValidityDateTime(gh, &vdate, &vtime, &sdate, &stime);
1020 DateTime datetime = { .date = vdate, .time = vtime };
1021
1022 VarScanKeys scanKeys = gribapiGetScanKeys(gh);
1023
1024 if (recID == 0)
1025 {
1026 datetime0 = datetime;
1027 gribapiGetDataDateTime(gh, &(taxis->rdate), &(taxis->rtime));
1028 fcast = gribapiTimeIsFC(gh);
1029 if (fcast) taxis->unit = gribapiGetTimeUnits(gh);
1030 taxis->fdate = taxis->rdate;
1031 taxis->ftime = taxis->rtime;
1032 taxis->sdate = sdate;
1033 taxis->stime = stime;
1034 taxis->vdate = vdate;
1035 taxis->vtime = vtime;
1036 }
1037
1038 if (CDI_Debug)
1039 {
1040 char paramstr[32];
1041 cdiParamToString(param, paramstr, sizeof(paramstr));
1042 Message("%4d name=%s id=%s ltype=%d lev1=%d lev2=%d vdate=%lld vtime=%d",
1043 recID + 1, varname, paramstr, leveltype1, level1, level2, vdate, vtime);
1044 }
1045
1046 var_tile_t *ptiles = memcmp(&tiles, &dummy_tiles, sizeof(var_tile_t)) ? &tiles : NULL;
1047 int recpos = 0;
1048 gribapiAddRecord(streamptr, param, gh, recsize, recpos, datatype, comptype, varname, leveltype1, leveltype2,
1049 lbounds, level1, level2, level_sf, level_unit, &scanKeys, ptiles, 1, fdbItemList[recID]);
1050 fdbItemList[recID] = NULL;
1051
1052 grib_handle_delete(gh);
1053 gh = NULL;
1054 }
1055
1056 if (gh) grib_handle_delete(gh);
1057
1058 cdi_generate_vars(streamptr);
1059
1060 taxis->type = fcast ? TAXIS_RELATIVE : TAXIS_ABSOLUTE;
1061 int taxisID = taxisCreate(taxis->type);
1062
1063 vlistDefTaxis(streamptr->vlistID, taxisID);
1064
1065 streamScanResizeRecords1(streamptr);
1066
1067 streamptr->record->buffer = gribbuffer;
1068 streamptr->record->buffersize = buffersize;
1069
1070 if (numTimesteps == 1) streamptr->ntsteps = 1;
1071 streamScanTimeConstAdjust(streamptr, taxis);
1072
1073 for (int tsID = 1; tsID < numTimesteps; tsID++)
1074 {
1075 const int recordOffset = timestepRecordOffset[tsID];
1076 const int vdate = recordInfoList[recordOffset].date;
1077 const int vtime = recordInfoList[recordOffset].time * 100;
1078 // printf("timestep=%d recOffset=%d date=%d time=%d\n", tsID + 1, recordOffset, vdate, vtime);
1079
1080 const int tsIDnext = tstepsNewEntry(streamptr);
1081 if (tsIDnext != tsID) Error("Internal error. tsID = %d", tsID);
1082
1083 streamptr->tsteps[tsID-1].next = true;
1084 streamptr->tsteps[tsID].position = 0;
1085
1086 taxis = &streamptr->tsteps[tsID].taxis;
1087
1088 cdi_create_records(streamptr, tsID);
1089 record_t *records = streamptr->tsteps[tsID].records;
1090
1091 const int nrecs = (tsID == 1) ? streamScanInitRecords2(streamptr) : streamScanInitRecords(streamptr, tsID);
1092 if (nrecs != numRecords) Error("Internal error. nrecs = %d", nrecs);
1093
1094 taxis->vdate = vdate;
1095 taxis->vtime = vtime;
1096
1097 int rindex = 0;
1098 for (int recID = 0; recID < numRecords; recID++)
1099 {
1100 records[recID].used = true;
1101 streamptr->tsteps[tsID].recIDs[rindex] = recID;
1102 rindex++;
1103
1104 records[recID].position = 0;
1105 records[recID].size = 0;
1106 records[recID].fdbItem = fdbItemList[recordOffset + recID];
1107 fdbItemList[recordOffset + recID] = NULL;
1108 }
1109
1110 if (tsID == 1) streamptr->tsteps[1].nrecs = numRecords;
1111 }
1112
1113 streamptr->rtsteps = numTimesteps;
1114 streamptr->ntsteps = numTimesteps;
1115
1116 for (int i = 0; i < numItems; ++i)
1117 if (fdbItemList[i]) free(fdbItemList[i]);
1118
1119 if (fdbItemList) free(fdbItemList);
1120 if (recordInfoList) free(recordInfoList);
1121 if (timestepRecordOffset) free(timestepRecordOffset);
1122 #endif
1123
1124 return 0;
1125 }
1126
1127
gribapiScanTimestep1(stream_t * streamptr)1128 int gribapiScanTimestep1(stream_t *streamptr)
1129 {
1130 off_t recpos = 0;
1131 void *gribbuffer = NULL;
1132 size_t buffersize = 0;
1133 DateTime datetime0 = { .date = 10101, .time = 0 };
1134 int nrecs_scanned = 0; // Only used for debug output.
1135 bool warn_time = true;
1136 int fcast = 0;
1137 grib_handle *gh = NULL;
1138
1139 streamptr->curTsID = 0;
1140
1141 const int tsID = tstepsNewEntry(streamptr);
1142 if (tsID != 0) Error("Internal problem! tstepsNewEntry returns %d", tsID);
1143
1144 taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
1145
1146 const int fileID = streamptr->fileID;
1147
1148 unsigned nrecs = 0;
1149 while (true)
1150 {
1151 const size_t recsize = gribGetSize(fileID);
1152 recpos = fileGetPos(fileID);
1153
1154 if (recsize == 0)
1155 {
1156 streamptr->ntsteps = 1;
1157 break;
1158 }
1159
1160 ensureBufferSize(recsize, &buffersize, &gribbuffer);
1161
1162 size_t readsize = recsize;
1163 // Search for next 'GRIB', read the following record, and position file offset after it.
1164 if (gribRead(fileID, gribbuffer, &readsize)) break;
1165
1166 int datatype, comptype = 0;
1167 gh = gribapiGetDiskRepresentation(recsize, &buffersize, &gribbuffer, &datatype, &comptype);
1168
1169 nrecs_scanned++;
1170 GRIB_CHECK(my_grib_set_double(gh, "missingValue", CDI_Default_Missval), 0);
1171
1172 const int param = gribapiGetParam(gh);
1173 int leveltype1 = -1, leveltype2 = -1, lbounds, level_sf, level_unit;
1174 var_tile_t tiles = dummy_tiles;
1175 int level1 = 0, level2 = 0;
1176 gribGetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles);
1177
1178 char varname[256];
1179 gribapiGetString(gh, "shortName", varname, sizeof(varname));
1180
1181 int64_t vdate, sdate;
1182 int vtime, stime;
1183 gribapiGetValidityDateTime(gh, &vdate, &vtime, &sdate, &stime);
1184 DateTime datetime = { .date = vdate, .time = vtime };
1185
1186 VarScanKeys scanKeys = gribapiGetScanKeys(gh);
1187
1188 if (nrecs == 0)
1189 {
1190 datetime0 = datetime;
1191 gribapiGetDataDateTime(gh, &(taxis->rdate), &(taxis->rtime));
1192 fcast = gribapiTimeIsFC(gh);
1193 if (fcast) taxis->unit = gribapiGetTimeUnits(gh);
1194 taxis->fdate = taxis->rdate;
1195 taxis->ftime = taxis->rtime;
1196 taxis->sdate = sdate;
1197 taxis->stime = stime;
1198 taxis->vdate = vdate;
1199 taxis->vtime = vtime;
1200 }
1201 else
1202 {
1203 const int tsteptype = gribapiGetTsteptype(gh);
1204 const size_t gridsize = gribapiGetGridsize(gh);
1205 checkTimeResult result = checkTime(streamptr,
1206 gribapiVarSet(param, level1, level2, leveltype1, tsteptype, gridsize, varname, scanKeys, tiles),
1207 &datetime, &datetime0);
1208 if (result == CHECKTIME_STOP)
1209 {
1210 break;
1211 }
1212 else if (result == CHECKTIME_SKIP)
1213 {
1214 gribWarning("Parameter already exist, skipped!", nrecs_scanned, tsID+1, varname, param, level1, level2);
1215 continue;
1216 }
1217 else if (result == CHECKTIME_INCONSISTENT && warn_time)
1218 {
1219 gribWarning("Inconsistent verification time!", nrecs_scanned, tsID+1, varname, param, level1, level2);
1220 warn_time = false;
1221 }
1222 assert(result == CHECKTIME_OK || result == CHECKTIME_INCONSISTENT);
1223 }
1224
1225 nrecs++;
1226
1227 if (CDI_Debug)
1228 {
1229 char paramstr[32];
1230 cdiParamToString(param, paramstr, sizeof(paramstr));
1231 Message("%4u %8d name=%s id=%s ltype=%d lev1=%d lev2=%d vdate=%lld vtime=%d",
1232 nrecs, (int)recpos, varname, paramstr, leveltype1, level1, level2, vdate, vtime);
1233 }
1234
1235 var_tile_t *ptiles = memcmp(&tiles, &dummy_tiles, sizeof(var_tile_t)) ? &tiles : NULL;
1236 gribapiAddRecord(streamptr, param, gh, recsize, recpos, datatype, comptype, varname, leveltype1, leveltype2,
1237 lbounds, level1, level2, level_sf, level_unit, &scanKeys, ptiles, 1, NULL);
1238
1239 grib_handle_delete(gh);
1240 gh = NULL;
1241 }
1242
1243 if (gh) grib_handle_delete(gh);
1244
1245 streamptr->rtsteps = 1;
1246
1247 if (nrecs == 0) return CDI_EUFSTRUCT;
1248
1249 cdi_generate_vars(streamptr);
1250
1251 taxis->type = fcast ? TAXIS_RELATIVE : TAXIS_ABSOLUTE;
1252 int taxisID = taxisCreate(taxis->type);
1253 // printf("1: %d %6d %d %6d %d %6d\n", taxis->rdate, taxis->rtime, taxis->sdate, taxis->stime, taxis->vdate, taxis->vtime);
1254
1255 const int vlistID = streamptr->vlistID;
1256 vlistDefTaxis(vlistID, taxisID);
1257
1258 streamScanResizeRecords1(streamptr);
1259
1260 streamptr->record->buffer = gribbuffer;
1261 streamptr->record->buffersize = buffersize;
1262
1263 streamScanTsFixNtsteps(streamptr, recpos);
1264 streamScanTimeConstAdjust(streamptr, taxis);
1265
1266 return 0;
1267 }
1268
1269
gribapiScanTimestep2(stream_t * streamptr)1270 int gribapiScanTimestep2(stream_t *streamptr)
1271 {
1272 int rstatus = 0;
1273 off_t recpos = 0;
1274 DateTime datetime0 = { LONG_MIN, LONG_MIN };
1275 // int gridID;
1276 int recID;
1277 grib_handle *gh = NULL;
1278
1279 streamptr->curTsID = 1;
1280
1281 const int fileID = streamptr->fileID;
1282 const int vlistID = streamptr->vlistID;
1283
1284 void *gribbuffer = streamptr->record->buffer;
1285 size_t buffersize = streamptr->record->buffersize;
1286
1287 int tsID = streamptr->rtsteps;
1288 if (tsID != 1) Error("Internal problem! unexpected timestep %d", tsID+1);
1289
1290 taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
1291
1292 fileSetPos(fileID, streamptr->tsteps[tsID].position, SEEK_SET);
1293
1294 cdi_create_records(streamptr, tsID);
1295 record_t *records = streamptr->tsteps[tsID].records;
1296
1297 const int nrecords = streamScanInitRecords2(streamptr);
1298
1299 int nrecs_scanned = nrecords; //Only used for debug output
1300 int rindex = 0;
1301 while (true)
1302 {
1303 if (rindex > nrecords) break;
1304
1305 const size_t recsize = gribGetSize(fileID);
1306 recpos = fileGetPos(fileID);
1307 if (recsize == 0)
1308 {
1309 streamptr->ntsteps = 2;
1310 break;
1311 }
1312
1313 ensureBufferSize(recsize, &buffersize, &gribbuffer);
1314
1315 size_t readsize = recsize;
1316 if (gribRead(fileID, gribbuffer, &readsize)) break;
1317
1318 grbDecompress(recsize, &buffersize, &gribbuffer);
1319
1320 nrecs_scanned++;
1321 gh = grib_handle_new_from_message(NULL, gribbuffer, recsize);
1322 GRIB_CHECK(my_grib_set_double(gh, "missingValue", CDI_Default_Missval), 0);
1323
1324 const int param = gribapiGetParam(gh);
1325 int level1 = 0, level2 = 0, leveltype1, leveltype2, lbounds, level_sf, level_unit;
1326 var_tile_t tiles = dummy_tiles;
1327 gribGetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles);
1328
1329 char varname[256];
1330 gribapiGetString(gh, "shortName", varname, sizeof(varname));
1331
1332 int64_t vdate, sdate;
1333 int vtime, stime;
1334 gribapiGetValidityDateTime(gh, &vdate, &vtime, &sdate, &stime);
1335 DateTime datetime = { .date = vdate, .time = vtime };
1336
1337 if (rindex == 0)
1338 {
1339 datetime0 = datetime;
1340 const int taxisID = vlistInqTaxis(vlistID);
1341 if (taxisInqType(taxisID) == TAXIS_RELATIVE)
1342 {
1343 taxis->type = TAXIS_RELATIVE;
1344 taxis->unit = gribapiGetTimeUnits(gh);
1345 gribapiGetDataDateTime(gh, &(taxis->rdate), &(taxis->rtime));
1346 }
1347 else
1348 {
1349 taxis->type = TAXIS_ABSOLUTE;
1350 }
1351 taxis->fdate = taxis->rdate;
1352 taxis->ftime = taxis->rtime;
1353 taxis->vdate = vdate;
1354 taxis->vtime = vtime;
1355 taxis->sdate = sdate;
1356 taxis->stime = stime;
1357 // printf("2: %d %6d %d %6d %d %6d\n", taxis->rdate, taxis->rtime, taxis->sdate, taxis->stime, taxis->vdate, taxis->vtime);
1358 }
1359
1360 VarScanKeys scanKeys = gribapiGetScanKeys(gh);
1361
1362 const int tsteptype = gribapiGetTsteptype(gh);
1363 const size_t gridsize = gribapiGetGridsize(gh);
1364 compvar2_t compVar = gribapiVarSet(param, level1, level2, leveltype1, tsteptype, gridsize, varname, scanKeys, tiles);
1365
1366 for (recID = 0; recID < nrecords; recID++)
1367 if (gribapiVarCompare(compVar, records[recID], 0) == 0) break;
1368
1369 if (recID == nrecords)
1370 {
1371 if (CDI_Inventory_Mode == 1)
1372 {
1373 gribWarning("Parameter not defined at timestep 1!", nrecs_scanned, tsID+1, varname, param, level1, level2);
1374 return CDI_EUFSTRUCT;
1375 }
1376 else
1377 {
1378 gribWarning("Parameter not defined at timestep 1, skipped!", nrecs_scanned, tsID+1, varname, param, level1, level2);
1379 continue;
1380 }
1381 }
1382
1383 if (records[recID].used)
1384 {
1385 if (CDI_Inventory_Mode == 1) break;
1386 else
1387 {
1388 if (datetimeDiffer(datetime, datetime0)) break;
1389
1390 gribWarning("Parameter already exist, skipped!", nrecs_scanned, tsID+1, varname, param, level1, level2);
1391 continue;
1392 }
1393 }
1394
1395 records[recID].used = true;
1396 streamptr->tsteps[tsID].recIDs[rindex] = recID;
1397
1398 if (CDI_Debug)
1399 {
1400 char paramstr[32];
1401 cdiParamToString(param, paramstr, sizeof(paramstr));
1402 Message("%4d %8d name=%s id=%s ltype=%d lev1=%d lev2=%d vdate=%lld vtime=%d",
1403 nrecs_scanned, (int)recpos, varname, paramstr, leveltype1, level1, level2, vdate, vtime);
1404 }
1405
1406 if (gribapiVarCompare(compVar, records[recID], 0) != 0)
1407 {
1408 Message("tsID = %d recID = %d param = %3d new %3d level = %3d new %3d",
1409 tsID, recID, records[recID].param, param, records[recID].ilevel, level1);
1410 return CDI_EUFSTRUCT;
1411 }
1412
1413 records[recID].position = recpos;
1414 records[recID].size = recsize;
1415
1416 const int varID = records[recID].varID;
1417
1418 if (tsteptype != vlistInqVarTsteptype(vlistID, varID))
1419 vlistDefVarTsteptype(vlistID, varID, tsteptype);
1420
1421 grib_handle_delete(gh);
1422 gh = NULL;
1423
1424 rindex++;
1425 }
1426
1427 if (gh) grib_handle_delete(gh);
1428
1429 int nrecs = 0;
1430 for (recID = 0; recID < nrecords; recID++)
1431 {
1432 if (! records[recID].used)
1433 {
1434 vlistDefVarTimetype(vlistID, records[recID].varID, TIME_CONSTANT);
1435 }
1436 else
1437 {
1438 nrecs++;
1439 }
1440 }
1441 streamptr->tsteps[tsID].nrecs = nrecs;
1442
1443 streamptr->rtsteps = 2;
1444
1445 streamScanTsFixNtsteps(streamptr, recpos);
1446
1447 streamptr->record->buffer = gribbuffer;
1448 streamptr->record->buffersize = buffersize;
1449
1450 return rstatus;
1451 }
1452
1453
gribapiScanTimestep(stream_t * streamptr)1454 int gribapiScanTimestep(stream_t *streamptr)
1455 {
1456 int vrecID, recID = -1;
1457 int nrecs = 0;
1458 int vlistID = streamptr->vlistID;
1459
1460 int tsID = streamptr->rtsteps;
1461 taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
1462
1463 if (streamptr->tsteps[tsID].recordSize == 0)
1464 {
1465 void *gribbuffer = streamptr->record->buffer;
1466 size_t buffersize = streamptr->record->buffersize;
1467
1468 cdi_create_records(streamptr, tsID);
1469 record_t *records = streamptr->tsteps[tsID].records;
1470
1471 nrecs = streamScanInitRecords(streamptr, tsID);
1472
1473 const int fileID = streamptr->fileID;
1474
1475 fileSetPos(fileID, streamptr->tsteps[tsID].position, SEEK_SET);
1476
1477 int nrecs_scanned = streamptr->tsteps[0].nallrecs + streamptr->tsteps[1].nrecs*(tsID-1); // Only used for debug output.
1478 int rindex = 0;
1479 off_t recpos = 0;
1480 DateTime datetime0 = { LONG_MIN, LONG_MIN };
1481 grib_handle *gh = NULL;
1482 char varname[256];
1483 while (true)
1484 {
1485 if (rindex > nrecs) break;
1486
1487 const size_t recsize = gribGetSize(fileID);
1488 recpos = fileGetPos(fileID);
1489 if (recsize == 0)
1490 {
1491 streamptr->ntsteps = streamptr->rtsteps + 1;
1492 break;
1493 }
1494
1495 if (rindex >= nrecs) break;
1496
1497 ensureBufferSize(recsize, &buffersize, &gribbuffer);
1498
1499 size_t readsize = recsize;
1500 if (gribRead(fileID, gribbuffer, &readsize))
1501 {
1502 Warning("Inconsistent timestep %d (GRIB record %d/%d)!", tsID+1, rindex+1,
1503 streamptr->tsteps[tsID].recordSize);
1504 break;
1505 }
1506
1507 grbDecompress(recsize, &buffersize, &gribbuffer);
1508
1509 nrecs_scanned++;
1510 gh = grib_handle_new_from_message(NULL, gribbuffer, recsize);
1511 GRIB_CHECK(my_grib_set_double(gh, "missingValue", CDI_Default_Missval), 0);
1512
1513 const int param = gribapiGetParam(gh);
1514 int level1 = 0, level2 = 0, leveltype1, leveltype2 = -1, lbounds, level_sf, level_unit;
1515 var_tile_t tiles = dummy_tiles;
1516 gribGetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles);
1517
1518 gribapiGetString(gh, "shortName", varname, sizeof(varname));
1519
1520 int64_t vdate, sdate;
1521 int vtime, stime;
1522 gribapiGetValidityDateTime(gh, &vdate, &vtime, &sdate, &stime);
1523 DateTime datetime = { .date = vdate, .time = vtime };
1524
1525 if (rindex == nrecs) break;
1526
1527 if (rindex == 0)
1528 {
1529 datetime0 = datetime;
1530 const int taxisID = vlistInqTaxis(vlistID);
1531 if (taxisInqType(taxisID) == TAXIS_RELATIVE)
1532 {
1533 taxis->type = TAXIS_RELATIVE;
1534 taxis->unit = gribapiGetTimeUnits(gh);
1535 gribapiGetDataDateTime(gh, &(taxis->rdate), &(taxis->rtime));
1536 }
1537 else
1538 {
1539 taxis->type = TAXIS_ABSOLUTE;
1540 }
1541 taxis->fdate = taxis->rdate;
1542 taxis->ftime = taxis->rtime;
1543 taxis->vdate = vdate;
1544 taxis->vtime = vtime;
1545 taxis->sdate = sdate;
1546 taxis->stime = stime;
1547 // printf("n: %d %6d %d %6d %d %6d\n", taxis->rdate, taxis->rtime, taxis->sdate, taxis->stime, taxis->vdate, taxis->vtime);
1548 }
1549
1550 VarScanKeys scanKeys = gribapiGetScanKeys(gh);
1551
1552 const int tsteptype = gribapiGetTsteptype(gh);
1553 const size_t gridsize = gribapiGetGridsize(gh);
1554 compvar2_t compVar = gribapiVarSet(param, level1, level2, leveltype1, tsteptype, gridsize, varname, scanKeys, tiles);
1555
1556 for (vrecID = 0; vrecID < nrecs; vrecID++)
1557 {
1558 recID = streamptr->tsteps[1].recIDs[vrecID];
1559 if (gribapiVarCompare(compVar, records[recID], 0) == 0) break;
1560 }
1561
1562 if (vrecID == nrecs)
1563 {
1564 if (CDI_Inventory_Mode == 1)
1565 {
1566 gribWarning("Parameter not defined at timestep 1!", nrecs_scanned, tsID+1, varname, param, level1, level2);
1567 return CDI_EUFSTRUCT;
1568 }
1569 else
1570 {
1571 gribWarning("Parameter not defined at timestep 1, skipped!", nrecs_scanned, tsID+1, varname, param, level1, level2);
1572 continue;
1573 }
1574 }
1575
1576 if (CDI_Inventory_Mode != 1)
1577 {
1578 if (records[recID].used)
1579 {
1580 if (datetimeDiffer(datetime, datetime0)) break;
1581
1582 if (CDI_Debug)
1583 gribWarning("Parameter already exist, skipped!", nrecs_scanned, tsID+1, varname, param, level1, level2);
1584
1585 continue;
1586 }
1587 }
1588
1589 records[recID].used = true;
1590 streamptr->tsteps[tsID].recIDs[rindex] = recID;
1591
1592 if (CDI_Debug)
1593 Message("%4d %8d %4d %8d %8lld %6d", rindex+1, (int)recpos, param, level1, vdate, vtime);
1594
1595 if (gribapiVarCompare(compVar, records[recID], 0) != 0)
1596 {
1597 Message("tsID = %d recID = %d param = %3d new %3d level = %3d new %3d",
1598 tsID, recID, records[recID].param, param, records[recID].ilevel, level1);
1599 Error("Invalid, unsupported or inconsistent record structure");
1600 }
1601
1602 records[recID].position = recpos;
1603 records[recID].size = recsize;
1604
1605 if (CDI_Debug)
1606 Message("%4d %8d %4d %8d %8lld %6d", rindex, (int)recpos, param, level1, vdate, vtime);
1607
1608 grib_handle_delete(gh);
1609 gh = NULL;
1610
1611 rindex++;
1612 }
1613
1614 if (gh) grib_handle_delete(gh);
1615
1616 for (vrecID = 0; vrecID < nrecs; vrecID++)
1617 {
1618 recID = streamptr->tsteps[tsID].recIDs[vrecID];
1619 if ( ! records[recID].used ) break;
1620 }
1621
1622 if (vrecID < nrecs)
1623 {
1624 gribWarning("Paramameter not found!", nrecs_scanned, tsID+1, varname, records[recID].param,
1625 records[recID].ilevel, records[recID].ilevel2);
1626 return CDI_EUFSTRUCT;
1627 }
1628
1629 streamptr->rtsteps++;
1630
1631 if (streamptr->ntsteps != streamptr->rtsteps)
1632 {
1633 tsID = tstepsNewEntry(streamptr);
1634 if (tsID != streamptr->rtsteps) Error("Internal error. tsID = %d", tsID);
1635
1636 streamptr->tsteps[tsID-1].next = true;
1637 streamptr->tsteps[tsID].position = recpos;
1638 }
1639
1640 fileSetPos(fileID, streamptr->tsteps[tsID].position, SEEK_SET);
1641 streamptr->tsteps[tsID].position = recpos;
1642
1643 streamptr->record->buffer = gribbuffer;
1644 streamptr->record->buffersize = buffersize;
1645 }
1646
1647 if (nrecs > 0 && nrecs < streamptr->tsteps[tsID].nrecs)
1648 {
1649 Warning("Incomplete timestep. Stop scanning at timestep %d.", tsID);
1650 streamptr->ntsteps = tsID;
1651 }
1652
1653 return streamptr->ntsteps;
1654 }
1655
1656 #ifdef gribWarning
1657 #undef gribWarning
1658 #endif
1659
gribapiDecode(void * gribbuffer,size_t gribsize,void * data,size_t gridsize,int unreduced,size_t * nmiss,double missval)1660 int gribapiDecode(void *gribbuffer, size_t gribsize, void *data, size_t gridsize,
1661 int unreduced, size_t *nmiss, double missval)
1662 {
1663 int status = 0;
1664
1665 if ( unreduced )
1666 {
1667 static bool lwarn = true;
1668 if ( lwarn )
1669 {
1670 lwarn = false;
1671 Warning("Conversion of gaussian reduced grids unsupported!");
1672 }
1673 }
1674
1675 size_t recsize = (size_t)gribsize;
1676 grib_handle *gh = grib_handle_new_from_message(NULL, gribbuffer, recsize);
1677 GRIB_CHECK(my_grib_set_double(gh, "missingValue", missval), 0);
1678
1679 // get the size of the values array
1680 size_t datasize;
1681 GRIB_CHECK(grib_get_size(gh, "values", &datasize), 0);
1682 // long numberOfPoints;
1683 // GRIB_CHECK(grib_get_long(gh, "numberOfPoints", &numberOfPoints), 0);
1684 // printf("values_size = %d numberOfPoints = %ld\n", datasize, numberOfPoints);
1685
1686 if ( gridsize != datasize )
1687 Error("Internal problem: gridsize(%zu) != datasize(%zu)!", gridsize, datasize);
1688 size_t dummy = datasize;
1689 GRIB_CHECK(grib_get_double_array(gh, "values", (double*)data, &dummy), 0);
1690
1691 long lpar;
1692 GRIB_CHECK(grib_get_long(gh, "gridDefinitionTemplateNumber", &lpar), 0);
1693 int gridtype = (int) lpar;
1694
1695 *nmiss = 0;
1696 if ( gridtype < 50 || gridtype > 53 )
1697 {
1698 GRIB_CHECK(grib_get_long(gh, "numberOfMissing", &lpar), 0);
1699 *nmiss = (int) lpar;
1700 // printf("gridtype %d, nmiss %d\n", gridtype, nmiss);
1701 }
1702
1703 grib_handle_delete(gh);
1704
1705 return status;
1706 }
1707
1708
1709 static
gribapiDefInstitut(grib_handle * gh,int vlistID,int varID)1710 void gribapiDefInstitut(grib_handle *gh, int vlistID, int varID)
1711 {
1712 int instID = (vlistInqInstitut(vlistID) != CDI_UNDEFID) ?
1713 vlistInqInstitut(vlistID) : vlistInqVarInstitut(vlistID, varID);
1714
1715 if ( instID != CDI_UNDEFID )
1716 {
1717 long center = institutInqCenter(instID);
1718 long subcenter = institutInqSubcenter(instID);
1719
1720 long center0, subcenter0;
1721 GRIB_CHECK(grib_get_long(gh, "centre", ¢er0), 0);
1722 GRIB_CHECK(grib_get_long(gh, "subCentre", &subcenter0), 0);
1723
1724 if ( center != center0 )
1725 GRIB_CHECK(my_grib_set_long(gh, "centre", center), 0);
1726 if ( subcenter != subcenter0 )
1727 GRIB_CHECK(my_grib_set_long(gh, "subCentre", subcenter), 0);
1728 }
1729
1730 int status;
1731 int centre, subCentre;
1732 status = cdiInqKeyInt(vlistID, CDI_GLOBAL, CDI_KEY_CENTRE, ¢re);
1733 if ( status == 0 ) grib_set_long(gh, "centre", centre);
1734 status = cdiInqKeyInt(vlistID, CDI_GLOBAL, CDI_KEY_SUBCENTRE, &subCentre);
1735 if ( status == 0 ) grib_set_long(gh, "subCentre", subCentre);
1736
1737 status = cdiInqKeyInt(vlistID, varID, CDI_KEY_CENTRE, ¢re);
1738 if ( status == 0 ) grib_set_long(gh, "centre", centre);
1739 status = cdiInqKeyInt(vlistID, varID, CDI_KEY_SUBCENTRE, &subCentre);
1740 if ( status == 0 ) grib_set_long(gh, "subCentre", subCentre);
1741 }
1742
1743 static
gribapiDefModel(grib_handle * gh,int vlistID,int varID)1744 void gribapiDefModel(grib_handle *gh, int vlistID, int varID)
1745 {
1746 int modelID = (vlistInqModel(vlistID) != CDI_UNDEFID) ?
1747 vlistInqModel(vlistID) : vlistInqVarModel(vlistID, varID);
1748
1749 if ( modelID != CDI_UNDEFID )
1750 GRIB_CHECK(my_grib_set_long(gh, "generatingProcessIdentifier", modelInqGribID(modelID)), 0);
1751 }
1752
1753 static
gribapiDefParam(int editionNumber,grib_handle * gh,int param,const char * name,const char * stdname)1754 void gribapiDefParam(int editionNumber, grib_handle *gh, int param, const char *name, const char *stdname)
1755 {
1756 bool ldefined = false;
1757
1758 int pdis, pcat, pnum;
1759 cdiDecodeParam(param, &pnum, &pcat, &pdis);
1760
1761 if ( pnum < 0 )
1762 {
1763 size_t len = strlen(stdname);
1764 if ( len )
1765 {
1766 int status = my_grib_set_string(gh, "cfName", stdname, &len);
1767 if ( status == 0 ) ldefined = true;
1768 else Warning("grib_api: No match for cfName=%s", stdname);
1769 }
1770
1771 if ( ldefined == false )
1772 {
1773 len = strlen(name);
1774 int status = my_grib_set_string(gh, "shortName", name, &len);
1775 if ( status == 0 ) ldefined = true;
1776 else Warning("grib_api: No match for shortName=%s", name);
1777 }
1778 }
1779
1780 if ( ldefined == false )
1781 {
1782 if ( pnum < 0 ) pnum = -pnum;
1783
1784 if ( pnum > 255 )
1785 {
1786 static bool lwarn_pnum = true;
1787 if ( lwarn_pnum )
1788 {
1789 Warning("Parameter number %d out of range (1-255), set to %d!", pnum, pnum%256);
1790 lwarn_pnum = false;
1791 }
1792 pnum = pnum%256;
1793 }
1794
1795 if ( editionNumber <= 1 )
1796 {
1797 static bool lwarn_pdis = true;
1798 if ( pdis != 255 && lwarn_pdis )
1799 {
1800 char paramstr[32];
1801 cdiParamToString(param, paramstr, sizeof(paramstr));
1802 Warning("Can't convert GRIB2 parameter ID (%s) to GRIB1, set to %d.%d!", paramstr, pnum, pcat);
1803 lwarn_pdis = false;
1804 }
1805
1806 GRIB_CHECK(my_grib_set_long(gh, "table2Version", pcat), 0);
1807 GRIB_CHECK(my_grib_set_long(gh, "indicatorOfParameter", pnum), 0);
1808 }
1809 else
1810 {
1811 GRIB_CHECK(my_grib_set_long(gh, "discipline", pdis), 0);
1812 GRIB_CHECK(my_grib_set_long(gh, "parameterCategory", pcat), 0);
1813 GRIB_CHECK(my_grib_set_long(gh, "parameterNumber", pnum), 0);
1814 }
1815 }
1816
1817 // printf("param: %d.%d.%d %s\n", pnum, pcat, pdis, name);
1818 }
1819
1820 static
getTimeunitFactor(const int timeunit)1821 int getTimeunitFactor(const int timeunit)
1822 {
1823 switch (timeunit)
1824 {
1825 case TUNIT_SECOND: return 1;
1826 case TUNIT_MINUTE: return 60;
1827 case TUNIT_HOUR: return 3600;
1828 case TUNIT_3HOURS: return 10800;
1829 case TUNIT_6HOURS: return 21600;
1830 case TUNIT_12HOURS: return 43200;
1831 case TUNIT_DAY: return 86400;
1832 }
1833
1834 return 3600;
1835 }
1836
1837 static
grib2ProDefTempHasStatisticalDef(const int proDefTempNum)1838 int grib2ProDefTempHasStatisticalDef(const int proDefTempNum)
1839 {
1840 switch (proDefTempNum)
1841 {
1842 case 8:
1843 case 9:
1844 case 10:
1845 case 11:
1846 case 12:
1847 case 13:
1848 case 14:
1849 case 34:
1850 case 42:
1851 case 43:
1852 case 46:
1853 case 47:
1854 case 61:
1855 case 91:
1856 case 1001:
1857 case 1101:
1858 case 40034:
1859 return 1;
1860 }
1861
1862 return 0;
1863 }
1864
1865 static
getUnitsOfTime(const int timeunit)1866 int getUnitsOfTime(const int timeunit)
1867 {
1868 switch (timeunit)
1869 {
1870 case TUNIT_SECOND: return 13;
1871 case TUNIT_MINUTE: return 0;
1872 case TUNIT_HOUR: return 1;
1873 case TUNIT_3HOURS: return 10;
1874 case TUNIT_6HOURS: return 11;
1875 case TUNIT_12HOURS: return 12;
1876 case TUNIT_DAY: return 2;
1877 }
1878
1879 return 1;
1880 }
1881
1882 static
gribapiDefStepUnits(int editionNumber,grib_handle * gh,int timeunit,int proDefTempNum,int gcinit)1883 void gribapiDefStepUnits(int editionNumber, grib_handle *gh, int timeunit, int proDefTempNum, int gcinit)
1884 {
1885 if ( !gcinit )
1886 {
1887 long unitsOfTime = getUnitsOfTime(timeunit);
1888
1889 grib_set_long(gh, "stepUnits", unitsOfTime);
1890 if ( editionNumber == 1 )
1891 {
1892 grib_set_long(gh, "unitOfTimeRange", unitsOfTime);
1893 }
1894 else if ( grib2ProDefTempHasStatisticalDef(proDefTempNum) )
1895 {
1896 grib_set_long(gh, "indicatorOfUnitForTimeRange", unitsOfTime);
1897 grib_set_long(gh, "indicatorOfUnitOfTimeRange", unitsOfTime);
1898 }
1899 else
1900 {
1901 // NOTE KNMI: HIRLAM model files LAMH_D11 are in grib1 and do NOT have key indicatorOfUnitForTimeRange
1902 // Watch out for compatibility issues.
1903 grib_set_long(gh, "indicatorOfUnitOfTimeRange", unitsOfTime);
1904 }
1905 }
1906 }
1907
1908 static
gribapiDefSteptype(int editionNumber,grib_handle * gh,int productDefinitionTemplate,int typeOfGeneratingProcess,int tsteptype,int gcinit)1909 int gribapiDefSteptype(int editionNumber, grib_handle *gh, int productDefinitionTemplate, int typeOfGeneratingProcess, int tsteptype, int gcinit)
1910 {
1911 const char *stepType = "instant";
1912 long proDefTempNum = 0;
1913
1914 if ( tsteptype >= TSTEP_INSTANT && tsteptype <= TSTEP_SUM )
1915 {
1916 stepType = cdiGribAPI_ts_str_map[tsteptype].sname;
1917 proDefTempNum = cdiGribAPI_ts_str_map[tsteptype].productionTemplate;
1918 }
1919
1920 if ( typeOfGeneratingProcess == 4 )
1921 {
1922 proDefTempNum = (proDefTempNum == 8) ? 11 : 1;
1923 }
1924
1925 if ( productDefinitionTemplate != -1 ) proDefTempNum = productDefinitionTemplate;
1926
1927 if ( !gcinit )
1928 {
1929 if ( editionNumber > 1 ) GRIB_CHECK(my_grib_set_long(gh, "productDefinitionTemplateNumber", proDefTempNum), 0);
1930 size_t len = strlen(stepType);
1931 int status = my_grib_set_string(gh, "stepType", stepType, &len);
1932 if (status != 0) GRIB_CHECK(my_grib_set_long(gh, "productDefinitionTemplateNumber", 0), 0);
1933 }
1934
1935 return (int)proDefTempNum;
1936 }
1937
1938 static
gribapiDefDateTimeAbs(int editionNumber,grib_handle * gh,int64_t date,int time,int productDefinitionTemplate,int typeOfGeneratingProcess,int tsteptype,int gcinit)1939 void gribapiDefDateTimeAbs(int editionNumber, grib_handle *gh, int64_t date, int time, int productDefinitionTemplate, int typeOfGeneratingProcess, int tsteptype, int gcinit)
1940 {
1941 (void ) gribapiDefSteptype(editionNumber, gh, productDefinitionTemplate, typeOfGeneratingProcess, tsteptype, gcinit);
1942
1943 if ( editionNumber > 1 ) GRIB_CHECK(my_grib_set_long(gh, "significanceOfReferenceTime", 0), 0);
1944 if ( editionNumber > 1 ) GRIB_CHECK(my_grib_set_long(gh, "stepRange", 0), 0);
1945
1946 if ( date == 0 ) date = 10101;
1947 gribapiSetDataDateTime(gh, date, time);
1948 }
1949
1950 static
gribapiDefDateTimeRel(int editionNumber,grib_handle * gh,int64_t fdate,int ftime,int64_t vdate,int vtime,int64_t sdate,int stime,int productDefinitionTemplate,int typeOfGeneratingProcess,int tsteptype,int timeunit,int calendar,int gcinit)1951 int gribapiDefDateTimeRel(int editionNumber, grib_handle *gh, int64_t fdate, int ftime, int64_t vdate, int vtime, int64_t sdate, int stime,
1952 int productDefinitionTemplate, int typeOfGeneratingProcess, int tsteptype, int timeunit, int calendar, int gcinit)
1953 {
1954 int status = -1;
1955 int year, month, day, hour, minute, second;
1956 int64_t julday1, julday2, days;
1957 int secofday1, secofday2, secs;
1958
1959 cdiDecodeDate(fdate, &year, &month, &day);
1960 cdiDecodeTime(ftime, &hour, &minute, &second);
1961 encode_juldaysec(calendar, year, month, day, hour, minute, second, &julday1, &secofday1);
1962
1963 if ( vdate == 0 && vtime == 0 ) { vdate = fdate; vtime = ftime; }
1964
1965 cdiDecodeDate(vdate, &year, &month, &day);
1966 cdiDecodeTime(vtime, &hour, &minute, &second);
1967 encode_juldaysec(calendar, year, month, day, hour, minute, second, &julday2, &secofday2);
1968
1969 (void) julday_sub(julday1, secofday1, julday2, secofday2, &days, &secs);
1970
1971 const int factor = getTimeunitFactor(timeunit);
1972
1973 if ( !(int)(fmod(days*86400.0 + secs, factor)))
1974 {
1975 const int proDefTempNum = gribapiDefSteptype(editionNumber, gh, productDefinitionTemplate, typeOfGeneratingProcess, tsteptype, gcinit);
1976 gribapiDefStepUnits(editionNumber, gh, timeunit, proDefTempNum, gcinit);
1977
1978 long startStep = 0;
1979 const double endStepF = (days*86400.0 + secs)/factor;
1980 const long maxStep = (editionNumber > 1) ? INT_MAX : 65000;
1981 if (endStepF > maxStep) return status;
1982 long endStep = (long) endStepF;
1983
1984 if (sdate != 0 && (tsteptype == TSTEP_RANGE || tsteptype == TSTEP_AVG || tsteptype == TSTEP_ACCUM || tsteptype == TSTEP_DIFF))
1985 {
1986 cdiDecodeDate(sdate, &year, &month, &day);
1987 cdiDecodeTime(stime, &hour, &minute, &second);
1988 encode_juldaysec(calendar, year, month, day, hour, minute, second, &julday2, &secofday2);
1989
1990 (void) julday_sub(julday1, secofday1, julday2, secofday2, &days, &secs);
1991
1992 startStep = (long) ((days*86400.0 + secs)/factor);
1993 }
1994
1995 if ( editionNumber > 1 ) GRIB_CHECK(my_grib_set_long(gh, "significanceOfReferenceTime", 1), 0);
1996 if ( editionNumber > 1 ) GRIB_CHECK(my_grib_set_long(gh, "stepRange", 0), 0);
1997
1998 if ( fdate == 0 ) fdate = 10101;
1999 gribapiSetDataDateTime(gh, fdate, ftime);
2000
2001 // printf(">>>>> tsteptype %d startStep %ld endStep %ld\n", tsteptype, startStep, endStep);
2002
2003 // Product Definition Template Number: defined in GRIB_API file 4.0.table point in time products:
2004 if ( (proDefTempNum >= 0 && proDefTempNum <= 7) ||
2005 proDefTempNum == 55 || proDefTempNum == 40055 ) // Tile
2006 startStep = endStep;
2007
2008 if (endStep < startStep) return status;
2009
2010 if ( editionNumber == 1 && (startStep > 255 || endStep > 255) )
2011 {
2012 startStep = 0;
2013 endStep = 0;
2014 }
2015 else
2016 {
2017 status = 0;
2018 }
2019
2020 if ( editionNumber > 1 ) GRIB_CHECK(my_grib_set_long(gh, "forecastTime", startStep), 0);
2021 //if ( editionNumber == 1 && startStep > 0) GRIB_CHECK(my_grib_set_long(gh, "startStep", startStep), 0);
2022 if ( editionNumber == 1 ) GRIB_CHECK(my_grib_set_long(gh, "startStep", startStep), 0);
2023 GRIB_CHECK(my_grib_set_long(gh, "endStep", endStep), 0);
2024 }
2025
2026 return status;
2027 }
2028
2029 static
gribapiDefTime(int editionNumber,int productDefinitionTemplate,int typeOfGeneratingProcess,grib_handle * gh,int64_t vdate,int vtime,int tsteptype,int numavg,int taxisID,int gcinit)2030 void gribapiDefTime(int editionNumber, int productDefinitionTemplate, int typeOfGeneratingProcess, grib_handle *gh,
2031 int64_t vdate, int vtime, int tsteptype, int numavg, int taxisID, int gcinit)
2032 {
2033 UNUSED(numavg);
2034
2035 int taxistype = (taxisID == -1) ? -1 : taxisInqType(taxisID);
2036
2037 if (typeOfGeneratingProcess == 196)
2038 {
2039 vdate = 10101;
2040 vtime = 0;
2041 taxistype = TAXIS_ABSOLUTE;
2042 }
2043
2044 if (taxistype == TAXIS_RELATIVE)
2045 {
2046 const int timeunit = taxisInqTunit(taxisID);
2047 const int calendar = taxisInqCalendar(taxisID);
2048
2049 int64_t fdate = taxisInqFdate(taxisID);
2050 int ftime = taxisInqFtime(taxisID);
2051 if (fdate == CDI_UNDEFID)
2052 {
2053 fdate = taxisInqRdate(taxisID);
2054 ftime = taxisInqRtime(taxisID);
2055 }
2056 if (vdate < fdate || (vdate == fdate && vtime < ftime))
2057 {
2058 fdate = vdate;
2059 ftime = vtime;
2060 }
2061
2062 const int64_t sdate = taxisInqSdate(taxisID);
2063 const int stime = taxisInqStime(taxisID);
2064
2065 int status = gribapiDefDateTimeRel(editionNumber, gh, fdate, ftime, vdate, vtime, sdate, stime, productDefinitionTemplate,
2066 typeOfGeneratingProcess, tsteptype, timeunit, calendar, gcinit);
2067 if (status != 0) taxistype = TAXIS_ABSOLUTE;
2068 }
2069
2070 if (taxistype == TAXIS_ABSOLUTE)
2071 gribapiDefDateTimeAbs(editionNumber, gh, vdate, vtime, productDefinitionTemplate, typeOfGeneratingProcess, tsteptype, gcinit);
2072 }
2073
2074 static
gribapiDefGridRegular(grib_handle * gh,int gridID,int gridtype,bool gridIsRotated,bool gridIsCurvilinear,int uvRelativeToGrid)2075 void gribapiDefGridRegular(grib_handle *gh, int gridID, int gridtype, bool gridIsRotated, bool gridIsCurvilinear, int uvRelativeToGrid)
2076 {
2077 const char *mesg;
2078 // clang-format off
2079 if (gridtype == GRID_GAUSSIAN) mesg = "regular_gg";
2080 else if (gridtype == GRID_GAUSSIAN_REDUCED) mesg = "reduced_gg";
2081 else if (gridIsRotated) mesg = "rotated_ll";
2082 else mesg = "regular_ll";
2083 // clang-format on
2084 size_t len = strlen(mesg);
2085 GRIB_CHECK(my_grib_set_string(gh, "gridType", mesg, &len), 0);
2086
2087 double xfirst = 0.0, xlast = 0.0, xinc = 0.0;
2088 double yfirst = 0.0, ylast = 0.0, yinc = 0.0;
2089
2090 size_t nlon = gridInqXsize(gridID);
2091 size_t nlat = gridInqYsize(gridID);
2092
2093 if ( gridtype == GRID_GAUSSIAN_REDUCED )
2094 {
2095 xfirst = (nlon == 2) ? gridInqXval(gridID, 0) : 0.0;
2096 xlast = (nlon == 2) ? gridInqXval(gridID, 1) : 360.0 - 360.0 * 0.5 / (double)nlat;
2097
2098 nlon = 0;
2099
2100 int *reducedPoints = (int *) Malloc(nlat*sizeof(int));
2101 long *pl = (long *) Malloc(nlat*sizeof(long));
2102 gridInqReducedPoints(gridID, reducedPoints);
2103 for ( size_t i = 0; i < nlat; ++i ) pl[i] = reducedPoints[i];
2104
2105 GRIB_CHECK(grib_set_long_array(gh, "pl", pl, nlat), 0);
2106
2107 Free(pl);
2108 Free(reducedPoints);
2109 }
2110 else
2111 {
2112 if ( nlon == 0 ) nlon = 1;
2113 else
2114 {
2115 xfirst = gridInqXval(gridID, 0);
2116 xlast = gridInqXval(gridID, (gridIsCurvilinear ? nlon*nlat : nlon) - 1);
2117 xinc = fabs(gridInqXinc(gridID));
2118 }
2119 }
2120
2121 if ( nlat == 0 ) nlat = 1;
2122 else
2123 {
2124 yfirst = gridInqYval(gridID, 0);
2125 ylast = gridInqYval(gridID, (gridIsCurvilinear ? nlon*nlat : nlat) - 1);
2126 yinc = fabs(gridInqYinc(gridID));
2127 }
2128
2129 double xfirsto = xfirst;
2130 double xlasto = xlast;
2131 while ( xfirsto > 360.0 ) xfirsto -= 360.0;
2132 while ( xlasto > 360.0 ) xlasto -= 360.0;
2133
2134 if ( gridtype != GRID_GAUSSIAN_REDUCED ) GRIB_CHECK(my_grib_set_long(gh, "Ni", nlon), 0);
2135 GRIB_CHECK(my_grib_set_double(gh, "longitudeOfFirstGridPointInDegrees", xfirsto), 0);
2136 GRIB_CHECK(my_grib_set_double(gh, "longitudeOfLastGridPointInDegrees", xlasto), 0);
2137 if ( gridtype != GRID_GAUSSIAN_REDUCED ) GRIB_CHECK(my_grib_set_double(gh, "iDirectionIncrementInDegrees", xinc), 0);
2138
2139 GRIB_CHECK(my_grib_set_long(gh, "Nj", (long)nlat), 0);
2140 GRIB_CHECK(my_grib_set_double(gh, "latitudeOfFirstGridPointInDegrees", yfirst), 0);
2141 GRIB_CHECK(my_grib_set_double(gh, "latitudeOfLastGridPointInDegrees", ylast), 0);
2142
2143 long xscan = xfirst > xlast;
2144 GRIB_CHECK(my_grib_set_long(gh, "iScansNegatively", xscan), 0);
2145 xscan = yfirst < ylast;
2146 GRIB_CHECK(my_grib_set_long(gh, "jScansPositively", xscan), 0);
2147
2148 if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_GAUSSIAN_REDUCED )
2149 {
2150 int np = gridInqNP(gridID);
2151 if ( np == 0 ) np = nlat/2;
2152 GRIB_CHECK(my_grib_set_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", np), 0);
2153 }
2154 else
2155 {
2156 GRIB_CHECK(my_grib_set_double(gh, "jDirectionIncrementInDegrees", yinc), 0);
2157 }
2158
2159 if ( gridIsRotated )
2160 {
2161 double xpole = 0.0, ypole = 0.0, angle = 0.0;
2162 gridInqParamRLL(gridID, &xpole, &ypole, &angle);
2163
2164 xpole += 180.0;
2165 if ( fabs(ypole) > 0.0 ) ypole = -ypole; // change from north to south pole
2166 if ( fabs(angle) > 0.0 ) angle = -angle;
2167 GRIB_CHECK(my_grib_set_double(gh, "latitudeOfSouthernPoleInDegrees", ypole), 0);
2168 GRIB_CHECK(my_grib_set_double(gh, "longitudeOfSouthernPoleInDegrees", xpole), 0);
2169 GRIB_CHECK(my_grib_set_double(gh, "angleOfRotation", angle), 0);
2170 }
2171
2172 if ( uvRelativeToGrid >= 0 ) GRIB_CHECK(my_grib_set_long(gh, "uvRelativeToGrid", uvRelativeToGrid), 0);
2173 }
2174
2175 static
radiusToShapeOfTheEarth(double radius)2176 long radiusToShapeOfTheEarth(double radius)
2177 {
2178 long shapeOfTheEarth = 0;
2179
2180 if (IS_EQUAL(radius, 6367470.0)) shapeOfTheEarth = 0;
2181 else if (IS_EQUAL(radius, 6378160.0)) shapeOfTheEarth = 2;
2182 else if (IS_EQUAL(radius, 6371229.0)) shapeOfTheEarth = 6;
2183 else if (IS_EQUAL(radius, 6371200.0)) shapeOfTheEarth = 8;
2184
2185 return shapeOfTheEarth;
2186 }
2187
2188 static
gribapiDefGridLCC(grib_handle * gh,int editionNumber,int gridID,int uvRelativeToGrid)2189 void gribapiDefGridLCC(grib_handle *gh, int editionNumber, int gridID, int uvRelativeToGrid)
2190 {
2191 const long xsize = (long) gridInqXsize(gridID);
2192 const long ysize = (long) gridInqYsize(gridID);
2193
2194 struct CDI_GridProjParams gpp;
2195 gridInqParamsLCC(gridID, &gpp);
2196 if (IS_EQUAL(gpp.x_0, gpp.mv) && IS_EQUAL(gpp.y_0, gpp.mv) && (IS_EQUAL(gpp.xval_0, gpp.mv) || IS_EQUAL(gpp.yval_0, gpp.mv)))
2197 {
2198 gpp.x_0 = gridInqXval(gridID, 0);
2199 gpp.y_0 = gridInqYval(gridID, 0);
2200 }
2201 gridVerifyProjParamsLCC(&gpp);
2202 if (gpp.xval_0 < 0.0) gpp.xval_0 += 360.0;
2203 if (gpp.lon_0 < 0.0) gpp.lon_0 += 360.0;
2204
2205 const bool lsouth = (gpp.lat_1 < 0.0);
2206 if (lsouth) { gpp.lat_1 = -gpp.lat_2; gpp.lat_2 = -gpp.lat_2; }
2207 int projflag = 0;
2208 if (lsouth) gribbyte_set_bit(&projflag, 1);
2209
2210 double xinc = gridInqXinc(gridID);
2211 double yinc = gridInqYinc(gridID);
2212 if (IS_EQUAL(xinc, 0.0)) xinc = gridInqXincInMeter(gridID);
2213 if (IS_EQUAL(yinc, 0.0)) yinc = gridInqYincInMeter(gridID);
2214
2215 static const char mesg[] = "lambert";
2216 size_t len = sizeof(mesg) -1;
2217 GRIB_CHECK(my_grib_set_string(gh, "gridType", mesg, &len), 0);
2218
2219 GRIB_CHECK(my_grib_set_long(gh, "Nx", xsize), 0);
2220 GRIB_CHECK(my_grib_set_long(gh, "Ny", ysize), 0);
2221 GRIB_CHECK(my_grib_set_long(gh, "DxInMetres", lround(fabs(xinc))), 0);
2222 GRIB_CHECK(my_grib_set_long(gh, "DyInMetres", lround(fabs(yinc))), 0);
2223 GRIB_CHECK(my_grib_set_double(gh, "longitudeOfFirstGridPointInDegrees", gpp.xval_0), 0);
2224 GRIB_CHECK(my_grib_set_double(gh, "latitudeOfFirstGridPointInDegrees", gpp.yval_0), 0);
2225 GRIB_CHECK(my_grib_set_double(gh, "LoVInDegrees", gpp.lon_0), 0);
2226 GRIB_CHECK(my_grib_set_double(gh, "Latin1InDegrees", gpp.lat_1), 0);
2227 GRIB_CHECK(my_grib_set_double(gh, "Latin2InDegrees", gpp.lat_2), 0);
2228 GRIB_CHECK(my_grib_set_long(gh, "projectionCentreFlag", projflag), 0);
2229
2230
2231 const long shapeOfTheEarth = radiusToShapeOfTheEarth(gpp.a);
2232 if (shapeOfTheEarth) GRIB_CHECK(my_grib_set_long(gh, "shapeOfTheEarth", shapeOfTheEarth), 0);
2233
2234 if (uvRelativeToGrid >= 0) GRIB_CHECK(my_grib_set_long(gh, "uvRelativeToGrid", uvRelativeToGrid), 0);
2235
2236 const long earthIsOblate = (IS_EQUAL(gpp.a, 6378160.0) && IS_EQUAL(gpp.rf, 297.0));
2237 if (earthIsOblate) GRIB_CHECK(my_grib_set_long(gh, "earthIsOblate", earthIsOblate), 0);
2238
2239 int scanflag = 0;
2240 gribbyte_set_bit(&scanflag, 2);
2241 if ( editionNumber <= 1 )
2242 GRIB_CHECK(my_grib_set_long(gh, "scanningMode", (long)scanflag), 0);
2243 }
2244
2245 static
gribapiDefGridSTERE(grib_handle * gh,int gridID)2246 void gribapiDefGridSTERE(grib_handle *gh, int gridID)
2247 {
2248 const long xsize = (long) gridInqXsize(gridID);
2249 const long ysize = (long) gridInqYsize(gridID);
2250
2251 struct CDI_GridProjParams gpp;
2252 gridInqParamsSTERE(gridID, &gpp);
2253 gridVerifyProjParamsSTERE(&gpp);
2254 if (gpp.xval_0 < 0.0) gpp.xval_0 += 360.0;
2255 int projflag = 0;
2256
2257 const double xinc = gridInqXinc(gridID);
2258 const double yinc = gridInqYinc(gridID);
2259
2260 static const char mesg[] = "polar_stereographic";
2261 size_t len = sizeof(mesg) -1;
2262 GRIB_CHECK(my_grib_set_string(gh, "gridType", mesg, &len), 0);
2263
2264 GRIB_CHECK(my_grib_set_long(gh, "Nx", xsize), 0);
2265 GRIB_CHECK(my_grib_set_long(gh, "Ny", ysize), 0);
2266 GRIB_CHECK(my_grib_set_long(gh, "DxInMetres", lround(xinc)), 0);
2267 GRIB_CHECK(my_grib_set_long(gh, "DyInMetres", lround(yinc)), 0);
2268 GRIB_CHECK(my_grib_set_double(gh, "longitudeOfFirstGridPointInDegrees", gpp.xval_0), 0);
2269 GRIB_CHECK(my_grib_set_double(gh, "latitudeOfFirstGridPointInDegrees", gpp.yval_0), 0);
2270 GRIB_CHECK(my_grib_set_double(gh, "LaDInDegrees", gpp.lat_1), 0);
2271 GRIB_CHECK(my_grib_set_double(gh, "orientationOfTheGridInDegrees", gpp.lon_0), 0);
2272 const long southPoleOnProjectionPlane = IS_EQUAL(gpp.lat_0, -90.0);
2273 GRIB_CHECK(my_grib_set_double(gh, "southPoleOnProjectionPlane", southPoleOnProjectionPlane), 0);
2274 GRIB_CHECK(my_grib_set_long(gh, "projectionCentreFlag", projflag), 0);
2275
2276 const long shapeOfTheEarth = radiusToShapeOfTheEarth(gpp.a);
2277 if (shapeOfTheEarth) GRIB_CHECK(my_grib_set_long(gh, "shapeOfTheEarth", shapeOfTheEarth), 0);
2278
2279 GRIB_CHECK(my_grib_set_long(gh, "resolutionAndComponentFlags", 8), 0);
2280 GRIB_CHECK(my_grib_set_long(gh, "iScansNegatively", 0), 0);
2281 GRIB_CHECK(my_grib_set_long(gh, "jScansPositively", 1), 0);
2282 }
2283
2284 static
gribapiDefGridGME(grib_handle * gh,int gridID,long gridsize)2285 void gribapiDefGridGME(grib_handle *gh, int gridID, long gridsize)
2286 {
2287 GRIB_CHECK(my_grib_set_long(gh, "gridDefinitionTemplateNumber", GRIB2_GTYPE_GME), 0);
2288
2289 int nd = 0, ni = 0, ni2 = 0, ni3 = 0;
2290 gridInqParamGME(gridID, &nd, &ni, &ni2, &ni3);
2291 GRIB_CHECK(my_grib_set_long(gh, "nd", nd), 0);
2292 GRIB_CHECK(my_grib_set_long(gh, "Ni", ni), 0);
2293 GRIB_CHECK(my_grib_set_long(gh, "n2", ni2), 0);
2294 GRIB_CHECK(my_grib_set_long(gh, "n3", ni3), 0);
2295 GRIB_CHECK(my_grib_set_long(gh, "latitudeOfThePolePoint", 90000000), 0);
2296 GRIB_CHECK(my_grib_set_long(gh, "longitudeOfThePolePoint", 0), 0);
2297
2298 GRIB_CHECK(my_grib_set_long(gh, "numberOfDataPoints", gridsize), 0);
2299 GRIB_CHECK(my_grib_set_long(gh, "totalNumberOfGridPoints", gridsize), 0);
2300 }
2301
2302 static
gribapiDefGridUnstructured(grib_handle * gh,int gridID)2303 void gribapiDefGridUnstructured(grib_handle *gh, int gridID)
2304 {
2305 static bool warning = true;
2306
2307 int status = my_grib_set_long(gh, "gridDefinitionTemplateNumber", GRIB2_GTYPE_UNSTRUCTURED);
2308 if (status != 0 && warning)
2309 {
2310 warning = false;
2311 Warning("Can't write reference grid!");
2312 Warning("gridDefinitionTemplateNumber %d not found (grib2/template.3.%d.def)!",
2313 GRIB2_GTYPE_UNSTRUCTURED, GRIB2_GTYPE_UNSTRUCTURED);
2314 }
2315 else
2316 {
2317 int errCount = 0;
2318 int number = 0;
2319 status = cdiInqKeyInt(gridID, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDUSED, &number);
2320 if (status < 0) errCount++;
2321 if (number < 0) number = 0;
2322 GRIB_CHECK(my_grib_set_long(gh, "numberOfGridUsed", number), 0);
2323
2324 int position = 0;
2325 status = cdiInqKeyInt(gridID, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDINREFERENCE, &position);
2326 if (status < 0) errCount++;
2327 if (position < 0) position = 0;
2328 GRIB_CHECK(my_grib_set_long(gh, "numberOfGridInReference", position), 0);
2329
2330 unsigned char uuid[CDI_UUID_SIZE];
2331 size_t len = CDI_UUID_SIZE;
2332 memset(uuid, 0, len);
2333 int length = CDI_UUID_SIZE;
2334 status = cdiInqKeyBytes(gridID, CDI_GLOBAL, CDI_KEY_UUID, uuid, &length);
2335 if (status < 0) errCount++;
2336 if (grib_set_bytes(gh, "uuidOfHGrid", uuid, &len) != 0) Warning("Can't write UUID!");
2337
2338 if (warning && errCount > 0)
2339 {
2340 warning = false;
2341 Warning("GRIB library doesn't support unstructured grids, grid information will be lost!");
2342 }
2343 }
2344 }
2345
2346 static
gribapiDefGridSpectral(grib_handle * gh,int gridID)2347 void gribapiDefGridSpectral(grib_handle *gh, int gridID)
2348 {
2349 const int trunc = gridInqTrunc(gridID);
2350 enum { numTruncAtt = 3 };
2351 static const char truncAttNames[numTruncAtt][2] = { "J", "K", "M" };
2352 for (size_t i = 0; i < numTruncAtt; ++i)
2353 GRIB_CHECK(my_grib_set_long(gh, truncAttNames[i], trunc), 0);
2354
2355 if ( gridInqComplexPacking(gridID) )
2356 {
2357 static const char truncAttNames2[numTruncAtt][3] = { "JS", "KS", "MS" };
2358 for (size_t i = 0; i < numTruncAtt; ++i)
2359 GRIB_CHECK(my_grib_set_long(gh, truncAttNames2[i], 20), 0);
2360 }
2361 }
2362
2363 static
gribapiDefPackingType(grib_handle * gh,bool lieee,bool lspectral,bool lcomplex,int comptype,size_t gridsize)2364 void gribapiDefPackingType(grib_handle *gh, bool lieee, bool lspectral, bool lcomplex, int comptype, size_t gridsize)
2365 {
2366 static const char mesg_spectral_complex[] = "spectral_complex";
2367 static const char mesg_spectral_simple[] = "spectral_simple";
2368 static const char mesg_grid_jpeg[] = "grid_jpeg";
2369 static const char mesg_grid_ccsds[] = "grid_ccsds";
2370 static const char mesg_ieee[] = "grid_ieee";
2371 static const char mesg_simple[] = "grid_simple";
2372 const char *mesg = mesg_simple;
2373
2374 if ( lspectral )
2375 {
2376 mesg = lcomplex ? mesg_spectral_complex : mesg_spectral_simple;
2377 }
2378 else if ( comptype == CDI_COMPRESS_JPEG && gridsize > 1 )
2379 {
2380 mesg = mesg_grid_jpeg;
2381 }
2382 else if ( (comptype == CDI_COMPRESS_SZIP || comptype == CDI_COMPRESS_AEC) && gridsize > 1 )
2383 {
2384 mesg = mesg_grid_ccsds;
2385 }
2386 else if ( lieee )
2387 {
2388 mesg = mesg_ieee;
2389 }
2390
2391 size_t len = strlen(mesg);
2392 GRIB_CHECK(my_grib_set_string(gh, "packingType", mesg, &len), 0);
2393 }
2394
2395 static
gribapiDefGrid(int editionNumber,grib_handle * gh,int gridID,int comptype,int datatype,int uvRelativeToGrid)2396 void gribapiDefGrid(int editionNumber, grib_handle *gh, int gridID, int comptype, int datatype, int uvRelativeToGrid)
2397 {
2398 const size_t gridsize = gridInqSize(gridID);
2399 bool gridIsRotated = false;
2400 bool gridIsCurvilinear = false;
2401 const int gridtype = grbGetGridtype(&gridID, gridsize, &gridIsRotated, &gridIsCurvilinear);
2402
2403 bool lieee = (editionNumber == 2 && (datatype == CDI_DATATYPE_FLT32 || datatype == CDI_DATATYPE_FLT64));
2404 bool lspectral = (gridtype == GRID_SPECTRAL);
2405 bool lcomplex = (lspectral && gridInqComplexPacking(gridID));
2406
2407 if ( lieee ) comptype = 0;
2408 if ( lspectral ) lieee = false;
2409
2410 if ( lspectral ) // gridType needs to be defined before packingType !!!
2411 {
2412 static const char mesg[] = "sh";
2413 size_t len = sizeof (mesg) -1;
2414 GRIB_CHECK(my_grib_set_string(gh, "gridType", mesg, &len), 0);
2415 }
2416
2417 gribapiDefPackingType(gh, lieee, lspectral, lcomplex, comptype, gridsize);
2418
2419 if ( lieee ) GRIB_CHECK(my_grib_set_long(gh, "precision", datatype == CDI_DATATYPE_FLT64 ? 2 : 1), 0);
2420
2421 if ( editionNumber == 2 ) GRIB_CHECK(my_grib_set_long(gh, "numberOfValues", (long)gridsize), 0);
2422
2423 switch (gridtype)
2424 {
2425 case GRID_LONLAT:
2426 case GRID_GAUSSIAN:
2427 case GRID_GAUSSIAN_REDUCED:
2428 case GRID_TRAJECTORY:
2429 {
2430 gribapiDefGridRegular(gh, gridID, gridtype, gridIsRotated, gridIsCurvilinear, uvRelativeToGrid);
2431 break;
2432 }
2433 case CDI_PROJ_LCC:
2434 {
2435 gribapiDefGridLCC(gh, editionNumber, gridID, uvRelativeToGrid);
2436 break;
2437 }
2438 case CDI_PROJ_STERE:
2439 {
2440 gribapiDefGridSTERE(gh, gridID);
2441 break;
2442 }
2443 case GRID_SPECTRAL:
2444 {
2445 gribapiDefGridSpectral(gh, gridID);
2446 break;
2447 }
2448 case GRID_GME:
2449 {
2450 if ( editionNumber <= 1 ) Error("GME grid can't be stored in GRIB edition %d!", editionNumber);
2451 gribapiDefGridGME(gh, gridID, (long) gridsize);
2452 break;
2453 }
2454 case GRID_UNSTRUCTURED:
2455 {
2456 if ( editionNumber <= 1 ) Error("Unstructured grid can't be stored in GRIB edition %d!", editionNumber);
2457 gribapiDefGridUnstructured(gh, gridID);
2458 break;
2459 }
2460 default:
2461 {
2462 Error("Unsupported grid type: %s", gridNamePtr(gridtype));
2463 break;
2464 }
2465 }
2466 }
2467
2468 static
getLevelFactor(double level,long * factor,long * out_scaled_value)2469 void getLevelFactor(double level, long *factor, long *out_scaled_value)
2470 {
2471 double scaled_value = level;
2472 long iscaled_value = lround(scaled_value);
2473 long i;
2474
2475 const double eps = 1.e-8;
2476 for ( i=0; (fabs(scaled_value - (double) iscaled_value) >= eps) && i < 7; i++ )
2477 {
2478 scaled_value *= 10.;
2479 iscaled_value = lround(scaled_value);
2480 }
2481
2482 (*factor) = i;
2483 (*out_scaled_value) = iscaled_value;
2484 }
2485
2486 static
gribapiDefLevelType(grib_handle * gh,int gcinit,const char * keyname,long leveltype)2487 void gribapiDefLevelType(grib_handle *gh, int gcinit, const char *keyname, long leveltype)
2488 {
2489 bool lset = false;
2490 if ( (leveltype == GRIB1_LTYPE_ISOBARIC_PA || leveltype == 99 || leveltype == 100) && gribEditionNumber(gh) == 1 )
2491 {
2492 if ( gribGetLong(gh, "indicatorOfTypeOfLevel") != leveltype ) lset = true;
2493 }
2494
2495 if ( !gcinit || lset ) GRIB_CHECK(my_grib_set_long(gh, keyname, leveltype), 0);
2496 }
2497
2498 static
grib1DefLevel(grib_handle * gh,int gcinit,long leveltype,bool hasBounds,double level,double dlevel1,double dlevel2)2499 void grib1DefLevel(grib_handle *gh, int gcinit, long leveltype, bool hasBounds, double level, double dlevel1, double dlevel2)
2500 {
2501 gribapiDefLevelType(gh, gcinit, "indicatorOfTypeOfLevel", leveltype);
2502
2503 if (hasBounds)
2504 {
2505 GRIB_CHECK(my_grib_set_long(gh, "topLevel", lround(dlevel1)), 0);
2506 GRIB_CHECK(my_grib_set_long(gh, "bottomLevel", lround(dlevel2)), 0);
2507 }
2508 else
2509 {
2510 GRIB_CHECK(my_grib_set_long(gh, "level", lround(level)), 0);
2511 }
2512 }
2513
2514 static
grib2DefLevel(grib_handle * gh,int gcinit,long leveltype1,long leveltype2,bool hasBounds,double level,double dlevel1,double dlevel2)2515 void grib2DefLevel(grib_handle *gh, int gcinit, long leveltype1, long leveltype2, bool hasBounds, double level, double dlevel1, double dlevel2)
2516 {
2517 gribapiDefLevelType(gh, gcinit, "typeOfFirstFixedSurface", leveltype1);
2518 if (hasBounds) gribapiDefLevelType(gh, gcinit, "typeOfSecondFixedSurface", leveltype2);
2519
2520 if (!hasBounds) dlevel1 = level;
2521
2522 long scaled_level, factor;
2523 getLevelFactor(dlevel1, &factor, &scaled_level);
2524 GRIB_CHECK(my_grib_set_long(gh, "scaleFactorOfFirstFixedSurface", factor), 0);
2525 GRIB_CHECK(my_grib_set_long(gh, "scaledValueOfFirstFixedSurface", scaled_level), 0);
2526
2527 if (hasBounds)
2528 {
2529 getLevelFactor(dlevel2, &factor, &scaled_level);
2530 GRIB_CHECK(my_grib_set_long(gh, "scaleFactorOfSecondFixedSurface", factor), 0);
2531 GRIB_CHECK(my_grib_set_long(gh, "scaledValueOfSecondFixedSurface", scaled_level), 0);
2532 }
2533 }
2534
2535 static
gribapiDefLevel(int editionNumber,grib_handle * gh,int zaxisID,int levelID,int gcinit,int proddef_template_num)2536 void gribapiDefLevel(int editionNumber, grib_handle *gh, int zaxisID, int levelID, int gcinit, int proddef_template_num)
2537 {
2538 int zaxistype = zaxisInqType(zaxisID);
2539 int ltype = 0, ltype2 = -1;
2540 cdiInqKeyInt(zaxisID, CDI_GLOBAL, CDI_KEY_TYPEOFFIRSTFIXEDSURFACE, <ype);
2541 cdiInqKeyInt(zaxisID, CDI_GLOBAL, CDI_KEY_TYPEOFSECONDFIXEDSURFACE, <ype2);
2542
2543 bool hasBounds = (zaxisInqLbounds(zaxisID, NULL) && zaxisInqUbounds(zaxisID, NULL));
2544 double level = zaxisInqLevels(zaxisID, NULL) ? zaxisInqLevel(zaxisID, levelID) : levelID + 1;
2545 double dlevel1 = hasBounds ? zaxisInqLbound(zaxisID, levelID) : level;
2546 double dlevel2 = hasBounds ? zaxisInqUbound(zaxisID, levelID) : 0.0;
2547
2548 if ( zaxistype == ZAXIS_GENERIC && ltype == 0 )
2549 {
2550 Message("Changed zaxis type from %s to %s", zaxisNamePtr(zaxistype), zaxisNamePtr(ZAXIS_PRESSURE));
2551 zaxistype = ZAXIS_PRESSURE;
2552 zaxisChangeType(zaxisID, zaxistype);
2553 cdiDefKeyString(zaxisID, CDI_GLOBAL, CDI_KEY_UNITS, "Pa");
2554 }
2555
2556 long grib_ltype = (editionNumber <= 1) ? zaxisTypeToGrib1ltype(zaxistype) : zaxisTypeToGrib2ltype(zaxistype);
2557 long grib_ltype2 = (ltype != ltype2 && ltype2 != -1) ? ltype2 : grib_ltype;
2558
2559 switch (zaxistype)
2560 {
2561 case ZAXIS_SURFACE:
2562 case ZAXIS_MEANSEA:
2563 case ZAXIS_HEIGHT:
2564 case ZAXIS_ALTITUDE:
2565 case ZAXIS_SIGMA:
2566 case ZAXIS_DEPTH_BELOW_SEA:
2567 case ZAXIS_ISENTROPIC:
2568 {
2569 if (zaxistype == ZAXIS_HEIGHT)
2570 {
2571 const double sf = zaxis_units_to_meter(zaxisID);
2572 level *= sf;
2573 dlevel1 *= sf;
2574 dlevel2 *= sf;
2575 }
2576
2577 if (editionNumber <= 1)
2578 {
2579 grib1DefLevel(gh, gcinit, grib_ltype, hasBounds, level, dlevel1, dlevel2);
2580 }
2581 else
2582 {
2583 /* PRODUCT DEFINITION TEMPLATE NUMBER 32:
2584
2585 "Analysis or forecast at a horizontal level or in a horizontal layer at a point
2586 in time for simulate (synthetic) satellite data"
2587
2588 The key/value pairs that are set in "grib2DefLevel" do not exist for this template.
2589 */
2590 if ( proddef_template_num != 32 )
2591 grib2DefLevel(gh, gcinit, grib_ltype, grib_ltype2, hasBounds, level, dlevel1, dlevel2);
2592 }
2593
2594 break;
2595 }
2596 case ZAXIS_CLOUD_BASE:
2597 case ZAXIS_CLOUD_TOP:
2598 case ZAXIS_ISOTHERM_ZERO:
2599 case ZAXIS_TROPOPAUSE:
2600 case ZAXIS_TOA:
2601 case ZAXIS_SEA_BOTTOM:
2602 case ZAXIS_LAKE_BOTTOM:
2603 case ZAXIS_SEDIMENT_BOTTOM:
2604 case ZAXIS_SEDIMENT_BOTTOM_TA:
2605 case ZAXIS_SEDIMENT_BOTTOM_TW:
2606 case ZAXIS_MIX_LAYER:
2607 case ZAXIS_ATMOSPHERE:
2608 {
2609 if (editionNumber <= 1)
2610 grib1DefLevel(gh, gcinit, grib_ltype, hasBounds, level, dlevel1, dlevel2);
2611 else
2612 grib2DefLevel(gh, gcinit, grib_ltype, grib_ltype2, hasBounds, level, dlevel1, dlevel2);
2613
2614 break;
2615 }
2616 case ZAXIS_HYBRID:
2617 case ZAXIS_HYBRID_HALF:
2618 {
2619 if (editionNumber <= 1)
2620 {
2621 grib_ltype = hasBounds ? GRIB1_LTYPE_HYBRID_LAYER : GRIB1_LTYPE_HYBRID;
2622 grib1DefLevel(gh, gcinit, grib_ltype, hasBounds, level, dlevel1, dlevel2);
2623 }
2624 else
2625 {
2626 grib2DefLevel(gh, gcinit, grib_ltype, grib_ltype, hasBounds, level, dlevel1, dlevel2);
2627 }
2628
2629 if (!gcinit)
2630 {
2631 int vctsize = zaxisInqVctSize(zaxisID);
2632 if (vctsize > 0)
2633 {
2634 GRIB_CHECK(my_grib_set_long(gh, "PVPresent", 1), 0);
2635 GRIB_CHECK(grib_set_double_array(gh, "pv", zaxisInqVctPtr(zaxisID), (size_t)vctsize), 0);
2636 }
2637 }
2638
2639 break;
2640 }
2641 case ZAXIS_PRESSURE:
2642 {
2643 if (level < 0) Warning("Pressure level of %f Pa is below zero!", level);
2644
2645 if (!zaxis_units_is_Pa(zaxisID))
2646 {
2647 level *= 100;
2648 dlevel1 *= 100;
2649 dlevel2 *= 100;
2650 }
2651
2652 if (editionNumber <= 1)
2653 {
2654 double dum;
2655 if (level < 32768 && (level < 100 || modf(level/100, &dum) > 0))
2656 grib_ltype = GRIB1_LTYPE_ISOBARIC_PA;
2657 else
2658 level /= 100;
2659
2660 grib1DefLevel(gh, gcinit, grib_ltype, hasBounds, level, dlevel1, dlevel2);
2661 }
2662 else
2663 {
2664 if (ltype2 == -1) ltype2 = GRIB2_LTYPE_ISOBARIC;
2665 grib2DefLevel(gh, gcinit, GRIB2_LTYPE_ISOBARIC, ltype2, hasBounds, level, dlevel1, dlevel2);
2666 }
2667
2668 break;
2669 }
2670 case ZAXIS_SNOW:
2671 {
2672 if (editionNumber <= 1)
2673 ; // not available
2674 else
2675 grib2DefLevel(gh, gcinit, grib_ltype, grib_ltype, hasBounds, level, dlevel1, dlevel2);
2676
2677 break;
2678 }
2679 case ZAXIS_DEPTH_BELOW_LAND:
2680 {
2681 if (editionNumber <= 1)
2682 {
2683 grib_ltype = hasBounds ? GRIB1_LTYPE_LANDDEPTH_LAYER : GRIB1_LTYPE_LANDDEPTH;
2684 double sf = zaxis_units_to_centimeter(zaxisID);
2685 grib1DefLevel(gh, gcinit, grib_ltype, hasBounds, level*sf, dlevel1*sf, dlevel2*sf);
2686 }
2687 else
2688 {
2689 double sf = zaxis_units_to_meter(zaxisID);
2690 grib2DefLevel(gh, gcinit, grib_ltype, grib_ltype, hasBounds, level*sf, dlevel1*sf, dlevel2*sf);
2691 }
2692
2693 break;
2694 }
2695 case ZAXIS_REFERENCE:
2696 {
2697 if (!gcinit) GRIB_CHECK(my_grib_set_long(gh, "genVertHeightCoords", 1), 0);
2698
2699 if (editionNumber <= 1)
2700 ; // not available
2701 else
2702 {
2703 if (hasBounds)
2704 {
2705 gribapiDefLevelType(gh, gcinit, "typeOfFirstFixedSurface", grib_ltype);
2706 gribapiDefLevelType(gh, gcinit, "typeOfSecondFixedSurface", grib_ltype2);
2707 GRIB_CHECK(my_grib_set_long(gh, "topLevel", (long) dlevel1), 0);
2708 GRIB_CHECK(my_grib_set_long(gh, "bottomLevel", (long) dlevel2), 0);
2709 }
2710 else
2711 {
2712 grib2DefLevel(gh, gcinit, grib_ltype, grib_ltype2, hasBounds, level, dlevel1, dlevel2);
2713 }
2714
2715 GRIB_CHECK(my_grib_set_long(gh, "NV", 6), 0);
2716 int number = 0;
2717 cdiInqKeyInt(zaxisID, CDI_GLOBAL, CDI_KEY_NUMBEROFVGRIDUSED, &number);
2718 GRIB_CHECK(my_grib_set_long(gh, "numberOfVGridUsed", number), 0);
2719 int nlev = 0;
2720 cdiInqKeyInt(zaxisID, CDI_GLOBAL, CDI_KEY_NLEV, &nlev);
2721 GRIB_CHECK(my_grib_set_long(gh, "nlev", nlev), 0);
2722 unsigned char uuid[CDI_UUID_SIZE];
2723 int length = CDI_UUID_SIZE;
2724 memset(uuid, 0, length);
2725 cdiInqKeyBytes(zaxisID, CDI_GLOBAL, CDI_KEY_UUID, uuid, &length);
2726 size_t len = CDI_UUID_SIZE;
2727 if ( grib_set_bytes(gh, "uuidOfVGrid", uuid, &len) != 0 ) Warning("Can't write UUID!");
2728 }
2729
2730 break;
2731 }
2732 case ZAXIS_GENERIC:
2733 {
2734 if (editionNumber <= 1)
2735 grib1DefLevel(gh, gcinit, ltype, hasBounds, level, dlevel1, dlevel2);
2736 else
2737 grib2DefLevel(gh, gcinit, ltype, ltype, hasBounds, level, dlevel1, dlevel2);
2738
2739 break;
2740 }
2741 default:
2742 {
2743 Error("Unsupported zaxis type: %s", zaxisNamePtr(zaxistype));
2744 break;
2745 }
2746 }
2747 }
2748
2749
gribapiGetScanningMode(grib_handle * gh)2750 int gribapiGetScanningMode(grib_handle *gh)
2751 {
2752 long iScansNegatively, jScansPositively, jPointsAreConsecutive;
2753 GRIB_CHECK(grib_get_long(gh, "iScansNegatively", &iScansNegatively), 0);
2754 GRIB_CHECK(grib_get_long(gh, "jScansPositively", &jScansPositively), 0);
2755 GRIB_CHECK(grib_get_long(gh, "jPointsAreConsecutive", &jPointsAreConsecutive), 0);
2756 int scanningMode
2757 = 128*(bool)iScansNegatively
2758 + 64 *(bool)jScansPositively
2759 + 32 *(bool)jPointsAreConsecutive;
2760 if (cdiDebugExt>=30)
2761 printf("gribapiGetScanningMode(): Scanning mode = %02d (%1d%1d%1d)*32; \n",\
2762 scanningMode,(int)jPointsAreConsecutive,(int)jScansPositively,(int)iScansNegatively);
2763
2764 return scanningMode;
2765 }
2766
2767
gribapiSetScanningMode(grib_handle * gh,int scanningMode)2768 void gribapiSetScanningMode(grib_handle *gh, int scanningMode)
2769 {
2770 // 127: reserved for testing; generated test data will be in 64 scanning mode
2771 //if (scanningMode== 127) scanningMode = 64;
2772
2773 const long iScansNegatively = (scanningMode & 128)/128;
2774 const long jScansPositively = (scanningMode & 64)/64;
2775 const long jPointsAreConsecutive = (scanningMode & 32)/32;
2776
2777 if (cdiDebugExt>=30)
2778 {
2779 long paramId, levelTypeId, levelId, uvRelativeToGrid;
2780 GRIB_CHECK(grib_get_long(gh, "uvRelativeToGrid", &uvRelativeToGrid), 0);
2781 GRIB_CHECK(grib_get_long(gh, "indicatorOfParameter", ¶mId), 0);
2782 GRIB_CHECK(grib_get_long(gh, "indicatorOfTypeOfLevel", &levelTypeId), 0);
2783 GRIB_CHECK(grib_get_long(gh, "level", &levelId), 0);
2784 printf("gribapiSetScanningMode(): (param,ltype,level) = (%3d,%3d,%4d); Scanning mode = %02d (%1d%1d%1d)*32; uvRelativeToGrid = %02d\n",\
2785 (int)paramId, (int)levelTypeId, (int)levelId,
2786 scanningMode,(int)jPointsAreConsecutive,(int)jScansPositively,(int)iScansNegatively,
2787 (int)uvRelativeToGrid);
2788 }
2789
2790 GRIB_CHECK(my_grib_set_long(gh, "iScansNegatively", iScansNegatively), 0);
2791 GRIB_CHECK(my_grib_set_long(gh, "jScansPositively", jScansPositively), 0);
2792 GRIB_CHECK(my_grib_set_long(gh, "jPointsAreConsecutive", jPointsAreConsecutive), 0);
2793 }
2794
2795
2796 /*
2797 TABLE 8. SCANNING MODE FLAG
2798
2799 (GDS Octet 28)
2800 BIT VALUE MEANING
2801 1 0 Points scan in +i direction
2802 1 Points scan in -i direction
2803 2 0 Points scan in -j direction
2804 1 Points scan in +j direction
2805 3 0 Adjacent points in i direction are consecutive
2806 (FORTRAN: (I,J))
2807 1 Adjacent points in j direction are consecutive
2808 (FORTRAN: (J,I))
2809
2810 => Scanning Mode 0 0 0 0 0 0 0 0 (00 dec) +i, -j; i direction consecutive (row-major order West->East & North->South)
2811 => Scanning Mode 0 1 0 0 0 0 0 0 (64 dec) +i, +j; i direction consecutive (row-major order West->East & South->North )
2812 => Scanning Mode 1 1 0 0 0 0 0 0 (96 dec) +i, +j; j direction consecutive (column-major order South->North & West->East )
2813
2814 NOTE: South->North - As if you would plot the data as image on the screen
2815 where [0,0] of the data is the top-left pixel.
2816
2817 grib2ppm LAMH_D11_201302150000_00000_oro | display ppm:-
2818 ImageMagick (display): [0,0] of an image belongs to the top-left pixel
2819 [DEFAULT] : 64 dec
2820
2821 iScansNegatively = 0;
2822 jScansPositively = 1;
2823 jPointsAreConsecutive = 0; => Scanning Mode 64
2824
2825 cdo selindexbox,1,726,100,550 LAMH_D11_201302150000_00000_oro LAMH_D11_201302150000_00000_oro_cropped
2826 grib2ppm LAMH_D11_201302150000_00000_oro_cropped | /usr/bin/display ppm:- &
2827 # ^^^ this image will be missing the souther parts of data
2828
2829 grib2ppm LAMH_D11_201302150000_00000_oro | /usr/bin/display ppm:- &
2830 # ^ full domain data
2831 */
2832
2833 #ifdef HIRLAM_EXTENSIONS
2834 static void
verticallyFlipGridDefinitionWhenScanningModeChanged(grib_handle * gh,double yfirst,double ylast,double yinc)2835 verticallyFlipGridDefinitionWhenScanningModeChanged(grib_handle *gh, double yfirst, double ylast, double yinc )
2836 {
2837 /*
2838 Nj = 550;
2839 latitudeOfFirstGridPointInDegrees = -30.8;
2840 latitudeOfLastGridPointInDegrees = 24.1;
2841 iScansNegatively = 0;
2842 jScansPositively = 0;
2843 jPointsAreConsecutive = 0;
2844 jDirectionIncrementInDegrees = 0.1;
2845
2846 When switching from scanning mode 0 <=> 64
2847 yfirst = -30.8 + (550-1)*0.1
2848
2849 yfirst = yfirst + (ysize-1) * yinc
2850 yinc = -1.0*yinc
2851 */
2852
2853 //long jDim=0;
2854 //GRIB_CHECK(grib_get_long(gh, "Nj", &jDim), 0);
2855
2856 double latitudeOfFirstGridPointInDegrees;
2857 double latitudeOfLastGridPointInDegrees;
2858 double jDirectionIncrementInDegrees;
2859
2860 //GRIB_CHECK(grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &latitudeOfFirstGridPointInDegrees), 0); // yfirst
2861 //GRIB_CHECK(grib_get_double(gh, "latitudeOfLastGridPointInDegrees", &latitudeOfLastGridPointInDegrees), 0); // ylast
2862 //GRIB_CHECK(grib_get_double(gh, "jDirectionIncrementInDegrees", &jDirectionIncrementInDegrees), 0); // yinc
2863
2864 if (cdiDebugExt>=10)
2865 {
2866 Message(" BEFORE: yfirst = %f; ylast = %f; yinc = %f; ", yfirst,ylast, yinc);
2867 }
2868
2869 GRIB_CHECK(my_grib_set_double(gh, "latitudeOfFirstGridPointInDegrees", ylast), 0);
2870 GRIB_CHECK(my_grib_set_double(gh, "latitudeOfLastGridPointInDegrees", yfirst), 0);
2871 //yinc *= -1.0; // don't set yinc here ...
2872 //GRIB_CHECK(my_grib_set_double(gh, "jDirectionIncrementInDegrees", yinc), 0);
2873
2874 if (cdiDebugExt>=10)
2875 {
2876 GRIB_CHECK(grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &latitudeOfFirstGridPointInDegrees), 0); // yfirst
2877 GRIB_CHECK(grib_get_double(gh, "latitudeOfLastGridPointInDegrees", &latitudeOfLastGridPointInDegrees), 0); // ylast
2878 GRIB_CHECK(grib_get_double(gh, "jDirectionIncrementInDegrees", &jDirectionIncrementInDegrees), 0); // yinc
2879 Message("CHANGED INTO: yfirst = %f, ylast = %f, yinc = %f",latitudeOfFirstGridPointInDegrees,latitudeOfLastGridPointInDegrees, jDirectionIncrementInDegrees);
2880 }
2881 }
2882
2883 static void
convertDataScanningMode(int scanModeIN,int scanModeOUT,double * data,size_t gridsize,size_t iDim,size_t jDim)2884 convertDataScanningMode(int scanModeIN, int scanModeOUT, double *data,
2885 size_t gridsize, size_t iDim, size_t jDim)
2886 {
2887 size_t i,j;
2888 size_t idxIN, idxOUT;
2889
2890 // 127: reserved for testing; it will generate test data in 64 scanning mode
2891 if (scanModeOUT== 127) // fill with testdata ...
2892 {
2893 scanModeOUT = 64;
2894 if (cdiDebugExt>=30) printf("convertDataScanningMode(): Generating test data in 64 scanning mode..\n");
2895 for (j=0; j<jDim; j++)
2896 {
2897 const size_t jXiDim = j*iDim;
2898 for (i=0; i<iDim; i++)
2899 {
2900 idxIN = i + jXiDim;
2901 data[idxIN] = (double) (100.0*j +i);
2902 }
2903 }
2904 }
2905
2906 if ( (iDim*jDim)!= gridsize)
2907 {
2908 if (cdiDebugExt>=30) printf("convertDataScanningMode(): ERROR: (iDim*jDim)!= gridsize; (%zu * %zu) != %zu\n", iDim,jDim, gridsize);
2909 return;
2910 }
2911 if (cdiDebugExt>=30) printf("convertDataScanningMode(): scanModeIN=%02d => scanModeOUT=%02d ; where: (iDim * jDim == gridsize) (%zu*%zu == %zu)\n",scanModeIN, scanModeOUT, iDim,jDim, gridsize);
2912
2913 if (cdiDebugExt>=100)
2914 {
2915 printf("convertDataScanningMode(): data IN:\n");
2916 for (j=0; j<jDim; j++)
2917 {
2918 size_t jXiDim = j*iDim;
2919 for (i=0; i<iDim; i++)
2920 {
2921 idxIN = i + jXiDim;
2922 printf("%03.0f, ",data[idxIN]);
2923 }
2924 printf("\n");
2925 }
2926 }
2927
2928 if (scanModeIN==scanModeOUT)
2929 {
2930 if (cdiDebugExt>=30) printf("convertDataScanningMode(): INFO: Nothing to do; scanModeIN==scanModeOUT..\n");
2931 return;
2932 }
2933
2934 if (0)
2935 {
2936 return;
2937 if (scanModeOUT==00)
2938 {
2939 if (cdiDebugExt>0) printf("convertDataScanningMode(): Leave data unchaged BUT set scanModeOUT=00.\n");
2940 // CHECK: Looks like that GRIB-API provide (no matter what) data in the scannning mode 00, even it is store in the gribfile as 64 !!
2941 return;
2942 }
2943 }
2944 double *dataCopy = (double *) malloc(gridsize*sizeof(double));
2945 memcpy((void*)dataCopy,(void*) data, gridsize*sizeof(double));
2946
2947 if (scanModeIN==64) // Scanning Mode (00 dec) +i, -j; i direction consecutive (row-major order West->East & South->North )
2948 { // Scanning Mode (64 dec) +i, +j; i direction consecutive (row-major order West->East & North->South )
2949 // Scanning Mode (96 dec) +i, +j; j direction consecutive (column-major order North->South & West->East )
2950 if (scanModeOUT==00)
2951 // CHECK: Looks like that GRIB-API provide (no matter what) data in the scannning mode 00, even it is store in the gribfile as 64 !!
2952 #define VERTICAL_FLIP
2953 #ifdef VERTICAL_FLIP
2954 { // flip the data vertically ..
2955 idxIN= 0; idxOUT= (jDim-1)*iDim;
2956 if (cdiDebugExt>=30) printf("convertDataScanningMode(): copying rows nr. (%04d : %04zu)\n",0,jDim-1);
2957 for (j=0; j<jDim; j++)
2958 {
2959 memcpy((void*)&data[idxOUT], (void*)&dataCopy[idxIN], iDim*sizeof(double));
2960 idxIN += iDim; idxOUT -= iDim;
2961 }
2962 } // end if (scanModeOUT==00)*/
2963 #endif
2964 #ifdef HORIZONTAL_FLIP
2965 { // flip data horizontally ...
2966 if (1)
2967 {
2968 if (cdiDebugExt>=30) printf("convertDataScanningMode(): copying columns nr. (%04d : %04d);\n", 0, iDim-1);
2969 for (i=0; i<iDim; i++)
2970 {
2971 for (j=0; j<jDim; j++)
2972 {
2973 size_t jXiDim = j*iDim;
2974 idxIN = i + jXiDim;
2975 //data[idxIN] = (double) (100.0*j +i); // just some testdata ..
2976 idxOUT = iDim - i -1 + jXiDim;
2977 //printf("[%03d=>%03d] = %f;",idxIN,idxOUT,dataCopy[idxIN]);
2978 data[idxOUT] = dataCopy[idxIN];
2979 }
2980 }
2981 }
2982 } // end if (scanModeOUT==00)
2983 #endif
2984
2985 if (scanModeOUT==96)
2986 { // transpose the data
2987 if (cdiDebugExt>=30) printf("convertDataScanningMode(): transpose data rows=>columns nr. (%04d : %04zu) => (%04d : %04zu);\n", 0, iDim-1, 0, jDim-1);
2988 for (j=0; j<jDim; j++)
2989 {
2990 size_t jXiDim = j*iDim;
2991 for (i=0; i<iDim; i++)
2992 {
2993 idxIN = i + jXiDim;
2994 idxOUT = j + i*jDim;
2995 //printf("[%03d=>%03d] = %f;",idxIN,idxOUT,dataCopy[idxIN]);
2996 data[idxOUT] = dataCopy[idxIN];
2997 }
2998 //printf(".\n");
2999 }
3000 } // end if (scanModeOUT==96)
3001 } // end if (scanModeIN==64)
3002
3003 if (scanModeIN==00) // Scanning Mode (00 dec) +i, -j; i direction consecutive (row-major order West->East & South->North )
3004 { // Scanning Mode (64 dec) +i, +j; i direction consecutive (row-major order West->East & North->South )
3005 // Scanning Mode (96 dec) +i, +j; j direction consecutive (column-major order North->South & West->East )
3006 if (scanModeOUT==64)
3007 { // flip the data vertically ..
3008 idxIN= 0; idxOUT= (jDim-1)*iDim;
3009 for (j=0; j<jDim; j++)
3010 {
3011 if (cdiDebugExt>=25) printf("convertDataScanningMode(): copying row nr. %04zu; [idxIN=%08zu] => [idxOUT=%08zu]\n",j, idxIN, idxOUT);
3012 memcpy((void*)&data[idxOUT], (void*)&dataCopy[idxIN], iDim*sizeof(double));
3013 idxIN += iDim; idxOUT -= iDim;
3014 }
3015 } // end if (scanModeOUT==64)
3016
3017 if (scanModeOUT==96)
3018 { // transpose the data
3019 size_t jInv;
3020 for (j=0; j<jDim; j++)
3021 {
3022 if (cdiDebugExt>=30) printf("convertDataScanningMode(): processing row nr. %04zu;\n", j);
3023 jInv = (jDim-1) -j;
3024 for (i=0; i<iDim; i++)
3025 data[j + i*jDim] = dataCopy[i + jInv*iDim]; // source data has -j
3026 }
3027 } // end if (scanModeOUT==96)
3028 } // end if (scanModeIN==00)
3029
3030 if (scanModeIN==96) // Scanning Mode (00 dec) +i, -j; i direction consecutive (row-major order West->East & South->North )
3031 { // Scanning Mode (64 dec) +i, +j; i direction consecutive (row-major order West->East & North->South )
3032 // Scanning Mode (96 dec) +i, +j; j direction consecutive (column-major order North->South & West->East )
3033 if (scanModeOUT==64)
3034 { // transpose the data
3035 for (j=0; j<jDim; j++)
3036 {
3037 if (cdiDebugExt>=30) printf("convertDataScanningMode(): processing row nr. %04zu;\n", j);
3038 size_t jXiDim = j*iDim;
3039 for (i=0; i<iDim; i++)
3040 //data[j + i*jDim] = dataCopy[i + j*iDim];
3041 data[i + jXiDim] = dataCopy[j + i*jDim];
3042 }
3043 } // end if (scanModeOUT==64)
3044
3045 if (scanModeOUT==00)
3046 { // transpose the data
3047 idxIN= 0; idxOUT= 0;
3048 size_t jInv;
3049 for (j=0; j<jDim; j++)
3050 {
3051 if (cdiDebugExt>=30) printf("convertDataScanningMode(): processing row nr. %04zu;\n", j);
3052 jInv = (jDim-1) -j;
3053 size_t jXiDim = j*iDim;
3054 for (i=0; i<iDim; i++)
3055 //data[jInv + iXjDim] = dataCopy[i + jXiDim]; // target data has -j
3056 data[i + jXiDim] = dataCopy[jInv + i*jDim]; // target data has -j
3057 }
3058 } // end if (scanModeOUT==00)
3059 } // end if (scanModeIN==96)
3060
3061 if (cdiDebugExt>=100)
3062 {
3063 printf("convertDataScanningMode(): data OUT (new scanning mode):\n");
3064 for (j=0; j<jDim; j++)
3065 {
3066 size_t jXiDim = j*iDim;
3067 for (i=0; i<iDim; i++)
3068 {
3069 idxIN = i + jXiDim;
3070 printf("%03.0f, ",data[idxIN]);
3071 }
3072 printf("\n");
3073 }
3074 }
3075
3076 free(dataCopy); return;
3077 }
3078 #endif //HIRLAM_EXTENSIONS
3079
3080 static
gribapiSetExtMode(grib_handle * gh,int gridID,size_t datasize,const void * data)3081 void gribapiSetExtMode(grib_handle *gh, int gridID, size_t datasize, const void *data)
3082 {
3083 #ifndef HIRLAM_EXTENSIONS
3084 (void)data;
3085 (void)datasize;
3086 #endif
3087 int gridtype = gridInqType(gridID);
3088 if ( gridtype == GRID_PROJECTION )
3089 {
3090 int projtype = gridInqProjType(gridID);
3091 if ( projtype == CDI_PROJ_RLL ) gridtype = GRID_LONLAT;
3092 else if ( projtype == CDI_PROJ_LCC ) gridtype = CDI_PROJ_LCC;
3093 else if ( projtype == CDI_PROJ_STERE ) gridtype = CDI_PROJ_STERE;
3094 }
3095
3096 if ( gridtype == GRID_GENERIC || gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN ||
3097 gridtype == GRID_GAUSSIAN_REDUCED || gridtype == CDI_PROJ_LCC )
3098 {
3099 #ifdef HIRLAM_EXTENSIONS
3100 int scanModeIN = 0;
3101 cdiInqKeyInt(gridID, CDI_GLOBAL, CDI_KEY_SCANNINGMODE, &scanModeIN);
3102
3103 if (cdiDebugExt>=100)
3104 {
3105 size_t gridsize = gridInqSize(gridID);
3106 Message("(scanModeIN=%d; gridsize=%zu", scanModeIN, gridsize);
3107 }
3108
3109 if ( cdiGribDataScanningMode.active ) // allowed modes: <0, 64, 96>; Default is 64
3110 {
3111 size_t iDim = gridInqXsize(gridID);
3112 size_t jDim = gridInqYsize(gridID);
3113
3114 double yfirst = gridInqYval(gridID, 0);
3115 double ylast = gridInqYval(gridID, jDim-1);
3116 double yinc = gridInqYinc(gridID);
3117
3118 int scanModeOUT = cdiGribDataScanningMode.value;
3119 convertDataScanningMode(scanModeIN, scanModeOUT, (double*)data, datasize, iDim, jDim);
3120 // This will overrule the old scanning mode of the given grid
3121 if (cdiDebugExt>=10) Message("Set GribDataScanningMode (%d) => (%d)", scanModeIN, cdiGribDataScanningMode.value);
3122 gribapiSetScanningMode(gh, cdiGribDataScanningMode.value);
3123
3124 if (((scanModeIN==00) && (cdiGribDataScanningMode.value==64)) ||
3125 ((scanModeIN==64) && (cdiGribDataScanningMode.value==00)) )
3126 verticallyFlipGridDefinitionWhenScanningModeChanged(gh, yfirst, ylast, yinc);
3127 }
3128 else
3129 {
3130 if (cdiDebugExt>=100) Message("Set GribDataScanningMode => (%d) based on used grid", scanModeIN);
3131 gribapiSetScanningMode(gh, scanModeIN);
3132 }
3133 #endif
3134 }
3135 }
3136
3137 // #define GRIBAPIENCODETEST 1
3138
gribapiEncode(int varID,int levelID,int vlistID,int gridID,int zaxisID,int64_t vdate,int vtime,int tsteptype,int numavg,size_t datasize,const void * data,size_t nmiss,void ** gribbuffer,size_t * gribbuffersize,int comptype,void * gribContainer)3139 size_t gribapiEncode(int varID, int levelID, int vlistID, int gridID, int zaxisID,
3140 int64_t vdate, int vtime, int tsteptype, int numavg,
3141 size_t datasize, const void *data, size_t nmiss, void **gribbuffer, size_t *gribbuffersize,
3142 int comptype, void *gribContainer)
3143 {
3144 void *dummy = NULL;
3145 long editionNumber = 2;
3146
3147 // extern unsigned char _grib_template_GRIB2[];
3148 cdi_check_gridsize_int_limit("GRIB", datasize);
3149
3150 int param = vlistInqVarParam(vlistID, varID);
3151 int datatype = vlistInqVarDatatype(vlistID, varID);
3152 int typeOfGeneratingProcess = 0;
3153 cdiInqKeyInt(vlistID, varID, CDI_KEY_TYPEOFGENERATINGPROCESS, &typeOfGeneratingProcess);
3154 int productDefinitionTemplate = 0;
3155 cdiInqKeyInt(vlistID, varID, CDI_KEY_PRODUCTDEFINITIONTEMPLATE, &productDefinitionTemplate);
3156
3157 int uvRelativeToGrid = -1;
3158 cdiInqKeyInt(vlistID, varID, CDI_KEY_UVRELATIVETOGRID, &uvRelativeToGrid);
3159
3160 char name[256], stdname[256];
3161 vlistInqVarName(vlistID, varID, name);
3162 vlistInqVarStdname(vlistID, varID, stdname);
3163
3164 #ifdef GRIBAPIENCODETEST
3165 grib_handle *gh = (grib_handle *) gribHandleNew(editionNumber);
3166 #else
3167 gribContainer_t *gc = (gribContainer_t *) gribContainer;
3168 assert(gc != NULL);
3169 grib_handle *gh = (struct grib_handle *) gc->gribHandle;
3170 #endif
3171
3172 GRIB_CHECK(grib_get_long(gh, "editionNumber", &editionNumber), 0);
3173
3174 if (editionNumber == 2)
3175 {
3176 if (! gc->init)
3177 {
3178 int backgroundProcess = 0;
3179 cdiInqKeyInt(vlistID, varID, CDI_KEY_BACKGROUNDPROCESS, &backgroundProcess);
3180 GRIB_CHECK(my_grib_set_long(gh, "typeOfGeneratingProcess", typeOfGeneratingProcess), 0);
3181 GRIB_CHECK(my_grib_set_long(gh, "backgroundProcess", backgroundProcess), 0);
3182 int status, tablesVersion, localTablesVersion;
3183 status = cdiInqKeyInt(vlistID, varID, CDI_KEY_TABLESVERSION, &tablesVersion);
3184 if ( status == 0 ) GRIB_CHECK(my_grib_set_long(gh, "tablesVersion", (long)tablesVersion), 0);
3185 status = cdiInqKeyInt(vlistID, varID, CDI_KEY_LOCALTABLESVERSION, &localTablesVersion);
3186 if ( status == 0 ) GRIB_CHECK(my_grib_set_long(gh, "localTablesVersion", (long)localTablesVersion), 0);
3187 int typeOfProcessedData = 0;
3188 status = cdiInqKeyInt(vlistID, varID, CDI_KEY_TYPEOFPROCESSEDDATA, &typeOfProcessedData);
3189 if ( status == 0 ) GRIB_CHECK(my_grib_set_long(gh, "typeOfProcessedData", (long)typeOfProcessedData), 0);
3190 /*
3191 int constituentType = 0;
3192 status = cdiInqKeyInt(vlistID, varID, CDI_KEY_CONSTITUENTTYPE, &constituentType);
3193 if ( status == 0 ) GRIB_CHECK(my_grib_set_long(gh, "constituentType", (long)constituentType), 0);
3194 */
3195 }
3196 }
3197
3198 gribapiDefTime((int)editionNumber, productDefinitionTemplate, typeOfGeneratingProcess,
3199 gh, vdate, vtime, tsteptype, numavg, vlistInqTaxis(vlistID), gc->init);
3200
3201 {
3202 int typeOfTimeIncrement = 0;
3203 int status = cdiInqKeyInt(vlistID, varID, CDI_KEY_TYPEOFTIMEINCREMENT, &typeOfTimeIncrement);
3204 if ( status == 0 ) grib_set_long(gh, "typeOfTimeIncrement", (long)typeOfTimeIncrement);
3205 }
3206
3207 {
3208 int status, perturbationNumber, numberOfForecastsInEnsemble, typeOfEnsembleForecast;
3209 status = cdiInqKeyInt(vlistID, varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, &typeOfEnsembleForecast);
3210 if (status == 0) grib_set_long(gh, "typeOfEnsembleForecast", typeOfEnsembleForecast);
3211 status = cdiInqKeyInt(vlistID, varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, &numberOfForecastsInEnsemble);
3212 if (status == 0) grib_set_long(gh, "numberOfForecastsInEnsemble", numberOfForecastsInEnsemble);
3213 status = cdiInqKeyInt(vlistID, varID, CDI_KEY_PERTURBATIONNUMBER, &perturbationNumber);
3214 if (status == 0) grib_set_long(gh, "perturbationNumber", perturbationNumber);
3215 }
3216
3217 if (! gc->init) gribapiDefInstitut(gh, vlistID, varID);
3218 if (! gc->init) gribapiDefModel(gh, vlistID, varID);
3219
3220 if (! gc->init) gribapiDefParam((int)editionNumber, gh, param, name, stdname);
3221
3222 if (! gc->init && editionNumber == 2)
3223 {
3224 int shapeOfTheEarth = 0;
3225 cdiInqKeyInt(vlistID, varID, CDI_KEY_SHAPEOFTHEEARTH, &shapeOfTheEarth);
3226 GRIB_CHECK(my_grib_set_long(gh, "shapeOfTheEarth", (long)shapeOfTheEarth), 0);
3227
3228 int grib2LocalSectionNumber, section2PaddingLength;
3229 int mpimType, mpimClass, mpimUser;
3230 if ( cdiInqKeyInt(vlistID, varID, CDI_KEY_MPIMTYPE, &mpimType) == CDI_NOERR &&
3231 cdiInqKeyInt(vlistID, varID, CDI_KEY_MPIMCLASS, &mpimClass) == CDI_NOERR &&
3232 cdiInqKeyInt(vlistID, varID, CDI_KEY_MPIMUSER, &mpimUser) == CDI_NOERR )
3233 {
3234 grib_set_long(gh, "grib2LocalSectionPresent", 1);
3235 grib_set_long(gh, "grib2LocalSectionNumber", 1);
3236 grib_set_long(gh, "mpimType", mpimType);
3237 grib_set_long(gh, "mpimClass", mpimClass);
3238 grib_set_long(gh, "mpimUser", mpimUser);
3239
3240 int revNumLen = 20;
3241 unsigned char revNumber[revNumLen];
3242 if ( cdiInqKeyBytes(vlistID, varID, CDI_KEY_REVNUMBER, revNumber, &revNumLen) == CDI_NOERR )
3243 {
3244 size_t revNumLenS = revNumLen;
3245 grib_set_bytes(gh, "revNumber", revNumber, &revNumLenS);
3246 }
3247 int revStatus;
3248 if ( cdiInqKeyInt(vlistID, varID, CDI_KEY_REVSTATUS, &revStatus) == CDI_NOERR )
3249 grib_set_long(gh, "revStatus", revStatus);
3250 }
3251 else if ( cdiInqKeyInt(vlistID, varID, CDI_KEY_GRIB2LOCALSECTIONNUMBER, &grib2LocalSectionNumber) == CDI_NOERR &&
3252 cdiInqKeyInt(vlistID, varID, CDI_KEY_SECTION2PADDINGLENGTH, §ion2PaddingLength) == CDI_NOERR )
3253 {
3254 grib_set_long(gh, "grib2LocalSectionPresent", 1);
3255 grib_set_long(gh, "grib2LocalSectionNumber", grib2LocalSectionNumber);
3256 unsigned char *section2Padding = (unsigned char*) Malloc(section2PaddingLength);
3257 cdiInqKeyBytes(vlistID, varID, CDI_KEY_SECTION2PADDING, section2Padding, §ion2PaddingLength);
3258 size_t len = section2PaddingLength;
3259 grib_set_bytes(gh, "section2Padding", section2Padding, &len);
3260 Free(section2Padding);
3261 }
3262 }
3263
3264 // bitsPerValue have to be defined first (complex packing)
3265 GRIB_CHECK(my_grib_set_long(gh, "bitsPerValue", (long)grbBitsPerValue(datatype)), 0);
3266
3267 if ( ! gc->init ) gribapiDefGrid((int)editionNumber, gh, gridID, comptype, datatype, uvRelativeToGrid);
3268
3269 gribapiDefLevel((int)editionNumber, gh, zaxisID, levelID, gc->init, productDefinitionTemplate);
3270
3271 vlist_t *vlistptr = vlist_to_pointer(vlistID);
3272 //if (!gc->init)
3273 {
3274 int ret = 0;
3275
3276 // NOTE: Optional key/value pairs: Note that we do not distinguish between tiles here!
3277
3278 for ( int i=0; i<vlistptr->vars[varID].opt_grib_nentries; i++ )
3279 {
3280 if ( vlistptr->vars[varID].opt_grib_kvpair[i].update )
3281 {
3282 // DR: Fix for multi-level fields (otherwise only the 1st level is correct)
3283 if ( zaxisInqSize(zaxisID)==(levelID+1) )
3284 vlistptr->vars[varID].opt_grib_kvpair[i].update = false;
3285
3286 if (vlistptr->vars[varID].opt_grib_kvpair[i].data_type == t_double)
3287 {
3288 if ( CDI_Debug )
3289 Message("key \"%s\" : double value = %g\n",
3290 vlistptr->vars[varID].opt_grib_kvpair[i].keyword,
3291 vlistptr->vars[varID].opt_grib_kvpair[i].dbl_val);
3292 my_grib_set_double(gh, vlistptr->vars[varID].opt_grib_kvpair[i].keyword,
3293 vlistptr->vars[varID].opt_grib_kvpair[i].dbl_val);
3294 GRIB_CHECK(ret, 0);
3295 }
3296 if (vlistptr->vars[varID].opt_grib_kvpair[i].data_type == t_int)
3297 {
3298 if ( CDI_Debug )
3299 Message("key \"%s\" : integer value = %d\n",
3300 vlistptr->vars[varID].opt_grib_kvpair[i].keyword,
3301 vlistptr->vars[varID].opt_grib_kvpair[i].int_val);
3302 my_grib_set_long(gh, vlistptr->vars[varID].opt_grib_kvpair[i].keyword,
3303 (long) vlistptr->vars[varID].opt_grib_kvpair[i].int_val);
3304 GRIB_CHECK(ret, 0);
3305 }
3306 }
3307 }
3308 }
3309
3310 if ( nmiss > 0 )
3311 {
3312 GRIB_CHECK(my_grib_set_long(gh, "bitmapPresent", 1), 0);
3313 GRIB_CHECK(my_grib_set_double(gh, "missingValue", vlistInqVarMissval(vlistID, varID)), 0);
3314 }
3315
3316 gribapiSetExtMode(gh, gridID, datasize, data);
3317
3318 GRIB_CHECK(grib_set_double_array(gh, "values", (const double*)data, datasize), 0);
3319
3320 // get the size of coded message
3321 size_t recsize = 0;
3322 GRIB_CHECK(grib_get_message(gh, (const void **)&dummy, &recsize), 0);
3323 recsize += 512; // add some space for possible filling
3324 *gribbuffersize = recsize;
3325 *gribbuffer = Malloc(*gribbuffersize);
3326
3327 if ( ! gc->init && editionNumber == 2 )
3328 {
3329 long pdis;
3330 grib_get_long(gh, "discipline", &pdis);
3331 if ( pdis != 255 )
3332 {
3333 char cdi_name[CDI_MAX_NAME]; cdi_name[0] = 0;
3334 vlistInqVarName(vlistID, varID, cdi_name);
3335 char grb_name[256];
3336 gribapiGetString(gh, "shortName", grb_name, sizeof(grb_name));
3337 strToLower(cdi_name);
3338 strToLower(grb_name);
3339 bool checkName = (!grb_name[0] && strncmp(cdi_name, "param", 5) == 0) ? false : true;
3340 if ( checkName && ((strlen(cdi_name) != strlen(grb_name)) || !strStartsWith(cdi_name, grb_name)) )
3341 Message("*** GRIB2 shortName does not correspond to chosen variable name: \"%s\" (\"%s\").",
3342 grb_name[0]?grb_name:"unknown", cdi_name);
3343 }
3344 }
3345
3346 // get a copy of the coded message
3347 GRIB_CHECK(grib_get_message_copy(gh, *gribbuffer, &recsize), 0);
3348
3349 #ifdef GRIBAPIENCODETEST
3350 gribHandleDelete(gh);
3351 #endif
3352
3353 gc->init = true;
3354
3355 return recsize;
3356 }
3357
3358
gribapiChangeParameterIdentification(grib_handle * gh,int code,int ltype,int level)3359 void gribapiChangeParameterIdentification(grib_handle *gh, int code, int ltype, int level)
3360 {
3361 // timeRangeIndicator: could be included later
3362 if (code!=-1) GRIB_CHECK(my_grib_set_long(gh, "indicatorOfParameter", code), 0);
3363 if (ltype!=-1) GRIB_CHECK(my_grib_set_long(gh, "indicatorOfTypeOfLevel", ltype), 0);
3364 if (level!=-1) GRIB_CHECK(my_grib_set_long(gh, "level", level), 0);
3365 }
3366
3367 #endif
3368
3369 /*
3370 * Local Variables:
3371 * c-file-style: "Java"
3372 * c-basic-offset: 2
3373 * indent-tabs-mode: nil
3374 * show-trailing-whitespace: t
3375 * require-trailing-newline: t
3376 * End:
3377 */
3378