21 #include "HDFCFUtil.h"
23 #include "HDF4RequestHandler.h"
33 const char *HDFEOS2::File::_geodim_x_names[] = {
"XDim",
"LonDim",
"nlon"};
36 const char *HDFEOS2::File::_geodim_y_names[] = {
"YDim",
"LatDim",
"nlat"};
39 const char *HDFEOS2::File::_latfield_names[] = {
40 "Latitude",
"Lat",
"YDim",
"LatCenter"
44 const char *HDFEOS2::File::_lonfield_names[] = {
45 "Longitude",
"Lon",
"XDim",
"LonCenter"
50 const char *HDFEOS2::File::_geogrid_names[] = {
"location"};
52 using namespace HDFEOS2;
56 template<
typename T,
typename U,
typename V,
typename W,
typename X>
57 static void _throw5(
const char *fname,
int line,
int numarg,
58 const T &a1,
const U &a2,
const V &a3,
const W &a4,
62 ss << fname <<
":" << line <<
":";
63 for (
int i = 0; i < numarg; ++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;
72 ss<<
" Argument number is beyond 5 ";
75 throw Exception(ss.str());
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)
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))
93 template<
typename T>
void operator()(T *ptr)
103 for (vector<GridDataset *>::const_iterator i = grids.begin();
104 i != grids.end(); ++i){
111 for (vector<SwathDataset *>::const_iterator i = swaths.begin();
112 i != swaths.end(); ++i){
118 for (vector<PointDataset *>::const_iterator i = points.begin();
119 i != points.end(); ++i){
126 File * File::Read(
const char *path, int32 mygridfd, int32 myswathfd)
throw(Exception)
129 File *file =
new File(path);
131 throw1(
"Memory allocation for file class failed. ");
133 file->gridfd = mygridfd;
134 file->swathfd = myswathfd;
138 if ((file->gridfd = GDopen(
const_cast<char *
>(file->path.c_str()),
139 DFACC_READ)) == -1) {
141 throw2(
"grid open", path);
145 vector<string> gridlist;
146 if (!Utility::ReadNamelist(file->path.c_str(), GDinqgrid, gridlist)) {
148 throw1(
"Grid ReadNamelist failed.");
152 for (vector<string>::const_iterator i = gridlist.begin();
153 i != gridlist.end(); ++i)
154 file->grids.push_back(GridDataset::Read(file->gridfd, *i));
158 throw1(
"GridDataset Read failed");
163 if ((file->swathfd = SWopen(
const_cast<char *
>(file->path.c_str()),
166 throw2(
"swath open", path);
170 vector<string> swathlist;
171 if (!Utility::ReadNamelist(file->path.c_str(), SWinqswath, swathlist)){
173 throw1(
"Swath ReadNamelist failed.");
177 for (vector<string>::const_iterator i = swathlist.begin();
178 i != swathlist.end(); ++i)
179 file->swaths.push_back(SwathDataset::Read(file->swathfd, *i));
183 throw1(
"SwathDataset Read failed.");
190 vector<string> pointlist;
191 if (!Utility::ReadNamelist(file->path.c_str(), PTinqpoint, pointlist)){
193 throw1(
"Point ReadNamelist failed.");
196 for (vector<string>::const_iterator i = pointlist.begin();
197 i != pointlist.end(); ++i)
198 file->points.push_back(PointDataset::Read(-1, *i));
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);
215 string File::get_geodim_x_name()
217 if(!_geodim_x_name.empty())
218 return _geodim_x_name;
219 _find_geodim_names();
220 return _geodim_x_name;
226 string File::get_geodim_y_name()
228 if(!_geodim_y_name.empty())
229 return _geodim_y_name;
230 _find_geodim_names();
231 return _geodim_y_name;
241 string File::get_latfield_name()
243 if(!_latfield_name.empty())
244 return _latfield_name;
245 _find_latlonfield_names();
246 return _latfield_name;
249 string File::get_lonfield_name()
251 if(!_lonfield_name.empty())
252 return _lonfield_name;
253 _find_latlonfield_names();
254 return _lonfield_name;
261 string File::get_geogrid_name()
263 if(!_geogrid_name.empty())
264 return _geogrid_name;
265 _find_geogrid_name();
266 return _geogrid_name;
274 void File::_find_geodim_names()
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]);
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]);
287 const size_t gs = grids.size();
288 const size_t ss = swaths.size();
289 for(
size_t i=0; ;i++)
291 Dataset *dataset=NULL;
293 dataset =
static_cast<Dataset*
>(grids[i]);
295 dataset =
static_cast<Dataset*
>(swaths[i-gs]);
299 const vector<Dimension *>& dims = dataset->getDimensions();
300 for(vector<Dimension*>::const_iterator it = dims.begin();
301 it != dims.end(); ++it)
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();
313 const size_t gs = grids.size();
317 Dataset *dataset=NULL;
318 dataset =
static_cast<Dataset*
>(grids[0]);
320 const vector<Dimension *>& dims = dataset->getDimensions();
321 for(vector<Dimension*>::const_iterator it = dims.begin();
322 it != dims.end(); ++it)
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();
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];
347 void File::_find_latlonfield_names()
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]);
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]);
357 const size_t gs = grids.size();
358 const size_t ss = swaths.size();
362 for(
size_t i=0;i<1 ;i++)
364 Dataset *dataset = NULL;
365 SwathDataset *sw = NULL;
367 dataset =
static_cast<Dataset*
>(grids[i]);
371 dataset =
static_cast<Dataset*
>(sw);
376 const vector<Field *>& fields = dataset->getDataFields();
377 for(vector<Field*>::const_iterator it = fields.begin();
378 it != fields.end(); ++it)
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();
388 const vector<Field *>& geofields = dataset->getDataFields();
389 for(vector<Field*>::const_iterator it = geofields.begin();
390 it != geofields.end(); ++it)
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();
402 if(_latfield_name.empty())
403 _latfield_name = _latfield_names[0];
404 if(_lonfield_name.empty())
405 _lonfield_name = _lonfield_names[0];
413 void File::_find_geogrid_name()
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]);
419 const size_t gs = grids.size();
420 const size_t ss = swaths.size();
421 for(
size_t i=0; ;i++)
425 dataset =
static_cast<Dataset*
>(grids[i]);
427 dataset =
static_cast<Dataset*
>(swaths[i-gs]);
431 if(geogrid_name_set.find(dataset->getName()) != geogrid_name_set.end())
432 _geogrid_name = dataset->getName();
434 if(_geogrid_name.empty())
435 _geogrid_name =
"location";
439 void File::check_onelatlon_grids() {
442 string LATFIELDNAME = this->get_latfield_name();
443 string LONFIELDNAME = this->get_lonfield_name();
446 string GEOGRIDNAME =
"location";
455 for(vector<GridDataset *>::const_iterator i = this->grids.begin();
456 i != this->grids.end(); ++i){
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){
467 if((*j)->getName()==LONFIELDNAME){
475 if(((*j)->getName()==LATFIELDNAME)||((*j)->getName()==LONFIELDNAME)){
476 (*i)->ownllflag =
true;
484 if(morellcount ==0 && onellcount ==2)
485 this->onelatlon =
true;
489 void File::handle_one_grid_zdim(GridDataset* gdset) {
492 string DIMXNAME = this->get_geodim_x_name();
493 string DIMYNAME = this->get_geodim_y_name();
495 bool missingfield_unlim_flag =
false;
496 Field *missingfield_unlim = NULL;
503 set<string> tempdimlist;
504 pair<set<string>::iterator,
bool> tempdimret;
506 for (vector<Field *>::const_iterator j =
507 gdset->getDataFields().begin();
508 j != gdset->getDataFields().end(); ++j) {
511 if ((*j)->getRank()==1){
515 if(((*j)->getDimensions())[0]->getName()!=DIMXNAME && ((*j)->getDimensions())[0]->getName()!=DIMYNAME){
517 tempdimret = tempdimlist.insert(((*j)->getDimensions())[0]->getName());
522 if(tempdimret.second ==
true) {
527 if((*j)->getName() ==
"Time")
543 for (vector<Dimension *>::const_iterator j =
544 gdset->getDimensions().begin(); j!= gdset->getDimensions().end();++j){
547 if((*j)->getName()!=DIMXNAME && (*j)->getName()!=DIMYNAME){
550 if((tempdimlist.find((*j)->getName())) == tempdimlist.end()){
553 Field *missingfield =
new Field();
554 missingfield->name = (*j)->getName();
555 missingfield->rank = 1;
558 missingfield->type = DFNT_INT32;
561 Dimension *dim =
new Dimension((*j)->getName(),(*j)->getSize());
564 missingfield->dims.push_back(dim);
574 if(0 == (*j)->getSize()) {
575 missingfield_unlim_flag =
true;
576 missingfield_unlim = missingfield;
581 missingfield->fieldtype = 4;
591 gdset->datafields.push_back(missingfield);
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) {
606 for (vector<Dimension *>::const_iterator j =
607 gdset->getDimensions().begin(); j!= gdset->getDimensions().end();++j){
609 if((*j)->getName() == (missingfield_unlim->getDimensions())[0]->getName()) {
610 if((*j)->getSize()!= 0) {
611 Dimension *dim = missingfield_unlim->getDimensions()[0];
614 dim->dimsize = (*j)->getSize();
615 missingfield_unlim_flag =
false;
621 if(
false == missingfield_unlim_flag)
629 void File::handle_one_grid_latlon(GridDataset* gdset)
throw(Exception)
633 string DIMXNAME = this->get_geodim_x_name();
634 string DIMYNAME = this->get_geodim_y_name();
636 string LATFIELDNAME = this->get_latfield_name();
637 string LONFIELDNAME = this->get_lonfield_name();
641 if(gdset->ownllflag) {
644 for (vector<Field *>::const_iterator j =
645 gdset->getDataFields().begin();
646 j != gdset->getDataFields().end(); ++j) {
648 if((*j)->getName() == LATFIELDNAME) {
659 if((*j)->getRank() > 2)
660 throw3(
"The rank of latitude is greater than 2",
661 gdset->getName(),(*j)->getName());
663 if((*j)->getRank() != 1) {
667 (*j)->ydimmajor = gdset->getCalculated().isYDimMajor();
673 int32 projectioncode = gdset->getProjection().getCode();
674 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
675 (*j)->condenseddim =
true;
685 if((*j)->condenseddim) {
689 for (vector<Dimension *>::const_iterator k =
690 (*j)->getDimensions().begin(); k!= (*j)->getDimensions().end();++k){
691 if((*k)->getName() == DIMYNAME) {
699 for (vector<Dimension *>::const_iterator k =
700 (*j)->getDimensions().begin(); k!= (*j)->getDimensions().end();++k){
701 if((*k)->getName() == DIMYNAME) {
713 else if ((*j)->getName() == LONFIELDNAME) {
722 if((*j)->getRank() >2)
723 throw3(
"The rank of Longitude is greater than 2",gdset->getName(),(*j)->getName());
725 if((*j)->getRank() != 1) {
728 (*j)->ydimmajor = gdset->getCalculated().isYDimMajor();
733 int32 projectioncode = gdset->getProjection().getCode();
734 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
735 (*j)->condenseddim =
true;
744 if((*j)->condenseddim) {
748 for (vector<Dimension *>::const_iterator k =
749 (*j)->getDimensions().begin(); k!= (*j)->getDimensions().end();++k){
750 if((*k)->getName() == DIMXNAME) {
757 for (vector<Dimension *>::const_iterator k =
758 (*j)->getDimensions().begin(); k!= (*j)->getDimensions().end();++k){
759 if((*k)->getName() == DIMXNAME) {
776 Field *latfield =
new Field();
777 Field *lonfield =
new Field();
779 latfield->name = LATFIELDNAME;
780 lonfield->name = LONFIELDNAME;
787 latfield->type = DFNT_FLOAT64;
788 lonfield->type = DFNT_FLOAT64;
791 latfield->fieldtype = 1;
794 lonfield->fieldtype = 2;
798 latfield->ydimmajor = gdset->getCalculated().isYDimMajor();
799 lonfield->ydimmajor = latfield->ydimmajor;
802 int xdimsize = gdset->getInfo().getX();
803 int ydimsize = gdset->getInfo().getY();
808 bool dmajor=(gdset->getProjection().getCode()==GCTP_LAMAZ)? gdset->getCalculated().DetectFieldMajorDimension()
809 : latfield->ydimmajor;
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);
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);
835 upleft =
const_cast<float64 *
>(gdset->getInfo().getUpLeft());
836 lowright =
const_cast<float64 *
>(gdset->getInfo().getLowRight());
839 int32 projectioncode = gdset->getProjection().getCode();
840 if(((
int)lowright[0]>180000000) && ((
int)upleft[0]>-1)) {
843 if(projectioncode == GCTP_GEO) {
844 lonfield->speciallon =
true;
850 if(HDF4RequestHandler::get_enable_eosgeo_cachefile() ==
true)
851 latfield->speciallon =
true;
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;
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;
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;
902 if(GCTP_SOM == projectioncode) {
903 lonfield->specialformat = 4;
904 latfield->specialformat = 4;
910 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
911 lonfield->condenseddim =
true;
912 latfield->condenseddim =
true;
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) {
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());
936 cache_fname +=HDFCFUtil::get_int_str(ydimsize);
937 cache_fname +=HDFCFUtil::get_int_str(xdimsize);
940 cache_fname +=HDFCFUtil::get_int_str(xdimsize);
941 cache_fname +=HDFCFUtil::get_int_str(ydimsize);
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);
954 params =
const_cast<float64 *
>(gdset->getProjection().getParam());
957 for(
int ipar = 0; ipar<13;ipar++)
958 cache_fname+=HDFCFUtil::get_double_str(params[ipar],17,6);
967 string cache_fpath =
"/tmp/"+cache_fname;
968 int result = stat(cache_fpath.c_str(), &st);
970 int actual_file_size = st.st_size;
971 cerr<<
"HDF-EOS2 actual_file_size is "<<actual_file_size <<endl;
972 int expected_file_size = 0;
973 if(gdset->getProjection().getCode() == GCTP_SOM)
974 expected_file_size = xdimsize*ydimsize*2*
sizeof(double)*NBLOCK;
975 else if(gdset->getProjection().getCode() == GCTP_CEA ||
976 gdset->getProjection().getCode() == GCTP_GEO)
977 expected_file_size = (xdimsize+ydimsize)*
sizeof(double);
979 expected_file_size = xdimsize*ydimsize*2*
sizeof(double);
981 cerr<<
"HDF-EOS2 expected_file_size is "<<expected_file_size <<endl;
982 if(actual_file_size != expected_file_size){
983 cerr<<
"field_cache is 1 "<<endl;
984 lonfield->field_cache = 1;
985 latfield->field_cache = 1;
988 cerr<<
"field cache is 2 "<<endl;
989 lonfield->field_cache = 2;
990 latfield->field_cache = 2;
1006 gdset->datafields.push_back(latfield);
1007 gdset->datafields.push_back(lonfield);
1015 if(lonfield->condenseddim) {
1019 for (vector<Dimension *>::const_iterator j =
1020 lonfield->getDimensions().begin(); j!= lonfield->getDimensions().end();++j){
1022 if((*j)->getName() == DIMXNAME) {
1026 if((*j)->getName() == DIMYNAME) {
1032 for (vector<Dimension *>::const_iterator j =
1033 lonfield->getDimensions().begin(); j!= lonfield->getDimensions().end();++j){
1035 if((*j)->getName() == DIMXNAME){
1039 if((*j)->getName() == DIMYNAME){
1050 void File::handle_onelatlon_grids() throw(Exception) {
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();
1060 string GEOGRIDNAME =
"location";
1063 map<string,string>temponelatlondimcvarlist;
1066 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1067 i != this->grids.end(); ++i){
1071 (*i)->setDimxName(DIMXNAME);
1072 (*i)->setDimyName(DIMYNAME);
1075 if((*i)->getName()==GEOGRIDNAME) {
1080 (*i)->lonfield->fieldtype = 2;
1081 (*i)->latfield->fieldtype = 1;
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());
1090 if((*i)->lonfield->rank != 1) {
1094 (*i)->lonfield->ydimmajor = (*i)->getCalculated().isYDimMajor();
1095 (*i)->latfield->ydimmajor = (*i)->lonfield->ydimmajor;
1100 int32 projectioncode = (*i)->getProjection().getCode();
1101 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
1102 (*i)->lonfield->condenseddim =
true;
1103 (*i)->latfield->condenseddim =
true;
1110 if((*i)->lonfield->condenseddim) {
1116 for (vector<Dimension *>::const_iterator j =
1117 (*i)->lonfield->getDimensions().begin(); j!= (*i)->lonfield->getDimensions().end();++j){
1118 if((*j)->getName() == DIMXNAME) {
1120 (*i)->lonfield->getName());
1122 if((*j)->getName() == DIMYNAME) {
1124 (*i)->latfield->getName());
1131 for (vector<Dimension *>::const_iterator j =
1132 (*i)->lonfield->getDimensions().begin(); j!= (*i)->lonfield->getDimensions().end();++j){
1133 if((*j)->getName() == DIMXNAME) {
1135 (*i)->lonfield->getName());
1137 if((*j)->getName() == DIMYNAME) {
1139 (*i)->latfield->getName());
1147 (*i)->lonfield->getName());
1149 (*i)->latfield->getName());
1151 temponelatlondimcvarlist = (*i)->dimcvarlist;
1159 for(vector<GridDataset *>::const_iterator i = this->grids.begin();
1160 i != this->grids.end(); ++i){
1162 string templatlonname1;
1163 string templatlonname2;
1165 if((*i)->getName() != GEOGRIDNAME) {
1167 map<string,string>::iterator tempmapit;
1170 tempmapit = temponelatlondimcvarlist.find(DIMXNAME);
1171 if(tempmapit != temponelatlondimcvarlist.end())
1172 templatlonname1= tempmapit->second;
1174 throw2(
"cannot find the dimension field of XDim", (*i)->getName());
1179 tempmapit = temponelatlondimcvarlist.find(DIMYNAME);
1180 if(tempmapit != temponelatlondimcvarlist.end())
1181 templatlonname2= tempmapit->second;
1183 throw2(
"cannot find the dimension field of YDim", (*i)->getName());
1191 void File::handle_grid_dim_cvar_maps() throw(Exception) {
1194 string DIMXNAME = this->get_geodim_x_name();
1196 string DIMYNAME = this->get_geodim_y_name();
1198 string LATFIELDNAME = this->get_latfield_name();
1200 string LONFIELDNAME = this->get_lonfield_name();
1205 string GEOGRIDNAME =
"location";
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) {
1233 map<string,string>tempncvarnamelist;
1234 string tempcorrectedlatname, tempcorrectedlonname;
1236 int total_fcounter = 0;
1238 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1239 i != this->grids.end(); ++i){
1244 for (vector<Field *>::const_iterator j =
1245 (*i)->getDataFields().begin();
1246 j != (*i)->getDataFields().end(); ++j)
1248 (*j)->newname = tempfieldnamelist[total_fcounter];
1252 if((*j)->fieldtype!=0) {
1254 tempncvarnamelist.insert(make_pair((*j)->getName(), (*j)->newname));
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;
1266 (*i)->ncvarnamelist = tempncvarnamelist;
1267 tempncvarnamelist.clear();
1273 if(this->onelatlon) {
1274 for(vector<GridDataset *>::const_iterator i = this->grids.begin();
1275 i != this->grids.end(); ++i){
1277 if((*i)->getName()!=GEOGRIDNAME){
1285 map<string,string>tempndimnamelist;
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)
1299 int total_dcounter = 0;
1300 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1301 i != this->grids.end(); ++i){
1303 for (map<string,string>::const_iterator j =
1304 (*i)->dimcvarlist.begin(); j!= (*i)->dimcvarlist.end();++j){
1307 if((DIMXNAME == (*j).first || DIMYNAME == (*j).first) && (
true==(this->onelatlon)))
1314 (*i)->ndimnamelist = tempndimnamelist;
1315 tempndimnamelist.clear();
1320 void File::handle_grid_coards() throw(Exception) {
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();
1330 string GEOGRIDNAME =
"location";
1334 vector<Dimension*> correcteddims;
1335 string tempcorrecteddimname;
1338 map<string,string> tempnewxdimnamelist;
1341 map<string,string> tempnewydimnamelist;
1344 Dimension *correcteddim;
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) {
1355 if((*j)->getName()==LATFIELDNAME && (*j)->getRank()==2 &&(*j)->condenseddim) {
1357 string templatdimname;
1358 map<string,string>::iterator tempmapit;
1361 tempmapit = (*i)->ncvarnamelist.find(LATFIELDNAME);
1362 if(tempmapit != (*i)->ncvarnamelist.end())
1363 templatdimname= tempmapit->second;
1365 throw2(
"cannot find the corrected field of Latitude", (*i)->getName());
1367 for(vector<Dimension *>::const_iterator k =(*j)->getDimensions().begin();
1368 k!=(*j)->getDimensions().end();++k){
1372 if((*k)->getName()==DIMYNAME) {
1373 correcteddim =
new Dimension(templatdimname,(*k)->getSize());
1374 correcteddims.push_back(correcteddim);
1375 (*j)->setCorrectedDimensions(correcteddims);
1379 (*j)->iscoard =
true;
1380 (*i)->iscoard =
true;
1382 this->iscoard =
true;
1385 correcteddims.clear();
1389 else if((*j)->getName()==LONFIELDNAME && (*j)->getRank()==2 &&(*j)->condenseddim){
1391 string templondimname;
1392 map<string,string>::iterator tempmapit;
1395 tempmapit = (*i)->ncvarnamelist.find(LONFIELDNAME);
1396 if(tempmapit != (*i)->ncvarnamelist.end())
1397 templondimname= tempmapit->second;
1399 throw2(
"cannot find the corrected field of Longitude", (*i)->getName());
1401 for(vector<Dimension *>::const_iterator k =(*j)->getDimensions().begin();
1402 k!=(*j)->getDimensions().end();++k){
1406 if((*k)->getName()==DIMXNAME) {
1407 correcteddim =
new Dimension(templondimname,(*k)->getSize());
1408 correcteddims.push_back(correcteddim);
1409 (*j)->setCorrectedDimensions(correcteddims);
1414 (*j)->iscoard =
true;
1415 (*i)->iscoard =
true;
1417 this->iscoard =
true;
1418 correcteddims.clear();
1422 else if(((*j)->getRank()==1) &&((*j)->getName()==LONFIELDNAME) ) {
1424 string templondimname;
1425 map<string,string>::iterator tempmapit;
1428 tempmapit = (*i)->ncvarnamelist.find(LONFIELDNAME);
1429 if(tempmapit != (*i)->ncvarnamelist.end())
1430 templondimname= tempmapit->second;
1432 throw2(
"cannot find the corrected field of Longitude", (*i)->getName());
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;
1440 this->iscoard =
true;
1441 correcteddims.clear();
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());
1448 if((((*j)->getDimensions())[0]->getName())==DIMXNAME) {
1457 else if(((*j)->getRank()==1) &&((*j)->getName()==LATFIELDNAME) ) {
1459 string templatdimname;
1460 map<string,string>::iterator tempmapit;
1463 tempmapit = (*i)->ncvarnamelist.find(LATFIELDNAME);
1464 if(tempmapit != (*i)->ncvarnamelist.end())
1465 templatdimname= tempmapit->second;
1467 throw2(
"cannot find the corrected field of Latitude", (*i)->getName());
1469 correcteddim =
new Dimension(templatdimname,((*j)->getDimensions())[0]->getSize());
1470 correcteddims.push_back(correcteddim);
1471 (*j)->setCorrectedDimensions(correcteddims);
1473 (*j)->iscoard =
true;
1474 (*i)->iscoard =
true;
1476 this->iscoard =
true;
1477 correcteddims.clear();
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){
1497 if(
true == this->onelatlon){
1500 if(
true == this->iscoard){
1503 string tempcorrectedxdimname;
1504 string tempcorrectedydimname;
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");
1511 map<string,string>::iterator tempdimmapit = tempnewxdimnamelist.begin();
1512 tempcorrectedxdimname = tempdimmapit->second;
1513 tempdimmapit = tempnewydimnamelist.begin();
1514 tempcorrectedydimname = tempdimmapit->second;
1516 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1517 i != this->grids.end(); ++i){
1520 map<string,string>::iterator tempmapit;
1521 tempmapit = (*i)->ndimnamelist.find(DIMXNAME);
1522 if(tempmapit != (*i)->ndimnamelist.end()) {
1526 throw2(
"cannot find the corrected dimension name", (*i)->getName());
1527 tempmapit = (*i)->ndimnamelist.find(DIMYNAME);
1528 if(tempmapit != (*i)->ndimnamelist.end()) {
1532 throw2(
"cannot find the corrected dimension name", (*i)->getName());
1537 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1538 i != this->grids.end(); ++i){
1541 string tempcorrectedxdimname;
1542 string tempcorrectedydimname;
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;
1551 throw2(
"cannot find the corrected COARD XDim dimension name", (*i)->getName());
1552 tempmapit = (*i)->ndimnamelist.find(DIMXNAME);
1553 if(tempmapit != (*i)->ndimnamelist.end()) {
1557 throw2(
"cannot find the corrected dimension name", (*i)->getName());
1559 tempdimmapit = tempnewydimnamelist.find((*i)->getName());
1560 if(tempdimmapit != tempnewydimnamelist.end())
1561 tempcorrectedydimname = tempdimmapit->second;
1563 throw2(
"cannot find the corrected COARD YDim dimension name", (*i)->getName());
1565 tempmapit = (*i)->ndimnamelist.find(DIMYNAME);
1566 if(tempmapit != (*i)->ndimnamelist.end()) {
1570 throw2(
"cannot find the corrected dimension name", (*i)->getName());
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){
1585 if((this->iscoard||(*i)->iscoard) && (*j).first !=DIMXNAME && (*j).first !=DIMYNAME) {
1586 string tempnewdimname;
1587 map<string,string>::iterator tempmapit;
1590 tempmapit = (*i)->ncvarnamelist.find((*j).second);
1591 if(tempmapit != (*i)->ncvarnamelist.end())
1592 tempnewdimname= tempmapit->second;
1594 throw3(
"cannot find the corrected field of ", (*j).second,(*i)->getName());
1597 tempmapit =(*i)->ndimnamelist.find((*j).first);
1598 if(tempmapit != (*i)->ndimnamelist.end())
1601 throw3(
"cannot find the corrected dimension name of ", (*j).first,(*i)->getName());
1608 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1609 i != this->grids.end(); ++i){
1611 for (vector<Field *>::const_iterator j =
1612 (*i)->getDataFields().begin();
1613 j != (*i)->getDataFields().end(); ++j) {
1616 if((*j)->iscoard ==
false) {
1619 for(vector<Dimension *>::const_iterator k=(*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k){
1622 map<string,string>::iterator tempmapit;
1625 tempmapit = (*i)->ndimnamelist.find((*k)->getName());
1626 if(tempmapit != (*i)->ndimnamelist.end())
1627 tempcorrecteddimname= tempmapit->second;
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);
1633 (*j)->setCorrectedDimensions(correcteddims);
1634 correcteddims.clear();
1643 void File::update_grid_field_corrected_dims() throw(Exception) {
1646 vector<Dimension*> correcteddims;
1647 string tempcorrecteddimname;
1649 Dimension *correcteddim;
1652 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1653 i != this->grids.end(); ++i){
1655 for (vector<Field *>::const_iterator j =
1656 (*i)->getDataFields().begin();
1657 j != (*i)->getDataFields().end(); ++j) {
1660 if((*j)->iscoard ==
false) {
1663 for(vector<Dimension *>::const_iterator k=(*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k){
1665 map<string,string>::iterator tempmapit;
1668 tempmapit = (*i)->ndimnamelist.find((*k)->getName());
1669 if(tempmapit != (*i)->ndimnamelist.end())
1670 tempcorrecteddimname= tempmapit->second;
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);
1676 (*j)->setCorrectedDimensions(correcteddims);
1677 correcteddims.clear();
1684 void File::handle_grid_cf_attrs() throw(Exception) {
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) {
1698 if((*j)->fieldtype == 0) {
1699 string tempcoordinates=
"";
1700 string tempfieldname=
"";
1701 string tempcorrectedfieldname=
"";
1703 for(vector<Dimension *>::const_iterator k=(*j)->getDimensions().begin();
1704 k!=(*j)->getDimensions().end();++k){
1707 map<string,string>::iterator tempmapit;
1708 map<string,string>::iterator tempmapit2;
1711 tempmapit = ((*i)->dimcvarlist).find((*k)->getName());
1712 if(tempmapit != ((*i)->dimcvarlist).end())
1713 tempfieldname = tempmapit->second;
1715 throw4(
"cannot find the dimension field name",
1716 (*i)->getName(),(*j)->getName(),(*k)->getName());
1719 tempmapit2 = ((*i)->ncvarnamelist).find(tempfieldname);
1720 if(tempmapit2 != ((*i)->ncvarnamelist).end())
1721 tempcorrectedfieldname = tempmapit2->second;
1723 throw4(
"cannot find the corrected dimension field name",
1724 (*i)->getName(),(*j)->getName(),(*k)->getName());
1727 tempcoordinates= tempcorrectedfieldname;
1729 tempcoordinates = tempcoordinates +
" "+tempcorrectedfieldname;
1735 (*j)->setCoordinates(tempcoordinates);
1739 if((*j)->fieldtype == 1) {
1740 string tempunits =
"degrees_north";
1741 (*j)->setUnits(tempunits);
1743 if((*j)->fieldtype == 2) {
1744 string tempunits =
"degrees_east";
1745 (*j)->setUnits(tempunits);
1752 if((*j)->fieldtype == 4) {
1753 string tempunits =
"level";
1754 (*j)->setUnits(tempunits);
1758 if((*j)->fieldtype == 5) {
1759 string tempunits =
"days since 1900-01-01 00:00:00";
1760 (*j)->setUnits(tempunits);
1768 if (
true == (*i)->addfvalueattr) {
1769 if((((*j)->getFillValue()).empty()) && ((*j)->getType()==DFNT_FLOAT32 )) {
1770 float tempfillvalue = FLT_MAX;
1771 (*j)->addFillValue(tempfillvalue);
1772 (*j)->setAddedFillValue(
true);
1780 void File::handle_grid_SOM_projection() throw(Exception) {
1789 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1790 i != this->grids.end(); ++i){
1791 if (GCTP_SOM == (*i)->getProjection().getCode()) {
1797 for(vector<Dimension *>::const_iterator j=(*i)->getDimensions().begin();
1798 j!=(*i)->getDimensions().end();++j){
1801 if(NBLOCK == (*j)->getSize()) {
1805 if ((*j)->getName().compare(0,3,
"SOM") == 0) {
1806 som_dimname = (*j)->getName();
1812 if(
""== som_dimname)
1813 throw4(
"Wrong number of block: The number of block of MISR SOM Grid ",
1814 (*i)->getName(),
" is not ",NBLOCK);
1816 map<string,string>::iterator tempmapit;
1819 string cor_som_dimname;
1820 tempmapit = (*i)->ndimnamelist.find(som_dimname);
1821 if(tempmapit != (*i)->ndimnamelist.end())
1822 cor_som_dimname = tempmapit->second;
1824 throw2(
"cannot find the corrected dimension name for ", som_dimname);
1827 string cor_som_cvname;
1830 for (vector<Field *>::iterator j = (*i)->datafields.begin();
1831 j != (*i)->datafields.end(); ) {
1835 if (1 == (*j)->fieldtype || 2 == (*j)->fieldtype) {
1837 Dimension *newdim =
new Dimension(som_dimname,NBLOCK);
1838 Dimension *newcor_dim =
new Dimension(cor_som_dimname,NBLOCK);
1839 vector<Dimension *>::iterator it_d;
1841 it_d = (*j)->dims.begin();
1842 (*j)->dims.insert(it_d,newdim);
1844 it_d = (*j)->correcteddims.begin();
1845 (*j)->correcteddims.insert(it_d,newcor_dim);
1854 if ( 4 == (*j)->fieldtype) {
1855 cor_som_cvname = (*j)->newname;
1857 j = (*i)->datafields.erase(j);
1872 for (vector<Field *>::const_iterator j = (*i)->getDataFields().begin();
1873 j != (*i)->getDataFields().end(); ++j) {
1875 if ( 0 == (*j)->fieldtype) {
1877 string temp_coordinates = (*j)->coordinates;
1880 found = temp_coordinates.find(cor_som_cvname);
1884 if (temp_coordinates.size() >cor_som_cvname.size())
1885 temp_coordinates.erase(found,cor_som_cvname.size()+1);
1887 temp_coordinates.erase(found,cor_som_cvname.size());
1889 else if (found != string::npos)
1890 temp_coordinates.erase(found-1,cor_som_cvname.size()+1);
1892 throw4(
"cannot find the coordinate variable ",cor_som_cvname,
1893 "from ",temp_coordinates);
1895 (*j)->setCoordinates(temp_coordinates);
1905 void File::check_swath_dimmap(
int numswath)
throw(Exception) {
1907 if(HDF4RequestHandler::get_disable_swath_dim_map() ==
true)
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) {
1925 if(tempnumdm != 0 && odd_num_map ==
false)
1926 handle_swath_dimmap =
true;
1935 bool fakedimmap =
false;
1939 if((this->swaths[0]->getName()).find(
"atml2")!=string::npos){
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;
1960 for (vector<Field *>::const_iterator j =
1961 this->swaths[0]->getDataFields().begin();
1962 j != this->swaths[0]->getDataFields().end(); ++j) {
1977 if(((*j)->getName()).find(
"Latitude") != string::npos){
1979 if ((*j)->getType() == DFNT_UINT16 ||
1980 (*j)->getType() == DFNT_INT16)
1981 (*j)->type = DFNT_FLOAT32;
1983 (*j)->fieldtype = 1;
1986 if((*j)->getRank() != 2)
1987 throw2(
"The lat/lon rank must be 2 for Java clients to work",
1990 (((*j)->getDimensions())[0])->getName(),(*j)->getName());
1994 if(((*j)->getName()).find(
"Longitude")!= string::npos) {
1996 if((*j)->getType() == DFNT_UINT16 ||
1997 (*j)->getType() == DFNT_INT16)
1998 (*j)->type = DFNT_FLOAT32;
2000 (*j)->fieldtype = 2;
2001 if((*j)->getRank() != 2)
2002 throw2(
"The lat/lon rank must be 2 for Java clients to work",
2005 (((*j)->getDimensions())[1])->getName(), (*j)->getName());
2017 if(
true == fakedimmap)
2018 handle_swath_dimmap =
false;
2025 void File::check_swath_dimmap_bk_compat(
int numswath){
2027 if(
true == handle_swath_dimmap) {
2029 if(numswath == 1 && (((this->swaths)[0])->name==
"MODIS_SWATH_Type_L1B"))
2030 backward_handle_swath_dimmap =
true;
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;
2044 if (
true == all_2_dimmaps_no_geodim)
2045 backward_handle_swath_dimmap =
true;
2053 void File::create_swath_latlon_dim_cvar_map() throw(Exception){
2055 vector<Field*> ori_lats;
2056 vector<Field*> ori_lons;
2057 if(handle_swath_dimmap ==
true && backward_handle_swath_dimmap ==
false) {
2062 multi_dimmap =
true;
2064 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2065 i != this->swaths.end(); ++i){
2067 bool has_cf_lat =
false;
2068 bool has_cf_lon =
false;
2070 for (vector<Field *>::const_iterator j =
2071 (*i)->getGeoFields().begin();
2072 j != (*i)->getGeoFields().end(); ++j) {
2077 if((*j)->getName()==
"Latitude" && (*j)->getRank() == 2){
2079 ori_lats.push_back(*j);
2081 else if((*j)->getName()==
"Longitude" && (*j)->getRank() == 2){
2083 ori_lons.push_back(*j);
2085 if(has_cf_lat ==
true && has_cf_lon ==
true)
2088 if(has_cf_lat ==
false || has_cf_lon ==
false) {
2089 multi_dimmap =
false;
2098 if(
true == multi_dimmap) {
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]);
2124 bool lat_in_geofields =
false;
2125 bool lon_in_geofields =
false;
2127 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2128 i != this->swaths.end(); ++i){
2130 int tempgeocount = 0;
2131 for (vector<Field *>::const_iterator j =
2132 (*i)->getGeoFields().begin();
2133 j != (*i)->getGeoFields().end(); ++j) {
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",
2142 lat_in_geofields =
true;
2153 if(handle_swath_dimmap ==
true) {
2156 if(
true == backward_handle_swath_dimmap) {
2159 for(vector<SwathDataset::DimensionMap *>::const_iterator
2160 l=(*i)->getDimensionMaps().begin(); l!=(*i)->getDimensionMaps().end();++l){
2164 if(((*j)->getDimensions()[0])->getName() == (*l)->getGeoDimension()) {
2172 (*j)->fieldtype = 1;
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",
2183 lon_in_geofields =
true;
2184 if((*j)->getRank() == 1) {
2194 (((*j)->getDimensions())[1])->getName(),
"Longitude");
2195 if(handle_swath_dimmap ==
true) {
2196 if(
true == backward_handle_swath_dimmap) {
2199 for(vector<SwathDataset::DimensionMap *>::const_iterator
2200 l=(*i)->getDimensionMaps().begin(); l!=(*i)->getDimensionMaps().end();++l){
2205 if(((*j)->getDimensions()[1])->getName() ==
2206 (*l)->getGeoDimension()) {
2208 (*l)->getDataDimension(),
"Longitude");
2214 (*j)->fieldtype = 2;
2217 if(tempgeocount == 2)
2224 if (lat_in_geofields!=lon_in_geofields)
2225 throw1(
"Latitude and longitude must be both under Geolocation fields or Data fields");
2228 if (!lat_in_geofields && !lon_in_geofields) {
2230 bool lat_in_datafields =
false;
2231 bool lon_in_datafields =
false;
2233 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2234 i != this->swaths.end(); ++i){
2236 int tempgeocount = 0;
2237 for (vector<Field *>::const_iterator j =
2238 (*i)->getDataFields().begin();
2239 j != (*i)->getDataFields().end(); ++j) {
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",
2250 lat_in_datafields =
true;
2259 (((*j)->getDimensions())[0])->getName(),
"Latitude");
2261 if(handle_swath_dimmap ==
true) {
2262 if(
true == backward_handle_swath_dimmap) {
2264 for(vector<SwathDataset::DimensionMap *>::const_iterator
2265 l=(*i)->getDimensionMaps().begin(); l!=(*i)->getDimensionMaps().end();++l){
2269 if(((*j)->getDimensions()[0])->getName() == (*l)->getGeoDimension()) {
2276 (*j)->fieldtype = 1;
2280 if((*j)->getName()==
"Longitude"){
2282 if((*j)->getRank() > 2) {
2283 throw2(
"Currently the lat/lon rank must be 1 or 2 for Java clients to work",
2289 lon_in_datafields =
true;
2290 if((*j)->getRank() == 1) {
2300 (((*j)->getDimensions())[1])->getName(),
"Longitude");
2301 if(handle_swath_dimmap ==
true) {
2302 if(
true == backward_handle_swath_dimmap) {
2304 for(vector<SwathDataset::DimensionMap *>::const_iterator
2305 l=(*i)->getDimensionMaps().begin(); l!=(*i)->getDimensionMaps().end();++l){
2308 if(((*j)->getDimensions()[1])->getName() == (*l)->getGeoDimension()) {
2310 (*l)->getDataDimension(),
"Longitude");
2316 (*j)->fieldtype = 2;
2319 if(tempgeocount == 2)
2326 if (lat_in_datafields!=lon_in_datafields)
2327 throw1(
"Latitude and longitude must be both under Geolocation fields or Data fields");
2335 void File:: create_swath_nonll_dim_cvar_map() throw(Exception)
2338 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2339 i != this->swaths.end(); ++i){
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);
2366 for (vector<Field *>::const_iterator j =
2367 (*i)->getGeoFields().begin();
2368 j != (*i)->getGeoFields().end(); ++j) {
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;
2384 if(((((*j)->getDimensions())[0])->getName())==(*j)->getName()){
2385 (*j)->oriname = (*j)->getName();
2388 (*j)->specialcoard =
true;
2392 (*j)->fieldtype = 3;
2402 for (vector<Field *>::const_iterator j =
2403 (*i)->getDataFields().begin();
2404 j != (*i)->getDataFields().end(); ++j) {
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;
2417 if(((((*j)->getDimensions())[0])->getName())==(*j)->getName()){
2418 (*j)->oriname = (*j)->getName();
2420 (*j)->specialcoard =
true;
2424 (*j)->fieldtype = 3;
2434 bool missingfield_unlim_flag =
false;
2435 Field *missingfield_unlim = NULL;
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()){
2442 Field *missingfield =
new Field();
2455 if(
true == multi_dimmap && (this->swaths.size() != 1)) {
2456 missingfield->name = (*j)->getName()+
"_"+(*i)->name;
2457 dim =
new Dimension(missingfield->name,(*j)->getSize());
2460 missingfield->name = (*j)->getName();
2461 dim =
new Dimension((*j)->getName(),(*j)->getSize());
2463 missingfield->rank = 1;
2464 missingfield->type = DFNT_INT32;
2467 missingfield->dims.push_back(dim);
2474 if(0 == (*j)->getSize()) {
2475 missingfield_unlim_flag =
true;
2476 missingfield_unlim = missingfield;
2480 missingfield->fieldtype = 4;
2482 (*i)->geofields.push_back(missingfield);
2484 (missingfield->getDimensions())[0]->getName(), missingfield->name);
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) {
2500 for (vector<Dimension *>::const_iterator k =
2501 (*j)->getDimensions().begin(); k!= (*j)->getDimensions().end();++k){
2503 if((*k)->getName() == (missingfield_unlim->getDimensions())[0]->getName()) {
2504 if((*k)->getSize()!= 0) {
2505 Dimension *dim = missingfield_unlim->getDimensions()[0];
2507 dim->dimsize = (*k)->getSize();
2508 missingfield_unlim_flag =
false;
2514 if(
false == missingfield_unlim_flag)
2519 (*i)->nonmisscvdimlist.clear();
2527 void File::handle_swath_dim_cvar_maps() throw(Exception) {
2530 vector <string> tempfieldnamelist;
2531 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2532 i != this->swaths.end(); ++i) {
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;
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){
2555 string new_field_name = (*j)->name+
"_"+(*i)->name;
2565 int total_fcounter = 0;
2570 map<string,string>tempncvarnamelist;
2571 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2572 i != this->swaths.end(); ++i){
2575 for (vector<Field *>::const_iterator j =
2576 (*i)->getGeoFields().begin();
2577 j != (*i)->getGeoFields().end(); ++j)
2580 (*j)->newname = tempfieldnamelist[total_fcounter];
2584 if((*j)->fieldtype!=0) {
2589 for (vector<Field *>::const_iterator j =
2590 (*i)->getDataFields().begin();
2591 j != (*i)->getDataFields().end(); ++j)
2593 (*j)->newname = tempfieldnamelist[total_fcounter];
2597 if((*j)->fieldtype!=0) {
2605 vector <string>tempalldimnamelist;
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)
2616 int total_dcounter = 0;
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){
2628 vector<Dimension*> correcteddims;
2629 string tempcorrecteddimname;
2630 Dimension *correcteddim;
2632 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2633 i != this->swaths.end(); ++i){
2636 for (vector<Field *>::const_iterator j =
2637 (*i)->getGeoFields().begin();
2638 j != (*i)->getGeoFields().end(); ++j) {
2640 for(vector<Dimension *>::const_iterator k=(*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k){
2642 map<string,string>::iterator tempmapit;
2645 if(handle_swath_dimmap ==
false || multi_dimmap ==
true) {
2648 tempmapit = (*i)->ndimnamelist.find((*k)->getName());
2649 if(tempmapit != (*i)->ndimnamelist.end())
2650 tempcorrecteddimname= tempmapit->second;
2652 throw4(
"cannot find the corrected dimension name",
2653 (*i)->getName(),(*j)->getName(),(*k)->getName());
2655 correcteddim =
new Dimension(tempcorrecteddimname,(*k)->getSize());
2659 bool isdimmapname =
false;
2662 for(vector<SwathDataset::DimensionMap *>::const_iterator
2663 l=(*i)->getDimensionMaps().begin(); l!=(*i)->getDimensionMaps().end();++l){
2666 if((*k)->getName() == (*l)->getGeoDimension()) {
2668 isdimmapname =
true;
2670 string temprepdimname = (*l)->getDataDimension();
2673 tempmapit = (*i)->ndimnamelist.find(temprepdimname);
2674 if(tempmapit != (*i)->ndimnamelist.end())
2675 tempcorrecteddimname= tempmapit->second;
2677 throw4(
"cannot find the corrected dimension name", (*i)->getName(),
2678 (*j)->getName(),(*k)->getName());
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) {
2687 correcteddim =
new Dimension(tempcorrecteddimname,(*m)->getSize());
2693 throw4(
"cannot find the corrected dimension size", (*i)->getName(),
2694 (*j)->getName(),(*k)->getName());
2698 if(
false == isdimmapname) {
2700 tempmapit = (*i)->ndimnamelist.find((*k)->getName());
2701 if(tempmapit != (*i)->ndimnamelist.end())
2702 tempcorrecteddimname= tempmapit->second;
2704 throw4(
"cannot find the corrected dimension name",
2705 (*i)->getName(),(*j)->getName(),(*k)->getName());
2707 correcteddim =
new Dimension(tempcorrecteddimname,(*k)->getSize());
2712 correcteddims.push_back(correcteddim);
2714 (*j)->setCorrectedDimensions(correcteddims);
2715 correcteddims.clear();
2719 for (vector<Field *>::const_iterator j =
2720 (*i)->getDataFields().begin();
2721 j != (*i)->getDataFields().end(); ++j) {
2723 for(vector<Dimension *>::const_iterator k=
2724 (*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k){
2726 if((handle_swath_dimmap ==
false) || multi_dimmap ==
true) {
2729 map<string,string>::iterator tempmapit;
2731 tempmapit = (*i)->ndimnamelist.find((*k)->getName());
2732 if(tempmapit != (*i)->ndimnamelist.end())
2733 tempcorrecteddimname= tempmapit->second;
2735 throw4(
"cannot find the corrected dimension name", (*i)->getName(),
2736 (*j)->getName(),(*k)->getName());
2738 correcteddim =
new Dimension(tempcorrecteddimname,(*k)->getSize());
2741 map<string,string>::iterator tempmapit;
2742 bool isdimmapname =
false;
2744 for(vector<SwathDataset::DimensionMap *>::const_iterator l=
2745 (*i)->getDimensionMaps().begin(); l!=(*i)->getDimensionMaps().end();++l){
2748 if((*k)->getName() == (*l)->getGeoDimension()) {
2749 isdimmapname =
true;
2751 string temprepdimname = (*l)->getDataDimension();
2754 tempmapit = (*i)->ndimnamelist.find(temprepdimname);
2755 if(tempmapit != (*i)->ndimnamelist.end())
2756 tempcorrecteddimname= tempmapit->second;
2758 throw4(
"cannot find the corrected dimension name",
2759 (*i)->getName(),(*j)->getName(),(*k)->getName());
2763 bool ddimsflag =
false;
2764 for(vector<Dimension *>::const_iterator m=
2765 (*i)->getDimensions().begin();m!=(*i)->getDimensions().end();++m) {
2768 if((*m)->getName() == temprepdimname) {
2769 correcteddim =
new Dimension(tempcorrecteddimname,(*m)->getSize());
2775 throw4(
"cannot find the corrected dimension size",
2776 (*i)->getName(),(*j)->getName(),(*k)->getName());
2784 tempmapit = (*i)->ndimnamelist.find((*k)->getName());
2785 if(tempmapit != (*i)->ndimnamelist.end())
2786 tempcorrecteddimname= tempmapit->second;
2788 throw4(
"cannot find the corrected dimension name",
2789 (*i)->getName(),(*j)->getName(),(*k)->getName());
2791 correcteddim =
new Dimension(tempcorrecteddimname,(*k)->getSize());
2795 correcteddims.push_back(correcteddim);
2797 (*j)->setCorrectedDimensions(correcteddims);
2798 correcteddims.clear();
2805 void File::handle_swath_cf_attrs() throw(Exception) {
2815 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2816 i != this->swaths.end(); ++i){
2819 for (vector<Field *>::const_iterator j =
2820 (*i)->getGeoFields().begin();
2821 j != (*i)->getGeoFields().end(); ++j) {
2824 if((*j)->fieldtype == 0) {
2825 string tempcoordinates=
"";
2826 string tempfieldname=
"";
2827 string tempcorrectedfieldname=
"";
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;
2836 for(vector<Dimension *>::const_iterator
2837 k=(*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k){
2840 map<string,string>::iterator tempmapit;
2841 map<string,string>::iterator tempmapit2;
2844 tempmapit = ((*i)->dimcvarlist).find((*k)->getName());
2845 if(tempmapit != ((*i)->dimcvarlist).end())
2846 tempfieldname = tempmapit->second;
2848 throw4(
"cannot find the dimension field name",(*i)->getName(),
2849 (*j)->getName(),(*k)->getName());
2852 tempmapit2 = ((*i)->ncvarnamelist).find(tempfieldname);
2853 if(tempmapit2 != ((*i)->ncvarnamelist).end())
2854 tempcorrectedfieldname = tempmapit2->second;
2856 throw4(
"cannot find the corrected dimension field name",
2857 (*i)->getName(),(*j)->getName(),(*k)->getName());
2859 if(
false == has_ll_coord)
2860 has_ll_coord= check_ll_in_coords(tempcorrectedfieldname);
2863 tempcoordinates= tempcorrectedfieldname;
2865 tempcoordinates = tempcoordinates +
" "+tempcorrectedfieldname;
2868 if(
true == has_ll_coord)
2869 (*j)->setCoordinates(tempcoordinates);
2874 if((*j)->fieldtype == 1) {
2875 string tempunits =
"degrees_north";
2876 (*j)->setUnits(tempunits);
2880 if((*j)->fieldtype == 2) {
2881 string tempunits =
"degrees_east";
2882 (*j)->setUnits(tempunits);
2889 if((*j)->fieldtype == 4) {
2890 string tempunits =
"level";
2891 (*j)->setUnits(tempunits);
2896 if((*j)->fieldtype == 5) {
2897 string tempunits =
"days since 1900-01-01 00:00:00";
2898 (*j)->setUnits(tempunits);
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);
2913 for (vector<Field *>::const_iterator j =
2914 (*i)->getDataFields().begin();
2915 j != (*i)->getDataFields().end(); ++j) {
2918 if((*j)->fieldtype == 0) {
2919 string tempcoordinates=
"";
2920 string tempfieldname=
"";
2921 string tempcorrectedfieldname=
"";
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;
2930 for(vector<Dimension *>::const_iterator k
2931 =(*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k){
2934 map<string,string>::iterator tempmapit;
2935 map<string,string>::iterator tempmapit2;
2938 tempmapit = ((*i)->dimcvarlist).find((*k)->getName());
2939 if(tempmapit != ((*i)->dimcvarlist).end())
2940 tempfieldname = tempmapit->second;
2942 throw4(
"cannot find the dimension field name",(*i)->getName(),
2943 (*j)->getName(),(*k)->getName());
2946 tempmapit2 = ((*i)->ncvarnamelist).find(tempfieldname);
2947 if(tempmapit2 != ((*i)->ncvarnamelist).end())
2948 tempcorrectedfieldname = tempmapit2->second;
2950 throw4(
"cannot find the corrected dimension field name",
2951 (*i)->getName(),(*j)->getName(),(*k)->getName());
2953 if(
false == has_ll_coord)
2954 has_ll_coord= check_ll_in_coords(tempcorrectedfieldname);
2957 tempcoordinates= tempcorrectedfieldname;
2959 tempcoordinates = tempcoordinates +
" "+tempcorrectedfieldname;
2962 if(
true == has_ll_coord)
2963 (*j)->setCoordinates(tempcoordinates);
2966 if(((*j)->fieldtype == 3)||((*j)->fieldtype == 4)) {
2967 string tempunits =
"level";
2968 (*j)->setUnits(tempunits);
2973 if((*j)->fieldtype == 5) {
2974 string tempunits =
"days since 1900-01-01 00:00:00";
2975 (*j)->setUnits(tempunits);
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);
2993 bool File::find_dim_in_dims(
const std::vector<Dimension*>&dims,
const std::string &dim_name) {
2995 bool ret_value =
false;
2996 for (
int i = 0; i <dims.size(); i++) {
2997 if((dims[i])->name == dim_name) {
3006 void File::check_dm_geo_dims_in_vars() {
3008 if(handle_swath_dimmap ==
false)
3010 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
3011 i != this->swaths.end(); ++i){
3014 if((*i)->get_num_map() > 0) {
3016 for (vector<Field *>::const_iterator j =
3017 (*i)->getDataFields().begin();
3018 j != (*i)->getDataFields().end(); ++j) {
3022 if((*j)->rank >=2) {
3023 for(vector<Dimension *>::const_iterator k=
3024 (*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k){
3028 bool not_match_geo_dim =
true;
3029 for(vector<SwathDataset::DimensionMap *>::const_iterator l=
3030 (*i)->getDimensionMaps().begin(); l!=(*i)->getDimensionMaps().end();++l){
3032 if(((*k)->getName() == (*l)->getGeoDimension()) && not_match_geo_dim){
3034 not_match_geo_dim =
false;
3040 if(match_dims == 2) {
3041 (*i)->GeoDim_in_vars =
true;
3046 if((*i)->GeoDim_in_vars ==
false) {
3047 for (vector<Field *>::const_iterator j =
3048 (*i)->getGeoFields().begin();
3049 j != (*i)->getGeoFields().end(); ++j) {
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){
3059 bool not_match_geo_dim =
true;
3061 for(vector<SwathDataset::DimensionMap *>::const_iterator l=
3062 (*i)->getDimensionMaps().begin(); l!=(*i)->getDimensionMaps().end();++l){
3064 if(((*k)->getName() == (*l)->getGeoDimension()) && not_match_geo_dim){
3066 not_match_geo_dim =
false;
3072 if(match_dims == 2){
3073 (*i)->GeoDim_in_vars =
true;
3086 bool SwathDataset::obtain_dmap_offset_inc(
const string& ori_dimname,
const string & mapped_dimname,
int &offset,
int&inc) {
3087 bool ret_value =
false;
3088 for(vector<DimensionMap *>::const_iterator
3089 i=this->dimmaps.begin(); i!=this->dimmaps.end();++i){
3090 if((*i)->geodim==ori_dimname && (*i)->datadim == mapped_dimname){
3091 offset = (*i)->offset;
3092 inc = (*i)->increment;
3104 void File::create_geo_varnames_list(vector<string> & geo_varnames,
const string & swathname,
3105 const string & fieldname,
int extra_ll_pairs,
bool oneswath) {
3107 if(
true == oneswath)
3108 geo_varnames.push_back(fieldname);
3110 string nfieldname = fieldname+
"_"+swathname;
3111 geo_varnames.push_back(nfieldname);
3113 for (
int i = 0; i <extra_ll_pairs;i++) {
3117 if(
true == oneswath)
3118 nfieldname = fieldname+
"_"+si.str();
3120 nfieldname = fieldname+
"_"+swathname+
"_"+si.str();
3121 geo_varnames.push_back(nfieldname);
3124 cerr<<
"ll_pairs is "<<extra_ll_pairs <<endl;
3125 for(
int i =0;i<geo_varnames.size();i++)
3126 cerr<<
"geo_varnames["<<i<<
"]= " <<geo_varnames[i] <<endl;
3132 void File::create_geo_dim_var_maps(SwathDataset*sd, Field*fd,
const vector<string>& lat_names,
3133 const vector<string>& lon_names,vector<Dimension*>& geo_var_dim1,
3134 vector<Dimension*>& geo_var_dim2) {
3135 string field_lat_dim1_name =(fd->dims)[0]->name;
3136 string field_lat_dim2_name =(fd->dims)[1]->name;
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;
3144 geo_var_dim1.push_back((fd->dims)[0]);
3145 geo_var_dim2.push_back((fd->dims)[1]);
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;
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]);
3171 Dimension *lat_dim =
new Dimension(data_dim1_name,dim1_size);
3172 geo_var_dim1.push_back(lat_dim);
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]);
3187 Dimension *lon_dim =
new Dimension(data_dim2_name,dim2_size);
3188 geo_var_dim2.push_back(lon_dim);
3195 for(
int i = 0; i<lat_names.size();i++) {
3205 void File::create_geo_vars(SwathDataset* sd,Field *orig_lat,Field*orig_lon,
3206 const vector<string>& lat_names,
const vector<string>& lon_names,
3207 vector<Dimension*>&geo_var_dim1,vector<Dimension*>&geo_var_dim2)
throw(Exception){
3218 for (vector<Field *>::iterator i = sd->geofields.begin();
3219 i!=sd->geofields.end();++i) {
3220 if((*i)->name ==
"Latitude")
3222 else if((*i)->name ==
"Longitude")
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;
3232 if(sd->GeoDim_in_vars ==
false) {
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;
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());
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);
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;
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);
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;
3274 cerr<<
"orig_lon "<<orig_lon->name <<endl;
3275 cerr<<
"orig_lon dim1 inc "<<orig_lon->ll_dim1_inc<<endl;
3276 cerr<<
"orig_lon dim1 offset "<<orig_lon->ll_dim1_offset<<endl;
3277 cerr<<
"orig_lon dim0 inc "<<orig_lon->ll_dim0_inc<<endl;
3278 cerr<<
"orig_lon dim0 offset "<<orig_lon->ll_dim0_offset<<endl;
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;
3293 if((this->swaths).size()>1) {
3294 orig_lat->name = lat_names[0];
3295 orig_lon->name = lon_names[0];
3299 orig_lat->fieldtype = 1;
3300 orig_lon->fieldtype = 2;
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;
3310 newlat->type = orig_lat->type;
3311 Field * newlon =
new Field();
3312 newlon->name = lon_names[i];
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;
3323 newlon->type = orig_lon->type;
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());
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);
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);
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;
3351 cerr<<
"newlon "<<newlon->name <<endl;
3352 cerr<<
"newlon dim1 inc "<<newlon->ll_dim1_inc<<endl;
3353 cerr<<
"newlon dim1 offset "<<newlon->ll_dim1_offset<<endl;
3354 cerr<<
"newlon dim0 inc "<<newlon->ll_dim0_inc<<endl;
3355 cerr<<
"newlon dim0 offset "<<newlon->ll_dim0_offset<<endl;
3357 sd->geofields.push_back(newlat);
3358 sd->geofields.push_back(newlon);
3362 for (vector<Field *>::const_iterator j =
3363 sd->getGeoFields().begin();
3364 j != sd->getGeoFields().end(); ++j) {
3365 cerr<<
"Field name: "<<(*j)->name <<endl;
3366 cerr<<
"Dimension name and size: "<<endl;
3367 for(vector<Dimension *>::const_iterator k=
3368 (*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k)
3369 cerr<<(*k)->getName() <<
": "<<(*k)->getSize() <<endl;
3378 void File::update_swath_dims_for_dimmap(SwathDataset* sd,
const std::vector<Dimension*>&geo_var_dim1,
3379 const std::vector<Dimension*>&geo_var_dim2) {
3384 for (vector<Field *>::const_iterator j = sd->getGeoFields().begin();
3385 j != sd->getGeoFields().end(); ++j) {
3387 if((*j)->fieldtype == 1 || (*j)->fieldtype == 2)
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;
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;
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;
3425 void File::create_swath_latlon_dim_cvar_map_for_dimmap(SwathDataset* sd, Field* ori_lat, Field* ori_lon)
throw(Exception) {
3427 bool one_swath = ((this->swaths).size() == 1);
3429 int num_extra_lat_lon_pairs = 0;
3432 if(sd->GeoDim_in_vars ==
true)
3433 cerr<<
" swath name is "<<sd->name <<endl;
3437 if(sd->GeoDim_in_vars ==
false)
3438 num_extra_lat_lon_pairs--;
3440 num_extra_lat_lon_pairs += (sd->num_map)/2;
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;
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);
3455 if((this->swaths).size() >1)
3456 update_swath_dims_for_dimmap(sd,geo_var_dim1,geo_var_dim2);
3463 void File::Prepare(
const char *eosfile_path)
throw(Exception)
3471 int numgrid = this->grids.size();
3472 int numswath = this->swaths.size();
3475 throw2(
"the number of grid is less than 0", eosfile_path);
3481 string DIMXNAME = this->get_geodim_x_name();
3483 string DIMYNAME = this->get_geodim_y_name();
3485 string LATFIELDNAME = this->get_latfield_name();
3487 string LONFIELDNAME = this->get_lonfield_name();
3490 string GEOGRIDNAME =
"location";
3495 check_onelatlon_grids();
3498 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
3499 i != this->grids.end(); ++i) {
3500 handle_one_grid_zdim(*i);
3504 if (
true == this->onelatlon)
3505 handle_onelatlon_grids();
3507 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
3508 i != this->grids.end(); ++i) {
3512 (*i)->setDimxName(DIMXNAME);
3513 (*i)->setDimyName(DIMYNAME);
3516 handle_one_grid_latlon(*i);
3521 handle_grid_dim_cvar_maps();
3524 handle_grid_coards();
3527 update_grid_field_corrected_dims();
3530 handle_grid_cf_attrs();
3533 handle_grid_SOM_projection();
3539 for(vector<GridDataset *>::const_iterator i = this->grids.begin();
3540 i != this->grids.end(); ++i){
3541 (*i)->SetScaleType((*i)->name);
3550 check_swath_dimmap(numswath);
3553 check_dm_geo_dims_in_vars();
3556 check_swath_dimmap_bk_compat(numswath);
3559 create_swath_latlon_dim_cvar_map();
3562 create_swath_nonll_dim_cvar_map();
3565 handle_swath_dim_cvar_maps();
3569 handle_swath_cf_attrs();
3572 for(vector<SwathDataset *>::const_iterator i = this->swaths.begin();
3573 i != this->swaths.end(); ++i)
3574 (*i)->SetScaleType((*i)->name);
3584 void correct_unlimited_missing_zdim(GridDataset* gdset)
throw(Exception) {
3586 for (vector<Field *>::const_iterator j =
3587 gdset->getDataFields().begin();
3588 j != gdset->getDataFields().end(); ++j) {
3591 if ((*j)->getRank()==1 && (*j)->){
3601 bool File::check_special_1d_grid() throw(Exception) {
3603 int numgrid = this->grids.size();
3604 int numswath = this->swaths.size();
3606 if (numgrid != 1 || numswath != 0)
3610 string DIMXNAME = this->get_geodim_x_name();
3611 string DIMYNAME = this->get_geodim_y_name();
3613 if(DIMXNAME !=
"XDim" || DIMYNAME !=
"YDim")
3616 int var_dimx_size = 0;
3617 int var_dimy_size = 0;
3619 GridDataset *mygrid = (this->grids)[0];
3621 int field_xydim_flag = 0;
3622 for (vector<Field *>::const_iterator i = mygrid->getDataFields().begin();
3623 i!= mygrid->getDataFields().end(); ++i) {
3625 if((*i)->name ==
"XDim"){
3627 var_dimx_size = ((*i)->getDimensions())[0]->getSize();
3629 if((*i)->name ==
"YDim"){
3631 var_dimy_size = ((*i)->getDimensions())[0]->getSize();
3634 if(2==field_xydim_flag)
3638 if(field_xydim_flag !=2)
3642 int xdimsize = mygrid->getInfo().getX();
3643 int ydimsize = mygrid->getInfo().getY();
3645 if(var_dimx_size != xdimsize || var_dimy_size != ydimsize)
3653 bool File::check_ll_in_coords(
const string& vname)
throw(Exception) {
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) {
3662 if((*j)->fieldtype == 1 || (*j)->fieldtype == 2) {
3663 if((*j)->getNewName() == vname) {
3671 for (vector<Field *>::const_iterator j =
3672 (*i)->getDataFields().begin();
3673 j != (*i)->getDataFields().end(); ++j) {
3676 if((*j)->fieldtype == 1 || (*j)->fieldtype == 2) {
3677 if((*j)->getNewName() == vname) {
3698 void Dataset::SetScaleType(
const string EOS2ObjName)
throw(Exception) {
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");
3721 vector<string> modis_div_scale_type;
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");
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";
3737 string modis_divequ_scale_group =
"MODIS_Grid";
3738 string modis_div_scale_group =
"MOD_Grid";
3742 string modis_equ_scale_group =
"MODIS_Grid_1km_2D";
3744 if(EOS2ObjName==
"mod05" || EOS2ObjName==
"mod06" || EOS2ObjName==
"mod07"
3745 || EOS2ObjName==
"mod08" || EOS2ObjName==
"atml2")
3747 scaletype = MODIS_MUL_SCALE;
3761 if(EOS2ObjName.find(
"MOD")==0 || EOS2ObjName.find(
"mod")==0)
3763 size_t pos = EOS2ObjName.rfind(modis_eq_scale_type);
3764 if(pos != string::npos &&
3765 (pos== (EOS2ObjName.length()-modis_eq_scale_type.length())))
3767 scaletype = MODIS_EQ_SCALE;
3771 for(
unsigned int k=0; k<modis_multi_scale_type.size(); k++)
3773 pos = EOS2ObjName.rfind(modis_multi_scale_type[k]);
3774 if (pos !=string::npos &&
3775 (pos== (EOS2ObjName.length()-modis_multi_scale_type[k].length())))
3777 scaletype = MODIS_MUL_SCALE;
3782 for(
unsigned int k=0; k<modis_div_scale_type.size(); k++)
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;
3792 if (EOS2ObjName !=
"MODIS_Grid_1km_2D")
3799 pos = EOS2ObjName.find(modis_divequ_scale_group);
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;
3811 scaletype = MODIS_DIV_SCALE;
3814 size_t div_scale_pos = EOS2ObjName.find(modis_div_scale_group);
3817 if ( 0 == div_scale_pos)
3818 scaletype = MODIS_DIV_SCALE;
3824 if (EOS2ObjName ==
"VIP_CMG_GRID")
3825 scaletype = MODIS_DIV_SCALE;
3828 int Dataset::obtain_dimsize_with_dimname(
const string & dimname) {
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;
3843 for_each(this->dims.begin(), this->dims.end(), delete_elem());
3844 for_each(this->correcteddims.begin(), this->correcteddims.end(), delete_elem());
3850 for_each(this->dims.begin(), this->dims.end(), delete_elem());
3851 for_each(this->datafields.begin(), this->datafields.end(),
3853 for_each(this->attrs.begin(), this->attrs.end(), delete_elem());
3857 void Dataset::ReadDimensions(int32 (*entries)(int32, int32, int32 *),
3858 int32 (*inq)(int32,
char *, int32 *),
3859 vector<Dimension *> &d_dims)
throw(Exception)
3869 if ((numdims = entries(this->datasetid, HDFE_NENTDIM, &bufsize)) == -1)
3870 throw2(
"dimension entry", this->name);
3874 vector<char> namelist;
3875 vector<int32> dimsize;
3877 namelist.resize(bufsize + 1);
3878 dimsize.resize(numdims);
3881 if (inq(this->datasetid, &namelist[0], &dimsize[0]) == -1)
3882 throw2(
"inquire dimension", this->name);
3884 vector<string> dimnames;
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);
3900 void Dataset::ReadFields(int32 (*entries)(int32, int32, int32 *),
3901 int32 (*inq)(int32,
char *, int32 *, int32 *),
3903 (int32,
char *, int32 *, int32 *, int32 *,
char *),
3905 (int32,
char *, int32 *, int32 *, int32 *, VOIDP),
3906 intn (*getfill)(int32,
char *, VOIDP),
3908 vector<Field *> &fields)
throw(Exception)
3912 int32 numfields = 0;
3919 if ((numfields = entries(this->datasetid, geofield ?
3920 HDFE_NENTGFLD : HDFE_NENTDFLD, &bufsize)) == -1)
3921 throw2(
"field entry", this->name);
3925 if (numfields > 0) {
3926 vector<char> namelist;
3928 namelist.resize(bufsize + 1);
3931 if (inq(this->datasetid, &namelist[0], NULL, NULL) == -1)
3932 throw2(
"inquire field", this->name);
3934 vector<string> fieldnames;
3939 for (vector<string>::const_iterator i = fieldnames.begin();
3940 i != fieldnames.end(); ++i) {
3942 Field *field =
new Field();
3944 throw1(
"The field is NULL");
3947 bool throw_error =
false;
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;
3961 err_msg =
"Obtain field info error for field name "+fieldname_for_eh;
3964 if(
false == throw_error) {
3966 vector<string> dimnames;
3970 if ((
int)dimnames.size() != field->rank) {
3972 err_msg =
"Dimension names size is not consistent with field rank. ";
3973 err_msg +=
"Field name is "+field->name;
3976 for (
int k = 0; k < field->rank; ++k) {
3977 Dimension *dim =
new Dimension(dimnames[k], dimsize[k]);
3978 field->dims.push_back(dim);
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();
3989 fields.push_back(field);
3993 if(
true == throw_error) {
4003 void Dataset::ReadAttributes(int32 (*inq)(int32,
char *, int32 *),
4004 intn (*attrinfo)(int32,
char *, int32 *, int32 *),
4005 intn (*readattr)(int32,
char *, VOIDP),
4006 vector<Attribute *> &obj_attrs)
throw(Exception)
4015 if ((numattrs = inq(this->datasetid, NULL, &bufsize)) == -1)
4016 throw2(
"inquire attribute", this->name);
4020 vector<char> namelist;
4022 namelist.resize(bufsize + 1);
4025 if (inq(this->datasetid, &namelist[0], &bufsize) == -1)
4026 throw2(
"inquire attribute", this->name);
4028 vector<string> attrnames;
4033 for (vector<string>::const_iterator i = attrnames.begin();
4034 i != attrnames.end(); ++i) {
4035 Attribute *attr =
new Attribute();
4042 if (attrinfo(this->datasetid,
const_cast<char *
>
4043 (attr->name.c_str()), &attr->type, &count) == -1) {
4045 throw3(
"attribute info", this->name, attr->name);
4048 attr->count = count/DFKNTsize(attr->type);
4049 attr->value.resize(count);
4056 if (readattr(this->datasetid,
4057 const_cast<char *
>(attr->name.c_str()),
4058 &attr->value[0]) == -1) {
4060 throw3(
"read attribute", this->name, attr->name);
4064 obj_attrs.push_back(attr);
4070 GridDataset::~GridDataset()
4072 if (this->datasetid != -1){
4073 GDdetach(this->datasetid);
4078 GridDataset * GridDataset::Read(int32 fd,
const string &gridname)
4082 bool GD_fun_err =
false;
4083 GridDataset *grid =
new GridDataset(gridname);
4086 if ((grid->datasetid = GDattach(fd,
const_cast<char *
>(gridname.c_str())))
4088 err_msg =
"attach grid";
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";
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";
4114 if (GDpixreginfo(grid->datasetid, &proj.pix) == -1) {
4115 err_msg =
"pixreg info";
4119 if (GDorigininfo(grid->datasetid, &proj.origin) == -1){
4120 err_msg =
"origin info";
4126 if(
true == GD_fun_err){
4128 throw2(err_msg,gridname);
4133 grid->ReadDimensions(GDnentries, GDinqdims, grid->dims);
4136 grid->ReadFields(GDnentries, GDinqfields, GDfieldinfo, GDreadfield,
4137 GDgetfillvalue,
false, grid->datafields);
4140 grid->ReadAttributes(GDinqattrs, GDattrinfo, GDreadattr, grid->attrs);
4150 GridDataset::Calculated & GridDataset::getCalculated()
const
4152 if (this->calculated.grid !=
this)
4153 this->calculated.grid =
this;
4154 return this->calculated;
4157 bool GridDataset::Calculated::isYDimMajor() throw(Exception)
4159 this->DetectMajorDimension();
4160 return this->ydimmajor;
4164 bool GridDataset::Calculated::isOrthogonal() throw(Exception)
4167 this->ReadLongitudeLatitude();
4168 return this->orthogonal;
4172 int GridDataset::Calculated::DetectFieldMajorDimension() throw(Exception)
4177 for (vector<Field *>::const_iterator i =
4178 this->grid->getDataFields().begin();
4179 i != this->grid->getDataFields().end(); ++i) {
4181 int xdimindex = -1, ydimindex = -1, index = 0;
4184 for (vector<Dimension *>::const_iterator j =
4185 (*i)->getDimensions().begin();
4186 j != (*i)->getDimensions().end(); ++j) {
4187 if ((*j)->getName() == this->grid->dimxname)
4189 else if ((*j)->getName() == this->grid->dimyname)
4193 if (xdimindex == -1 || ydimindex == -1)
4196 int major = ydimindex < xdimindex ? 1 : 0;
4207 else if (ym != major)
4208 throw2(
"inconsistent XDim, YDim order", this->grid->getName());
4214 throw2(
"lack of data fields", this->grid->getName());
4219 void GridDataset::Calculated::DetectMajorDimension() throw(Exception)
4227 for (vector<Field *>::const_iterator i =
4228 this->grid->getDataFields().begin();
4229 i != this->grid->getDataFields().end(); ++i) {
4231 int xdimindex = -1, ydimindex = -1, index = 0;
4234 for (vector<Dimension *>::const_iterator j =
4235 (*i)->getDimensions().begin();
4236 j != (*i)->getDimensions().end(); ++j) {
4237 if ((*j)->getName() == this->grid->dimxname)
4239 else if ((*j)->getName() == this->grid->dimyname)
4243 if (xdimindex == -1 || ydimindex == -1)
4248 if(this->grid->getProjection().getCode() == GCTP_LAMAZ)
4251 major = ydimindex < xdimindex ? 1 : 0;
4262 else if (ym != major)
4263 throw2(
"inconsistent XDim, YDim order", this->grid->getName());
4268 throw2(
"lack of data fields", this->grid->getName());
4269 this->ydimmajor = ym != 0;
4278 static bool IsDisjoint(
const vector<Field *> &l,
4279 const vector<Field *> &r)
4281 for (vector<Field *>::const_iterator i = l.begin(); i != l.end(); ++i)
4283 if (find(r.begin(), r.end(), *i) != r.end())
4291 static bool IsDisjoint(vector<pair<Field *, string> > &l,
const vector<Field *> &r)
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())
4302 static bool IsSubset(
const vector<Field *> &s,
const vector<Field *> &b)
4304 for (vector<Field *>::const_iterator i = s.begin(); i != s.end(); ++i)
4306 if (find(b.begin(), b.end(), *i) == b.end())
4313 static bool IsSubset(vector<pair<Field *, string> > &s,
const vector<Field *> &b)
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())
4325 SwathDataset::~SwathDataset()
4327 if (this->datasetid != -1) {
4328 SWdetach(this->datasetid);
4331 for_each(this->dimmaps.begin(), this->dimmaps.end(), delete_elem());
4332 for_each(this->indexmaps.begin(), this->indexmaps.end(),
4335 for_each(this->geofields.begin(), this->geofields.end(),
4342 SwathDataset * SwathDataset::Read(int32 fd,
const string &swathname)
4345 SwathDataset *swath =
new SwathDataset(swathname);
4347 throw1(
"Cannot allocate HDF5 Swath object");
4350 if ((swath->datasetid = SWattach(fd,
4351 const_cast<char *
>(swathname.c_str())))
4354 throw2(
"attach swath", swathname);
4362 swath->ReadDimensions(SWnentries, SWinqdims, swath->dims);
4365 swath->ReadFields(SWnentries, SWinqdatafields, SWfieldinfo, SWreadfield,
4366 SWgetfillvalue,
false, swath->datafields);
4369 swath->ReadFields(SWnentries, SWinqgeofields, SWfieldinfo, SWreadfield,
4370 SWgetfillvalue,
true, swath->geofields);
4373 swath->ReadAttributes(SWinqattrs, SWattrinfo, SWreadattr, swath->attrs);
4376 swath->set_num_map(swath->ReadDimensionMaps(swath->dimmaps));
4379 swath->ReadIndexMaps(swath->indexmaps);
4391 int SwathDataset::ReadDimensionMaps(vector<DimensionMap *> &swath_dimmaps)
4394 int32 nummaps, bufsize;
4397 if ((nummaps = SWnentries(this->datasetid, HDFE_NENTMAP, &bufsize)) == -1)
4398 throw2(
"dimmap entry", this->name);
4409 vector<char> namelist;
4410 vector<int32> offset, increment;
4412 namelist.resize(bufsize + 1);
4413 offset.resize(nummaps);
4414 increment.resize(nummaps);
4415 if (SWinqmaps(this->datasetid, &namelist[0], &offset[0], &increment[0])
4417 throw2(
"inquire dimmap", this->name);
4419 vector<string> mapnames;
4422 for (vector<string>::const_iterator i = mapnames.begin();
4423 i != mapnames.end(); ++i) {
4424 vector<string> parts;
4426 if (parts.size() != 2)
4427 throw3(
"dimmap part", parts.size(),
4430 DimensionMap *dimmap =
new DimensionMap(parts[0], parts[1],
4433 swath_dimmaps.push_back(dimmap);
4441 void SwathDataset::ReadIndexMaps(vector<IndexMap *> &swath_indexmaps)
4444 int32 numindices, bufsize;
4446 if ((numindices = SWnentries(this->datasetid, HDFE_NENTIMAP, &bufsize)) ==
4448 throw2(
"indexmap entry", this->name);
4449 if (numindices > 0) {
4451 vector<char> namelist;
4453 namelist.resize(bufsize + 1);
4454 if (SWinqidxmaps(this->datasetid, &namelist[0], NULL) == -1)
4455 throw2(
"inquire indexmap", this->name);
4457 vector<string> mapnames;
4459 for (vector<string>::const_iterator i = mapnames.begin();
4460 i != mapnames.end(); ++i) {
4461 IndexMap *indexmap =
new IndexMap();
4462 vector<string> parts;
4464 if (parts.size() != 2)
4465 throw3(
"indexmap part", parts.size(),
4467 indexmap->geo = parts[0];
4468 indexmap->data = parts[1];
4473 SWdiminfo(this->datasetid,
4474 const_cast<char *
>(indexmap->geo.c_str())))
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,
4486 swath_indexmaps.push_back(indexmap);
4492 PointDataset::~PointDataset()
4496 PointDataset * PointDataset::Read(int32 ,
const string &pointname)
4499 PointDataset *point =
new PointDataset(pointname);
4505 bool Utility::ReadNamelist(
const char *path,
4506 int32 (*inq)(
char *,
char *, int32 *),
4507 vector<string> &names)
4509 char *fname =
const_cast<char *
>(path);
4513 if ((numobjs = inq(fname, NULL, &bufsize)) == -1)
return false;
4515 vector<char> buffer;
4516 buffer.resize(bufsize + 1);
4517 if (inq(fname, &buffer[0], &bufsize) == -1)
return false;
#define throw1(a1)
The followings are convenient functions to throw exceptions with different.
static bool insert_map(std::map< std::string, std::string > &m, std::string key, std::string val)
static void Handle_NameClashing(std::vector< std::string > &newobjnamelist)
General routines to handle name clashings.
static std::string get_CF_string(std::string s)
Change special characters to "_".
static void Split(const char *s, int len, char sep, std::vector< std::string > &names)