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