1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Author: Uwe Schulzweida
5
6 */
7
8 #ifdef HAVE_CONFIG_H
9 #include "config.h"
10 #endif
11
12 #define H5_USE_16_API
13
14 #ifdef HAVE_LIBHDF5
15 #include "hdf5.h"
16 #endif
17
18 #include <cdi.h>
19 #include "cdo_options.h"
20
21 #include "cdo_output.h"
22 #include "griddes.h"
23
24 #ifdef HAVE_LIBHDF5
25 static herr_t
obj_info(hid_t loc_id,const char * name,void * objname)26 obj_info(hid_t loc_id, const char *name, void *objname)
27 {
28 herr_t lexist = 0;
29
30 H5G_stat_t statbuf;
31 H5Gget_objinfo(loc_id, name, false, &statbuf);
32
33 if (strcmp(name, (char *) objname) == 0)
34 {
35 lexist = 1;
36
37 H5G_obj_t obj_type = statbuf.type;
38
39 switch (obj_type)
40 {
41 case H5G_GROUP:
42 if (Options::cdoVerbose) cdo_print("HDF5 object '%s' is a group", name);
43 break;
44 case H5G_DATASET:
45 if (Options::cdoVerbose) cdo_print("HDF5 object '%s' is a dataset", name);
46 break;
47 case H5G_TYPE:
48 if (Options::cdoVerbose) cdo_print("HDF5 object '%s' is a named datatype", name);
49 break;
50 default:
51 /*cdo_abort(" Unable to identify an object %s", name);*/
52 break;
53 }
54 }
55
56 return lexist;
57 }
58 #endif
59
60 #ifdef HAVE_LIBHDF5
61 static int
h5find_object(hid_t file_id,const char * name)62 h5find_object(hid_t file_id, const char *name)
63 {
64 return (int) H5Giterate(file_id, "/", nullptr, obj_info, (void *) name);
65 }
66 #endif
67
68 #ifdef HAVE_LIBHDF5
69 static void
fill_gridvals(size_t xsize,size_t ysize,double * xvals,double * yvals)70 fill_gridvals(size_t xsize, size_t ysize, double *xvals, double *yvals)
71 {
72 size_t i, j, ii, jj;
73 size_t index, index2;
74
75 double xmin = -180;
76 double xmax = 180;
77 double ymin = -90;
78 double ymax = 90;
79
80 for (ii = 0; ii < xsize / 2; ++ii)
81 {
82 index2 = ysize / 2 * xsize + ii;
83 if (xvals[index2] > -180 && xvals[index2] < 360)
84 {
85 xmin = xvals[index2];
86 break;
87 }
88 }
89
90 for (ii = xsize - 1; ii > xsize / 2; --ii)
91 {
92 index2 = ysize / 2 * xsize + ii;
93 if (xvals[index2] > -180 && xvals[index2] < 360)
94 {
95 xmax = xvals[index2];
96 break;
97 }
98 }
99 /*
100 for ( jj = 0; jj < ysize; ++jj )
101 {
102 index2 = jj*xsize + xsize/2;
103 if ( xvals[index2] < -180 || xvals[index2] > 360 ) xvals[index2] = 0;
104 index2 = jj*xsize + xsize/2-1;
105 if ( xvals[index2] < -180 || xvals[index2] > 360 ) xvals[index2] = 0;
106 }
107 */
108 for (jj = 0; jj < ysize / 2; ++jj)
109 {
110 index2 = jj * xsize + xsize / 2;
111 if (yvals[index2] > -90 && yvals[index2] < 90)
112 {
113 ymax = yvals[index2];
114 break;
115 }
116 }
117
118 for (jj = ysize - 1; jj > ysize / 2; --jj)
119 {
120 index2 = jj * xsize + xsize / 2;
121 if (yvals[index2] > -90 && yvals[index2] < 90)
122 {
123 ymin = yvals[index2];
124 break;
125 }
126 }
127
128 /* printf("xmin %g, xmax %g, ymin %g, ymax %g\n", xmin, xmax, ymin, ymax); */
129
130 for (i = 0; i < xsize * ysize; ++i)
131 {
132 if (xvals[i] > -180 && xvals[i] < 360)
133 {
134 if (xvals[i] < xmin) xmin = xvals[i];
135 if (xvals[i] > xmax) xmax = xvals[i];
136 }
137
138 if (yvals[i] > -90 && yvals[i] < 90)
139 {
140 if (yvals[i] < ymin) ymin = yvals[i];
141 if (yvals[i] > ymax) ymax = yvals[i];
142 }
143 }
144
145 for (j = 0; j < ysize; ++j)
146 for (i = 0; i < xsize; ++i)
147 {
148 index = j * xsize + i;
149
150 if (xvals[index] < -180 || xvals[index] > 360)
151 {
152 if (i < xsize / 2)
153 xvals[index] = xmin;
154 else
155 xvals[index] = xmax;
156 /*
157 if ( j < ysize/2 )
158 for ( jj = j+1; jj < ysize/2; ++jj )
159 {
160 index2 = jj*xsize + i;
161 if ( xvals[index2] > -180 && xvals[index2] < 360 )
162 {
163 xvals[index] = xvals[index2];
164 break;
165 }
166 }
167 else
168 for ( jj = j-1; jj > ysize/2; --jj )
169 {
170 index2 = jj*xsize + i;
171 if ( xvals[index2] > -180 && xvals[index2] < 360 )
172 {
173 xvals[index] = xvals[index2];
174 break;
175 }
176 }
177 */
178 /*
179 if ( i < xsize/2 )
180 {
181 xvals[index] = xmin;
182 for ( ii = i+1; ii < xsize/2; ++ii )
183 {
184 index2 = j*xsize + ii;
185 if ( xvals[index2] > -180 && xvals[index2] < 360 )
186 {
187 xvals[index] = (xmin*(ii-i) + xvals[index2]*(i))/ii;
188 break;
189 }
190 }
191 }
192 else
193 {
194 for ( ii = i-1; ii >= xsize/2; --ii )
195 {
196 index2 = j*xsize + ii;
197 if ( xvals[index2] > -180 && xvals[index2] < 360 )
198 {
199 xvals[index] = (xmax*(i-ii) + xvals[index2]*((xsize-1)-i))/(xsize-1-ii); break;
200 }
201 }
202 }
203 */
204 }
205
206 if (yvals[index] < -90 || yvals[index] > 90)
207 {
208 if (j < ysize / 2)
209 yvals[index] = ymax;
210 else
211 yvals[index] = ymin;
212
213 if (i < xsize / 2)
214 for (ii = i + 1; ii < xsize / 2; ++ii)
215 {
216 index2 = j * xsize + ii;
217 if (yvals[index2] > -90 && yvals[index2] < 90)
218 {
219 yvals[index] = yvals[index2];
220 break;
221 }
222 }
223 else
224 for (ii = i - 1; ii > xsize / 2; --ii)
225 {
226 index2 = j * xsize + ii;
227 if (yvals[index2] > -90 && yvals[index2] < 90)
228 {
229 yvals[index] = yvals[index2];
230 break;
231 }
232 }
233 }
234 }
235 }
236
237 int
grid_from_h5_file(const char * gridfile)238 grid_from_h5_file(const char *gridfile)
239 {
240 int gridID = -1;
241 hid_t lon_id = -1; /* Dataset ID */
242 hid_t lat_id = -1; /* Dataset ID */
243 hid_t att_id;
244 hid_t dataspace;
245 hsize_t dims_out[9]; /* dataset dimensions */
246 herr_t status = 0; /* Generic return value */
247 int rank;
248 GridDesciption grid;
249
250 hid_t fapl_id = H5Pcreate(H5P_FILE_ACCESS);
251 H5Pset_fclose_degree(fapl_id, H5F_CLOSE_STRONG);
252
253 /* Open an existing file. */
254 hid_t file_id = H5Fopen(gridfile, H5F_ACC_RDONLY, fapl_id);
255
256 H5Pclose(fapl_id);
257
258 if (file_id < 0) return gridID;
259
260 if (h5find_object(file_id, "lon") > 0 && h5find_object(file_id, "lat") > 0)
261 {
262 lon_id = H5Dopen(file_id, "/lon");
263 lat_id = H5Dopen(file_id, "/lat");
264 }
265 else if (h5find_object(file_id, "Longitude") > 0 && h5find_object(file_id, "Latitude") > 0)
266 {
267 lon_id = H5Dopen(file_id, "/Longitude");
268 lat_id = H5Dopen(file_id, "/Latitude");
269 }
270
271 if (lon_id >= 0 && lat_id >= 0)
272 {
273 dataspace = H5Dget_space(lon_id); /* dataspace handle */
274 rank = H5Sget_simple_extent_ndims(dataspace);
275 status = H5Sget_simple_extent_dims(dataspace, dims_out, nullptr);
276
277 if (rank != 2)
278 {
279 // if ( Options::cdoVerbose ) cdo_warning("Unexpected rank = %d!", rank);
280 goto RETURN;
281 }
282
283 // check for netcdf4 attribute
284 if (H5Aexists(lon_id, "DIMENSION_LIST")) goto RETURN;
285 if (H5Aexists(lat_id, "DIMENSION_LIST")) goto RETURN;
286
287 if (H5Aexists(lon_id, "bounds")) goto RETURN;
288 if (H5Aexists(lat_id, "bounds")) goto RETURN;
289
290 /*
291 printf("\nRank: %d\nDimensions: %lu x %lu \n", rank,
292 (unsigned long)(dims_out[1]), (unsigned long)(dims_out[0]));
293 */
294
295 hid_t type_id = H5Dget_type(lon_id); /* get datatype*/
296
297 hid_t native_type = H5Tget_native_type(type_id, H5T_DIR_ASCEND);
298 int ftype = 0;
299 if (H5Tequal(native_type, H5T_NATIVE_SCHAR) > 0)
300 {
301 ftype = 0;
302 }
303 else if (H5Tequal(native_type, H5T_NATIVE_UCHAR) > 0)
304 {
305 ftype = 0;
306 }
307 else if (H5Tequal(native_type, H5T_NATIVE_SHORT) > 0)
308 {
309 ftype = 0;
310 }
311 else if (H5Tequal(native_type, H5T_NATIVE_USHORT) > 0)
312 {
313 ftype = 0;
314 }
315 else if (H5Tequal(native_type, H5T_NATIVE_INT) > 0)
316 {
317 ftype = 0;
318 }
319 else if (H5Tequal(native_type, H5T_NATIVE_UINT) > 0)
320 {
321 ftype = 0;
322 }
323 else if (H5Tequal(native_type, H5T_NATIVE_FLOAT) > 0)
324 {
325 ftype = 1;
326 }
327 else if (H5Tequal(native_type, H5T_NATIVE_DOUBLE) > 0)
328 {
329 ftype = 1;
330 }
331 else
332 {
333 cdo_warning("Grid has unsupported native datatype!");
334 goto RETURN;
335 }
336 H5Tclose(native_type);
337
338 grid.xsize = dims_out[1];
339 grid.ysize = dims_out[0];
340 grid.size = grid.xsize * grid.ysize;
341
342 grid.xvals.resize(grid.size);
343 grid.yvals.resize(grid.size);
344
345 if (ftype)
346 {
347 status = H5Dread(lon_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, grid.xvals.data());
348 status = H5Dread(lat_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, grid.yvals.data());
349 }
350 else
351 {
352 std::vector<int> iarray(grid.size);
353 status = H5Dread(lon_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, iarray.data());
354 for (size_t i = 0; i < grid.size; ++i) grid.xvals[i] = iarray[i];
355 status = H5Dread(lat_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, iarray.data());
356 for (size_t i = 0; i < grid.size; ++i) grid.yvals[i] = iarray[i];
357 }
358
359 status = H5Sclose(dataspace);
360
361 /* Close the dataset. */
362 status = H5Dclose(lon_id);
363 status = H5Dclose(lat_id);
364
365 fill_gridvals(grid.xsize, grid.ysize, grid.xvals.data(), grid.yvals.data());
366
367 grid.type = GRID_CURVILINEAR;
368 grid.datatype = CDI_DATATYPE_FLT32;
369
370 gridID = grid_define(grid);
371 }
372 else if (h5find_object(file_id, "where") > 0)
373 {
374 double xscale = 1, yscale = 1;
375 double xoffset = 0, yoffset = 0;
376 hid_t grp_id;
377
378 grp_id = H5Gopen(file_id, "/where/lon/what");
379 if (grp_id >= 0)
380 {
381 att_id = H5Aopen_name(grp_id, "gain");
382 if (att_id >= 0)
383 {
384 status = H5Aread(att_id, H5T_NATIVE_DOUBLE, &xscale);
385 H5Aclose(att_id);
386 }
387
388 att_id = H5Aopen_name(grp_id, "offset");
389 if (att_id >= 0)
390 {
391 status = H5Aread(att_id, H5T_NATIVE_DOUBLE, &xoffset);
392 H5Aclose(att_id);
393 }
394
395 H5Gclose(grp_id);
396 }
397
398 grp_id = H5Gopen(file_id, "/where/lat/what");
399 if (grp_id >= 0)
400 {
401 att_id = H5Aopen_name(grp_id, "gain");
402 if (att_id >= 0)
403 {
404 status = H5Aread(att_id, H5T_NATIVE_DOUBLE, &yscale);
405 H5Aclose(att_id);
406 }
407
408 att_id = H5Aopen_name(grp_id, "offset");
409 if (att_id >= 0)
410 {
411 status = H5Aread(att_id, H5T_NATIVE_DOUBLE, &yoffset);
412 H5Aclose(att_id);
413 }
414
415 H5Gclose(grp_id);
416 }
417
418 /* Open an existing dataset. */
419 lon_id = H5Dopen(file_id, "/where/lon/data");
420 if (lon_id >= 0) lat_id = H5Dopen(file_id, "/where/lat/data");
421
422 if (lon_id >= 0 && lat_id >= 0)
423 {
424 dataspace = H5Dget_space(lon_id); /* dataspace handle */
425 rank = H5Sget_simple_extent_ndims(dataspace);
426 status = H5Sget_simple_extent_dims(dataspace, dims_out, nullptr);
427
428 if (rank != 2)
429 {
430 // if ( Options::cdoVerbose ) cdo_warning("Unexpected rank = %d!", rank);
431 goto RETURN;
432 }
433 /*
434 printf("\nRank: %d\nDimensions: %lu x %lu \n", rank, (unsigned long)(dims_out[1]), (unsigned long)(dims_out[0]));
435 */
436
437 hid_t type_id = H5Dget_type(lon_id); /* get datatype*/
438
439 hid_t native_type = H5Tget_native_type(type_id, H5T_DIR_ASCEND);
440 int ftype = 0;
441 if (H5Tequal(native_type, H5T_NATIVE_SCHAR) > 0)
442 {
443 ftype = 0;
444 }
445 else if (H5Tequal(native_type, H5T_NATIVE_UCHAR) > 0)
446 {
447 ftype = 0;
448 }
449 else if (H5Tequal(native_type, H5T_NATIVE_SHORT) > 0)
450 {
451 ftype = 0;
452 }
453 else if (H5Tequal(native_type, H5T_NATIVE_USHORT) > 0)
454 {
455 ftype = 0;
456 }
457 else if (H5Tequal(native_type, H5T_NATIVE_INT) > 0)
458 {
459 ftype = 0;
460 }
461 else if (H5Tequal(native_type, H5T_NATIVE_UINT) > 0)
462 {
463 ftype = 0;
464 }
465 else if (H5Tequal(native_type, H5T_NATIVE_FLOAT) > 0)
466 {
467 ftype = 1;
468 }
469 else if (H5Tequal(native_type, H5T_NATIVE_DOUBLE) > 0)
470 {
471 ftype = 1;
472 }
473 else
474 {
475 cdo_warning("Grid has unsupported native datatype!");
476 goto RETURN;
477 }
478 H5Tclose(native_type);
479
480 grid.xsize = dims_out[1];
481 grid.ysize = dims_out[0];
482 grid.size = grid.xsize * grid.ysize;
483
484 grid.xvals.resize(grid.size);
485 grid.yvals.resize(grid.size);
486
487 if (ftype)
488 {
489 status = H5Dread(lon_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, grid.xvals.data());
490 status = H5Dread(lat_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, grid.yvals.data());
491 }
492 else
493 {
494 std::vector<int> iarray(grid.size);
495 status = H5Dread(lon_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, iarray.data());
496 for (size_t i = 0; i < grid.size; ++i) grid.xvals[i] = iarray[i];
497 status = H5Dread(lat_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, iarray.data());
498 for (size_t i = 0; i < grid.size; ++i) grid.yvals[i] = iarray[i];
499 }
500
501 status = H5Sclose(dataspace);
502
503 /* Close the dataset. */
504 status = H5Dclose(lon_id);
505 status = H5Dclose(lat_id);
506
507 for (size_t i = 0; i < grid.size; ++i) grid.xvals[i] = grid.xvals[i] * xscale + xoffset;
508 for (size_t i = 0; i < grid.size; ++i) grid.yvals[i] = grid.yvals[i] * yscale + yoffset;
509
510 grid.type = GRID_CURVILINEAR;
511 grid.datatype = CDI_DATATYPE_FLT32;
512
513 gridID = grid_define(grid);
514 }
515 }
516
517 RETURN:
518
519 /* Close file */
520 if (file_id >= 0) status = H5Fclose(file_id);
521 (void) status;
522
523 if (gridID != -1 && Options::cdoVerbose) cdo_print("%s: grid created.", __func__);
524
525 return gridID;
526 }
527 #else
528 int
grid_from_h5_file(const char * gridfile)529 grid_from_h5_file(const char *gridfile)
530 {
531 (void) gridfile;
532 cdo_warning("HDF5 support not compiled in!");
533 return -1;
534 }
535 #endif
536