20 #include "HDFCFUtil.h"
22 #include "HDF4RequestHandler.h"
32 const char *HDFEOS2::File::_geodim_x_names[] = {
"XDim",
"LonDim",
"nlon"};
35 const char *HDFEOS2::File::_geodim_y_names[] = {
"YDim",
"LatDim",
"nlat"};
38 const char *HDFEOS2::File::_latfield_names[] = {
39 "Latitude",
"Lat",
"YDim",
"LatCenter"
43 const char *HDFEOS2::File::_lonfield_names[] = {
44 "Longitude",
"Lon",
"XDim",
"LonCenter"
49 const char *HDFEOS2::File::_geogrid_names[] = {
"location"};
51 using namespace HDFEOS2;
55 template<
typename T,
typename U,
typename V,
typename W,
typename X>
56 static void _throw5(
const char *fname,
int line,
int numarg,
57 const T &a1,
const U &a2,
const V &a3,
const W &a4,
61 ss << fname <<
":" << line <<
":";
62 for (
int i = 0; i < numarg; ++i) {
65 case 0: ss << a1;
break;
66 case 1: ss << a2;
break;
67 case 2: ss << a3;
break;
68 case 3: ss << a4;
break;
69 case 4: ss << a5;
break;
71 ss<<
" Argument number is beyond 5 ";
74 throw Exception(ss.str());
80 #define throw1(a1) _throw5(__FILE__, __LINE__, 1, a1, 0, 0, 0, 0)
81 #define throw2(a1, a2) _throw5(__FILE__, __LINE__, 2, a1, a2, 0, 0, 0)
82 #define throw3(a1, a2, a3) _throw5(__FILE__, __LINE__, 3, a1, a2, a3, 0, 0)
83 #define throw4(a1, a2, a3, a4) _throw5(__FILE__, __LINE__, 4, a1, a2, a3, a4, 0)
84 #define throw5(a1, a2, a3, a4, a5) _throw5(__FILE__, __LINE__, 5, a1, a2, a3, a4, a5)
86 #define assert_throw0(e) do { if (!(e)) throw1("assertion failure"); } while (false)
87 #define assert_range_throw0(e, ge, l) assert_throw0((ge) <= (e) && (e) < (l))
92 template<
typename T>
void operator()(T *ptr)
102 for (vector<GridDataset *>::const_iterator i = grids.begin();
103 i != grids.end(); ++i){
110 for (vector<SwathDataset *>::const_iterator i = swaths.begin();
111 i != swaths.end(); ++i){
117 for (vector<PointDataset *>::const_iterator i = points.begin();
118 i != points.end(); ++i){
125 File * File::Read(
const char *path, int32 mygridfd, int32 myswathfd)
throw(Exception)
128 File *file =
new File(path);
130 throw1(
"Memory allocation for file class failed. ");
132 file->gridfd = mygridfd;
133 file->swathfd = myswathfd;
137 if ((file->gridfd = GDopen(
const_cast<char *
>(file->path.c_str()),
138 DFACC_READ)) == -1) {
140 throw2(
"grid open", path);
144 vector<string> gridlist;
145 if (!Utility::ReadNamelist(file->path.c_str(), GDinqgrid, gridlist)) {
147 throw1(
"Grid ReadNamelist failed.");
151 for (vector<string>::const_iterator i = gridlist.begin();
152 i != gridlist.end(); ++i)
153 file->grids.push_back(GridDataset::Read(file->gridfd, *i));
157 throw1(
"GridDataset Read failed");
162 if ((file->swathfd = SWopen(
const_cast<char *
>(file->path.c_str()),
165 throw2(
"swath open", path);
169 vector<string> swathlist;
170 if (!Utility::ReadNamelist(file->path.c_str(), SWinqswath, swathlist)){
172 throw1(
"Swath ReadNamelist failed.");
176 for (vector<string>::const_iterator i = swathlist.begin();
177 i != swathlist.end(); ++i)
178 file->swaths.push_back(SwathDataset::Read(file->swathfd, *i));
182 throw1(
"SwathDataset Read failed.");
189 vector<string> pointlist;
190 if (!Utility::ReadNamelist(file->path.c_str(), PTinqpoint, pointlist)){
192 throw1(
"Point ReadNamelist failed.");
195 for (vector<string>::const_iterator i = pointlist.begin();
196 i != pointlist.end(); ++i)
197 file->points.push_back(PointDataset::Read(-1, *i));
200 if(file->grids.size() == 0 && file->swaths.size() == 0
201 && file->points.size() == 0) {
202 Exception e(
"Not an HDF-EOS2 file");
203 e.setFileType(
false);
214 string File::get_geodim_x_name()
216 if(!_geodim_x_name.empty())
217 return _geodim_x_name;
218 _find_geodim_names();
219 return _geodim_x_name;
225 string File::get_geodim_y_name()
227 if(!_geodim_y_name.empty())
228 return _geodim_y_name;
229 _find_geodim_names();
230 return _geodim_y_name;
240 string File::get_latfield_name()
242 if(!_latfield_name.empty())
243 return _latfield_name;
244 _find_latlonfield_names();
245 return _latfield_name;
248 string File::get_lonfield_name()
250 if(!_lonfield_name.empty())
251 return _lonfield_name;
252 _find_latlonfield_names();
253 return _lonfield_name;
260 string File::get_geogrid_name()
262 if(!_geogrid_name.empty())
263 return _geogrid_name;
264 _find_geogrid_name();
265 return _geogrid_name;
273 void File::_find_geodim_names()
275 set<string> geodim_x_name_set;
276 for(
size_t i = 0; i<
sizeof(_geodim_x_names) /
sizeof(
const char *); i++)
277 geodim_x_name_set.insert(_geodim_x_names[i]);
279 set<string> geodim_y_name_set;
280 for(
size_t i = 0; i<
sizeof(_geodim_y_names) /
sizeof(
const char *); i++)
281 geodim_y_name_set.insert(_geodim_y_names[i]);
286 const size_t gs = grids.size();
287 const size_t ss = swaths.size();
288 for(
size_t i=0; ;i++)
290 Dataset *dataset=NULL;
292 dataset =
static_cast<Dataset*
>(grids[i]);
294 dataset =
static_cast<Dataset*
>(swaths[i-gs]);
298 const vector<Dimension *>& dims = dataset->getDimensions();
299 for(vector<Dimension*>::const_iterator it = dims.begin();
300 it != dims.end(); ++it)
302 if(geodim_x_name_set.find((*it)->getName()) != geodim_x_name_set.end())
303 _geodim_x_name = (*it)->getName();
304 else if(geodim_y_name_set.find((*it)->getName()) != geodim_y_name_set.end())
305 _geodim_y_name = (*it)->getName();
312 const size_t gs = grids.size();
316 Dataset *dataset=NULL;
317 dataset =
static_cast<Dataset*
>(grids[0]);
319 const vector<Dimension *>& dims = dataset->getDimensions();
320 for(vector<Dimension*>::const_iterator it = dims.begin();
321 it != dims.end(); ++it)
329 if(geodim_x_name_set.find((*it)->getName()) != geodim_x_name_set.end())
330 _geodim_x_name = (*it)->getName();
331 else if(geodim_y_name_set.find((*it)->getName()) != geodim_y_name_set.end())
332 _geodim_y_name = (*it)->getName();
335 if(_geodim_x_name.empty())
336 _geodim_x_name = _geodim_x_names[0];
337 if(_geodim_y_name.empty())
338 _geodim_y_name = _geodim_y_names[0];
346 void File::_find_latlonfield_names()
348 set<string> latfield_name_set;
349 for(
size_t i = 0; i<
sizeof(_latfield_names) /
sizeof(
const char *); i++)
350 latfield_name_set.insert(_latfield_names[i]);
352 set<string> lonfield_name_set;
353 for(
size_t i = 0; i<
sizeof(_lonfield_names) /
sizeof(
const char *); i++)
354 lonfield_name_set.insert(_lonfield_names[i]);
356 const size_t gs = grids.size();
357 const size_t ss = swaths.size();
361 for(
size_t i=0;i<1 ;i++)
363 Dataset *dataset = NULL;
364 SwathDataset *sw = NULL;
366 dataset =
static_cast<Dataset*
>(grids[i]);
370 dataset =
static_cast<Dataset*
>(sw);
375 const vector<Field *>& fields = dataset->getDataFields();
376 for(vector<Field*>::const_iterator it = fields.begin();
377 it != fields.end(); ++it)
379 if(latfield_name_set.find((*it)->getName()) != latfield_name_set.end())
380 _latfield_name = (*it)->getName();
381 else if(lonfield_name_set.find((*it)->getName()) != lonfield_name_set.end())
382 _lonfield_name = (*it)->getName();
387 const vector<Field *>& geofields = dataset->getDataFields();
388 for(vector<Field*>::const_iterator it = geofields.begin();
389 it != geofields.end(); ++it)
391 if(latfield_name_set.find((*it)->getName()) != latfield_name_set.end())
392 _latfield_name = (*it)->getName();
393 else if(lonfield_name_set.find((*it)->getName()) != lonfield_name_set.end())
394 _lonfield_name = (*it)->getName();
401 if(_latfield_name.empty())
402 _latfield_name = _latfield_names[0];
403 if(_lonfield_name.empty())
404 _lonfield_name = _lonfield_names[0];
412 void File::_find_geogrid_name()
414 set<string> geogrid_name_set;
415 for(
size_t i = 0; i<
sizeof(_geogrid_names) /
sizeof(
const char *); i++)
416 geogrid_name_set.insert(_geogrid_names[i]);
418 const size_t gs = grids.size();
419 const size_t ss = swaths.size();
420 for(
size_t i=0; ;i++)
424 dataset =
static_cast<Dataset*
>(grids[i]);
426 dataset =
static_cast<Dataset*
>(swaths[i-gs]);
430 if(geogrid_name_set.find(dataset->getName()) != geogrid_name_set.end())
431 _geogrid_name = dataset->getName();
433 if(_geogrid_name.empty())
434 _geogrid_name =
"location";
438 void File::check_onelatlon_grids() {
441 string LATFIELDNAME = this->get_latfield_name();
442 string LONFIELDNAME = this->get_lonfield_name();
445 string GEOGRIDNAME =
"location";
454 for(vector<GridDataset *>::const_iterator i = this->grids.begin();
455 i != this->grids.end(); ++i){
458 for(vector<Field *>::const_iterator j =
459 (*i)->getDataFields().begin();
460 j != (*i)->getDataFields().end(); ++j) {
461 if((*i)->getName()==GEOGRIDNAME){
462 if((*j)->getName()==LATFIELDNAME){
466 if((*j)->getName()==LONFIELDNAME){
474 if(((*j)->getName()==LATFIELDNAME)||((*j)->getName()==LONFIELDNAME)){
475 (*i)->ownllflag =
true;
483 if(morellcount ==0 && onellcount ==2)
484 this->onelatlon =
true;
488 void File::handle_one_grid_zdim(GridDataset* gdset) {
491 string DIMXNAME = this->get_geodim_x_name();
492 string DIMYNAME = this->get_geodim_y_name();
494 bool missingfield_unlim_flag =
false;
495 Field *missingfield_unlim = NULL;
502 set<string> tempdimlist;
503 pair<set<string>::iterator,
bool> tempdimret;
505 for (vector<Field *>::const_iterator j =
506 gdset->getDataFields().begin();
507 j != gdset->getDataFields().end(); ++j) {
510 if ((*j)->getRank()==1){
514 if(((*j)->getDimensions())[0]->getName()!=DIMXNAME && ((*j)->getDimensions())[0]->getName()!=DIMYNAME){
516 tempdimret = tempdimlist.insert(((*j)->getDimensions())[0]->getName());
521 if(tempdimret.second ==
true) {
526 if((*j)->getName() ==
"Time")
542 for (vector<Dimension *>::const_iterator j =
543 gdset->getDimensions().begin(); j!= gdset->getDimensions().end();++j){
546 if((*j)->getName()!=DIMXNAME && (*j)->getName()!=DIMYNAME){
549 if((tempdimlist.find((*j)->getName())) == tempdimlist.end()){
552 Field *missingfield =
new Field();
553 missingfield->name = (*j)->getName();
554 missingfield->rank = 1;
557 missingfield->type = DFNT_INT32;
560 Dimension *dim =
new Dimension((*j)->getName(),(*j)->getSize());
563 missingfield->dims.push_back(dim);
569 int missingdimsize[1];
570 missingdimsize[0]= (*j)->getSize();
571 if(0 == (*j)->getSize()) {
572 missingfield_unlim_flag =
true;
573 missingfield_unlim = missingfield;
578 missingfield->fieldtype = 4;
588 gdset->datafields.push_back(missingfield);
597 bool temp_missingfield_unlim_flag = missingfield_unlim_flag;
598 if(
true == temp_missingfield_unlim_flag) {
599 for (vector<Field *>::const_iterator i =
600 gdset->getDataFields().begin();
601 i != gdset->getDataFields().end(); ++i) {
603 for (vector<Dimension *>::const_iterator j =
604 gdset->getDimensions().begin(); j!= gdset->getDimensions().end();++j){
606 if((*j)->getName() == (missingfield_unlim->getDimensions())[0]->getName()) {
607 if((*j)->getSize()!= 0) {
608 Dimension *dim = missingfield_unlim->getDimensions()[0];
611 dim->dimsize = (*j)->getSize();
612 missingfield_unlim_flag =
false;
618 if(
false == missingfield_unlim_flag)
626 void File::handle_one_grid_latlon(GridDataset* gdset)
throw(Exception)
630 string DIMXNAME = this->get_geodim_x_name();
631 string DIMYNAME = this->get_geodim_y_name();
633 string LATFIELDNAME = this->get_latfield_name();
634 string LONFIELDNAME = this->get_lonfield_name();
638 if(gdset->ownllflag) {
641 for (vector<Field *>::const_iterator j =
642 gdset->getDataFields().begin();
643 j != gdset->getDataFields().end(); ++j) {
645 if((*j)->getName() == LATFIELDNAME) {
656 if((*j)->getRank() > 2)
657 throw3(
"The rank of latitude is greater than 2",
658 gdset->getName(),(*j)->getName());
660 if((*j)->getRank() != 1) {
664 (*j)->ydimmajor = gdset->getCalculated().isYDimMajor();
670 int32 projectioncode = gdset->getProjection().getCode();
671 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
672 (*j)->condenseddim =
true;
682 if((*j)->condenseddim) {
686 for (vector<Dimension *>::const_iterator k =
687 (*j)->getDimensions().begin(); k!= (*j)->getDimensions().end();++k){
688 if((*k)->getName() == DIMYNAME) {
696 for (vector<Dimension *>::const_iterator k =
697 (*j)->getDimensions().begin(); k!= (*j)->getDimensions().end();++k){
698 if((*k)->getName() == DIMYNAME) {
710 else if ((*j)->getName() == LONFIELDNAME) {
719 if((*j)->getRank() >2)
720 throw3(
"The rank of Longitude is greater than 2",gdset->getName(),(*j)->getName());
722 if((*j)->getRank() != 1) {
725 (*j)->ydimmajor = gdset->getCalculated().isYDimMajor();
730 int32 projectioncode = gdset->getProjection().getCode();
731 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
732 (*j)->condenseddim =
true;
741 if((*j)->condenseddim) {
745 for (vector<Dimension *>::const_iterator k =
746 (*j)->getDimensions().begin(); k!= (*j)->getDimensions().end();++k){
747 if((*k)->getName() == DIMXNAME) {
754 for (vector<Dimension *>::const_iterator k =
755 (*j)->getDimensions().begin(); k!= (*j)->getDimensions().end();++k){
756 if((*k)->getName() == DIMXNAME) {
773 Field *latfield =
new Field();
774 Field *lonfield =
new Field();
776 latfield->name = LATFIELDNAME;
777 lonfield->name = LONFIELDNAME;
784 latfield->type = DFNT_FLOAT64;
785 lonfield->type = DFNT_FLOAT64;
788 latfield->fieldtype = 1;
791 lonfield->fieldtype = 2;
795 latfield->ydimmajor = gdset->getCalculated().isYDimMajor();
796 lonfield->ydimmajor = latfield->ydimmajor;
799 int xdimsize = gdset->getInfo().getX();
800 int ydimsize = gdset->getInfo().getY();
805 bool dmajor=(gdset->getProjection().getCode()==GCTP_LAMAZ)? gdset->getCalculated().DetectFieldMajorDimension()
806 : latfield->ydimmajor;
809 Dimension *dimlaty =
new Dimension(DIMYNAME,ydimsize);
810 latfield->dims.push_back(dimlaty);
811 Dimension *dimlony =
new Dimension(DIMYNAME,ydimsize);
812 lonfield->dims.push_back(dimlony);
813 Dimension *dimlatx =
new Dimension(DIMXNAME,xdimsize);
814 latfield->dims.push_back(dimlatx);
815 Dimension *dimlonx =
new Dimension(DIMXNAME,xdimsize);
816 lonfield->dims.push_back(dimlonx);
819 Dimension *dimlatx =
new Dimension(DIMXNAME,xdimsize);
820 latfield->dims.push_back(dimlatx);
821 Dimension *dimlonx =
new Dimension(DIMXNAME,xdimsize);
822 lonfield->dims.push_back(dimlonx);
823 Dimension *dimlaty =
new Dimension(DIMYNAME,ydimsize);
824 latfield->dims.push_back(dimlaty);
825 Dimension *dimlony =
new Dimension(DIMYNAME,ydimsize);
826 lonfield->dims.push_back(dimlony);
832 upleft =
const_cast<float64 *
>(gdset->getInfo().getUpLeft());
833 lowright =
const_cast<float64 *
>(gdset->getInfo().getLowRight());
836 int32 projectioncode = gdset->getProjection().getCode();
837 if(((
int)lowright[0]>180000000) && ((
int)upleft[0]>-1)) {
840 if(projectioncode == GCTP_GEO) {
841 lonfield->speciallon =
true;
847 if(HDF4RequestHandler::get_enable_eosgeo_cachefile() ==
true)
848 latfield->speciallon =
true;
860 if(((
int)(lowright[0]/1000)==0) &&((int)(upleft[0]/1000)==0)
861 && ((
int)(upleft[1]/1000)==0) && ((int)(lowright[1]/1000)==0)) {
862 if(projectioncode == GCTP_GEO){
863 lonfield->specialformat = 1;
864 latfield->specialformat = 1;
873 if(((
int)(lowright[0])==0) &&((int)(upleft[0])==0)
874 && ((
int)(upleft[1])==0) && ((int)(lowright[1])==0)) {
875 if(projectioncode == GCTP_GEO){
876 lonfield->specialformat = 2;
877 latfield->specialformat = 2;
878 gdset->addfvalueattr =
true;
889 if(((
int)(lowright[0])==-1) &&((int)(upleft[0])==-1)
890 && ((
int)(upleft[1])==-1) && ((int)(lowright[1])==-1)) {
891 lonfield->specialformat = 3;
892 latfield->specialformat = 3;
893 lonfield->condenseddim =
true;
894 latfield->condenseddim =
true;
899 if(GCTP_SOM == projectioncode) {
900 lonfield->specialformat = 4;
901 latfield->specialformat = 4;
907 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
908 lonfield->condenseddim =
true;
909 latfield->condenseddim =
true;
916 string check_eos_geo_cache_key =
"H4.EnableEOSGeoCacheFile";
917 bool enable_eos_geo_cache_key =
false;
918 enable_eos_geo_cache_key = HDFCFUtil::check_beskeys(check_eos_geo_cache_key);
919 if(
true == enable_eos_geo_cache_key) {
924 cache_fname =HDFCFUtil::get_int_str(gdset->getProjection().getCode());
925 cache_fname +=HDFCFUtil::get_int_str(gdset->getProjection().getZone());
926 cache_fname +=HDFCFUtil::get_int_str(gdset->getProjection().getSphere());
927 cache_fname +=HDFCFUtil::get_int_str(gdset->getProjection().getPix());
928 cache_fname +=HDFCFUtil::get_int_str(gdset->getProjection().getOrigin());
933 cache_fname +=HDFCFUtil::get_int_str(ydimsize);
934 cache_fname +=HDFCFUtil::get_int_str(xdimsize);
937 cache_fname +=HDFCFUtil::get_int_str(xdimsize);
938 cache_fname +=HDFCFUtil::get_int_str(ydimsize);
943 cache_fname +=HDFCFUtil::get_double_str(upleft[0],17,6);
944 cache_fname +=HDFCFUtil::get_double_str(upleft[1],17,6);
945 cache_fname +=HDFCFUtil::get_double_str(lowright[0],17,6);
946 cache_fname +=HDFCFUtil::get_double_str(lowright[1],17,6);
951 params =
const_cast<float64 *
>(gdset->getProjection().getParam());
954 for(
int ipar = 0; ipar<13;ipar++)
955 cache_fname+=HDFCFUtil::get_double_str(params[ipar],17,6);
964 string cache_fpath =
"/tmp/"+cache_fname;
965 int result = stat(cache_fpath.c_str(), &st);
967 int actual_file_size = st.st_size;
968 cerr<<
"HDF-EOS2 actual_file_size is "<<actual_file_size <<endl;
969 int expected_file_size = 0;
970 if(gdset->getProjection().getCode() == GCTP_SOM)
971 expected_file_size = xdimsize*ydimsize*2*
sizeof(double)*NBLOCK;
972 else if(gdset->getProjection().getCode() == GCTP_CEA ||
973 gdset->getProjection().getCode() == GCTP_GEO)
974 expected_file_size = (xdimsize+ydimsize)*
sizeof(double);
976 expected_file_size = xdimsize*ydimsize*2*
sizeof(double);
978 cerr<<
"HDF-EOS2 expected_file_size is "<<expected_file_size <<endl;
979 if(actual_file_size != expected_file_size){
980 cerr<<
"field_cache is 1 "<<endl;
981 lonfield->field_cache = 1;
982 latfield->field_cache = 1;
985 cerr<<
"field cache is 2 "<<endl;
986 lonfield->field_cache = 2;
987 latfield->field_cache = 2;
1003 gdset->datafields.push_back(latfield);
1004 gdset->datafields.push_back(lonfield);
1012 if(lonfield->condenseddim) {
1016 for (vector<Dimension *>::const_iterator j =
1017 lonfield->getDimensions().begin(); j!= lonfield->getDimensions().end();++j){
1019 if((*j)->getName() == DIMXNAME) {
1023 if((*j)->getName() == DIMYNAME) {
1029 for (vector<Dimension *>::const_iterator j =
1030 lonfield->getDimensions().begin(); j!= lonfield->getDimensions().end();++j){
1032 if((*j)->getName() == DIMXNAME){
1036 if((*j)->getName() == DIMYNAME){
1047 void File::handle_onelatlon_grids() throw(Exception) {
1050 string DIMXNAME = this->get_geodim_x_name();
1051 string DIMYNAME = this->get_geodim_y_name();
1052 string LATFIELDNAME = this->get_latfield_name();
1053 string LONFIELDNAME = this->get_lonfield_name();
1057 string GEOGRIDNAME =
"location";
1060 map<string,string>temponelatlondimcvarlist;
1063 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1064 i != this->grids.end(); ++i){
1068 (*i)->setDimxName(DIMXNAME);
1069 (*i)->setDimyName(DIMYNAME);
1072 if((*i)->getName()==GEOGRIDNAME) {
1077 (*i)->lonfield->fieldtype = 2;
1078 (*i)->latfield->fieldtype = 1;
1081 if((*i)->lonfield->rank >2 || (*i)->latfield->rank >2)
1082 throw2(
"Either the rank of lat or the lon is greater than 2",(*i)->getName());
1083 if((*i)->lonfield->rank !=(*i)->latfield->rank)
1084 throw2(
"The rank of the latitude is not the same as the rank of the longitude",(*i)->getName());
1087 if((*i)->lonfield->rank != 1) {
1091 (*i)->lonfield->ydimmajor = (*i)->getCalculated().isYDimMajor();
1092 (*i)->latfield->ydimmajor = (*i)->lonfield->ydimmajor;
1097 int32 projectioncode = (*i)->getProjection().getCode();
1098 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
1099 (*i)->lonfield->condenseddim =
true;
1100 (*i)->latfield->condenseddim =
true;
1107 if((*i)->lonfield->condenseddim) {
1113 for (vector<Dimension *>::const_iterator j =
1114 (*i)->lonfield->getDimensions().begin(); j!= (*i)->lonfield->getDimensions().end();++j){
1115 if((*j)->getName() == DIMXNAME) {
1117 (*i)->lonfield->getName());
1119 if((*j)->getName() == DIMYNAME) {
1121 (*i)->latfield->getName());
1128 for (vector<Dimension *>::const_iterator j =
1129 (*i)->lonfield->getDimensions().begin(); j!= (*i)->lonfield->getDimensions().end();++j){
1130 if((*j)->getName() == DIMXNAME) {
1132 (*i)->lonfield->getName());
1134 if((*j)->getName() == DIMYNAME) {
1136 (*i)->latfield->getName());
1144 (*i)->lonfield->getName());
1146 (*i)->latfield->getName());
1148 temponelatlondimcvarlist = (*i)->dimcvarlist;
1156 for(vector<GridDataset *>::const_iterator i = this->grids.begin();
1157 i != this->grids.end(); ++i){
1159 string templatlonname1;
1160 string templatlonname2;
1162 if((*i)->getName() != GEOGRIDNAME) {
1164 map<string,string>::iterator tempmapit;
1167 tempmapit = temponelatlondimcvarlist.find(DIMXNAME);
1168 if(tempmapit != temponelatlondimcvarlist.end())
1169 templatlonname1= tempmapit->second;
1171 throw2(
"cannot find the dimension field of XDim", (*i)->getName());
1176 tempmapit = temponelatlondimcvarlist.find(DIMYNAME);
1177 if(tempmapit != temponelatlondimcvarlist.end())
1178 templatlonname2= tempmapit->second;
1180 throw2(
"cannot find the dimension field of YDim", (*i)->getName());
1188 void File::handle_grid_dim_cvar_maps() throw(Exception) {
1191 string DIMXNAME = this->get_geodim_x_name();
1193 string DIMYNAME = this->get_geodim_y_name();
1195 string LATFIELDNAME = this->get_latfield_name();
1197 string LONFIELDNAME = this->get_lonfield_name();
1202 string GEOGRIDNAME =
"location";
1216 vector <string> tempfieldnamelist;
1217 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1218 i != this->grids.end(); ++i) {
1219 for (vector<Field *>::const_iterator j = (*i)->getDataFields().begin();
1220 j!= (*i)->getDataFields().end(); ++j) {
1230 map<string,string>tempncvarnamelist;
1231 string tempcorrectedlatname, tempcorrectedlonname;
1233 int total_fcounter = 0;
1235 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1236 i != this->grids.end(); ++i){
1241 for (vector<Field *>::const_iterator j =
1242 (*i)->getDataFields().begin();
1243 j != (*i)->getDataFields().end(); ++j)
1245 (*j)->newname = tempfieldnamelist[total_fcounter];
1249 if((*j)->fieldtype!=0) {
1251 tempncvarnamelist.insert(make_pair((*j)->getName(), (*j)->newname));
1254 if((this->onelatlon)&&(((*i)->getName())==GEOGRIDNAME)) {
1255 if((*j)->getName()==LATFIELDNAME)
1256 tempcorrectedlatname = (*j)->newname;
1257 if((*j)->getName()==LONFIELDNAME)
1258 tempcorrectedlonname = (*j)->newname;
1263 (*i)->ncvarnamelist = tempncvarnamelist;
1264 tempncvarnamelist.clear();
1270 if(this->onelatlon) {
1271 for(vector<GridDataset *>::const_iterator i = this->grids.begin();
1272 i != this->grids.end(); ++i){
1274 if((*i)->getName()!=GEOGRIDNAME){
1282 map<string,string>tempndimnamelist;
1285 vector <string>tempalldimnamelist;
1286 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1287 i != this->grids.end(); ++i)
1288 for (map<string,string>::const_iterator j =
1289 (*i)->dimcvarlist.begin(); j!= (*i)->dimcvarlist.end();++j)
1296 int total_dcounter = 0;
1297 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1298 i != this->grids.end(); ++i){
1300 for (map<string,string>::const_iterator j =
1301 (*i)->dimcvarlist.begin(); j!= (*i)->dimcvarlist.end();++j){
1304 if((DIMXNAME == (*j).first || DIMYNAME == (*j).first) && (
true==(this->onelatlon)))
1311 (*i)->ndimnamelist = tempndimnamelist;
1312 tempndimnamelist.clear();
1317 void File::handle_grid_coards() throw(Exception) {
1320 string DIMXNAME = this->get_geodim_x_name();
1321 string DIMYNAME = this->get_geodim_y_name();
1322 string LATFIELDNAME = this->get_latfield_name();
1323 string LONFIELDNAME = this->get_lonfield_name();
1327 string GEOGRIDNAME =
"location";
1331 vector<Dimension*> correcteddims;
1332 string tempcorrecteddimname;
1335 map<string,string> tempnewxdimnamelist;
1338 map<string,string> tempnewydimnamelist;
1341 Dimension *correcteddim;
1343 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1344 i != this->grids.end(); ++i){
1345 for (vector<Field *>::const_iterator j =
1346 (*i)->getDataFields().begin();
1347 j != (*i)->getDataFields().end(); ++j) {
1352 if((*j)->getName()==LATFIELDNAME && (*j)->getRank()==2 &&(*j)->condenseddim) {
1354 string templatdimname;
1355 map<string,string>::iterator tempmapit;
1358 tempmapit = (*i)->ncvarnamelist.find(LATFIELDNAME);
1359 if(tempmapit != (*i)->ncvarnamelist.end())
1360 templatdimname= tempmapit->second;
1362 throw2(
"cannot find the corrected field of Latitude", (*i)->getName());
1364 for(vector<Dimension *>::const_iterator k =(*j)->getDimensions().begin();
1365 k!=(*j)->getDimensions().end();++k){
1369 if((*k)->getName()==DIMYNAME) {
1370 correcteddim =
new Dimension(templatdimname,(*k)->getSize());
1371 correcteddims.push_back(correcteddim);
1372 (*j)->setCorrectedDimensions(correcteddims);
1376 (*j)->iscoard =
true;
1377 (*i)->iscoard =
true;
1379 this->iscoard =
true;
1382 correcteddims.clear();
1386 else if((*j)->getName()==LONFIELDNAME && (*j)->getRank()==2 &&(*j)->condenseddim){
1388 string templondimname;
1389 map<string,string>::iterator tempmapit;
1392 tempmapit = (*i)->ncvarnamelist.find(LONFIELDNAME);
1393 if(tempmapit != (*i)->ncvarnamelist.end())
1394 templondimname= tempmapit->second;
1396 throw2(
"cannot find the corrected field of Longitude", (*i)->getName());
1398 for(vector<Dimension *>::const_iterator k =(*j)->getDimensions().begin();
1399 k!=(*j)->getDimensions().end();++k){
1403 if((*k)->getName()==DIMXNAME) {
1404 correcteddim =
new Dimension(templondimname,(*k)->getSize());
1405 correcteddims.push_back(correcteddim);
1406 (*j)->setCorrectedDimensions(correcteddims);
1411 (*j)->iscoard =
true;
1412 (*i)->iscoard =
true;
1414 this->iscoard =
true;
1415 correcteddims.clear();
1419 else if(((*j)->getRank()==1) &&((*j)->getName()==LONFIELDNAME) ) {
1421 string templondimname;
1422 map<string,string>::iterator tempmapit;
1425 tempmapit = (*i)->ncvarnamelist.find(LONFIELDNAME);
1426 if(tempmapit != (*i)->ncvarnamelist.end())
1427 templondimname= tempmapit->second;
1429 throw2(
"cannot find the corrected field of Longitude", (*i)->getName());
1431 correcteddim =
new Dimension(templondimname,((*j)->getDimensions())[0]->getSize());
1432 correcteddims.push_back(correcteddim);
1433 (*j)->setCorrectedDimensions(correcteddims);
1434 (*j)->iscoard =
true;
1435 (*i)->iscoard =
true;
1437 this->iscoard =
true;
1438 correcteddims.clear();
1440 if((((*j)->getDimensions())[0]->getName()!=DIMXNAME)
1441 &&((((*j)->getDimensions())[0]->getName())!=DIMYNAME)){
1442 throw3(
"the dimension name of longitude should not be ",
1443 ((*j)->getDimensions())[0]->getName(),(*i)->getName());
1445 if((((*j)->getDimensions())[0]->getName())==DIMXNAME) {
1454 else if(((*j)->getRank()==1) &&((*j)->getName()==LATFIELDNAME) ) {
1456 string templatdimname;
1457 map<string,string>::iterator tempmapit;
1460 tempmapit = (*i)->ncvarnamelist.find(LATFIELDNAME);
1461 if(tempmapit != (*i)->ncvarnamelist.end())
1462 templatdimname= tempmapit->second;
1464 throw2(
"cannot find the corrected field of Latitude", (*i)->getName());
1466 correcteddim =
new Dimension(templatdimname,((*j)->getDimensions())[0]->getSize());
1467 correcteddims.push_back(correcteddim);
1468 (*j)->setCorrectedDimensions(correcteddims);
1470 (*j)->iscoard =
true;
1471 (*i)->iscoard =
true;
1473 this->iscoard =
true;
1474 correcteddims.clear();
1476 if(((((*j)->getDimensions())[0]->getName())!=DIMXNAME)
1477 &&((((*j)->getDimensions())[0]->getName())!=DIMYNAME))
1478 throw3(
"the dimension name of latitude should not be ",
1479 ((*j)->getDimensions())[0]->getName(),(*i)->getName());
1480 if((((*j)->getDimensions())[0]->getName())==DIMXNAME){
1494 if(
true == this->onelatlon){
1497 if(
true == this->iscoard){
1500 string tempcorrectedxdimname;
1501 string tempcorrectedydimname;
1503 if((
int)(tempnewxdimnamelist.size())!= 1)
1504 throw1(
"the corrected dimension name should have only one pair");
1505 if((
int)(tempnewydimnamelist.size())!= 1)
1506 throw1(
"the corrected dimension name should have only one pair");
1508 map<string,string>::iterator tempdimmapit = tempnewxdimnamelist.begin();
1509 tempcorrectedxdimname = tempdimmapit->second;
1510 tempdimmapit = tempnewydimnamelist.begin();
1511 tempcorrectedydimname = tempdimmapit->second;
1513 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1514 i != this->grids.end(); ++i){
1517 map<string,string>::iterator tempmapit;
1518 tempmapit = (*i)->ndimnamelist.find(DIMXNAME);
1519 if(tempmapit != (*i)->ndimnamelist.end()) {
1523 throw2(
"cannot find the corrected dimension name", (*i)->getName());
1524 tempmapit = (*i)->ndimnamelist.find(DIMYNAME);
1525 if(tempmapit != (*i)->ndimnamelist.end()) {
1529 throw2(
"cannot find the corrected dimension name", (*i)->getName());
1534 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1535 i != this->grids.end(); ++i){
1538 string tempcorrectedxdimname;
1539 string tempcorrectedydimname;
1542 map<string,string>::iterator tempdimmapit;
1543 map<string,string>::iterator tempmapit;
1544 tempdimmapit = tempnewxdimnamelist.find((*i)->getName());
1545 if(tempdimmapit != tempnewxdimnamelist.end())
1546 tempcorrectedxdimname = tempdimmapit->second;
1548 throw2(
"cannot find the corrected COARD XDim dimension name", (*i)->getName());
1549 tempmapit = (*i)->ndimnamelist.find(DIMXNAME);
1550 if(tempmapit != (*i)->ndimnamelist.end()) {
1554 throw2(
"cannot find the corrected dimension name", (*i)->getName());
1556 tempdimmapit = tempnewydimnamelist.find((*i)->getName());
1557 if(tempdimmapit != tempnewydimnamelist.end())
1558 tempcorrectedydimname = tempdimmapit->second;
1560 throw2(
"cannot find the corrected COARD YDim dimension name", (*i)->getName());
1562 tempmapit = (*i)->ndimnamelist.find(DIMYNAME);
1563 if(tempmapit != (*i)->ndimnamelist.end()) {
1567 throw2(
"cannot find the corrected dimension name", (*i)->getName());
1575 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1576 i != this->grids.end(); ++i){
1577 for (map<string,string>::const_iterator j =
1578 (*i)->dimcvarlist.begin(); j!= (*i)->dimcvarlist.end();++j){
1582 if((this->iscoard||(*i)->iscoard) && (*j).first !=DIMXNAME && (*j).first !=DIMYNAME) {
1583 string tempnewdimname;
1584 map<string,string>::iterator tempmapit;
1587 tempmapit = (*i)->ncvarnamelist.find((*j).second);
1588 if(tempmapit != (*i)->ncvarnamelist.end())
1589 tempnewdimname= tempmapit->second;
1591 throw3(
"cannot find the corrected field of ", (*j).second,(*i)->getName());
1594 tempmapit =(*i)->ndimnamelist.find((*j).first);
1595 if(tempmapit != (*i)->ndimnamelist.end())
1598 throw3(
"cannot find the corrected dimension name of ", (*j).first,(*i)->getName());
1605 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1606 i != this->grids.end(); ++i){
1608 for (vector<Field *>::const_iterator j =
1609 (*i)->getDataFields().begin();
1610 j != (*i)->getDataFields().end(); ++j) {
1613 if((*j)->iscoard ==
false) {
1616 for(vector<Dimension *>::const_iterator k=(*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k){
1619 map<string,string>::iterator tempmapit;
1622 tempmapit = (*i)->ndimnamelist.find((*k)->getName());
1623 if(tempmapit != (*i)->ndimnamelist.end())
1624 tempcorrecteddimname= tempmapit->second;
1626 throw4(
"cannot find the corrected dimension name", (*i)->getName(),(*j)->getName(),(*k)->getName());
1627 correcteddim =
new Dimension(tempcorrecteddimname,(*k)->getSize());
1628 correcteddims.push_back(correcteddim);
1630 (*j)->setCorrectedDimensions(correcteddims);
1631 correcteddims.clear();
1640 void File::update_grid_field_corrected_dims() throw(Exception) {
1643 vector<Dimension*> correcteddims;
1644 string tempcorrecteddimname;
1646 Dimension *correcteddim;
1649 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1650 i != this->grids.end(); ++i){
1652 for (vector<Field *>::const_iterator j =
1653 (*i)->getDataFields().begin();
1654 j != (*i)->getDataFields().end(); ++j) {
1657 if((*j)->iscoard ==
false) {
1660 for(vector<Dimension *>::const_iterator k=(*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k){
1662 map<string,string>::iterator tempmapit;
1665 tempmapit = (*i)->ndimnamelist.find((*k)->getName());
1666 if(tempmapit != (*i)->ndimnamelist.end())
1667 tempcorrecteddimname= tempmapit->second;
1669 throw4(
"cannot find the corrected dimension name", (*i)->getName(),(*j)->getName(),(*k)->getName());
1670 correcteddim =
new Dimension(tempcorrecteddimname,(*k)->getSize());
1671 correcteddims.push_back(correcteddim);
1673 (*j)->setCorrectedDimensions(correcteddims);
1674 correcteddims.clear();
1681 void File::handle_grid_cf_attrs() throw(Exception) {
1688 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1689 i != this->grids.end(); ++i){
1690 for (vector<Field *>::const_iterator j =
1691 (*i)->getDataFields().begin();
1692 j != (*i)->getDataFields().end(); ++j) {
1695 if((*j)->fieldtype == 0) {
1696 string tempcoordinates=
"";
1697 string tempfieldname=
"";
1698 string tempcorrectedfieldname=
"";
1700 for(vector<Dimension *>::const_iterator k=(*j)->getDimensions().begin();
1701 k!=(*j)->getDimensions().end();++k){
1704 map<string,string>::iterator tempmapit;
1705 map<string,string>::iterator tempmapit2;
1708 tempmapit = ((*i)->dimcvarlist).find((*k)->getName());
1709 if(tempmapit != ((*i)->dimcvarlist).end())
1710 tempfieldname = tempmapit->second;
1712 throw4(
"cannot find the dimension field name",
1713 (*i)->getName(),(*j)->getName(),(*k)->getName());
1716 tempmapit2 = ((*i)->ncvarnamelist).find(tempfieldname);
1717 if(tempmapit2 != ((*i)->ncvarnamelist).end())
1718 tempcorrectedfieldname = tempmapit2->second;
1720 throw4(
"cannot find the corrected dimension field name",
1721 (*i)->getName(),(*j)->getName(),(*k)->getName());
1724 tempcoordinates= tempcorrectedfieldname;
1726 tempcoordinates = tempcoordinates +
" "+tempcorrectedfieldname;
1729 (*j)->setCoordinates(tempcoordinates);
1733 if((*j)->fieldtype == 1) {
1734 string tempunits =
"degrees_north";
1735 (*j)->setUnits(tempunits);
1737 if((*j)->fieldtype == 2) {
1738 string tempunits =
"degrees_east";
1739 (*j)->setUnits(tempunits);
1746 if((*j)->fieldtype == 4) {
1747 string tempunits =
"level";
1748 (*j)->setUnits(tempunits);
1752 if((*j)->fieldtype == 5) {
1753 string tempunits =
"days since 1900-01-01 00:00:00";
1754 (*j)->setUnits(tempunits);
1762 if (
true == (*i)->addfvalueattr) {
1763 if((((*j)->getFillValue()).empty()) && ((*j)->getType()==DFNT_FLOAT32 )) {
1764 float tempfillvalue = HUGE;
1765 (*j)->addFillValue(tempfillvalue);
1766 (*j)->setAddedFillValue(
true);
1774 void File::handle_grid_SOM_projection() throw(Exception) {
1783 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
1784 i != this->grids.end(); ++i){
1785 if (GCTP_SOM == (*i)->getProjection().getCode()) {
1791 for(vector<Dimension *>::const_iterator j=(*i)->getDimensions().begin();
1792 j!=(*i)->getDimensions().end();++j){
1795 if(NBLOCK == (*j)->getSize()) {
1799 if ((*j)->getName().compare(0,3,
"SOM") == 0) {
1800 som_dimname = (*j)->getName();
1806 if(
""== som_dimname)
1807 throw4(
"Wrong number of block: The number of block of MISR SOM Grid ",
1808 (*i)->getName(),
" is not ",NBLOCK);
1810 map<string,string>::iterator tempmapit;
1813 string cor_som_dimname;
1814 tempmapit = (*i)->ndimnamelist.find(som_dimname);
1815 if(tempmapit != (*i)->ndimnamelist.end())
1816 cor_som_dimname = tempmapit->second;
1818 throw2(
"cannot find the corrected dimension name for ", som_dimname);
1821 string cor_som_cvname;
1824 for (vector<Field *>::iterator j = (*i)->datafields.begin();
1825 j != (*i)->datafields.end(); ) {
1829 if (1 == (*j)->fieldtype || 2 == (*j)->fieldtype) {
1831 Dimension *newdim =
new Dimension(som_dimname,NBLOCK);
1832 Dimension *newcor_dim =
new Dimension(cor_som_dimname,NBLOCK);
1833 vector<Dimension *>::iterator it_d;
1835 it_d = (*j)->dims.begin();
1836 (*j)->dims.insert(it_d,newdim);
1838 it_d = (*j)->correcteddims.begin();
1839 (*j)->correcteddims.insert(it_d,newcor_dim);
1848 if ( 4 == (*j)->fieldtype) {
1849 cor_som_cvname = (*j)->newname;
1851 j = (*i)->datafields.erase(j);
1866 for (vector<Field *>::const_iterator j = (*i)->getDataFields().begin();
1867 j != (*i)->getDataFields().end(); ++j) {
1869 if ( 0 == (*j)->fieldtype) {
1871 string temp_coordinates = (*j)->coordinates;
1874 found = temp_coordinates.find(cor_som_cvname);
1878 if (temp_coordinates.size() >cor_som_cvname.size())
1879 temp_coordinates.erase(found,cor_som_cvname.size()+1);
1881 temp_coordinates.erase(found,cor_som_cvname.size());
1883 else if (found != string::npos)
1884 temp_coordinates.erase(found-1,cor_som_cvname.size()+1);
1886 throw4(
"cannot find the coordinate variable ",cor_som_cvname,
1887 "from ",temp_coordinates);
1889 (*j)->setCoordinates(temp_coordinates);
1898 int File::obtain_dimmap_num(
int numswath)
throw(Exception) {
1902 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
1903 i != this->swaths.end(); ++i){
1904 tempnumdm += (*i)->get_num_map();
1916 bool fakedimmap =
false;
1920 if((this->swaths[0]->getName()).find(
"atml2")!=string::npos){
1926 for (vector<Field *>::const_iterator j =
1927 this->swaths[0]->getGeoFields().begin();
1928 j != this->swaths[0]->getGeoFields().end(); ++j) {
1929 if((*j)->getName() ==
"Latitude" || (*j)->getName() ==
"Longitude") {
1930 if ((*j)->getType() == DFNT_UINT16 ||
1931 (*j)->getType() == DFNT_INT16)
1932 (*j)->type = DFNT_FLOAT32;
1941 for (vector<Field *>::const_iterator j =
1942 this->swaths[0]->getDataFields().begin();
1943 j != this->swaths[0]->getDataFields().end(); ++j) {
1958 if(((*j)->getName()).find(
"Latitude") != string::npos){
1960 if ((*j)->getType() == DFNT_UINT16 ||
1961 (*j)->getType() == DFNT_INT16)
1962 (*j)->type = DFNT_FLOAT32;
1964 (*j)->fieldtype = 1;
1967 if((*j)->getRank() != 2)
1968 throw2(
"The lat/lon rank must be 2 for Java clients to work",
1971 (((*j)->getDimensions())[0])->getName(),(*j)->getName());
1975 if(((*j)->getName()).find(
"Longitude")!= string::npos) {
1977 if((*j)->getType() == DFNT_UINT16 ||
1978 (*j)->getType() == DFNT_INT16)
1979 (*j)->type = DFNT_FLOAT32;
1981 (*j)->fieldtype = 2;
1982 if((*j)->getRank() != 2)
1983 throw2(
"The lat/lon rank must be 2 for Java clients to work",
1986 (((*j)->getDimensions())[1])->getName(), (*j)->getName());
1998 if(
true == fakedimmap)
2007 void File::create_swath_latlon_dim_cvar_map(
int numdm)
throw(Exception){
2022 bool lat_in_geofields =
false;
2023 bool lon_in_geofields =
false;
2025 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2026 i != this->swaths.end(); ++i){
2028 int tempgeocount = 0;
2029 for (vector<Field *>::const_iterator j =
2030 (*i)->getGeoFields().begin();
2031 j != (*i)->getGeoFields().end(); ++j) {
2035 if((*j)->getName()==
"Latitude" ){
2036 if((*j)->getRank() > 2)
2037 throw2(
"Currently the lat/lon rank must be 1 or 2 for Java clients to work",
2040 lat_in_geofields =
true;
2054 for(vector<SwathDataset::DimensionMap *>::const_iterator
2055 l=(*i)->getDimensionMaps().begin(); l!=(*i)->getDimensionMaps().end();++l){
2059 if(((*j)->getDimensions()[0])->getName() == (*l)->getGeoDimension()) {
2065 (*j)->fieldtype = 1;
2069 if((*j)->getName()==
"Longitude"){
2070 if((*j)->getRank() > 2)
2071 throw2(
"Currently the lat/lon rank must be 1 or 2 for Java clients to work",
2076 lon_in_geofields =
true;
2077 if((*j)->getRank() == 1) {
2087 (((*j)->getDimensions())[1])->getName(),
"Longitude");
2091 for(vector<SwathDataset::DimensionMap *>::const_iterator
2092 l=(*i)->getDimensionMaps().begin(); l!=(*i)->getDimensionMaps().end();++l){
2097 if(((*j)->getDimensions()[1])->getName() ==
2098 (*l)->getGeoDimension()) {
2100 (*l)->getDataDimension(),
"Longitude");
2105 (*j)->fieldtype = 2;
2108 if(tempgeocount == 2)
2115 if (lat_in_geofields!=lon_in_geofields)
2116 throw1(
"Latitude and longitude must be both under Geolocation fields or Data fields");
2119 if (!lat_in_geofields && !lon_in_geofields) {
2121 bool lat_in_datafields =
false;
2122 bool lon_in_datafields =
false;
2124 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2125 i != this->swaths.end(); ++i){
2127 int tempgeocount = 0;
2128 for (vector<Field *>::const_iterator j =
2129 (*i)->getDataFields().begin();
2130 j != (*i)->getDataFields().end(); ++j) {
2136 if((*j)->getName()==
"Latitude" ){
2137 if((*j)->getRank() > 2) {
2138 throw2(
"Currently the lat/lon rank must be 1 or 2 for Java clients to work",
2141 lat_in_datafields =
true;
2150 (((*j)->getDimensions())[0])->getName(),
"Latitude");
2156 for(vector<SwathDataset::DimensionMap *>::const_iterator
2157 l=(*i)->getDimensionMaps().begin(); l!=(*i)->getDimensionMaps().end();++l){
2161 if(((*j)->getDimensions()[0])->getName() == (*l)->getGeoDimension()) {
2167 (*j)->fieldtype = 1;
2171 if((*j)->getName()==
"Longitude"){
2173 if((*j)->getRank() > 2) {
2174 throw2(
"Currently the lat/lon rank must be 1 or 2 for Java clients to work",
2180 lon_in_datafields =
true;
2181 if((*j)->getRank() == 1) {
2191 (((*j)->getDimensions())[1])->getName(),
"Longitude");
2194 for(vector<SwathDataset::DimensionMap *>::const_iterator
2195 l=(*i)->getDimensionMaps().begin(); l!=(*i)->getDimensionMaps().end();++l){
2198 if(((*j)->getDimensions()[1])->getName() == (*l)->getGeoDimension()) {
2200 (*l)->getDataDimension(),
"Longitude");
2205 (*j)->fieldtype = 2;
2208 if(tempgeocount == 2)
2215 if (lat_in_datafields!=lon_in_datafields)
2216 throw1(
"Latitude and longitude must be both under Geolocation fields or Data fields");
2224 void File:: create_swath_nonll_dim_cvar_map() throw(Exception)
2227 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2228 i != this->swaths.end(); ++i){
2245 pair<set<string>::iterator,
bool> tempdimret;
2246 for(map<string,string>::const_iterator j = (*i)->dimcvarlist.begin();
2247 j!= (*i)->dimcvarlist.end();++j){
2248 tempdimret = (*i)->nonmisscvdimlist.insert((*j).first);
2255 for (vector<Field *>::const_iterator j =
2256 (*i)->getGeoFields().begin();
2257 j != (*i)->getGeoFields().end(); ++j) {
2259 if((*j)->getRank()==1) {
2260 if((*i)->nonmisscvdimlist.find((((*j)->getDimensions())[0])->getName()) == (*i)->nonmisscvdimlist.end()){
2261 tempdimret = (*i)->nonmisscvdimlist.insert((((*j)->getDimensions())[0])->getName());
2262 if((*j)->getName() ==
"Time")
2263 (*j)->fieldtype = 5;
2273 if(((((*j)->getDimensions())[0])->getName())==(*j)->getName()){
2274 (*j)->oriname = (*j)->getName();
2277 (*j)->specialcoard =
true;
2281 (*j)->fieldtype = 3;
2291 for (vector<Field *>::const_iterator j =
2292 (*i)->getDataFields().begin();
2293 j != (*i)->getDataFields().end(); ++j) {
2295 if((*j)->getRank()==1) {
2296 if((*i)->nonmisscvdimlist.find((((*j)->getDimensions())[0])->getName()) == (*i)->nonmisscvdimlist.end()){
2297 tempdimret = (*i)->nonmisscvdimlist.insert((((*j)->getDimensions())[0])->getName());
2298 if((*j)->getName() ==
"Time")
2299 (*j)->fieldtype = 5;
2306 if(((((*j)->getDimensions())[0])->getName())==(*j)->getName()){
2307 (*j)->oriname = (*j)->getName();
2309 (*j)->specialcoard =
true;
2313 (*j)->fieldtype = 3;
2323 bool missingfield_unlim_flag =
false;
2324 Field *missingfield_unlim = NULL;
2326 for (vector<Dimension *>::const_iterator j =
2327 (*i)->getDimensions().begin(); j!= (*i)->getDimensions().end();++j){
2329 if(((*i)->nonmisscvdimlist.find((*j)->getName())) == (*i)->nonmisscvdimlist.end()){
2332 Field *missingfield =
new Field();
2341 missingfield->name = (*j)->getName();
2342 missingfield->rank = 1;
2343 missingfield->type = DFNT_INT32;
2344 Dimension *dim =
new Dimension((*j)->getName(),(*j)->getSize());
2347 missingfield->dims.push_back(dim);
2351 int missingdimsize[1];
2352 missingdimsize[0]= (*j)->getSize();
2354 if(0 == (*j)->getSize()) {
2355 missingfield_unlim_flag =
true;
2356 missingfield_unlim = missingfield;
2360 missingfield->fieldtype = 4;
2362 (*i)->geofields.push_back(missingfield);
2364 (missingfield->getDimensions())[0]->getName(), missingfield->name);
2374 bool temp_missingfield_unlim_flag = missingfield_unlim_flag;
2375 if(
true == temp_missingfield_unlim_flag) {
2376 for (vector<Field *>::const_iterator j =
2377 (*i)->getDataFields().begin();
2378 j != (*i)->getDataFields().end(); ++j) {
2380 for (vector<Dimension *>::const_iterator k =
2381 (*j)->getDimensions().begin(); k!= (*j)->getDimensions().end();++k){
2383 if((*k)->getName() == (missingfield_unlim->getDimensions())[0]->getName()) {
2384 if((*k)->getSize()!= 0) {
2385 Dimension *dim = missingfield_unlim->getDimensions()[0];
2387 dim->dimsize = (*k)->getSize();
2388 missingfield_unlim_flag =
false;
2394 if(
false == missingfield_unlim_flag)
2399 (*i)->nonmisscvdimlist.clear();
2407 void File::handle_swath_dim_cvar_maps(
int tempnumdm)
throw(Exception) {
2410 vector <string> tempfieldnamelist;
2411 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2412 i != this->swaths.end(); ++i) {
2415 for (vector<Field *>::const_iterator j =
2416 (*i)->getGeoFields().begin();
2417 j != (*i)->getGeoFields().end(); ++j) {
2421 for (vector<Field *>::const_iterator j = (*i)->getDataFields().begin();
2422 j!= (*i)->getDataFields().end(); ++j) {
2429 int total_fcounter = 0;
2434 map<string,string>tempncvarnamelist;
2435 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2436 i != this->swaths.end(); ++i){
2439 for (vector<Field *>::const_iterator j =
2440 (*i)->getGeoFields().begin();
2441 j != (*i)->getGeoFields().end(); ++j)
2444 (*j)->newname = tempfieldnamelist[total_fcounter];
2448 if((*j)->fieldtype!=0) {
2453 for (vector<Field *>::const_iterator j =
2454 (*i)->getDataFields().begin();
2455 j != (*i)->getDataFields().end(); ++j)
2457 (*j)->newname = tempfieldnamelist[total_fcounter];
2461 if((*j)->fieldtype!=0) {
2469 vector <string>tempalldimnamelist;
2471 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2472 i != this->swaths.end(); ++i)
2473 for (map<string,string>::const_iterator j =
2474 (*i)->dimcvarlist.begin(); j!= (*i)->dimcvarlist.end();++j)
2480 int total_dcounter = 0;
2482 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2483 i != this->swaths.end(); ++i){
2484 for (map<string,string>::const_iterator j =
2485 (*i)->dimcvarlist.begin(); j!= (*i)->dimcvarlist.end();++j){
2492 vector<Dimension*> correcteddims;
2493 string tempcorrecteddimname;
2494 Dimension *correcteddim;
2496 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2497 i != this->swaths.end(); ++i){
2500 for (vector<Field *>::const_iterator j =
2501 (*i)->getGeoFields().begin();
2502 j != (*i)->getGeoFields().end(); ++j) {
2504 for(vector<Dimension *>::const_iterator k=(*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k){
2506 map<string,string>::iterator tempmapit;
2508 if(tempnumdm == 0) {
2511 tempmapit = (*i)->ndimnamelist.find((*k)->getName());
2512 if(tempmapit != (*i)->ndimnamelist.end())
2513 tempcorrecteddimname= tempmapit->second;
2515 throw4(
"cannot find the corrected dimension name",
2516 (*i)->getName(),(*j)->getName(),(*k)->getName());
2518 correcteddim =
new Dimension(tempcorrecteddimname,(*k)->getSize());
2522 bool isdimmapname =
false;
2525 for(vector<SwathDataset::DimensionMap *>::const_iterator
2526 l=(*i)->getDimensionMaps().begin(); l!=(*i)->getDimensionMaps().end();++l){
2529 if((*k)->getName() == (*l)->getGeoDimension()) {
2531 isdimmapname =
true;
2533 string temprepdimname = (*l)->getDataDimension();
2536 tempmapit = (*i)->ndimnamelist.find(temprepdimname);
2537 if(tempmapit != (*i)->ndimnamelist.end())
2538 tempcorrecteddimname= tempmapit->second;
2540 throw4(
"cannot find the corrected dimension name", (*i)->getName(),
2541 (*j)->getName(),(*k)->getName());
2545 bool ddimsflag =
false;
2546 for(vector<Dimension *>::const_iterator m=(*i)->getDimensions().begin();
2547 m!=(*i)->getDimensions().end();++m) {
2548 if((*m)->getName() == temprepdimname) {
2550 correcteddim =
new Dimension(tempcorrecteddimname,(*m)->getSize());
2556 throw4(
"cannot find the corrected dimension size", (*i)->getName(),
2557 (*j)->getName(),(*k)->getName());
2561 if(
false == isdimmapname) {
2563 tempmapit = (*i)->ndimnamelist.find((*k)->getName());
2564 if(tempmapit != (*i)->ndimnamelist.end())
2565 tempcorrecteddimname= tempmapit->second;
2567 throw4(
"cannot find the corrected dimension name",
2568 (*i)->getName(),(*j)->getName(),(*k)->getName());
2570 correcteddim =
new Dimension(tempcorrecteddimname,(*k)->getSize());
2575 correcteddims.push_back(correcteddim);
2577 (*j)->setCorrectedDimensions(correcteddims);
2578 correcteddims.clear();
2582 for (vector<Field *>::const_iterator j =
2583 (*i)->getDataFields().begin();
2584 j != (*i)->getDataFields().end(); ++j) {
2586 for(vector<Dimension *>::const_iterator k=
2587 (*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k){
2589 if(tempnumdm == 0) {
2591 map<string,string>::iterator tempmapit;
2593 tempmapit = (*i)->ndimnamelist.find((*k)->getName());
2594 if(tempmapit != (*i)->ndimnamelist.end())
2595 tempcorrecteddimname= tempmapit->second;
2597 throw4(
"cannot find the corrected dimension name", (*i)->getName(),
2598 (*j)->getName(),(*k)->getName());
2600 correcteddim =
new Dimension(tempcorrecteddimname,(*k)->getSize());
2603 map<string,string>::iterator tempmapit;
2604 bool isdimmapname =
false;
2606 for(vector<SwathDataset::DimensionMap *>::const_iterator l=
2607 (*i)->getDimensionMaps().begin(); l!=(*i)->getDimensionMaps().end();++l){
2610 if((*k)->getName() == (*l)->getGeoDimension()) {
2611 isdimmapname =
true;
2613 string temprepdimname = (*l)->getDataDimension();
2616 tempmapit = (*i)->ndimnamelist.find(temprepdimname);
2617 if(tempmapit != (*i)->ndimnamelist.end())
2618 tempcorrecteddimname= tempmapit->second;
2620 throw4(
"cannot find the corrected dimension name",
2621 (*i)->getName(),(*j)->getName(),(*k)->getName());
2625 bool ddimsflag =
false;
2626 for(vector<Dimension *>::const_iterator m=
2627 (*i)->getDimensions().begin();m!=(*i)->getDimensions().end();++m) {
2630 if((*m)->getName() == temprepdimname) {
2631 correcteddim =
new Dimension(tempcorrecteddimname,(*m)->getSize());
2637 throw4(
"cannot find the corrected dimension size",
2638 (*i)->getName(),(*j)->getName(),(*k)->getName());
2646 tempmapit = (*i)->ndimnamelist.find((*k)->getName());
2647 if(tempmapit != (*i)->ndimnamelist.end())
2648 tempcorrecteddimname= tempmapit->second;
2650 throw4(
"cannot find the corrected dimension name",
2651 (*i)->getName(),(*j)->getName(),(*k)->getName());
2653 correcteddim =
new Dimension(tempcorrecteddimname,(*k)->getSize());
2657 correcteddims.push_back(correcteddim);
2659 (*j)->setCorrectedDimensions(correcteddims);
2660 correcteddims.clear();
2667 void File::handle_swath_cf_attrs() throw(Exception) {
2677 for (vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2678 i != this->swaths.end(); ++i){
2681 for (vector<Field *>::const_iterator j =
2682 (*i)->getGeoFields().begin();
2683 j != (*i)->getGeoFields().end(); ++j) {
2686 if((*j)->fieldtype == 0) {
2687 string tempcoordinates=
"";
2688 string tempfieldname=
"";
2689 string tempcorrectedfieldname=
"";
2691 for(vector<Dimension *>::const_iterator
2692 k=(*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k){
2695 map<string,string>::iterator tempmapit;
2696 map<string,string>::iterator tempmapit2;
2699 tempmapit = ((*i)->dimcvarlist).find((*k)->getName());
2700 if(tempmapit != ((*i)->dimcvarlist).end())
2701 tempfieldname = tempmapit->second;
2703 throw4(
"cannot find the dimension field name",(*i)->getName(),
2704 (*j)->getName(),(*k)->getName());
2707 tempmapit2 = ((*i)->ncvarnamelist).find(tempfieldname);
2708 if(tempmapit2 != ((*i)->ncvarnamelist).end())
2709 tempcorrectedfieldname = tempmapit2->second;
2711 throw4(
"cannot find the corrected dimension field name",
2712 (*i)->getName(),(*j)->getName(),(*k)->getName());
2715 tempcoordinates= tempcorrectedfieldname;
2717 tempcoordinates = tempcoordinates +
" "+tempcorrectedfieldname;
2720 (*j)->setCoordinates(tempcoordinates);
2725 if((*j)->fieldtype == 1) {
2726 string tempunits =
"degrees_north";
2727 (*j)->setUnits(tempunits);
2731 if((*j)->fieldtype == 2) {
2732 string tempunits =
"degrees_east";
2733 (*j)->setUnits(tempunits);
2740 if((*j)->fieldtype == 4) {
2741 string tempunits =
"level";
2742 (*j)->setUnits(tempunits);
2747 if((*j)->fieldtype == 5) {
2748 string tempunits =
"days since 1900-01-01 00:00:00";
2749 (*j)->setUnits(tempunits);
2755 if((((*j)->getFillValue()).empty()) &&
2756 ((*j)->getType()==DFNT_FLOAT32 || (*j)->getType()==DFNT_FLOAT64)) {
2757 float tempfillvalue = -9999.0;
2758 (*j)->addFillValue(tempfillvalue);
2759 (*j)->setAddedFillValue(
true);
2764 for (vector<Field *>::const_iterator j =
2765 (*i)->getDataFields().begin();
2766 j != (*i)->getDataFields().end(); ++j) {
2769 if((*j)->fieldtype == 0) {
2770 string tempcoordinates=
"";
2771 string tempfieldname=
"";
2772 string tempcorrectedfieldname=
"";
2774 for(vector<Dimension *>::const_iterator k
2775 =(*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k){
2778 map<string,string>::iterator tempmapit;
2779 map<string,string>::iterator tempmapit2;
2782 tempmapit = ((*i)->dimcvarlist).find((*k)->getName());
2783 if(tempmapit != ((*i)->dimcvarlist).end())
2784 tempfieldname = tempmapit->second;
2786 throw4(
"cannot find the dimension field name",(*i)->getName(),
2787 (*j)->getName(),(*k)->getName());
2790 tempmapit2 = ((*i)->ncvarnamelist).find(tempfieldname);
2791 if(tempmapit2 != ((*i)->ncvarnamelist).end())
2792 tempcorrectedfieldname = tempmapit2->second;
2794 throw4(
"cannot find the corrected dimension field name",
2795 (*i)->getName(),(*j)->getName(),(*k)->getName());
2798 tempcoordinates= tempcorrectedfieldname;
2800 tempcoordinates = tempcoordinates +
" "+tempcorrectedfieldname;
2803 (*j)->setCoordinates(tempcoordinates);
2806 if(((*j)->fieldtype == 3)||((*j)->fieldtype == 4)) {
2807 string tempunits =
"level";
2808 (*j)->setUnits(tempunits);
2813 if((*j)->fieldtype == 5) {
2814 string tempunits =
"days since 1900-01-01 00:00:00";
2815 (*j)->setUnits(tempunits);
2822 if((((*j)->getFillValue()).empty()) &&
2823 ((*j)->getType()==DFNT_FLOAT32 || (*j)->getType()==DFNT_FLOAT64)) {
2824 float tempfillvalue = -9999.0;
2825 (*j)->addFillValue(tempfillvalue);
2826 (*j)->setAddedFillValue(
true);
2835 void File::Prepare(
const char *eosfile_path)
throw(Exception)
2843 int numgrid = this->grids.size();
2844 int numswath = this->swaths.size();
2847 throw2(
"the number of grid is less than 0", eosfile_path);
2853 string DIMXNAME = this->get_geodim_x_name();
2855 string DIMYNAME = this->get_geodim_y_name();
2857 string LATFIELDNAME = this->get_latfield_name();
2859 string LONFIELDNAME = this->get_lonfield_name();
2862 string GEOGRIDNAME =
"location";
2867 check_onelatlon_grids();
2870 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
2871 i != this->grids.end(); ++i) {
2872 handle_one_grid_zdim(*i);
2876 if (
true == this->onelatlon)
2877 handle_onelatlon_grids();
2879 for (vector<GridDataset *>::const_iterator i = this->grids.begin();
2880 i != this->grids.end(); ++i) {
2884 (*i)->setDimxName(DIMXNAME);
2885 (*i)->setDimyName(DIMYNAME);
2888 handle_one_grid_latlon(*i);
2893 handle_grid_dim_cvar_maps();
2896 handle_grid_coards();
2899 update_grid_field_corrected_dims();
2902 handle_grid_cf_attrs();
2905 handle_grid_SOM_projection();
2911 for(vector<GridDataset *>::const_iterator i = this->grids.begin();
2912 i != this->grids.end(); ++i){
2913 (*i)->SetScaleType((*i)->name);
2922 int tempnumdm = obtain_dimmap_num(numswath);
2925 create_swath_latlon_dim_cvar_map(tempnumdm);
2928 create_swath_nonll_dim_cvar_map();
2931 handle_swath_dim_cvar_maps(tempnumdm);
2935 handle_swath_cf_attrs();
2938 for(vector<SwathDataset *>::const_iterator i = this->swaths.begin();
2939 i != this->swaths.end(); ++i)
2940 (*i)->SetScaleType((*i)->name);
2948 void correct_unlimited_missing_zdim(GridDataset* gdset)
throw(Exception) {
2950 for (vector<Field *>::const_iterator j =
2951 gdset->getDataFields().begin();
2952 j != gdset->getDataFields().end(); ++j) {
2955 if ((*j)->getRank()==1 && (*j)->){
2965 bool File::check_special_1d_grid() throw(Exception) {
2967 int numgrid = this->grids.size();
2968 int numswath = this->swaths.size();
2970 if (numgrid != 1 || numswath != 0)
2974 string DIMXNAME = this->get_geodim_x_name();
2975 string DIMYNAME = this->get_geodim_y_name();
2977 if(DIMXNAME !=
"XDim" || DIMYNAME !=
"YDim")
2980 int var_dimx_size = 0;
2981 int var_dimy_size = 0;
2983 GridDataset *mygrid = (this->grids)[0];
2985 int field_xydim_flag = 0;
2986 for (vector<Field *>::const_iterator i = mygrid->getDataFields().begin();
2987 i!= mygrid->getDataFields().end(); ++i) {
2989 if((*i)->name ==
"XDim"){
2991 var_dimx_size = ((*i)->getDimensions())[0]->getSize();
2993 if((*i)->name ==
"YDim"){
2995 var_dimy_size = ((*i)->getDimensions())[0]->getSize();
2998 if(2==field_xydim_flag)
3002 if(field_xydim_flag !=2)
3006 int xdimsize = mygrid->getInfo().getX();
3007 int ydimsize = mygrid->getInfo().getY();
3009 if(var_dimx_size != xdimsize || var_dimy_size != ydimsize)
3028 void Dataset::SetScaleType(
const string EOS2ObjName)
throw(Exception) {
3036 vector<string> modis_multi_scale_type;
3037 modis_multi_scale_type.push_back(
"L1B");
3038 modis_multi_scale_type.push_back(
"GEO");
3039 modis_multi_scale_type.push_back(
"BRDF");
3040 modis_multi_scale_type.push_back(
"0.05Deg");
3041 modis_multi_scale_type.push_back(
"Reflectance");
3042 modis_multi_scale_type.push_back(
"MOD17A2");
3043 modis_multi_scale_type.push_back(
"North");
3044 modis_multi_scale_type.push_back(
"South");
3045 modis_multi_scale_type.push_back(
"MOD_Swath_Sea_Ice");
3046 modis_multi_scale_type.push_back(
"MOD_Grid_MOD15A2");
3047 modis_multi_scale_type.push_back(
"MOD_Grid_MOD16A2");
3048 modis_multi_scale_type.push_back(
"MOD_Grid_MOD16A3");
3049 modis_multi_scale_type.push_back(
"MODIS_NACP_LAI");
3051 vector<string> modis_div_scale_type;
3054 modis_div_scale_type.push_back(
"VI");
3055 modis_div_scale_type.push_back(
"1km_2D");
3056 modis_div_scale_type.push_back(
"L2g_2d");
3057 modis_div_scale_type.push_back(
"CMG");
3058 modis_div_scale_type.push_back(
"MODIS SWATH TYPE L2");
3063 string modis_eq_scale_type =
"LST";
3064 string modis_equ_scale_lst_group1=
"MODIS_Grid_8Day_1km_LST21";
3065 string modis_equ_scale_lst_group2=
"MODIS_Grid_Daily_1km_LST21";
3067 string modis_divequ_scale_group =
"MODIS_Grid";
3068 string modis_div_scale_group =
"MOD_Grid";
3072 string modis_equ_scale_group =
"MODIS_Grid_1km_2D";
3074 if(EOS2ObjName==
"mod05" || EOS2ObjName==
"mod06" || EOS2ObjName==
"mod07"
3075 || EOS2ObjName==
"mod08" || EOS2ObjName==
"atml2")
3077 scaletype = MODIS_MUL_SCALE;
3091 if(EOS2ObjName.find(
"MOD")==0 || EOS2ObjName.find(
"mod")==0)
3093 size_t pos = EOS2ObjName.rfind(modis_eq_scale_type);
3094 if(pos != string::npos &&
3095 (pos== (EOS2ObjName.length()-modis_eq_scale_type.length())))
3097 scaletype = MODIS_EQ_SCALE;
3101 for(
unsigned int k=0; k<modis_multi_scale_type.size(); k++)
3103 pos = EOS2ObjName.rfind(modis_multi_scale_type[k]);
3104 if (pos !=string::npos &&
3105 (pos== (EOS2ObjName.length()-modis_multi_scale_type[k].length())))
3107 scaletype = MODIS_MUL_SCALE;
3112 for(
unsigned int k=0; k<modis_div_scale_type.size(); k++)
3114 pos = EOS2ObjName.rfind(modis_div_scale_type[k]);
3115 if (pos != string::npos &&
3116 (pos==(EOS2ObjName.length()-modis_div_scale_type[k].length()))){
3117 scaletype = MODIS_DIV_SCALE;
3122 if (EOS2ObjName !=
"MODIS_Grid_1km_2D")
3129 pos = EOS2ObjName.find(modis_divequ_scale_group);
3135 size_t eq_scale_pos = EOS2ObjName.find(modis_equ_scale_group)
3136 *EOS2ObjName.find(modis_equ_scale_lst_group1)
3137 *EOS2ObjName.find(modis_equ_scale_lst_group2);
3138 if (0 == eq_scale_pos)
3139 scaletype = MODIS_EQ_SCALE;
3141 scaletype = MODIS_DIV_SCALE;
3144 size_t div_scale_pos = EOS2ObjName.find(modis_div_scale_group);
3147 if ( 0 == div_scale_pos)
3148 scaletype = MODIS_DIV_SCALE;
3154 if (EOS2ObjName ==
"VIP_CMG_GRID")
3155 scaletype = MODIS_DIV_SCALE;
3163 for_each(this->dims.begin(), this->dims.end(), delete_elem());
3164 for_each(this->correcteddims.begin(), this->correcteddims.end(), delete_elem());
3170 for_each(this->dims.begin(), this->dims.end(), delete_elem());
3171 for_each(this->datafields.begin(), this->datafields.end(),
3173 for_each(this->attrs.begin(), this->attrs.end(), delete_elem());
3177 void Dataset::ReadDimensions(int32 (*entries)(int32, int32, int32 *),
3178 int32 (*inq)(int32,
char *, int32 *),
3179 vector<Dimension *> &d_dims)
throw(Exception)
3189 if ((numdims = entries(this->datasetid, HDFE_NENTDIM, &bufsize)) == -1)
3190 throw2(
"dimension entry", this->name);
3194 vector<char> namelist;
3195 vector<int32> dimsize;
3197 namelist.resize(bufsize + 1);
3198 dimsize.resize(numdims);
3201 if (inq(this->datasetid, &namelist[0], &dimsize[0]) == -1)
3202 throw2(
"inquire dimension", this->name);
3204 vector<string> dimnames;
3210 for (vector<string>::const_iterator i = dimnames.begin();
3211 i != dimnames.end(); ++i) {
3212 Dimension *dim =
new Dimension(*i, dimsize[count]);
3213 d_dims.push_back(dim);
3220 void Dataset::ReadFields(int32 (*entries)(int32, int32, int32 *),
3221 int32 (*inq)(int32,
char *, int32 *, int32 *),
3223 (int32,
char *, int32 *, int32 *, int32 *,
char *),
3225 (int32,
char *, int32 *, int32 *, int32 *, VOIDP),
3226 intn (*getfill)(int32,
char *, VOIDP),
3228 vector<Field *> &fields)
throw(Exception)
3232 int32 numfields = 0;
3239 if ((numfields = entries(this->datasetid, geofield ?
3240 HDFE_NENTGFLD : HDFE_NENTDFLD, &bufsize)) == -1)
3241 throw2(
"field entry", this->name);
3245 if (numfields > 0) {
3246 vector<char> namelist;
3248 namelist.resize(bufsize + 1);
3251 if (inq(this->datasetid, &namelist[0], NULL, NULL) == -1)
3252 throw2(
"inquire field", this->name);
3254 vector<string> fieldnames;
3259 for (vector<string>::const_iterator i = fieldnames.begin();
3260 i != fieldnames.end(); ++i) {
3262 Field *field =
new Field();
3264 throw1(
"The field is NULL");
3267 bool throw_error =
false;
3276 if ((fldinfo(this->datasetid,
3277 const_cast<char *
>(field->name.c_str()),
3278 &field->rank, dimsize, &field->type, dimlist)) == -1){
3279 string fieldname_for_eh = field->name;
3281 err_msg =
"Obtain field info error for field name "+fieldname_for_eh;
3284 if(
false == throw_error) {
3286 vector<string> dimnames;
3290 if ((
int)dimnames.size() != field->rank) {
3292 err_msg =
"Dimension names size is not consistent with field rank. ";
3293 err_msg +=
"Field name is "+field->name;
3296 for (
int k = 0; k < field->rank; ++k) {
3297 Dimension *dim =
new Dimension(dimnames[k], dimsize[k]);
3298 field->dims.push_back(dim);
3302 field->filler.resize(DFKNTsize(field->type));
3303 if (getfill(this->datasetid,
3304 const_cast<char *
>(field->name.c_str()),
3305 &field->filler[0]) == -1)
3306 field->filler.clear();
3309 fields.push_back(field);
3313 if(
true == throw_error) {
3323 void Dataset::ReadAttributes(int32 (*inq)(int32,
char *, int32 *),
3324 intn (*attrinfo)(int32,
char *, int32 *, int32 *),
3325 intn (*readattr)(int32,
char *, VOIDP),
3326 vector<Attribute *> &obj_attrs)
throw(Exception)
3335 if ((numattrs = inq(this->datasetid, NULL, &bufsize)) == -1)
3336 throw2(
"inquire attribute", this->name);
3340 vector<char> namelist;
3342 namelist.resize(bufsize + 1);
3345 if (inq(this->datasetid, &namelist[0], &bufsize) == -1)
3346 throw2(
"inquire attribute", this->name);
3348 vector<string> attrnames;
3353 for (vector<string>::const_iterator i = attrnames.begin();
3354 i != attrnames.end(); ++i) {
3355 Attribute *attr =
new Attribute();
3362 if (attrinfo(this->datasetid,
const_cast<char *
>
3363 (attr->name.c_str()), &attr->type, &count) == -1) {
3365 throw3(
"attribute info", this->name, attr->name);
3368 attr->count = count/DFKNTsize(attr->type);
3369 attr->value.resize(count);
3376 if (readattr(this->datasetid,
3377 const_cast<char *
>(attr->name.c_str()),
3378 &attr->value[0]) == -1) {
3380 throw3(
"read attribute", this->name, attr->name);
3384 obj_attrs.push_back(attr);
3390 GridDataset::~GridDataset()
3392 if (this->datasetid != -1){
3393 GDdetach(this->datasetid);
3398 GridDataset * GridDataset::Read(int32 fd,
const string &gridname)
3402 bool GD_fun_err =
false;
3403 GridDataset *grid =
new GridDataset(gridname);
3406 if ((grid->datasetid = GDattach(fd,
const_cast<char *
>(gridname.c_str())))
3408 err_msg =
"attach grid";
3416 Info &info = grid->info;
3417 if (GDgridinfo(grid->datasetid, &info.xdim, &info.ydim, info.upleft,
3418 info.lowright) == -1) {
3419 err_msg =
"grid info";
3427 Projection &proj = grid->proj;
3428 if (GDprojinfo(grid->datasetid, &proj.code, &proj.zone, &proj.sphere,
3429 proj.param) == -1) {
3430 err_msg =
"projection info";
3434 if (GDpixreginfo(grid->datasetid, &proj.pix) == -1) {
3435 err_msg =
"pixreg info";
3439 if (GDorigininfo(grid->datasetid, &proj.origin) == -1){
3440 err_msg =
"origin info";
3446 if(
true == GD_fun_err){
3448 throw2(err_msg,gridname);
3453 grid->ReadDimensions(GDnentries, GDinqdims, grid->dims);
3456 grid->ReadFields(GDnentries, GDinqfields, GDfieldinfo, GDreadfield,
3457 GDgetfillvalue,
false, grid->datafields);
3460 grid->ReadAttributes(GDinqattrs, GDattrinfo, GDreadattr, grid->attrs);
3470 GridDataset::Calculated & GridDataset::getCalculated()
const
3472 if (this->calculated.grid !=
this)
3473 this->calculated.grid =
this;
3474 return this->calculated;
3477 bool GridDataset::Calculated::isYDimMajor() throw(Exception)
3479 this->DetectMajorDimension();
3480 return this->ydimmajor;
3484 bool GridDataset::Calculated::isOrthogonal() throw(Exception)
3487 this->ReadLongitudeLatitude();
3488 return this->orthogonal;
3492 int GridDataset::Calculated::DetectFieldMajorDimension() throw(Exception)
3497 for (vector<Field *>::const_iterator i =
3498 this->grid->getDataFields().begin();
3499 i != this->grid->getDataFields().end(); ++i) {
3501 int xdimindex = -1, ydimindex = -1, index = 0;
3504 for (vector<Dimension *>::const_iterator j =
3505 (*i)->getDimensions().begin();
3506 j != (*i)->getDimensions().end(); ++j) {
3507 if ((*j)->getName() == this->grid->dimxname)
3509 else if ((*j)->getName() == this->grid->dimyname)
3513 if (xdimindex == -1 || ydimindex == -1)
3516 int major = ydimindex < xdimindex ? 1 : 0;
3527 else if (ym != major)
3528 throw2(
"inconsistent XDim, YDim order", this->grid->getName());
3534 throw2(
"lack of data fields", this->grid->getName());
3539 void GridDataset::Calculated::DetectMajorDimension() throw(Exception)
3547 for (vector<Field *>::const_iterator i =
3548 this->grid->getDataFields().begin();
3549 i != this->grid->getDataFields().end(); ++i) {
3551 int xdimindex = -1, ydimindex = -1, index = 0;
3554 for (vector<Dimension *>::const_iterator j =
3555 (*i)->getDimensions().begin();
3556 j != (*i)->getDimensions().end(); ++j) {
3557 if ((*j)->getName() == this->grid->dimxname)
3559 else if ((*j)->getName() == this->grid->dimyname)
3563 if (xdimindex == -1 || ydimindex == -1)
3568 if(this->grid->getProjection().getCode() == GCTP_LAMAZ)
3571 major = ydimindex < xdimindex ? 1 : 0;
3582 else if (ym != major)
3583 throw2(
"inconsistent XDim, YDim order", this->grid->getName());
3588 throw2(
"lack of data fields", this->grid->getName());
3589 this->ydimmajor = ym != 0;
3598 static bool IsDisjoint(
const vector<Field *> &l,
3599 const vector<Field *> &r)
3601 for (vector<Field *>::const_iterator i = l.begin(); i != l.end(); ++i)
3603 if (find(r.begin(), r.end(), *i) != r.end())
3611 static bool IsDisjoint(vector<pair<Field *, string> > &l,
const vector<Field *> &r)
3613 for (vector<pair<Field *, string> >::const_iterator i =
3614 l.begin(); i != l.end(); ++i) {
3615 if (find(r.begin(), r.end(), i->first) != r.end())
3622 static bool IsSubset(
const vector<Field *> &s,
const vector<Field *> &b)
3624 for (vector<Field *>::const_iterator i = s.begin(); i != s.end(); ++i)
3626 if (find(b.begin(), b.end(), *i) == b.end())
3633 static bool IsSubset(vector<pair<Field *, string> > &s,
const vector<Field *> &b)
3635 for (vector<pair<Field *, string> >::const_iterator i
3636 = s.begin(); i != s.end(); ++i) {
3637 if (find(b.begin(), b.end(), i->first) == b.end())
3645 SwathDataset::~SwathDataset()
3647 if (this->datasetid != -1) {
3648 SWdetach(this->datasetid);
3651 for_each(this->dimmaps.begin(), this->dimmaps.end(), delete_elem());
3652 for_each(this->indexmaps.begin(), this->indexmaps.end(),
3655 for_each(this->geofields.begin(), this->geofields.end(),
3662 SwathDataset * SwathDataset::Read(int32 fd,
const string &swathname)
3665 SwathDataset *swath =
new SwathDataset(swathname);
3667 throw1(
"Cannot allocate HDF5 Swath object");
3670 if ((swath->datasetid = SWattach(fd,
3671 const_cast<char *
>(swathname.c_str())))
3674 throw2(
"attach swath", swathname);
3682 swath->ReadDimensions(SWnentries, SWinqdims, swath->dims);
3685 swath->ReadFields(SWnentries, SWinqdatafields, SWfieldinfo, SWreadfield,
3686 SWgetfillvalue,
false, swath->datafields);
3689 swath->ReadFields(SWnentries, SWinqgeofields, SWfieldinfo, SWreadfield,
3690 SWgetfillvalue,
true, swath->geofields);
3693 swath->ReadAttributes(SWinqattrs, SWattrinfo, SWreadattr, swath->attrs);
3696 swath->set_num_map(swath->ReadDimensionMaps(swath->dimmaps));
3699 swath->ReadIndexMaps(swath->indexmaps);
3711 int SwathDataset::ReadDimensionMaps(vector<DimensionMap *> &swath_dimmaps)
3714 int32 nummaps, bufsize;
3717 if ((nummaps = SWnentries(this->datasetid, HDFE_NENTMAP, &bufsize)) == -1)
3718 throw2(
"dimmap entry", this->name);
3729 vector<char> namelist;
3730 vector<int32> offset, increment;
3732 namelist.resize(bufsize + 1);
3733 offset.resize(nummaps);
3734 increment.resize(nummaps);
3735 if (SWinqmaps(this->datasetid, &namelist[0], &offset[0], &increment[0])
3737 throw2(
"inquire dimmap", this->name);
3739 vector<string> mapnames;
3742 for (vector<string>::const_iterator i = mapnames.begin();
3743 i != mapnames.end(); ++i) {
3744 vector<string> parts;
3746 if (parts.size() != 2)
3747 throw3(
"dimmap part", parts.size(),
3750 DimensionMap *dimmap =
new DimensionMap(parts[0], parts[1],
3753 swath_dimmaps.push_back(dimmap);
3761 void SwathDataset::ReadIndexMaps(vector<IndexMap *> &swath_indexmaps)
3764 int32 numindices, bufsize;
3766 if ((numindices = SWnentries(this->datasetid, HDFE_NENTIMAP, &bufsize)) ==
3768 throw2(
"indexmap entry", this->name);
3769 if (numindices > 0) {
3771 vector<char> namelist;
3773 namelist.resize(bufsize + 1);
3774 if (SWinqidxmaps(this->datasetid, &namelist[0], NULL) == -1)
3775 throw2(
"inquire indexmap", this->name);
3777 vector<string> mapnames;
3779 for (vector<string>::const_iterator i = mapnames.begin();
3780 i != mapnames.end(); ++i) {
3781 IndexMap *indexmap =
new IndexMap();
3782 vector<string> parts;
3784 if (parts.size() != 2)
3785 throw3(
"indexmap part", parts.size(),
3787 indexmap->geo = parts[0];
3788 indexmap->data = parts[1];
3793 SWdiminfo(this->datasetid,
3794 const_cast<char *
>(indexmap->geo.c_str())))
3796 throw3(
"dimension info", this->name, indexmap->geo);
3797 indexmap->indices.resize(dimsize);
3798 if (SWidxmapinfo(this->datasetid,
3799 const_cast<char *
>(indexmap->geo.c_str()),
3800 const_cast<char *
>(indexmap->data.c_str()),
3801 &indexmap->indices[0]) == -1)
3802 throw4(
"indexmap info", this->name, indexmap->geo,
3806 swath_indexmaps.push_back(indexmap);
3812 PointDataset::~PointDataset()
3816 PointDataset * PointDataset::Read(int32 ,
const string &pointname)
3819 PointDataset *point =
new PointDataset(pointname);
3825 bool Utility::ReadNamelist(
const char *path,
3826 int32 (*inq)(
char *,
char *, int32 *),
3827 vector<string> &names)
3829 char *fname =
const_cast<char *
>(path);
3833 if ((numobjs = inq(fname, NULL, &bufsize)) == -1)
return false;
3835 vector<char> buffer;
3836 buffer.resize(bufsize + 1);
3837 if (inq(fname, &buffer[0], &bufsize) == -1)
return false;