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