bes  Updated for version 3.20.13
HDFCFUtil.cc
1 #include "config.h"
2 #include "config_hdf.h"
3 
4 #include "HDFCFUtil.h"
5 #include <BESDebug.h>
6 #include <BESLog.h>
7 #include <math.h>
8 #include"dodsutil.h"
9 #include"HDFSPArray_RealField.h"
10 #include"HDFSPArrayMissField.h"
11 #include"HDFEOS2GeoCFProj.h"
12 #include"HDFEOS2GeoCF1D.h"
13 
14 #include <libdap/debug.h>
15 
16 #define SIGNED_BYTE_TO_INT32 1
17 
18 // HDF datatype headers for both the default and the CF options
19 #include "HDFByte.h"
20 #include "HDFInt16.h"
21 #include "HDFUInt16.h"
22 #include "HDFInt32.h"
23 #include "HDFUInt32.h"
24 #include "HDFFloat32.h"
25 #include "HDFFloat64.h"
26 #include "HDFStr.h"
27 #include "HDF4RequestHandler.h"
28 //
29 using namespace std;
30 using namespace libdap;
31 
32 #define ERR_LOC1(x) #x
33 #define ERR_LOC2(x) ERR_LOC1(x)
34 #define ERR_LOC __FILE__ " : " ERR_LOC2(__LINE__)
35 
36 // Check the BES key.
37 // This function will check a BES key specified at the file h4.conf.in.
38 // If the key's value is either true or yes. The handler claims to find
39 // a key and will do some operations. Otherwise, will do different operations.
40 // For example, One may find a line H4.EnableCF=true at h4.conf.in.
41 // That means, the HDF4 handler will handle the HDF4 files by following CF conventions.
42 #if 0
43 bool
44 HDFCFUtil::check_beskeys(const string& key) {
45 
46  bool found = false;
47  string doset ="";
48  const string dosettrue ="true";
49  const string dosetyes = "yes";
50 
51  TheBESKeys::TheKeys()->get_value( key, doset, found ) ;
52  if( true == found ) {
53  doset = BESUtil::lowercase( doset ) ;
54  if( dosettrue == doset || dosetyes == doset )
55  return true;
56  }
57  return false;
58 
59 }
60 #endif
61 
66 static void
67 split_helper(vector<string> &tokens, const string &text, const char sep)
68 {
69  string::size_type start = 0;
70  string::size_type end = 0;
71 
72  while ((end = text.find(sep, start)) != string::npos) {
73  tokens.push_back(text.substr(start, end - start));
74  start = end + 1;
75  }
76  tokens.push_back(text.substr(start));
77 }
78 
79 // From a string separated by a separator to a list of string,
80 // for example, split "ab,c" to {"ab","c"}
81 void
82 HDFCFUtil::Split(const char *s, int len, char sep, std::vector<std::string> &names)
83 {
84  names.clear();
85  split_helper(names, string(s, len), sep);
86 #if 0
87  // Replaced with the above since valgrind reports errors
88  // with this code. jhrg 6/22/15
89  for (int i = 0, j = 0; j <= len; ++j) {
90  if ((j == len && len) || s[j] == sep) {
91  string elem(s + i, j - i);
92  names.push_back(elem);
93  i = j + 1;
94  continue;
95  }
96  }
97 #endif
98 }
99 
100 // Assume sz is Null terminated string.
101 void
102 HDFCFUtil::Split(const char *sz, char sep, std::vector<std::string> &names)
103 {
104  names.clear();
105 
106  // Split() was showing up in some valgrind runs as having an off-by-one
107  // error. I added this help track it down.
108  DBG(std::cerr << "HDFCFUtil::Split: sz: <" << sz << ">, sep: <" << sep << ">" << std::endl);
109 
110  split_helper(names, string(sz), sep);
111 #if 0
112  // Replaced with a direct call to the new helper code.
113  // jhrg 6/22/15
114  Split(sz, (int)strlen(sz), sep, names);
115 #endif
116 }
117 
118 #if 0
119 // From a string separated by a separator to a list of string,
120 // for example, split "ab,c" to {"ab","c"}
121 void
122 HDFCFUtil::Split(const char *s, int len, char sep, std::vector<std::string> &names)
123 {
124  names.clear();
125  for (int i = 0, j = 0; j <= len; ++j) {
126  if ((j == len && len) || s[j] == sep) {
127  string elem(s + i, j - i);
128  names.push_back(elem);
129  i = j + 1;
130  continue;
131  }
132  }
133 }
134 
135 // Assume sz is Null terminated string.
136 void
137 HDFCFUtil::Split(const char *sz, char sep, std::vector<std::string> &names)
138 {
139  Split(sz, (int)strlen(sz), sep, names);
140 }
141 
142 #endif
143 
144 // This is a safer way to insert and update a c++ map value.
145 // Otherwise, the local testsuite at The HDF Group will fail for HDF-EOS2 data
146 // under iMac machine platform.
147 // The implementation replaces the element even if the key exists.
148 // This function is equivalent to map[key]=value
149 bool
150 HDFCFUtil::insert_map(std::map<std::string,std::string>& m, string key, string val)
151 {
152  pair<map<string,string>::iterator, bool> ret;
153  ret = m.insert(make_pair(key, val));
154  if(ret.second == false){
155  m.erase(key);
156  ret = m.insert(make_pair(key, val));
157  if(ret.second == false){
158  BESDEBUG("h4","insert_map():insertion failed on Key=" << key << " Val=" << val << endl);
159  }
160  }
161  return ret.second;
162 }
163 
164 // Obtain CF string
165 string
167 {
168 
169  if(""==s)
170  return s;
171  string insertString(1,'_');
172 
173  // Always start with _ if the first character is not a letter
174  if (true == isdigit(s[0]))
175  s.insert(0,insertString);
176 
177  // New name conventions drop the first '/' from the path.
178  if ('/' ==s[0])
179  s.erase(0,1);
180 
181  for(auto &si: s)
182  if((false == isalnum(si)) && (si!='_'))
183  si='_';
184 
185  return s;
186 
187 }
188 
189 // Obtain the unique name for the clashed names and save it to set namelist.
190 // This is a recursive call. A unique name list is represented as a set.
191 // If we find that a name already exists in the nameset, we will add a number
192 // at the end of the name to form a new name. If the new name still exists
193 // in the nameset, we will increase the index number and check again until
194 // a unique name is generated.
195 void
196 HDFCFUtil::gen_unique_name(string &str,set<string>& namelist, int&clash_index) {
197 
198  pair<set<string>::iterator,bool> ret;
199  string newstr = "";
200  stringstream sclash_index;
201  sclash_index << clash_index;
202  newstr = str + sclash_index.str();
203 
204  ret = namelist.insert(newstr);
205  if (false == ret.second) {
206  clash_index++;
207  gen_unique_name(str,namelist,clash_index);
208  }
209  else
210  str = newstr;
211 }
212 
213 // General routine to handle the name clashing
214 // The input parameters include:
215 // name vector -> newobjnamelist(The name vector is changed to a unique name list
216 // a pre-allocated object name set ->objnameset(can be used to determine if a name exists)
217 void
218 HDFCFUtil::Handle_NameClashing(vector<string>&newobjnamelist,set<string>&objnameset) {
219 
220  pair<set<string>::iterator,bool> setret;
221  set<string>::iterator iss;
222 
223  vector<string> clashnamelist;
224 
225  // clash index to original index mapping
226  map<int,int> cl_to_ol;
227  int ol_index = 0;
228  int cl_index = 0;
229 
230 
231  for (const auto &newobjname:newobjnamelist) {
232  setret = objnameset.insert(newobjname);
233  if (false == setret.second ) {
234  clashnamelist.insert(clashnamelist.end(),newobjname);
235  cl_to_ol[cl_index] = ol_index;
236  cl_index++;
237  }
238  ol_index++;
239  }
240 
241  // Now change the clashed elements to unique elements,
242  // Generate the set which has the same size as the original vector.
243  for (auto &clashname:clashnamelist) {
244  int clash_index = 1;
245  string temp_clashname = clashname +'_';
246  HDFCFUtil::gen_unique_name(temp_clashname,objnameset,clash_index);
247  clashname = temp_clashname;
248  }
249 
250  // Now go back to the original vector, make it unique.
251  for (unsigned int i =0; i <clashnamelist.size(); i++)
252  newobjnamelist[cl_to_ol[i]] = clashnamelist[i];
253 
254 }
255 
256 // General routine to handle the name clashing
257 // The input parameter just includes:
258 // name vector -> newobjnamelist(The name vector is changed to a unique name list
259 void
260 HDFCFUtil::Handle_NameClashing(vector<string>&newobjnamelist) {
261 
262  set<string> objnameset;
263  Handle_NameClashing(newobjnamelist,objnameset);
264 }
265 
266 // Borrowed codes from ncdas.cc in netcdf_handle4 module.
267 string
268 HDFCFUtil::print_attr(int32 type, int loc, void *vals)
269 {
270  ostringstream rep;
271 
272  union {
273  char *cp;
274  unsigned char *ucp;
275  short *sp;
276  unsigned short *usp;
277  int32 /*nclong*/ *lp;
278  unsigned int *ui;
279  float *fp;
280  double *dp;
281  } gp;
282 
283  switch (type) {
284 
285  // Mapping both DFNT_UINT8 and DFNT_INT8 to unsigned char
286  // may cause overflow. Documented at jira ticket HFRHANDLER-169.
287  case DFNT_UINT8:
288  {
289  unsigned char uc;
290  gp.ucp = (unsigned char *) vals;
291 
292  uc = *(gp.ucp+loc);
293  rep << (int)uc;
294  return rep.str();
295  }
296 
297  case DFNT_INT8:
298  {
299  char c;
300  gp.cp = (char *) vals;
301 
302  c = *(gp.cp+loc);
303  rep << (int)c;
304  return rep.str();
305  }
306 
307  case DFNT_UCHAR:
308  case DFNT_CHAR:
309  {
310  // Use the customized escattr function. Don't escape \n,\t and \r. KY 2013-10-14
311  return escattr(static_cast<const char*>(vals));
312  }
313 
314  case DFNT_INT16:
315  {
316  gp.sp = (short *) vals;
317  rep << *(gp.sp+loc);
318  return rep.str();
319  }
320 
321  case DFNT_UINT16:
322  {
323  gp.usp = (unsigned short *) vals;
324  rep << *(gp.usp+loc);
325  return rep.str();
326  }
327 
328  case DFNT_INT32:
329  {
330  gp.lp = (int32 *) vals;
331  rep << *(gp.lp+loc);
332  return rep.str();
333  }
334 
335  case DFNT_UINT32:
336  {
337  gp.ui = (unsigned int *) vals;
338  rep << *(gp.ui+loc);
339  return rep.str();
340  }
341 
342  case DFNT_FLOAT:
343  {
344  float attr_val = *(float*)vals;
345  bool is_a_fin = isfinite(attr_val);
346  gp.fp = (float *) vals;
347  rep << showpoint;
348  // setprecision seeme to cause the one-bit error when
349  // converting from float to string. Watch whether this
350  // is an isue.
351  rep << setprecision(10);
352  rep << *(gp.fp+loc);
353  string tmp_rep_str = rep.str();
354  if (tmp_rep_str.find('.') == string::npos
355  && tmp_rep_str.find('e') == string::npos
356  && tmp_rep_str.find('E') == string::npos) {
357  if(true == is_a_fin)
358  rep << ".";
359  }
360  return rep.str();
361  }
362 
363  case DFNT_DOUBLE:
364  {
365 
366  double attr_val = *(double*)vals;
367  bool is_a_fin = isfinite(attr_val);
368  gp.dp = (double *) vals;
369  rep << std::showpoint;
370  rep << std::setprecision(17);
371  rep << *(gp.dp+loc);
372  string tmp_rep_str = rep.str();
373  if (tmp_rep_str.find('.') == string::npos
374  && tmp_rep_str.find('e') == string::npos
375  && tmp_rep_str.find('E') == string::npos) {
376  if(true == is_a_fin)
377  rep << ".";
378  }
379  return rep.str();
380  }
381  default:
382  return string("UNKNOWN");
383  }
384 
385 }
386 
387 // Print datatype in string. This is used to generate DAS.
388 string
390 {
391 
392  // I expanded the list based on libdap/AttrTable.h.
393  static const char UNKNOWN[]="Unknown";
394  static const char BYTE[]="Byte";
395  static const char INT16[]="Int16";
396  static const char UINT16[]="UInt16";
397  static const char INT32[]="Int32";
398  static const char UINT32[]="UInt32";
399  static const char FLOAT32[]="Float32";
400  static const char FLOAT64[]="Float64";
401  static const char STRING[]="String";
402 
403  // I got different cases from hntdefs.h.
404  switch (type) {
405 
406  case DFNT_CHAR:
407  return STRING;
408 
409  case DFNT_UCHAR8:
410  return STRING;
411 
412  case DFNT_UINT8:
413  return BYTE;
414 
415  case DFNT_INT8:
416 // ADD the macro
417  {
418 #ifndef SIGNED_BYTE_TO_INT32
419  return BYTE;
420 #else
421  return INT32;
422 #endif
423  }
424  case DFNT_UINT16:
425  return UINT16;
426 
427  case DFNT_INT16:
428  return INT16;
429 
430  case DFNT_INT32:
431  return INT32;
432 
433  case DFNT_UINT32:
434  return UINT32;
435 
436  case DFNT_FLOAT:
437  return FLOAT32;
438 
439  case DFNT_DOUBLE:
440  return FLOAT64;
441 
442  default:
443  return UNKNOWN;
444  }
445 
446 }
447 
448 // Obtain HDF4 datatype size.
449 short
451 {
452 
453  // .
454  switch (type) {
455 
456  case DFNT_CHAR:
457  return sizeof(char);
458 
459  case DFNT_UCHAR8:
460  return sizeof(unsigned char);
461 
462  case DFNT_UINT8:
463  return sizeof(unsigned char);
464 
465  case DFNT_INT8:
466 // ADD the macro
467  {
468 #ifndef SIGNED_BYTE_TO_INT32
469  return sizeof(char);
470 #else
471  return sizeof(int);
472 #endif
473  }
474  case DFNT_UINT16:
475  return sizeof(unsigned short);
476 
477  case DFNT_INT16:
478  return sizeof(short);
479 
480  case DFNT_INT32:
481  return sizeof(int);
482 
483  case DFNT_UINT32:
484  return sizeof(unsigned int);
485 
486  case DFNT_FLOAT:
487  return sizeof(float);
488 
489  case DFNT_DOUBLE:
490  return sizeof(double);
491 
492  default:
493  return -1;
494  }
495 
496 }
497 // Subset of latitude and longitude to follow the parameters from the DAP expression constraint
498 // Somehow this function doesn't work. Now it is not used. Still keep it here for the future investigation.
499 template < typename T >
500 void HDFCFUtil::LatLon2DSubset (T * outlatlon,
501  int majordim,
502  int minordim,
503  T * latlon,
504  const int32 * offset,
505  const int32 * count,
506  const int32 * step)
507 {
508 
509 #if 0
510  // float64 templatlon[majordim][minordim];
511  // --std=c++11 on OSX causes 'typeof' to fail. This is a GNU gcc-specific
512  // keyword. jhrg 3/28/19
513  T (*templatlonptr)[majordim][minordim] = (typeof templatlonptr) latlon;
514 #endif
515  T (*templatlonptr)[majordim][minordim] = (T *) latlon;
516  int i;
517  int j;
518  int k;
519 
520  // do subsetting
521  // Find the correct index
522  int dim0count = count[0];
523  int dim1count = count[1];
524  int dim0index[dim0count], dim1index[dim1count];
525 
526  for (i = 0; i < count[0]; i++) // count[0] is the least changing dimension
527  dim0index[i] = offset[0] + i * step[0];
528 
529 
530  for (j = 0; j < count[1]; j++)
531  dim1index[j] = offset[1] + j * step[1];
532 
533  // Now assign the subsetting data
534  k = 0;
535 
536  for (i = 0; i < count[0]; i++) {
537  for (j = 0; j < count[1]; j++) {
538  outlatlon[k] = (*templatlonptr)[dim0index[i]][dim1index[j]];
539  k++;
540 
541  }
542  }
543 }
544 
545 // CF requires the _FillValue attribute datatype is the same as the corresponding field datatype.
546 // For some NASA files, this is not true.
547 // So we need to check if the _FillValue's datatype is the same as the attribute's.
548 // If not, we need to correct them.
549 void HDFCFUtil::correct_fvalue_type(AttrTable *at,int32 dtype) {
550 
551  AttrTable::Attr_iter it = at->attr_begin();
552  bool find_fvalue = false;
553  while (it!=at->attr_end() && false==find_fvalue) {
554  if (at->get_name(it) =="_FillValue")
555  {
556  find_fvalue = true;
557  string fillvalue ="";
558  string fillvalue_type ="";
559  if((*at->get_attr_vector(it)).size() !=1)
560  throw InternalErr(__FILE__,__LINE__,"The number of _FillValue must be 1.");
561  fillvalue = (*at->get_attr_vector(it)->begin());
562  fillvalue_type = at->get_type(it);
563  string var_type = HDFCFUtil::print_type(dtype);
564 
565  if(fillvalue_type != var_type){
566 
567  at->del_attr("_FillValue");
568 
569  if (fillvalue_type == "String") {
570 
571  // String fillvalue is always represented as /+octal numbers when its type is forced to
572  // change to string(check HFRHANDLER-223). So we have to retrieve it back.
573  if(fillvalue.size() >1) {
574 
575  long int fillvalue_int = 0;
576  vector<char> fillvalue_temp(fillvalue.size());
577  char *pEnd;
578  fillvalue_int = strtol((fillvalue.substr(1)).c_str(),&pEnd,8);
579  stringstream convert_str;
580  convert_str << fillvalue_int;
581  at->append_attr("_FillValue",var_type,convert_str.str());
582  }
583  else {
584 
585  // If somehow the fillvalue type is DFNT_CHAR or DFNT_UCHAR, and the size is 1,
586  // that means the fillvalue type is wrongly defined, we treat as a 8-bit integer number.
587  // Note, the code will only assume the value ranges from 0 to 128.(JIRA HFRHANDLER-248).
588  // KY 2014-04-24
589 
590  short fillvalue_int = fillvalue.at(0);
591 
592  stringstream convert_str;
593  convert_str << fillvalue_int;
594  if(fillvalue_int <0 || fillvalue_int >128)
595  throw InternalErr(__FILE__,__LINE__,
596  "If the fillvalue is a char type, the value must be between 0 and 128.");
597 
598 
599  at->append_attr("_FillValue",var_type,convert_str.str());
600  }
601  }
602 
603  else
604  at->append_attr("_FillValue",var_type,fillvalue);
605  }
606  }
607  it++;
608  }
609 
610 }
611 
612 // CF requires the scale_factor and add_offset attribute datatypes must be the same as the corresponding field datatype.
613 // For some NASA files, this is not true.
614 // So we need to check if the _FillValue's datatype is the same as the attribute's.
615 // If not, we need to correct them.
617 
618  AttrTable::Attr_iter it = at->attr_begin();
619  bool find_scale = false;
620  bool find_offset = false;
621 
622  // Declare scale_factor,add_offset, fillvalue and valid_range type in string format.
623  string scale_factor_type;
624  string add_offset_type;
625  string scale_factor_value;
626  string add_offset_value;
627 
628  while (it!=at->attr_end() &&((find_scale!=true) ||(find_offset!=true))) {
629  if (at->get_name(it) =="scale_factor")
630  {
631  find_scale = true;
632  scale_factor_value = (*at->get_attr_vector(it)->begin());
633  scale_factor_type = at->get_type(it);
634  }
635 
636  if(at->get_name(it)=="add_offset")
637  {
638  find_offset = true;
639  add_offset_value = (*at->get_attr_vector(it)->begin());
640  add_offset_type = at->get_type(it);
641  }
642 
643  it++;
644  }
645 
646  // Change offset type to the scale type
647  if((true==find_scale) && (true==find_offset)) {
648  if(scale_factor_type != add_offset_type) {
649  at->del_attr("add_offset");
650  at->append_attr("add_offset",scale_factor_type,add_offset_value);
651  }
652  }
653 }
654 
655 
656 #ifdef USE_HDFEOS2_LIB
657 
658 // For MODIS (confirmed by level 1B) products, values between 65500(MIN_NON_SCALE_SPECIAL_VALUE)
659 // and 65535(MAX_NON_SCALE_SPECIAL_VALUE) are treated as
660 // special values. These values represent non-physical data values caused by various failures.
661 // For example, 65533 represents "when Detector is saturated".
662 bool HDFCFUtil::is_special_value(int32 dtype, float fillvalue, float realvalue) {
663 
664  bool ret_value = false;
665 
666  if (DFNT_UINT16 == dtype) {
667 
668  auto fillvalue_int = (int)fillvalue;
669  if (MAX_NON_SCALE_SPECIAL_VALUE == fillvalue_int) {
670  auto realvalue_int = (int)realvalue;
671  if (realvalue_int <= MAX_NON_SCALE_SPECIAL_VALUE && realvalue_int >=MIN_NON_SCALE_SPECIAL_VALUE)
672  ret_value = true;
673  }
674  }
675 
676  return ret_value;
677 
678 }
679 
680 // Check if the MODIS file has dimension map and return the number of dimension maps
681 // Note: This routine is only applied to a MODIS geo-location file when the corresponding
682 // MODIS swath uses dimension maps and has the MODIS geo-location file under the same
683 // file directory. This is also restricted by turning on H4.EnableCheckMODISGeoFile to be true(file h4.conf.in).
684 // By default, this key is turned off. Also this function is only used once in one service. So it won't
685 // affect performance. KY 2014-02-18
686 int HDFCFUtil::check_geofile_dimmap(const string & geofilename) {
687 
688  int32 fileid = SWopen(const_cast<char*>(geofilename.c_str()),DFACC_READ);
689  if (fileid < 0)
690  return -1;
691  string swathname = "MODIS_Swath_Type_GEO";
692  int32 datasetid = SWattach(fileid,const_cast<char*>(swathname.c_str()));
693  if (datasetid < 0) {
694  SWclose(fileid);
695  return -1;
696  }
697 
698  int32 nummaps = 0;
699  int32 bufsize;
700 
701  // Obtain number of dimension maps and the buffer size.
702  if ((nummaps = SWnentries(datasetid, HDFE_NENTMAP, &bufsize)) == -1) {
703  SWdetach(datasetid);
704  SWclose(fileid);
705  return -1;
706  }
707 
708  SWdetach(datasetid);
709  SWclose(fileid);
710  return nummaps;
711 
712 }
713 
714 // Check if we need to change the datatype for MODIS fields. The datatype needs to be changed
715 // mainly because of non-CF scale and offset rules. To avoid violating CF conventions, we apply
716 // the non-CF MODIS scale and offset rule to MODIS data. So the final data type may be different
717 // than the original one due to this operation. For example, the original datatype may be int16.
718 // After applying the scale/offset rule, the datatype may become float32.
719 // The following are useful notes about MODIS SCALE OFFSET HANDLING
720 // Note: MODIS Scale and offset handling needs to re-organized. But it may take big efforts.
721 // Instead, I remove the global variable mtype, and _das; move the old calculate_dtype code
722 // back to HDFEOS2.cc. The code is a little better organized. If possible, we may think to overhaul
723 // the handling of MODIS scale-offset part. KY 2012-6-19
724 //
725 bool HDFCFUtil::change_data_type(DAS & das, SOType scaletype, const string &new_field_name)
726 {
727 
728  AttrTable *at = das.get_table(new_field_name);
729 
730  // The following codes do these:
731  // For MODIS level 1B(using the swath name), check radiance,reflectance etc.
732  // If the DISABLESCALE key is true, no need to check scale and offset for other types.
733  // Otherwise, continue checking the scale and offset names.
734  // KY 2013-12-13
735 
736  if(scaletype!=DEFAULT_CF_EQU && at!=nullptr)
737  {
738  AttrTable::Attr_iter it = at->attr_begin();
739  string scale_factor_value="";
740  string add_offset_value="0";
741  string radiance_scales_value="";
742  string radiance_offsets_value="";
743  string reflectance_scales_value="";
744  string reflectance_offsets_value="";
745  string scale_factor_type;
746  string add_offset_type;
747 
748  while (it!=at->attr_end())
749  {
750  if(at->get_name(it)=="radiance_scales")
751  radiance_scales_value = *(at->get_attr_vector(it)->begin());
752  if(at->get_name(it)=="radiance_offsets")
753  radiance_offsets_value = *(at->get_attr_vector(it)->begin());
754  if(at->get_name(it)=="reflectance_scales")
755  reflectance_scales_value = *(at->get_attr_vector(it)->begin());
756  if(at->get_name(it)=="reflectance_offsets")
757  reflectance_offsets_value = *(at->get_attr_vector(it)->begin());
758 
759  // Ideally may just check if the attribute name is scale_factor.
760  // But don't know if some products use attribut name like "scale_factor "
761  // So now just filter out the attribute scale_factor_err. KY 2012-09-20
762  if(at->get_name(it).find("scale_factor")!=string::npos){
763  string temp_attr_name = at->get_name(it);
764  if (temp_attr_name != "scale_factor_err") {
765  scale_factor_value = *(at->get_attr_vector(it)->begin());
766  scale_factor_type = at->get_type(it);
767  }
768  }
769  if(at->get_name(it).find("add_offset")!=string::npos)
770  {
771  string temp_attr_name = at->get_name(it);
772  if (temp_attr_name !="add_offset_err") {
773  add_offset_value = *(at->get_attr_vector(it)->begin());
774  add_offset_type = at->get_type(it);
775  }
776  }
777  it++;
778  }
779 
780  if((radiance_scales_value.length()!=0 && radiance_offsets_value.length()!=0)
781  || (reflectance_scales_value.length()!=0 && reflectance_offsets_value.length()!=0))
782  return true;
783 
784  if(scale_factor_value.length()!=0)
785  {
786  if(!(atof(scale_factor_value.c_str())==1 && atof(add_offset_value.c_str())==0))
787  return true;
788  }
789  }
790 
791  return false;
792 }
793 
794 // Obtain the MODIS swath dimension map info.
795 void HDFCFUtil::obtain_dimmap_info(const string& filename,HDFEOS2::Dataset*dataset,
796  vector<struct dimmap_entry> & dimmaps,
797  string & modis_geofilename, bool& geofile_has_dimmap) {
798 
799 
800  auto sw = static_cast<HDFEOS2::SwathDataset *>(dataset);
801  const vector<HDFEOS2::SwathDataset::DimensionMap*>& origdimmaps = sw->getDimensionMaps();
802  struct dimmap_entry tempdimmap;
803 
804  // if having dimension maps, we need to retrieve the dimension map info.
805  for(const auto & origdimmap:origdimmaps){
806  tempdimmap.geodim = origdimmap->getGeoDimension();
807  tempdimmap.datadim = origdimmap->getDataDimension();
808  tempdimmap.offset = origdimmap->getOffset();
809  tempdimmap.inc = origdimmap->getIncrement();
810  dimmaps.push_back(tempdimmap);
811  }
812 
813 #if 0
814  string check_modis_geofile_key ="H4.EnableCheckMODISGeoFile";
815  bool check_geofile_key = false;
816  check_geofile_key = HDFCFUtil::check_beskeys(check_modis_geofile_key);
817 #endif
818 
819  // Only when there is dimension map, we need to consider the additional MODIS geolocation files.
820  // Will check if the check modis_geo_location file key is turned on.
821  if((origdimmaps.empty() != false) && (true == HDF4RequestHandler::get_enable_check_modis_geo_file()) ) {
822 
823  // Has to use C-style since basename and dirname are not C++ routines.
824  char*tempcstr;
825  tempcstr = new char [filename.size()+1];
826  strncpy (tempcstr,filename.c_str(),filename.size());
827  string basefilename = basename(tempcstr);
828  string dirfilename = dirname(tempcstr);
829  delete [] tempcstr;
830 
831  // If the current file is a MOD03 or a MYD03 file, we don't need to check the extra MODIS geolocation file at all.
832  bool is_modis_geofile = false;
833  if(basefilename.size() >5) {
834  if((0 == basefilename.compare(0,5,"MOD03")) || (0 == basefilename.compare(0,5,"MYD03")))
835  is_modis_geofile = true;
836  }
837 
838  // This part is implemented specifically for supporting MODIS dimension map data.
839  // MODIS Aqua Swath dimension map geolocation file always starts with MYD03
840  // MODIS Terra Swath dimension map geolocation file always starts with MOD03
841  // We will check the first three characters to see if the dimension map geolocation file exists.
842  // An example MODIS swath filename is MOD05_L2.A2008120.0000.005.2008121182723.hdf
843  // An example MODIS geo-location file name is MOD03.A2008120.0000.005.2010003235220.hdf
844  // Notice that the "A2008120.0000" in the middle of the name is the "Acquistion Date" of the data
845  // So the geo-location file name should have exactly the same string. We will use this
846  // string to identify if a MODIS geo-location file exists or not.
847  // Note the string size is 14(.A2008120.0000).
848  // More information of naming convention, check http://modis-atmos.gsfc.nasa.gov/products_filename.html
849  // KY 2010-5-10
850 
851 
852  // Obtain string "MYD" or "MOD"
853  // Here we need to consider when MOD03 or MYD03 use the dimension map.
854  if ((false == is_modis_geofile) && (basefilename.size() >3)) {
855 
856  string fnameprefix = basefilename.substr(0,3);
857 
858  if(fnameprefix == "MYD" || fnameprefix =="MOD") {
859  size_t fnamemidpos = basefilename.find(".A");
860  if(fnamemidpos != string::npos) {
861  string fnamemiddle = basefilename.substr(fnamemidpos,14);
862  if(fnamemiddle.size()==14) {
863  string geofnameprefix = fnameprefix+"03";
864  // geofnamefp will be something like "MOD03.A2008120.0000"
865  string geofnamefp = geofnameprefix + fnamemiddle;
866  DIR *dirp;
867  struct dirent* dirs;
868 
869  // Go through the directory to see if we have the corresponding MODIS geolocation file
870  dirp = opendir(dirfilename.c_str());
871  if (nullptr == dirp)
872  throw InternalErr(__FILE__,__LINE__,"opendir fails.");
873 
874  while ((dirs = readdir(dirp))!= nullptr){
875  if(strncmp(dirs->d_name,geofnamefp.c_str(),geofnamefp.size())==0){
876  modis_geofilename = dirfilename + "/"+ dirs->d_name;
877  int num_dimmap = HDFCFUtil::check_geofile_dimmap(modis_geofilename);
878  if (num_dimmap < 0) {
879  closedir(dirp);
880  throw InternalErr(__FILE__,__LINE__,"this file is not a MODIS geolocation file.");
881  }
882  geofile_has_dimmap = (num_dimmap >0)?true:false;
883  break;
884  }
885  }
886  closedir(dirp);
887  }
888  }
889  }
890  }
891  }
892 }
893 
894 // This is for the case that the separate MODIS geo-location file is used.
895 // Some geolocation names at the MODIS data file are not consistent with
896 // the names in the MODIS geo-location file. So need to correct them.
897 bool HDFCFUtil::is_modis_dimmap_nonll_field(string & fieldname) {
898 
899  bool modis_dimmap_nonll_field = false;
900  vector<string> modis_dimmap_nonll_fieldlist;
901 
902  modis_dimmap_nonll_fieldlist.emplace_back("Height");
903  modis_dimmap_nonll_fieldlist.emplace_back("SensorZenith");
904  modis_dimmap_nonll_fieldlist.emplace_back("SensorAzimuth");
905  modis_dimmap_nonll_fieldlist.emplace_back("Range");
906  modis_dimmap_nonll_fieldlist.emplace_back("SolarZenith");
907  modis_dimmap_nonll_fieldlist.emplace_back("SolarAzimuth");
908  modis_dimmap_nonll_fieldlist.emplace_back("Land/SeaMask");
909  modis_dimmap_nonll_fieldlist.emplace_back("gflags");
910  modis_dimmap_nonll_fieldlist.emplace_back("Solar_Zenith");
911  modis_dimmap_nonll_fieldlist.emplace_back("Solar_Azimuth");
912  modis_dimmap_nonll_fieldlist.emplace_back("Sensor_Azimuth");
913  modis_dimmap_nonll_fieldlist.emplace_back("Sensor_Zenith");
914 
915  map<string,string>modis_field_to_geofile_field;
916  map<string,string>::iterator itmap;
917  modis_field_to_geofile_field["Solar_Zenith"] = "SolarZenith";
918  modis_field_to_geofile_field["Solar_Azimuth"] = "SolarAzimuth";
919  modis_field_to_geofile_field["Sensor_Zenith"] = "SensorZenith";
920  modis_field_to_geofile_field["Solar_Azimuth"] = "SolarAzimuth";
921 
922  for (const auto & modis_dimmap_nonll_f:modis_dimmap_nonll_fieldlist) {
923 
924  if (fieldname == modis_dimmap_nonll_f) {
925  itmap = modis_field_to_geofile_field.find(fieldname);
926  if (itmap !=modis_field_to_geofile_field.end())
927  fieldname = itmap->second;
928  modis_dimmap_nonll_field = true;
929  break;
930  }
931  }
932 
933  return modis_dimmap_nonll_field;
934 }
935 
936 void HDFCFUtil::add_scale_offset_attrs(AttrTable*at,const std::string& s_type, float svalue_f, double svalue_d, bool add_offset_found,
937  const std::string& o_type, float ovalue_f, double ovalue_d) {
938  at->del_attr("scale_factor");
939  string print_rep;
940  if(s_type!="Float64") {
941  print_rep = HDFCFUtil::print_attr(DFNT_FLOAT32,0,(void*)(&svalue_f));
942  at->append_attr("scale_factor", "Float32", print_rep);
943  }
944  else {
945  print_rep = HDFCFUtil::print_attr(DFNT_FLOAT64,0,(void*)(&svalue_d));
946  at->append_attr("scale_factor", "Float64", print_rep);
947  }
948 
949  if (true == add_offset_found) {
950  at->del_attr("add_offset");
951  if(o_type!="Float64") {
952  print_rep = HDFCFUtil::print_attr(DFNT_FLOAT32,0,(void*)(&ovalue_f));
953  at->append_attr("add_offset", "Float32", print_rep);
954  }
955  else {
956  print_rep = HDFCFUtil::print_attr(DFNT_FLOAT64,0,(void*)(&ovalue_d));
957  at->append_attr("add_offset", "Float64", print_rep);
958  }
959  }
960 }
961 
962 void HDFCFUtil::add_scale_str_offset_attrs(AttrTable*at,const std::string& s_type, const std::string& s_value_str, bool add_offset_found,
963  const std::string& o_type, float ovalue_f, double ovalue_d) {
964  at->del_attr("scale_factor");
965  string print_rep;
966  if(s_type!="Float64")
967  at->append_attr("scale_factor", "Float32", s_value_str);
968  else
969  at->append_attr("scale_factor", "Float64", s_value_str);
970 
971  if (true == add_offset_found) {
972  at->del_attr("add_offset");
973  if(o_type!="Float64") {
974  print_rep = HDFCFUtil::print_attr(DFNT_FLOAT32,0,(void*)(&ovalue_f));
975  at->append_attr("add_offset", "Float32", print_rep);
976  }
977  else {
978  print_rep = HDFCFUtil::print_attr(DFNT_FLOAT64,0,(void*)(&ovalue_d));
979  at->append_attr("add_offset", "Float64", print_rep);
980  }
981  }
982 }
984 void HDFCFUtil::handle_modis_special_attrs_disable_scale_comp(AttrTable *at,
985  const string &filename,
986  bool is_grid,
987  const string & newfname,
988  SOType sotype){
989 
990  // Declare scale_factor,add_offset, fillvalue and valid_range type in string format.
991  string scale_factor_type;
992  string add_offset_type;
993 
994  // Scale and offset values
995  string scale_factor_value="";
996  float orig_scale_value_float = 1;
997  double orig_scale_value_double = 1;
998  string add_offset_value="0";
999  float orig_offset_value_float = 0;
1000  double orig_offset_value_double = 0;
1001  bool add_offset_found = false;
1002 
1003 
1004  // Go through all attributes to find scale_factor,add_offset,_FillValue,valid_range
1005  // Number_Type information
1006  AttrTable::Attr_iter it = at->attr_begin();
1007  while (it!=at->attr_end())
1008  {
1009  if(at->get_name(it)=="scale_factor")
1010  {
1011  scale_factor_value = (*at->get_attr_vector(it)->begin());
1012  scale_factor_type = at->get_type(it);
1013  if(scale_factor_type =="Float64")
1014  orig_scale_value_double=atof(scale_factor_value.c_str());
1015  else
1016  orig_scale_value_float = (float)(atof(scale_factor_value.c_str()));
1017  }
1018 
1019  if(at->get_name(it)=="add_offset")
1020  {
1021  add_offset_value = (*at->get_attr_vector(it)->begin());
1022  add_offset_type = at->get_type(it);
1023 
1024  if(add_offset_type == "Float64")
1025  orig_offset_value_double = atof(add_offset_value.c_str());
1026  else
1027  orig_offset_value_float = (float)(atof(add_offset_value.c_str()));
1028  add_offset_found = true;
1029  }
1030 
1031  it++;
1032  }
1033 
1034  // According to our observations, it seems that MODIS products ALWAYS use the "scale" factor to
1035  // make bigger values smaller.
1036  // So for MODIS_MUL_SCALE products, if the scale of some variable is greater than 1,
1037  // it means that for this variable, the MODIS type for this variable may be MODIS_DIV_SCALE.
1038  // For the similar logic, we may need to change MODIS_DIV_SCALE to MODIS_MUL_SCALE and MODIS_EQ_SCALE
1039  // to MODIS_DIV_SCALE.
1040  // We indeed find such a case. HDF-EOS2 Grid MODIS_Grid_1km_2D of MOD(or MYD)09GA is a MODIS_EQ_SCALE.
1041  // However,
1042  // the scale_factor of the variable Range_1 in the MOD09GA product is 25. According to our observation,
1043  // this variable should be MODIS_DIV_SCALE.We verify this is true according to MODIS 09 product document
1044  // http://modis-sr.ltdri.org/products/MOD09_UserGuide_v1_3.pdf.
1045  // Since this conclusion is based on our observation, we would like to add a BESlog to detect if we find
1046  // the similar cases so that we can verify with the corresponding product documents to see if this is true.
1047  // More information,
1048  // We just verified with the data producer, the scale_factor for Range_1 and Range_c is 25 but the
1049  // equation is still multiplication instead of division. So we have to make this as a special case that
1050  // we don't want to change the scale and offset settings.
1051  // KY 2014-01-13
1052 
1053 
1054  if(scale_factor_value.length()!=0) {
1055  if (MODIS_EQ_SCALE == sotype || MODIS_MUL_SCALE == sotype) {
1056  if (orig_scale_value_float > 1 || orig_scale_value_double >1) {
1057 
1058  bool need_change_scale = true;
1059  if(true == is_grid) {
1060  if ((filename.size() >5) && ((filename.compare(0,5,"MOD09") == 0)|| (filename.compare(0,5,"MYD09")==0))) {
1061  if ((newfname.size() >5) && newfname.find("Range") != string::npos)
1062  need_change_scale = false;
1063  }
1064  else if((filename.size() >7)&&((filename.compare(0,7,"MOD16A2") == 0)|| (filename.compare(0,7,"MYD16A2")==0) ||
1065  (filename.compare(0,7,"MOD16A3")==0) || (filename.compare(0,7,"MYD16A3")==0)))
1066  need_change_scale = false;
1067  }
1068  if(true == need_change_scale) {
1069  sotype = MODIS_DIV_SCALE;
1070  (*BESLog::TheLog())<< "The field " << newfname << " scale factor is "<< scale_factor_value << endl
1071  << " But the original scale factor type is MODIS_MUL_SCALE or MODIS_EQ_SCALE. " << endl
1072  << " Now change it to MODIS_DIV_SCALE. "<<endl;
1073  }
1074  }
1075  }
1076 
1077  if (MODIS_DIV_SCALE == sotype) {
1078  if (orig_scale_value_float < 1 || orig_scale_value_double<1) {
1079  sotype = MODIS_MUL_SCALE;
1080  (*BESLog::TheLog())<< "The field " << newfname << " scale factor is "<< scale_factor_value << endl
1081  << " But the original scale factor type is MODIS_DIV_SCALE. " << endl
1082  << " Now change it to MODIS_MUL_SCALE. "<<endl;
1083  }
1084  }
1085 
1086 
1087  if ((MODIS_MUL_SCALE == sotype) &&(true == add_offset_found)) {
1088 
1089  float new_offset_value_float=0;
1090  double new_offset_value_double=0;
1091  if(add_offset_type!="Float64")
1092  new_offset_value_float = (orig_offset_value_float ==0)?0:(-1 * orig_offset_value_float *orig_scale_value_float);
1093  else
1094  new_offset_value_double = (orig_offset_value_double ==0)?0:(-1 * orig_offset_value_double *orig_scale_value_double);
1095 
1096  // May need to use another function to avoid the rounding error by atof.
1097  add_scale_str_offset_attrs(at,scale_factor_type,scale_factor_value,add_offset_found,
1098  add_offset_type,new_offset_value_float,new_offset_value_double);
1099 
1100  }
1101 
1102  if (MODIS_DIV_SCALE == sotype) {
1103 
1104  float new_scale_value_float=1;
1105  double new_scale_value_double=1;
1106  float new_offset_value_float=0;
1107  double new_offset_value_double=0;
1108 
1109 
1110  if(scale_factor_type !="Float64") {
1111  new_scale_value_float = 1.0/orig_scale_value_float;
1112  if (true == add_offset_found) {
1113  if(add_offset_type !="Float64")
1114  new_offset_value_float = (orig_offset_value_float==0)?0:(-1 * orig_offset_value_float *new_scale_value_float);
1115  else
1116  new_offset_value_double = (orig_offset_value_double==0)?0:(-1 * orig_offset_value_double *new_scale_value_float);
1117  }
1118  }
1119 
1120  else {
1121  new_scale_value_double = 1.0/orig_scale_value_double;
1122  if (true == add_offset_found) {
1123  if(add_offset_type !="Float64")
1124  new_offset_value_float = (orig_offset_value_float==0)?0:(-1 * orig_offset_value_float *new_scale_value_double);
1125  else
1126  new_offset_value_double = (orig_offset_value_double==0)?0:(-1 * orig_offset_value_double *new_scale_value_double);
1127  }
1128  }
1129 
1130  add_scale_offset_attrs(at,scale_factor_type,new_scale_value_float,new_scale_value_double,add_offset_found,
1131  add_offset_type,new_offset_value_float,new_offset_value_double);
1132 
1133  }
1134 
1135  }
1136 
1137 }
1138 
1139 // These routines will handle scale_factor,add_offset,valid_min,valid_max and other attributes
1140 // such as Number_Type to make sure the CF is followed.
1141 // For example, For the case that the scale and offset rule doesn't follow CF, the scale_factor and add_offset attributes are renamed
1142 // to orig_scale_factor and orig_add_offset to keep the original field info.
1143 // Note: This routine should only be called when MODIS Scale and offset rules don't follow CF.
1144 // Input parameters:
1145 // AttrTable at - DAS attribute table
1146 // string newfname - field name: this is used to identify special MODIS level 1B emissive and refSB fields
1147 // SOType sotype - MODIS scale and offset type
1148 // bool gridname_change_valid_range - Flag to check if this is the special MODIS VIP product
1149 // bool changedtype - Flag to check if the datatype of this field needs to be changed
1150 // bool change_fvtype - Flag to check if the attribute _FillValue needs to change to data type
1151 
1153 // Divide this function into smaller functions:
1154 //
1155 void HDFCFUtil::handle_modis_special_attrs(AttrTable *at, const string & filename,
1156  bool is_grid,const string & newfname,
1157  SOType sotype, bool gridname_change_valid_range,
1158  bool changedtype, bool & change_fvtype){
1159 
1160  // Declare scale_factor,add_offset, fillvalue and valid_range type in string format.
1161  string scale_factor_type;
1162  string add_offset_type;
1163  string fillvalue_type;
1164  string valid_range_type;
1165 
1166 
1167  // Scale and offset values
1168  string scale_factor_value="";
1169  float orig_scale_value = 1;
1170  string add_offset_value="0";
1171  float orig_offset_value = 0;
1172  bool add_offset_found = false;
1173 
1174  // Fillvalue
1175  string fillvalue="";
1176 
1177  // Valid range value
1178  string valid_range_value="";
1179 
1180  bool has_valid_range = false;
1181 
1182  // We may need to change valid_range to valid_min and valid_max. Initialize them here.
1183  float orig_valid_min = 0;
1184  float orig_valid_max = 0;
1185 
1186  // Number_Type also needs to be adjusted when datatype is changed
1187  string number_type_value="";
1188  string number_type_dap_type="";
1189 
1190  // Go through all attributes to find scale_factor,add_offset,_FillValue,valid_range
1191  // Number_Type information
1192  AttrTable::Attr_iter it = at->attr_begin();
1193  while (it!=at->attr_end())
1194  {
1195  if(at->get_name(it)=="scale_factor")
1196  {
1197  scale_factor_value = (*at->get_attr_vector(it)->begin());
1198  orig_scale_value = (float)(atof(scale_factor_value.c_str()));
1199  scale_factor_type = at->get_type(it);
1200  }
1201 
1202  if(at->get_name(it)=="add_offset")
1203  {
1204  add_offset_value = (*at->get_attr_vector(it)->begin());
1205  orig_offset_value = (float)(atof(add_offset_value.c_str()));
1206  add_offset_type = at->get_type(it);
1207  add_offset_found = true;
1208  }
1209 
1210  if(at->get_name(it)=="_FillValue")
1211  {
1212  fillvalue = (*at->get_attr_vector(it)->begin());
1213  fillvalue_type = at->get_type(it);
1214  }
1215 
1216  if(at->get_name(it)=="valid_range")
1217  {
1218  vector<string> *avalue = at->get_attr_vector(it);
1219  auto ait = avalue->begin();
1220  while(ait!=avalue->end())
1221  {
1222  valid_range_value += *ait;
1223  ait++;
1224  if(ait!=avalue->end())
1225  valid_range_value += ", ";
1226  }
1227  valid_range_type = at->get_type(it);
1228  if (false == gridname_change_valid_range) {
1229  orig_valid_min = (float)(atof((avalue->at(0)).c_str()));
1230  orig_valid_max = (float)(atof((avalue->at(1)).c_str()));
1231  }
1232  has_valid_range = true;
1233  }
1234 
1235  if(true == changedtype && (at->get_name(it)=="Number_Type"))
1236  {
1237  number_type_value = (*at->get_attr_vector(it)->begin());
1238  number_type_dap_type= at->get_type(it);
1239  }
1240 
1241  it++;
1242  }
1243 
1244  // Rename scale_factor and add_offset attribute names. Otherwise, they will be
1245  // misused by CF tools to generate wrong data values based on the CF scale and offset rule.
1246  if(scale_factor_value.length()!=0)
1247  {
1248  if(!(atof(scale_factor_value.c_str())==1 && atof(add_offset_value.c_str())==0)) //Rename them.
1249  {
1250  at->del_attr("scale_factor");
1251  at->append_attr("orig_scale_factor", scale_factor_type, scale_factor_value);
1252  if(add_offset_found)
1253  {
1254  at->del_attr("add_offset");
1255  at->append_attr("orig_add_offset", add_offset_type, add_offset_value);
1256  }
1257  }
1258  }
1259 
1260  // Change _FillValue datatype
1261  if(true == changedtype && fillvalue.length()!=0 && fillvalue_type!="Float32" && fillvalue_type!="Float64")
1262  {
1263  change_fvtype = true;
1264  at->del_attr("_FillValue");
1265  at->append_attr("_FillValue", "Float32", fillvalue);
1266  }
1267 
1268  float valid_max = 0;
1269  float valid_min = 0;
1270 
1271  it = at->attr_begin();
1272  bool handle_modis_l1b = false;
1273 
1274  // MODIS level 1B's Emissive and RefSB fields' scale_factor and add_offset
1275  // get changed according to different vertical levels.
1276  // So we need to handle them specifically.
1277  if (sotype == MODIS_MUL_SCALE && true ==changedtype) {
1278 
1279  // Emissive should be at the end of the field name such as "..._Emissive"
1280  string emissive_str = "Emissive";
1281  string RefSB_str = "RefSB";
1282  bool is_emissive_field = false;
1283  bool is_refsb_field = false;
1284  if(newfname.find(emissive_str)!=string::npos) {
1285  if(0 == newfname.compare(newfname.size()-emissive_str.size(),emissive_str.size(),emissive_str))
1286  is_emissive_field = true;
1287  }
1288 
1289  if(newfname.find(RefSB_str)!=string::npos) {
1290  if(0 == newfname.compare(newfname.size()-RefSB_str.size(),RefSB_str.size(),RefSB_str))
1291  is_refsb_field = true;
1292  }
1293 
1294  // We will calculate the maximum and minimum scales and offsets within all the vertical levels.
1295  if ((true == is_emissive_field) || (true== is_refsb_field)){
1296 
1297  float scale_max = 0;
1298  float scale_min = 100000;
1299 
1300  float offset_max = 0;
1301  float offset_min = 0;
1302 
1303  float temp_var_val = 0;
1304 
1305  string orig_long_name_value;
1306  string modify_long_name_value;
1307  string str_removed_from_long_name=" Scaled Integers";
1308  string radiance_units_value;
1309 
1310  while (it!=at->attr_end())
1311  {
1312  // Radiance field(Emissive field)
1313  if(true == is_emissive_field) {
1314 
1315  if ("radiance_scales" == (at->get_name(it))) {
1316  vector<string> *avalue = at->get_attr_vector(it);
1317  for (vector<string>::const_iterator ait = avalue->begin();ait !=avalue->end();++ait) {
1318  temp_var_val = (float)(atof((*ait).c_str()));
1319 #if 0
1320  // Check if this works in the future.
1321  for (const auto &avalue_ele:*avalue) {
1322  temp_var_val = (float)(atof((avalue_ele).c_str()));
1323 #endif
1324  if (temp_var_val > scale_max)
1325  scale_max = temp_var_val;
1326  if (temp_var_val < scale_min)
1327  scale_min = temp_var_val;
1328  }
1329  }
1330 
1331  if ("radiance_offsets" == (at->get_name(it))) {
1332  vector<string> *avalue = at->get_attr_vector(it);
1333  for (vector<string>::const_iterator ait = avalue->begin();ait !=avalue->end();++ait) {
1334  temp_var_val = (float)(atof((*ait).c_str()));
1335  if (temp_var_val > offset_max)
1336  offset_max = temp_var_val;
1337  if (temp_var_val < scale_min)
1338  offset_min = temp_var_val;
1339  }
1340  }
1341 
1342  if ("long_name" == (at->get_name(it))) {
1343  orig_long_name_value = (*at->get_attr_vector(it)->begin());
1344  if (orig_long_name_value.find(str_removed_from_long_name)!=string::npos) {
1345  if(0 == orig_long_name_value.compare(orig_long_name_value.size()-str_removed_from_long_name.size(),
1346  str_removed_from_long_name.size(),str_removed_from_long_name)) {
1347 
1348  modify_long_name_value =
1349  orig_long_name_value.substr(0,orig_long_name_value.size()-str_removed_from_long_name.size());
1350  at->del_attr("long_name");
1351  at->append_attr("long_name","String",modify_long_name_value);
1352  at->append_attr("orig_long_name","String",orig_long_name_value);
1353  }
1354  }
1355  }
1356 
1357  if ("radiance_units" == (at->get_name(it)))
1358  radiance_units_value = (*at->get_attr_vector(it)->begin());
1359  }
1360 
1361  if (true == is_refsb_field) {
1362  if ("reflectance_scales" == (at->get_name(it))) {
1363  vector<string> *avalue = at->get_attr_vector(it);
1364  for (vector<string>::const_iterator ait = avalue->begin();ait !=avalue->end();++ait) {
1365  temp_var_val = (float)(atof((*ait).c_str()));
1366  if (temp_var_val > scale_max)
1367  scale_max = temp_var_val;
1368  if (temp_var_val < scale_min)
1369  scale_min = temp_var_val;
1370  }
1371  }
1372 
1373  if ("reflectance_offsets" == (at->get_name(it))) {
1374 
1375  vector<string> *avalue = at->get_attr_vector(it);
1376  for (vector<string>::const_iterator ait = avalue->begin();ait !=avalue->end();++ait) {
1377  temp_var_val = (float)(atof((*ait).c_str()));
1378  if (temp_var_val > offset_max)
1379  offset_max = temp_var_val;
1380  if (temp_var_val < scale_min)
1381  offset_min = temp_var_val;
1382  }
1383  }
1384 
1385  if ("long_name" == (at->get_name(it))) {
1386  orig_long_name_value = (*at->get_attr_vector(it)->begin());
1387  if (orig_long_name_value.find(str_removed_from_long_name)!=string::npos) {
1388  if(0 == orig_long_name_value.compare(orig_long_name_value.size()-str_removed_from_long_name.size(),
1389  str_removed_from_long_name.size(),str_removed_from_long_name)) {
1390 
1391  modify_long_name_value =
1392  orig_long_name_value.substr(0,orig_long_name_value.size()-str_removed_from_long_name.size());
1393  at->del_attr("long_name");
1394  at->append_attr("long_name","String",modify_long_name_value);
1395  at->append_attr("orig_long_name","String",orig_long_name_value);
1396  }
1397  }
1398  }
1399  }
1400  it++;
1401  }
1402 
1403  // Calculate the final valid_max and valid_min.
1404  if (scale_min <= 0)
1405  throw InternalErr(__FILE__,__LINE__,"the scale factor should always be greater than 0.");
1406 
1407  if (orig_valid_max > offset_min)
1408  valid_max = (orig_valid_max-offset_min)*scale_max;
1409  else
1410  valid_max = (orig_valid_max-offset_min)*scale_min;
1411 
1412  if (orig_valid_min > offset_max)
1413  valid_min = (orig_valid_min-offset_max)*scale_min;
1414  else
1415  valid_min = (orig_valid_min -offset_max)*scale_max;
1416 
1417  // These physical variables should be greater than 0
1418  if (valid_min < 0)
1419  valid_min = 0;
1420 
1421  string print_rep = HDFCFUtil::print_attr(DFNT_FLOAT32,0,(void*)(&valid_min));
1422  at->append_attr("valid_min","Float32",print_rep);
1423  print_rep = HDFCFUtil::print_attr(DFNT_FLOAT32,0,(void*)(&valid_max));
1424  at->append_attr("valid_max","Float32",print_rep);
1425  at->del_attr("valid_range");
1426  handle_modis_l1b = true;
1427 
1428  // Change the units for the emissive field
1429  if (true == is_emissive_field && radiance_units_value.size() >0) {
1430  at->del_attr("units");
1431  at->append_attr("units","String",radiance_units_value);
1432  }
1433  }
1434  }
1435 
1436  // Handle other valid_range attributes
1437  if(true == changedtype && true == has_valid_range && false == handle_modis_l1b) {
1438 
1439  // If the gridname_change_valid_range is true, call a special function to handle this.
1440  if (true == gridname_change_valid_range)
1441  HDFCFUtil::handle_modis_vip_special_attrs(valid_range_value,scale_factor_value,valid_min,valid_max);
1442  else if(scale_factor_value.length()!=0) {
1443 
1444  // We found MODIS products always scale to a smaller value. If somehow the original scale factor
1445  // is smaller than 1, the scale/offset should be the multiplication rule.(new_data =scale*(old_data-offset))
1446  // If the original scale factor is greater than 1, the scale/offset rule should be division rule.
1447  // New_data = (old_data-offset)/scale.
1448  // We indeed find such a case. HDF-EOS2 Grid MODIS_Grid_1km_2D of MOD(or MYD)09GA is a MODIS_EQ_SCALE.
1449  // However,
1450  // the scale_factor of the variable Range_1 in the MOD09GA product is 25. According to our observation,
1451  // this variable should be MODIS_DIV_SCALE.We verify this is true according to MODIS 09 product document
1452  // http://modis-sr.ltdri.org/products/MOD09_UserGuide_v1_3.pdf.
1453  // Since this conclusion is based on our observation, we would like to add a BESlog to detect if we find
1454  // the similar cases so that we can verify with the corresponding product documents to see if this is true.
1455  // More information,
1456  // We just verified with the data producer, the scale_factor for Range_1 and Range_c is 25 but the
1457  // equation is still multiplication instead of division. So we have to make this as a special case that
1458  // we don't want to change the scale and offset settings.
1459  // KY 2014-01-13
1460 
1461  if (MODIS_EQ_SCALE == sotype || MODIS_MUL_SCALE == sotype) {
1462  if (orig_scale_value > 1) {
1463 
1464  bool need_change_scale = true;
1465  if(true == is_grid) {
1466  if ((filename.size() >5) && ((filename.compare(0,5,"MOD09") == 0)|| (filename.compare(0,5,"MYD09")==0))) {
1467  if ((newfname.size() >5) && newfname.find("Range") != string::npos)
1468  need_change_scale = false;
1469  }
1470 
1471  else if((filename.size() >7)&&((filename.compare(0,7,"MOD16A2") == 0)|| (filename.compare(0,7,"MYD16A2")==0) ||
1472  (filename.compare(0,7,"MOD16A3")==0) || (filename.compare(0,7,"MYD16A3")==0)))
1473  need_change_scale = false;
1474  }
1475  if(true == need_change_scale) {
1476  sotype = MODIS_DIV_SCALE;
1477  (*BESLog::TheLog())<< "The field " << newfname << " scale factor is "<< orig_scale_value << endl
1478  << " But the original scale factor type is MODIS_MUL_SCALE or MODIS_EQ_SCALE. " << endl
1479  << " Now change it to MODIS_DIV_SCALE. "<<endl;
1480  }
1481  }
1482  }
1483 
1484  if (MODIS_DIV_SCALE == sotype) {
1485  if (orig_scale_value < 1) {
1486  sotype = MODIS_MUL_SCALE;
1487  (*BESLog::TheLog())<< "The field " << newfname << " scale factor is "<< orig_scale_value << endl
1488  << " But the original scale factor type is MODIS_DIV_SCALE. " << endl
1489  << " Now change it to MODIS_MUL_SCALE. "<<endl;
1490  }
1491  }
1492 
1493  if(sotype == MODIS_MUL_SCALE) {
1494  valid_min = (orig_valid_min - orig_offset_value)*orig_scale_value;
1495  valid_max = (orig_valid_max - orig_offset_value)*orig_scale_value;
1496  }
1497  else if (sotype == MODIS_DIV_SCALE) {
1498  valid_min = (orig_valid_min - orig_offset_value)/orig_scale_value;
1499  valid_max = (orig_valid_max - orig_offset_value)/orig_scale_value;
1500  }
1501  else if (sotype == MODIS_EQ_SCALE) {
1502  valid_min = orig_valid_min * orig_scale_value + orig_offset_value;
1503  valid_max = orig_valid_max * orig_scale_value + orig_offset_value;
1504  }
1505  }
1506  else {// This is our current assumption.
1507  if (MODIS_MUL_SCALE == sotype || MODIS_DIV_SCALE == sotype) {
1508  valid_min = orig_valid_min - orig_offset_value;
1509  valid_max = orig_valid_max - orig_offset_value;
1510  }
1511  else if (MODIS_EQ_SCALE == sotype) {
1512  valid_min = orig_valid_min + orig_offset_value;
1513  valid_max = orig_valid_max + orig_offset_value;
1514  }
1515  }
1516 
1517  // Append the corrected valid_min and valid_max. (valid_min,valid_max) is equivalent to valid_range.
1518  string print_rep = HDFCFUtil::print_attr(DFNT_FLOAT32,0,(void*)(&valid_min));
1519  at->append_attr("valid_min","Float32",print_rep);
1520  print_rep = HDFCFUtil::print_attr(DFNT_FLOAT32,0,(void*)(&valid_max));
1521  at->append_attr("valid_max","Float32",print_rep);
1522  at->del_attr("valid_range");
1523  }
1524 
1525  // We also find that there is an attribute called "Number_Type" that will stores the original attribute
1526  // datatype. If the datatype gets changed, the attribute is confusion. here we can change the attribute
1527  // name to "Number_Type_Orig"
1528  if(true == changedtype && number_type_dap_type !="" ) {
1529  at->del_attr("Number_Type");
1530  at->append_attr("Number_Type_Orig",number_type_dap_type,number_type_value);
1531  }
1532 }
1533 
1534  // This routine makes the MeaSUREs VIP attributes follow CF.
1535 // valid_range attribute uses char type instead of the raw data's datatype.
1536 void HDFCFUtil::handle_modis_vip_special_attrs(const std::string& valid_range_value,
1537  const std::string& scale_factor_value,
1538  float& valid_min, float &valid_max) {
1539 
1540  int16 vip_orig_valid_min = 0;
1541  int16 vip_orig_valid_max = 0;
1542 
1543  size_t found = valid_range_value.find_first_of(",");
1544  size_t found_from_end = valid_range_value.find_last_of(",");
1545  if (string::npos == found)
1546  throw InternalErr(__FILE__,__LINE__,"should find the separator ,");
1547  if (found != found_from_end)
1548  throw InternalErr(__FILE__,__LINE__,"There should be only one separator.");
1549 
1550 #if 0
1551  //istringstream(valid_range_value.substr(0,found))>>orig_valid_min;
1552  //istringstream(valid_range_value.substr(found+1))>>orig_valid_max;
1553 #endif
1554 
1555  vip_orig_valid_min = (short) (atoi((valid_range_value.substr(0,found)).c_str()));
1556  vip_orig_valid_max = (short) (atoi((valid_range_value.substr(found+1)).c_str()));
1557 
1558  int16 scale_factor_number = 1;
1559 
1560  scale_factor_number = (short)(atoi(scale_factor_value.c_str()));
1561 
1562  if(scale_factor_number !=0) {
1563  valid_min = (float)(vip_orig_valid_min/scale_factor_number);
1564  valid_max = (float)(vip_orig_valid_max/scale_factor_number);
1565  }
1566  else
1567  throw InternalErr(__FILE__,__LINE__,"The scale_factor_number should not be zero.");
1568 }
1569 
1570 // Make AMSR-E attributes follow CF.
1571 // Change SCALE_FACTOR and OFFSET to CF names: scale_factor and add_offset.
1572 void HDFCFUtil::handle_amsr_attrs(AttrTable *at) {
1573 
1574  AttrTable::Attr_iter it = at->attr_begin();
1575 
1576  string scale_factor_value="";
1577  string add_offset_value="0";
1578 
1579  string scale_factor_type;
1580  string add_offset_type;
1581 
1582  bool OFFSET_found = false;
1583  bool Scale_found = false;
1584  bool SCALE_FACTOR_found = false;
1585 
1586  while (it!=at->attr_end())
1587  {
1588  if(at->get_name(it)=="SCALE_FACTOR")
1589  {
1590  scale_factor_value = (*at->get_attr_vector(it)->begin());
1591  scale_factor_type = at->get_type(it);
1592  SCALE_FACTOR_found = true;
1593  }
1594 
1595  if(at->get_name(it)=="Scale")
1596  {
1597  scale_factor_value = (*at->get_attr_vector(it)->begin());
1598  scale_factor_type = at->get_type(it);
1599  Scale_found = true;
1600  }
1601 
1602  if(at->get_name(it)=="OFFSET")
1603  {
1604  add_offset_value = (*at->get_attr_vector(it)->begin());
1605  add_offset_type = at->get_type(it);
1606  OFFSET_found = true;
1607  }
1608  it++;
1609  }
1610 
1611  if (true == SCALE_FACTOR_found) {
1612  at->del_attr("SCALE_FACTOR");
1613  at->append_attr("scale_factor",scale_factor_type,scale_factor_value);
1614  }
1615 
1616  if (true == Scale_found) {
1617  at->del_attr("Scale");
1618  at->append_attr("scale_factor",scale_factor_type,scale_factor_value);
1619  }
1620 
1621  if (true == OFFSET_found) {
1622  at->del_attr("OFFSET");
1623  at->append_attr("add_offset",add_offset_type,add_offset_value);
1624  }
1625 
1626 }
1627 
1628 //This function obtains the latitude and longitude dimension info. of an
1629 //HDF-EOS2 grid after the handler translates the HDF-EOS to CF.
1630 // Dimension info. includes dimension name and dimension size.
1631 void HDFCFUtil::obtain_grid_latlon_dim_info(const HDFEOS2::GridDataset* gdset,
1632  string & dim0name,
1633  int32 & dim0size,
1634  string & dim1name,
1635  int32& dim1size){
1636 
1637  const vector<HDFEOS2::Field*>gfields = gdset->getDataFields();
1638  for (const auto &gf:gfields) {
1639 
1640  // Check the dimensions for Latitude
1641  if(1 == gf->getFieldType()) {
1642  const vector<HDFEOS2::Dimension*>& dims= gf->getCorrectedDimensions();
1643 
1644  //2-D latitude
1645  if(2 == dims.size()) {
1646  // Most time, it is YDim Major. We will see if we have a real case to
1647  // check if the handling is right for the XDim Major case.
1648  if(true == gf->getYDimMajor()) {
1649  dim0name = dims[0]->getName();
1650  dim0size = dims[0]->getSize();
1651  dim1name = dims[1]->getName();
1652  dim1size = dims[1]->getSize();
1653  }
1654  else {
1655  dim0name = dims[1]->getName();
1656  dim0size = dims[1]->getSize();
1657  dim1name = dims[0]->getName();
1658  dim1size = dims[0]->getSize();
1659  }
1660  break;
1661  }
1662 
1663  //1-D latitude
1664  if( 1 == dims.size()) {
1665  dim0name = dims[0]->getName();
1666  dim0size = dims[0]->getSize();
1667  }
1668  }
1669 
1670  // Longitude, if longitude is checked first, it goes here.
1671  if(2 == gf->getFieldType()) {
1672  const vector<HDFEOS2::Dimension*>& dims= gf->getCorrectedDimensions();
1673  if(2 == dims.size()) {
1674 
1675  // Most time, it is YDim Major. We will see if we have a real case to
1676  // check if the handling is right for the XDim Major case.
1677  if(true == gf->getYDimMajor()) {
1678  dim0name = dims[0]->getName();
1679  dim0size = dims[0]->getSize();
1680  dim1name = dims[1]->getName();
1681  dim1size = dims[1]->getSize();
1682  }
1683  else {
1684  dim0name = dims[1]->getName();
1685  dim0size = dims[1]->getSize();
1686  dim1name = dims[0]->getName();
1687  dim1size = dims[0]->getSize();
1688  }
1689  break;
1690  }
1691  if( 1 == dims.size()) {
1692  dim1name = dims[0]->getName();
1693  dim1size = dims[0]->getSize();
1694 
1695  }
1696  }
1697  }
1698  if(""==dim0name || ""==dim1name || dim0size<0 || dim1size <0)
1699  throw InternalErr (__FILE__, __LINE__,"Fail to obtain grid lat/lon dimension info.");
1700 
1701 }
1702 
1703 // This function adds the 1-D cf grid projection mapping attribute to data variables
1704 // it is called by the function add_cf_grid_attrs.
1705 void HDFCFUtil::add_cf_grid_mapping_attr(DAS &das, const HDFEOS2::GridDataset*gdset,const string& cf_projection,
1706  const string & dim0name,int32 dim0size,const string &dim1name,int32 dim1size) {
1707 
1708  // Check >=2-D fields, check if they hold the dim0name,dim0size etc., yes, add the attribute cf_projection.
1709  const vector<HDFEOS2::Field*>gfields = gdset->getDataFields();
1710 
1711  for (const auto &gf:gfields) {
1712  if(0 == gf->getFieldType() && gf->getRank() >1) {
1713  bool has_dim0 = false;
1714  bool has_dim1 = false;
1715  const vector<HDFEOS2::Dimension*>& dims= gf->getCorrectedDimensions();
1716  for (const auto &dim:dims) {
1717  if(dim->getName()== dim0name && dim->getSize() == dim0size)
1718  has_dim0 = true;
1719  else if(dim->getName()== dim1name && dim->getSize() == dim1size)
1720  has_dim1 = true;
1721 
1722  }
1723  if(true == has_dim0 && true == has_dim1) {// Need to add the grid_mapping attribute
1724  AttrTable *at = das.get_table(gf->getNewName());
1725  if (!at)
1726  at = das.add_table(gf->getNewName(), new AttrTable);
1727 
1728  // The dummy projection name is the value of the grid_mapping attribute
1729  at->append_attr("grid_mapping","String",cf_projection);
1730  }
1731  }
1732  }
1733 }
1734 
1735 //This function adds 1D grid mapping CF attributes to CV and data variables.
1736 void HDFCFUtil::add_cf_grid_cv_attrs(DAS & das, const HDFEOS2::GridDataset *gdset) {
1737 
1738  //1. Check the projection information, now, we only handle sinusoidal now
1739  if(GCTP_SNSOID == gdset->getProjection().getCode()) {
1740 
1741  //2. Obtain the dimension information from latitude and longitude(fieldtype =1 or fieldtype =2)
1742  string dim0name;
1743  string dim1name;
1744  int32 dim0size = -1;
1745  int32 dim1size = -1;
1746 
1747  HDFCFUtil::obtain_grid_latlon_dim_info(gdset,dim0name,dim0size,dim1name,dim1size);
1748 
1749  //3. Add 1D CF attributes to the 1-D CV variables and the dummy projection variable
1750  AttrTable *at = das.get_table(dim0name);
1751  if (!at)
1752  at = das.add_table(dim0name, new AttrTable);
1753  at->append_attr("standard_name","String","projection_y_coordinate");
1754 
1755  string long_name="y coordinate of projection for grid "+ gdset->getName();
1756  at->append_attr("long_name","String",long_name);
1757  // Change this to meter.
1758  at->append_attr("units","string","meter");
1759 
1760  at->append_attr("_CoordinateAxisType","string","GeoY");
1761 
1762  at = das.get_table(dim1name);
1763  if (!at)
1764  at = das.add_table(dim1name, new AttrTable);
1765 
1766  at->append_attr("standard_name","String","projection_x_coordinate");
1767  long_name="x coordinate of projection for grid "+ gdset->getName();
1768  at->append_attr("long_name","String",long_name);
1769 
1770  // change this to meter.
1771  at->append_attr("units","string","meter");
1772  at->append_attr("_CoordinateAxisType","string","GeoX");
1773 
1774  // Add the attributes for the dummy projection variable.
1775  string cf_projection_base = "eos_cf_projection";
1776  string cf_projection = HDFCFUtil::get_CF_string(gdset->getName()) +"_"+cf_projection_base;
1777  at = das.get_table(cf_projection);
1778  if (!at)
1779  at = das.add_table(cf_projection, new AttrTable);
1780 
1781  at->append_attr("grid_mapping_name","String","sinusoidal");
1782  at->append_attr("longitude_of_central_meridian","Float64","0.0");
1783  at->append_attr("earth_radius","Float64","6371007.181");
1784 
1785  at->append_attr("_CoordinateAxisTypes","string","GeoX GeoY");
1786 
1787  // Fill in the data fields that contains the dim0name and dim1name dimensions with the grid_mapping
1788  // We only apply to >=2D data fields.
1789  HDFCFUtil::add_cf_grid_mapping_attr(das,gdset,cf_projection,dim0name,dim0size,dim1name,dim1size);
1790  }
1791 
1792 }
1793 
1794 // This function adds the 1-D horizontal coordinate variables as well as the dummy projection variable to the grid.
1795 //Note: Since we don't add these artifical CF variables to our main engineering at HDFEOS2.cc, the information
1796 // to handle DAS won't pass to DDS by the file pointer, we need to re-call the routines to check projection
1797 // and dimension. The time to retrieve these information is trivial compared with the whole translation.
1798 void HDFCFUtil::add_cf_grid_cvs(DDS & dds, const HDFEOS2::GridDataset *gdset) {
1799 
1800  //1. Check the projection information, now, we only handle sinusoidal now
1801  if(GCTP_SNSOID == gdset->getProjection().getCode()) {
1802 
1803  //2. Obtain the dimension information from latitude and longitude(fieldtype =1 or fieldtype =2)
1804  string dim0name;
1805  string dim1name;
1806  int32 dim0size = -1;
1807  int32 dim1size = -1;
1808  HDFCFUtil::obtain_grid_latlon_dim_info(gdset,dim0name,dim0size,dim1name,dim1size);
1809 
1810  //3. Add the 1-D CV variables and the dummy projection variable
1811  // Note: we just need to pass the parameters that calculate 1-D cv to the data reading function,
1812  // in that way, we save the open cost of HDF-EOS2.
1813  BaseType *bt_dim0 = nullptr;
1814  BaseType *bt_dim1 = nullptr;
1815 
1816  HDFEOS2GeoCF1D * ar_dim0 = nullptr;
1817  HDFEOS2GeoCF1D * ar_dim1 = nullptr;
1818 
1819  const float64 *upleft = nullptr;
1820  const float64 *lowright = nullptr;
1821 
1822  try {
1823 
1824  bt_dim0 = new(HDFFloat64)(dim0name,gdset->getName());
1825  bt_dim1 = new(HDFFloat64)(dim1name,gdset->getName());
1826 
1827  // Obtain the upleft and lowright coordinates
1828  upleft = gdset->getInfo().getUpLeft();
1829  lowright = gdset->getInfo().getLowRight();
1830 
1831  // Note ar_dim0 is y, ar_dim1 is x.
1832  ar_dim0 = new HDFEOS2GeoCF1D(GCTP_SNSOID,
1833  upleft[1],lowright[1],dim0size,dim0name,bt_dim0);
1834  ar_dim0->append_dim(dim0size,dim0name);
1835 
1836  ar_dim1 = new HDFEOS2GeoCF1D(GCTP_SNSOID,
1837  upleft[0],lowright[0],dim1size,dim1name,bt_dim1);
1838  ar_dim1->append_dim(dim1size,dim1name);
1839  dds.add_var(ar_dim0);
1840  dds.add_var(ar_dim1);
1841 
1842  }
1843  catch(...) {
1844  if(bt_dim0)
1845  delete bt_dim0;
1846  if(bt_dim1)
1847  delete bt_dim1;
1848  if(ar_dim0)
1849  delete ar_dim0;
1850  if(ar_dim1)
1851  delete ar_dim1;
1852  throw InternalErr(__FILE__,__LINE__,"Unable to allocate the HDFEOS2GeoCF1D instance.");
1853  }
1854 
1855  if(bt_dim0)
1856  delete bt_dim0;
1857  if(bt_dim1)
1858  delete bt_dim1;
1859  if(ar_dim0)
1860  delete ar_dim0;
1861  if(ar_dim1)
1862  delete ar_dim1;
1863 
1864  // Also need to add the dummy projection variable.
1865  string cf_projection_base = "eos_cf_projection";
1866 
1867  // To handle multi-grid cases, we need to add the grid name.
1868  string cf_projection = HDFCFUtil::get_CF_string(gdset->getName()) +"_"+cf_projection_base;
1869 
1870  HDFEOS2GeoCFProj * dummy_proj_cf = new HDFEOS2GeoCFProj(cf_projection,gdset->getName());
1871  dds.add_var(dummy_proj_cf);
1872  if(dummy_proj_cf)
1873  delete dummy_proj_cf;
1874 
1875  }
1876 
1877 }
1878 #endif
1879 
1880  // Check OBPG attributes. Specifically, check if slope and intercept can be obtained from the file level.
1881  // If having global slope and intercept, obtain OBPG scaling, slope and intercept values.
1882 void HDFCFUtil::check_obpg_global_attrs(HDFSP::File *f, std::string & scaling,
1883  float & slope,bool &global_slope_flag,
1884  float & intercept, bool & global_intercept_flag){
1885 
1886  HDFSP::SD* spsd = f->getSD();
1887 
1888  for (const auto &attr:spsd->getAttributes()) {
1889 
1890  //We want to add two new attributes, "scale_factor" and "add_offset" to data fields if the scaling equation is linear.
1891  // OBPG products use "Slope" instead of "scale_factor", "intercept" instead of "add_offset". "Scaling" describes if the equation is linear.
1892  // Their values will be copied directly from File attributes. This addition is for OBPG L3 only.
1893  // We also add OBPG L2 support since all OBPG level 2 products(CZCS,MODISA,MODIST,OCTS,SeaWiFS) we evaluate use Slope,intercept and linear equation
1894  // for the final data. KY 2012-09-06
1895  if(f->getSPType()==OBPGL3 || f->getSPType() == OBPGL2)
1896  {
1897  if(attr->getName()=="Scaling")
1898  {
1899  string tmpstring(attr->getValue().begin(), attr->getValue().end());
1900  scaling = tmpstring;
1901  }
1902  if(attr->getName()=="Slope" || attr->getName()=="slope")
1903  {
1904  global_slope_flag = true;
1905 
1906  switch(attr->getType())
1907  {
1908 #define GET_SLOPE(TYPE, CAST) \
1909  case DFNT_##TYPE: \
1910  { \
1911  CAST tmpvalue = *(CAST*)&(attr->getValue()[0]); \
1912  slope = (float)tmpvalue; \
1913  } \
1914  break;
1915  GET_SLOPE(INT16, int16)
1916  GET_SLOPE(INT32, int32)
1917  GET_SLOPE(FLOAT32, float)
1918  GET_SLOPE(FLOAT64, double)
1919  default:
1920  throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
1921 #undef GET_SLOPE
1922  }
1923 #if 0
1924  //slope = *(float*)&((*i)->getValue()[0]);
1925 #endif
1926  }
1927  if(attr->getName()=="Intercept" || attr->getName()=="intercept")
1928  {
1929  global_intercept_flag = true;
1930  switch(attr->getType())
1931  {
1932 #define GET_INTERCEPT(TYPE, CAST) \
1933  case DFNT_##TYPE: \
1934  { \
1935  CAST tmpvalue = *(CAST*)&(attr->getValue()[0]); \
1936  intercept = (float)tmpvalue; \
1937  } \
1938  break;
1939  GET_INTERCEPT(INT16, int16)
1940  GET_INTERCEPT(INT32, int32)
1941  GET_INTERCEPT(FLOAT32, float)
1942  GET_INTERCEPT(FLOAT64, double)
1943  default:
1944  throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
1945 #undef GET_INTERCEPT
1946  }
1947 #if 0
1948  //intercept = *(float*)&((*i)->getValue()[0]);
1949 #endif
1950  }
1951  }
1952  }
1953 }
1954 
1955 // For some OBPG files that only provide slope and intercept at the file level,
1956 // global slope and intercept are needed to add to all fields and their names are needed to be changed to scale_factor and add_offset.
1957 // For OBPG files that provide slope and intercept at the field level, slope and intercept are needed to rename to scale_factor and add_offset.
1958 void HDFCFUtil::add_obpg_special_attrs(const HDFSP::File*f,DAS &das,
1959  const HDFSP::SDField *onespsds,
1960  const string& scaling, float& slope,
1961  const bool& global_slope_flag,
1962  float& intercept,
1963  const bool & global_intercept_flag) {
1964 
1965  AttrTable *at = das.get_table(onespsds->getNewName());
1966  if (!at)
1967  at = das.add_table(onespsds->getNewName(), new AttrTable);
1968 
1969  //For OBPG L2 and L3 only.Some OBPG products put "slope" and "Intercept" etc. in the field attributes
1970  // We still need to handle the scale and offset here.
1971  bool scale_factor_flag = false;
1972  bool add_offset_flag = false;
1973  bool slope_flag = false;
1974  bool intercept_flag = false;
1975 
1976  if(f->getSPType()==OBPGL3 || f->getSPType() == OBPGL2) {// Begin OPBG CF attribute handling(Checking "slope" and "intercept")
1977 #if 0
1978  for(vector<HDFSP::Attribute *>::const_iterator i=onespsds->getAttributes().begin();
1979  i!=onespsds->getAttributes().end();i++) {
1980 #endif
1981  for (const auto &attr:onespsds->getAttributes()) {
1982 
1983  if(global_slope_flag != true && (attr->getName()=="Slope" || attr->getName()=="slope"))
1984  {
1985  slope_flag = true;
1986 
1987  switch(attr->getType())
1988  {
1989 #define GET_SLOPE(TYPE, CAST) \
1990  case DFNT_##TYPE: \
1991  { \
1992  CAST tmpvalue = *(CAST*)&(attr->getValue()[0]); \
1993  slope = (float)tmpvalue; \
1994  } \
1995  break;
1996 
1997  GET_SLOPE(INT16, int16)
1998  GET_SLOPE(INT32, int32)
1999  GET_SLOPE(FLOAT32, float)
2000  GET_SLOPE(FLOAT64, double)
2001  default:
2002  throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
2003 
2004 #undef GET_SLOPE
2005  }
2006 #if 0
2007  //slope = *(float*)&((*i)->getValue()[0]);
2008 #endif
2009  }
2010  if(global_intercept_flag != true && (attr->getName()=="Intercept" || attr->getName()=="intercept"))
2011  {
2012  intercept_flag = true;
2013  switch(attr->getType())
2014  {
2015 #define GET_INTERCEPT(TYPE, CAST) \
2016  case DFNT_##TYPE: \
2017  { \
2018  CAST tmpvalue = *(CAST*)&(attr->getValue()[0]); \
2019  intercept = (float)tmpvalue; \
2020  } \
2021  break;
2022  GET_INTERCEPT(INT16, int16)
2023  GET_INTERCEPT(INT32, int32)
2024  GET_INTERCEPT(FLOAT32, float)
2025  GET_INTERCEPT(FLOAT64, double)
2026  default:
2027  throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
2028 
2029 #undef GET_INTERCEPT
2030  }
2031 #if 0
2032  //intercept = *(float*)&((*i)->getValue()[0]);
2033 #endif
2034  }
2035  }
2036  } // End of checking "slope" and "intercept"
2037 
2038  // Checking if OBPG has "scale_factor" ,"add_offset", generally checking for "long_name" attributes.
2039  for (const auto& attr:onespsds->getAttributes()) {
2040 
2041  if((f->getSPType()==OBPGL3 || f->getSPType() == OBPGL2) && attr->getName()=="scale_factor")
2042  scale_factor_flag = true;
2043 
2044  if((f->getSPType()==OBPGL3 || f->getSPType() == OBPGL2) && attr->getName()=="add_offset")
2045  add_offset_flag = true;
2046  }
2047 
2048  // Checking if the scale and offset equation is linear, this is only for OBPG level 3.
2049  // Also checking if having the fill value and adding fill value if not having one for some OBPG data type
2050  if((f->getSPType() == OBPGL2 || f->getSPType()==OBPGL3 )&& onespsds->getFieldType()==0)
2051  {
2052 
2053  if((OBPGL3 == f->getSPType() && (scaling.find("linear")!=string::npos)) || OBPGL2 == f->getSPType() ) {
2054 
2055  if(false == scale_factor_flag && (true == slope_flag || true == global_slope_flag))
2056  {
2057  string print_rep = HDFCFUtil::print_attr(DFNT_FLOAT32, 0, (void*)&slope);
2058  at->append_attr("scale_factor", HDFCFUtil::print_type(DFNT_FLOAT32), print_rep);
2059  }
2060 
2061  if(false == add_offset_flag && (true == intercept_flag || true == global_intercept_flag))
2062  {
2063  string print_rep = HDFCFUtil::print_attr(DFNT_FLOAT32, 0, (void*)&intercept);
2064  at->append_attr("add_offset", HDFCFUtil::print_type(DFNT_FLOAT32), print_rep);
2065  }
2066  }
2067 
2068  bool has_fill_value = false;
2069  for(const auto &attr:onespsds->getAttributes()) {
2070  if ("_FillValue" == attr->getNewName()){
2071  has_fill_value = true;
2072  break;
2073  }
2074  }
2075 
2076 
2077  // Add fill value to some type of OBPG data. 16-bit integer, fill value = -32767, unsigned 16-bit integer fill value = 65535
2078  // This is based on the evaluation of the example files. KY 2012-09-06
2079  if ((false == has_fill_value) &&(DFNT_INT16 == onespsds->getType())) {
2080  short fill_value = -32767;
2081  string print_rep = HDFCFUtil::print_attr(DFNT_INT16,0,(void*)&fill_value);
2082  at->append_attr("_FillValue",HDFCFUtil::print_type(DFNT_INT16),print_rep);
2083  }
2084 
2085  if ((false == has_fill_value) &&(DFNT_UINT16 == onespsds->getType())) {
2086  unsigned short fill_value = 65535;
2087  string print_rep = HDFCFUtil::print_attr(DFNT_UINT16,0,(void*)&fill_value);
2088  at->append_attr("_FillValue",HDFCFUtil::print_type(DFNT_UINT16),print_rep);
2089  }
2090 
2091  }// Finish OBPG handling
2092 
2093 }
2094 
2095 // Handle HDF4 OTHERHDF products that follow SDS dimension scale model.
2096 // The special handling of AVHRR data is also included.
2097 void HDFCFUtil::handle_otherhdf_special_attrs(const HDFSP::File*f,DAS &das) {
2098 
2099  // For some HDF4 files that follow HDF4 dimension scales, P.O. DAAC's AVHRR files.
2100  // The "otherHDF" category can almost make AVHRR files work, except
2101  // that AVHRR uses the attribute name "unit" instead of "units" for latitude and longitude,
2102  // I have to correct the name to follow CF conventions(using "units"). I won't check
2103  // the latitude and longitude values since latitude and longitude values for some files(LISO files)
2104  // are not in the standard range(0-360 for lon and 0-180 for lat). KY 2011-3-3
2105  const vector<HDFSP::SDField *>& spsds = f->getSD()->getFields();
2106 
2107  if(f->getSPType() == OTHERHDF) {
2108 
2109  bool latflag = false;
2110  bool latunitsflag = false; //CF conventions use "units" instead of "unit"
2111  bool lonflag = false;
2112  bool lonunitsflag = false; // CF conventions use "units" instead of "unit"
2113  int llcheckoverflag = 0;
2114 
2115  // Here I try to condense the code within just two for loops
2116  // The outer loop: Loop through all SDS objects
2117  // The inner loop: Loop through all attributes of this SDS
2118  // Inside two inner loops(since "units" are common for other fields),
2119  // inner loop 1: when an attribute's long name value is L(l)atitude,mark it.
2120  // inner loop 2: when an attribute's name is units, mark it.
2121  // Outside the inner loop, when latflag is true, and latunitsflag is false,
2122  // adding a new attribute called units and the value should be "degrees_north".
2123  // doing the same thing for longitude.
2124 
2125  for (const auto &fd:spsds){
2126 
2127  // Ignore ALL coordinate variables if this is "OTHERHDF" case and some dimensions
2128  // don't have dimension scale data.
2129  if ( true == f->Has_Dim_NoScale_Field() && (fd->getFieldType() !=0) && (fd->IsDimScale() == false))
2130  continue;
2131 
2132  // Ignore the empty(no data) dimension variable.
2133  if (OTHERHDF == f->getSPType() && true == fd->IsDimNoScale())
2134  continue;
2135 
2136  AttrTable *at = das.get_table(fd->getNewName());
2137  if (!at)
2138  at = das.add_table(fd->getNewName(), new AttrTable);
2139 
2140  for (const auto& attr:fd->getAttributes()) {
2141  if(attr->getType()==DFNT_UCHAR || attr->getType() == DFNT_CHAR){
2142 
2143  if(attr->getName() == "long_name") {
2144  string tempstring2(attr->getValue().begin(),attr->getValue().end());
2145  string tempfinalstr= string(tempstring2.c_str());// This may remove some garbage characters
2146  if(tempfinalstr=="latitude" || tempfinalstr == "Latitude") // Find long_name latitude
2147  latflag = true;
2148  if(tempfinalstr=="longitude" || tempfinalstr == "Longitude") // Find long_name latitude
2149  lonflag = true;
2150  }
2151  }
2152  }
2153 
2154  if(latflag) {
2155  for (const auto& attr:fd->getAttributes()) {
2156  if (attr->getName() == "units")
2157  latunitsflag = true;
2158  }
2159  }
2160 
2161  if(lonflag) {
2162  for(const auto& attr:fd->getAttributes()) {
2163  if(attr->getName() == "units")
2164  lonunitsflag = true;
2165  }
2166  }
2167  if(latflag && !latunitsflag){ // No "units" for latitude, add "units"
2168  at->append_attr("units","String","degrees_north");
2169  latflag = false;
2170  latunitsflag = false;
2171  llcheckoverflag++;
2172  }
2173 
2174  if(lonflag && !lonunitsflag){ // No "units" for latitude, add "units"
2175  at->append_attr("units","String","degrees_east");
2176  lonflag = false;
2177  lonunitsflag = false;
2178  llcheckoverflag++;
2179  }
2180  if(llcheckoverflag ==2) break;
2181 
2182  }
2183 
2184  }
2185 
2186 }
2187 
2188 // Add missing CF attributes for non-CV varibles
2189 void
2190 HDFCFUtil::add_missing_cf_attrs(const HDFSP::File*f,DAS &das) {
2191 
2192  const vector<HDFSP::SDField *>& spsds = f->getSD()->getFields();
2193 
2194  // TRMM level 3 grid
2195  if(TRMML3A_V6== f->getSPType() || TRMML3C_V6==f->getSPType() || TRMML3S_V7 == f->getSPType() || TRMML3M_V7 == f->getSPType()) {
2196 
2197  for(const auto &fd:spsds){
2198  if(fd->getFieldType() == 0 && fd->getType()==DFNT_FLOAT32) {
2199 
2200  AttrTable *at = das.get_table(fd->getNewName());
2201  if (!at)
2202  at = das.add_table(fd->getNewName(), new AttrTable);
2203  string print_rep = "-9999.9";
2204  at->append_attr("_FillValue",HDFCFUtil::print_type(DFNT_FLOAT32),print_rep);
2205 
2206  }
2207  }
2208 
2209  for(const auto &fd:spsds){
2210  if(fd->getFieldType() == 0 && ((fd->getType()==DFNT_INT32) || (fd->getType()==DFNT_INT16))) {
2211 
2212  AttrTable *at = das.get_table(fd->getNewName());
2213  if (!at)
2214  at = das.add_table(fd->getNewName(), new AttrTable);
2215  string print_rep = "-9999";
2216  if(fd->getType()==DFNT_INT32)
2217  at->append_attr("_FillValue",HDFCFUtil::print_type(DFNT_INT32),print_rep);
2218  else
2219  at->append_attr("_FillValue",HDFCFUtil::print_type(DFNT_INT16),print_rep);
2220 
2221  }
2222  }
2223 
2224  // nlayer for TRMM single grid version 7, the units should be "km"
2225  if(TRMML3S_V7 == f->getSPType()) {
2226  for(const auto &fd:spsds){
2227  if(fd->getFieldType() == 6 && fd->getNewName()=="nlayer") {
2228 
2229  AttrTable *at = das.get_table(fd->getNewName());
2230  if (!at)
2231  at = das.add_table(fd->getNewName(), new AttrTable);
2232  at->append_attr("units","String","km");
2233 
2234  }
2235  else if(fd->getFieldType() == 4) {
2236 
2237  if (fd->getNewName()=="nh3" ||
2238  fd->getNewName()=="ncat3" ||
2239  fd->getNewName()=="nthrshZO" ||
2240  fd->getNewName()=="nthrshHB" ||
2241  fd->getNewName()=="nthrshSRT")
2242  {
2243 
2244  string references =
2245  "http://pps.gsfc.nasa.gov/Documents/filespec.TRMM.V7.pdf";
2246  string comment;
2247 
2248  AttrTable *at = das.get_table(fd->getNewName());
2249  if (!at)
2250  at = das.add_table(fd->getNewName(), new AttrTable);
2251 
2252  if(fd->getNewName()=="nh3") {
2253  comment="Index number to represent the fixed heights above the earth ellipsoid,";
2254  comment= comment + " at 2, 4, 6 km plus one for path-average.";
2255  }
2256 
2257  else if(fd->getNewName()=="ncat3") {
2258  comment="Index number to represent catgories for probability distribution functions.";
2259  comment=comment + "Check more information from the references.";
2260  }
2261 
2262  else if(fd->getNewName()=="nthrshZO")
2263  comment="Q-thresholds for Zero order used for probability distribution functions.";
2264 
2265  else if(fd->getNewName()=="nthrshHB")
2266  comment="Q-thresholds for HB used for probability distribution functions.";
2267 
2268  else if(fd->getNewName()=="nthrshSRT")
2269  comment="Q-thresholds for SRT used for probability distribution functions.";
2270 
2271  at->append_attr("comment","String",comment);
2272  at->append_attr("references","String",references);
2273 
2274  }
2275 
2276  }
2277 
2278  }
2279 
2280  // 3A26 use special values such as -666, -777,-999 in their fields.
2281  // Although the document doesn't provide range for some fields, the meaning of those fields should be greater than 0.
2282  // So add valid_min = 0 and fill_value = -999 .
2283  string base_filename;
2284  size_t last_slash_pos = f->getPath().find_last_of("/");
2285  if(last_slash_pos != string::npos)
2286  base_filename = f->getPath().substr(last_slash_pos+1);
2287  if(""==base_filename)
2288  base_filename = f->getPath();
2289  bool t3a26_flag = ((base_filename.find("3A26")!=string::npos)?true:false);
2290 
2291  if(true == t3a26_flag) {
2292 
2293  for(const auto &fd:spsds){
2294 
2295  if(fd->getFieldType() == 0 && (fd->getType()==DFNT_FLOAT32)) {
2296  AttrTable *at = das.get_table(fd->getNewName());
2297  if (!at)
2298  at = das.add_table(fd->getNewName(), new AttrTable);
2299  at->del_attr("_FillValue");
2300  at->append_attr("_FillValue","Float32","-999");
2301  at->append_attr("valid_min","Float32","0");
2302 
2303  }
2304  }
2305  }
2306 
2307  }
2308 
2309  // nlayer for TRMM single grid version 7, the units should be "km"
2310  if(TRMML3M_V7 == f->getSPType()) {
2311  for(const auto &fd:spsds){
2312 
2313  if(fd->getFieldType() == 4 ) {
2314 
2315  string references ="http://pps.gsfc.nasa.gov/Documents/filespec.TRMM.V7.pdf";
2316  if (fd->getNewName()=="nh1") {
2317 
2318  AttrTable *at = das.get_table(fd->getNewName());
2319  if (!at)
2320  at = das.add_table(fd->getNewName(), new AttrTable);
2321 
2322  string comment="Number of fixed heights above the earth ellipsoid,";
2323  comment= comment + " at 2, 4, 6, 10, and 15 km plus one for path-average.";
2324 
2325  at->append_attr("comment","String",comment);
2326  at->append_attr("references","String",references);
2327 
2328  }
2329  if (fd->getNewName()=="nh3") {
2330 
2331  AttrTable *at = das.get_table(fd->getNewName());
2332  if (!at)
2333  at = das.add_table(fd->getNewName(), new AttrTable);
2334 
2335  string comment="Number of fixed heights above the earth ellipsoid,";
2336  comment= comment + " at 2, 4, 6 km plus one for path-average.";
2337 
2338  at->append_attr("comment","String",comment);
2339  at->append_attr("references","String",references);
2340 
2341  }
2342 
2343  if (fd->getNewName()=="nang") {
2344 
2345  AttrTable *at = das.get_table(fd->getNewName());
2346  if (!at)
2347  at = das.add_table(fd->getNewName(), new AttrTable);
2348 
2349  string comment="Number of fixed incidence angles, at 0, 5, 10 and 15 degree and all angles.";
2350  references = "http://pps.gsfc.nasa.gov/Documents/ICSVol4.pdf";
2351 
2352  at->append_attr("comment","String",comment);
2353  at->append_attr("references","String",references);
2354 
2355  }
2356 
2357  if (fd->getNewName()=="ncat2") {
2358 
2359  AttrTable *at = das.get_table(fd->getNewName());
2360  if (!at)
2361  at = das.add_table(fd->getNewName(), new AttrTable);
2362 
2363  string comment="Second number of categories for histograms (30). ";
2364  comment=comment + "Check more information from the references.";
2365 
2366  at->append_attr("comment","String",comment);
2367  at->append_attr("references","String",references);
2368 
2369  }
2370 
2371  }
2372  }
2373 
2374  }
2375 
2376  }
2377 
2378  // TRMM level 2 swath
2379  else if(TRMML2_V7== f->getSPType()) {
2380 
2381  string base_filename;
2382  size_t last_slash_pos = f->getPath().find_last_of("/");
2383  if(last_slash_pos != string::npos)
2384  base_filename = f->getPath().substr(last_slash_pos+1);
2385  if(""==base_filename)
2386  base_filename = f->getPath();
2387  bool t2b31_flag = ((base_filename.find("2B31")!=string::npos)?true:false);
2388  bool t2a21_flag = ((base_filename.find("2A21")!=string::npos)?true:false);
2389  bool t2a12_flag = ((base_filename.find("2A12")!=string::npos)?true:false);
2390  // 2A23 is temporarily not supported perhaps due to special fill values
2391 #if 0
2392  //bool t2a23_flag = ((base_filename.find("2A23")!=string::npos)?true:false);
2393 #endif
2394  bool t2a25_flag = ((base_filename.find("2A25")!=string::npos)?true:false);
2395  bool t1c21_flag = ((base_filename.find("1C21")!=string::npos)?true:false);
2396  bool t1b21_flag = ((base_filename.find("1B21")!=string::npos)?true:false);
2397  bool t1b11_flag = ((base_filename.find("1B11")!=string::npos)?true:false);
2398  bool t1b01_flag = ((base_filename.find("1B01")!=string::npos)?true:false);
2399 
2400  // Handle scale and offset
2401 
2402  // group 1: 2B31,2A12,2A21
2403  if(t2b31_flag || t2a12_flag || t2a21_flag) {
2404 
2405  // special for 2B31
2406  if(t2b31_flag) {
2407 
2408  for (const auto &fd:spsds){
2409 
2410  if(fd->getFieldType() == 0 && fd->getType()==DFNT_INT16) {
2411 
2412  AttrTable *at = das.get_table(fd->getNewName());
2413  if (!at)
2414  at = das.add_table(fd->getNewName(), new AttrTable);
2415 
2416  AttrTable::Attr_iter it = at->attr_begin();
2417  while (it!=at->attr_end()) {
2418  if(at->get_name(it)=="scale_factor")
2419  {
2420  // Scale type and value
2421  string scale_factor_value="";
2422  string scale_factor_type;
2423 
2424  scale_factor_value = (*at->get_attr_vector(it)->begin());
2425  scale_factor_type = at->get_type(it);
2426 
2427  if(scale_factor_type == "Float64") {
2428  double new_scale = 1.0/strtod(scale_factor_value.c_str(),nullptr);
2429  at->del_attr("scale_factor");
2430  string print_rep = HDFCFUtil::print_attr(DFNT_FLOAT64,0,(void*)(&new_scale));
2431  at->append_attr("scale_factor", scale_factor_type,print_rep);
2432 
2433  }
2434 
2435  if(scale_factor_type == "Float32") {
2436  float new_scale = 1.0f/strtof(scale_factor_value.c_str(),nullptr);
2437  at->del_attr("scale_factor");
2438  string print_rep = HDFCFUtil::print_attr(DFNT_FLOAT32,0,(void*)(&new_scale));
2439  at->append_attr("scale_factor", scale_factor_type,print_rep);
2440 
2441  }
2442 
2443  break;
2444  }
2445  ++it;
2446  }
2447  }
2448  }
2449  }
2450 
2451  // Special for 2A12
2452  if (t2a12_flag==true) {
2453 
2454  for (const auto &fd:spsds){
2455 
2456  if (fd->getFieldType() == 6 && fd->getNewName()=="nlayer") {
2457 
2458  AttrTable *at = das.get_table(fd->getNewName());
2459  if (!at)
2460  at = das.add_table(fd->getNewName(), new AttrTable);
2461  at->append_attr("units","String","km");
2462 
2463  }
2464 
2465  // signed char maps to int32, so use int32 for the fillvalue.
2466  if (fd->getFieldType() == 0 && fd->getType()==DFNT_INT8) {
2467 
2468  AttrTable *at = das.get_table(fd->getNewName());
2469  if (!at)
2470  at = das.add_table(fd->getNewName(), new AttrTable);
2471  at->append_attr("_FillValue","Int32","-99");
2472 
2473  }
2474 
2475  }
2476  }
2477 
2478 
2479  // for all 2A12,2A21 and 2B31
2480  // Add fillvalues for float32 and int32.
2481  for (const auto & fd:spsds){
2482  if (fd->getFieldType() == 0 && fd->getType()==DFNT_FLOAT32) {
2483 
2484  AttrTable *at = das.get_table(fd->getNewName());
2485  if (!at)
2486  at = das.add_table(fd->getNewName(), new AttrTable);
2487  string print_rep = "-9999.9";
2488  at->append_attr("_FillValue",HDFCFUtil::print_type(DFNT_FLOAT32),print_rep);
2489 
2490  }
2491  }
2492 
2493  for (const auto &fd:spsds){
2494 
2495 
2496  if(fd->getFieldType() == 0 && fd->getType()==DFNT_INT16) {
2497 
2498  AttrTable *at = das.get_table(fd->getNewName());
2499  if (!at)
2500  at = das.add_table(fd->getNewName(), new AttrTable);
2501 
2502  string print_rep = "-9999";
2503  at->append_attr("_FillValue",HDFCFUtil::print_type(DFNT_INT32),print_rep);
2504 
2505  }
2506  }
2507 
2508  }
2509 
2510  // group 2: 2A21 and 2A25.
2511  else if(t2a21_flag == true || t2a25_flag == true) {
2512 
2513  // 2A25: handle reflectivity and rain rate scales
2514  if (t2a25_flag == true) {
2515 
2516  unsigned char handle_scale = 0;
2517 
2518  for (const auto &fd:spsds){
2519 
2520  if (fd->getFieldType() == 0 && fd->getType()==DFNT_INT16) {
2521  bool has_dBZ = false;
2522  bool has_rainrate = false;
2523  bool has_scale = false;
2524  string scale_factor_value;
2525  string scale_factor_type;
2526 
2527  AttrTable *at = das.get_table(fd->getNewName());
2528  if (!at)
2529  at = das.add_table(fd->getNewName(), new AttrTable);
2530  AttrTable::Attr_iter it = at->attr_begin();
2531  while (it!=at->attr_end()) {
2532  if(at->get_name(it)=="units"){
2533  string units_value = *at->get_attr_vector(it)->begin();
2534  if("dBZ" == units_value) {
2535  has_dBZ = true;
2536  }
2537 
2538  else if("mm/hr" == units_value){
2539  has_rainrate = true;
2540  }
2541  }
2542  if(at->get_name(it)=="scale_factor")
2543  {
2544  scale_factor_value = *at->get_attr_vector(it)->begin();
2545  scale_factor_type = at->get_type(it);
2546  has_scale = true;
2547  }
2548  ++it;
2549 
2550  }
2551 
2552  if((true == has_rainrate || true == has_dBZ) && true == has_scale) {
2553 
2554  handle_scale++;
2555  short valid_min = 0;
2556  short valid_max = 0;
2557 
2558  // Here just use 32-bit floating-point for the scale_factor, should be okay.
2559  if(true == has_rainrate)
2560  valid_max = (short)(300*strtof(scale_factor_value.c_str(),nullptr));
2561  else if(true == has_dBZ)
2562  valid_max = (short)(80*strtof(scale_factor_value.c_str(),nullptr));
2563 
2564  string print_rep1 = HDFCFUtil::print_attr(DFNT_INT16,0,(void*)(&valid_min));
2565  at->append_attr("valid_min","Int16",print_rep1);
2566  print_rep1 = HDFCFUtil::print_attr(DFNT_INT16,0,(void*)(&valid_max));
2567  at->append_attr("valid_max","Int16",print_rep1);
2568 
2569  at->del_attr("scale_factor");
2570  if(scale_factor_type == "Float64") {
2571  double new_scale = 1.0/strtod(scale_factor_value.c_str(),nullptr);
2572  string print_rep2 = HDFCFUtil::print_attr(DFNT_FLOAT64,0,(void*)(&new_scale));
2573  at->append_attr("scale_factor", scale_factor_type,print_rep2);
2574 
2575  }
2576 
2577  if(scale_factor_type == "Float32") {
2578  float new_scale = 1.0/strtof(scale_factor_value.c_str(),nullptr);
2579  string print_rep3 = HDFCFUtil::print_attr(DFNT_FLOAT32,0,(void*)(&new_scale));
2580  at->append_attr("scale_factor", scale_factor_type,print_rep3);
2581 
2582  }
2583 
2584 
2585  }
2586 
2587  if(2 == handle_scale)
2588  break;
2589 
2590  }
2591  }
2592  }
2593  }
2594 
2595  // 1B21,1C21 and 1B11
2596  else if (t1b21_flag || t1c21_flag || t1b11_flag) {
2597 
2598  // 1B21,1C21 scale_factor to CF and valid_range for dBm and dBZ.
2599  if (t1b21_flag || t1c21_flag) {
2600 
2601  for (const auto &fd:spsds){
2602 
2603  if(fd->getFieldType() == 0 && fd->getType()==DFNT_INT16) {
2604 
2605  bool has_dBm = false;
2606  bool has_dBZ = false;
2607 
2608  AttrTable *at = das.get_table(fd->getNewName());
2609  if (!at)
2610  at = das.add_table(fd->getNewName(), new AttrTable);
2611  AttrTable::Attr_iter it = at->attr_begin();
2612 
2613  while (it!=at->attr_end()) {
2614  if(at->get_name(it)=="units"){
2615 
2616  string units_value = *at->get_attr_vector(it)->begin();
2617  if("dBm" == units_value) {
2618  has_dBm = true;
2619  break;
2620  }
2621 
2622  else if("dBZ" == units_value){
2623  has_dBZ = true;
2624  break;
2625  }
2626  }
2627  ++it;
2628  }
2629 
2630  if(has_dBm == true || has_dBZ == true) {
2631  it = at->attr_begin();
2632  while (it!=at->attr_end()) {
2633  if(at->get_name(it)=="scale_factor")
2634  {
2635 
2636  string scale_value = *at->get_attr_vector(it)->begin();
2637 
2638  if(true == has_dBm) {
2639  short valid_min = (short)(-120 *strtof(scale_value.c_str(),nullptr));
2640  short valid_max = (short)(-20 *strtof(scale_value.c_str(),nullptr));
2641  string print_rep = HDFCFUtil::print_attr(DFNT_INT16,0,(void*)(&valid_min));
2642  at->append_attr("valid_min","Int16",print_rep);
2643  print_rep = HDFCFUtil::print_attr(DFNT_INT16,0,(void*)(&valid_max));
2644  at->append_attr("valid_max","Int16",print_rep);
2645  break;
2646 
2647  }
2648 
2649  else if(true == has_dBZ){
2650  short valid_min = (short)(-20 *strtof(scale_value.c_str(),nullptr));
2651  short valid_max = (short)(80 *strtof(scale_value.c_str(),nullptr));
2652  string print_rep = HDFCFUtil::print_attr(DFNT_INT16,0,(void*)(&valid_min));
2653  at->append_attr("valid_min","Int16",print_rep);
2654  print_rep = HDFCFUtil::print_attr(DFNT_INT16,0,(void*)(&valid_max));
2655  at->append_attr("valid_max","Int16",print_rep);
2656  break;
2657 
2658  }
2659  }
2660  ++it;
2661  }
2662 
2663  }
2664  }
2665  }
2666  }
2667 
2668  // For all 1B21,1C21 and 1B11 int16-bit products,change scale to follow CF
2669  // I find that one 1B21 variable binStormHeight has fillvalue -9999,
2670  // so add _FillValue -9999 for int16-bit variables.
2671  for(const auto &fd:spsds){
2672 
2673  if(fd->getFieldType() == 0 && fd->getType()==DFNT_INT16) {
2674 
2675  AttrTable *at = das.get_table(fd->getNewName());
2676  if (!at)
2677  at = das.add_table(fd->getNewName(), new AttrTable);
2678  AttrTable::Attr_iter it = at->attr_begin();
2679 
2680 
2681  while (it!=at->attr_end()) {
2682 
2683  if(at->get_name(it)=="scale_factor")
2684  {
2685  // Scale type and value
2686  string scale_factor_value="";
2687  string scale_factor_type;
2688 
2689  scale_factor_value = (*at->get_attr_vector(it)->begin());
2690  scale_factor_type = at->get_type(it);
2691 
2692  if(scale_factor_type == "Float64") {
2693  double new_scale = 1.0/strtod(scale_factor_value.c_str(),nullptr);
2694  at->del_attr("scale_factor");
2695  string print_rep = HDFCFUtil::print_attr(DFNT_FLOAT64,0,(void*)(&new_scale));
2696  at->append_attr("scale_factor", scale_factor_type,print_rep);
2697 
2698  }
2699 
2700  if(scale_factor_type == "Float32") {
2701  float new_scale = 1.0f/strtof(scale_factor_value.c_str(),nullptr);
2702  at->del_attr("scale_factor");
2703  string print_rep = HDFCFUtil::print_attr(DFNT_FLOAT32,0,(void*)(&new_scale));
2704  at->append_attr("scale_factor", scale_factor_type,print_rep);
2705 
2706  }
2707 
2708  break;
2709 
2710  }
2711  ++it;
2712 
2713  }
2714 
2715  at->append_attr("_FillValue","Int16","-9999");
2716 
2717  }
2718  }
2719  }
2720 
2721  // For 1B01 product, just add the fillvalue.
2722  else if (t1b01_flag == true) {
2723 
2724  for (const auto &fd:spsds){
2725 
2726  if(fd->getFieldType() == 0 && fd->getType()==DFNT_FLOAT32) {
2727 
2728  AttrTable *at = das.get_table(fd->getNewName());
2729  if (!at)
2730  at = das.add_table(fd->getNewName(), new AttrTable);
2731 
2732  at->append_attr("_FillValue","Float32","-9999.9");
2733  }
2734  }
2735  }
2736 
2737  AttrTable *at = das.get_table("HDF_GLOBAL");
2738  if (!at)
2739  at = das.add_table("HDF_GLOBAL", new AttrTable);
2740  string references ="http://pps.gsfc.nasa.gov/Documents/filespec.TRMM.V7.pdf";
2741  string comment="The HDF4 OPeNDAP handler adds _FillValue, valid_min and valid_max for some TRMM level 1 and level 2 products.";
2742  comment= comment + " It also changes scale_factor to follow CF conventions. ";
2743 
2744  at->append_attr("comment","String",comment);
2745  at->append_attr("references","String",references);
2746 
2747  }
2748 
2749 }
2750 
2751 //
2752 // Many CERES products compose of multiple groups
2753 // There are many fields in CERES data(a few hundred) and the full name(with the additional path)
2754 // is very long. It causes Java clients choken since Java clients append names in the URL
2755 // To improve the performance and to make Java clients access the data, we will ignore
2756 // the path of these fields. Users can turn off this feature by commenting out the line: H4.EnableCERESMERRAShortName=true
2757 // or can set the H4.EnableCERESMERRAShortName=false
2758 // We still preserve the full path as an attribute in case users need to check them.
2759 // Kent 2012-6-29
2760 
2761 void HDFCFUtil::handle_merra_ceres_attrs_with_bes_keys(const HDFSP::File *f, DAS &das,const string& filename) {
2762 
2763  string base_filename = filename.substr(filename.find_last_of("/")+1);
2764 
2765 #if 0
2766  string check_ceres_merra_short_name_key="H4.EnableCERESMERRAShortName";
2767  bool turn_on_ceres_merra_short_name_key= false;
2768 
2769  turn_on_ceres_merra_short_name_key = HDFCFUtil::check_beskeys(check_ceres_merra_short_name_key);
2770 #endif
2771 
2772  bool merra_is_eos2 = false;
2773  if(0== (base_filename.compare(0,5,"MERRA"))) {
2774 
2775 #if 0
2776  for (vector < HDFSP::Attribute * >::const_iterator i =
2777  f->getSD()->getAttributes ().begin ();
2778  i != f->getSD()->getAttributes ().end (); ++i) {
2779 #endif
2780  for (const auto & fd:f->getSD()->getAttributes()) {
2781 
2782  // Check if this MERRA file is an HDF-EOS2 or not.
2783  if((fd->getName().compare(0, 14, "StructMetadata" )== 0) ||
2784  (fd->getName().compare(0, 14, "structmetadata" )== 0)) {
2785  merra_is_eos2 = true;
2786  break;
2787  }
2788 
2789  }
2790  }
2791 
2792  if (true == HDF4RequestHandler::get_enable_ceres_merra_short_name() && (CER_ES4 == f->getSPType() || CER_SRB == f->getSPType()
2793  || CER_CDAY == f->getSPType() || CER_CGEO == f->getSPType()
2794  || CER_SYN == f->getSPType() || CER_ZAVG == f->getSPType()
2795  || CER_AVG == f->getSPType() || (true == merra_is_eos2))) {
2796 
2797  const vector<HDFSP::SDField *>& spsds = f->getSD()->getFields();
2798  for (const auto & fd:spsds){
2799 
2800  AttrTable *at = das.get_table(fd->getNewName());
2801  if (!at)
2802  at = das.add_table(fd->getNewName(), new AttrTable);
2803 
2804  at->append_attr("fullpath","String",fd->getSpecFullPath());
2805 
2806  }
2807 
2808  }
2809 
2810 }
2811 
2812 
2813 // Handle the attributes when the BES key EnableVdataDescAttr is enabled..
2814 void HDFCFUtil::handle_vdata_attrs_with_desc_key(const HDFSP::File*f,libdap::DAS &das) {
2815 
2816  // Check the EnableVdataDescAttr key. If this key is turned on, the handler-added attribute VDdescname and
2817  // the attributes of vdata and vdata fields will be outputed to DAS. Otherwise, these attributes will
2818  // not outputed to DAS. The key will be turned off by default to shorten the DAP output. KY 2012-09-18
2819 
2820 #if 0
2821  string check_vdata_desc_key="H4.EnableVdataDescAttr";
2822  bool turn_on_vdata_desc_key= false;
2823 
2824  turn_on_vdata_desc_key = HDFCFUtil::check_beskeys(check_vdata_desc_key);
2825 #endif
2826 
2827  string VDdescname = "hdf4_vd_desc";
2828  string VDdescvalue = "This is an HDF4 Vdata.";
2829  string VDfieldprefix = "Vdata_field_";
2830  string VDattrprefix = "Vdata_attr_";
2831  string VDfieldattrprefix ="Vdata_field_attr_";
2832 
2833  // To speed up the performance for handling CERES data, we turn off some CERES vdata fields, this should be resumed in the future version with BESKeys.
2834 #if 0
2835  string check_ceres_vdata_key="H4.EnableCERESVdata";
2836  bool turn_on_ceres_vdata_key= false;
2837  turn_on_ceres_vdata_key = HDFCFUtil::check_beskeys(check_ceres_vdata_key);
2838 #endif
2839 
2840  bool output_vdata_flag = true;
2841  if (false == HDF4RequestHandler::get_enable_ceres_vdata() &&
2842  (CER_AVG == f->getSPType() ||
2843  CER_ES4 == f->getSPType() ||
2844  CER_SRB == f->getSPType() ||
2845  CER_ZAVG == f->getSPType()))
2846  output_vdata_flag = false;
2847 
2848  if (true == output_vdata_flag) {
2849 
2850  for (const auto &vd:f->getVDATAs()) {
2851 
2852  AttrTable *at = das.get_table(vd->getNewName());
2853  if(!at)
2854  at = das.add_table(vd->getNewName(),new AttrTable);
2855 
2856  if (true == HDF4RequestHandler::get_enable_vdata_desc_attr()) {
2857  // Add special vdata attributes
2858  bool emptyvddasflag = true;
2859  if (!(vd->getAttributes().empty()))
2860  emptyvddasflag = false;
2861  if (vd->getTreatAsAttrFlag())
2862  emptyvddasflag = false;
2863  else {
2864  for (const auto &vfd:vd->getFields()) {
2865  if(!(vfd->getAttributes().empty())) {
2866  emptyvddasflag = false;
2867  break;
2868  }
2869  }
2870  }
2871 
2872  if(emptyvddasflag)
2873  continue;
2874  at->append_attr(VDdescname, "String" , VDdescvalue);
2875 
2876  for(const auto &va:vd->getAttributes()) {
2877 
2878  if(va->getType()==DFNT_UCHAR || va->getType() == DFNT_CHAR){
2879 
2880  string tempstring2(va->getValue().begin(),va->getValue().end());
2881  string tempfinalstr= string(tempstring2.c_str());
2882  at->append_attr(VDattrprefix+va->getNewName(), "String" , HDFCFUtil::escattr(tempfinalstr));
2883  }
2884  else {
2885  for (int loc=0; loc < va->getCount() ; loc++) {
2886  string print_rep = HDFCFUtil::print_attr(va->getType(), loc, (void*) &(va->getValue()[0]));
2887  at->append_attr(VDattrprefix+va->getNewName(), HDFCFUtil::print_type(va->getType()), print_rep);
2888  }
2889  }
2890  }
2891  }
2892 
2893  if(false == (vd->getTreatAsAttrFlag())){
2894 
2895  if (true == HDF4RequestHandler::get_enable_vdata_desc_attr()) {
2896 
2897  //NOTE: for vdata field, we assume that no special characters are found. We need to escape the special characters when the data type is char.
2898  // We need to create a DAS container for each field so that the attributes can be put inside.
2899  for (const auto &vdf:vd->getFields()) {
2900 
2901  // This vdata field will NOT be treated as attributes, only save the field attribute as the attribute
2902  // First check if the field has attributes, if it doesn't have attributes, no need to create a container.
2903 
2904  if (vdf->getAttributes().empty() ==false) {
2905 
2906  AttrTable *at_v = das.get_table(vdf->getNewName());
2907  if(!at_v)
2908  at_v = das.add_table(vdf->getNewName(),new AttrTable);
2909 
2910  for (const auto &va:vdf->getAttributes()) {
2911 
2912  if(va->getType()==DFNT_UCHAR || va->getType() == DFNT_CHAR){
2913 
2914  string tempstring2(va->getValue().begin(),va->getValue().end());
2915  string tempfinalstr= string(tempstring2.c_str());
2916  at_v->append_attr(va->getNewName(), "String" , HDFCFUtil::escattr(tempfinalstr));
2917  }
2918  else {
2919  for (int loc=0; loc < va->getCount() ; loc++) {
2920  string print_rep = HDFCFUtil::print_attr(va->getType(), loc, (void*) &(va->getValue()[0]));
2921  at_v->append_attr(va->getNewName(), HDFCFUtil::print_type(va->getType()), print_rep);
2922  }
2923  }
2924 
2925  }
2926  }
2927  }
2928  }
2929 
2930  }
2931 
2932  else {
2933 
2934  for(const auto & vdf:vd->getFields()) {
2935 
2936  if(vdf->getFieldOrder() == 1) {
2937  if(vdf->getType()==DFNT_UCHAR || vdf->getType() == DFNT_CHAR){
2938  string tempfinalstr;
2939  tempfinalstr.resize(vdf->getValue().size());
2940  copy(vdf->getValue().begin(),vdf->getValue().end(),tempfinalstr.begin());
2941  at->append_attr(VDfieldprefix+vdf->getNewName(), "String" , HDFCFUtil::escattr(tempfinalstr));
2942  }
2943  else {
2944  for ( int loc=0; loc < vdf->getNumRec(); loc++) {
2945  string print_rep = HDFCFUtil::print_attr(vdf->getType(), loc, (void*) &(vdf->getValue()[0]));
2946  at->append_attr(VDfieldprefix+vdf->getNewName(), HDFCFUtil::print_type(vdf->getType()), print_rep);
2947  }
2948  }
2949  }
2950  else {//When field order is greater than 1,we want to print each record in group with single quote,'0 1 2','3 4 5', etc.
2951 
2952  if (vdf->getValue().size() != (unsigned int)(DFKNTsize(vdf->getType())*(vdf->getFieldOrder())*(vdf->getNumRec()))){
2953  throw InternalErr(__FILE__,__LINE__,"the vdata field size doesn't match the vector value");
2954  }
2955 
2956  if(vdf->getNumRec()==1){
2957  if(vdf->getType()==DFNT_UCHAR || vdf->getType() == DFNT_CHAR){
2958  string tempstring2(vdf->getValue().begin(),vdf->getValue().end());
2959  string tempfinalstr= string(tempstring2.c_str());
2960  at->append_attr(VDfieldprefix+vdf->getNewName(),"String",HDFCFUtil::escattr(tempfinalstr));
2961  }
2962  else {
2963  for (int loc=0; loc < vdf->getFieldOrder(); loc++) {
2964  string print_rep = HDFCFUtil::print_attr(vdf->getType(), loc, (void*) &(vdf->getValue()[0]));
2965  at->append_attr(VDfieldprefix+vdf->getNewName(), HDFCFUtil::print_type(vdf->getType()), print_rep);
2966  }
2967  }
2968 
2969  }
2970 
2971  else {
2972  if(vdf->getType()==DFNT_UCHAR || vdf->getType() == DFNT_CHAR){
2973 
2974  for(int tempcount = 0; tempcount < vdf->getNumRec()*DFKNTsize(vdf->getType());tempcount ++) {
2975  vector<char>::const_iterator tempit;
2976  tempit = vdf->getValue().begin()+tempcount*(vdf->getFieldOrder());
2977  string tempstring2(tempit,tempit+vdf->getFieldOrder());
2978  string tempfinalstr= string(tempstring2.c_str());
2979  string tempoutstring = "'"+tempfinalstr+"'";
2980  at->append_attr(VDfieldprefix+vdf->getNewName(),"String",HDFCFUtil::escattr(tempoutstring));
2981  }
2982  }
2983 
2984  else {
2985  for(int tempcount = 0; tempcount < vdf->getNumRec();tempcount ++) {
2986  at->append_attr(VDfieldprefix+vdf->getNewName(),HDFCFUtil::print_type(vdf->getType()),"'");
2987  for (int loc=0; loc < vdf->getFieldOrder(); loc++) {
2988  string print_rep = HDFCFUtil::print_attr(vdf->getType(), loc, (void*) &(vdf->getValue()[tempcount*(vdf->getFieldOrder())]));
2989  at->append_attr(VDfieldprefix+vdf->getNewName(), HDFCFUtil::print_type(vdf->getType()), print_rep);
2990  }
2991  at->append_attr(VDfieldprefix+vdf->getNewName(),HDFCFUtil::print_type(vdf->getType()),"'");
2992  }
2993  }
2994  }
2995  }
2996 
2997 
2998  if (true == HDF4RequestHandler::get_enable_vdata_desc_attr()) {
2999  for(const auto &va:vdf->getAttributes()) {
3000 
3001  if(va->getType()==DFNT_UCHAR || va->getType() == DFNT_CHAR){
3002 
3003  string tempstring2(va->getValue().begin(),va->getValue().end());
3004  string tempfinalstr= string(tempstring2.c_str());
3005  at->append_attr(VDfieldattrprefix+va->getNewName(), "String" , HDFCFUtil::escattr(tempfinalstr));
3006  }
3007  else {
3008  for (int loc=0; loc < va->getCount() ; loc++) {
3009  string print_rep = HDFCFUtil::print_attr(va->getType(), loc, (void*) &(va->getValue()[0]));
3010  at->append_attr(VDfieldattrprefix+va->getNewName(), HDFCFUtil::print_type(va->getType()), print_rep);
3011  }
3012  }
3013  }
3014  }
3015  }
3016  }
3017 
3018  }
3019  }
3020 
3021 }
3022 
3023 void HDFCFUtil::map_eos2_objects_attrs(libdap::DAS &das,const string &filename) {
3024 
3025  intn status_n =-1;
3026  int32 status_32 = -1;
3027  int32 file_id = -1;
3028  int32 vgroup_id = -1;
3029  int32 num_of_lones = 0;
3030  uint16 name_len = 0;
3031 
3032  file_id = Hopen (filename.c_str(), DFACC_READ, 0);
3033  if(file_id == FAIL)
3034  throw InternalErr(__FILE__,__LINE__,"Hopen failed.");
3035 
3036  status_n = Vstart (file_id);
3037  if(status_n == FAIL) {
3038  Hclose(file_id);
3039  throw InternalErr(__FILE__,__LINE__,"Vstart failed.");
3040  }
3041 
3042  string err_msg;
3043  bool unexpected_fail = false;
3044  //Get and print the names and class names of all the lone vgroups.
3045  // First, call Vlone with num_of_lones set to 0 to get the number of
3046  // lone vgroups in the file, but not to get their reference numbers.
3047  num_of_lones = Vlone (file_id, nullptr, num_of_lones );
3048 
3049  //
3050  // Then, if there are any lone vgroups,
3051  if (num_of_lones > 0)
3052  {
3053  // Use the num_of_lones returned to allocate sufficient space for the
3054  // buffer ref_array to hold the reference numbers of all lone vgroups,
3055  vector<int32> ref_array;
3056  ref_array.resize(num_of_lones);
3057 
3058 
3059  // and call Vlone again to retrieve the reference numbers into
3060  // the buffer ref_array.
3061 
3062  num_of_lones = Vlone (file_id, ref_array.data(), num_of_lones);
3063 
3064  // Loop the name and class of each lone vgroup.
3065  for (int lone_vg_number = 0; lone_vg_number < num_of_lones;
3066  lone_vg_number++)
3067  {
3068 
3069  // Attach to the current vgroup then get and display its
3070  // name and class. Note: the current vgroup must be detached before
3071  // moving to the next.
3072  vgroup_id = Vattach(file_id, ref_array[lone_vg_number], "r");
3073  if(vgroup_id == FAIL) {
3074  unexpected_fail = true;
3075  err_msg = string(ERR_LOC) + " Vattach failed. ";
3076  goto cleanFail;
3077  }
3078 
3079  status_32 = Vgetnamelen(vgroup_id, &name_len);
3080  if(status_32 == FAIL) {
3081  unexpected_fail = true;
3082  Vdetach(vgroup_id);
3083  err_msg = string(ERR_LOC) + " Vgetnamelen failed. ";
3084  goto cleanFail;
3085  }
3086 
3087  vector<char> vgroup_name;
3088  vgroup_name.resize(name_len+1);
3089  status_32 = Vgetname (vgroup_id, vgroup_name.data());
3090  if(status_32 == FAIL) {
3091  unexpected_fail = true;
3092  Vdetach(vgroup_id);
3093  err_msg = string(ERR_LOC) + " Vgetname failed. ";
3094  goto cleanFail;
3095  }
3096 
3097  status_32 = Vgetclassnamelen(vgroup_id, &name_len);
3098  if(status_32 == FAIL) {
3099  unexpected_fail = true;
3100  Vdetach(vgroup_id);
3101  err_msg = string(ERR_LOC) + " Vgetclassnamelen failed. ";
3102  goto cleanFail;
3103  }
3104 
3105  vector<char>vgroup_class;
3106  vgroup_class.resize(name_len+1);
3107  status_32 = Vgetclass (vgroup_id, vgroup_class.data());
3108  if(status_32 == FAIL) {
3109  unexpected_fail = true;
3110  Vdetach(vgroup_id);
3111  err_msg = string(ERR_LOC) + " Vgetclass failed. ";
3112  goto cleanFail;
3113  }
3114 
3115  string vgroup_name_str(vgroup_name.begin(),vgroup_name.end());
3116  vgroup_name_str = vgroup_name_str.substr(0,vgroup_name_str.size()-1);
3117 
3118  string vgroup_class_str(vgroup_class.begin(),vgroup_class.end());
3119  vgroup_class_str = vgroup_class_str.substr(0,vgroup_class_str.size()-1);
3120  try {
3121  if(vgroup_class_str =="GRID")
3122  map_eos2_one_object_attrs_wrapper(das,file_id,vgroup_id,vgroup_name_str,true);
3123  else if(vgroup_class_str =="SWATH")
3124  map_eos2_one_object_attrs_wrapper(das,file_id,vgroup_id,vgroup_name_str,false);
3125  }
3126  catch(...) {
3127  Vdetach(vgroup_id);
3128  Vend(file_id);
3129  Hclose(file_id);
3130  throw InternalErr(__FILE__,__LINE__,"map_eos2_one_object_attrs_wrapper failed.");
3131  }
3132  Vdetach (vgroup_id);
3133  }// for
3134  }// if
3135 
3136  //Terminate access to the V interface and close the file.
3137 cleanFail:
3138  Vend (file_id);
3139  Hclose (file_id);
3140  if(true == unexpected_fail)
3141  throw InternalErr(__FILE__,__LINE__,err_msg);
3142 
3143  return;
3144 
3145 }
3146 
3147 void HDFCFUtil::map_eos2_one_object_attrs_wrapper(libdap:: DAS &das,int32 file_id,int32 vgroup_id, const string& vgroup_name,bool is_grid) {
3148 
3149  int32 num_gobjects = Vntagrefs (vgroup_id);
3150  if(num_gobjects < 0)
3151  throw InternalErr(__FILE__,__LINE__,"Cannot obtain the number of objects under a vgroup.");
3152 
3153  for(int i = 0; i<num_gobjects;i++) {
3154 
3155  int32 obj_tag;
3156  int32 obj_ref;
3157  if (Vgettagref (vgroup_id, i, &obj_tag, &obj_ref) == FAIL)
3158  throw InternalErr(__FILE__,__LINE__,"Failed to obtain the tag and reference of an object under a vgroup.");
3159 
3160  if (Visvg (vgroup_id, obj_ref) == TRUE) {
3161 
3162  int32 object_attr_vgroup = Vattach(file_id,obj_ref,"r");
3163  if(object_attr_vgroup == FAIL)
3164  throw InternalErr(__FILE__,__LINE__,"Failed to attach an EOS2 vgroup.");
3165 
3166  uint16 name_len = 0;
3167  int32 status_32 = Vgetnamelen(object_attr_vgroup, &name_len);
3168  if(status_32 == FAIL) {
3169  Vdetach(object_attr_vgroup);
3170  throw InternalErr(__FILE__,__LINE__,"Failed to obtain an EOS2 vgroup name length.");
3171  }
3172  vector<char> attr_vgroup_name;
3173  attr_vgroup_name.resize(name_len+1);
3174  status_32 = Vgetname (object_attr_vgroup, attr_vgroup_name.data());
3175  if(status_32 == FAIL) {
3176  Vdetach(object_attr_vgroup);
3177  throw InternalErr(__FILE__,__LINE__,"Failed to obtain an EOS2 vgroup name. ");
3178  }
3179 
3180  string attr_vgroup_name_str(attr_vgroup_name.begin(),attr_vgroup_name.end());
3181  attr_vgroup_name_str = attr_vgroup_name_str.substr(0,attr_vgroup_name_str.size()-1);
3182 
3183  try {
3184  if(true == is_grid && attr_vgroup_name_str=="Grid Attributes"){
3185  map_eos2_one_object_attrs(das,file_id,object_attr_vgroup,vgroup_name);
3186  Vdetach(object_attr_vgroup);
3187  break;
3188  }
3189  else if(false == is_grid && attr_vgroup_name_str=="Swath Attributes") {
3190  map_eos2_one_object_attrs(das,file_id,object_attr_vgroup,vgroup_name);
3191  Vdetach(object_attr_vgroup);
3192  break;
3193  }
3194  }
3195  catch(...) {
3196  Vdetach(object_attr_vgroup);
3197  throw InternalErr(__FILE__,__LINE__,"Cannot map eos2 object attributes to DAP2.");
3198  }
3199  Vdetach(object_attr_vgroup);
3200 
3201  }
3202 
3203  }
3204 }
3205 
3206 void HDFCFUtil::map_eos2_one_object_attrs(libdap:: DAS &das,int32 file_id, int32 obj_attr_group_id, const string& vgroup_name) {
3207 
3208  int32 num_gobjects = Vntagrefs(obj_attr_group_id);
3209  if(num_gobjects < 0)
3210  throw InternalErr(__FILE__,__LINE__,"Cannot obtain the number of objects under a vgroup.");
3211 
3212  string vgroup_cf_name = HDFCFUtil::get_CF_string(vgroup_name);
3213  AttrTable *at = das.get_table(vgroup_cf_name);
3214  if(!at)
3215  at = das.add_table(vgroup_cf_name,new AttrTable);
3216 
3217  for(int i = 0; i<num_gobjects;i++) {
3218 
3219  int32 obj_tag, obj_ref;
3220  if (Vgettagref(obj_attr_group_id, i, &obj_tag, &obj_ref) == FAIL)
3221  throw InternalErr(__FILE__,__LINE__,"Failed to obtain the tag and reference of an object under a vgroup.");
3222 
3223  if(Visvs(obj_attr_group_id,obj_ref)) {
3224 
3225  int32 vdata_id = VSattach(file_id,obj_ref,"r");
3226  if(vdata_id == FAIL)
3227  throw InternalErr(__FILE__,__LINE__,"Failed to attach a vdata.");
3228 
3229  // EOS2 object vdatas are actually EOS2 object attributes.
3230  if(VSisattr(vdata_id)) {
3231 
3232  // According to EOS2 library, EOS2 number of field and record must be 1.
3233  int32 num_field = VFnfields(vdata_id);
3234  if(num_field !=1) {
3235  VSdetach(vdata_id);
3236  throw InternalErr(__FILE__,__LINE__,"Number of vdata field for an EOS2 object must be 1.");
3237  }
3238  int32 num_record = VSelts(vdata_id);
3239  if(num_record !=1){
3240  VSdetach(vdata_id);
3241  throw InternalErr(__FILE__,__LINE__,"Number of vdata record for an EOS2 object must be 1.");
3242  }
3243  char vdata_name[VSNAMELENMAX];
3244  if(VSQueryname(vdata_id,vdata_name) == FAIL) {
3245  VSdetach(vdata_id);
3246  throw InternalErr(__FILE__,__LINE__,"Failed to obtain EOS2 object vdata name.");
3247  }
3248  string vdatanamestr(vdata_name);
3249  string vdataname_cfstr = HDFCFUtil::get_CF_string(vdatanamestr);
3250  int32 fieldsize = VFfieldesize(vdata_id,0);
3251  if(fieldsize == FAIL) {
3252  VSdetach(vdata_id);
3253  throw InternalErr(__FILE__,__LINE__,"Failed to obtain EOS2 object vdata field size.");
3254  }
3255 
3256  const char* fieldname = VFfieldname(vdata_id,0);
3257  if(fieldname == nullptr) {
3258  VSdetach(vdata_id);
3259  throw InternalErr(__FILE__,__LINE__,"Failed to obtain EOS2 object vdata field name.");
3260  }
3261  int32 fieldtype = VFfieldtype(vdata_id,0);
3262  if(fieldtype == FAIL) {
3263  VSdetach(vdata_id);
3264  throw InternalErr(__FILE__,__LINE__,"Failed to obtain EOS2 object vdata field type.");
3265  }
3266 
3267  if(VSsetfields(vdata_id,fieldname) == FAIL) {
3268  VSdetach(vdata_id);
3269  throw InternalErr(__FILE__,__LINE__,"EOS2 object vdata: VSsetfields failed.");
3270  }
3271 
3272  vector<char> vdata_value;
3273  vdata_value.resize(fieldsize);
3274  if(VSread(vdata_id,(uint8*)vdata_value.data(),1,FULL_INTERLACE) == FAIL) {
3275  VSdetach(vdata_id);
3276  throw InternalErr(__FILE__,__LINE__,"EOS2 object vdata: VSread failed.");
3277  }
3278 
3279  // Map the attributes to DAP2.
3280  if(fieldtype == DFNT_UCHAR || fieldtype == DFNT_CHAR){
3281  string tempstring(vdata_value.begin(),vdata_value.end());
3282  // Remove the nullptr term
3283  auto tempstring2 = string(tempstring.c_str());
3284  at->append_attr(vdataname_cfstr,"String",HDFCFUtil::escattr(tempstring2));
3285  }
3286  else {
3287  string print_rep = HDFCFUtil::print_attr(fieldtype, 0, (void*) vdata_value.data());
3288  at->append_attr(vdataname_cfstr, HDFCFUtil::print_type(fieldtype), print_rep);
3289  }
3290 
3291  }
3292  VSdetach(vdata_id);
3293 
3294  }
3295  }
3296 
3297  return;
3298 }
3299 
3300 // Part of a large fix for attributes. Escaping the values of the attributes
3301 // may have been a bad idea. It breaks using JSON, for example. If this is a
3302 // bad idea - to turn of escaping - then we'll have to figure out how to store
3303 // 'serialized JSON' in attributes because it's being used in netcdf/hdf files.
3304 // If we stick with this, there's clearly a more performant solution - eliminate
3305 // the calls to this code.
3306 // jhrg 6/25/21
3307 #define ESCAPE_STRING_ATTRIBUTES 0
3308 
3309 string HDFCFUtil::escattr(string s)
3310 {
3311  const string printable = " ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789~`!@#$%^&*()_-+={[}]|\\:;<,>.?/'\"\n\t\r";
3312  const string ESC = "\\";
3313 #if ESCAPE_STRING_ATTRIBUTES
3314  const string DOUBLE_ESC = ESC + ESC;
3315  const string QUOTE = "\"";
3316  const string ESCQUOTE = ESC + QUOTE;
3317 
3318  // escape \ with a second backslash
3319  size_t ind = 0;
3320  while ((ind = s.find(ESC, ind)) != string::npos) {
3321  s.replace(ind, 1, DOUBLE_ESC);
3322  ind += DOUBLE_ESC.length();
3323  }
3324 
3325  // escape " with backslash
3326  ind = 0;
3327  while ((ind = s.find(QUOTE, ind)) != string::npos) {
3328  //comment out the following line since it wastes the CPU operation.
3329  s.replace(ind, 1, ESCQUOTE);
3330  ind += ESCQUOTE.length();
3331  }
3332 #endif
3333 
3334  size_t ind = 0;
3335  while ((ind = s.find_first_not_of(printable, ind)) != string::npos) {
3336  s.replace(ind, 1, ESC + octstring(s[ind]));
3337  }
3338 
3339  return s;
3340 }
3341 
3342 
3343 void HDFCFUtil::parser_trmm_v7_gridheader(const vector<char>& value,
3344  int& latsize, int&lonsize,
3345  float& lat_start, float& lon_start,
3346  float& lat_res, float& lon_res,
3347  bool check_reg_orig ){
3348 
3349  float lat_north = 0.;
3350  float lat_south = 0.;
3351  float lon_east = 0.;
3352  float lon_west = 0.;
3353 
3354  vector<string> ind_elems;
3355  char sep='\n';
3356  HDFCFUtil::Split(value.data(),sep,ind_elems);
3357  // The number of elements in the GridHeader is 9.
3358  //The string vector will add a leftover. So the size should be 10.
3359  // For the MacOS clang compiler, the string vector size may become 11.
3360  // So we change the condition to be "<9" is wrong.
3361  if(ind_elems.size()<9)
3362  throw InternalErr(__FILE__,__LINE__,"The number of elements in the TRMM level 3 GridHeader is not right.");
3363 
3364  if(false == check_reg_orig) {
3365  if (0 != ind_elems[1].find("Registration=CENTER"))
3366  throw InternalErr(__FILE__,__LINE__,"The TRMM grid registration is not center.");
3367  }
3368 
3369  if (0 == ind_elems[2].find("LatitudeResolution")){
3370 
3371  size_t equal_pos = ind_elems[2].find_first_of('=');
3372  if(string::npos == equal_pos)
3373  throw InternalErr(__FILE__,__LINE__,"Cannot find latitude resolution for TRMM level 3 products");
3374 
3375  size_t scolon_pos = ind_elems[2].find_first_of(';');
3376  if(string::npos == scolon_pos)
3377  throw InternalErr(__FILE__,__LINE__,"Cannot find latitude resolution for TRMM level 3 products");
3378  if (equal_pos < scolon_pos){
3379 
3380  string latres_str = ind_elems[2].substr(equal_pos+1,scolon_pos-equal_pos-1);
3381  lat_res = strtof(latres_str.c_str(),nullptr);
3382  }
3383  else
3384  throw InternalErr(__FILE__,__LINE__,"latitude resolution is not right for TRMM level 3 products");
3385  }
3386  else
3387  throw InternalErr(__FILE__,__LINE__,"The TRMM grid LatitudeResolution doesn't exist.");
3388 
3389  if (0 == ind_elems[3].find("LongitudeResolution")){
3390 
3391  size_t equal_pos = ind_elems[3].find_first_of('=');
3392  if(string::npos == equal_pos)
3393  throw InternalErr(__FILE__,__LINE__,"Cannot find longitude resolution for TRMM level 3 products");
3394 
3395  size_t scolon_pos = ind_elems[3].find_first_of(';');
3396  if(string::npos == scolon_pos)
3397  throw InternalErr(__FILE__,__LINE__,"Cannot find longitude resolution for TRMM level 3 products");
3398  if (equal_pos < scolon_pos){
3399  string lonres_str = ind_elems[3].substr(equal_pos+1,scolon_pos-equal_pos-1);
3400  lon_res = strtof(lonres_str.c_str(),nullptr);
3401  }
3402  else
3403  throw InternalErr(__FILE__,__LINE__,"longitude resolution is not right for TRMM level 3 products");
3404  }
3405  else
3406  throw InternalErr(__FILE__,__LINE__,"The TRMM grid LongitudeResolution doesn't exist.");
3407 
3408  if (0 == ind_elems[4].find("NorthBoundingCoordinate")){
3409 
3410  size_t equal_pos = ind_elems[4].find_first_of('=');
3411  if(string::npos == equal_pos)
3412  throw InternalErr(__FILE__,__LINE__,"Cannot find latitude resolution for TRMM level 3 products");
3413 
3414  size_t scolon_pos = ind_elems[4].find_first_of(';');
3415  if(string::npos == scolon_pos)
3416  throw InternalErr(__FILE__,__LINE__,"Cannot find latitude resolution for TRMM level 3 products");
3417  if (equal_pos < scolon_pos){
3418  string north_bounding_str = ind_elems[4].substr(equal_pos+1,scolon_pos-equal_pos-1);
3419  lat_north = strtof(north_bounding_str.c_str(),nullptr);
3420  }
3421  else
3422  throw InternalErr(__FILE__,__LINE__,"NorthBoundingCoordinate is not right for TRMM level 3 products");
3423 
3424  }
3425  else
3426  throw InternalErr(__FILE__,__LINE__,"The TRMM grid NorthBoundingCoordinate doesn't exist.");
3427 
3428  if (0 == ind_elems[5].find("SouthBoundingCoordinate")){
3429 
3430  size_t equal_pos = ind_elems[5].find_first_of('=');
3431  if(string::npos == equal_pos)
3432  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for TRMM level 3 products");
3433 
3434  size_t scolon_pos = ind_elems[5].find_first_of(';');
3435  if(string::npos == scolon_pos)
3436  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for TRMM level 3 products");
3437  if (equal_pos < scolon_pos){
3438  string lat_south_str = ind_elems[5].substr(equal_pos+1,scolon_pos-equal_pos-1);
3439  lat_south = strtof(lat_south_str.c_str(),nullptr);
3440  }
3441  else
3442  throw InternalErr(__FILE__,__LINE__,"south bound coordinate is not right for TRMM level 3 products");
3443 
3444  }
3445  else
3446  throw InternalErr(__FILE__,__LINE__,"The TRMM grid SouthBoundingCoordinate doesn't exist.");
3447 
3448  if (0 == ind_elems[6].find("EastBoundingCoordinate")){
3449 
3450  size_t equal_pos = ind_elems[6].find_first_of('=');
3451  if(string::npos == equal_pos)
3452  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for TRMM level 3 products");
3453 
3454  size_t scolon_pos = ind_elems[6].find_first_of(';');
3455  if(string::npos == scolon_pos)
3456  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for TRMM level 3 products");
3457  if (equal_pos < scolon_pos){
3458  string lon_east_str = ind_elems[6].substr(equal_pos+1,scolon_pos-equal_pos-1);
3459  lon_east = strtof(lon_east_str.c_str(),nullptr);
3460  }
3461  else
3462  throw InternalErr(__FILE__,__LINE__,"south bound coordinate is not right for TRMM level 3 products");
3463 
3464  }
3465  else
3466  throw InternalErr(__FILE__,__LINE__,"The TRMM grid EastBoundingCoordinate doesn't exist.");
3467 
3468  if (0 == ind_elems[7].find("WestBoundingCoordinate")){
3469 
3470  size_t equal_pos = ind_elems[7].find_first_of('=');
3471  if(string::npos == equal_pos)
3472  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for TRMM level 3 products");
3473 
3474  size_t scolon_pos = ind_elems[7].find_first_of(';');
3475  if(string::npos == scolon_pos)
3476  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for TRMM level 3 products");
3477  if (equal_pos < scolon_pos){
3478  string lon_west_str = ind_elems[7].substr(equal_pos+1,scolon_pos-equal_pos-1);
3479  lon_west = strtof(lon_west_str.c_str(),nullptr);
3480  }
3481  else
3482  throw InternalErr(__FILE__,__LINE__,"south bound coordinate is not right for TRMM level 3 products");
3483 
3484  }
3485  else
3486  throw InternalErr(__FILE__,__LINE__,"The TRMM grid WestBoundingCoordinate doesn't exist.");
3487 
3488  if (false == check_reg_orig) {
3489  if (0 != ind_elems[8].find("Origin=SOUTHWEST"))
3490  throw InternalErr(__FILE__,__LINE__,"The TRMM grid origin is not SOUTHWEST.");
3491  }
3492 
3493  // Since we only treat the case when the Registration is center, so the size should be the
3494  // regular number size - 1.
3495  latsize =(int)((lat_north-lat_south)/lat_res);
3496  lonsize =(int)((lon_east-lon_west)/lon_res);
3497  lat_start = lat_south;
3498  lon_start = lon_west;
3499 }
3500 
3501 // Somehow the conversion of double to c++ string with sprintf causes the memory error in
3502 // the testing code. I used the following code to do the conversion. Most part of the code
3503 // in reverse, int_to_str and dtoa are adapted from geeksforgeeks.org
3504 
3505 // reverses a string 'str' of length 'len'
3506 void HDFCFUtil::rev_str(char *str, int len)
3507 {
3508  int i=0;
3509  int j=len-1;
3510  char temp = 0;
3511  while (i<j)
3512  {
3513  temp = str[i];
3514  str[i] = str[j];
3515  str[j] = temp;
3516  i++;
3517  j--;
3518  }
3519 }
3520 
3521 // Converts a given integer x to string str[]. d is the number
3522 // of digits required in output. If d is more than the number
3523 // of digits in x, then 0s are added at the beginning.
3524 int HDFCFUtil::int_to_str(int x, char str[], int d)
3525 {
3526  int i = 0;
3527  while (x)
3528  {
3529  str[i++] = (x%10) + '0';
3530  x = x/10;
3531  }
3532 
3533  // If number of digits required is more, then
3534  // add 0s at the beginning
3535  while (i < d)
3536  str[i++] = '0';
3537 
3538  rev_str(str, i);
3539  str[i] = '\0';
3540  return i;
3541 }
3542 
3543 // Converts a double floating point number to string.
3544 void HDFCFUtil::dtoa(double n, char *res, int afterpoint)
3545 {
3546  // Extract integer part
3547  auto ipart = (int)n;
3548 
3549  // Extract the double part
3550  double fpart = n - (double)ipart;
3551 
3552  // convert integer part to string
3553  int i = int_to_str(ipart, res, 0);
3554 
3555  // check for display option after point
3556  if (afterpoint != 0)
3557  {
3558  res[i] = '.'; // add dot
3559 
3560  // Get the value of fraction part upto given no.
3561  // of points after dot. The third parameter is needed
3562  // to handle cases like 233.007
3563  fpart = fpart * pow(10, afterpoint);
3564 
3565  // A round-error of 1 is found when casting to the integer for some numbers.
3566  // We need to correct it.
3567  auto final_fpart = (int)fpart;
3568  if(fpart -(int)fpart >0.5)
3569  final_fpart = (int)fpart +1;
3570  int_to_str(final_fpart, res + i + 1, afterpoint);
3571  }
3572 }
3573 
3574 
3575 string HDFCFUtil::get_double_str(double x,int total_digit,int after_point) {
3576 
3577  string str;
3578  if(x!=0) {
3579  vector<char> res;
3580  res.resize(total_digit);
3581  for(int i = 0; i<total_digit;i++)
3582  res[i] = '\0';
3583  if (x<0) {
3584  str.push_back('-');
3585  dtoa(-x,res.data(),after_point);
3586  for(int i = 0; i<total_digit;i++) {
3587  if(res[i] != '\0')
3588  str.push_back(res[i]);
3589  }
3590  }
3591  else {
3592  dtoa(x, res.data(), after_point);
3593  for(int i = 0; i<total_digit;i++) {
3594  if(res[i] != '\0')
3595  str.push_back(res[i]);
3596  }
3597  }
3598 
3599  }
3600  else
3601  str.push_back('0');
3602 
3603  return str;
3604 
3605 }
3606 
3607 string HDFCFUtil::get_int_str(int x) {
3608 
3609  string str;
3610  if(x > 0 && x <10)
3611  str.push_back((char)x+'0');
3612 
3613  else if (x >10 && x<100) {
3614  str.push_back((char)(x/10)+'0');
3615  str.push_back((char)(x%10)+'0');
3616  }
3617  else {
3618  int num_digit = 0;
3619  int abs_x = (x<0)?-x:x;
3620  while(abs_x/=10)
3621  num_digit++;
3622  if(x<=0)
3623  num_digit++;
3624  vector<char> buf;
3625  buf.resize(num_digit);
3626  snprintf(buf.data(),num_digit,"%d",x);
3627  str.assign(buf.data());
3628 
3629  }
3630 
3631  return str;
3632 
3633 }
3634 
3635 ssize_t HDFCFUtil::read_vector_from_file(int fd, vector<double> &val, size_t dtypesize) {
3636 
3637  ssize_t ret_val;
3638  ret_val = read(fd,val.data(),val.size()*dtypesize);
3639 
3640  return ret_val;
3641 }
3642 // Need to wrap a 'read buffer' from a pure file call here since read() is also a DAP function to read DAP data.
3643 ssize_t HDFCFUtil::read_buffer_from_file(int fd, void*buf, size_t total_read) {
3644 
3645  ssize_t ret_val;
3646  ret_val = read(fd,buf,total_read);
3647 
3648  return ret_val;
3649 }
3650 void HDFCFUtil::close_fileid(int32 sdfd, int32 fileid,int32 gridfd, int32 swathfd,bool pass_fileid) {
3651 
3652  if(false == pass_fileid) {
3653  if(sdfd != -1)
3654  SDend(sdfd);
3655  if(fileid != -1)
3656  Hclose(fileid);
3657 #ifdef USE_HDFEOS2_LIB
3658  if(gridfd != -1)
3659  GDclose(gridfd);
3660  if(swathfd != -1)
3661  SWclose(swathfd);
3662 
3663 #endif
3664  }
3665 
3666 }
3667 
3668 // Obtain the cache name. Since AIRS version 6 level 3 all share the same latitude and longitude,
3669 // we provide one set of latitude and longitude cache files for all AIRS level 3 version 6 products.
3670 string HDFCFUtil::obtain_cache_fname(const string & fprefix, const string &fname, const string &vname) {
3671 
3672  string cache_fname = fprefix;
3673  string base_file_name = basename(fname);
3674  // Identify this file from product name: AIRS, product level: .L3. or .L2. and version .v6.
3675  if((base_file_name.size() >12)
3676  && (base_file_name.compare(0,4,"AIRS") == 0)
3677  && (base_file_name.find(".L3.")!=string::npos)
3678  && (base_file_name.find(".v6.")!=string::npos)
3679  && ((vname == "Latitude") ||(vname == "Longitude")))
3680  cache_fname = cache_fname +"AIRS"+".L3."+".v6."+vname;
3681  else
3682  cache_fname = cache_fname + base_file_name +"_"+vname;
3683 
3684  return cache_fname;
3685 }
3686 
3687 // The current DDS cache is only for some products of which object information
3688 // can be retrieved via HDF4 SDS APIs. Currently only AIRS version 6 products are supported.
3689 size_t HDFCFUtil::obtain_dds_cache_size(const HDFSP::File*spf) {
3690 
3691  size_t total_bytes_written = 0;
3692  const vector<HDFSP::SDField *>& spsds = spf->getSD()->getFields();
3693  for (const auto &fd:spsds){
3694 
3695  // We will not handle when the SDS datatype is DFNT_CHAR now.
3696  if(DFNT_CHAR == fd->getType()) {
3697  total_bytes_written = 0;
3698  break;
3699  }
3700  else {
3701  // We need to store dimension names and variable names.
3702  int temp_rank = fd->getRank();
3703  for (const auto & dim:fd->getDimensions())
3704  total_bytes_written += (dim->getName()).size()+1;
3705 
3706  total_bytes_written +=(fd->getName()).size()+1;
3707 
3708  // Many a time variable new name is the same as variable name,
3709  // so we may just save one byte('\0') for such as a case.
3710  if((fd->getName()) != (fd->getNewName()))
3711  total_bytes_written +=(fd->getNewName()).size()+1;
3712  else
3713  total_bytes_written +=1;
3714 
3715  // We need to store 4 integers: rank, variable datatype, SDS reference number, fieldtype
3716  total_bytes_written +=(temp_rank+4)*sizeof(int);
3717  }
3718  }
3719 
3720  if(total_bytes_written != 0)
3721  total_bytes_written +=1;
3722 
3723  return total_bytes_written;
3724 
3725 }
3726 
3727 // Write the DDS of the special SDS-only HDF to a cache.
3728 void HDFCFUtil::write_sp_sds_dds_cache(const HDFSP::File* spf,FILE*dds_file,size_t total_bytes_dds_cache,const string &dds_filename) {
3729 
3730  BESDEBUG("h4"," Coming to write SDS DDS to a cache" << endl);
3731  char delimiter = '\0';
3732  char cend = '\n';
3733  size_t total_written_bytes_count = 0;
3734 
3735  // The buffer to hold information to write to a DDS cache file.
3736  vector<char>temp_buf;
3737  temp_buf.resize(total_bytes_dds_cache);
3738  char* temp_pointer = temp_buf.data();
3739 
3740  const vector<HDFSP::SDField *>& spsds = spf->getSD()->getFields();
3741 
3742  //Read SDS
3743  for(const auto &fd:spsds){
3744 
3745  // First, rank, fieldtype, SDS reference number, variable datatype, dimsize(rank)
3746  int sds_rank = fd->getRank();
3747  int sds_ref = fd->getFieldRef();
3748  int sds_dtype = fd->getType();
3749  int sds_ftype = fd->getFieldType();
3750 
3751  vector<int32>dimsizes;
3752  dimsizes.resize(sds_rank);
3753 
3754  // Size for each dimension
3755  const vector<HDFSP::Dimension*>& dims= fd->getDimensions();
3756  for(int i = 0; i < sds_rank; i++)
3757  dimsizes[i] = dims[i]->getSize();
3758 
3759  memcpy((void*)temp_pointer,(void*)&sds_rank,sizeof(int));
3760  temp_pointer +=sizeof(int);
3761  memcpy((void*)temp_pointer,(void*)&sds_ref,sizeof(int));
3762  temp_pointer +=sizeof(int);
3763  memcpy((void*)temp_pointer,(void*)&sds_dtype,sizeof(int));
3764  temp_pointer +=sizeof(int);
3765  memcpy((void*)temp_pointer,(void*)&sds_ftype,sizeof(int));
3766  temp_pointer +=sizeof(int);
3767 
3768  memcpy((void*)temp_pointer,(void*)dimsizes.data(),sds_rank*sizeof(int));
3769  temp_pointer +=sds_rank*sizeof(int);
3770 
3771  // total written bytes so far
3772  total_written_bytes_count +=(sds_rank+4)*sizeof(int);
3773 
3774  // Second, variable name,variable new name and SDS dim name(rank)
3775  // we need to a delimiter to distinguish each name.
3776  string temp_varname = fd->getName();
3777  vector<char>temp_val1(temp_varname.begin(),temp_varname.end());
3778  memcpy((void*)temp_pointer,(void*)temp_val1.data(),temp_varname.size());
3779  temp_pointer +=temp_varname.size();
3780  memcpy((void*)temp_pointer,&delimiter,1);
3781  temp_pointer +=1;
3782 
3783  total_written_bytes_count =total_written_bytes_count +(1+temp_varname.size());
3784 
3785  // To save the dds cache size and the speed to retrieve variable new name
3786  // we only save variable cf name when the variable cf name is not the
3787  // same as the variable name. When they are the same, a delimiter is saved for
3788  // variable cf name.
3789  if(fd->getName() == fd->getNewName()){
3790  memcpy((void*)temp_pointer,&delimiter,1);
3791  temp_pointer +=1;
3792  total_written_bytes_count +=1;
3793  }
3794  else {
3795  string temp_varnewname = fd->getNewName();
3796  vector<char>temp_val2(temp_varnewname.begin(),temp_varnewname.end());
3797  memcpy((void*)temp_pointer,(void*)temp_val2.data(),temp_varnewname.size());
3798  temp_pointer +=temp_varnewname.size();
3799  memcpy((void*)temp_pointer,&delimiter,1);
3800  temp_pointer +=1;
3801  total_written_bytes_count =total_written_bytes_count +(1+temp_varnewname.size());
3802  }
3803 
3804  // Now the name for each dimensions.
3805  for(int i = 0; i < sds_rank; i++) {
3806  string temp_dimname = dims[i]->getName();
3807  vector<char>temp_val3(temp_dimname.begin(),temp_dimname.end());
3808  memcpy((void*)temp_pointer,(void*)temp_val3.data(),temp_dimname.size());
3809  temp_pointer +=temp_dimname.size();
3810  memcpy((void*)temp_pointer,&delimiter,1);
3811  temp_pointer +=1;
3812  total_written_bytes_count =total_written_bytes_count +(1+temp_dimname.size());
3813  }
3814  }
3815 
3816  memcpy((void*)temp_pointer,&cend,1);
3817  total_written_bytes_count +=1;
3818 
3819  if(total_written_bytes_count != total_bytes_dds_cache) {
3820  stringstream s_total_written_count;
3821  s_total_written_count << total_written_bytes_count;
3822  stringstream s_total_bytes_dds_cache;
3823  s_total_bytes_dds_cache << total_bytes_dds_cache;
3824  string msg = "DDs cached file "+ dds_filename +" buffer size should be " + s_total_bytes_dds_cache.str() ;
3825  msg = msg + ". But the real size written in the buffer is " + s_total_written_count.str();
3826  throw InternalErr (__FILE__, __LINE__,msg);
3827  }
3828 
3829  size_t bytes_really_written = fwrite((const void*)temp_buf.data(),1,total_bytes_dds_cache,dds_file);
3830 
3831  if(bytes_really_written != total_bytes_dds_cache) {
3832  stringstream s_expected_bytes;
3833  s_expected_bytes << total_bytes_dds_cache;
3834  stringstream s_really_written_bytes;
3835  s_really_written_bytes << bytes_really_written;
3836  string msg = "DDs cached file "+ dds_filename +" size should be " + s_expected_bytes.str() ;
3837  msg += ". But the real size written to the file is " + s_really_written_bytes.str();
3838  throw InternalErr (__FILE__, __LINE__,msg);
3839  }
3840 
3841 }
3842 
3843 // Read DDS of a special SDS-only HDF file from a cache.
3844 void HDFCFUtil::read_sp_sds_dds_cache(FILE* dds_file,libdap::DDS * dds_ptr,
3845  const std::string &cache_filename, const std::string &hdf4_filename) {
3846 
3847  BESDEBUG("h4"," Coming to read SDS DDS from a cache" << endl);
3848 
3849  // Check the cache file size.
3850  struct stat sb;
3851  if(stat(cache_filename.c_str(),&sb)!=0) {
3852  string err_mesg="The DDS cache file " + cache_filename;
3853  err_mesg = err_mesg + " doesn't exist. ";
3854  throw InternalErr(__FILE__,__LINE__,err_mesg);
3855  }
3856 
3857  auto bytes_expected_read = (size_t)sb.st_size;
3858 
3859  // Allocate the buffer size based on the file size.
3860  vector<char> temp_buf;
3861  temp_buf.resize(bytes_expected_read);
3862  size_t bytes_really_read = fread((void*)temp_buf.data(),1,bytes_expected_read,dds_file);
3863 
3864  // Now bytes_really_read should be the same as bytes_expected_read if the element size is 1.
3865  if(bytes_really_read != bytes_expected_read) {
3866  stringstream s_bytes_really_read;
3867  s_bytes_really_read << bytes_really_read ;
3868  stringstream s_bytes_expected_read;
3869  s_bytes_expected_read << bytes_expected_read;
3870  string msg = "The expected bytes to read from DDS cache file " + cache_filename +" is " + s_bytes_expected_read.str();
3871  msg = msg + ". But the real read size from the buffer is " + s_bytes_really_read.str();
3872  throw InternalErr (__FILE__, __LINE__,msg);
3873  }
3874  char* temp_pointer = temp_buf.data();
3875 
3876  char delimiter = '\0';
3877  // The end of the whole string.
3878  char cend = '\n';
3879  bool end_file_flag = false;
3880 
3881 
3882  do {
3883  int sds_rank = *((int *)(temp_pointer));
3884  temp_pointer = temp_pointer + sizeof(int);
3885 
3886  int sds_ref = *((int *)(temp_pointer));
3887  temp_pointer = temp_pointer + sizeof(int);
3888 
3889  int sds_dtype = *((int *)(temp_pointer));
3890  temp_pointer = temp_pointer + sizeof(int);
3891 
3892  int sds_ftype = *((int *)(temp_pointer));
3893  temp_pointer = temp_pointer + sizeof(int);
3894 
3895  vector<int32>dimsizes;
3896  if(sds_rank <=0)
3897  throw InternalErr (__FILE__, __LINE__,"SDS rank must be >0");
3898 
3899  dimsizes.resize(sds_rank);
3900  for (int i = 0; i <sds_rank;i++) {
3901  dimsizes[i] = *((int *)(temp_pointer));
3902  temp_pointer = temp_pointer + sizeof(int);
3903  }
3904 
3905  vector<string>dimnames;
3906  dimnames.resize(sds_rank);
3907  string varname;
3908  string varnewname;
3909  for (int i = 0; i <sds_rank+2;i++) {
3910  vector<char> temp_vchar;
3911  char temp_char = *temp_pointer;
3912 
3913  // Only apply when varnewname is stored as the delimiter.
3914  if(temp_char == delimiter)
3915  temp_vchar.push_back(temp_char);
3916  while(temp_char !=delimiter) {
3917  temp_vchar.push_back(temp_char);
3918  temp_pointer++;
3919  temp_char = *temp_pointer;
3920  //temp_char = *(++temp_pointer); work
3921  //temp_char = *(temp_pointer++); not working
3922  }
3923  string temp_string(temp_vchar.begin(),temp_vchar.end());
3924  if(i == 0)
3925  varname = temp_string;
3926  else if( i == 1)
3927  varnewname = temp_string;
3928  else
3929  dimnames[i-2] = temp_string;
3930  temp_pointer++;
3931  }
3932 
3933  // If varnewname is only the delimiter, varname and varnewname is the same.
3934  if(varnewname[0] == delimiter)
3935  varnewname = varname;
3936  // Assemble DDS.
3937  // 1. Create a basetype
3938  BaseType *bt = nullptr;
3939  switch(sds_dtype) {
3940 #define HANDLE_CASE(tid, type) \
3941  case tid: \
3942  bt = new (type)(varnewname,hdf4_filename); \
3943  break;
3944  HANDLE_CASE(DFNT_FLOAT32, HDFFloat32)
3945  HANDLE_CASE(DFNT_FLOAT64, HDFFloat64)
3946  HANDLE_CASE(DFNT_CHAR, HDFStr)
3947 #ifndef SIGNED_BYTE_TO_INT32
3948  HANDLE_CASE(DFNT_INT8, HDFByte)
3949 #else
3950  HANDLE_CASE(DFNT_INT8,HDFInt32)
3951 #endif
3952  HANDLE_CASE(DFNT_UINT8, HDFByte)
3953  HANDLE_CASE(DFNT_INT16, HDFInt16)
3954  HANDLE_CASE(DFNT_UINT16, HDFUInt16)
3955  HANDLE_CASE(DFNT_INT32, HDFInt32)
3956  HANDLE_CASE(DFNT_UINT32, HDFUInt32)
3957  HANDLE_CASE(DFNT_UCHAR8, HDFByte)
3958  default:
3959  throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
3960 #undef HANDLE_CASE
3961  }
3962 
3963  if(nullptr == bt)
3964  throw InternalErr(__FILE__,__LINE__,"Cannot create the basetype when creating DDS from a cache file.");
3965 
3966  SPType sptype = OTHERHDF;
3967 
3968  // sds_ftype indicates if this is a general data field or geolocation field.
3969  // 4 indicates this is a missing non-latlon geo-location fields.
3970  if(sds_ftype != 4){
3971  HDFSPArray_RealField *ar = nullptr;
3972  try {
3973  // pass sds id as 0 since the sds id will be retrieved from SDStart if necessary.
3974  ar = new HDFSPArray_RealField(
3975  sds_rank,
3976  hdf4_filename,
3977  0,
3978  sds_ref,
3979  sds_dtype,
3980  sptype,
3981  varname,
3982  dimsizes,
3983  varnewname,
3984  bt);
3985  }
3986  catch(...) {
3987  delete bt;
3988  throw InternalErr(__FILE__,__LINE__,"Unable to allocate the HDFSPArray_RealField instance.");
3989  }
3990 
3991  for(int i = 0; i <sds_rank; i++)
3992  ar->append_dim(dimsizes[i],dimnames[i]);
3993  dds_ptr->add_var(ar);
3994  delete bt;
3995  delete ar;
3996  }
3997  else {
3998  HDFSPArrayMissGeoField *ar = nullptr;
3999  if(sds_rank == 1) {
4000  try {
4001  ar = new HDFSPArrayMissGeoField(
4002  sds_rank,
4003  dimsizes[0],
4004  varnewname,
4005  bt);
4006  }
4007  catch(...) {
4008  delete bt;
4009  throw InternalErr(__FILE__,__LINE__,"Unable to allocate the HDFSPArray_RealField instance.");
4010  }
4011 
4012  ar->append_dim(dimsizes[0],dimnames[0]);
4013  dds_ptr->add_var(ar);
4014  delete bt;
4015  delete ar;
4016  }
4017  else
4018  throw InternalErr(__FILE__,__LINE__,"SDS rank must be 1 for the missing coordinate.");
4019  }
4020 
4021  if(*temp_pointer == cend)
4022  end_file_flag = true;
4023  if((temp_pointer - temp_buf.data()) > (int)bytes_expected_read) {
4024  string msg = cache_filename + " doesn't have the end-line character at the end. The file may be corrupted.";
4025  throw InternalErr (__FILE__, __LINE__,msg);
4026  }
4027  } while(false == end_file_flag);
4028 
4029  dds_ptr->set_dataset_name(basename(hdf4_filename));
4030 }
4031 
4032 
4033 #if 0
4034 void HDFCFUtil::close_fileid(int32 sdfd, int32 fileid,int32 gridfd, int32 swathfd) {
4035 
4036  if(sdfd != -1)
4037  SDend(sdfd);
4038  if(fileid != -1)
4039  Hclose(fileid);
4040  if(gridfd != -1)
4041  GDclose(gridfd);
4042  if(swathfd != -1)
4043  SWclose(swathfd);
4044 
4045 }
4046 
4047 void HDFCFUtil::reset_fileid(int& sdfd, int& fileid,int& gridfd, int& swathfd) {
4048 
4049  sdfd = -1;
4050  fileid = -1;
4051  gridfd = -1;
4052  swathfd = -1;
4053 
4054 }
4055 #endif
static std::string lowercase(const std::string &s)
Definition: BESUtil.cc:254
const std::vector< Attribute * > & getAttributes() const
Get the attributes of this field.
Definition: HDFSP.h:306
int32 getType() const
Get the data type of this field.
Definition: HDFSP.h:300
const std::string & getNewName() const
Get the CF name(special characters replaced by underscores) of this field.
Definition: HDFSP.h:288
bool Has_Dim_NoScale_Field() const
This file has a field that is a SDS dimension but no dimension scale.
Definition: HDFSP.h:748
SD * getSD() const
Public interface to Obtain SD.
Definition: HDFSP.h:763
const std::vector< VDATA * > & getVDATAs() const
Public interface to Obtain Vdata.
Definition: HDFSP.h:769
SPType getSPType() const
Obtain special HDF4 product type.
Definition: HDFSP.h:741
const std::string & getPath() const
Obtain the path of the file.
Definition: HDFSP.h:757
One instance of this class represents one SDS object.
Definition: HDFSP.h:337
This class retrieves all SDS objects and SD file attributes.
Definition: HDFSP.h:542
const std::vector< SDField * > & getFields() const
Redundant member function.
Definition: HDFSP.h:561
const std::vector< Attribute * > & getAttributes() const
Public interface to obtain the SD(file) attributes.
Definition: HDFSP.h:567
Definition: HDFStr.h:51
void get_value(const std::string &s, std::string &val, bool &found)
Retrieve the value of a given key, if set.
Definition: TheBESKeys.cc:340
static TheBESKeys * TheKeys()
Definition: TheBESKeys.cc:71
static short obtain_type_size(int32)
Obtain datatype size.
Definition: HDFCFUtil.cc:450
static bool insert_map(std::map< std::string, std::string > &m, std::string key, std::string val)
Definition: HDFCFUtil.cc:150
static std::string print_attr(int32, int, void *)
Print attribute values in string.
Definition: HDFCFUtil.cc:268
static std::string print_type(int32)
Print datatype in string.
Definition: HDFCFUtil.cc:389
static void gen_unique_name(std::string &str, std::set< std::string > &namelist, int &clash_index)
Obtain the unique name for the clashed names and save it to set namelist.
Definition: HDFCFUtil.cc:196
static void Handle_NameClashing(std::vector< std::string > &newobjnamelist)
General routines to handle name clashings.
Definition: HDFCFUtil.cc:260
static void correct_scale_offset_type(libdap::AttrTable *at)
Definition: HDFCFUtil.cc:616
static std::string get_CF_string(std::string s)
Change special characters to "_".
Definition: HDFCFUtil.cc:166
static void Split(const char *s, int len, char sep, std::vector< std::string > &names)
Definition: HDFCFUtil.cc:82
static std::string escattr(std::string s)
Definition: HDFCFUtil.cc:3309
static void correct_fvalue_type(libdap::AttrTable *at, int32 dtype)
Definition: HDFCFUtil.cc:549
static void close_fileid(int32 sdfd, int32 file_id, int32 gridfd, int32 swathfd, bool pass_fileid_key)
Definition: HDFCFUtil.cc:3650
Definition: HDFCFUtil.h:52