ocssw  1.0
/disk01/web/ocssw/build/src/libhdfutils/hdf5util.cpp (r8102/r7537)
Go to the documentation of this file.
00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <time.h>
00004 #include <string.h>
00005 #include <math.h>
00006 
00007 #include "hdf5util.h"
00008 #include "hdf5_Aquarius.h"
00009 #include "passthebuck.h"
00010 
00011 #define FAIL -1
00012 
00013 #define DATACENTER "NASA/GSFC Aquarius Data Processing Center"
00014 #define MISSIONCHARACTERISTICS "Nominal orbit: inclination=98.0 (Sun-synchronous); node=6PM (ascending); eccentricity=<0.002; altitude=657 km; ground speed=6.825 km/sec"
00015 
00016 #define SENSORCHARACTERISTICS "Number of beams=3; channels per receiver=4; frequency 1.413 GHz; bits per sample=16; instatntaneous field of view=6.5 degrees; science data block period=1.44 sec."
00017 
00018 // Millisecs since start of day
00019 int32_t get_millisec( string *ydhmsf_str) {
00020   istringstream istr;
00021   int32_t itemp;
00022 
00023   istr.clear(); istr.str( ydhmsf_str->substr(7,2)); istr >> itemp;
00024   int32_t millisec = itemp * 3600000;
00025   istr.clear(); istr.str( ydhmsf_str->substr(9,2)); istr >> itemp;
00026   millisec += itemp * 60000;
00027   istr.clear(); istr.str( ydhmsf_str->substr(11,5)); istr >> itemp;
00028   millisec += itemp;
00029 
00030   return millisec;
00031 }
00032 
00033 namespace Hdf {
00034 
00035   int h5a_set(hid_t dataset, const char *nam, hid_t typ,
00036           hid_t cnt, VOIDP data) {
00037 
00038     hid_t attr, atype, wr_typ, aid;
00039     hsize_t fdim[1];
00040     herr_t  ret; 
00041 
00042     /* Save old error handler */
00043     H5E_auto_t old_func;
00044     void *old_client_data;
00045     H5Eget_auto(H5E_DEFAULT , &old_func, &old_client_data);
00046 
00047     /* Turn off error handling */
00048     H5Eset_auto( H5E_DEFAULT, NULL, NULL);
00049 
00050     /* Probe. Likely to fail, but that's okay */
00051     attr = H5Aopen_name(dataset, nam);
00052 
00053     /* Restore previous error handler */
00054     H5Eset_auto( H5E_DEFAULT, old_func, old_client_data);
00055     
00056     // Write attribute if it exists and return
00057     if (attr > 0) {
00058       if(H5Awrite(attr, typ, data)) { 
00059     fprintf(stderr,"-E- %s line %d: ",__FILE__,__LINE__);
00060     fprintf(stderr,"H5Awrite(%ld,\"%s\",%ld,%ld,data) failed.  ",
00061         (long) dataset, nam, (long) typ, (long) cnt);
00062     return(HDF_FUNCTION_ERROR);
00063       }
00064       ret = H5Aclose(attr);
00065       return(LIFE_IS_GOOD); 
00066     }
00067 
00068 
00069     if (typ == H5T_STRING) {
00070       /*
00071        * Create string attribute.
00072        */
00073       aid  = H5Screate(H5S_SCALAR);
00074 
00075       atype = H5Tcopy(H5T_C_S1);
00076       H5Tset_size(atype, cnt);
00077 
00078       attr = H5Acreate1(dataset, nam, atype, aid, H5P_DEFAULT);
00079       wr_typ = atype;
00080 
00081     } else if (cnt == 1) {
00082       /*
00083        * Create scalar attribute.
00084        */
00085       aid  = H5Screate(H5S_SCALAR);
00086       attr = H5Acreate1(dataset, nam, typ, aid, H5P_DEFAULT);
00087       wr_typ = typ;
00088     
00089     } else {
00090       aid = H5Screate(H5S_SIMPLE);
00091       fdim[0] = cnt;
00092       ret  = H5Sset_extent_simple(aid, 1, fdim, NULL);
00093     
00094       /*
00095        * Create array attribute.
00096        */
00097       attr = H5Acreate1(dataset, nam, typ, aid, H5P_DEFAULT);
00098       wr_typ = typ;
00099     }
00100 
00101     /*
00102      * Write attribute.
00103      */
00104     if(H5Awrite(attr, wr_typ, data)) { 
00105       fprintf(stderr,"-E- %s line %d: ",__FILE__,__LINE__);
00106       fprintf(stderr,"H5Awrite(%ld,\"%s\",%ld,%ld,data) failed.  ",
00107           (long) dataset, nam, (long) typ, (long) cnt);
00108       return(HDF_FUNCTION_ERROR);
00109     }
00110 
00111     /*
00112      * Close attribute and file dataspaces.
00113      */
00114     ret = H5Sclose(aid);
00115     
00116     /*
00117      * Close the attributes.
00118      */ 
00119     ret = H5Aclose(attr);
00120     
00121     return(LIFE_IS_GOOD);
00122   }
00123 
00124 
00125   hid_t h5d_create(hid_t grp, const char *nam, hid_t typ,
00126            int rank, 
00127            hsize_t d0, hsize_t d1, hsize_t d2, 
00128            hsize_t d3, hsize_t d4, hsize_t d5, 
00129            hid_t *dataset, 
00130            hid_t *dataspace,
00131            hid_t plist) {
00132 
00133     hsize_t dimsizes[6];
00134     herr_t ret;
00135 
00136     if(rank > 6){
00137       fprintf(stderr,"-E- %s line %d: ",__FILE__,__LINE__);
00138       fprintf(stderr,"sd_create() expects to be passed a rank <= 6.  ");
00139       return(PROGRAMMER_BOOBOO);
00140     }
00141     dimsizes[0]=d0; dimsizes[1]=d1; dimsizes[2]=d2;
00142     dimsizes[3]=d3; dimsizes[4]=d4; dimsizes[5]=d5;
00143 
00144     *dataspace = H5Screate(H5S_SIMPLE);
00145     ret = H5Sset_extent_simple(*dataspace, rank, dimsizes, NULL);
00146 
00147     if((*dataset = H5Dcreate1(grp, nam, typ, 
00148                  *dataspace, plist)) == FAIL){
00149       fprintf(stderr,"-E- %s line %d: ",__FILE__,__LINE__);
00150       fprintf(stderr,"H5Dcreate(%ld,\"%s\",%ld,%ld,[%ld,%ld,%ld]) failed.  ",
00151           (long) grp, nam, (long) typ, (long) rank, 
00152           (long) d0, (long) d1, (long) d2);
00153       return(HDF_FUNCTION_ERROR);
00154     }
00155     
00156     return(LIFE_IS_GOOD);
00157   }
00158 
00159 
00160   int SetScalarH5A(hid_t id, const char *name, hid_t type, const void *value){
00161     if (type == H5T_STRING)
00162       h5a_set(id, name, type, strlen((char *) value)+1, (VOIDP) value);
00163     else
00164       h5a_set(id, name, type, 1, (VOIDP) value);
00165 
00166     return(LIFE_IS_GOOD);
00167   }
00168 
00169 
00170 
00171   /* ----------------------------------------------------- */
00172   /* Create an H5D using the wrappers for the HDF routines */
00173   /* ----------------------------------------------------- */
00174   int CreateH5D(
00175         hid_t grp,      /* group id */
00176         const char *sname,    /* short name */
00177         const char *lname,    /* long name */
00178         const char *units,    /* units (not set if passed NULL or "") */
00179         double low,     /* low end of valid range */
00180         double high,    /* high end of range (no range set if low >= high) */
00181         float slope,    /* scale factor (not set if 0)  */
00182         float offset,   /* scaling offset (not set if 0)  */
00183         hid_t nt,       /* HDF number type */
00184         int   rank,     /* number of dimensions (must be <= 3) */
00185         int32_t  d0,       /* size of 1st dimension */
00186         int32_t  d1,       /* size of 2nd dimension */
00187         int32_t  d2,       /* size of 3rd dimension */
00188         int32_t  d3,       /* size of 4th dimension */
00189         int32_t  d4,       /* size of 5th dimension */
00190         int32_t  d5,       /* size of 6th dimension */
00191         const char *dn0,      /* name of 1st dimension */
00192         const char *dn1,      /* name of 2nd dimension */
00193         const char *dn2,      /* name of 3rd dimension */
00194         const char *dn3,      /* name of 4th dimension */
00195         const char *dn4,      /* name of 5th dimension */
00196         const char *dn5,      /* name of 6th dimension */
00197         hid_t plist     /* Dataset property list */
00198         ) {
00199     
00200     hid_t dataset,dataspace;
00201 
00202     /* Create the SDS */
00203     PTB( h5d_create(grp, sname, nt, rank, d0, d1, d2, d3, d4, d5,
00204             &dataset, &dataspace, plist) );
00205 
00206     /* Add a "long_name" attribute */
00207     PTB( SetScalarH5A (dataset, "long_name", (hid_t) H5T_STRING, lname)   );
00208     /* Add a "valid_range" attribute if one is specified */
00209     if (nt == H5T_NATIVE_UCHAR) {
00210       unsigned char vr[2];
00211       vr[0] = (unsigned char)low;
00212       vr[1] = (unsigned char)high;
00213       PTB(h5a_set(dataset, "valid_min", H5T_NATIVE_UCHAR, 1, (VOIDP) &vr[0]));
00214       PTB(h5a_set(dataset, "valid_max", H5T_NATIVE_UCHAR, 1, (VOIDP) &vr[1]));
00215     } else if (nt == H5T_NATIVE_USHORT) {
00216       short vr[2];
00217       vr[0] = (short)low;
00218       vr[1] = (short)high;
00219       PTB(h5a_set(dataset, "valid_min", H5T_NATIVE_USHORT, 1, (VOIDP) &vr[0]));
00220       PTB(h5a_set(dataset, "valid_max", H5T_NATIVE_USHORT, 1, (VOIDP) &vr[1]));
00221     } else if (nt == H5T_STD_I32LE) {
00222       int32_t vr[2];
00223       vr[0] = (int32_t)low;
00224       vr[1] = (int32_t)high;
00225       PTB(h5a_set(dataset, "valid_min", H5T_STD_I32LE, 1, (VOIDP) &vr[0]));
00226       PTB(h5a_set(dataset, "valid_max", H5T_STD_I32LE, 1, (VOIDP) &vr[1]));
00227     } else if (nt == H5T_STD_U32LE) {
00228       int32_t vr[2];
00229       vr[0] = (int32_t)low;
00230       vr[1] = (int32_t)high;
00231       PTB(h5a_set(dataset, "valid_min", H5T_STD_U32LE, 1, (VOIDP) &vr[0]));
00232       PTB(h5a_set(dataset, "valid_max", H5T_STD_U32LE, 1, (VOIDP) &vr[1]));
00233     } else if (nt == H5T_NATIVE_FLOAT) {
00234       float vr[2];
00235       vr[0] = (float)low;
00236       vr[1] = (float)high;
00237       PTB(h5a_set(dataset, "valid_min", H5T_NATIVE_FLOAT, 1, (VOIDP) &vr[0]));
00238       PTB(h5a_set(dataset, "valid_max", H5T_NATIVE_FLOAT, 1, (VOIDP) &vr[1]));
00239     } else if (nt == H5T_NATIVE_DOUBLE) {
00240       double vr[2];
00241       vr[0] = (float)low;
00242       vr[1] = (float)high;
00243       PTB(h5a_set(dataset, "valid_min", H5T_NATIVE_DOUBLE, 1, (VOIDP) &vr[0]));
00244       PTB(h5a_set(dataset, "valid_max", H5T_NATIVE_DOUBLE, 1, (VOIDP) &vr[1]));
00245     } else {
00246       fprintf(stderr,"-E- %s line %d: ",__FILE__,__LINE__);
00247       fprintf(stderr,"Got unsupported number type (%ld) ", (long) nt);
00248       fprintf(stderr,"while trying to create SDS, \"%s\", ",sname);
00249       exit(1);
00250     }
00251 
00252 
00253     /* Add a "slope" attribute if one is specified, and also
00254        an intercept attribute */
00255     if (slope != 0) {
00256       float vr[1], intercept[1];
00257       vr[0] = slope;
00258       intercept[0] = 0.0;
00259       
00260       PTB( h5a_set(dataset, "slope", H5T_NATIVE_FLOAT, 1, (VOIDP) vr) );
00261     }
00262     if (offset != 0) {
00263       float intercept[1];
00264       intercept[0] = offset;
00265       PTB( h5a_set(dataset, "intercept", H5T_NATIVE_FLOAT, 1, 
00266            (VOIDP) intercept) );
00267     }
00268 
00269     /* Add a "units" attribute if one is specified */
00270     if(units != NULL && *units != 0) {
00271       PTB( SetScalarH5A (dataset, "units", (hid_t) H5T_STRING, units) );
00272     }
00273 
00274 
00275     /* Release this SDS */
00276     PTB( H5Sclose(dataspace) );
00277     PTB( H5Dclose(dataset) );
00278 
00279     return(LIFE_IS_GOOD);
00280   }
00281 
00282 
00283 
00284   /* ----------------------------------------------------- */
00285   /* Create an H5D using the wrappers for the HDF routines */
00286   /* ----------------------------------------------------- */
00287   int CreateH5D(
00288         hid_t grp,      /* group id */
00289         const char *sname,    /* short name */
00290         const char *lname,    /* long name */
00291         const char *stdname,    /* standard name */
00292         const char *units,    /* units (not set if passed NULL or "") */
00293         double low,     /* low end of valid range */
00294         double high,    /* high end of range (no range set if low >= high) */
00295         float slope,    /* scale factor (not set if 0)  */
00296         float offset,   /* scaling offset (not set if 0)  */
00297         float fillvalue,/* fill value */
00298         hid_t nt,       /* HDF number type */
00299         int   rank,     /* number of dimensions (must be <= 3) */
00300         int32_t  d0,       /* size of 1st dimension */
00301         int32_t  d1,       /* size of 2nd dimension */
00302         int32_t  d2,       /* size of 3rd dimension */
00303         int32_t  d3,       /* size of 4th dimension */
00304         int32_t  d4,       /* size of 5th dimension */
00305         int32_t  d5,       /* size of 6th dimension */
00306         const char *dn0,      /* name of 1st dimension */
00307         const char *dn1,      /* name of 2nd dimension */
00308         const char *dn2,      /* name of 3rd dimension */
00309         const char *dn3,      /* name of 4th dimension */
00310         const char *dn4,      /* name of 5th dimension */
00311         const char *dn5,      /* name of 6th dimension */
00312         hid_t plist     /* Dataset property list */
00313         ) {
00314     
00315     hid_t dataset;
00316 
00317     CreateH5D( grp, sname, lname, units, low, high, slope, offset,
00318            nt, rank, d0, d1, d2, d3, d4, d5, 
00319            dn0, dn1, dn2, dn3, dn4, dn5, plist);
00320 
00321     dataset = H5Dopen1( grp, sname);
00322 
00323     if (nt == H5T_NATIVE_SHORT) {
00324       short _FillValue;
00325       _FillValue = (short) fillvalue;
00326       PTB(h5a_set(dataset, " _FillValue", H5T_NATIVE_SHORT, 1, (VOIDP) &_FillValue));
00327     } else if (nt == H5T_STD_I32LE) {
00328       int32_t _FillValue;
00329       _FillValue = (int32_t) fillvalue;
00330       PTB(h5a_set(dataset, " _FillValue", H5T_STD_I32LE, 1, (VOIDP) &_FillValue));
00331     } else if (nt == H5T_NATIVE_FLOAT) {
00332       float _FillValue;
00333       _FillValue = (float) fillvalue;
00334       PTB(h5a_set(dataset, " _FillValue", H5T_NATIVE_FLOAT, 1, (VOIDP) &_FillValue));
00335     } else if (nt == H5T_NATIVE_DOUBLE) {
00336       double _FillValue;
00337       _FillValue = (float) fillvalue;
00338       PTB(h5a_set(dataset, " _FillValue", H5T_NATIVE_DOUBLE, 1, (VOIDP) &_FillValue));
00339     } else {
00340       fprintf(stderr,"-E- %s line %d: ",__FILE__,__LINE__);
00341       fprintf(stderr,"Got unsupported number type (%ld) ", (long) nt);
00342       fprintf(stderr,"while trying to create SDS, \"%s\", ",sname);
00343       return(PROGRAMMER_BOOBOO);
00344     }
00345 
00346 
00347     // Standard Name
00348     if ( strlen(stdname) != 0) {
00349       PTB( SetScalarH5A (dataset, "standard_name", (hid_t) H5T_STRING, stdname)   );
00350     }
00351 
00352     /* Release this SDS */
00353     PTB( H5Dclose(dataset) );
00354 
00355     return(LIFE_IS_GOOD);
00356   }
00357 
00358 
00359 
00360   /*---------------------- */
00361   /* Write to HDF5 dataset */
00362   /* --------------------- */
00363   herr_t h5d_write(hid_t id, const char *name, VOIDP data, hsize_t rank,
00364            hsize_t s[6], hsize_t e[6])
00365   {
00366 
00367     hid_t dataset, dataspace, filespace, datatype;
00368     hsize_t start[6], edge[6], dims[6];
00369     herr_t status;
00370 
00371     for (size_t i=0; i<6; i++) {
00372       start[i] = s[i];
00373       edge[i]  = e[i];
00374     }
00375 
00376 
00377     dataset = H5Dopen1(id, name);
00378     if (dataset == -1) {
00379       printf("Dataset: %s not found (%d)\n", name, id);
00380       exit(1);
00381     }
00382 
00383     /* Get datatype */
00384     /* ------------ */
00385     datatype = H5Dget_type(dataset);
00386 
00387     /*
00388      * Select a hyperslab.
00389      */
00390     filespace = H5Dget_space(dataset);
00391 
00392 
00393     // Check for out-of-bounds write
00394     status = H5Sget_simple_extent_dims(filespace, dims, NULL);
00395     for (size_t i=0; i<rank; i++) {
00396       if (start[i]+edge[i] > dims[i]) {
00397     printf("Write to dim (0-based): %d of \"%s\" out of bounds (%d %d)\n", 
00398            (int) i, name, (int) (start[i]+edge[i]), (int) dims[i]);
00399     exit (110);
00400       }
00401     }
00402 
00403     status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, start, NULL,
00404                  edge, NULL);  
00405 
00406     dataspace = H5Screate_simple(rank, edge, NULL); 
00407 
00408     /*
00409      * Write the data to the hyperslab.
00410      */
00411     status = H5Dwrite(dataset, datatype, dataspace, filespace,
00412               H5P_DEFAULT, data);
00413 
00414     H5Dclose(dataset);
00415     H5Sclose(dataspace);
00416     H5Sclose(filespace);
00417 
00418     return(LIFE_IS_GOOD);
00419   }
00420 
00421 
00422 
00423   /*----------------------- */
00424   /* Read from HDF5 dataset */
00425   /* ---------------------- */
00426   herr_t h5d_read(hid_t id, const char *name, VOIDP data, hsize_t rank,
00427           hsize_t s[6], hsize_t e[6])
00428   {
00429 
00430     hid_t dataset, dataspace, filespace, datatype;
00431     hsize_t start[6], edge[6], dims[6];
00432     herr_t status;
00433 
00434     dataset = H5Dopen1(id, name);
00435     if (dataset == -1) {
00436       printf("Dataset: %s not found (%d)\n", name, id);
00437       exit(1);
00438     }
00439 
00440     // Get datatype
00441     datatype = H5Dget_type(dataset);
00442 
00443     // Select hyperslab
00444     filespace = H5Dget_space(dataset);
00445 
00446     // Return dimensions if data = NULL
00447     if (data == NULL) {
00448       for (size_t i=0; i<6; i++) e[i] = 0;
00449       status = H5Sget_simple_extent_dims(filespace, e, NULL);
00450       H5Dclose(dataset);
00451       H5Sclose(filespace);
00452       return(LIFE_IS_GOOD);
00453     }
00454 
00455     // Copy start and edge
00456     for (size_t i=0; i<6; i++) {
00457       start[i] = s[i];
00458       edge[i]  = e[i];
00459     }
00460 
00461     // Check for out-of-bounds read
00462     status = H5Sget_simple_extent_dims(filespace, dims, NULL);
00463     for (size_t i=0; i<rank; i++) {
00464       if (start[i]+edge[i] > dims[i]) {
00465     printf("Read to dim (0-based): %d of \"%s\" out of bounds (%d %d)\n", 
00466            (int) i, name, (int) (start[i]+edge[i]), (int) dims[i]);
00467     exit (110);
00468       }
00469     }
00470 
00471     status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, start, NULL,
00472                  edge, NULL);  
00473 
00474     dataspace = H5Screate_simple(rank, edge, NULL); 
00475 
00476     /*
00477      * Read the data from the hyperslab.
00478      */
00479     status = H5Dread(dataset, datatype, dataspace, filespace,
00480              H5P_DEFAULT, data);
00481 
00482     H5Dclose(dataset);
00483     H5Sclose(dataspace);
00484     H5Sclose(filespace);
00485 
00486     return(LIFE_IS_GOOD);
00487   }
00488   
00489 } // namespace Hdf