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//
29using namespace std;
30using 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
43bool
44HDFCFUtil::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
66static void
67split_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"}
81void
82HDFCFUtil::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.
101void
102HDFCFUtil::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"}
121void
122HDFCFUtil::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.
136void
137HDFCFUtil::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
149bool
150HDFCFUtil::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
165string
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.
195void
196HDFCFUtil::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)
217void
218HDFCFUtil::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
259void
260HDFCFUtil::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.
267string
268HDFCFUtil::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.
388string
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.
449short
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.
499template < typename T >
500void 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.
549void 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".
662bool 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
686int 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//
725bool 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.
795void 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.
897bool 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
936void 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
962void 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}
984void 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//
1155void 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.
1536void 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.
1572void 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.
1631void 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.
1705void 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.
1736void 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.
1798void 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.
1882void 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.
1958void 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.
2097void 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
2189void
2190HDFCFUtil::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
2761void 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..
2814void 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
3023void 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.
3137cleanFail:
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
3147void 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
3206void 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
3309string 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
3343void 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'
3506void 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.
3524int 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.
3544void 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
3575string 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
3607string 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
3635ssize_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.
3643ssize_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}
3650void 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.
3670string 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.
3689size_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.
3728void 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.
3844void 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
4034void 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
4047void 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
const std::vector< VDATA * > & getVDATAs() const
Public interface to Obtain Vdata.
Definition: HDFSP.h:769
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
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< Attribute * > & getAttributes() const
Public interface to obtain the SD(file) attributes.
Definition: HDFSP.h:567
const std::vector< SDField * > & getFields() const
Redundant member function.
Definition: HDFSP.h:561
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