bes Updated for version 3.20.13
HDFEOS2.cc
1
2// This file is part of the hdf4 data handler for the OPeNDAP data server.
3//
4// Authors: MuQun Yang <myang6@hdfgroup.org> Choonghwan Lee
5// Copyright (c) 2009 The HDF Group
7
8#ifdef USE_HDFEOS2_LIB
9
10//#include <BESDEBUG.h> // Should provide BESDebug info. later
11#include <sstream>
12#include <algorithm>
13#include <functional>
14#include <vector>
15#include <map>
16#include <set>
17#include <cfloat>
18#include <math.h>
19#include <sys/stat.h>
20
21#include "HDFCFUtil.h"
22#include "HDFEOS2.h"
23#include "HDF4RequestHandler.h"
24
25// Names to be searched by
26// get_geodim_x_name()
27// get_geodim_y_name()
28// get_latfield_name()
29// get_lonfield_name()
30// get_geogrid_name()
31
32// Possible XDim names
33const char *HDFEOS2::File::_geodim_x_names[] = {"XDim", "LonDim","nlon"};
34
35// Possible YDim names.
36const char *HDFEOS2::File::_geodim_y_names[] = {"YDim", "LatDim","nlat"};
37
38// Possible latitude names.
39const char *HDFEOS2::File::_latfield_names[] = {
40 "Latitude", "Lat","YDim", "LatCenter"
41};
42
43// Possible longitude names.
44const char *HDFEOS2::File::_lonfield_names[] = {
45 "Longitude", "Lon","XDim", "LonCenter"
46};
47
48// For some grid products, latitude and longitude are put under a special geogrid.
49// One possible name is "location".
50const char *HDFEOS2::File::_geogrid_names[] = {"location"};
51
52using namespace HDFEOS2;
53using namespace std;
54
55// The following is for generating exception messages.
56template<typename T, typename U, typename V, typename W, typename X>
57static void _throw5(const char *fname, int line, int numarg,
58 const T &a1, const U &a2, const V &a3, const W &a4,
59 const X &a5)
60{
61 ostringstream ss;
62 ss << fname << ":" << line << ":";
63 for (int i = 0; i < numarg; ++i) {
64 ss << " ";
65 switch (i) {
66 case 0: ss << a1; break;
67 case 1: ss << a2; break;
68 case 2: ss << a3; break;
69 case 3: ss << a4; break;
70 case 4: ss << a5; break;
71 default:
72 ss<<" Argument number is beyond 5 ";
73 }
74 }
75 throw Exception(ss.str());
76}
77
79// number of arguments.
81#define throw1(a1) _throw5(__FILE__, __LINE__, 1, a1, 0, 0, 0, 0)
82#define throw2(a1, a2) _throw5(__FILE__, __LINE__, 2, a1, a2, 0, 0, 0)
83#define throw3(a1, a2, a3) _throw5(__FILE__, __LINE__, 3, a1, a2, a3, 0, 0)
84#define throw4(a1, a2, a3, a4) _throw5(__FILE__, __LINE__, 4, a1, a2, a3, a4, 0)
85#define throw5(a1, a2, a3, a4, a5) _throw5(__FILE__, __LINE__, 5, a1, a2, a3, a4, a5)
86
87#define assert_throw0(e) do { if (!(e)) throw1("assertion failure"); } while (false)
88#define assert_range_throw0(e, ge, l) assert_throw0((ge) <= (e) && (e) < (l))
89
90// Convenient function used in destructors.
91struct delete_elem
92{
93 template<typename T> void operator()(T *ptr)
94 {
95 delete ptr;
96 }
97};
98
99// Destructor for class File.
100File::~File()
101{
102 if(gridfd !=-1) {
103 for (vector<GridDataset *>::const_iterator i = grids.begin();
104 i != grids.end(); ++i){
105 delete *i;
106 }
107 // Grid file IDs will be closed in HDF4RequestHandler.cc.
108 }
109
110 if(swathfd !=-1) {
111 for (vector<SwathDataset *>::const_iterator i = swaths.begin();
112 i != swaths.end(); ++i){
113 delete *i;
114 }
115
116 }
117
118 for (vector<PointDataset *>::const_iterator i = points.begin();
119 i != points.end(); ++i){
120 delete *i;
121 }
122
123}
124
126File * File::Read(const char *path, int32 mygridfd, int32 myswathfd) throw(Exception)
127{
128
129 File *file = new File(path);
130 if(file == nullptr)
131 throw1("Memory allocation for file class failed. ");
132
133 file->gridfd = mygridfd;
134 file->swathfd = myswathfd;
135
136#if 0
137 // Read information of all Grid objects in this file.
138 if ((file->gridfd = GDopen(const_cast<char *>(file->path.c_str()),
139 DFACC_READ)) == -1) {
140 delete file;
141 throw2("grid open", path);
142 }
143#endif
144
145 vector<string> gridlist;
146 if (!Utility::ReadNamelist(file->path.c_str(), GDinqgrid, gridlist)) {
147 delete file;
148 throw1("Grid ReadNamelist failed.");
149 }
150
151 try {
152 for (const auto &grid: gridlist)
153 file->grids.push_back(GridDataset::Read(file->gridfd, grid));
154 }
155 catch(...) {
156 delete file;
157 throw1("GridDataset Read failed");
158 }
159
160#if 0
161 // Read information of all Swath objects in this file
162 if ((file->swathfd = SWopen(const_cast<char *>(file->path.c_str()),
163 DFACC_READ)) == -1){
164 delete file;
165 throw2("swath open", path);
166 }
167#endif
168
169 vector<string> swathlist;
170 if (!Utility::ReadNamelist(file->path.c_str(), SWinqswath, swathlist)){
171 delete file;
172 throw1("Swath ReadNamelist failed.");
173 }
174
175 try {
176 for (const auto &swath:swathlist)
177 file->swaths.push_back(SwathDataset::Read(file->swathfd, swath));
178 }
179 catch(...) {
180 delete file;
181 throw1("SwathDataset Read failed.");
182 }
183
184
185 // We only obtain the name list of point objects but not don't provide
186 // other information of these objects.
187 // The client will only get the name list of point objects.
188 vector<string> pointlist;
189 if (!Utility::ReadNamelist(file->path.c_str(), PTinqpoint, pointlist)){
190 delete file;
191 throw1("Point ReadNamelist failed.");
192 }
193 //See if I can make coverity happy because it doesn't understand throw macro.
194 for (const auto&point: pointlist)
195 file->points.push_back(PointDataset::Read(-1, point));
196
197 // If this is not an HDF-EOS2 file, returns exception as false.
198 if (file->grids.empty() && file->swaths.empty()
199 && file->points.empty()) {
200 Exception e("Not an HDF-EOS2 file");
201 e.setFileType(false);
202 delete file;
203 throw e;
204 }
205 return file;
206}
207
208
209// A grid's X-dimension can have different names: XDim, LatDim, etc.
210// This function returns the name of X-dimension which is used in the given file.
211// For better performance, we check the first grid or swath only.
212string File::get_geodim_x_name()
213{
214 if(!_geodim_x_name.empty())
215 return _geodim_x_name;
216 _find_geodim_names();
217 return _geodim_x_name;
218}
219
220// A grid's Y-dimension can have different names: YDim, LonDim, etc.
221// This function returns the name of Y-dimension which is used in the given file.
222// For better performance, we check the first grid or swath only.
223string File::get_geodim_y_name()
224{
225 if(!_geodim_y_name.empty())
226 return _geodim_y_name;
227 _find_geodim_names();
228 return _geodim_y_name;
229}
230
231// In some cases, values of latitude and longitude are stored in data fields.
232// Since the latitude field and longitude field do not have unique names
233// (e.g., latitude field can be "latitude", "Lat", ...),
234// we need to retrieve the field name.
235// The following two functions does this job.
236// For better performance, we check the first grid or swath only.
237
238string File::get_latfield_name()
239{
240 if(!_latfield_name.empty())
241 return _latfield_name;
242 _find_latlonfield_names();
243 return _latfield_name;
244}
245
246string File::get_lonfield_name()
247{
248 if(!_lonfield_name.empty())
249 return _lonfield_name;
250 _find_latlonfield_names();
251 return _lonfield_name;
252}
253
254// In some cases, a dedicated grid is used to store the values of
255// latitude and longitude. The following function finds the name
256// of the geo grid.
257
258string File::get_geogrid_name()
259{
260 if(!_geogrid_name.empty())
261 return _geogrid_name;
262 _find_geogrid_name();
263 return _geogrid_name;
264}
265
266// Internal function used by
267// get_geodim_x_name and get_geodim_y_name functions.
268// This function is not intended to be used outside the
269// get_geodim_x_name and get_geodim_y_name functions.
270
271void File::_find_geodim_names()
272{
273 set<string> geodim_x_name_set;
274 for(size_t i = 0; i<sizeof(_geodim_x_names) / sizeof(const char *); i++)
275 geodim_x_name_set.insert(_geodim_x_names[i]);
276
277 set<string> geodim_y_name_set;
278 for(size_t i = 0; i<sizeof(_geodim_y_names) / sizeof(const char *); i++)
279 geodim_y_name_set.insert(_geodim_y_names[i]);
280
281#if 0
282 // The following code is only used for grid. It also causes the coverity unhappy
283 // for the code block for(size_t i=0; ;i++), so simplify it after this code block.
284 const size_t gs = grids.size();
285 const size_t ss = swaths.size();
286 for(size_t i=0; ;i++)
287 {
288 Dataset *dataset=nullptr;
289 if(i<gs)
290 dataset = static_cast<Dataset*>(grids[i]);
291 else if(i < gs + ss)
292 dataset = static_cast<Dataset*>(swaths[i-gs]);
293 else
294 break;
295
296 const vector<Dimension *>& dims = dataset->getDimensions();
297 for(vector<Dimension*>::const_iterator it = dims.begin();
298 it != dims.end(); ++it)
299 {
300 if(geodim_x_name_set.find((*it)->getName()) != geodim_x_name_set.end())
301 _geodim_x_name = (*it)->getName();
302 else if(geodim_y_name_set.find((*it)->getName()) != geodim_y_name_set.end())
303 _geodim_y_name = (*it)->getName();
304 }
305 // For performance, we're checking this for the first grid or swath
306 break;
307 }
308#endif
309
310 const size_t gs = grids.size();
311 // For performance, we're checking this for the first grid
312 if(gs >0)
313 {
314 const Dataset *dataset=nullptr;
315 dataset = static_cast<Dataset*>(grids[0]);
316
317 const vector<Dimension *>& dims = dataset->getDimensions();
318 for (const auto &dim:dims)
319 {
320 // Essentially this code will grab any dimension names that is
321 // NOT predefined "XDim","LonDim","nlon" for geodim_x_name;
322 // any dimension names that is NOT predefined "YDim","LatDim","nlat"
323 // for geodim_y_name. This is in theory not right. Given the
324 // fact that this works with the current HDF-EOS2 products and there
325 // will be no more HDF-EOS2 products. We will leave the code this way.
326 if(geodim_x_name_set.find(dim->getName()) != geodim_x_name_set.end())
327 _geodim_x_name = dim->getName();
328 else if(geodim_y_name_set.find(dim->getName()) != geodim_y_name_set.end())
329 _geodim_y_name = dim->getName();
330 }
331 }
332 if(_geodim_x_name.empty())
333 _geodim_x_name = _geodim_x_names[0];
334 if(_geodim_y_name.empty())
335 _geodim_y_name = _geodim_y_names[0];
336}
337
338// Internal function used by
339// get_latfield_name and get_lonfield_name functions.
340// This function is not intended to be used outside
341// the get_latfield_name and get_lonfield_name functions.
342
343void File::_find_latlonfield_names()
344{
345 set<string> latfield_name_set;
346 for(size_t i = 0; i<sizeof(_latfield_names) / sizeof(const char *); i++)
347 latfield_name_set.insert(_latfield_names[i]);
348
349 set<string> lonfield_name_set;
350 for(size_t i = 0; i<sizeof(_lonfield_names) / sizeof(const char *); i++)
351 lonfield_name_set.insert(_lonfield_names[i]);
352
353 const size_t gs = grids.size();
354 const size_t ss = swaths.size();
355 // KY: converity structurally dead code i++ is never reached
356 // i++ is unreachable,so comment out this one
357#if 0
358 //for(size_t i=0; ;i++)
359#endif
360 for(size_t i=0;i<1 ;i++)
361 {
362 const Dataset *dataset = nullptr;
363 SwathDataset *sw = nullptr;
364 if(i<gs)
365 dataset = static_cast<Dataset*>(grids[i]);
366 else if(i < gs + ss)
367 {
368 sw = swaths[i-gs];
369 dataset = static_cast<Dataset*>(sw);
370 }
371 else
372 break;
373
374 const vector<Field *>& fields = dataset->getDataFields();
375 for (const auto &field:fields)
376 {
377 if(latfield_name_set.find(field->getName()) != latfield_name_set.end())
378 _latfield_name = field->getName();
379 else if(lonfield_name_set.find(field->getName()) != lonfield_name_set.end())
380 _lonfield_name = field->getName();
381 }
382
383 if(sw)
384 {
385 const vector<Field *>& geofields = dataset->getDataFields();
386 for(const auto &gfield:geofields)
387 {
388 if(latfield_name_set.find(gfield->getName()) != latfield_name_set.end())
389 _latfield_name = gfield->getName();
390 else if(lonfield_name_set.find(gfield->getName()) != lonfield_name_set.end())
391 _lonfield_name = gfield->getName();
392 }
393 }
394 // For performance, we're checking this for the first grid or swath
395 // comment out to fix coverity issues
396 //break;
397 }
398 if(_latfield_name.empty())
399 _latfield_name = _latfield_names[0];
400 if(_lonfield_name.empty())
401 _lonfield_name = _lonfield_names[0];
402
403}
404
405// Internal function used by
406// the get_geogrid_name function.
407// This function is not intended to be used outside the get_geogrid_name function.
408
409void File::_find_geogrid_name()
410{
411 set<string> geogrid_name_set;
412 for(size_t i = 0; i<sizeof(_geogrid_names) / sizeof(const char *); i++)
413 geogrid_name_set.insert(_geogrid_names[i]);
414
415 const size_t gs = grids.size();
416 const size_t ss = swaths.size();
417 for(size_t i=0; ;i++)
418 {
419 const Dataset *dataset = nullptr;
420 if(i<gs)
421 dataset = static_cast<Dataset*>(grids[i]);
422 else if(i < gs + ss)
423 dataset = static_cast<Dataset*>(swaths[i-gs]);
424 else
425 break;
426
427 if(geogrid_name_set.find(dataset->getName()) != geogrid_name_set.end())
428 _geogrid_name = dataset->getName();
429 }
430 if(_geogrid_name.empty())
431 _geogrid_name = "location";
432}
433
434// Check if we have the dedicated lat/lon grid.
435void File::check_onelatlon_grids() {
436
437 // 0. obtain "Latitude","Longitude" and "location" set.
438 string LATFIELDNAME = this->get_latfield_name();
439 string LONFIELDNAME = this->get_lonfield_name();
440
441 // Now only there is only one geo grid name "location", so don't call it know.
442 string GEOGRIDNAME = "location";
443
444 //only one lat/lon and it is under GEOGRIDNAME
445 int onellcount = 0;
446
447 // Check if lat/lon is found under other grids.
448 int morellcount = 0;
449
450 // Loop through all grids
451 for (const auto &grid:grids){
452
453 // Loop through all fields
454 for (const auto &field:grid->getDataFields()) {
455 if(grid->getName()==GEOGRIDNAME){
456 if(field->getName()==LATFIELDNAME){
457 onellcount++;
458 grid->latfield = field;
459 }
460 if(field->getName()==LONFIELDNAME){
461 onellcount++;
462 grid->lonfield = field;
463 }
464 if(onellcount == 2)
465 break;//Finish this grid
466 }
467 else {// Here we assume that lat and lon are always in pairs.
468 if((field->getName()==LATFIELDNAME)||(field->getName()==LONFIELDNAME)){
469 grid->ownllflag = true;
470 morellcount++;
471 break;
472 }
473 }
474 }
475 }
476
477 if(morellcount ==0 && onellcount ==2)
478 this->onelatlon = true;
479}
480
481// For one grid, need to handle the third-dimension(both existing and missing) coordinate variables
482void File::handle_one_grid_zdim(GridDataset* gdset) {
483
484 // Obtain "XDim","YDim"
485 string DIMXNAME = this->get_geodim_x_name();
486 string DIMYNAME = this->get_geodim_y_name();
487
488 bool missingfield_unlim_flag = false;
489 const Field *missingfield_unlim = nullptr;
490
491 // This is a big assumption, it may be wrong since not every 1-D field
492 // with the third dimension(name and size) is a coordinate
493 // variable. We have to watch the products we've supported. KY 2012-6-13
494
495 // Unique 1-D field's dimension name list.
496 set<string> tempdimlist;
497 pair<set<string>::iterator,bool> tempdimret;
498
499 for (const auto &field:gdset->getDataFields()) {
500 //We only need to search those 1-D fields
501
502 if (field->getRank()==1){
503
504 // DIMXNAME and DIMYNAME correspond to latitude and longitude.
505 // They should NOT be treated as dimension names missing fields. It will be handled differently.
506 if((field->getDimensions())[0]->getName()!=DIMXNAME && (field->getDimensions())[0]->getName()!=DIMYNAME){
507
508 tempdimret = tempdimlist.insert((field->getDimensions())[0]->getName());
509
510 // Kent: The following implementation may not be always right. This essentially is the flaw of the
511 // data product if a file encounters this case. Only unique 1-D third-dimension field should be provided.
512 // Only pick up the first 1-D field that the third-dimension
513 if(tempdimret.second == true) {
514
515 HDFCFUtil::insert_map(gdset->dimcvarlist, (field->getDimensions())[0]->getName(),
516 field->getName());
517 field->fieldtype = 3;
518 if(field->getName() == "Time")
519 field->fieldtype = 5;// IDV can handle 4-D fields when the 4th dim is Time.
520 }
521 }
522 }
523 }
524
525 // Add the missing Z-dimension field.
526 // Some dimension name's corresponding fields are missing,
527 // so add the missing Z-dimension fields based on the dimension names. When the real
528 // data is read, nature number 1,2,3,.... will be filled!
529 // NOTE: The latitude and longitude dim names are not handled yet.
530
531 // The above handling is also based on a big assumption. This is the best the
532 // handler can do without having the extra information outside the HDF-EOS2 file. KY 2012-6-12
533 // Loop through all dimensions of this grid.
534 for (const auto &gdim:gdset->getDimensions()) {
535
536 // Don't handle DIMXNAME and DIMYNAME yet.
537 if(gdim->getName()!=DIMXNAME && gdim->getName()!=DIMYNAME){
538
539 // This dimension needs a field
540 if((tempdimlist.find(gdim->getName())) == tempdimlist.end()){
541
542 // Need to create a new data field vector element with the name and dimension as above.
543 auto missingfield = new Field();
544 missingfield->name = gdim->getName();
545 missingfield->rank = 1;
546
547 //This is an HDF constant.the data type is always integer.
548 missingfield->type = DFNT_INT32;
549
550 // Dimension of the missing field
551 auto dim = new Dimension(gdim->getName(),gdim->getSize());
552
553 // only 1 dimension
554 missingfield->dims.push_back(dim);
555
556 if (0 == gdim->getSize()) {
557 missingfield_unlim_flag = true;
558 missingfield_unlim = missingfield;
559 }
560
561 // Provide information for the missing data, since we need to calculate the data, so
562 // the information is different than a normal field.
563#if 0
564 // The following code is for debugging purpose.
565 //int missingdatarank =1;
566 //int missingdatatypesize = 4;
567
568 //int missingdimsize[1]; //unused variable. SBL 2/7/20
569 //missingdimsize[0]= gdim->getSize(); //no purpose
570#endif
571
572
573 // added Z-dimension coordinate variable with nature number
574 missingfield->fieldtype = 4;
575
576 // input data is empty now. We need to review this approach in the future.
577 // The data will be retrieved in HDFEOS2ArrayMissGeoField.cc. KY 2013-06-14
578#if 0
579// LightVector<char>inputdata;
580// missingfield->data = nullptr;
581 //missingfield->data = new MissingFieldData(missingdatarank,missingdatatypesize,missingdimsize,inputdata);
582 // The data will be handled separately, we don't need to provide data.
583#endif
584 gdset->datafields.push_back(missingfield);
585 HDFCFUtil::insert_map(gdset->dimcvarlist, (missingfield->getDimensions())[0]->getName(),
586 missingfield->name);
587
588 }
589 }
590 }
591
592 //Correct the unlimited dimension size.
593 bool temp_missingfield_unlim_flag = missingfield_unlim_flag;
594 if(true == temp_missingfield_unlim_flag) {
595 for (unsigned int i =0; i<gdset->getDataFields().size(); i++) {
596
597 for (const auto &gdim:gdset->getDimensions()) {
598
599 if(gdim->getName() == (missingfield_unlim->getDimensions())[0]->getName()) {
600 if(gdim->getSize()!= 0) {
601 Dimension *dim = missingfield_unlim->getDimensions()[0];
602
603 // The unlimited dimension size is updated.
604 dim->dimsize = gdim->getSize();
605 missingfield_unlim_flag = false;
606 break;
607 }
608 }
609
610 }
611 if(false == missingfield_unlim_flag)
612 break;
613 }
614 }
615
616}
617
618// For one grid, need to handle lat/lon(both existing lat/lon and calculated lat/lon from EOS2 APIs)
619void File::handle_one_grid_latlon(GridDataset* gdset) throw(Exception)
620{
621
622 // Obtain "XDim","YDim","Latitude","Longitude"
623 string DIMXNAME = this->get_geodim_x_name();
624 string DIMYNAME = this->get_geodim_y_name();
625
626 string LATFIELDNAME = this->get_latfield_name();
627 string LONFIELDNAME = this->get_lonfield_name();
628
629
630 // This grid has its own latitude/longitude
631 if(gdset->ownllflag) {
632
633 // Searching the lat/lon field from the grid.
634 for (const auto &field:gdset->getDataFields()) {
635
636 if(field->getName() == LATFIELDNAME) {
637
638 // set the flag to tell if this is lat or lon.
639 // The unit will be different for lat and lon.
640 field->fieldtype = 1;
641
642 // Latitude rank should not be greater than 2.
643 // Here I don't check if the rank of latitude is the same as the longitude.
644 // Hopefully it never happens for HDF-EOS2 cases.
645 // We are still investigating if Java clients work
646 // when the rank of latitude and longitude is greater than 2.
647 if(field->getRank() > 2)
648 throw3("The rank of latitude is greater than 2",
649 gdset->getName(),field->getName());
650
651 if(field->getRank() != 1) {
652
653 // Obtain the major dim. For most cases, it is YDim Major.
654 // But still need to watch out the rare cases.
655 field->ydimmajor = gdset->getCalculated().isYDimMajor();
656
657 // If the 2-D lat/lon can be condensed to 1-D.
658 // For current HDF-EOS2 files, only GEO and CEA can be condensed.
659 // To gain performance,
660 // we don't check the real latitude values.
661 int32 projectioncode = gdset->getProjection().getCode();
662 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
663 field->condenseddim = true;
664 }
665
666 // Now we want to handle the dim and the var lists.
667 // If the lat/lon can be condensed to 1-D array,
668 // COARD convention needs to be followed.
669 // Since COARD requires that the dimension name of lat/lon is the same as the field name of lat/lon,
670 // we still need to handle this case in the later step(in function handle_grid_coards).
671 // Regardless of dimension major, always lat->YDim, lon->XDim;
672 // We don't need to adjust the dimension rank.
673 for (const auto &dim:field->getDimensions()) {
674 if (dim->getName() == DIMYNAME)
675 HDFCFUtil::insert_map(gdset->dimcvarlist, dim->getName(), field->getName());
676 }
677
678#if 0
679 // Don't know why this if/else block is necessary. It runs exactly the same code.
680 // Perhaps it tries to fix something and forgets to clean up. Comment out for a while. KY 2022-06-16
681 // If the 2-D array can be condensed to a 1-D array.
682 if(field->condenseddim) {
683
684 // Regardless of dimension major, always lat->YDim, lon->XDim;
685 // We don't need to adjust the dimension rank.
686 for (vector<Dimension *>::const_iterator k =
687 field->getDimensions().begin(); k!= field->getDimensions().end();++k){
688 if((*k)->getName() == DIMYNAME) {
689 HDFCFUtil::insert_map(gdset->dimcvarlist, dim->getName(), field->getName());
690 }
691 }
692 }
693
694 // 2-D lat/lon case. Since dimension order doesn't matter, so we always assume lon->XDim, lat->YDim.
695 else {
696 for (vector<Dimension *>::const_iterator k =
697 field->getDimensions().begin(); k!= field->getDimensions().end();++k){
698 if((*k)->getName() == DIMYNAME) {
699 HDFCFUtil::insert_map(gdset->dimcvarlist, (*k)->getName(), field->getName());
700 }
701 }
702 }
703#endif
704 }
705 // This is the 1-D case, just inserting the dimension, field pair.
706 else {
707 HDFCFUtil::insert_map(gdset->dimcvarlist, ((field->getDimensions())[0])->getName(),
708 field->getName());
709 }
710 }
711 else if (field->getName() == LONFIELDNAME) {
712
713 // set the flag to tell if this is lat or lon. The unit will be different for lat and lon.
714 field->fieldtype = 2;
715
716 // longitude rank should not be greater than 2.
717 // Here I don't check if the rank of latitude and longitude is the same.
718 // Hopefully it never happens for HDF-EOS2 cases.
719 // We are still investigating if Java clients work when the rank of latitude and longitude is greater than 2.
720 if(field->getRank() >2)
721 throw3("The rank of Longitude is greater than 2",gdset->getName(),field->getName());
722
723 if(field->getRank() != 1) {
724
725 // Obtain the major dim. For most cases, it is YDim Major. But still need to check for rare cases.
726 field->ydimmajor = gdset->getCalculated().isYDimMajor();
727
728 // If the 2-D lat/lon can be condensed to 1-D.
729 //For current HDF-EOS2 files, only GEO and CEA can be condensed. To gain performance,
730 // we don't check with real values.
731 int32 projectioncode = gdset->getProjection().getCode();
732 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
733 field->condenseddim = true;
734 }
735
736 // Now we want to handle the dim and the var lists.
737 // If the lat/lon can be condensed to 1-D array, COARD convention needs to be followed.
738 // Since COARD requires that the dimension name of lat/lon is the same as the field name of lat/lon,
739 // we still need to handle this case at last.
740 // When the field can be condensed, regardless of dimension major, the EOS convention is always lat->YDim, lon->XDim;
741 // We don't need to adjust the dimension rank.
742 // For 2-D lat/lon case: since dimension order doesn't matter, so we always assume lon->XDim, lat->YDim.
743 for (const auto &dim:field->getDimensions()) {
744 if(dim->getName() == DIMXNAME)
745 HDFCFUtil::insert_map(gdset->dimcvarlist, dim->getName(), field->getName());
746 }
747
748#if 0
749 // Leave this block commented out for a while, may delete it in the future. KY 2022-06-16
750 // Can be condensed to 1-D array.
751 if(field->condenseddim) {
752
753 // Regardless of dimension major, the EOS convention is always lat->YDim, lon->XDim;
754 // We don't need to adjust the dimension rank.
755 for (vector<Dimension *>::const_iterator k =
756 field->getDimensions().begin(); k!= field->getDimensions().end();++k){
757 if((*k)->getName() == DIMXNAME) {
758 HDFCFUtil::insert_map(gdset->dimcvarlist, (*k)->getName(), field->getName());
759 }
760 }
761 }
762 // 2-D lat/lon case. Since dimension order doesn't matter, so we always assume lon->XDim, lat->YDim.
763 else {
764 for (vector<Dimension *>::const_iterator k =
765 field->getDimensions().begin(); k!= field->getDimensions().end();++k){
766 if((*k)->getName() == DIMXNAME) {
767 HDFCFUtil::insert_map(gdset->dimcvarlist, (*k)->getName(), field->getName());
768 }
769 }
770 }
771#endif
772 }
773 // This is the 1-D case, just inserting the dimension, field pair.
774 else {
775 HDFCFUtil::insert_map(gdset->dimcvarlist, ((field->getDimensions())[0])->getName(),
776 field->getName());
777 }
778 }
779 }
780 }
781 else { // this grid's lat/lon has to be calculated.
782
783 // Latitude and Longitude
784 auto latfield = new Field();
785 auto lonfield = new Field();
786
787 latfield->name = LATFIELDNAME;
788 lonfield->name = LONFIELDNAME;
789
790 // lat/lon is a 2-D array
791 latfield->rank = 2;
792 lonfield->rank = 2;
793
794 // The data type is always float64. DFNT_FLOAT64 is the equivalent float64 type.
795 latfield->type = DFNT_FLOAT64;
796 lonfield->type = DFNT_FLOAT64;
797
798 // Latitude's fieldtype is 1
799 latfield->fieldtype = 1;
800
801 // Longitude's fieldtype is 2
802 lonfield->fieldtype = 2;
803
804 // Check if YDim is the major order.
805 // Obtain the major dim. For most cases, it is YDim Major. But some cases may be not. Still need to check.
806 latfield->ydimmajor = gdset->getCalculated().isYDimMajor();
807 lonfield->ydimmajor = latfield->ydimmajor;
808
809 // Obtain XDim and YDim size.
810 int xdimsize = gdset->getInfo().getX();
811 int ydimsize = gdset->getInfo().getY();
812
813 // Add dimensions. If it is YDim major,the dimension name list is "YDim XDim", otherwise, it is "XDim YDim".
814 // For LAMAZ projection, Y dimension is always supposed to be major for calculating lat/lon,
815 // but for dimension order, it should be consistent with data fields. (LD -2012/01/16
816 bool dmajor=(gdset->getProjection().getCode()==GCTP_LAMAZ)? gdset->getCalculated().DetectFieldMajorDimension()
817 : latfield->ydimmajor;
818
819 if(dmajor) {
820 auto dimlaty = new Dimension(DIMYNAME,ydimsize);
821 latfield->dims.push_back(dimlaty);
822 auto dimlony = new Dimension(DIMYNAME,ydimsize);
823 lonfield->dims.push_back(dimlony);
824 auto dimlatx = new Dimension(DIMXNAME,xdimsize);
825 latfield->dims.push_back(dimlatx);
826 auto dimlonx = new Dimension(DIMXNAME,xdimsize);
827 lonfield->dims.push_back(dimlonx);
828 }
829 else {
830 auto dimlatx = new Dimension(DIMXNAME,xdimsize);
831 latfield->dims.push_back(dimlatx);
832 auto dimlonx = new Dimension(DIMXNAME,xdimsize);
833 lonfield->dims.push_back(dimlonx);
834 auto dimlaty = new Dimension(DIMYNAME,ydimsize);
835 latfield->dims.push_back(dimlaty);
836 auto dimlony = new Dimension(DIMYNAME,ydimsize);
837 lonfield->dims.push_back(dimlony);
838 }
839
840 // Obtain info upleft and lower right for special longitude.
841#if 0
842 float64* upleft;
843 float64* lowright;
844 upleft = const_cast<float64 *>(gdset->getInfo().getUpLeft());
845 lowright = const_cast<float64 *>(gdset->getInfo().getLowRight());
846#endif
847
848 const float64* upleft = gdset->getInfo().getUpLeft();
849 const float64* lowright = gdset->getInfo().getLowRight();
850
851 // SOme special longitude is from 0 to 360.We need to check this case.
852 int32 projectioncode = gdset->getProjection().getCode();
853 if(((int)lowright[0]>180000000) && ((int)upleft[0]>-1)) {
854 // We can only handle geographic projection now.
855 // This is the only case we can handle.
856 if(projectioncode == GCTP_GEO) {// Will handle when data is read.
857 lonfield->speciallon = true;
858 // When HDF-EOS2 cache is involved, we have to also set the
859 // speciallon flag for the latfield since the cache file
860 // includes both lat and lon fields, and even the request
861 // is only to generate the lat field, the lon field also needs to
862 // be updated to write the proper cache. KY 2016-03-16
863 if(HDF4RequestHandler::get_enable_eosgeo_cachefile() == true)
864 latfield->speciallon = true;
865 }
866 }
867
868 // Some MODIS MCD files don't follow standard format for lat/lon (DDDMMMSSS);
869 // they simply represent lat/lon as -180.0000000 or -90.000000.
870 // HDF-EOS2 library won't give the correct value based on these value.
871 // These need to be remembered and resumed to the correct format when retrieving the data.
872 // Since so far we haven't found region of satellite files is within 0.1666 degree(1 minute)
873 // so, we divide the corner coordinate by 1000 and see if the integral part is 0.
874 // If it is 0, we know this file uses special lat/lon coordinate.
875
876 if(((int)(lowright[0]/1000)==0) &&((int)(upleft[0]/1000)==0)
877 && ((int)(upleft[1]/1000)==0) && ((int)(lowright[1]/1000)==0)) {
878 if(projectioncode == GCTP_GEO){
879 lonfield->specialformat = 1;
880 latfield->specialformat = 1;
881 }
882 }
883
884 // Some TRMM CERES Grid Data have "default" to be set for the corner coordinate,
885 // which they really mean for the whole globe(-180 - 180 lon and -90 - 90 lat).
886 // We will remember the information and change
887 // those values when we read the lat and lon.
888
889 if(((int)(lowright[0])==0) &&((int)(upleft[0])==0)
890 && ((int)(upleft[1])==0) && ((int)(lowright[1])==0)) {
891 if(projectioncode == GCTP_GEO){
892 lonfield->specialformat = 2;
893 latfield->specialformat = 2;
894 gdset->addfvalueattr = true;
895 }
896 }
897
898 //One MOD13C2 file doesn't provide projection code
899 // The upperleft and lowerright coordinates are all -1
900 // We have to calculate lat/lon by ourselves.
901 // Since it doesn't provide the project code, we double check their information
902 // and find that it covers the whole globe with 0.05 degree resolution.
903 // Lat. is from 90 to -90 and Lon is from -180 to 180.
904
905 if(((int)(lowright[0])==-1) &&((int)(upleft[0])==-1)
906 && ((int)(upleft[1])==-1) && ((int)(lowright[1])==-1)) {
907 lonfield->specialformat = 3;
908 latfield->specialformat = 3;
909 lonfield->condenseddim = true;
910 latfield->condenseddim = true;
911 }
912
913
914 // We need to handle SOM projection in a different way.
915 if(GCTP_SOM == projectioncode) {
916 lonfield->specialformat = 4;
917 latfield->specialformat = 4;
918 }
919
920 // Check if the 2-D lat/lon can be condensed to 1-D.
921 //For current HDF-EOS2 files, only GEO and CEA can be condensed. To gain performance,
922 // we just check the projection code, don't check with real values.
923 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
924 lonfield->condenseddim = true;
925 latfield->condenseddim = true;
926 }
927
928 // Add latitude and longitude fields to the field list.
929 gdset->datafields.push_back(latfield);
930 gdset->datafields.push_back(lonfield);
931
932 // Now we want to handle the dim and the var lists.
933 // If the lat/lon can be condensed to 1-D array, COARD convention needs to be followed.
934 // Since COARD requires that the dimension name of lat/lon is the same as the field name of lat/lon,
935 // we still need to handle this case later(function handle_grid_coards).
936
937 // There are two cases,
938 // 1) lat/lon can be condensed to 1-D array. lat to YDim, Lon to XDim, we don't need to adjust the rank.
939 // 2) 2-D lat./lon. The dimension order doesn't matter. So always assume lon to XDim, lat to YDim.
940 // So we can handle them with one loop.
941 for (const auto &dim:lonfield->getDimensions()) {
942 if(dim->getName() == DIMXNAME)
943 HDFCFUtil::insert_map(gdset->dimcvarlist, dim->getName(), lonfield->getName());
944
945 if(dim->getName() == DIMYNAME)
946 HDFCFUtil::insert_map(gdset->dimcvarlist, dim->getName(), latfield->getName());
947 }
948 }
949
950}
951
952// For the case of which all grids have one dedicated lat/lon grid,
953// this function shows how to handle lat/lon fields.
954void File::handle_onelatlon_grids() throw(Exception) {
955
956 // Obtain "XDim","YDim","Latitude","Longitude" and "location" set.
957 string DIMXNAME = this->get_geodim_x_name();
958 string DIMYNAME = this->get_geodim_y_name();
959 string LATFIELDNAME = this->get_latfield_name();
960 string LONFIELDNAME = this->get_lonfield_name();
961
962 // Now only there is only one geo grid name "location", so don't call it now.
963 // string GEOGRIDNAME = this->get_geogrid_name();
964 string GEOGRIDNAME = "location";
965
966 //Dimension name and the corresponding field name when only one lat/lon is used for all grids.
967 map<string,string>temponelatlondimcvarlist;
968
969 // First we need to obtain dimcvarlist for the grid that contains lat/lon.
970 for (const auto &grid:this->grids) {
971
972 // Set the horizontal dimension name "dimxname" and "dimyname"
973 // This will be used to detect the dimension major order.
974 grid->setDimxName(DIMXNAME);
975 grid->setDimyName(DIMYNAME);
976
977 // Handle lat/lon. Note that other grids need to point to this lat/lon.
978 if(grid->getName()==GEOGRIDNAME) {
979
980 // Figure out dimension order,2D or 1D for lat/lon
981 // if lat/lon field's pointed value is changed, the value of the lat/lon field is also changed.
982 // set the flag to tell if this is lat or lon. The unit will be different for lat and lon.
983 grid->lonfield->fieldtype = 2;
984 grid->latfield->fieldtype = 1;
985
986 // latitude and longitude rank must be equal and should not be greater than 2.
987 if(grid->lonfield->rank >2 || grid->latfield->rank >2)
988 throw2("Either the rank of lat or the lon is greater than 2",grid->getName());
989 if(grid->lonfield->rank !=grid->latfield->rank)
990 throw2("The rank of the latitude is not the same as the rank of the longitude",grid->getName());
991
992 // For 2-D lat/lon arrays
993 if(grid->lonfield->rank != 1) {
994
995 // Obtain the major dim. For most cases, it is YDim Major.
996 //But for some cases it is not. So still need to check.
997 grid->lonfield->ydimmajor = grid->getCalculated().isYDimMajor();
998 grid->latfield->ydimmajor = grid->lonfield->ydimmajor;
999
1000 // Check if the 2-D lat/lon can be condensed to 1-D.
1001 //For current HDF-EOS2 files, only GEO and CEA can be condensed. To gain performance,
1002 // we just check the projection code, don't check the real values.
1003 int32 projectioncode = grid->getProjection().getCode();
1004 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
1005 grid->lonfield->condenseddim = true;
1006 grid->latfield->condenseddim = true;
1007 }
1008
1009 // Now we want to handle the dim and the var lists.
1010 // If the lat/lon can be condensed to 1-D array, COARD convention needs to be followed.
1011 // Since COARD requires that the dimension name of lat/lon is the same as the field name of lat/lon,
1012 // we still need to handle this case later(function handle_grid_coards). Now we do the first step.
1013
1014 // There are two cases,
1015 // 1) lat/lon can be condensed to 1-D array. lat to YDim, Lon to XDim, we don't need to adjust the rank.
1016 // 2) 2-D lat./lon. The dimension order doesn't matter. So always assume lon to XDim, lat to YDim.
1017 // So we can handle them with one loop.
1018
1019 for (const auto &dim:grid->lonfield->getDimensions()) {
1020 if(dim->getName() == DIMXNAME) {
1021 HDFCFUtil::insert_map(grid->dimcvarlist, dim->getName(),
1022 grid->lonfield->getName());
1023 }
1024 if(dim->getName() == DIMYNAME) {
1025 HDFCFUtil::insert_map(grid->dimcvarlist, dim->getName(),
1026 grid->latfield->getName());
1027 }
1028 }
1029 }
1030 else { // This is the 1-D case, just inserting the dimension, field pair.
1031 HDFCFUtil::insert_map(grid->dimcvarlist, (grid->lonfield->getDimensions())[0]->getName(),
1032 grid->lonfield->getName());
1033 HDFCFUtil::insert_map(grid->dimcvarlist, (grid->latfield->getDimensions())[0]->getName(),
1034 grid->latfield->getName());
1035 }
1036 temponelatlondimcvarlist = grid->dimcvarlist;
1037 break;
1038
1039 }
1040
1041 }
1042
1043 // Now we need to assign the dim->cvar relation for lat/lon(xdim->lon,ydim->lat) to grids that don't contain lat/lon
1044 for (const auto &grid:this->grids) {
1045
1046 string templatlonname1;
1047 string templatlonname2;
1048
1049 if(grid->getName() != GEOGRIDNAME) {
1050
1051 map<string,string>::iterator tempmapit;
1052
1053 // Find DIMXNAME field
1054 tempmapit = temponelatlondimcvarlist.find(DIMXNAME);
1055 if(tempmapit != temponelatlondimcvarlist.end())
1056 templatlonname1= tempmapit->second;
1057 else
1058 throw2("cannot find the dimension field of XDim", grid->getName());
1059
1060 HDFCFUtil::insert_map(grid->dimcvarlist, DIMXNAME, templatlonname1);
1061
1062 // Find DIMYNAME field
1063 tempmapit = temponelatlondimcvarlist.find(DIMYNAME);
1064 if(tempmapit != temponelatlondimcvarlist.end())
1065 templatlonname2= tempmapit->second;
1066 else
1067 throw2("cannot find the dimension field of YDim", grid->getName());
1068 HDFCFUtil::insert_map(grid->dimcvarlist, DIMYNAME, templatlonname2);
1069 }
1070 }
1071
1072}
1073
1074// Handle the dimension name to coordinate variable map for grid.
1075void File::handle_grid_dim_cvar_maps() throw(Exception) {
1076
1077 // Obtain "XDim","YDim","Latitude","Longitude" and "location" set.
1078 string DIMXNAME = this->get_geodim_x_name();
1079
1080 string DIMYNAME = this->get_geodim_y_name();
1081
1082 string LATFIELDNAME = this->get_latfield_name();
1083
1084 string LONFIELDNAME = this->get_lonfield_name();
1085
1086
1087 // Now only there is only one geo grid name "location", so don't call it know.
1088 // string GEOGRIDNAME = this->get_geogrid_name();
1089 string GEOGRIDNAME = "location";
1090
1094
1095 // 1. Handle name clashings
1096 // 1.1 build up a temp. name list
1097 // Note here: we don't include grid and swath names(simply (*j)->name) due to the products we observe
1098 // Adding the grid/swath names makes the names artificially long. Will check user's feedback
1099 // and may change them later. KY 2012-6-25
1100 // The above assumption is purely for practical reason. Field names for all NASA multi-grid/swath products
1101 // (AIRS, AMSR-E, some MODIS, MISR) can all be distinguished regardless of grid/swath names. However,
1102 // this needs to be carefully watched out. KY 2013-07-08
1103 vector <string> tempfieldnamelist;
1104 for (const auto &grid:this->grids) {
1105 for (const auto &field:grid->getDataFields())
1106 tempfieldnamelist.push_back(HDFCFUtil::get_CF_string(field->name));
1107 }
1108 HDFCFUtil::Handle_NameClashing(tempfieldnamelist);
1109
1110 // 2. Create a map for dimension field name <original field name, corrected field name>
1111 // Also assure the uniqueness of all field names,save the new field names.
1112
1113 //the original dimension field name to the corrected dimension field name
1114 map<string,string>tempncvarnamelist;
1115 string tempcorrectedlatname, tempcorrectedlonname;
1116
1117 int total_fcounter = 0;
1118
1119 for (const auto &grid:this->grids) {
1120
1121 // Here we can't use getDataFields call since for no lat/lon fields
1122 // are created for one global lat/lon case. We have to use the dimcvarnamelist
1123 // map we just created.
1124 for (const auto &field:grid->getDataFields())
1125 {
1126 field->newname = tempfieldnamelist[total_fcounter];
1127 total_fcounter++;
1128
1129 // If this field is a dimension field, save the name/new name pair.
1130 if(field->fieldtype!=0) {
1131
1132 tempncvarnamelist.insert(make_pair(field->getName(), field->newname));
1133
1134 // For one latlon case, remember the corrected latitude and longitude field names.
1135 if((this->onelatlon)&&((grid->getName())==GEOGRIDNAME)) {
1136 if(field->getName()==LATFIELDNAME)
1137 tempcorrectedlatname = field->newname;
1138 if(field->getName()==LONFIELDNAME)
1139 tempcorrectedlonname = field->newname;
1140 }
1141 }
1142 }
1143
1144 grid->ncvarnamelist = tempncvarnamelist;
1145 tempncvarnamelist.clear();
1146 }
1147
1148 // For one lat/lon case, we have to add the lat/lon field name to other grids.
1149 // We know the original lat and lon names. So just retrieve the corrected lat/lon names from
1150 // the geo grid(GEOGRIDNAME).
1151 if(this->onelatlon) {
1152 for(const auto &grid:this->grids) {
1153 // Lat/lon names must be in this group.
1154 if(grid->getName()!=GEOGRIDNAME){
1155 HDFCFUtil::insert_map(grid->ncvarnamelist, LATFIELDNAME, tempcorrectedlatname);
1156 HDFCFUtil::insert_map(grid->ncvarnamelist, LONFIELDNAME, tempcorrectedlonname);
1157 }
1158 }
1159 }
1160
1161 // 3. Create a map for dimension name < original dimension name, corrected dimension name>
1162 map<string,string>tempndimnamelist;//the original dimension name to the corrected dimension name
1163
1165 vector <string>tempalldimnamelist;
1166 for (const auto &grid:this->grids) {
1167 for (map<string,string>::const_iterator j =
1168 grid->dimcvarlist.begin(); j!= grid->dimcvarlist.end();++j)
1169 tempalldimnamelist.push_back(HDFCFUtil::get_CF_string((*j).first));
1170 }
1171
1172 HDFCFUtil::Handle_NameClashing(tempalldimnamelist);
1173
1174 // Since DIMXNAME and DIMYNAME are not in the original dimension name list, we use the dimension name,field map
1175 // we just formed.
1176 int total_dcounter = 0;
1177 for (const auto &grid:this->grids) {
1178
1179 for (map<string,string>::const_iterator j =
1180 grid->dimcvarlist.begin(); j!= grid->dimcvarlist.end();++j){
1181
1182 // We have to handle DIMXNAME and DIMYNAME separately.
1183 if((DIMXNAME == (*j).first || DIMYNAME == (*j).first) && (true==(this->onelatlon)))
1184 HDFCFUtil::insert_map(tempndimnamelist, (*j).first,(*j).first);
1185 else
1186 HDFCFUtil::insert_map(tempndimnamelist, (*j).first, tempalldimnamelist[total_dcounter]);
1187 total_dcounter++;
1188 }
1189
1190 grid->ndimnamelist = tempndimnamelist;
1191 tempndimnamelist.clear();
1192 }
1193}
1194
1195// Follow COARDS for grids.
1196void File::handle_grid_coards() throw(Exception) {
1197
1198 // Obtain "XDim","YDim","Latitude","Longitude" and "location" set.
1199 string DIMXNAME = this->get_geodim_x_name();
1200 string DIMYNAME = this->get_geodim_y_name();
1201 string LATFIELDNAME = this->get_latfield_name();
1202 string LONFIELDNAME = this->get_lonfield_name();
1203
1204 // Now only there is only one geo grid name "location", so don't call it know.
1205 // string GEOGRIDNAME = this->get_geogrid_name();
1206 string GEOGRIDNAME = "location";
1207
1208
1209 // Revisit the lat/lon fields to check if 1-D COARD convention needs to be followed.
1210 vector<Dimension*> correcteddims;
1211 string tempcorrecteddimname;
1212
1213 // grid name to the corrected latitude field name
1214 map<string,string> tempnewxdimnamelist;
1215
1216 // grid name to the corrected longitude field name
1217 map<string,string> tempnewydimnamelist;
1218
1219 // temporary dimension pointer
1220 Dimension *correcteddim;
1221
1222 for (const auto &grid:this->grids) {
1223 for (const auto &field:grid->getDataFields()) {
1224
1225 // Now handling COARD cases, since latitude/longitude can be either 1-D or 2-D array.
1226 // So we need to correct both cases.
1227 // 2-D lat to 1-D COARD lat
1228 if(field->getName()==LATFIELDNAME && field->getRank()==2 &&field->condenseddim) {
1229
1230 string templatdimname;
1231 map<string,string>::iterator tempmapit;
1232
1233 // Find the new name of LATFIELDNAME
1234 tempmapit = grid->ncvarnamelist.find(LATFIELDNAME);
1235 if(tempmapit != grid->ncvarnamelist.end())
1236 templatdimname= tempmapit->second;
1237 else
1238 throw2("cannot find the corrected field of Latitude", grid->getName());
1239
1240 for (const auto &dim:field->getDimensions()) {
1241
1242 // Since hhis is the latitude, we create the corrected dimension with the corrected latitude field name
1243 // latitude[YDIM]->latitude[latitude]
1244 if(dim->getName()==DIMYNAME) {
1245 correcteddim = new Dimension(templatdimname,dim->getSize());
1246 correcteddims.push_back(correcteddim);
1247 field->setCorrectedDimensions(correcteddims);
1248 HDFCFUtil::insert_map(tempnewydimnamelist, grid->getName(), templatdimname);
1249 }
1250 }
1251 field->iscoard = true;
1252 grid->iscoard = true;
1253 if (this->onelatlon)
1254 this->iscoard = true;
1255
1256 // Clear the local temporary vector
1257 correcteddims.clear();
1258 }
1259
1260 // 2-D lon to 1-D COARD lon
1261 else if(field->getName()==LONFIELDNAME && field->getRank()==2 &&field->condenseddim){
1262
1263 string templondimname;
1264 map<string,string>::iterator tempmapit;
1265
1266 // Find the new name of LONFIELDNAME
1267 tempmapit = grid->ncvarnamelist.find(LONFIELDNAME);
1268 if(tempmapit != grid->ncvarnamelist.end())
1269 templondimname= tempmapit->second;
1270 else
1271 throw2("cannot find the corrected field of Longitude", grid->getName());
1272
1273 for (const auto &dim:field->getDimensions()) {
1274
1275 // Since this is the longitude, we create the corrected dimension with the corrected longitude field name
1276 // longitude[XDIM]->longitude[longitude]
1277 if(dim->getName()==DIMXNAME) {
1278 correcteddim = new Dimension(templondimname,dim->getSize());
1279 correcteddims.push_back(correcteddim);
1280 field->setCorrectedDimensions(correcteddims);
1281 HDFCFUtil::insert_map(tempnewxdimnamelist, grid->getName(), templondimname);
1282 }
1283 }
1284
1285 field->iscoard = true;
1286 grid->iscoard = true;
1287 if(this->onelatlon)
1288 this->iscoard = true;
1289 correcteddims.clear();
1290 }
1291 // 1-D lon to 1-D COARD lon
1292 // (this code can be combined with the 2-D lon to 1-D lon case, should handle this later, KY 2013-07-10).
1293 else if((field->getRank()==1) &&(field->getName()==LONFIELDNAME) ) {
1294
1295 string templondimname;
1296 map<string,string>::iterator tempmapit;
1297
1298 // Find the new name of LONFIELDNAME
1299 tempmapit = grid->ncvarnamelist.find(LONFIELDNAME);
1300 if(tempmapit != grid->ncvarnamelist.end())
1301 templondimname= tempmapit->second;
1302 else
1303 throw2("cannot find the corrected field of Longitude", grid->getName());
1304
1305 correcteddim = new Dimension(templondimname,(field->getDimensions())[0]->getSize());
1306 correcteddims.push_back(correcteddim);
1307 field->setCorrectedDimensions(correcteddims);
1308 field->iscoard = true;
1309 grid->iscoard = true;
1310 if(this->onelatlon)
1311 this->iscoard = true;
1312 correcteddims.clear();
1313
1314 if(((field->getDimensions())[0]->getName()!=DIMXNAME)
1315 &&(((field->getDimensions())[0]->getName())!=DIMYNAME)){
1316 throw3("the dimension name of longitude should not be ",
1317 (field->getDimensions())[0]->getName(),grid->getName());
1318 }
1319 if(((field->getDimensions())[0]->getName())==DIMXNAME) {
1320 HDFCFUtil::insert_map(tempnewxdimnamelist, grid->getName(), templondimname);
1321 }
1322 else {
1323 HDFCFUtil::insert_map(tempnewydimnamelist, grid->getName(), templondimname);
1324 }
1325 }
1326 // 1-D lat to 1-D COARD lat
1327 // (this case can be combined with the 2-D lat to 1-D lat case, should handle this later. KY 2013-7-10).
1328 else if((field->getRank()==1) &&(field->getName()==LATFIELDNAME) ) {
1329
1330 string templatdimname;
1331 map<string,string>::iterator tempmapit;
1332
1333 // Find the new name of LATFIELDNAME
1334 tempmapit = grid->ncvarnamelist.find(LATFIELDNAME);
1335 if(tempmapit != grid->ncvarnamelist.end())
1336 templatdimname= tempmapit->second;
1337 else
1338 throw2("cannot find the corrected field of Latitude", grid->getName());
1339
1340 correcteddim = new Dimension(templatdimname,(field->getDimensions())[0]->getSize());
1341 correcteddims.push_back(correcteddim);
1342 field->setCorrectedDimensions(correcteddims);
1343
1344 field->iscoard = true;
1345 grid->iscoard = true;
1346 if(this->onelatlon)
1347 this->iscoard = true;
1348 correcteddims.clear();
1349
1350 if((((field->getDimensions())[0]->getName())!=DIMXNAME)
1351 &&(((field->getDimensions())[0]->getName())!=DIMYNAME))
1352 throw3("the dimension name of latitude should not be ",
1353 (field->getDimensions())[0]->getName(),grid->getName());
1354 if(((field->getDimensions())[0]->getName())==DIMXNAME){
1355 HDFCFUtil::insert_map(tempnewxdimnamelist, grid->getName(), templatdimname);
1356 }
1357 else {
1358 HDFCFUtil::insert_map(tempnewydimnamelist, grid->getName(), templatdimname);
1359 }
1360 }
1361 }
1362 }
1363
1364 // If COARDS follows, apply the new DIMXNAME and DIMYNAME name to the ndimnamelist
1365 // One lat/lon for all grids.
1366 if(true == this->onelatlon){
1367
1368 // COARDS is followed.
1369 if(true == this->iscoard){
1370
1371 // For this case, only one pair of corrected XDim and YDim for all grids.
1372 string tempcorrectedxdimname;
1373 string tempcorrectedydimname;
1374
1375 if((int)(tempnewxdimnamelist.size())!= 1)
1376 throw1("the corrected dimension name should have only one pair");
1377 if((int)(tempnewydimnamelist.size())!= 1)
1378 throw1("the corrected dimension name should have only one pair");
1379
1380 map<string,string>::iterator tempdimmapit = tempnewxdimnamelist.begin();
1381 tempcorrectedxdimname = tempdimmapit->second;
1382 tempdimmapit = tempnewydimnamelist.begin();
1383 tempcorrectedydimname = tempdimmapit->second;
1384
1385 for (const auto &grid:this->grids) {
1386
1387 // Find the DIMXNAME and DIMYNAME in the dimension name list.
1388 map<string,string>::iterator tempmapit;
1389 tempmapit = grid->ndimnamelist.find(DIMXNAME);
1390 if(tempmapit != grid->ndimnamelist.end()) {
1391 HDFCFUtil::insert_map(grid->ndimnamelist, DIMXNAME, tempcorrectedxdimname);
1392 }
1393 else
1394 throw2("cannot find the corrected dimension name", grid->getName());
1395 tempmapit = grid->ndimnamelist.find(DIMYNAME);
1396 if(tempmapit != grid->ndimnamelist.end()) {
1397 HDFCFUtil::insert_map(grid->ndimnamelist, DIMYNAME, tempcorrectedydimname);
1398 }
1399 else
1400 throw2("cannot find the corrected dimension name", grid->getName());
1401 }
1402 }
1403 }
1404 else {// We have to search each grid
1405 for (const auto &grid:this->grids) {
1406
1407 if(grid->iscoard){
1408
1409 string tempcorrectedxdimname;
1410 string tempcorrectedydimname;
1411
1412 // Find the DIMXNAME and DIMYNAME in the dimension name list.
1413 map<string,string>::iterator tempdimmapit;
1414 map<string,string>::iterator tempmapit;
1415 tempdimmapit = tempnewxdimnamelist.find(grid->getName());
1416 if(tempdimmapit != tempnewxdimnamelist.end())
1417 tempcorrectedxdimname = tempdimmapit->second;
1418 else
1419 throw2("cannot find the corrected COARD XDim dimension name", grid->getName());
1420 tempmapit = grid->ndimnamelist.find(DIMXNAME);
1421 if(tempmapit != grid->ndimnamelist.end()) {
1422 HDFCFUtil::insert_map(grid->ndimnamelist, DIMXNAME, tempcorrectedxdimname);
1423 }
1424 else
1425 throw2("cannot find the corrected dimension name", grid->getName());
1426
1427 tempdimmapit = tempnewydimnamelist.find(grid->getName());
1428 if(tempdimmapit != tempnewydimnamelist.end())
1429 tempcorrectedydimname = tempdimmapit->second;
1430 else
1431 throw2("cannot find the corrected COARD YDim dimension name", grid->getName());
1432
1433 tempmapit = grid->ndimnamelist.find(DIMYNAME);
1434 if(tempmapit != grid->ndimnamelist.end()) {
1435 HDFCFUtil::insert_map(grid->ndimnamelist, DIMYNAME, tempcorrectedydimname);
1436 }
1437 else
1438 throw2("cannot find the corrected dimension name", grid->getName());
1439 }
1440 }
1441 }
1442
1443
1444 // For 1-D lat/lon cases, Make the third (other than lat/lon coordinate variable) dimension to follow COARD conventions.
1445
1446 for (const auto &grid:this->grids){
1447 for (map<string,string>::const_iterator j =
1448 grid->dimcvarlist.begin(); j!= grid->dimcvarlist.end();++j){
1449
1450 // It seems that the condition for onelatlon case is if(this->iscoard) is true instead if
1451 // this->onelatlon is true.So change it. KY 2010-7-4
1452 if((this->iscoard||grid->iscoard) && (*j).first !=DIMXNAME && (*j).first !=DIMYNAME) {
1453 string tempnewdimname;
1454 map<string,string>::iterator tempmapit;
1455
1456 // Find the new field name of the corresponding dimennsion name
1457 tempmapit = grid->ncvarnamelist.find((*j).second);
1458 if(tempmapit != grid->ncvarnamelist.end())
1459 tempnewdimname= tempmapit->second;
1460 else
1461 throw3("cannot find the corrected field of ", (*j).second,grid->getName());
1462
1463 // Make the new field name to the correponding dimension name
1464 tempmapit =grid->ndimnamelist.find((*j).first);
1465 if(tempmapit != grid->ndimnamelist.end())
1466 HDFCFUtil::insert_map(grid->ndimnamelist, (*j).first, tempnewdimname);
1467 else
1468 throw3("cannot find the corrected dimension name of ", (*j).first,grid->getName());
1469
1470 }
1471 }
1472 }
1473}
1474
1475// Create the corrected dimension vector for each field when COARDS is not followed.
1476void File::update_grid_field_corrected_dims() throw(Exception) {
1477
1478 // Revisit the lat/lon fields to check if 1-D COARD convention needs to be followed.
1479 vector<Dimension*> correcteddims;
1480 string tempcorrecteddimname;
1481 // temporary dimension pointer
1482 Dimension *correcteddim;
1483
1484 for (const auto &grid:this->grids) {
1485
1486 for (const auto &field:grid->getDataFields()) {
1487
1488 // When the corrected dimension name of lat/lon has been updated,
1489 if (field->iscoard == false) {
1490
1491 // Just obtain the corrected dim names and save the corrected dimensions for each field.
1492 for (const auto &dim:field->getDimensions()){
1493
1494 map<string,string>::iterator tempmapit;
1495
1496 // Find the new name of this field
1497 tempmapit = grid->ndimnamelist.find(dim->getName());
1498 if(tempmapit != grid->ndimnamelist.end())
1499 tempcorrecteddimname= tempmapit->second;
1500 else
1501 throw4("cannot find the corrected dimension name", grid->getName(),field->getName(),dim->getName());
1502 correcteddim = new Dimension(tempcorrecteddimname,dim->getSize());
1503 correcteddims.push_back(correcteddim);
1504 }
1505 field->setCorrectedDimensions(correcteddims);
1506 correcteddims.clear();
1507 }
1508 }
1509 }
1510
1511}
1512
1513void File::handle_grid_cf_attrs() throw(Exception) {
1514
1515 // Create "coordinates" ,"units" attributes. The "units" attributes only apply to latitude and longitude.
1516 // This is the last round of looping through everything,
1517 // we will match dimension name list to the corresponding dimension field name
1518 // list for every field.
1519
1520 for (const auto &grid:this->grids) {
1521 for (const auto &field:grid->getDataFields()) {
1522
1523 // Real fields: adding coordinate attributesinate attributes
1524 if(field->fieldtype == 0) {
1525 string tempcoordinates="";
1526 string tempfieldname="";
1527 string tempcorrectedfieldname="";
1528 int tempcount = 0;
1529 for (const auto &dim:field->getDimensions()) {
1530
1531 // Handle coordinates attributes
1532 map<string,string>::iterator tempmapit;
1533 map<string,string>::iterator tempmapit2;
1534
1535 // Find the dimension field name
1536 tempmapit = (grid->dimcvarlist).find(dim->getName());
1537 if(tempmapit != (grid->dimcvarlist).end())
1538 tempfieldname = tempmapit->second;
1539 else
1540 throw4("cannot find the dimension field name",
1541 grid->getName(),field->getName(),dim->getName());
1542
1543 // Find the corrected dimension field name
1544 tempmapit2 = (grid->ncvarnamelist).find(tempfieldname);
1545 if(tempmapit2 != (grid->ncvarnamelist).end())
1546 tempcorrectedfieldname = tempmapit2->second;
1547 else
1548 throw4("cannot find the corrected dimension field name",
1549 grid->getName(),field->getName(),dim->getName());
1550
1551 if(tempcount == 0)
1552 tempcoordinates= tempcorrectedfieldname;
1553 else
1554 tempcoordinates = tempcoordinates +" "+tempcorrectedfieldname;
1555 tempcount++;
1556 }
1557 field->setCoordinates(tempcoordinates);
1558 }
1559
1560 // Add units for latitude and longitude
1561 if(field->fieldtype == 1) {// latitude,adding the "units" degrees_north.
1562 string tempunits = "degrees_north";
1563 field->setUnits(tempunits);
1564 }
1565 if(field->fieldtype == 2) { // longitude, adding the units degrees_east.
1566 string tempunits = "degrees_east";
1567 field->setUnits(tempunits);
1568 }
1569
1570 // Add units for Z-dimension, now it is always "level"
1571 // This also needs to be corrected since the Z-dimension may not always be "level".
1572 // KY 2012-6-13
1573 // We decide not to touch "units" when the Z-dimension is an existing field(fieldtype =3).
1574 if(field->fieldtype == 4) {
1575 string tempunits ="level";
1576 field->setUnits(tempunits);
1577 }
1578
1579 // The units of the time is not right. KY 2012-6-13(documented at jira HFRHANDLER-167)
1580 if(field->fieldtype == 5) {
1581 string tempunits ="days since 1900-01-01 00:00:00";
1582 field->setUnits(tempunits);
1583 }
1584
1585 // We meet a really special case for CERES TRMM data. We attribute it to the specialformat 2 case
1586 // since the corner coordinate is set to default in HDF-EOS2 structmetadata. We also find that there are
1587 // values such as 3.4028235E38 that is the maximum single precision floating point value. This value
1588 // is a fill value but the fillvalue attribute is not set. So we add the fillvalue attribute for this case.
1589 // We may find such cases for other products and will tackle them also.
1590 if (true == grid->addfvalueattr) {
1591 if(((field->getFillValue()).empty()) && (field->getType()==DFNT_FLOAT32 )) {
1592 float tempfillvalue = FLT_MAX; // Replaced HUGE with FLT_MAX. jhrg 12/3/20
1593 field->addFillValue(tempfillvalue);
1594 field->setAddedFillValue(true);
1595 }
1596 }
1597 }
1598 }
1599}
1600
1601// Special handling SOM(Space Oblique Mercator) projection files
1602void File::handle_grid_SOM_projection() throw(Exception) {
1603
1604 // since the latitude and longitude of the SOM projection are 3-D, so we need to handle this projection in a special way.
1605 // Based on our current understanding, the third dimension size is always 180.
1606 // If the size is not 180, the latitude and longitude will not be calculated correctly.
1607 // This is according to the MISR Lat/lon calculation document
1608 // at http://eosweb.larc.nasa.gov/PRODOCS/misr/DPS/DPS_v50_RevS.pdf
1609 // KY 2012-6-12
1610
1611 for (const auto &grid:this->grids) {
1612 if (GCTP_SOM == grid->getProjection().getCode()) {
1613
1614 // 0. Getting the SOM dimension for latitude and longitude.
1615
1616 // Obtain SOM's dimension name.
1617 string som_dimname;
1618 for (const auto &dim:grid->getDimensions()) {
1619
1620 // NBLOCK is from misrproj.h. It is the number of block that MISR team support for the SOM projection.
1621 if(NBLOCK == dim->getSize()) {
1622
1623 // To make sure we catch the right dimension, check the first three characters of the dim. name
1624 // It should be SOM
1625 if (dim->getName().compare(0,3,"SOM") == 0) {
1626 som_dimname = dim->getName();
1627 break;
1628 }
1629 }
1630 }
1631
1632 if(""== som_dimname)
1633 throw4("Wrong number of block: The number of block of MISR SOM Grid ",
1634 grid->getName()," is not ",NBLOCK);
1635
1636 map<string,string>::iterator tempmapit;
1637
1638 // Find the corrected (CF) dimension name
1639 string cor_som_dimname;
1640 tempmapit = grid->ndimnamelist.find(som_dimname);
1641 if(tempmapit != grid->ndimnamelist.end())
1642 cor_som_dimname = tempmapit->second;
1643 else
1644 throw2("cannot find the corrected dimension name for ", som_dimname);
1645
1646 // Find the corrected(CF) name of this field
1647 string cor_som_cvname;
1648
1649 // Here we cannot use getDataFields() since the returned elements cannot be modified. KY 2012-6-12
1650 // Here we cannot simply change the vector with for range loop since we need to remove an element. KY 2022-06-17
1651 for (vector<Field *>::iterator j = grid->datafields.begin();
1652 j != grid->datafields.end(); ) {
1653
1654 // Only 6-7 fields, so just loop through
1655 // 1. Set the SOM dimension for latitude and longitude
1656 if (1 == (*j)->fieldtype || 2 == (*j)->fieldtype) {
1657
1658 auto newdim = new Dimension(som_dimname,NBLOCK);
1659 auto newcor_dim = new Dimension(cor_som_dimname,NBLOCK);
1660 vector<Dimension *>::iterator it_d;
1661
1662 it_d = (*j)->dims.begin();
1663 (*j)->dims.insert(it_d,newdim);
1664
1665 it_d = (*j)->correcteddims.begin();
1666 (*j)->correcteddims.insert(it_d,newcor_dim);
1667
1668
1669 }
1670
1671 // 2. Remove the added coordinate variable for the SOM dimension
1672 // The added variable is a variable with the nature number
1673 // Why removing it? Since we cannot follow the general rule to create coordinate variables for MISR products.
1674 // The third-dimension belongs to lat/lon rather than a missing coordinate variable.
1675 if ( 4 == (*j)->fieldtype) {
1676 cor_som_cvname = (*j)->newname;
1677 delete (*j);
1678 j = grid->datafields.erase(j);
1679 }
1680 else {
1681 ++j;
1682 }
1683 }
1684
1685 // 3. Fix the "coordinates" attribute: remove the SOM CV name from the coordinate attribute.
1686 // Notice this is a little inefficient. Since we only have a few fields and non-SOM projection products
1687 // won't be affected, and more importantly, to keep the SOM projection handling in a central place,
1688 // I handle the adjustment of "coordinates" attribute here. KY 2012-6-12
1689
1690 // MISR data cannot be visualized by Panoply and IDV. So the coordinates attribute
1691 // created here reflects the coordinates of this variable more accurately. KY 2012-6-13
1692 for (const auto &field:grid->getDataFields()) {
1693
1694 if ( 0 == field->fieldtype) {
1695
1696 string temp_coordinates = field->coordinates;
1697
1698 size_t found;
1699 found = temp_coordinates.find(cor_som_cvname);
1700
1701 if (0 == found) {
1702 // Need also to remove the space after the SOM CV name.
1703 if (temp_coordinates.size() >cor_som_cvname.size())
1704 temp_coordinates.erase(found,cor_som_cvname.size()+1);
1705 else
1706 temp_coordinates.erase(found,cor_som_cvname.size());
1707 }
1708 else if (found != string::npos)
1709 temp_coordinates.erase(found-1,cor_som_cvname.size()+1);
1710 else
1711 throw4("cannot find the coordinate variable ",cor_som_cvname,
1712 "from ",temp_coordinates);
1713
1714 field->setCoordinates(temp_coordinates);
1715
1716 }
1717 }
1718 }
1719 }
1720}
1721
1722// Check if we need to handle dim. map and set handle_swath_dimmap if necessary.
1723// The input parameter is the number of swath.
1724void File::check_swath_dimmap(int numswath) throw(Exception) {
1725
1726 if(HDF4RequestHandler::get_disable_swath_dim_map() == true)
1727 return;
1728
1729 // Check if there are dimension maps and if the num of dim. maps is odd in this case.
1730 int tempnumdm = 0;
1731 int temp_num_map = 0;
1732 bool odd_num_map = false;
1733 for (const auto &swath:this->swaths) {
1734 temp_num_map = swath->get_num_map();
1735 tempnumdm += temp_num_map;
1736 if(temp_num_map%2!=0) {
1737 odd_num_map =true;
1738 break;
1739 }
1740 }
1741
1742 // We only handle even number of dimension maps like MODIS(2-D lat/lon)
1743 if(tempnumdm != 0 && odd_num_map == false)
1744 handle_swath_dimmap = true;
1745
1746 // MODATML2 and MYDATML2 in year 2010 include dimension maps. But the dimension map
1747 // is not used. Furthermore, they provide additional latitude/longtiude
1748 // for 10 KM under the data field. So we have to handle this differently.
1749 // MODATML2 in year 2000 version doesn't include dimension map, so we
1750 // have to consider both with dimension map and without dimension map cases.
1751 // The swath name is atml2.
1752
1753 bool fakedimmap = false;
1754
1755 if(numswath == 1) {// Start special atml2-like handling
1756
1757 if((this->swaths[0]->getName()).find("atml2")!=string::npos){
1758
1759 if(tempnumdm >0)
1760 fakedimmap = true;
1761 int templlflag = 0;
1762
1763 for (const auto &gfield:this->swaths[0]->getGeoFields()) {
1764 if(gfield->getName() == "Latitude" || gfield->getName() == "Longitude") {
1765 if (gfield->getType() == DFNT_UINT16 ||
1766 gfield->getType() == DFNT_INT16)
1767 gfield->type = DFNT_FLOAT32;
1768 templlflag ++;
1769 if(templlflag == 2)
1770 break;
1771 }
1772 }
1773
1774 templlflag = 0;
1775
1776 for (const auto &dfield:this->swaths[0]->getDataFields()) {
1777
1778 // We meet a very speical MODIS case.
1779 // The latitude and longitude types are int16.
1780 // The number are in thousand. The scale factor
1781 // attribute is 0.01. This attribute cannot be
1782 // retrieved by EOS2 APIs. So we have to hard code this.
1783 // We can only use the swath name to
1784 // identify this case.
1785 // The swath name is atml2. It has only one swath.
1786 // We have to change lat and lon to float type array;
1787 // since after applying the scale factor, the array becomes
1788 // float data.
1789 // KY-2010-7-12
1790
1791 if((dfield->getName()).find("Latitude") != string::npos){
1792
1793 if (dfield->getType() == DFNT_UINT16 ||
1794 dfield->getType() == DFNT_INT16)
1795 dfield->type = DFNT_FLOAT32;
1796
1797 dfield->fieldtype = 1;
1798
1799 // Also need to link the dimension to the coordinate variable list
1800 if(dfield->getRank() != 2)
1801 throw2("The lat/lon rank must be 2 for Java clients to work",
1802 dfield->getRank());
1803 HDFCFUtil::insert_map(this->swaths[0]->dimcvarlist,
1804 ((dfield->getDimensions())[0])->getName(),dfield->getName());
1805 templlflag ++;
1806 }
1807
1808 if((dfield->getName()).find("Longitude")!= string::npos) {
1809
1810 if(dfield->getType() == DFNT_UINT16 ||
1811 dfield->getType() == DFNT_INT16)
1812 dfield->type = DFNT_FLOAT32;
1813
1814 dfield->fieldtype = 2;
1815 if(dfield->getRank() != 2)
1816 throw2("The lat/lon rank must be 2 for Java clients to work",
1817 dfield->getRank());
1818 HDFCFUtil::insert_map(this->swaths[0]->dimcvarlist,
1819 ((dfield->getDimensions())[1])->getName(), dfield->getName());
1820 templlflag ++;
1821 }
1822
1823 if(templlflag == 2)
1824 break;
1825 }
1826 }
1827 }// End of special atml2 handling
1828
1829 // Although this file includes dimension maps, it doesn't use it at all. So set
1830 // handle_swath_dimmap to 0.
1831 if(true == fakedimmap)
1832 handle_swath_dimmap = false;
1833 return;
1834
1835}
1836
1837// If dim. map needs to be handled, we need to check if we fall into the case
1838// that backward compatibility of MODIS Level 1B etc. should be supported.
1839void File::check_swath_dimmap_bk_compat(int numswath){
1840
1841 if(true == handle_swath_dimmap) {
1842
1843 if(numswath == 1 && (((this->swaths)[0])->name== "MODIS_SWATH_Type_L1B"))
1844 backward_handle_swath_dimmap = true;
1845 else {
1846 // If the number of dimmaps is 2 for every swath
1847 // and latitude/longitude need to be interpolated,
1848 // this also falls back to the backward compatibility case.
1849 // GeoDim_in_vars needs to be checked first.
1850 bool all_2_dimmaps_no_geodim = true;
1851 for (const auto &swath:this->swaths) {
1852 if (swath->get_num_map() !=2 || swath->GeoDim_in_vars == true) {
1853 all_2_dimmaps_no_geodim = false;
1854 break;
1855 }
1856 }
1857 if (true == all_2_dimmaps_no_geodim)
1858 backward_handle_swath_dimmap = true;
1859 }
1860 }
1861 return;
1862}
1863
1864// Create the dimension name to coordinate variable name map for lat/lon.
1865void File::create_swath_latlon_dim_cvar_map() throw(Exception){
1866
1867 vector<Field*> ori_lats;
1868 vector<Field*> ori_lons;
1869 if(handle_swath_dimmap == true && backward_handle_swath_dimmap == false) {
1870
1871 // We need to check if "Latitude and Longitude" both exist in all swaths under GeoFields.
1872 // The latitude and longitude must be 2-D arrays.
1873 // This is the basic requirement to handle our defined multiple dimension map case.
1874 multi_dimmap = true;
1875
1876 for (const auto &swath:this->swaths) {
1877
1878 bool has_cf_lat = false;
1879 bool has_cf_lon = false;
1880
1881 for (const auto &gfield:swath->getGeoFields()) {
1882
1883 // Here we assume it is always lat[f0][f1] and lon [f0][f1].
1884 // lat[f0][f1] and lon[f1][f0] should not occur.
1885 // So far only "Latitude" and "Longitude" are used as standard names of lat and lon for swath.
1886 if(gfield->getName()=="Latitude" && gfield->getRank() == 2){
1887 has_cf_lat = true;
1888 ori_lats.push_back(gfield);
1889 }
1890 else if(gfield->getName()=="Longitude" && gfield->getRank() == 2){
1891 has_cf_lon = true;
1892 ori_lons.push_back(gfield);
1893 }
1894 if(has_cf_lat == true && has_cf_lon == true)
1895 break;
1896 }
1897 if(has_cf_lat == false || has_cf_lon == false) {
1898 multi_dimmap = false;
1899 break;
1900 }
1901 }
1902 }
1903
1904 // By our best knowledge so far, we know we come to a multiple dimension map case
1905 // that we can handle. We will create dim to coordinate variable map for lat and lon
1906 // with the following block and finish this function.
1907 if (true == multi_dimmap) {
1908
1909 int ll_count = 0;
1910 for (const auto &swath:this->swaths) {
1911 create_swath_latlon_dim_cvar_map_for_dimmap(swath,ori_lats[ll_count],ori_lons[ll_count]);
1912 ll_count++;
1913 }
1914 return;
1915
1916 }
1917
1918 // For the cases that multi_dimmap is not true, do the following:
1919 // 1. Prepare the right dimension name and the dimension field list for each swath.
1920 // The assumption is that within a swath, the dimension name is unique.
1921 // The dimension field name(even with the added Z-like field) is unique.
1922 // A map <dimension name, dimension field name> will be created.
1923 // The name clashing handling for multiple swaths will not be done in this step.
1924
1925 // 1.1 Obtain the dimension names corresponding to the latitude and longitude,
1926 // save them to the <dimname, dimfield> map.
1927
1928 // We found a special MODIS product: the Latitude and Longitude are put under the Data fields
1929 // rather than GeoLocation fields.
1930 // So we need to go to the "Data Fields" to grab the "Latitude and Longitude".
1931
1932 bool lat_in_geofields = false;
1933 bool lon_in_geofields = false;
1934
1935 for (const auto &swath:this->swaths) {
1936
1937 int tempgeocount = 0;
1938 for (const auto &gfield:swath->getGeoFields()) {
1939
1940 // Here we assume it is always lat[f0][f1] and lon [f0][f1]. No lat[f0][f1] and lon[f1][f0] occur.
1941 // So far only "Latitude" and "Longitude" are used as standard names of lat and lon for swath.
1942 if(gfield->getName()=="Latitude" ){
1943 if(gfield->getRank() > 2)
1944 throw2("Currently the lat/lon rank must be 1 or 2 for Java clients to work",
1945 gfield->getRank());
1946
1947 lat_in_geofields = true;
1948
1949 // Since under our assumption, lat/lon are always 2-D for a swath and
1950 // dimension order doesn't matter for Java clients,
1951 // so we always map Latitude to the first dimension and longitude to the second dimension.
1952 // Save this information in the coordinate variable name and field map.
1953 // For rank =1 case, we only handle the cross-section along the same
1954 // longitude line. So Latitude should be the dimension name.
1955 HDFCFUtil::insert_map(swath->dimcvarlist, ((gfield->getDimensions())[0])->getName(), "Latitude");
1956
1957 // Have dimension map, we want to remember the dimension and remove it from the list.
1958 if(handle_swath_dimmap == true) {
1959
1960 // We need to keep the backward compatibility when handling MODIS level 1B etc.
1961 if(true == backward_handle_swath_dimmap) {
1962
1963 // We have to loop through the dimension map
1964 for (const auto &dmap:swath->getDimensionMaps()) {
1965
1966 // This dimension name will be replaced by the mapped dimension name,
1967 // the mapped dimension name can be obtained from the getDataDimension() method.
1968 if((gfield->getDimensions()[0])->getName() == dmap->getGeoDimension()) {
1969 HDFCFUtil::insert_map(swath->dimcvarlist, dmap->getDataDimension(), "Latitude");
1970 break;
1971 }
1972 }
1973 }
1974 }
1975
1976 gfield->fieldtype = 1;
1977 tempgeocount ++;
1978 }
1979
1980 if(gfield->getName()=="Longitude"){
1981 if(gfield->getRank() > 2)
1982 throw2("Currently the lat/lon rank must be 1 or 2 for Java clients to work",
1983 gfield->getRank());
1984
1985 // Only lat-level cross-section(for Panoply)is supported
1986 // when longitude/latitude is 1-D, so ignore the longitude as the dimension field.
1987 lon_in_geofields = true;
1988 if(gfield->getRank() == 1) {
1989 tempgeocount++;
1990 continue;
1991 }
1992
1993 // Since under our assumption, lat/lon are almost always 2-D for
1994 // a swath and dimension order doesn't matter for Java clients,
1995 // we always map Latitude to the first dimension and longitude to the second dimension.
1996 // Save this information in the dimensiion name and coordinate variable map.
1997 HDFCFUtil::insert_map(swath->dimcvarlist,
1998 ((gfield->getDimensions())[1])->getName(), "Longitude");
1999 if(handle_swath_dimmap == true) {
2000 if(true == backward_handle_swath_dimmap) {
2001
2002 // We have to loop through the dimension map
2003 for (const auto &dmap:swath->getDimensionMaps()) {
2004
2005 // This dimension name will be replaced by the mapped dimension name,
2006 // This name can be obtained by getDataDimension() fuction of
2007 // dimension map class.
2008 if((gfield->getDimensions()[1])->getName() ==
2009 dmap->getGeoDimension()) {
2010 HDFCFUtil::insert_map(swath->dimcvarlist,
2011 dmap->getDataDimension(), "Longitude");
2012 break;
2013 }
2014 }
2015 }
2016 }
2017 gfield->fieldtype = 2;
2018 tempgeocount++;
2019 }
2020 if(tempgeocount == 2)
2021 break;
2022 }
2023 }// end of creating the <dimname,dimfield> map.
2024
2025 // If lat and lon are not together, throw an error.
2026 //if (lat_in_geofields ^ lon_in_geofields)
2027 if (lat_in_geofields!=lon_in_geofields)
2028 throw1("Latitude and longitude must be both under Geolocation fields or Data fields");
2029
2030 // Check if they are under data fields(The code may be combined with the above, see HFRHANDLER-166)
2031 if (!lat_in_geofields && !lon_in_geofields) {
2032
2033 bool lat_in_datafields = false;
2034 bool lon_in_datafields = false;
2035
2036 for (const auto &swath:this->swaths) {
2037
2038 int tempgeocount = 0;
2039 for (const auto &dfield:swath->getDataFields()) {
2040
2041 // Here we assume it is always lat[f0][f1] and lon [f0][f1].
2042 // No lat[f0][f1] and lon[f1][f0] occur.
2043 // So far only "Latitude" and "Longitude" are used as
2044 // standard names of Lat and lon for swath.
2045 if (dfield->getName()=="Latitude" ){
2046 if (dfield->getRank() > 2) {
2047 throw2("Currently the lat/lon rank must be 1 or 2 for Java clients to work",
2048 dfield->getRank());
2049 }
2050 lat_in_datafields = true;
2051
2052 // Since under our assumption, lat/lon are always 2-D
2053 // for a swath and dimension order doesn't matter for Java clients,
2054 // we always map Latitude the first dimension and longitude the second dimension.
2055 // Save this information in the coordinate variable name and field map.
2056 // For rank =1 case, we only handle the cross-section along the same longitude line.
2057 // So Latitude should be the dimension name.
2058 HDFCFUtil::insert_map(swath->dimcvarlist,
2059 ((dfield->getDimensions())[0])->getName(), "Latitude");
2060
2061 if (handle_swath_dimmap == true) {
2062 if (true == backward_handle_swath_dimmap) {
2063 // We have to loop through the dimension map
2064 for (const auto &dmap:swath->getDimensionMaps()) {
2065 // This dimension name will be replaced by the mapped dimension name,
2066 // the mapped dimension name can be obtained from the getDataDimension() method.
2067 if((dfield->getDimensions()[0])->getName() == dmap->getGeoDimension()) {
2068 HDFCFUtil::insert_map(swath->dimcvarlist, dmap->getDataDimension(), "Latitude");
2069 break;
2070 }
2071 }
2072 }
2073 }
2074 dfield->fieldtype = 1;
2075 tempgeocount ++;
2076 }
2077
2078 if(dfield->getName()=="Longitude"){
2079
2080 if(dfield->getRank() > 2) {
2081 throw2("Currently the lat/lon rank must be 1 or 2 for Java clients to work",
2082 dfield->getRank());
2083 }
2084
2085 // Only lat-level cross-section(for Panoply)is supported when
2086 // longitude/latitude is 1-D, so ignore the longitude as the dimension field.
2087 lon_in_datafields = true;
2088 if(dfield->getRank() == 1) {
2089 tempgeocount++;
2090 continue;
2091 }
2092
2093 // Since under our assumption,
2094 // lat/lon are almost always 2-D for a swath and dimension order doesn't matter for Java clients,
2095 // we always map Latitude the first dimension and longitude the second dimension.
2096 // Save this information in the dimensiion name and coordinate variable map.
2097 HDFCFUtil::insert_map(swath->dimcvarlist,
2098 ((dfield->getDimensions())[1])->getName(), "Longitude");
2099 if(handle_swath_dimmap == true) {
2100 if(true == backward_handle_swath_dimmap) {
2101 // We have to loop through the dimension map
2102 for (const auto &dmap:swath->getDimensionMaps()) {
2103 // This dimension name will be replaced by the mapped dimension name,
2104 // This name can be obtained by getDataDimension() fuction of dimension map class.
2105 if((dfield->getDimensions()[1])->getName() == dmap->getGeoDimension()) {
2106 HDFCFUtil::insert_map(swath->dimcvarlist,
2107 dmap->getDataDimension(), "Longitude");
2108 break;
2109 }
2110 }
2111 }
2112 }
2113 dfield->fieldtype = 2;
2114 tempgeocount++;
2115 }
2116 if(tempgeocount == 2)
2117 break;
2118 }
2119 }// end of creating the <dimname,dimfield> map.
2120
2121 // If lat and lon are not together, throw an error.
2122 //if (lat_in_datafields ^ lon_in_datafields)
2123 if (lat_in_datafields!=lon_in_datafields)
2124 throw1("Latitude and longitude must be both under Geolocation fields or Data fields");
2125
2126 }
2127
2128}
2129
2130// Create the dimension name to coordinate variable name map for coordinate variables that are not lat/lon.
2131// The input parameter is the number of dimension maps in this file.
2132void File:: create_swath_nonll_dim_cvar_map() throw(Exception)
2133{
2134 // Handle existing and missing fields
2135 for (const auto &swath:this->swaths) {
2136
2137 // Since we find multiple 1-D fields with the same dimension names for some Swath files(AIRS level 1B),
2138 // we currently always treat the third dimension field as a missing field, this may be corrected later.
2139 // Corrections for the above: MODIS data do include the unique existing Z fields, so we have to search
2140 // the existing Z field. KY 2010-8-11
2141 // Our correction is to search all 1-D fields with the same dimension name within one swath,
2142 // if only one field is found, we use this field as the third dimension.
2143 // 1.1 Add the missing Z-dimension field.
2144 // Some dimension name's corresponding fields are missing,
2145 // so add the missing Z-dimension fields based on the dimension name. When the real
2146 // data is read, nature number 1,2,3,.... will be filled!
2147 // NOTE: The latitude and longitude dim names are not handled yet.
2148
2149 // Build a unique 1-D dimension name list.
2150 // Now the list only includes dimension names of "latitude" and "longitude".
2151
2152 pair<set<string>::iterator,bool> tempdimret;
2153 for (map<string,string>::const_iterator j = swath->dimcvarlist.begin();
2154 j!= swath->dimcvarlist.end();++j){
2155 tempdimret = swath->nonmisscvdimlist.insert((*j).first);
2156 }
2157
2158 // Search the geofield group and see if there are any existing 1-D Z dimension data.
2159 // If 1-D field data with the same dimension name is found under GeoField,
2160 // we still search if that 1-D field is the dimension
2161 // field of a dimension name.
2162 for (const auto &gfield:swath->getGeoFields()) {
2163
2164 if(gfield->getRank()==1) {
2165 if(swath->nonmisscvdimlist.find(((gfield->getDimensions())[0])->getName()) == swath->nonmisscvdimlist.end()){
2166 tempdimret = swath->nonmisscvdimlist.insert(((gfield->getDimensions())[0])->getName());
2167 if(gfield->getName() =="Time")
2168 gfield->fieldtype = 5;// This is for IDV.
2169
2170 // This is for temporarily COARD fix.
2171 // For 2-D lat/lon, the third dimension should NOT follow
2172 // COARD conventions. It will cause Panoply and IDV failed.
2173 // KY 2010-7-21
2174 // It turns out that we need to keep the original field name of the third dimension.
2175 // So assign the flag and save the original name.
2176 // KY 2010-9-9
2177#if 0
2178 if((((gfield->getDimensions())[0])->getName())==gfield->getName()){
2179 gfield->oriname = gfield->getName();
2180 // netCDF-Java fixes the problem, now goes back to COARDS.
2181 //gfield->name = gfield->getName() +"_d";
2182 gfield->specialcoard = true;
2183 }
2184#endif
2185 HDFCFUtil::insert_map(swath->dimcvarlist, ((gfield->getDimensions())[0])->getName(), gfield->getName());
2186 gfield->fieldtype = 3;
2187
2188 }
2189 }
2190 }
2191
2192 // We will also check the third dimension inside DataFields
2193 // This may cause potential problems for AIRS data
2194 // We will double CHECK KY 2010-6-26
2195 // So far the tests seem okay. KY 2010-8-11
2196 for (const auto &dfield:swath->getDataFields()) {
2197
2198 if(dfield->getRank()==1) {
2199 if(swath->nonmisscvdimlist.find(((dfield->getDimensions())[0])->getName()) == swath->nonmisscvdimlist.end()){
2200 tempdimret = swath->nonmisscvdimlist.insert(((dfield->getDimensions())[0])->getName());
2201 if(dfield->getName() =="Time")
2202 dfield->fieldtype = 5;// This is for IDV.
2203
2204 // This is for temporarily COARD fix.
2205 // For 2-D lat/lon, the third dimension should NOT follow
2206 // COARD conventions. It will cause Panoply and IDV failed.
2207 // KY 2010-7-21
2208#if 0
2209 if((((dfield->getDimensions())[0])->getName())==dfield->getName()){
2210 dfield->oriname = dfield->getName();
2211 //dfield->name = dfield->getName() +"_d";
2212 dfield->specialcoard = true;
2213 }
2214#endif
2215 HDFCFUtil::insert_map(swath->dimcvarlist, ((dfield->getDimensions())[0])->getName(), dfield->getName());
2216 dfield->fieldtype = 3;
2217
2218 }
2219 }
2220 }
2221
2222
2223 // S1.2.3 Handle the missing fields
2224 // Loop through all dimensions of this swath to search the missing fields
2225 //
2226 bool missingfield_unlim_flag = false;
2227 Field *missingfield_unlim = nullptr;
2228
2229 for (const auto &sdim:swath->getDimensions()) {
2230
2231 if((swath->nonmisscvdimlist.find(sdim->getName())) == swath->nonmisscvdimlist.end()){// This dimension needs a field
2232
2233 // Need to create a new data field vector element with the name and dimension as above.
2234 auto missingfield = new Field();
2235
2236 // This is for temporarily COARD fix.
2237 // For 2-D lat/lon, the third dimension should NOT follow
2238 // COARD conventions. It will cause Panoply and IDV failed.
2239 // Since Swath is always 2-D lat/lon, so we are okay here. Add a "_d" for each field name.
2240 // KY 2010-7-21
2241 // netCDF-Java now first follows COARDS, change back
2242 // missingfield->name = sdim->getName()+"_d";
2243 Dimension *dim;
2244 // When we can handle multiple dimension maps and the
2245 // number of swath is >1, we add the swath name as suffix to
2246 // avoid the name clashing.
2247 if(true == multi_dimmap && (this->swaths.size() != 1)) {
2248 missingfield->name = sdim->getName()+"_"+swath->name;
2249 dim = new Dimension(missingfield->name,sdim->getSize());
2250 }
2251 else {
2252 missingfield->name = sdim->getName();
2253 dim = new Dimension(sdim->getName(),sdim->getSize());
2254 }
2255 missingfield->rank = 1;
2256 missingfield->type = DFNT_INT32;//This is an HDF constant.the data type is always integer.
2257
2258 // only 1 dimension
2259 missingfield->dims.push_back(dim);
2260
2261 // Provide information for the missing data, since we need to calculate the data, so
2262 // the information is different than a normal field.
2263 // int missingdimsize[1]; //unused variable. SBL 2/7/20
2264 // missingdimsize[0]= sdim->getSize();
2265
2266 if(0 == sdim->getSize()) {
2267 missingfield_unlim_flag = true;
2268 missingfield_unlim = missingfield;
2269 }
2270
2271 //added Z-dimension coordinate variable with nature number
2272 missingfield->fieldtype = 4;
2273
2274 swath->geofields.push_back(missingfield);
2275 HDFCFUtil::insert_map(swath->dimcvarlist,
2276 (missingfield->getDimensions())[0]->getName(), missingfield->name);
2277 }
2278 }
2279
2280 //Correct the unlimited dimension size.
2281 // The code on the following is ok.
2282 // However, coverity is picky about changing the missingfield_unlim_flag in the middle.
2283 // use a temporary variable for the if block.
2284 // The following code correct the dimension size of unlimited dimension.
2285
2286 bool temp_missingfield_unlim_flag = missingfield_unlim_flag;
2287 if(true == temp_missingfield_unlim_flag) {
2288 for (const auto &dfield:swath->getDataFields()) {
2289
2290 for (const auto &fdim:dfield->getDimensions()) {
2291
2292 if(fdim->getName() == (missingfield_unlim->getDimensions())[0]->getName()) {
2293 if(fdim->getSize()!= 0) {
2294 Dimension *dim = missingfield_unlim->getDimensions()[0];
2295 // Correct the dimension size.
2296 dim->dimsize = fdim->getSize();
2297 missingfield_unlim_flag = false;
2298 break;
2299 }
2300 }
2301
2302 }
2303 if(false == missingfield_unlim_flag)
2304 break;
2305 }
2306 }
2307
2308 swath->nonmisscvdimlist.clear();// clear this set.
2309
2310 }// End of handling non-latlon cv
2311
2312}
2313
2314// Handle swath dimension name to coordinate variable name maps.
2315void File::handle_swath_dim_cvar_maps() throw(Exception) {
2316
2317 // Start handling name clashing
2318 vector <string> tempfieldnamelist;
2319 for (const auto &swath:this->swaths) {
2320
2321 // First handle geofield, all dimension fields are under the geofield group.
2322 for (const auto &gfield:swath->getGeoFields()) {
2323 if(gfield->fieldtype == 0 && (this->swaths.size() !=1) &&
2324 (true == handle_swath_dimmap) &&
2325 (backward_handle_swath_dimmap == false)){
2326 string new_field_name = gfield->name+"_"+swath->name;
2327 tempfieldnamelist.push_back(HDFCFUtil::get_CF_string(new_field_name));
2328 }
2329 else
2330 tempfieldnamelist.push_back(HDFCFUtil::get_CF_string(gfield->name));
2331 }
2332
2333 for (const auto &dfield:swath->getDataFields()) {
2334 if(dfield->fieldtype == 0 && (this->swaths.size() !=1) &&
2335 true == multi_dimmap){
2336 // If we can handle multi dim. maps fro multi swaths, we
2337 // create the field name with the swath name as suffix to
2338 // avoid name clashing.
2339 string new_field_name = dfield->name+"_"+swath->name;
2340 tempfieldnamelist.push_back(HDFCFUtil::get_CF_string(new_field_name));
2341 }
2342 else
2343 tempfieldnamelist.push_back(HDFCFUtil::get_CF_string(dfield->name));
2344 }
2345 }
2346
2347 HDFCFUtil::Handle_NameClashing(tempfieldnamelist);
2348
2349 int total_fcounter = 0;
2350
2351 // Create a map for dimension field name <original field name, corrected field name>
2352 // Also assure the uniqueness of all field names,save the new field names.
2353 //the original dimension field name to the corrected dimension field name
2354 map<string,string>tempncvarnamelist;
2355 for (const auto &swath:this->swaths) {
2356
2357 // First handle geofield, all dimension fields are under the geofield group.
2358 for (const auto &gfield:swath->getGeoFields())
2359 {
2360
2361 gfield->newname = tempfieldnamelist[total_fcounter];
2362 total_fcounter++;
2363
2364 // If this field is a dimension field, save the name/new name pair.
2365 if(gfield->fieldtype!=0) {
2366 HDFCFUtil::insert_map(swath->ncvarnamelist, gfield->getName(), gfield->newname);
2367 }
2368 }
2369
2370 for (const auto &dfield:swath->getDataFields())
2371 {
2372 dfield->newname = tempfieldnamelist[total_fcounter];
2373 total_fcounter++;
2374
2375 // If this field is a dimension field, save the name/new name pair.
2376 if(dfield->fieldtype!=0) {
2377 HDFCFUtil::insert_map(swath->ncvarnamelist, dfield->getName(), dfield->newname);
2378 }
2379 }
2380 } // end of creating a map for dimension field name <original field name, corrected field name>
2381
2382 // Create a map for dimension name < original dimension name, corrected dimension name>
2383
2384 vector <string>tempalldimnamelist;
2385
2386 for (const auto &swath:this->swaths) {
2387 for (map<string,string>::const_iterator j =
2388 swath->dimcvarlist.begin(); j!= swath->dimcvarlist.end();++j)
2389 tempalldimnamelist.push_back(HDFCFUtil::get_CF_string((*j).first));
2390 }
2391
2392 // Handle name clashing will make the corrected dimension name follow CF
2393 HDFCFUtil::Handle_NameClashing(tempalldimnamelist);
2394
2395 int total_dcounter = 0;
2396
2397 for (const auto &swath:this->swaths) {
2398 for (map<string,string>::const_iterator j =
2399 swath->dimcvarlist.begin(); j!= swath->dimcvarlist.end();++j){
2400 HDFCFUtil::insert_map(swath->ndimnamelist, (*j).first, tempalldimnamelist[total_dcounter]);
2401 total_dcounter++;
2402 }
2403 }
2404
2405 // Create corrected dimension vectors.
2406 vector<Dimension*> correcteddims;
2407 string tempcorrecteddimname;
2408 Dimension *correcteddim;
2409
2410 for (const auto &swath:this->swaths) {
2411
2412 // First the geofield.
2413 for (const auto &gfield:swath->getGeoFields()) {
2414
2415 for (const auto &gdim:gfield->getDimensions()) {
2416
2417 map<string,string>::iterator tempmapit;
2418
2419 // No dimension map or dimension names were handled. just obtain the new dimension name.
2420 if(handle_swath_dimmap == false || multi_dimmap == true) {
2421
2422 // Find the new name of this field
2423 tempmapit = swath->ndimnamelist.find(gdim->getName());
2424 if(tempmapit != swath->ndimnamelist.end())
2425 tempcorrecteddimname= tempmapit->second;
2426 else
2427 throw4("cannot find the corrected dimension name",
2428 swath->getName(),gfield->getName(),gdim->getName());
2429
2430 correcteddim = new Dimension(tempcorrecteddimname,gdim->getSize());
2431 }
2432 else {
2433 // have dimension map, use the datadim and datadim size to replace the geodim and geodim size.
2434 bool isdimmapname = false;
2435
2436 // We have to loop through the dimension map
2437 for (const auto &sdmap:swath->getDimensionMaps()) {
2438
2439 // This dimension name is the geo dimension name in the dimension map,
2440 // replace the name with data dimension name.
2441 if(gdim->getName() == sdmap->getGeoDimension()) {
2442
2443 isdimmapname = true;
2444 gfield->dmap = true;
2445 string temprepdimname = sdmap->getDataDimension();
2446
2447 // Find the new name of this data dimension name
2448 tempmapit = swath->ndimnamelist.find(temprepdimname);
2449 if(tempmapit != swath->ndimnamelist.end())
2450 tempcorrecteddimname= tempmapit->second;
2451 else
2452 throw4("cannot find the corrected dimension name", swath->getName(),
2453 gfield->getName(),gdim->getName());
2454
2455 // Find the size of this data dimension name
2456 // We have to loop through the Dimensions of this swath
2457 bool ddimsflag = false;
2458 for (const auto &sdim:swath->getDimensions()) {
2459 if(sdim->getName() == temprepdimname) {
2460 // Find the dimension size, create the correcteddim
2461 correcteddim = new Dimension(tempcorrecteddimname,sdim->getSize());
2462 ddimsflag = true;
2463 break;
2464 }
2465 }
2466 if(!ddimsflag)
2467 throw4("cannot find the corrected dimension size", swath->getName(),
2468 gfield->getName(),gdim->getName());
2469 break;
2470 }
2471 }
2472 if(false == isdimmapname) { // Still need to assign the corrected dimensions.
2473 // Find the new name of this field
2474 tempmapit = swath->ndimnamelist.find(gdim->getName());
2475 if(tempmapit != swath->ndimnamelist.end())
2476 tempcorrecteddimname= tempmapit->second;
2477 else
2478 throw4("cannot find the corrected dimension name",
2479 swath->getName(),gfield->getName(),gdim->getName());
2480
2481 correcteddim = new Dimension(tempcorrecteddimname,gdim->getSize());
2482
2483 }
2484 }
2485
2486 correcteddims.push_back(correcteddim);
2487 }
2488 gfield->setCorrectedDimensions(correcteddims);
2489 correcteddims.clear();
2490 }// End of creating the corrected dimension vectors for GeoFields.
2491
2492 // Then the data field.
2493 for (const auto &dfield:swath->getDataFields()) {
2494
2495 for (const auto &fdim:dfield->getDimensions()) {
2496
2497 if((handle_swath_dimmap == false) || multi_dimmap == true) {
2498
2499 map<string,string>::iterator tempmapit;
2500 // Find the new name of this field
2501 tempmapit = swath->ndimnamelist.find(fdim->getName());
2502 if(tempmapit != swath->ndimnamelist.end())
2503 tempcorrecteddimname= tempmapit->second;
2504 else
2505 throw4("cannot find the corrected dimension name", swath->getName(),
2506 dfield->getName(),fdim->getName());
2507
2508 correcteddim = new Dimension(tempcorrecteddimname,fdim->getSize());
2509 }
2510 else {
2511 map<string,string>::iterator tempmapit;
2512 bool isdimmapname = false;
2513
2514 // We have to loop through dimension map
2515 for (const auto &smap:swath->getDimensionMaps()) {
2516 // This dimension name is the geo dimension name in the dimension map,
2517 // replace the name with data dimension name.
2518 if(fdim->getName() == smap->getGeoDimension()) {
2519 isdimmapname = true;
2520 dfield->dmap = true;
2521 string temprepdimname = smap->getDataDimension();
2522
2523 // Find the new name of this data dimension name
2524 tempmapit = swath->ndimnamelist.find(temprepdimname);
2525 if(tempmapit != swath->ndimnamelist.end())
2526 tempcorrecteddimname= tempmapit->second;
2527 else
2528 throw4("cannot find the corrected dimension name",
2529 swath->getName(),dfield->getName(),fdim->getName());
2530
2531 // Find the size of this data dimension name
2532 // We have to loop through the Dimensions of this swath
2533 bool ddimsflag = false;
2534 for(const auto &sdim:swath->getDimensions()) {
2535
2536 // Find the dimension size, create the correcteddim
2537 if(sdim->getName() == temprepdimname) {
2538 correcteddim = new Dimension(tempcorrecteddimname,sdim->getSize());
2539 ddimsflag = true;
2540 break;
2541 }
2542 }
2543 if(!ddimsflag)
2544 throw4("cannot find the corrected dimension size",
2545 swath->getName(),dfield->getName(),fdim->getName());
2546 break;
2547 }
2548 }
2549 // Not a dimension with dimension map; Still need to assign the corrected dimensions.
2550 if(!isdimmapname) {
2551
2552 // Find the new name of this field
2553 tempmapit = swath->ndimnamelist.find(fdim->getName());
2554 if(tempmapit != swath->ndimnamelist.end())
2555 tempcorrecteddimname= tempmapit->second;
2556 else
2557 throw4("cannot find the corrected dimension name",
2558 swath->getName(),dfield->getName(),fdim->getName());
2559
2560 correcteddim = new Dimension(tempcorrecteddimname,fdim->getSize());
2561 }
2562
2563 }
2564 correcteddims.push_back(correcteddim);
2565 }
2566 dfield->setCorrectedDimensions(correcteddims);
2567 correcteddims.clear();
2568 }// End of creating the dimensions for data fields.
2569 }
2570}
2571
2572// Handle CF attributes for swaths.
2573// The CF attributes include "coordinates", "units" for coordinate variables and "_FillValue".
2574void File::handle_swath_cf_attrs() throw(Exception) {
2575
2576 // Create "coordinates" ,"units" attributes. The "units" attributes only apply to latitude and longitude.
2577 // This is the last round of looping through everything,
2578 // we will match dimension name list to the corresponding dimension field name
2579 // list for every field.
2580 // Since we find some swath files don't specify fillvalue when -9999.0 is found in the real data,
2581 // we specify fillvalue for those fields. This is entirely
2582 // artifical and we will evaluate this approach. KY 2010-3-3
2583
2584 for (const auto &swath:this->swaths) {
2585
2586 // Handle GeoField first.
2587 for (const auto &gfield:swath->getGeoFields()) {
2588
2589 // Real fields: adding the coordinate attribute
2590 if(gfield->fieldtype == 0) {// currently it is always true.
2591 string tempcoordinates="";
2592 string tempfieldname="";
2593 string tempcorrectedfieldname="";
2594 int tempcount = 0;
2595 bool has_ll_coord = false;
2596 if(swath->get_num_map() == 0)
2597 has_ll_coord = true;
2598 else if(handle_swath_dimmap == true) {
2599 if(backward_handle_swath_dimmap == true || multi_dimmap == true)
2600 has_ll_coord = true;
2601 }
2602 for (const auto &dim:gfield->getDimensions()) {
2603
2604 // handle coordinates attributes
2605 map<string,string>::iterator tempmapit;
2606 map<string,string>::iterator tempmapit2;
2607
2608 // Find the dimension field name
2609 tempmapit = (swath->dimcvarlist).find(dim->getName());
2610 if(tempmapit != (swath->dimcvarlist).end())
2611 tempfieldname = tempmapit->second;
2612 else
2613 throw4("cannot find the dimension field name",swath->getName(),
2614 gfield->getName(),dim->getName());
2615
2616 // Find the corrected dimension field name
2617 tempmapit2 = (swath->ncvarnamelist).find(tempfieldname);
2618 if(tempmapit2 != (swath->ncvarnamelist).end())
2619 tempcorrectedfieldname = tempmapit2->second;
2620 else
2621 throw4("cannot find the corrected dimension field name",
2622 swath->getName(),gfield->getName(),dim->getName());
2623
2624 if(false == has_ll_coord)
2625 has_ll_coord= check_ll_in_coords(tempcorrectedfieldname);
2626
2627 if(tempcount == 0)
2628 tempcoordinates= tempcorrectedfieldname;
2629 else
2630 tempcoordinates = tempcoordinates +" "+tempcorrectedfieldname;
2631 tempcount++;
2632 }
2633 if(true == has_ll_coord)
2634 gfield->setCoordinates(tempcoordinates);
2635 }
2636
2637 // Add units for latitude and longitude
2638 // latitude,adding the CF units degrees_north.
2639 if(gfield->fieldtype == 1) {
2640 string tempunits = "degrees_north";
2641 gfield->setUnits(tempunits);
2642 }
2643
2644 // longitude, adding the CF units degrees_east
2645 if(gfield->fieldtype == 2) {
2646 string tempunits = "degrees_east";
2647 gfield->setUnits(tempunits);
2648 }
2649
2650 // Add units for Z-dimension, now it is always "level"
2651 // We decide not touch the units if the third-dimension CV exists(fieldtype =3)
2652 // KY 2013-02-15
2653 //if((gfield->fieldtype == 3)||(gfield->fieldtype == 4))
2654 if(gfield->fieldtype == 4) {
2655 string tempunits ="level";
2656 gfield->setUnits(tempunits);
2657 }
2658
2659 // Add units for "Time",
2660 // Be aware that it is always "days since 1900-01-01 00:00:00"(JIRA HFRHANDLER-167)
2661 if(gfield->fieldtype == 5) {
2662 string tempunits = "days since 1900-01-01 00:00:00";
2663 gfield->setUnits(tempunits);
2664 }
2665 // Set the fill value for floating type data that doesn't have the fill value.
2666 // We found _FillValue attribute is missing from some swath data.
2667 // To cover the most cases, an attribute called _FillValue(the value is -9999.0)
2668 // is added to the data whose type is float32 or float64.
2669 if(((gfield->getFillValue()).empty()) &&
2670 (gfield->getType()==DFNT_FLOAT32 || gfield->getType()==DFNT_FLOAT64)) {
2671 float tempfillvalue = -9999.0;
2672 gfield->addFillValue(tempfillvalue);
2673 gfield->setAddedFillValue(true);
2674 }
2675 }
2676
2677 // Data fields
2678 for (const auto &dfield:swath->getDataFields()) {
2679
2680 // Real fields: adding coordinate attributes
2681 if(dfield->fieldtype == 0) {// currently it is always true.
2682 string tempcoordinates="";
2683 string tempfieldname="";
2684 string tempcorrectedfieldname="";
2685 int tempcount = 0;
2686 bool has_ll_coord = false;
2687 if(swath->get_num_map() == 0)
2688 has_ll_coord = true;
2689 else if(handle_swath_dimmap == true) {
2690 if(backward_handle_swath_dimmap == true || multi_dimmap == true)
2691 has_ll_coord = true;
2692 }
2693 for (const auto &dim:dfield->getDimensions()) {
2694
2695 // handle coordinates attributes
2696 map<string,string>::iterator tempmapit;
2697 map<string,string>::iterator tempmapit2;
2698
2699 // Find the dimension field name
2700 tempmapit = (swath->dimcvarlist).find(dim->getName());
2701 if(tempmapit != (swath->dimcvarlist).end())
2702 tempfieldname = tempmapit->second;
2703 else
2704 throw4("cannot find the dimension field name",swath->getName(),
2705 dfield->getName(),dim->getName());
2706
2707 // Find the corrected dimension field name
2708 tempmapit2 = (swath->ncvarnamelist).find(tempfieldname);
2709 if(tempmapit2 != (swath->ncvarnamelist).end())
2710 tempcorrectedfieldname = tempmapit2->second;
2711 else
2712 throw4("cannot find the corrected dimension field name",
2713 swath->getName(),dfield->getName(),dim->getName());
2714
2715 if(false == has_ll_coord)
2716 has_ll_coord= check_ll_in_coords(tempcorrectedfieldname);
2717
2718 if(tempcount == 0)
2719 tempcoordinates= tempcorrectedfieldname;
2720 else
2721 tempcoordinates = tempcoordinates +" "+tempcorrectedfieldname;
2722 tempcount++;
2723 }
2724 if(true == has_ll_coord)
2725 dfield->setCoordinates(tempcoordinates);
2726 }
2727 // Add units for Z-dimension, now it is always "level"
2728 if((dfield->fieldtype == 3)||(dfield->fieldtype == 4)) {
2729 string tempunits ="level";
2730 dfield->setUnits(tempunits);
2731 }
2732
2733 // Add units for "Time", Be aware that it is always "days since 1900-01-01 00:00:00"
2734 // documented at JIRA (HFRHANDLER-167)
2735 if(dfield->fieldtype == 5) {
2736 string tempunits = "days since 1900-01-01 00:00:00";
2737 dfield->setUnits(tempunits);
2738 }
2739
2740 // Set the fill value for floating type data that doesn't have the fill value.
2741 // We found _FillValue attribute is missing from some swath data.
2742 // To cover the most cases, an attribute called _FillValue(the value is -9999.0)
2743 // is added to the data whose type is float32 or float64.
2744 if(((dfield->getFillValue()).empty()) &&
2745 (dfield->getType()==DFNT_FLOAT32 || dfield->getType()==DFNT_FLOAT64)) {
2746 float tempfillvalue = -9999.0;
2747 dfield->addFillValue(tempfillvalue);
2748 dfield->setAddedFillValue(true);
2749 }
2750 }
2751 }
2752}
2753
2754// Find dimension that has the dimension name.
2755bool File::find_dim_in_dims(const std::vector<Dimension*>&dims,const std::string &dim_name) const {
2756
2757 bool ret_value = false;
2758 for (const auto &dim:dims) {
2759 if (dim->name == dim_name) {
2760 ret_value = true;
2761 break;
2762 }
2763 }
2764 return ret_value;
2765}
2766
2767// Check if the original dimension names in Lat/lon that holds the dimension maps are used by data fields.
2768void File::check_dm_geo_dims_in_vars() {
2769
2770 if(handle_swath_dimmap == false)
2771 return;
2772
2773 for (const auto &swath:this->swaths) {
2774
2775 // Currently we only support swath that has 2-D lat/lon(MODIS).
2776 if(swath->get_num_map() > 0) {
2777
2778 for (const auto &dfield:swath->getDataFields()) {
2779
2780 int match_dims = 0;
2781 // We will only check the variables >=2D since lat/lon are 2D.
2782 if (dfield->rank >=2) {
2783 for (const auto &dim:dfield->getDimensions()) {
2784
2785 // There may be multiple dimension maps that hold the same geo-dimension.
2786 // We should not count this duplicately.
2787 bool not_match_geo_dim = true;
2788 for (const auto &sdmap:swath->getDimensionMaps()) {
2789
2790 if((dim->getName() == sdmap->getGeoDimension()) && not_match_geo_dim){
2791 match_dims++;
2792 not_match_geo_dim = false;
2793 }
2794 }
2795 }
2796 }
2797 // This variable holds the GeoDimensions,this swath
2798 if(match_dims == 2) {
2799 swath->GeoDim_in_vars = true;
2800 break;
2801 }
2802 }
2803
2804 if(swath->GeoDim_in_vars == false) {
2805
2806 for (const auto &gfield:swath->getGeoFields()) {
2807
2808 int match_dims = 0;
2809 // We will only check the variables >=2D since lat/lon are 2D.
2810 if(gfield->rank >=2 && (gfield->name != "Latitude" && gfield->name != "Longitude")) {
2811 for (const auto &dim:gfield->getDimensions()) {
2812
2813 // There may be multiple dimension maps that hold the same geo-dimension.
2814 // We should not count this duplicately.
2815 bool not_match_geo_dim = true;
2816
2817 for (const auto &sdmap:swath->getDimensionMaps()) {
2818 if((dim->getName() == sdmap->getGeoDimension()) && not_match_geo_dim){
2819 match_dims++;
2820 not_match_geo_dim = false;
2821 }
2822 }
2823 }
2824 }
2825 // This variable holds the GeoDimensions,this swath
2826 if(match_dims == 2){
2827 swath->GeoDim_in_vars = true;
2828 break;
2829 }
2830 }
2831 }
2832 }
2833 }
2834
2835 return;
2836}
2837
2838// Based on the dimension name and the mapped dimension name,obtain the offset and increment.
2839// return false if there is no match.
2840bool SwathDataset::obtain_dmap_offset_inc(const string& ori_dimname, const string & mapped_dimname,int &offset,int&inc) {
2841 bool ret_value = false;
2842 for (const auto &sdmap:this->dimmaps) {
2843 if(sdmap->geodim==ori_dimname && sdmap->datadim == mapped_dimname){
2844 offset = sdmap->offset;
2845 inc = sdmap->increment;
2846 ret_value = true;
2847 break;
2848 }
2849 }
2850 return ret_value;
2851}
2852
2853// For the multi-dimension map case, generate all the lat/lon names. The original lat/lon
2854// names should be used.
2855// For one swath, we don't need to provide the swath name. Yes, swath name is missed. However. this is
2856// what users want. For multi swaths, swath names are added.
2857void File::create_geo_varnames_list(vector<string> & geo_varnames,const string & swathname,
2858 const string & fieldname,int extra_ll_pairs,bool oneswath) {
2859 // We will always keep Latitude and Longitude
2860 if(true == oneswath)
2861 geo_varnames.push_back(fieldname);
2862 else {
2863 string nfieldname = fieldname+"_"+swathname;
2864 geo_varnames.push_back(nfieldname);
2865 }
2866 for (int i = 0; i <extra_ll_pairs;i++) {
2867 string nfieldname;
2868 stringstream si;
2869 si << (i+1);
2870 if( true == oneswath) // No swath name is needed.
2871 nfieldname = fieldname+"_"+si.str();
2872 else
2873 nfieldname = fieldname+"_"+swathname+"_"+si.str();
2874 geo_varnames.push_back(nfieldname);
2875 }
2876
2877#if 0
2878cerr<<"ll_pairs is "<<extra_ll_pairs <<endl;
2879for(int i =0;i<geo_varnames.size();i++)
2880 cerr<<"geo_varnames["<<i<<"]= " <<geo_varnames[i] <<endl;
2881#endif
2882}
2883
2884// Make just one routine for both latitude and longtitude dimmaps.
2885// In longitude part, we just ignore ..
2886void File::create_geo_dim_var_maps(SwathDataset*sd, Field*fd,const vector<string>& lat_names,
2887 const vector<string>& lon_names,vector<Dimension*>& geo_var_dim1,
2888 vector<Dimension*>& geo_var_dim2) {
2889 string field_lat_dim1_name =(fd->dims)[0]->name;
2890 string field_lat_dim2_name =(fd->dims)[1]->name;
2891
2892 // Keep the original Latitude/Longitude and the dimensions when GeoDim_in_vars is true.
2893 if(sd->GeoDim_in_vars == true) {
2894 if((this->swaths).size() >1) {
2895 (fd->dims)[0]->name = field_lat_dim1_name+"_"+sd->name;
2896 (fd->dims)[1]->name = field_lat_dim2_name+"_"+sd->name;
2897 }
2898 geo_var_dim1.push_back((fd->dims)[0]);
2899 geo_var_dim2.push_back((fd->dims)[1]);
2900 }
2901
2902 // Create dimension list for the lats and lons.
2903 // Consider the multi-swath case and if the dimension names of orig. lat/lon
2904 // are used.
2905 // We also need to consider one geo-dim can map to multiple data dim.
2906 // One caveat for the current approach is that we don't consider
2907 // two dimension maps are not created in order. HDFEOS2 API implies
2908 // the dimension maps are created in order.
2909 short dim1_map_count = 0;
2910 short dim2_map_count = 0;
2911 for (const auto &sdmap:sd->getDimensionMaps()) {
2912 if(sdmap->getGeoDimension()==field_lat_dim1_name){
2913 string data_dim1_name = sdmap->getDataDimension();
2914 int dim1_size = sd->obtain_dimsize_with_dimname(data_dim1_name);
2915 if((this->swaths).size() > 1)
2916 data_dim1_name = data_dim1_name+"_"+sd->name;
2917
2918 if(sd->GeoDim_in_vars == false && dim1_map_count == 0) {
2919 (fd->dims)[0]->name = data_dim1_name;
2920 (fd->dims)[0]->dimsize = dim1_size;
2921 geo_var_dim1.push_back((fd->dims)[0]);
2922 }
2923 else {
2924 auto lat_dim = new Dimension(data_dim1_name,dim1_size);
2925 geo_var_dim1.push_back(lat_dim);
2926 }
2927 dim1_map_count++;
2928 }
2929 else if(sdmap->getGeoDimension()==field_lat_dim2_name){
2930 string data_dim2_name = sdmap->getDataDimension();
2931 int dim2_size = sd->obtain_dimsize_with_dimname(data_dim2_name);
2932 if((this->swaths).size() > 1)
2933 data_dim2_name = data_dim2_name+"_"+sd->name;
2934 if(sd->GeoDim_in_vars == false && dim2_map_count == 0) {
2935 (fd->dims)[1]->name = data_dim2_name;
2936 (fd->dims)[1]->dimsize = dim2_size;
2937 geo_var_dim2.push_back((fd->dims)[1]);
2938 }
2939 else {
2940 auto lon_dim = new Dimension(data_dim2_name,dim2_size);
2941 geo_var_dim2.push_back(lon_dim);
2942 }
2943 dim2_map_count++;
2944 }
2945 }
2946
2947 // Build up dimension names to coordinate var lists.
2948 for(int i = 0; i<lat_names.size();i++) {
2949 HDFCFUtil::insert_map(sd->dimcvarlist, (geo_var_dim1[i])->name,lat_names[i]);
2950 HDFCFUtil::insert_map(sd->dimcvarlist, (geo_var_dim2[i])->name,lon_names[i]);
2951 }
2952
2953 return;
2954}
2955
2956// Generate lat/lon variables for the multi-dimension map case for this swath.
2957// Original lat/lon variable information is provided.
2958void File::create_geo_vars(SwathDataset* sd,Field *orig_lat,Field*orig_lon,
2959 const vector<string>& lat_names,const vector<string>& lon_names,
2960 vector<Dimension*>&geo_var_dim1,vector<Dimension*>&geo_var_dim2) throw(Exception){
2961
2962#if 0
2963 // Handle existing latitude and longitude.
2964 // If we don't need to keep GeoDim in the latitude and longitude,
2965 // dimensions of latitude and longitude need to be updated.
2966 Field* orig_lat;
2967 Field* orig_lon;
2968
2969 // Here we don't need to search the data fields,
2970 // we only support the standard Swath:lat/lon under /geolocation fields.
2971 for (vector<Field *>::iterator i = sd->geofields.begin();
2972 i!=sd->geofields.end();++i) {
2973 if((*i)->name == "Latitude")
2974 orig_lat =(*i);
2975 else if((*i)->name == "Longitude")
2976 orig_lon =(*i);
2977 }
2978#endif
2979
2980 // Need to have ll dimension names to obtain the dimension maps
2981 string ll_ori_dim0_name = (orig_lon->dims)[0]->name;
2982 string ll_ori_dim1_name = (orig_lon->dims)[1]->name;
2983 int dmap_offset = 0;
2984 int dmap_inc = 0;
2985 if(sd->GeoDim_in_vars == false) {
2986
2987
2988#if 0
2989 (orig_lat->dims)[0]->name = geo_var_dim1[0]->name;
2990 (orig_lat->dims)[0]->dimsize = geo_var_dim1[0]->dimsize;
2991 (orig_lat->dims)[1]->name = geo_var_dim2[0]->name;
2992 (orig_lat->dims)[1]->dimsize = geo_var_dim2[0]->dimsize;
2993#endif
2994 // The original lat's dim has been updated in create_geo_dim_var_maps
2995 // In theory , it should be reasonable to do it here. Later.
2996 (orig_lon->dims)[0]->name = geo_var_dim1[0]->name;
2997 (orig_lon->dims)[0]->dimsize = geo_var_dim1[0]->dimsize;
2998 (orig_lon->dims)[1]->name = geo_var_dim2[0]->name;
2999 (orig_lon->dims)[1]->dimsize = geo_var_dim2[0]->dimsize;
3000 string ll_datadim0_name = geo_var_dim1[0]->name;
3001 string ll_datadim1_name = geo_var_dim2[0]->name;
3002 if(this->swaths.size() >1) {
3003 string prefix_remove = "_"+sd->name;
3004 ll_datadim0_name = ll_datadim0_name.substr(0,ll_datadim0_name.size()-prefix_remove.size());
3005 ll_datadim1_name = ll_datadim1_name.substr(0,ll_datadim1_name.size()-prefix_remove.size());
3006 }
3007
3008 // dimension map offset and inc should be retrieved.
3009 if(false == sd->obtain_dmap_offset_inc(ll_ori_dim0_name,ll_datadim0_name,dmap_offset,dmap_inc)){
3010 throw5("Cannot retrieve dimension map offset and inc ",sd->name,
3011 orig_lon->name,ll_ori_dim0_name,ll_datadim0_name);
3012 }
3013 orig_lon->ll_dim0_inc = dmap_inc;
3014 orig_lon->ll_dim0_offset = dmap_offset;
3015 orig_lat->ll_dim0_inc = dmap_inc;
3016 orig_lat->ll_dim0_offset = dmap_offset;
3017
3018 if(false == sd->obtain_dmap_offset_inc(ll_ori_dim1_name,ll_datadim1_name,dmap_offset,dmap_inc)){
3019 throw5("Cannot retrieve dimension map offset and inc ",sd->name,
3020 orig_lon->name,ll_ori_dim1_name,ll_datadim1_name);
3021 }
3022 orig_lon->ll_dim1_inc = dmap_inc;
3023 orig_lon->ll_dim1_offset = dmap_offset;
3024 orig_lat->ll_dim1_inc = dmap_inc;
3025 orig_lat->ll_dim1_offset = dmap_offset;
3026#if 0
3027cerr<<"orig_lon "<<orig_lon->name <<endl;
3028cerr<<"orig_lon dim1 inc "<<orig_lon->ll_dim1_inc<<endl;
3029cerr<<"orig_lon dim1 offset "<<orig_lon->ll_dim1_offset<<endl;
3030cerr<<"orig_lon dim0 inc "<<orig_lon->ll_dim0_inc<<endl;
3031cerr<<"orig_lon dim0 offset "<<orig_lon->ll_dim0_offset<<endl;
3032#endif
3033
3034
3035 }
3036 else {
3037 // if GeoDim is used, we still need to update longitude's dimension names
3038 // if multiple swaths. Latitude was done in create_geo_dim_var_maps().
3039 if((this->swaths).size() >1) {
3040 (orig_lon->dims)[0]->name = (orig_lon->dims)[0]->name + "_" + sd->name;
3041 (orig_lon->dims)[1]->name = (orig_lon->dims)[1]->name + "_" + sd->name;
3042 }
3043 }
3044
3045 // We also need to update the latitude and longitude names when num_swath is not 1.
3046 if((this->swaths).size()>1) {
3047 orig_lat->name = lat_names[0];
3048 orig_lon->name = lon_names[0];
3049 }
3050
3051 // Mark the original lat/lon as coordinate variables.
3052 orig_lat->fieldtype = 1;
3053 orig_lon->fieldtype = 2;
3054
3055 // The added fields.
3056 for (int i = 1; i <lat_names.size();i++) {
3057 auto newlat = new Field();
3058 newlat->name = lat_names[i];
3059 (newlat->dims).push_back(geo_var_dim1[i]);
3060 (newlat->dims).push_back(geo_var_dim2[i]);
3061 newlat->fieldtype = 1;
3062 newlat->rank = 2;
3063 newlat->type = orig_lat->type;
3064 auto newlon = new Field();
3065 newlon->name = lon_names[i];
3066 // Here we need to create new Dimensions
3067 // for Longitude.
3068 auto lon_dim1=
3069 new Dimension(geo_var_dim1[i]->name,geo_var_dim1[i]->dimsize);
3070 auto lon_dim2=
3071 new Dimension(geo_var_dim2[i]->name,geo_var_dim2[i]->dimsize);
3072 (newlon->dims).push_back(lon_dim1);
3073 (newlon->dims).push_back(lon_dim2);
3074 newlon->fieldtype = 2;
3075 newlon->rank = 2;
3076 newlon->type = orig_lon->type;
3077
3078 string ll_datadim0_name = geo_var_dim1[i]->name;
3079 string ll_datadim1_name = geo_var_dim2[i]->name;
3080 if(this->swaths.size() >1) {
3081 string prefix_remove = "_"+sd->name;
3082 ll_datadim0_name = ll_datadim0_name.substr(0,ll_datadim0_name.size()-prefix_remove.size());
3083 ll_datadim1_name = ll_datadim1_name.substr(0,ll_datadim1_name.size()-prefix_remove.size());
3084 }
3085
3086 // Obtain dimension map offset and inc for the new lat/lon.
3087 if(false == sd->obtain_dmap_offset_inc(ll_ori_dim0_name,ll_datadim0_name,dmap_offset,dmap_inc)){
3088 throw5("Cannot retrieve dimension map offset and inc ",sd->name,
3089 newlon->name,ll_ori_dim0_name,ll_datadim0_name);
3090 }
3091 newlon->ll_dim0_inc = dmap_inc;
3092 newlon->ll_dim0_offset = dmap_offset;
3093 newlat->ll_dim0_inc = dmap_inc;
3094 newlat->ll_dim0_offset = dmap_offset;
3095 if(false == sd->obtain_dmap_offset_inc(ll_ori_dim1_name,ll_datadim1_name,dmap_offset,dmap_inc)){
3096 throw5("Cannot retrieve dimension map offset and inc ",sd->name,
3097 newlon->name,ll_ori_dim0_name,ll_datadim1_name);
3098 }
3099 newlon->ll_dim1_inc = dmap_inc;
3100 newlon->ll_dim1_offset = dmap_offset;
3101 newlat->ll_dim1_inc = dmap_inc;
3102 newlat->ll_dim1_offset = dmap_offset;
3103#if 0
3104cerr<<"newlon "<<newlon->name <<endl;
3105cerr<<"newlon dim1 inc "<<newlon->ll_dim1_inc<<endl;
3106cerr<<"newlon dim1 offset "<<newlon->ll_dim1_offset<<endl;
3107cerr<<"newlon dim0 inc "<<newlon->ll_dim0_inc<<endl;
3108cerr<<"newlon dim0 offset "<<newlon->ll_dim0_offset<<endl;
3109#endif
3110 sd->geofields.push_back(newlat);
3111 sd->geofields.push_back(newlon);
3112 }
3113#if 0
3114//cerr<<"Latlon under swath: "<<sd->name<<endl;
3115for (vector<Field *>::const_iterator j =
3116 sd->getGeoFields().begin();
3117 j != sd->getGeoFields().end(); ++j) {
3118cerr<<"Field name: "<<(*j)->name <<endl;
3119cerr<<"Dimension name and size: "<<endl;
3120 for(vector<Dimension *>::const_iterator k=
3121 (*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k)
3122cerr<<(*k)->getName() <<": "<<(*k)->getSize() <<endl;
3123}
3124#endif
3125
3126}
3127
3128// We need to update the dimensions for all the swath and all the fields when
3129// we can handle the multi-dimension map case. This only applies to >1 swath.
3130// For one swath, dimension names are not needed to be updated.
3131void File::update_swath_dims_for_dimmap(SwathDataset* sd,const std::vector<Dimension*>&geo_var_dim1,
3132 const std::vector<Dimension*>&geo_var_dim2) {
3133
3134 // Loop through each field under geofields and data fields. update dimensions.
3135 // Obtain each dimension name + _+swath_name, if match with geo_var_dim1 or geo_var_dim2;
3136 // Update the dimension names with the matched one.
3137 for (const auto &gfield:sd->getGeoFields()) {
3138 // No need to update latitude/longitude
3139 if(gfield->fieldtype == 1 || gfield->fieldtype == 2)
3140 continue;
3141 for (const auto &dim:gfield->getDimensions()) {
3142 string new_dim_name = dim->name +"_"+sd->name;
3143 if (find_dim_in_dims(geo_var_dim1,new_dim_name) ||
3144 find_dim_in_dims(geo_var_dim2,new_dim_name))
3145 dim->name = new_dim_name;
3146 }
3147 }
3148
3149 for (const auto &dfield:sd->getDataFields()){
3150 for (const auto &dim:dfield->getDimensions()) {
3151 string new_dim_name = dim->name +"_"+sd->name;
3152 if(find_dim_in_dims(geo_var_dim1,new_dim_name) ||
3153 find_dim_in_dims(geo_var_dim2,new_dim_name))
3154 dim->name = new_dim_name;
3155 }
3156 }
3157
3158 // We also need to update the dimension name of this swath.
3159 for (const auto &dim:sd->getDimensions()) {
3160 string new_dim_name = dim->name +"_"+sd->name;
3161 if(find_dim_in_dims(geo_var_dim1,new_dim_name) ||
3162 find_dim_in_dims(geo_var_dim2,new_dim_name))
3163 dim->name = new_dim_name;
3164 }
3165
3166 return;
3167
3168}
3169
3170// This is the main function to handle the multi-dimension map case.
3171// It creates the lat/lon lists, handle dimension names and then
3172// provide the dimension name to coordinate variable map.
3173void File::create_swath_latlon_dim_cvar_map_for_dimmap(SwathDataset* sd, Field* ori_lat, Field* ori_lon) throw(Exception) {
3174
3175 bool one_swath = ((this->swaths).size() == 1);
3176
3177 int num_extra_lat_lon_pairs = 0;
3178
3179#if 0
3180if(sd->GeoDim_in_vars == true)
3181 cerr<<" swath name is "<<sd->name <<endl;
3182#endif
3183
3184 // Since the original lat/lon will be kept for lat/lon with the first dimension map..
3185 if(sd->GeoDim_in_vars == false)
3186 num_extra_lat_lon_pairs--;
3187
3188 num_extra_lat_lon_pairs += (sd->num_map)/2;
3189
3190 vector<string> lat_names;
3191 create_geo_varnames_list(lat_names,sd->name,ori_lat->name,num_extra_lat_lon_pairs,one_swath);
3192 vector<string>lon_names;
3193 create_geo_varnames_list(lon_names,sd->name,ori_lon->name,num_extra_lat_lon_pairs,one_swath);
3194 vector<Dimension*> geo_var_dim1;
3195 vector<Dimension*> geo_var_dim2;
3196
3197 // Define dimensions or obtain dimensions for new field.
3198 create_geo_dim_var_maps(sd, ori_lat,lat_names,lon_names,geo_var_dim1,geo_var_dim2);
3199 create_geo_vars(sd,ori_lat,ori_lon,lat_names,lon_names,geo_var_dim1,geo_var_dim2);
3200
3201 // Update dims for vars,this is only necessary when there are multiple swaths
3202 // Dimension names need to be updated to include swath names.
3203 if((this->swaths).size() >1)
3204 update_swath_dims_for_dimmap(sd,geo_var_dim1,geo_var_dim2);
3205
3206}
3207
3211void File::Prepare(const char *eosfile_path) throw(Exception)
3212{
3213
3214 // check if this is a special HDF-EOS2 grid(MOD08_M3) that have all dimension scales
3215 // added by the HDF4 library but the relation between the dimension scale and the dimension is not
3216 // specified. If the return value is true, we will specify
3217
3218 // Obtain the number of swaths and the number of grids
3219 auto numgrid = (int)(this->grids.size());
3220 auto numswath = (int)(this->swaths.size());
3221
3222 if(numgrid < 0)
3223 throw2("the number of grid is less than 0", eosfile_path);
3224
3225 // First handle grids
3226 if (numgrid > 0) {
3227
3228 // Obtain "XDim","YDim","Latitude","Longitude" and "location" set.
3229 string DIMXNAME = this->get_geodim_x_name();
3230
3231 string DIMYNAME = this->get_geodim_y_name();
3232
3233 string LATFIELDNAME = this->get_latfield_name();
3234
3235 string LONFIELDNAME = this->get_lonfield_name();
3236
3237 // Now only there is only one geo grid name "location", so don't call it know.
3238 string GEOGRIDNAME = "location";
3239
3240 // Check global lat/lon for multiple grids.
3241 // We want to check if there is a global lat/lon for multiple grids.
3242 // AIRS level 3 data provides lat/lon under the GEOGRIDNAME grid.
3243 check_onelatlon_grids();
3244
3245 // Handle the third-dimension(both existing and missing) coordinate variables
3246 for (const auto &grid:this->grids)
3247 handle_one_grid_zdim(grid);
3248
3249
3250 // Handle lat/lon fields for the case of which all grids have one dedicated lat/lon grid.
3251 if (true == this->onelatlon)
3252 handle_onelatlon_grids();
3253 else {
3254 for (const auto &grid:this->grids) {
3255
3256 // Set the horizontal dimension name "dimxname" and "dimyname"
3257 // This will be used to detect the dimension major order.
3258 grid->setDimxName(DIMXNAME);
3259 grid->setDimyName(DIMYNAME);
3260
3261 // Handle lat/lon(both existing lat/lon and calculated lat/lon from EOS2 APIs)
3262 handle_one_grid_latlon(grid);
3263 }
3264 }
3265
3266 // Handle the dimension name to coordinate variable map for grid.
3267 handle_grid_dim_cvar_maps();
3268
3269 // Follow COARDS for grids.
3270 handle_grid_coards();
3271
3272 // Create the corrected dimension vector for each field when COARDS is not followed.
3273 update_grid_field_corrected_dims();
3274
3275 // Handle CF attributes for grids.
3276 handle_grid_cf_attrs();
3277
3278 // Special handling SOM(Space Oblique Mercator) projection files
3279 handle_grid_SOM_projection();
3280
3281
3282 }// End of handling grid
3283
3284 // Check and set the scale type
3285 for (const auto& grid:this->grids)
3286 grid->SetScaleType(grid->name);
3287
3288 if(numgrid==0) {
3289
3290 // Now we handle swath case.
3291 if (numswath > 0) {
3292
3293 // Check if we need to handle dimension maps in this file.
3294 check_swath_dimmap(numswath);
3295
3296 // Check if GeoDim is used by variables when dimension maps are present.
3297 check_dm_geo_dims_in_vars();
3298
3299 // If we need to handle dimension maps,check if we need to keep the old way.
3300 check_swath_dimmap_bk_compat(numswath);
3301
3302 // Create the dimension name to coordinate variable name map for lat/lon.
3303 create_swath_latlon_dim_cvar_map();
3304
3305 // Create the dimension name to coordinate variable name map for non lat/lon coordinate variables.
3306 create_swath_nonll_dim_cvar_map();
3307
3308 // Handle swath dimension name to coordinate variable name maps.
3309 handle_swath_dim_cvar_maps();
3310
3311 // Handle CF attributes for swaths.
3312 // The CF attributes include "coordinates", "units" for coordinate variables and "_FillValue".
3313 handle_swath_cf_attrs();
3314
3315 // Check and set the scale type
3316 for (const auto &swath:this->swaths)
3317 swath->SetScaleType(swath->name);
3318 }
3319
3320 }// End of handling swath
3321
3322}
3323
3324
3325
3326#if 0
3327void correct_unlimited_missing_zdim(GridDataset* gdset) throw(Exception) {
3328
3329 for (vector<Field *>::const_iterator j =
3330 gdset->getDataFields().begin();
3331 j != gdset->getDataFields().end(); ++j) {
3332
3333 //We only need to search those 1-D fields
3334 if ((*j)->getRank()==1 && (*j)->){
3335
3336
3337
3338 }
3339
3340 }
3341}
3342#endif
3343
3344bool File::check_special_1d_grid() throw(Exception) {
3345
3346 auto numgrid = (int)(this->grids.size());
3347 auto numswath = (int)(this->swaths.size());
3348
3349 if (numgrid != 1 || numswath != 0)
3350 return false;
3351
3352 // Obtain "XDim","YDim","Latitude","Longitude" and "location" set.
3353 string DIMXNAME = this->get_geodim_x_name();
3354 string DIMYNAME = this->get_geodim_y_name();
3355
3356 if(DIMXNAME != "XDim" || DIMYNAME != "YDim")
3357 return false;
3358
3359 int var_dimx_size = 0;
3360 int var_dimy_size = 0;
3361
3362 GridDataset *mygrid = (this->grids)[0];
3363
3364 int field_xydim_flag = 0;
3365 for (const auto &dfield:mygrid->getDataFields()) {
3366 if(1==dfield->rank) {
3367 if(dfield->name == "XDim"){
3368 field_xydim_flag++;
3369 var_dimx_size = (dfield->getDimensions())[0]->getSize();
3370 }
3371 if(dfield->name == "YDim"){
3372 field_xydim_flag++;
3373 var_dimy_size = (dfield->getDimensions())[0]->getSize();
3374 }
3375 }
3376 if(2==field_xydim_flag)
3377 break;
3378 }
3379
3380 if(field_xydim_flag !=2)
3381 return false;
3382
3383 // Obtain XDim and YDim size.
3384 int xdimsize = mygrid->getInfo().getX();
3385 int ydimsize = mygrid->getInfo().getY();
3386
3387 if(var_dimx_size != xdimsize || var_dimy_size != ydimsize)
3388 return false;
3389
3390 return true;
3391
3392}
3393
3394
3395bool File::check_ll_in_coords(const string& vname) throw(Exception) {
3396
3397 bool ret_val = false;
3398 for (const auto &swath:this->swaths) {
3399 for (const auto &gfield:swath->getGeoFields()) {
3400 // Real fields: adding the coordinate attribute
3401 if(gfield->fieldtype == 1 || gfield->fieldtype == 2) {// currently it is always true.
3402 if(gfield->getNewName() == vname) {
3403 ret_val = true;
3404 break;
3405 }
3406 }
3407 }
3408 if(true == ret_val)
3409 break;
3410 for (const auto &dfield:swath->getDataFields()) {
3411
3412 // Real fields: adding the coordinate attribute
3413 if(dfield->fieldtype == 1 || dfield->fieldtype == 2) {// currently it is always true.
3414 if(dfield->getNewName() == vname) {
3415 ret_val = true;
3416 break;
3417 }
3418 }
3419 }
3420 if(true == ret_val)
3421 break;
3422
3423 }
3424 return ret_val;
3425}
3426
3427
3428// Set scale and offset type
3429// MODIS data has three scale and offset rules.
3430// They are
3431// MODIS_EQ_SCALE: raw_data = scale*data + offset
3432// MODIS_MUL_SCALE: raw_data = scale*(data -offset)
3433// MODIS_DIV_SCALE: raw_data = (data-offset)/scale
3434
3435void Dataset::SetScaleType(const string & EOS2ObjName) throw(Exception) {
3436
3437
3438 // Group features of MODIS products.
3439 // Using vector of strings instead of the following.
3440 // C++11 may allow the vector of string to be assigned as follows
3441 // string modis_type1[] = {"L1B", "GEO", "BRDF", "0.05Deg", "Reflectance", "MOD17A2", "North","South","MOD_Swath_Sea_Ice","MOD_Grid_MOD15A2","MODIS_NACP_LAI"};
3442
3443 vector<string> modis_multi_scale_type;
3444 modis_multi_scale_type.push_back("L1B");
3445 modis_multi_scale_type.push_back("GEO");
3446 modis_multi_scale_type.push_back("BRDF");
3447 modis_multi_scale_type.push_back("0.05Deg");
3448 modis_multi_scale_type.push_back("Reflectance");
3449 modis_multi_scale_type.push_back("MOD17A2");
3450 modis_multi_scale_type.push_back("North");
3451 modis_multi_scale_type.push_back("South");
3452 modis_multi_scale_type.push_back("MOD_Swath_Sea_Ice");
3453 modis_multi_scale_type.push_back("MOD_Grid_MOD15A2");
3454 modis_multi_scale_type.push_back("MOD_Grid_MOD16A2");
3455 modis_multi_scale_type.push_back("MOD_Grid_MOD16A3");
3456 modis_multi_scale_type.push_back("MODIS_NACP_LAI");
3457
3458 vector<string> modis_div_scale_type;
3459
3460 // Always start with MOD or mod.
3461 modis_div_scale_type.push_back("VI");
3462 modis_div_scale_type.push_back("1km_2D");
3463 modis_div_scale_type.push_back("L2g_2d");
3464 modis_div_scale_type.push_back("CMG");
3465 modis_div_scale_type.push_back("MODIS SWATH TYPE L2");
3466
3467 // This one doesn't start with "MOD" or "mod".
3468 //modis_div_scale_type.push_back("VIP_CMG_GRID");
3469
3470 string modis_eq_scale_type = "LST";
3471 string modis_equ_scale_lst_group1="MODIS_Grid_8Day_1km_LST21";
3472 string modis_equ_scale_lst_group2="MODIS_Grid_Daily_1km_LST21";
3473
3474 string modis_divequ_scale_group = "MODIS_Grid";
3475 string modis_div_scale_group = "MOD_Grid";
3476
3477 // The scale-offset rule for the following group may be MULTI but since add_offset is 0. So
3478 // the MULTI rule is equal to the EQU rule. KY 2013-01-25
3479 string modis_equ_scale_group = "MODIS_Grid_1km_2D";
3480
3481 if(EOS2ObjName=="mod05" || EOS2ObjName=="mod06" || EOS2ObjName=="mod07"
3482 || EOS2ObjName=="mod08" || EOS2ObjName=="atml2")
3483 {
3484 scaletype = MODIS_MUL_SCALE;
3485 return;
3486 }
3487
3488 // Find one MYD09GA2012.version005 file that
3489 // the grid names change to MODIS_Grid_500m_2D.
3490 // So add this one. KY 2012-11-20
3491
3492 // Find the grid name in one MCD43C1 file starts with "MCD_CMG_BRDF",
3493 // however, the offset of the data is 0.
3494 // So we may not handle this file here since it follows the CF conventions.
3495 // May need to double check them later. KY 2013-01-24
3496
3497
3498 if(EOS2ObjName.find("MOD")==0 || EOS2ObjName.find("mod")==0)
3499 {
3500 size_t pos = EOS2ObjName.rfind(modis_eq_scale_type);
3501 if(pos != string::npos &&
3502 (pos== (EOS2ObjName.length()-modis_eq_scale_type.length())))
3503 {
3504 scaletype = MODIS_EQ_SCALE;
3505 return;
3506 }
3507
3508 for(unsigned int k=0; k<modis_multi_scale_type.size(); k++)
3509 {
3510 pos = EOS2ObjName.rfind(modis_multi_scale_type[k]);
3511 if (pos !=string::npos &&
3512 (pos== (EOS2ObjName.length()-modis_multi_scale_type[k].length())))
3513 {
3514 scaletype = MODIS_MUL_SCALE;
3515 return;
3516 }
3517 }
3518
3519 for(unsigned int k=0; k<modis_div_scale_type.size(); k++)
3520 {
3521 pos = EOS2ObjName.rfind(modis_div_scale_type[k]);
3522 if (pos != string::npos &&
3523 (pos==(EOS2ObjName.length()-modis_div_scale_type[k].length()))){
3524 scaletype = MODIS_DIV_SCALE;
3525
3526 // We have a case that group
3527 // MODIS_Grid_1km_2D should apply the equal scale equation.
3528 // This will be handled after this loop.
3529 if (EOS2ObjName != "MODIS_Grid_1km_2D")
3530 return;
3531 }
3532 }
3533
3534 // Special handling for MOD_Grid and MODIS_Grid_500m_2D.
3535 // Check if the group name starts with the modis_divequ and modis_div_scale.
3536 pos = EOS2ObjName.find(modis_divequ_scale_group);
3537
3538 // Find the "MODIS_Grid???" group.
3539 // We have to separate MODIS_Grid_1km_2D(EQ),MODIS_Grid_8Day_1km_LST21
3540 // and MODIS_Grid_Daily_1km_LST21 from other grids(DIV).
3541 if (0 == pos) {
3542 size_t eq_scale_pos = EOS2ObjName.find(modis_equ_scale_group)
3543 *EOS2ObjName.find(modis_equ_scale_lst_group1)
3544 *EOS2ObjName.find(modis_equ_scale_lst_group2);
3545 if (0 == eq_scale_pos)
3546 scaletype = MODIS_EQ_SCALE;
3547 else
3548 scaletype = MODIS_DIV_SCALE;
3549 }
3550 else {
3551 size_t div_scale_pos = EOS2ObjName.find(modis_div_scale_group);
3552
3553 // Find the "MOD_Grid???" group.
3554 if ( 0 == div_scale_pos)
3555 scaletype = MODIS_DIV_SCALE;
3556 }
3557 }
3558
3559 // MEASuRES VIP files have the grid name VIP_CMG_GRID.
3560 // This applies to all VIP version 2 files. KY 2013-01-24
3561 if (EOS2ObjName =="VIP_CMG_GRID")
3562 scaletype = MODIS_DIV_SCALE;
3563}
3564
3565int Dataset::obtain_dimsize_with_dimname(const string & dimname) {
3566 int ret_value = -1;
3567 for(vector<Dimension *>::const_iterator k=
3568 this->getDimensions().begin();k!=this->getDimensions().end();++k){
3569 if((*k)->name == dimname){
3570 ret_value = (*k)->dimsize;
3571 break;
3572 }
3573 }
3574 return ret_value;
3575}
3576
3577// Release resources
3578Field::~Field()
3579{
3580 for_each(this->dims.begin(), this->dims.end(), delete_elem());
3581 for_each(this->correcteddims.begin(), this->correcteddims.end(), delete_elem());
3582}
3583
3584// Release resources
3585Dataset::~Dataset()
3586{
3587 for_each(this->dims.begin(), this->dims.end(), delete_elem());
3588 for_each(this->datafields.begin(), this->datafields.end(),
3589 delete_elem());
3590 for_each(this->attrs.begin(), this->attrs.end(), delete_elem());
3591}
3592
3593// Retrieve dimensions of grids or swaths
3594void Dataset::ReadDimensions(int32 (*entries)(int32, int32, int32 *),
3595 int32 (*inq)(int32, char *, int32 *),
3596 vector<Dimension *> &d_dims) throw(Exception)
3597{
3598 // Initialize number of dimensions
3599 int32 numdims = 0;
3600
3601 // Initialize buf size
3602 int32 bufsize = 0;
3603
3604 // Obtain the number of dimensions and buffer size of holding ","
3605 // separated dimension name list.
3606 if ((numdims = entries(this->datasetid, HDFE_NENTDIM, &bufsize)) == -1)
3607 throw2("dimension entry", this->name);
3608
3609 // Read all dimension information
3610 if (numdims > 0) {
3611 vector<char> namelist;
3612 vector<int32> dimsize;
3613
3614 namelist.resize(bufsize + 1);
3615 dimsize.resize(numdims);
3616
3617 // Inquiry dimension name lists and sizes for all dimensions
3618 if (inq(this->datasetid, namelist.data(), dimsize.data()) == -1)
3619 throw2("inquire dimension", this->name);
3620
3621 vector<string> dimnames;
3622
3623 // Make the "," separated name string to a string list without ",".
3624 // This split is for global dimension of a Swath or a Grid object.
3625 HDFCFUtil::Split(namelist.data(), bufsize, ',', dimnames);
3626 int count = 0;
3627 for (vector<string>::const_iterator i = dimnames.begin();
3628 i != dimnames.end(); ++i) {
3629 Dimension *dim = new Dimension(*i, dimsize[count]);
3630 d_dims.push_back(dim);
3631 ++count;
3632 }
3633 }
3634}
3635
3636// Retrieve data field info.
3637void Dataset::ReadFields(int32 (*entries)(int32, int32, int32 *),
3638 int32 (*inq)(int32, char *, int32 *, int32 *),
3639 intn (*fldinfo)
3640 (int32, char *, int32 *, int32 *, int32 *, char *),
3641 intn (*readfld) //unused SBL 2/7/20
3642 (int32, char *, int32 *, int32 *, int32 *, VOIDP),
3643 intn (*getfill)(int32, char *, VOIDP),
3644 bool geofield,
3645 vector<Field *> &fields) throw(Exception)
3646{
3647
3648 // Initalize the number of fields
3649 int32 numfields = 0;
3650
3651 // Initalize the buffer size
3652 int32 bufsize = 0;
3653
3654 // Obtain the number of fields and buffer size for the current Swath or
3655 // Grid.
3656 if ((numfields = entries(this->datasetid, geofield ?
3657 HDFE_NENTGFLD : HDFE_NENTDFLD, &bufsize)) == -1)
3658 throw2("field entry", this->name);
3659
3660 // Obtain the information of fields (either data fields or geo-location
3661 // fields of a Swath object)
3662 if (numfields > 0) {
3663 vector<char> namelist;
3664
3665 namelist.resize(bufsize + 1);
3666
3667 // Inquiry fieldname list of the current object
3668 if (inq(this->datasetid, namelist.data(), nullptr, nullptr) == -1)
3669 throw2("inquire field", this->name);
3670
3671 vector<string> fieldnames;
3672
3673 // Split the field namelist, make the "," separated name string to a
3674 // string list without ",".
3675 HDFCFUtil::Split(namelist.data(), bufsize, ',', fieldnames);
3676 for (vector<string>::const_iterator i = fieldnames.begin();
3677 i != fieldnames.end(); ++i) {
3678
3679 Field *field = new Field();
3680 if(field == nullptr)
3681 throw1("The field is nullptr");
3682 field->name = *i;
3683
3684 bool throw_error = false;
3685 string err_msg;
3686
3687 // XXX: We assume the maximum number of dimension for an EOS field
3688 // is 16.
3689 int32 dimsize[16];
3690 char dimlist[512]; // XXX: what an HDF-EOS2 developer recommended
3691
3692 // Obtain most information of a field such as rank, dimension etc.
3693 if ((fldinfo(this->datasetid,
3694 const_cast<char *>(field->name.c_str()),
3695 &field->rank, dimsize, &field->type, dimlist)) == -1){
3696 string fieldname_for_eh = field->name;
3697 throw_error = true;
3698 err_msg ="Obtain field info error for field name "+fieldname_for_eh;
3699 }
3700
3701 if(false == throw_error) {
3702
3703 vector<string> dimnames;
3704
3705 // Split the dimension name list for a field
3706 HDFCFUtil::Split(dimlist, ',', dimnames);
3707 if ((int)dimnames.size() != field->rank) {
3708 throw_error = true;
3709 err_msg = "Dimension names size is not consistent with field rank. ";
3710 err_msg += "Field name is "+field->name;
3711 }
3712 else {
3713 for (int k = 0; k < field->rank; ++k) {
3714 Dimension *dim = new Dimension(dimnames[k], dimsize[k]);
3715 field->dims.push_back(dim);
3716 }
3717
3718 // Get fill value of a field
3719 field->filler.resize(DFKNTsize(field->type));
3720 if (getfill(this->datasetid,
3721 const_cast<char *>(field->name.c_str()),
3722 &field->filler[0]) == -1)
3723 field->filler.clear();
3724
3725 // Append the field into the fields vector.
3726 fields.push_back(field);
3727 }
3728 }
3729
3730 if(true == throw_error) {
3731 delete field;
3732 throw1(err_msg);
3733
3734 }
3735 }
3736 }
3737}
3738
3739// Retrieve attribute info.
3740void Dataset::ReadAttributes(int32 (*inq)(int32, char *, int32 *),
3741 intn (*attrinfo)(int32, char *, int32 *, int32 *),
3742 intn (*readattr)(int32, char *, VOIDP),
3743 vector<Attribute *> &obj_attrs) throw(Exception)
3744{
3745 // Initalize the number of attributes to be 0
3746 int32 numattrs = 0;
3747
3748 // Initalize the buffer size to be 0
3749 int32 bufsize = 0;
3750
3751 // Obtain the number of attributes in a Grid or Swath
3752 if ((numattrs = inq(this->datasetid, nullptr, &bufsize)) == -1)
3753 throw2("inquire attribute", this->name);
3754
3755 // Obtain the list of "name, type, value" tuple
3756 if (numattrs > 0) {
3757 vector<char> namelist;
3758
3759 namelist.resize(bufsize + 1);
3760
3761 // inquiry namelist and buffer size
3762 if (inq(this->datasetid, namelist.data(), &bufsize) == -1)
3763 throw2("inquire attribute", this->name);
3764
3765 vector<string> attrnames;
3766
3767 // Split the attribute namelist, make the "," separated name string to
3768 // a string list without ",".
3769 HDFCFUtil::Split(namelist.data(), bufsize, ',', attrnames);
3770 for (vector<string>::const_iterator i = attrnames.begin();
3771 i != attrnames.end(); ++i) {
3772 Attribute *attr = new Attribute();
3773 attr->name = *i;
3774 attr->newname = HDFCFUtil::get_CF_string(attr->name);
3775
3776 int32 count = 0;
3777
3778 // Obtain the datatype and byte count of this attribute
3779 if (attrinfo(this->datasetid, const_cast<char *>
3780 (attr->name.c_str()), &attr->type, &count) == -1) {
3781 delete attr;
3782 throw3("attribute info", this->name, attr->name);
3783 }
3784
3785 attr->count = count/DFKNTsize(attr->type);
3786 attr->value.resize(count);
3787
3788
3789 // Obtain the attribute value. Note that this function just
3790 // provides a copy of all attribute values.
3791 //The client should properly interpret to obtain individual
3792 // attribute value.
3793 if (readattr(this->datasetid,
3794 const_cast<char *>(attr->name.c_str()),
3795 &attr->value[0]) == -1) {
3796 delete attr;
3797 throw3("read attribute", this->name, attr->name);
3798 }
3799
3800 // Append this attribute to attrs list.
3801 obj_attrs.push_back(attr);
3802 }
3803 }
3804}
3805
3806// Release grid dataset resources
3807GridDataset::~GridDataset()
3808{
3809 if (this->datasetid != -1){
3810 GDdetach(this->datasetid);
3811 }
3812}
3813
3814// Retrieve all information of the grid dataset
3815GridDataset * GridDataset::Read(int32 fd, const string &gridname)
3816 throw(Exception)
3817{
3818 string err_msg;
3819 bool GD_fun_err = false;
3820 GridDataset *grid = new GridDataset(gridname);
3821
3822 // Open this Grid object
3823 if ((grid->datasetid = GDattach(fd, const_cast<char *>(gridname.c_str())))
3824 == -1) {
3825 err_msg = "attach grid";
3826 GD_fun_err = true;
3827 goto cleanFail;
3828 }
3829
3830 // Obtain the size of XDim and YDim as well as latitude and longitude of
3831 // the upleft corner and the low right corner.
3832 {
3833 Info &info = grid->info;
3834 if (GDgridinfo(grid->datasetid, &info.xdim, &info.ydim, info.upleft,
3835 info.lowright) == -1) {
3836 err_msg = "grid info";
3837 GD_fun_err = true;
3838 goto cleanFail;
3839 }
3840 }
3841
3842 // Obtain projection information.
3843 {
3844 Projection &proj = grid->proj;
3845 if (GDprojinfo(grid->datasetid, &proj.code, &proj.zone, &proj.sphere,
3846 proj.param) == -1) {
3847 err_msg = "projection info";
3848 GD_fun_err = true;
3849 goto cleanFail;
3850 }
3851 if (GDpixreginfo(grid->datasetid, &proj.pix) == -1) {
3852 err_msg = "pixreg info";
3853 GD_fun_err = true;
3854 goto cleanFail;
3855 }
3856 if (GDorigininfo(grid->datasetid, &proj.origin) == -1){
3857 err_msg = "origin info";
3858 GD_fun_err = true;
3859 goto cleanFail;
3860 }
3861 }
3862cleanFail:
3863 if(true == GD_fun_err){
3864 delete grid;
3865 throw2(err_msg,gridname);
3866 }
3867
3868 try {
3869 // Read dimensions
3870 grid->ReadDimensions(GDnentries, GDinqdims, grid->dims);
3871
3872 // Read all fields of this Grid.
3873 grid->ReadFields(GDnentries, GDinqfields, GDfieldinfo, GDreadfield,
3874 GDgetfillvalue, false, grid->datafields);
3875
3876 // Read all attributes of this Grid.
3877 grid->ReadAttributes(GDinqattrs, GDattrinfo, GDreadattr, grid->attrs);
3878 }
3879 catch (...) {
3880 delete grid;
3881 throw;
3882 }
3883
3884 return grid;
3885}
3886
3887GridDataset::Calculated & GridDataset::getCalculated() const
3888{
3889 if (this->calculated.grid != this)
3890 this->calculated.grid = this;
3891 return this->calculated;
3892}
3893
3894bool GridDataset::Calculated::isYDimMajor() throw(Exception)
3895{
3896 this->DetectMajorDimension();
3897 return this->ydimmajor;
3898}
3899
3900#if 0
3901bool GridDataset::Calculated::isOrthogonal() throw(Exception)
3902{
3903 if (!this->valid)
3904 this->ReadLongitudeLatitude();
3905 return this->orthogonal;
3906}
3907#endif
3908
3909int GridDataset::Calculated::DetectFieldMajorDimension() throw(Exception)
3910{
3911 int ym = -1;
3912
3913 // Traverse all data fields
3914 for (vector<Field *>::const_iterator i =
3915 this->grid->getDataFields().begin();
3916 i != this->grid->getDataFields().end(); ++i) {
3917
3918 int xdimindex = -1, ydimindex = -1, index = 0;
3919
3920 // Traverse all dimensions in each data field
3921 for (vector<Dimension *>::const_iterator j =
3922 (*i)->getDimensions().begin();
3923 j != (*i)->getDimensions().end(); ++j) {
3924 if ((*j)->getName() == this->grid->dimxname)
3925 xdimindex = index;
3926 else if ((*j)->getName() == this->grid->dimyname)
3927 ydimindex = index;
3928 ++index;
3929 }
3930 if (xdimindex == -1 || ydimindex == -1)
3931 continue;
3932
3933 int major = ydimindex < xdimindex ? 1 : 0;
3934
3935 if (ym == -1)
3936 ym = major;
3937 break;
3938
3939 // TO gain performance, just check one field.
3940 // The dimension order for all data fields in a grid should be
3941 // consistent.
3942 // Kent adds this if 0 block 2012-09-19
3943#if 0
3944 else if (ym != major)
3945 throw2("inconsistent XDim, YDim order", this->grid->getName());
3946#endif
3947
3948 }
3949 // At least one data field should refer to XDim or YDim
3950 if (ym == -1)
3951 throw2("lack of data fields", this->grid->getName());
3952
3953 return ym;
3954}
3955
3956void GridDataset::Calculated::DetectMajorDimension() throw(Exception)
3957{
3958 int ym = -1;
3959 // ydimmajor := true if (YDim, XDim)
3960 // ydimmajor := false if (XDim, YDim)
3961
3962 // Traverse all data fields
3963
3964 for (vector<Field *>::const_iterator i =
3965 this->grid->getDataFields().begin();
3966 i != this->grid->getDataFields().end(); ++i) {
3967
3968 int xdimindex = -1, ydimindex = -1, index = 0;
3969
3970 // Traverse all dimensions in each data field
3971 for (vector<Dimension *>::const_iterator j =
3972 (*i)->getDimensions().begin();
3973 j != (*i)->getDimensions().end(); ++j) {
3974 if ((*j)->getName() == this->grid->dimxname)
3975 xdimindex = index;
3976 else if ((*j)->getName() == this->grid->dimyname)
3977 ydimindex = index;
3978 ++index;
3979 }
3980 if (xdimindex == -1 || ydimindex == -1)
3981 continue;
3982
3983 // Change the way of deciding the major dimesion (LD -2012/01/10).
3984 int major;
3985 if(this->grid->getProjection().getCode() == GCTP_LAMAZ)
3986 major = 1;
3987 else
3988 major = ydimindex < xdimindex ? 1 : 0;
3989
3990 if (ym == -1)
3991 ym = major;
3992 break;
3993
3994 // TO gain performance, just check one field.
3995 // The dimension order for all data fields in a grid should be
3996 // consistent.
3997 // Kent adds this if 0 block
3998#if 0
3999 else if (ym != major)
4000 throw2("inconsistent XDim, YDim order", this->grid->getName());
4001#endif
4002 }
4003 // At least one data field should refer to XDim or YDim
4004 if (ym == -1)
4005 throw2("lack of data fields", this->grid->getName());
4006 this->ydimmajor = ym != 0;
4007}
4008
4009// The following internal utilities are not used currently, will see if
4010// they are necessary in the future. KY 2012-09-19
4011// The internal utility method to check if two vectors have overlapped.
4012// If not, return true.
4013// Not used. Temporarily comment out to avoid the compiler warnings.
4014#if 0
4015static bool IsDisjoint(const vector<Field *> &l,
4016 const vector<Field *> &r)
4017{
4018 for (vector<Field *>::const_iterator i = l.begin(); i != l.end(); ++i)
4019 {
4020 if (find(r.begin(), r.end(), *i) != r.end())
4021 return false;
4022 }
4023 return true;
4024}
4025
4026// The internal utility method to check if two vectors have overlapped.
4027// If not, return true.
4028static bool IsDisjoint(vector<pair<Field *, string> > &l, const vector<Field *> &r)
4029{
4030 for (vector<pair<Field *, string> >::const_iterator i =
4031 l.begin(); i != l.end(); ++i) {
4032 if (find(r.begin(), r.end(), i->first) != r.end())
4033 return false;
4034 }
4035 return true;
4036}
4037
4038// The internal utility method to check if vector s is a subset of vector b.
4039static bool IsSubset(const vector<Field *> &s, const vector<Field *> &b)
4040{
4041 for (vector<Field *>::const_iterator i = s.begin(); i != s.end(); ++i)
4042 {
4043 if (find(b.begin(), b.end(), *i) == b.end())
4044 return false;
4045 }
4046 return true;
4047}
4048
4049// The internal utility method to check if vector s is a subset of vector b.
4050static bool IsSubset(vector<pair<Field *, string> > &s, const vector<Field *> &b)
4051{
4052 for (vector<pair<Field *, string> >::const_iterator i
4053 = s.begin(); i != s.end(); ++i) {
4054 if (find(b.begin(), b.end(), i->first) == b.end())
4055 return false;
4056 }
4057 return true;
4058}
4059
4060#endif
4061// Destructor, release resources
4062SwathDataset::~SwathDataset()
4063{
4064 if (this->datasetid != -1) {
4065 SWdetach(this->datasetid);
4066 }
4067
4068 for_each(this->dimmaps.begin(), this->dimmaps.end(), delete_elem());
4069 for_each(this->indexmaps.begin(), this->indexmaps.end(),
4070 delete_elem());
4071
4072 for_each(this->geofields.begin(), this->geofields.end(),
4073 delete_elem());
4074 return;
4075
4076}
4077
4078// Read all information of this swath
4079SwathDataset * SwathDataset::Read(int32 fd, const string &swathname)
4080 throw(Exception)
4081{
4082 SwathDataset *swath = new SwathDataset(swathname);
4083 if(swath == nullptr)
4084 throw1("Cannot allocate HDF5 Swath object");
4085
4086 // Open this Swath object
4087 if ((swath->datasetid = SWattach(fd,
4088 const_cast<char *>(swathname.c_str())))
4089 == -1) {
4090 delete swath;
4091 throw2("attach swath", swathname);
4092 }
4093
4094 //if(swath != nullptr) {// See if I can make coverity happy.coverity doesn't know I call throw already.
4095
4096 try {
4097
4098 // Read dimensions of this Swath
4099 swath->ReadDimensions(SWnentries, SWinqdims, swath->dims);
4100
4101 // Read all information related to data fields of this Swath
4102 swath->ReadFields(SWnentries, SWinqdatafields, SWfieldinfo, SWreadfield,
4103 SWgetfillvalue, false, swath->datafields);
4104
4105 // Read all information related to geo-location fields of this Swath
4106 swath->ReadFields(SWnentries, SWinqgeofields, SWfieldinfo, SWreadfield,
4107 SWgetfillvalue, true, swath->geofields);
4108
4109 // Read all attributes of this Swath
4110 swath->ReadAttributes(SWinqattrs, SWattrinfo, SWreadattr, swath->attrs);
4111
4112 // Read dimension map and save the number of dimension map for dim. subsetting
4113 swath->set_num_map(swath->ReadDimensionMaps(swath->dimmaps));
4114
4115 // Read index maps, we haven't found any files with the Index Maps.
4116 swath->ReadIndexMaps(swath->indexmaps);
4117 }
4118 catch (...) {
4119 delete swath;
4120 throw;
4121 }
4122 //}
4123
4124 return swath;
4125}
4126
4127// Read dimension map info.
4128int SwathDataset::ReadDimensionMaps(vector<DimensionMap *> &swath_dimmaps)
4129 throw(Exception)
4130{
4131 int32 nummaps, bufsize;
4132
4133 // Obtain number of dimension maps and the buffer size.
4134 if ((nummaps = SWnentries(this->datasetid, HDFE_NENTMAP, &bufsize)) == -1)
4135 throw2("dimmap entry", this->name);
4136
4137 // Will use Split utility method to generate a list of dimension map.
4138 // An example of original representation of a swath dimension map:(d1/d2,
4139 // d3/d4,...)
4140 // d1:the name of this dimension of the data field , d2: the name of the
4141 // corresponding dimension of the geo-location field
4142 // The list will be decomposed from (d1/d2,d3/d4,...) to {[d1,d2,0(offset),
4143 // 1(increment)],[d3,d4,0(offset),1(increment)],...}
4144
4145 if (nummaps > 0) {
4146 vector<char> namelist;
4147 vector<int32> offset, increment;
4148
4149 namelist.resize(bufsize + 1);
4150 offset.resize(nummaps);
4151 increment.resize(nummaps);
4152 if (SWinqmaps(this->datasetid, namelist.data(), offset.data(), increment.data())
4153 == -1)
4154 throw2("inquire dimmap", this->name);
4155
4156 vector<string> mapnames;
4157 HDFCFUtil::Split(namelist.data(), bufsize, ',', mapnames);
4158 int count = 0;
4159 for (vector<string>::const_iterator i = mapnames.begin();
4160 i != mapnames.end(); ++i) {
4161 vector<string> parts;
4162 HDFCFUtil::Split(i->c_str(), '/', parts);
4163 if (parts.size() != 2)
4164 throw3("dimmap part", parts.size(),
4165 this->name);
4166
4167 DimensionMap *dimmap = new DimensionMap(parts[0], parts[1],
4168 offset[count],
4169 increment[count]);
4170 swath_dimmaps.push_back(dimmap);
4171 ++count;
4172 }
4173 }
4174 return nummaps;
4175}
4176
4177// The following function is nevered tested and probably will never be used.
4178void SwathDataset::ReadIndexMaps(vector<IndexMap *> &swath_indexmaps)
4179 throw(Exception)
4180{
4181 int32 numindices, bufsize;
4182
4183 if ((numindices = SWnentries(this->datasetid, HDFE_NENTIMAP, &bufsize)) ==
4184 -1)
4185 throw2("indexmap entry", this->name);
4186 if (numindices > 0) {
4187 // TODO: I have never seen any EOS2 files that have index map.
4188 vector<char> namelist;
4189
4190 namelist.resize(bufsize + 1);
4191 if (SWinqidxmaps(this->datasetid, namelist.data(), nullptr) == -1)
4192 throw2("inquire indexmap", this->name);
4193
4194 vector<string> mapnames;
4195 HDFCFUtil::Split(namelist.data(), bufsize, ',', mapnames);
4196 for (vector<string>::const_iterator i = mapnames.begin();
4197 i != mapnames.end(); ++i) {
4198 IndexMap *indexmap = new IndexMap();
4199 vector<string> parts;
4200 HDFCFUtil::Split(i->c_str(), '/', parts);
4201 if (parts.size() != 2)
4202 throw3("indexmap part", parts.size(),
4203 this->name);
4204 indexmap->geo = parts[0];
4205 indexmap->data = parts[1];
4206
4207 {
4208 int32 dimsize;
4209 if ((dimsize =
4210 SWdiminfo(this->datasetid,
4211 const_cast<char *>(indexmap->geo.c_str())))
4212 == -1)
4213 throw3("dimension info", this->name, indexmap->geo);
4214 indexmap->indices.resize(dimsize);
4215 if (SWidxmapinfo(this->datasetid,
4216 const_cast<char *>(indexmap->geo.c_str()),
4217 const_cast<char *>(indexmap->data.c_str()),
4218 &indexmap->indices[0]) == -1)
4219 throw4("indexmap info", this->name, indexmap->geo,
4220 indexmap->data);
4221 }
4222
4223 swath_indexmaps.push_back(indexmap);
4224 }
4225 }
4226}
4227
4228
4229PointDataset::~PointDataset()
4230{
4231}
4232
4233PointDataset * PointDataset::Read(int32 /*fd*/, const string &pointname)
4234 throw(Exception)
4235{
4236 PointDataset *point = new PointDataset(pointname);
4237 return point;
4238}
4239
4240
4241// Read name list from the EOS2 APIs and store names into a vector
4242bool Utility::ReadNamelist(const char *path,
4243 int32 (*inq)(char *, char *, int32 *),
4244 vector<string> &names)
4245{
4246 char *fname = const_cast<char *>(path);
4247 int32 bufsize;
4248 int numobjs;
4249
4250 if ((numobjs = inq(fname, nullptr, &bufsize)) == -1) return false;
4251 if (numobjs > 0) {
4252 vector<char> buffer;
4253 buffer.resize(bufsize + 1);
4254 if (inq(fname, buffer.data(), &bufsize) == -1) return false;
4255 HDFCFUtil::Split(buffer.data(), bufsize, ',', names);
4256 }
4257 return true;
4258}
4259#endif
4260
4261
static bool insert_map(std::map< std::string, std::string > &m, std::string key, std::string val)
Definition: HDFCFUtil.cc:150
static void Handle_NameClashing(std::vector< std::string > &newobjnamelist)
General routines to handle name clashings.
Definition: HDFCFUtil.cc:260
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