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