23#include "HDF4RequestHandler.h"
33const char *HDFEOS2::File::_geodim_x_names[] = {
"XDim",
"LonDim",
"nlon"};
36const char *HDFEOS2::File::_geodim_y_names[] = {
"YDim",
"LatDim",
"nlat"};
39const char *HDFEOS2::File::_latfield_names[] = {
40 "Latitude",
"Lat",
"YDim",
"LatCenter"
44const char *HDFEOS2::File::_lonfield_names[] = {
45 "Longitude",
"Lon",
"XDim",
"LonCenter"
50const char *HDFEOS2::File::_geogrid_names[] = {
"location"};
52using namespace HDFEOS2;
56template<
typename T,
typename U,
typename V,
typename W,
typename X>
57static void _throw5(
const char *fname,
int line,
int numarg,
58 const T &a1,
const U &a2,
const V &a3,
const W &a4,
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){
126File * 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 (
const auto &grid: gridlist)
153 file->grids.push_back(GridDataset::Read(file->gridfd, grid));
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 (
const auto &swath:swathlist)
177 file->swaths.push_back(SwathDataset::Read(file->swathfd, swath));
181 throw1(
"SwathDataset Read failed.");
188 vector<string> pointlist;
189 if (!Utility::ReadNamelist(file->path.c_str(), PTinqpoint, pointlist)){
191 throw1(
"Point ReadNamelist failed.");
194 for (
const auto&point: pointlist)
195 file->points.push_back(PointDataset::Read(-1, point));
198 if (file->grids.empty() && file->swaths.empty()
199 && file->points.empty()) {
200 Exception e(
"Not an HDF-EOS2 file");
201 e.setFileType(
false);
212string File::get_geodim_x_name()
214 if(!_geodim_x_name.empty())
215 return _geodim_x_name;
216 _find_geodim_names();
217 return _geodim_x_name;
223string File::get_geodim_y_name()
225 if(!_geodim_y_name.empty())
226 return _geodim_y_name;
227 _find_geodim_names();
228 return _geodim_y_name;
238string File::get_latfield_name()
240 if(!_latfield_name.empty())
241 return _latfield_name;
242 _find_latlonfield_names();
243 return _latfield_name;
246string File::get_lonfield_name()
248 if(!_lonfield_name.empty())
249 return _lonfield_name;
250 _find_latlonfield_names();
251 return _lonfield_name;
258string File::get_geogrid_name()
260 if(!_geogrid_name.empty())
261 return _geogrid_name;
262 _find_geogrid_name();
263 return _geogrid_name;
271void File::_find_geodim_names()
273 set<string> geodim_x_name_set;
274 for(
size_t i = 0; i<
sizeof(_geodim_x_names) /
sizeof(
const char *); i++)
275 geodim_x_name_set.insert(_geodim_x_names[i]);
277 set<string> geodim_y_name_set;
278 for(
size_t i = 0; i<
sizeof(_geodim_y_names) /
sizeof(
const char *); i++)
279 geodim_y_name_set.insert(_geodim_y_names[i]);
284 const size_t gs = grids.size();
285 const size_t ss = swaths.size();
286 for(
size_t i=0; ;i++)
288 Dataset *dataset=
nullptr;
290 dataset =
static_cast<Dataset*
>(grids[i]);
292 dataset =
static_cast<Dataset*
>(swaths[i-gs]);
296 const vector<Dimension *>& dims = dataset->getDimensions();
297 for(vector<Dimension*>::const_iterator it = dims.begin();
298 it != dims.end(); ++it)
300 if(geodim_x_name_set.find((*it)->getName()) != geodim_x_name_set.end())
301 _geodim_x_name = (*it)->getName();
302 else if(geodim_y_name_set.find((*it)->getName()) != geodim_y_name_set.end())
303 _geodim_y_name = (*it)->getName();
310 const size_t gs = grids.size();
314 const Dataset *dataset=
nullptr;
315 dataset =
static_cast<Dataset*
>(grids[0]);
317 const vector<Dimension *>& dims = dataset->getDimensions();
318 for (
const auto &dim:dims)
326 if(geodim_x_name_set.find(dim->getName()) != geodim_x_name_set.end())
327 _geodim_x_name = dim->getName();
328 else if(geodim_y_name_set.find(dim->getName()) != geodim_y_name_set.end())
329 _geodim_y_name = dim->getName();
332 if(_geodim_x_name.empty())
333 _geodim_x_name = _geodim_x_names[0];
334 if(_geodim_y_name.empty())
335 _geodim_y_name = _geodim_y_names[0];
343void File::_find_latlonfield_names()
345 set<string> latfield_name_set;
346 for(
size_t i = 0; i<
sizeof(_latfield_names) /
sizeof(
const char *); i++)
347 latfield_name_set.insert(_latfield_names[i]);
349 set<string> lonfield_name_set;
350 for(
size_t i = 0; i<
sizeof(_lonfield_names) /
sizeof(
const char *); i++)
351 lonfield_name_set.insert(_lonfield_names[i]);
353 const size_t gs = grids.size();
354 const size_t ss = swaths.size();
360 for(
size_t i=0;i<1 ;i++)
362 const Dataset *dataset =
nullptr;
363 SwathDataset *sw =
nullptr;
365 dataset =
static_cast<Dataset*
>(grids[i]);
369 dataset =
static_cast<Dataset*
>(sw);
374 const vector<Field *>& fields = dataset->getDataFields();
375 for (
const auto &field:fields)
377 if(latfield_name_set.find(field->getName()) != latfield_name_set.end())
378 _latfield_name = field->getName();
379 else if(lonfield_name_set.find(field->getName()) != lonfield_name_set.end())
380 _lonfield_name = field->getName();
385 const vector<Field *>& geofields = dataset->getDataFields();
386 for(
const auto &gfield:geofields)
388 if(latfield_name_set.find(gfield->getName()) != latfield_name_set.end())
389 _latfield_name = gfield->getName();
390 else if(lonfield_name_set.find(gfield->getName()) != lonfield_name_set.end())
391 _lonfield_name = gfield->getName();
398 if(_latfield_name.empty())
399 _latfield_name = _latfield_names[0];
400 if(_lonfield_name.empty())
401 _lonfield_name = _lonfield_names[0];
409void File::_find_geogrid_name()
411 set<string> geogrid_name_set;
412 for(
size_t i = 0; i<
sizeof(_geogrid_names) /
sizeof(
const char *); i++)
413 geogrid_name_set.insert(_geogrid_names[i]);
415 const size_t gs = grids.size();
416 const size_t ss = swaths.size();
417 for(
size_t i=0; ;i++)
419 const Dataset *dataset =
nullptr;
421 dataset =
static_cast<Dataset*
>(grids[i]);
423 dataset =
static_cast<Dataset*
>(swaths[i-gs]);
427 if(geogrid_name_set.find(dataset->getName()) != geogrid_name_set.end())
428 _geogrid_name = dataset->getName();
430 if(_geogrid_name.empty())
431 _geogrid_name =
"location";
435void File::check_onelatlon_grids() {
438 string LATFIELDNAME = this->get_latfield_name();
439 string LONFIELDNAME = this->get_lonfield_name();
442 string GEOGRIDNAME =
"location";
451 for (
const auto &grid:grids){
454 for (
const auto &field:grid->getDataFields()) {
455 if(grid->getName()==GEOGRIDNAME){
456 if(field->getName()==LATFIELDNAME){
458 grid->latfield = field;
460 if(field->getName()==LONFIELDNAME){
462 grid->lonfield = field;
468 if((field->getName()==LATFIELDNAME)||(field->getName()==LONFIELDNAME)){
469 grid->ownllflag =
true;
477 if(morellcount ==0 && onellcount ==2)
478 this->onelatlon =
true;
482void File::handle_one_grid_zdim(GridDataset* gdset) {
485 string DIMXNAME = this->get_geodim_x_name();
486 string DIMYNAME = this->get_geodim_y_name();
488 bool missingfield_unlim_flag =
false;
489 const Field *missingfield_unlim =
nullptr;
496 set<string> tempdimlist;
497 pair<set<string>::iterator,
bool> tempdimret;
499 for (
const auto &field:gdset->getDataFields()) {
502 if (field->getRank()==1){
506 if((field->getDimensions())[0]->getName()!=DIMXNAME && (field->getDimensions())[0]->getName()!=DIMYNAME){
508 tempdimret = tempdimlist.insert((field->getDimensions())[0]->getName());
513 if(tempdimret.second ==
true) {
517 field->fieldtype = 3;
518 if(field->getName() ==
"Time")
519 field->fieldtype = 5;
534 for (
const auto &gdim:gdset->getDimensions()) {
537 if(gdim->getName()!=DIMXNAME && gdim->getName()!=DIMYNAME){
540 if((tempdimlist.find(gdim->getName())) == tempdimlist.end()){
543 auto missingfield =
new Field();
544 missingfield->name = gdim->getName();
545 missingfield->rank = 1;
548 missingfield->type = DFNT_INT32;
551 auto dim =
new Dimension(gdim->getName(),gdim->getSize());
554 missingfield->dims.push_back(dim);
556 if (0 == gdim->getSize()) {
557 missingfield_unlim_flag =
true;
558 missingfield_unlim = missingfield;
574 missingfield->fieldtype = 4;
584 gdset->datafields.push_back(missingfield);
593 bool temp_missingfield_unlim_flag = missingfield_unlim_flag;
594 if(
true == temp_missingfield_unlim_flag) {
595 for (
unsigned int i =0; i<gdset->getDataFields().size(); i++) {
597 for (
const auto &gdim:gdset->getDimensions()) {
599 if(gdim->getName() == (missingfield_unlim->getDimensions())[0]->getName()) {
600 if(gdim->getSize()!= 0) {
601 Dimension *dim = missingfield_unlim->getDimensions()[0];
604 dim->dimsize = gdim->getSize();
605 missingfield_unlim_flag =
false;
611 if(
false == missingfield_unlim_flag)
619void File::handle_one_grid_latlon(GridDataset* gdset)
throw(Exception)
623 string DIMXNAME = this->get_geodim_x_name();
624 string DIMYNAME = this->get_geodim_y_name();
626 string LATFIELDNAME = this->get_latfield_name();
627 string LONFIELDNAME = this->get_lonfield_name();
631 if(gdset->ownllflag) {
634 for (
const auto &field:gdset->getDataFields()) {
636 if(field->getName() == LATFIELDNAME) {
640 field->fieldtype = 1;
647 if(field->getRank() > 2)
648 throw3(
"The rank of latitude is greater than 2",
649 gdset->getName(),field->getName());
651 if(field->getRank() != 1) {
655 field->ydimmajor = gdset->getCalculated().isYDimMajor();
661 int32 projectioncode = gdset->getProjection().getCode();
662 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
663 field->condenseddim =
true;
673 for (
const auto &dim:field->getDimensions()) {
674 if (dim->getName() == DIMYNAME)
682 if(field->condenseddim) {
686 for (vector<Dimension *>::const_iterator k =
687 field->getDimensions().begin(); k!= field->getDimensions().end();++k){
688 if((*k)->getName() == DIMYNAME) {
696 for (vector<Dimension *>::const_iterator k =
697 field->getDimensions().begin(); k!= field->getDimensions().end();++k){
698 if((*k)->getName() == DIMYNAME) {
711 else if (field->getName() == LONFIELDNAME) {
714 field->fieldtype = 2;
720 if(field->getRank() >2)
721 throw3(
"The rank of Longitude is greater than 2",gdset->getName(),field->getName());
723 if(field->getRank() != 1) {
726 field->ydimmajor = gdset->getCalculated().isYDimMajor();
731 int32 projectioncode = gdset->getProjection().getCode();
732 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
733 field->condenseddim =
true;
743 for (
const auto &dim:field->getDimensions()) {
744 if(dim->getName() == DIMXNAME)
751 if(field->condenseddim) {
755 for (vector<Dimension *>::const_iterator k =
756 field->getDimensions().begin(); k!= field->getDimensions().end();++k){
757 if((*k)->getName() == DIMXNAME) {
764 for (vector<Dimension *>::const_iterator k =
765 field->getDimensions().begin(); k!= field->getDimensions().end();++k){
766 if((*k)->getName() == DIMXNAME) {
784 auto latfield =
new Field();
785 auto lonfield =
new Field();
787 latfield->name = LATFIELDNAME;
788 lonfield->name = LONFIELDNAME;
795 latfield->type = DFNT_FLOAT64;
796 lonfield->type = DFNT_FLOAT64;
799 latfield->fieldtype = 1;
802 lonfield->fieldtype = 2;
806 latfield->ydimmajor = gdset->getCalculated().isYDimMajor();
807 lonfield->ydimmajor = latfield->ydimmajor;
810 int xdimsize = gdset->getInfo().getX();
811 int ydimsize = gdset->getInfo().getY();
816 bool dmajor=(gdset->getProjection().getCode()==GCTP_LAMAZ)? gdset->getCalculated().DetectFieldMajorDimension()
817 : latfield->ydimmajor;
820 auto dimlaty =
new Dimension(DIMYNAME,ydimsize);
821 latfield->dims.push_back(dimlaty);
822 auto dimlony =
new Dimension(DIMYNAME,ydimsize);
823 lonfield->dims.push_back(dimlony);
824 auto dimlatx =
new Dimension(DIMXNAME,xdimsize);
825 latfield->dims.push_back(dimlatx);
826 auto dimlonx =
new Dimension(DIMXNAME,xdimsize);
827 lonfield->dims.push_back(dimlonx);
830 auto dimlatx =
new Dimension(DIMXNAME,xdimsize);
831 latfield->dims.push_back(dimlatx);
832 auto dimlonx =
new Dimension(DIMXNAME,xdimsize);
833 lonfield->dims.push_back(dimlonx);
834 auto dimlaty =
new Dimension(DIMYNAME,ydimsize);
835 latfield->dims.push_back(dimlaty);
836 auto dimlony =
new Dimension(DIMYNAME,ydimsize);
837 lonfield->dims.push_back(dimlony);
844 upleft =
const_cast<float64 *
>(gdset->getInfo().getUpLeft());
845 lowright =
const_cast<float64 *
>(gdset->getInfo().getLowRight());
848 const float64* upleft = gdset->getInfo().getUpLeft();
849 const float64* lowright = gdset->getInfo().getLowRight();
852 int32 projectioncode = gdset->getProjection().getCode();
853 if(((
int)lowright[0]>180000000) && ((
int)upleft[0]>-1)) {
856 if(projectioncode == GCTP_GEO) {
857 lonfield->speciallon =
true;
863 if(HDF4RequestHandler::get_enable_eosgeo_cachefile() ==
true)
864 latfield->speciallon =
true;
876 if(((
int)(lowright[0]/1000)==0) &&((
int)(upleft[0]/1000)==0)
877 && ((
int)(upleft[1]/1000)==0) && ((
int)(lowright[1]/1000)==0)) {
878 if(projectioncode == GCTP_GEO){
879 lonfield->specialformat = 1;
880 latfield->specialformat = 1;
889 if(((
int)(lowright[0])==0) &&((
int)(upleft[0])==0)
890 && ((
int)(upleft[1])==0) && ((
int)(lowright[1])==0)) {
891 if(projectioncode == GCTP_GEO){
892 lonfield->specialformat = 2;
893 latfield->specialformat = 2;
894 gdset->addfvalueattr =
true;
905 if(((
int)(lowright[0])==-1) &&((
int)(upleft[0])==-1)
906 && ((
int)(upleft[1])==-1) && ((
int)(lowright[1])==-1)) {
907 lonfield->specialformat = 3;
908 latfield->specialformat = 3;
909 lonfield->condenseddim =
true;
910 latfield->condenseddim =
true;
915 if(GCTP_SOM == projectioncode) {
916 lonfield->specialformat = 4;
917 latfield->specialformat = 4;
923 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
924 lonfield->condenseddim =
true;
925 latfield->condenseddim =
true;
929 gdset->datafields.push_back(latfield);
930 gdset->datafields.push_back(lonfield);
941 for (
const auto &dim:lonfield->getDimensions()) {
942 if(dim->getName() == DIMXNAME)
945 if(dim->getName() == DIMYNAME)
954void File::handle_onelatlon_grids() throw(Exception) {
957 string DIMXNAME = this->get_geodim_x_name();
958 string DIMYNAME = this->get_geodim_y_name();
959 string LATFIELDNAME = this->get_latfield_name();
960 string LONFIELDNAME = this->get_lonfield_name();
964 string GEOGRIDNAME =
"location";
967 map<string,string>temponelatlondimcvarlist;
970 for (
const auto &grid:this->grids) {
974 grid->setDimxName(DIMXNAME);
975 grid->setDimyName(DIMYNAME);
978 if(grid->getName()==GEOGRIDNAME) {
983 grid->lonfield->fieldtype = 2;
984 grid->latfield->fieldtype = 1;
987 if(grid->lonfield->rank >2 || grid->latfield->rank >2)
988 throw2(
"Either the rank of lat or the lon is greater than 2",grid->getName());
989 if(grid->lonfield->rank !=grid->latfield->rank)
990 throw2(
"The rank of the latitude is not the same as the rank of the longitude",grid->getName());
993 if(grid->lonfield->rank != 1) {
997 grid->lonfield->ydimmajor = grid->getCalculated().isYDimMajor();
998 grid->latfield->ydimmajor = grid->lonfield->ydimmajor;
1003 int32 projectioncode = grid->getProjection().getCode();
1004 if(projectioncode == GCTP_GEO || projectioncode ==GCTP_CEA) {
1005 grid->lonfield->condenseddim =
true;
1006 grid->latfield->condenseddim =
true;
1019 for (
const auto &dim:grid->lonfield->getDimensions()) {
1020 if(dim->getName() == DIMXNAME) {
1022 grid->lonfield->getName());
1024 if(dim->getName() == DIMYNAME) {
1026 grid->latfield->getName());
1032 grid->lonfield->getName());
1034 grid->latfield->getName());
1036 temponelatlondimcvarlist = grid->dimcvarlist;
1044 for (
const auto &grid:this->grids) {
1046 string templatlonname1;
1047 string templatlonname2;
1049 if(grid->getName() != GEOGRIDNAME) {
1051 map<string,string>::iterator tempmapit;
1054 tempmapit = temponelatlondimcvarlist.find(DIMXNAME);
1055 if(tempmapit != temponelatlondimcvarlist.end())
1056 templatlonname1= tempmapit->second;
1058 throw2(
"cannot find the dimension field of XDim", grid->getName());
1063 tempmapit = temponelatlondimcvarlist.find(DIMYNAME);
1064 if(tempmapit != temponelatlondimcvarlist.end())
1065 templatlonname2= tempmapit->second;
1067 throw2(
"cannot find the dimension field of YDim", grid->getName());
1075void File::handle_grid_dim_cvar_maps() throw(Exception) {
1078 string DIMXNAME = this->get_geodim_x_name();
1080 string DIMYNAME = this->get_geodim_y_name();
1082 string LATFIELDNAME = this->get_latfield_name();
1084 string LONFIELDNAME = this->get_lonfield_name();
1089 string GEOGRIDNAME =
"location";
1103 vector <string> tempfieldnamelist;
1104 for (
const auto &grid:this->grids) {
1105 for (
const auto &field:grid->getDataFields())
1114 map<string,string>tempncvarnamelist;
1115 string tempcorrectedlatname, tempcorrectedlonname;
1117 int total_fcounter = 0;
1119 for (
const auto &grid:this->grids) {
1124 for (
const auto &field:grid->getDataFields())
1126 field->newname = tempfieldnamelist[total_fcounter];
1130 if(field->fieldtype!=0) {
1132 tempncvarnamelist.insert(make_pair(field->getName(), field->newname));
1135 if((this->onelatlon)&&((grid->getName())==GEOGRIDNAME)) {
1136 if(field->getName()==LATFIELDNAME)
1137 tempcorrectedlatname = field->newname;
1138 if(field->getName()==LONFIELDNAME)
1139 tempcorrectedlonname = field->newname;
1144 grid->ncvarnamelist = tempncvarnamelist;
1145 tempncvarnamelist.clear();
1151 if(this->onelatlon) {
1152 for(
const auto &grid:this->grids) {
1154 if(grid->getName()!=GEOGRIDNAME){
1162 map<string,string>tempndimnamelist;
1165 vector <string>tempalldimnamelist;
1166 for (
const auto &grid:this->grids) {
1167 for (map<string,string>::const_iterator j =
1168 grid->dimcvarlist.begin(); j!= grid->dimcvarlist.end();++j)
1176 int total_dcounter = 0;
1177 for (
const auto &grid:this->grids) {
1179 for (map<string,string>::const_iterator j =
1180 grid->dimcvarlist.begin(); j!= grid->dimcvarlist.end();++j){
1183 if((DIMXNAME == (*j).first || DIMYNAME == (*j).first) && (
true==(this->onelatlon)))
1190 grid->ndimnamelist = tempndimnamelist;
1191 tempndimnamelist.clear();
1196void File::handle_grid_coards() throw(Exception) {
1199 string DIMXNAME = this->get_geodim_x_name();
1200 string DIMYNAME = this->get_geodim_y_name();
1201 string LATFIELDNAME = this->get_latfield_name();
1202 string LONFIELDNAME = this->get_lonfield_name();
1206 string GEOGRIDNAME =
"location";
1210 vector<Dimension*> correcteddims;
1211 string tempcorrecteddimname;
1214 map<string,string> tempnewxdimnamelist;
1217 map<string,string> tempnewydimnamelist;
1220 Dimension *correcteddim;
1222 for (
const auto &grid:this->grids) {
1223 for (
const auto &field:grid->getDataFields()) {
1228 if(field->getName()==LATFIELDNAME && field->getRank()==2 &&field->condenseddim) {
1230 string templatdimname;
1231 map<string,string>::iterator tempmapit;
1234 tempmapit = grid->ncvarnamelist.find(LATFIELDNAME);
1235 if(tempmapit != grid->ncvarnamelist.end())
1236 templatdimname= tempmapit->second;
1238 throw2(
"cannot find the corrected field of Latitude", grid->getName());
1240 for (
const auto &dim:field->getDimensions()) {
1244 if(dim->getName()==DIMYNAME) {
1245 correcteddim =
new Dimension(templatdimname,dim->getSize());
1246 correcteddims.push_back(correcteddim);
1247 field->setCorrectedDimensions(correcteddims);
1251 field->iscoard =
true;
1252 grid->iscoard =
true;
1253 if (this->onelatlon)
1254 this->iscoard =
true;
1257 correcteddims.clear();
1261 else if(field->getName()==LONFIELDNAME && field->getRank()==2 &&field->condenseddim){
1263 string templondimname;
1264 map<string,string>::iterator tempmapit;
1267 tempmapit = grid->ncvarnamelist.find(LONFIELDNAME);
1268 if(tempmapit != grid->ncvarnamelist.end())
1269 templondimname= tempmapit->second;
1271 throw2(
"cannot find the corrected field of Longitude", grid->getName());
1273 for (
const auto &dim:field->getDimensions()) {
1277 if(dim->getName()==DIMXNAME) {
1278 correcteddim =
new Dimension(templondimname,dim->getSize());
1279 correcteddims.push_back(correcteddim);
1280 field->setCorrectedDimensions(correcteddims);
1285 field->iscoard =
true;
1286 grid->iscoard =
true;
1288 this->iscoard =
true;
1289 correcteddims.clear();
1293 else if((field->getRank()==1) &&(field->getName()==LONFIELDNAME) ) {
1295 string templondimname;
1296 map<string,string>::iterator tempmapit;
1299 tempmapit = grid->ncvarnamelist.find(LONFIELDNAME);
1300 if(tempmapit != grid->ncvarnamelist.end())
1301 templondimname= tempmapit->second;
1303 throw2(
"cannot find the corrected field of Longitude", grid->getName());
1305 correcteddim =
new Dimension(templondimname,(field->getDimensions())[0]->getSize());
1306 correcteddims.push_back(correcteddim);
1307 field->setCorrectedDimensions(correcteddims);
1308 field->iscoard =
true;
1309 grid->iscoard =
true;
1311 this->iscoard =
true;
1312 correcteddims.clear();
1314 if(((field->getDimensions())[0]->getName()!=DIMXNAME)
1315 &&(((field->getDimensions())[0]->getName())!=DIMYNAME)){
1316 throw3(
"the dimension name of longitude should not be ",
1317 (field->getDimensions())[0]->getName(),grid->getName());
1319 if(((field->getDimensions())[0]->getName())==DIMXNAME) {
1328 else if((field->getRank()==1) &&(field->getName()==LATFIELDNAME) ) {
1330 string templatdimname;
1331 map<string,string>::iterator tempmapit;
1334 tempmapit = grid->ncvarnamelist.find(LATFIELDNAME);
1335 if(tempmapit != grid->ncvarnamelist.end())
1336 templatdimname= tempmapit->second;
1338 throw2(
"cannot find the corrected field of Latitude", grid->getName());
1340 correcteddim =
new Dimension(templatdimname,(field->getDimensions())[0]->getSize());
1341 correcteddims.push_back(correcteddim);
1342 field->setCorrectedDimensions(correcteddims);
1344 field->iscoard =
true;
1345 grid->iscoard =
true;
1347 this->iscoard =
true;
1348 correcteddims.clear();
1350 if((((field->getDimensions())[0]->getName())!=DIMXNAME)
1351 &&(((field->getDimensions())[0]->getName())!=DIMYNAME))
1352 throw3(
"the dimension name of latitude should not be ",
1353 (field->getDimensions())[0]->getName(),grid->getName());
1354 if(((field->getDimensions())[0]->getName())==DIMXNAME){
1366 if(
true == this->onelatlon){
1369 if(
true == this->iscoard){
1372 string tempcorrectedxdimname;
1373 string tempcorrectedydimname;
1375 if((
int)(tempnewxdimnamelist.size())!= 1)
1376 throw1(
"the corrected dimension name should have only one pair");
1377 if((
int)(tempnewydimnamelist.size())!= 1)
1378 throw1(
"the corrected dimension name should have only one pair");
1380 map<string,string>::iterator tempdimmapit = tempnewxdimnamelist.begin();
1381 tempcorrectedxdimname = tempdimmapit->second;
1382 tempdimmapit = tempnewydimnamelist.begin();
1383 tempcorrectedydimname = tempdimmapit->second;
1385 for (
const auto &grid:this->grids) {
1388 map<string,string>::iterator tempmapit;
1389 tempmapit = grid->ndimnamelist.find(DIMXNAME);
1390 if(tempmapit != grid->ndimnamelist.end()) {
1394 throw2(
"cannot find the corrected dimension name", grid->getName());
1395 tempmapit = grid->ndimnamelist.find(DIMYNAME);
1396 if(tempmapit != grid->ndimnamelist.end()) {
1400 throw2(
"cannot find the corrected dimension name", grid->getName());
1405 for (
const auto &grid:this->grids) {
1409 string tempcorrectedxdimname;
1410 string tempcorrectedydimname;
1413 map<string,string>::iterator tempdimmapit;
1414 map<string,string>::iterator tempmapit;
1415 tempdimmapit = tempnewxdimnamelist.find(grid->getName());
1416 if(tempdimmapit != tempnewxdimnamelist.end())
1417 tempcorrectedxdimname = tempdimmapit->second;
1419 throw2(
"cannot find the corrected COARD XDim dimension name", grid->getName());
1420 tempmapit = grid->ndimnamelist.find(DIMXNAME);
1421 if(tempmapit != grid->ndimnamelist.end()) {
1425 throw2(
"cannot find the corrected dimension name", grid->getName());
1427 tempdimmapit = tempnewydimnamelist.find(grid->getName());
1428 if(tempdimmapit != tempnewydimnamelist.end())
1429 tempcorrectedydimname = tempdimmapit->second;
1431 throw2(
"cannot find the corrected COARD YDim dimension name", grid->getName());
1433 tempmapit = grid->ndimnamelist.find(DIMYNAME);
1434 if(tempmapit != grid->ndimnamelist.end()) {
1438 throw2(
"cannot find the corrected dimension name", grid->getName());
1446 for (
const auto &grid:this->grids){
1447 for (map<string,string>::const_iterator j =
1448 grid->dimcvarlist.begin(); j!= grid->dimcvarlist.end();++j){
1452 if((this->iscoard||grid->iscoard) && (*j).first !=DIMXNAME && (*j).first !=DIMYNAME) {
1453 string tempnewdimname;
1454 map<string,string>::iterator tempmapit;
1457 tempmapit = grid->ncvarnamelist.find((*j).second);
1458 if(tempmapit != grid->ncvarnamelist.end())
1459 tempnewdimname= tempmapit->second;
1461 throw3(
"cannot find the corrected field of ", (*j).second,grid->getName());
1464 tempmapit =grid->ndimnamelist.find((*j).first);
1465 if(tempmapit != grid->ndimnamelist.end())
1468 throw3(
"cannot find the corrected dimension name of ", (*j).first,grid->getName());
1476void File::update_grid_field_corrected_dims() throw(Exception) {
1479 vector<Dimension*> correcteddims;
1480 string tempcorrecteddimname;
1482 Dimension *correcteddim;
1484 for (
const auto &grid:this->grids) {
1486 for (
const auto &field:grid->getDataFields()) {
1489 if (field->iscoard ==
false) {
1492 for (
const auto &dim:field->getDimensions()){
1494 map<string,string>::iterator tempmapit;
1497 tempmapit = grid->ndimnamelist.find(dim->getName());
1498 if(tempmapit != grid->ndimnamelist.end())
1499 tempcorrecteddimname= tempmapit->second;
1501 throw4(
"cannot find the corrected dimension name", grid->getName(),field->getName(),dim->getName());
1502 correcteddim =
new Dimension(tempcorrecteddimname,dim->getSize());
1503 correcteddims.push_back(correcteddim);
1505 field->setCorrectedDimensions(correcteddims);
1506 correcteddims.clear();
1513void File::handle_grid_cf_attrs() throw(Exception) {
1520 for (
const auto &grid:this->grids) {
1521 for (
const auto &field:grid->getDataFields()) {
1524 if(field->fieldtype == 0) {
1525 string tempcoordinates=
"";
1526 string tempfieldname=
"";
1527 string tempcorrectedfieldname=
"";
1529 for (
const auto &dim:field->getDimensions()) {
1532 map<string,string>::iterator tempmapit;
1533 map<string,string>::iterator tempmapit2;
1536 tempmapit = (grid->dimcvarlist).find(dim->getName());
1537 if(tempmapit != (grid->dimcvarlist).end())
1538 tempfieldname = tempmapit->second;
1540 throw4(
"cannot find the dimension field name",
1541 grid->getName(),field->getName(),dim->getName());
1544 tempmapit2 = (grid->ncvarnamelist).find(tempfieldname);
1545 if(tempmapit2 != (grid->ncvarnamelist).end())
1546 tempcorrectedfieldname = tempmapit2->second;
1548 throw4(
"cannot find the corrected dimension field name",
1549 grid->getName(),field->getName(),dim->getName());
1552 tempcoordinates= tempcorrectedfieldname;
1554 tempcoordinates = tempcoordinates +
" "+tempcorrectedfieldname;
1557 field->setCoordinates(tempcoordinates);
1561 if(field->fieldtype == 1) {
1562 string tempunits =
"degrees_north";
1563 field->setUnits(tempunits);
1565 if(field->fieldtype == 2) {
1566 string tempunits =
"degrees_east";
1567 field->setUnits(tempunits);
1574 if(field->fieldtype == 4) {
1575 string tempunits =
"level";
1576 field->setUnits(tempunits);
1580 if(field->fieldtype == 5) {
1581 string tempunits =
"days since 1900-01-01 00:00:00";
1582 field->setUnits(tempunits);
1590 if (
true == grid->addfvalueattr) {
1591 if(((field->getFillValue()).empty()) && (field->getType()==DFNT_FLOAT32 )) {
1592 float tempfillvalue = FLT_MAX;
1593 field->addFillValue(tempfillvalue);
1594 field->setAddedFillValue(
true);
1602void File::handle_grid_SOM_projection() throw(Exception) {
1611 for (
const auto &grid:this->grids) {
1612 if (GCTP_SOM == grid->getProjection().getCode()) {
1618 for (
const auto &dim:grid->getDimensions()) {
1621 if(NBLOCK == dim->getSize()) {
1625 if (dim->getName().compare(0,3,
"SOM") == 0) {
1626 som_dimname = dim->getName();
1632 if(
""== som_dimname)
1633 throw4(
"Wrong number of block: The number of block of MISR SOM Grid ",
1634 grid->getName(),
" is not ",NBLOCK);
1636 map<string,string>::iterator tempmapit;
1639 string cor_som_dimname;
1640 tempmapit = grid->ndimnamelist.find(som_dimname);
1641 if(tempmapit != grid->ndimnamelist.end())
1642 cor_som_dimname = tempmapit->second;
1644 throw2(
"cannot find the corrected dimension name for ", som_dimname);
1647 string cor_som_cvname;
1651 for (vector<Field *>::iterator j = grid->datafields.begin();
1652 j != grid->datafields.end(); ) {
1656 if (1 == (*j)->fieldtype || 2 == (*j)->fieldtype) {
1658 auto newdim =
new Dimension(som_dimname,NBLOCK);
1659 auto newcor_dim =
new Dimension(cor_som_dimname,NBLOCK);
1660 vector<Dimension *>::iterator it_d;
1662 it_d = (*j)->dims.begin();
1663 (*j)->dims.insert(it_d,newdim);
1665 it_d = (*j)->correcteddims.begin();
1666 (*j)->correcteddims.insert(it_d,newcor_dim);
1675 if ( 4 == (*j)->fieldtype) {
1676 cor_som_cvname = (*j)->newname;
1678 j = grid->datafields.erase(j);
1692 for (
const auto &field:grid->getDataFields()) {
1694 if ( 0 == field->fieldtype) {
1696 string temp_coordinates = field->coordinates;
1699 found = temp_coordinates.find(cor_som_cvname);
1703 if (temp_coordinates.size() >cor_som_cvname.size())
1704 temp_coordinates.erase(found,cor_som_cvname.size()+1);
1706 temp_coordinates.erase(found,cor_som_cvname.size());
1708 else if (found != string::npos)
1709 temp_coordinates.erase(found-1,cor_som_cvname.size()+1);
1711 throw4(
"cannot find the coordinate variable ",cor_som_cvname,
1712 "from ",temp_coordinates);
1714 field->setCoordinates(temp_coordinates);
1724void File::check_swath_dimmap(
int numswath)
throw(Exception) {
1726 if(HDF4RequestHandler::get_disable_swath_dim_map() ==
true)
1731 int temp_num_map = 0;
1732 bool odd_num_map =
false;
1733 for (
const auto &swath:this->swaths) {
1734 temp_num_map = swath->get_num_map();
1735 tempnumdm += temp_num_map;
1736 if(temp_num_map%2!=0) {
1743 if(tempnumdm != 0 && odd_num_map ==
false)
1744 handle_swath_dimmap =
true;
1753 bool fakedimmap =
false;
1757 if((this->swaths[0]->getName()).find(
"atml2")!=string::npos){
1763 for (
const auto &gfield:this->swaths[0]->getGeoFields()) {
1764 if(gfield->getName() ==
"Latitude" || gfield->getName() ==
"Longitude") {
1765 if (gfield->getType() == DFNT_UINT16 ||
1766 gfield->getType() == DFNT_INT16)
1767 gfield->type = DFNT_FLOAT32;
1776 for (
const auto &dfield:this->swaths[0]->getDataFields()) {
1791 if((dfield->getName()).find(
"Latitude") != string::npos){
1793 if (dfield->getType() == DFNT_UINT16 ||
1794 dfield->getType() == DFNT_INT16)
1795 dfield->type = DFNT_FLOAT32;
1797 dfield->fieldtype = 1;
1800 if(dfield->getRank() != 2)
1801 throw2(
"The lat/lon rank must be 2 for Java clients to work",
1804 ((dfield->getDimensions())[0])->getName(),dfield->getName());
1808 if((dfield->getName()).find(
"Longitude")!= string::npos) {
1810 if(dfield->getType() == DFNT_UINT16 ||
1811 dfield->getType() == DFNT_INT16)
1812 dfield->type = DFNT_FLOAT32;
1814 dfield->fieldtype = 2;
1815 if(dfield->getRank() != 2)
1816 throw2(
"The lat/lon rank must be 2 for Java clients to work",
1819 ((dfield->getDimensions())[1])->getName(), dfield->getName());
1831 if(
true == fakedimmap)
1832 handle_swath_dimmap =
false;
1839void File::check_swath_dimmap_bk_compat(
int numswath){
1841 if(
true == handle_swath_dimmap) {
1843 if(numswath == 1 && (((this->swaths)[0])->name==
"MODIS_SWATH_Type_L1B"))
1844 backward_handle_swath_dimmap =
true;
1850 bool all_2_dimmaps_no_geodim =
true;
1851 for (
const auto &swath:this->swaths) {
1852 if (swath->get_num_map() !=2 || swath->GeoDim_in_vars ==
true) {
1853 all_2_dimmaps_no_geodim =
false;
1857 if (
true == all_2_dimmaps_no_geodim)
1858 backward_handle_swath_dimmap =
true;
1865void File::create_swath_latlon_dim_cvar_map() throw(Exception){
1867 vector<Field*> ori_lats;
1868 vector<Field*> ori_lons;
1869 if(handle_swath_dimmap ==
true && backward_handle_swath_dimmap ==
false) {
1874 multi_dimmap =
true;
1876 for (
const auto &swath:this->swaths) {
1878 bool has_cf_lat =
false;
1879 bool has_cf_lon =
false;
1881 for (
const auto &gfield:swath->getGeoFields()) {
1886 if(gfield->getName()==
"Latitude" && gfield->getRank() == 2){
1888 ori_lats.push_back(gfield);
1890 else if(gfield->getName()==
"Longitude" && gfield->getRank() == 2){
1892 ori_lons.push_back(gfield);
1894 if(has_cf_lat ==
true && has_cf_lon ==
true)
1897 if(has_cf_lat ==
false || has_cf_lon ==
false) {
1898 multi_dimmap =
false;
1907 if (
true == multi_dimmap) {
1910 for (
const auto &swath:this->swaths) {
1911 create_swath_latlon_dim_cvar_map_for_dimmap(swath,ori_lats[ll_count],ori_lons[ll_count]);
1932 bool lat_in_geofields =
false;
1933 bool lon_in_geofields =
false;
1935 for (
const auto &swath:this->swaths) {
1937 int tempgeocount = 0;
1938 for (
const auto &gfield:swath->getGeoFields()) {
1942 if(gfield->getName()==
"Latitude" ){
1943 if(gfield->getRank() > 2)
1944 throw2(
"Currently the lat/lon rank must be 1 or 2 for Java clients to work",
1947 lat_in_geofields =
true;
1958 if(handle_swath_dimmap ==
true) {
1961 if(
true == backward_handle_swath_dimmap) {
1964 for (
const auto &dmap:swath->getDimensionMaps()) {
1968 if((gfield->getDimensions()[0])->getName() == dmap->getGeoDimension()) {
1976 gfield->fieldtype = 1;
1980 if(gfield->getName()==
"Longitude"){
1981 if(gfield->getRank() > 2)
1982 throw2(
"Currently the lat/lon rank must be 1 or 2 for Java clients to work",
1987 lon_in_geofields =
true;
1988 if(gfield->getRank() == 1) {
1998 ((gfield->getDimensions())[1])->getName(),
"Longitude");
1999 if(handle_swath_dimmap ==
true) {
2000 if(
true == backward_handle_swath_dimmap) {
2003 for (
const auto &dmap:swath->getDimensionMaps()) {
2008 if((gfield->getDimensions()[1])->getName() ==
2009 dmap->getGeoDimension()) {
2011 dmap->getDataDimension(),
"Longitude");
2017 gfield->fieldtype = 2;
2020 if(tempgeocount == 2)
2027 if (lat_in_geofields!=lon_in_geofields)
2028 throw1(
"Latitude and longitude must be both under Geolocation fields or Data fields");
2031 if (!lat_in_geofields && !lon_in_geofields) {
2033 bool lat_in_datafields =
false;
2034 bool lon_in_datafields =
false;
2036 for (
const auto &swath:this->swaths) {
2038 int tempgeocount = 0;
2039 for (
const auto &dfield:swath->getDataFields()) {
2045 if (dfield->getName()==
"Latitude" ){
2046 if (dfield->getRank() > 2) {
2047 throw2(
"Currently the lat/lon rank must be 1 or 2 for Java clients to work",
2050 lat_in_datafields =
true;
2059 ((dfield->getDimensions())[0])->getName(),
"Latitude");
2061 if (handle_swath_dimmap ==
true) {
2062 if (
true == backward_handle_swath_dimmap) {
2064 for (
const auto &dmap:swath->getDimensionMaps()) {
2067 if((dfield->getDimensions()[0])->getName() == dmap->getGeoDimension()) {
2074 dfield->fieldtype = 1;
2078 if(dfield->getName()==
"Longitude"){
2080 if(dfield->getRank() > 2) {
2081 throw2(
"Currently the lat/lon rank must be 1 or 2 for Java clients to work",
2087 lon_in_datafields =
true;
2088 if(dfield->getRank() == 1) {
2098 ((dfield->getDimensions())[1])->getName(),
"Longitude");
2099 if(handle_swath_dimmap ==
true) {
2100 if(
true == backward_handle_swath_dimmap) {
2102 for (
const auto &dmap:swath->getDimensionMaps()) {
2105 if((dfield->getDimensions()[1])->getName() == dmap->getGeoDimension()) {
2107 dmap->getDataDimension(),
"Longitude");
2113 dfield->fieldtype = 2;
2116 if(tempgeocount == 2)
2123 if (lat_in_datafields!=lon_in_datafields)
2124 throw1(
"Latitude and longitude must be both under Geolocation fields or Data fields");
2132void File:: create_swath_nonll_dim_cvar_map() throw(Exception)
2135 for (
const auto &swath:this->swaths) {
2152 pair<set<string>::iterator,
bool> tempdimret;
2153 for (map<string,string>::const_iterator j = swath->dimcvarlist.begin();
2154 j!= swath->dimcvarlist.end();++j){
2155 tempdimret = swath->nonmisscvdimlist.insert((*j).first);
2162 for (
const auto &gfield:swath->getGeoFields()) {
2164 if(gfield->getRank()==1) {
2165 if(swath->nonmisscvdimlist.find(((gfield->getDimensions())[0])->getName()) == swath->nonmisscvdimlist.end()){
2166 tempdimret = swath->nonmisscvdimlist.insert(((gfield->getDimensions())[0])->getName());
2167 if(gfield->getName() ==
"Time")
2168 gfield->fieldtype = 5;
2178 if((((gfield->getDimensions())[0])->getName())==gfield->getName()){
2179 gfield->oriname = gfield->getName();
2182 gfield->specialcoard =
true;
2185 HDFCFUtil::insert_map(swath->dimcvarlist, ((gfield->getDimensions())[0])->getName(), gfield->getName());
2186 gfield->fieldtype = 3;
2196 for (
const auto &dfield:swath->getDataFields()) {
2198 if(dfield->getRank()==1) {
2199 if(swath->nonmisscvdimlist.find(((dfield->getDimensions())[0])->getName()) == swath->nonmisscvdimlist.end()){
2200 tempdimret = swath->nonmisscvdimlist.insert(((dfield->getDimensions())[0])->getName());
2201 if(dfield->getName() ==
"Time")
2202 dfield->fieldtype = 5;
2209 if((((dfield->getDimensions())[0])->getName())==dfield->getName()){
2210 dfield->oriname = dfield->getName();
2212 dfield->specialcoard =
true;
2215 HDFCFUtil::insert_map(swath->dimcvarlist, ((dfield->getDimensions())[0])->getName(), dfield->getName());
2216 dfield->fieldtype = 3;
2226 bool missingfield_unlim_flag =
false;
2227 Field *missingfield_unlim =
nullptr;
2229 for (
const auto &sdim:swath->getDimensions()) {
2231 if((swath->nonmisscvdimlist.find(sdim->getName())) == swath->nonmisscvdimlist.end()){
2234 auto missingfield =
new Field();
2247 if(
true == multi_dimmap && (this->swaths.size() != 1)) {
2248 missingfield->name = sdim->getName()+
"_"+swath->name;
2249 dim =
new Dimension(missingfield->name,sdim->getSize());
2252 missingfield->name = sdim->getName();
2253 dim =
new Dimension(sdim->getName(),sdim->getSize());
2255 missingfield->rank = 1;
2256 missingfield->type = DFNT_INT32;
2259 missingfield->dims.push_back(dim);
2266 if(0 == sdim->getSize()) {
2267 missingfield_unlim_flag =
true;
2268 missingfield_unlim = missingfield;
2272 missingfield->fieldtype = 4;
2274 swath->geofields.push_back(missingfield);
2276 (missingfield->getDimensions())[0]->getName(), missingfield->name);
2286 bool temp_missingfield_unlim_flag = missingfield_unlim_flag;
2287 if(
true == temp_missingfield_unlim_flag) {
2288 for (
const auto &dfield:swath->getDataFields()) {
2290 for (
const auto &fdim:dfield->getDimensions()) {
2292 if(fdim->getName() == (missingfield_unlim->getDimensions())[0]->getName()) {
2293 if(fdim->getSize()!= 0) {
2294 Dimension *dim = missingfield_unlim->getDimensions()[0];
2296 dim->dimsize = fdim->getSize();
2297 missingfield_unlim_flag =
false;
2303 if(
false == missingfield_unlim_flag)
2308 swath->nonmisscvdimlist.clear();
2315void File::handle_swath_dim_cvar_maps() throw(Exception) {
2318 vector <string> tempfieldnamelist;
2319 for (
const auto &swath:this->swaths) {
2322 for (
const auto &gfield:swath->getGeoFields()) {
2323 if(gfield->fieldtype == 0 && (this->swaths.size() !=1) &&
2324 (
true == handle_swath_dimmap) &&
2325 (backward_handle_swath_dimmap ==
false)){
2326 string new_field_name = gfield->name+
"_"+swath->name;
2333 for (
const auto &dfield:swath->getDataFields()) {
2334 if(dfield->fieldtype == 0 && (this->swaths.size() !=1) &&
2335 true == multi_dimmap){
2339 string new_field_name = dfield->name+
"_"+swath->name;
2349 int total_fcounter = 0;
2354 map<string,string>tempncvarnamelist;
2355 for (
const auto &swath:this->swaths) {
2358 for (
const auto &gfield:swath->getGeoFields())
2361 gfield->newname = tempfieldnamelist[total_fcounter];
2365 if(gfield->fieldtype!=0) {
2370 for (
const auto &dfield:swath->getDataFields())
2372 dfield->newname = tempfieldnamelist[total_fcounter];
2376 if(dfield->fieldtype!=0) {
2384 vector <string>tempalldimnamelist;
2386 for (
const auto &swath:this->swaths) {
2387 for (map<string,string>::const_iterator j =
2388 swath->dimcvarlist.begin(); j!= swath->dimcvarlist.end();++j)
2395 int total_dcounter = 0;
2397 for (
const auto &swath:this->swaths) {
2398 for (map<string,string>::const_iterator j =
2399 swath->dimcvarlist.begin(); j!= swath->dimcvarlist.end();++j){
2406 vector<Dimension*> correcteddims;
2407 string tempcorrecteddimname;
2408 Dimension *correcteddim;
2410 for (
const auto &swath:this->swaths) {
2413 for (
const auto &gfield:swath->getGeoFields()) {
2415 for (
const auto &gdim:gfield->getDimensions()) {
2417 map<string,string>::iterator tempmapit;
2420 if(handle_swath_dimmap ==
false || multi_dimmap ==
true) {
2423 tempmapit = swath->ndimnamelist.find(gdim->getName());
2424 if(tempmapit != swath->ndimnamelist.end())
2425 tempcorrecteddimname= tempmapit->second;
2427 throw4(
"cannot find the corrected dimension name",
2428 swath->getName(),gfield->getName(),gdim->getName());
2430 correcteddim =
new Dimension(tempcorrecteddimname,gdim->getSize());
2434 bool isdimmapname =
false;
2437 for (
const auto &sdmap:swath->getDimensionMaps()) {
2441 if(gdim->getName() == sdmap->getGeoDimension()) {
2443 isdimmapname =
true;
2444 gfield->dmap =
true;
2445 string temprepdimname = sdmap->getDataDimension();
2448 tempmapit = swath->ndimnamelist.find(temprepdimname);
2449 if(tempmapit != swath->ndimnamelist.end())
2450 tempcorrecteddimname= tempmapit->second;
2452 throw4(
"cannot find the corrected dimension name", swath->getName(),
2453 gfield->getName(),gdim->getName());
2457 bool ddimsflag =
false;
2458 for (
const auto &sdim:swath->getDimensions()) {
2459 if(sdim->getName() == temprepdimname) {
2461 correcteddim =
new Dimension(tempcorrecteddimname,sdim->getSize());
2467 throw4(
"cannot find the corrected dimension size", swath->getName(),
2468 gfield->getName(),gdim->getName());
2472 if(
false == isdimmapname) {
2474 tempmapit = swath->ndimnamelist.find(gdim->getName());
2475 if(tempmapit != swath->ndimnamelist.end())
2476 tempcorrecteddimname= tempmapit->second;
2478 throw4(
"cannot find the corrected dimension name",
2479 swath->getName(),gfield->getName(),gdim->getName());
2481 correcteddim =
new Dimension(tempcorrecteddimname,gdim->getSize());
2486 correcteddims.push_back(correcteddim);
2488 gfield->setCorrectedDimensions(correcteddims);
2489 correcteddims.clear();
2493 for (
const auto &dfield:swath->getDataFields()) {
2495 for (
const auto &fdim:dfield->getDimensions()) {
2497 if((handle_swath_dimmap ==
false) || multi_dimmap ==
true) {
2499 map<string,string>::iterator tempmapit;
2501 tempmapit = swath->ndimnamelist.find(fdim->getName());
2502 if(tempmapit != swath->ndimnamelist.end())
2503 tempcorrecteddimname= tempmapit->second;
2505 throw4(
"cannot find the corrected dimension name", swath->getName(),
2506 dfield->getName(),fdim->getName());
2508 correcteddim =
new Dimension(tempcorrecteddimname,fdim->getSize());
2511 map<string,string>::iterator tempmapit;
2512 bool isdimmapname =
false;
2515 for (
const auto &smap:swath->getDimensionMaps()) {
2518 if(fdim->getName() == smap->getGeoDimension()) {
2519 isdimmapname =
true;
2520 dfield->dmap =
true;
2521 string temprepdimname = smap->getDataDimension();
2524 tempmapit = swath->ndimnamelist.find(temprepdimname);
2525 if(tempmapit != swath->ndimnamelist.end())
2526 tempcorrecteddimname= tempmapit->second;
2528 throw4(
"cannot find the corrected dimension name",
2529 swath->getName(),dfield->getName(),fdim->getName());
2533 bool ddimsflag =
false;
2534 for(
const auto &sdim:swath->getDimensions()) {
2537 if(sdim->getName() == temprepdimname) {
2538 correcteddim =
new Dimension(tempcorrecteddimname,sdim->getSize());
2544 throw4(
"cannot find the corrected dimension size",
2545 swath->getName(),dfield->getName(),fdim->getName());
2553 tempmapit = swath->ndimnamelist.find(fdim->getName());
2554 if(tempmapit != swath->ndimnamelist.end())
2555 tempcorrecteddimname= tempmapit->second;
2557 throw4(
"cannot find the corrected dimension name",
2558 swath->getName(),dfield->getName(),fdim->getName());
2560 correcteddim =
new Dimension(tempcorrecteddimname,fdim->getSize());
2564 correcteddims.push_back(correcteddim);
2566 dfield->setCorrectedDimensions(correcteddims);
2567 correcteddims.clear();
2574void File::handle_swath_cf_attrs() throw(Exception) {
2584 for (
const auto &swath:this->swaths) {
2587 for (
const auto &gfield:swath->getGeoFields()) {
2590 if(gfield->fieldtype == 0) {
2591 string tempcoordinates=
"";
2592 string tempfieldname=
"";
2593 string tempcorrectedfieldname=
"";
2595 bool has_ll_coord =
false;
2596 if(swath->get_num_map() == 0)
2597 has_ll_coord =
true;
2598 else if(handle_swath_dimmap ==
true) {
2599 if(backward_handle_swath_dimmap ==
true || multi_dimmap ==
true)
2600 has_ll_coord =
true;
2602 for (
const auto &dim:gfield->getDimensions()) {
2605 map<string,string>::iterator tempmapit;
2606 map<string,string>::iterator tempmapit2;
2609 tempmapit = (swath->dimcvarlist).find(dim->getName());
2610 if(tempmapit != (swath->dimcvarlist).end())
2611 tempfieldname = tempmapit->second;
2613 throw4(
"cannot find the dimension field name",swath->getName(),
2614 gfield->getName(),dim->getName());
2617 tempmapit2 = (swath->ncvarnamelist).find(tempfieldname);
2618 if(tempmapit2 != (swath->ncvarnamelist).end())
2619 tempcorrectedfieldname = tempmapit2->second;
2621 throw4(
"cannot find the corrected dimension field name",
2622 swath->getName(),gfield->getName(),dim->getName());
2624 if(
false == has_ll_coord)
2625 has_ll_coord= check_ll_in_coords(tempcorrectedfieldname);
2628 tempcoordinates= tempcorrectedfieldname;
2630 tempcoordinates = tempcoordinates +
" "+tempcorrectedfieldname;
2633 if(
true == has_ll_coord)
2634 gfield->setCoordinates(tempcoordinates);
2639 if(gfield->fieldtype == 1) {
2640 string tempunits =
"degrees_north";
2641 gfield->setUnits(tempunits);
2645 if(gfield->fieldtype == 2) {
2646 string tempunits =
"degrees_east";
2647 gfield->setUnits(tempunits);
2654 if(gfield->fieldtype == 4) {
2655 string tempunits =
"level";
2656 gfield->setUnits(tempunits);
2661 if(gfield->fieldtype == 5) {
2662 string tempunits =
"days since 1900-01-01 00:00:00";
2663 gfield->setUnits(tempunits);
2669 if(((gfield->getFillValue()).empty()) &&
2670 (gfield->getType()==DFNT_FLOAT32 || gfield->getType()==DFNT_FLOAT64)) {
2671 float tempfillvalue = -9999.0;
2672 gfield->addFillValue(tempfillvalue);
2673 gfield->setAddedFillValue(
true);
2678 for (
const auto &dfield:swath->getDataFields()) {
2681 if(dfield->fieldtype == 0) {
2682 string tempcoordinates=
"";
2683 string tempfieldname=
"";
2684 string tempcorrectedfieldname=
"";
2686 bool has_ll_coord =
false;
2687 if(swath->get_num_map() == 0)
2688 has_ll_coord =
true;
2689 else if(handle_swath_dimmap ==
true) {
2690 if(backward_handle_swath_dimmap ==
true || multi_dimmap ==
true)
2691 has_ll_coord =
true;
2693 for (
const auto &dim:dfield->getDimensions()) {
2696 map<string,string>::iterator tempmapit;
2697 map<string,string>::iterator tempmapit2;
2700 tempmapit = (swath->dimcvarlist).find(dim->getName());
2701 if(tempmapit != (swath->dimcvarlist).end())
2702 tempfieldname = tempmapit->second;
2704 throw4(
"cannot find the dimension field name",swath->getName(),
2705 dfield->getName(),dim->getName());
2708 tempmapit2 = (swath->ncvarnamelist).find(tempfieldname);
2709 if(tempmapit2 != (swath->ncvarnamelist).end())
2710 tempcorrectedfieldname = tempmapit2->second;
2712 throw4(
"cannot find the corrected dimension field name",
2713 swath->getName(),dfield->getName(),dim->getName());
2715 if(
false == has_ll_coord)
2716 has_ll_coord= check_ll_in_coords(tempcorrectedfieldname);
2719 tempcoordinates= tempcorrectedfieldname;
2721 tempcoordinates = tempcoordinates +
" "+tempcorrectedfieldname;
2724 if(
true == has_ll_coord)
2725 dfield->setCoordinates(tempcoordinates);
2728 if((dfield->fieldtype == 3)||(dfield->fieldtype == 4)) {
2729 string tempunits =
"level";
2730 dfield->setUnits(tempunits);
2735 if(dfield->fieldtype == 5) {
2736 string tempunits =
"days since 1900-01-01 00:00:00";
2737 dfield->setUnits(tempunits);
2744 if(((dfield->getFillValue()).empty()) &&
2745 (dfield->getType()==DFNT_FLOAT32 || dfield->getType()==DFNT_FLOAT64)) {
2746 float tempfillvalue = -9999.0;
2747 dfield->addFillValue(tempfillvalue);
2748 dfield->setAddedFillValue(
true);
2755bool File::find_dim_in_dims(
const std::vector<Dimension*>&dims,
const std::string &dim_name)
const {
2757 bool ret_value =
false;
2758 for (
const auto &dim:dims) {
2759 if (dim->name == dim_name) {
2768void File::check_dm_geo_dims_in_vars() {
2770 if(handle_swath_dimmap ==
false)
2773 for (
const auto &swath:this->swaths) {
2776 if(swath->get_num_map() > 0) {
2778 for (
const auto &dfield:swath->getDataFields()) {
2782 if (dfield->rank >=2) {
2783 for (
const auto &dim:dfield->getDimensions()) {
2787 bool not_match_geo_dim =
true;
2788 for (
const auto &sdmap:swath->getDimensionMaps()) {
2790 if((dim->getName() == sdmap->getGeoDimension()) && not_match_geo_dim){
2792 not_match_geo_dim =
false;
2798 if(match_dims == 2) {
2799 swath->GeoDim_in_vars =
true;
2804 if(swath->GeoDim_in_vars ==
false) {
2806 for (
const auto &gfield:swath->getGeoFields()) {
2810 if(gfield->rank >=2 && (gfield->name !=
"Latitude" && gfield->name !=
"Longitude")) {
2811 for (
const auto &dim:gfield->getDimensions()) {
2815 bool not_match_geo_dim =
true;
2817 for (
const auto &sdmap:swath->getDimensionMaps()) {
2818 if((dim->getName() == sdmap->getGeoDimension()) && not_match_geo_dim){
2820 not_match_geo_dim =
false;
2826 if(match_dims == 2){
2827 swath->GeoDim_in_vars =
true;
2840bool SwathDataset::obtain_dmap_offset_inc(
const string& ori_dimname,
const string & mapped_dimname,
int &offset,
int&inc) {
2841 bool ret_value =
false;
2842 for (
const auto &sdmap:this->dimmaps) {
2843 if(sdmap->geodim==ori_dimname && sdmap->datadim == mapped_dimname){
2844 offset = sdmap->offset;
2845 inc = sdmap->increment;
2857void File::create_geo_varnames_list(vector<string> & geo_varnames,
const string & swathname,
2858 const string & fieldname,
int extra_ll_pairs,
bool oneswath) {
2860 if(
true == oneswath)
2861 geo_varnames.push_back(fieldname);
2863 string nfieldname = fieldname+
"_"+swathname;
2864 geo_varnames.push_back(nfieldname);
2866 for (
int i = 0; i <extra_ll_pairs;i++) {
2870 if(
true == oneswath)
2871 nfieldname = fieldname+
"_"+si.str();
2873 nfieldname = fieldname+
"_"+swathname+
"_"+si.str();
2874 geo_varnames.push_back(nfieldname);
2878cerr<<
"ll_pairs is "<<extra_ll_pairs <<endl;
2879for(
int i =0;i<geo_varnames.size();i++)
2880 cerr<<
"geo_varnames["<<i<<
"]= " <<geo_varnames[i] <<endl;
2886void File::create_geo_dim_var_maps(SwathDataset*sd, Field*fd,
const vector<string>& lat_names,
2887 const vector<string>& lon_names,vector<Dimension*>& geo_var_dim1,
2888 vector<Dimension*>& geo_var_dim2) {
2889 string field_lat_dim1_name =(fd->dims)[0]->name;
2890 string field_lat_dim2_name =(fd->dims)[1]->name;
2893 if(sd->GeoDim_in_vars ==
true) {
2894 if((this->swaths).size() >1) {
2895 (fd->dims)[0]->name = field_lat_dim1_name+
"_"+sd->name;
2896 (fd->dims)[1]->name = field_lat_dim2_name+
"_"+sd->name;
2898 geo_var_dim1.push_back((fd->dims)[0]);
2899 geo_var_dim2.push_back((fd->dims)[1]);
2909 short dim1_map_count = 0;
2910 short dim2_map_count = 0;
2911 for (
const auto &sdmap:sd->getDimensionMaps()) {
2912 if(sdmap->getGeoDimension()==field_lat_dim1_name){
2913 string data_dim1_name = sdmap->getDataDimension();
2914 int dim1_size = sd->obtain_dimsize_with_dimname(data_dim1_name);
2915 if((this->swaths).size() > 1)
2916 data_dim1_name = data_dim1_name+
"_"+sd->name;
2918 if(sd->GeoDim_in_vars ==
false && dim1_map_count == 0) {
2919 (fd->dims)[0]->name = data_dim1_name;
2920 (fd->dims)[0]->dimsize = dim1_size;
2921 geo_var_dim1.push_back((fd->dims)[0]);
2924 auto lat_dim =
new Dimension(data_dim1_name,dim1_size);
2925 geo_var_dim1.push_back(lat_dim);
2929 else if(sdmap->getGeoDimension()==field_lat_dim2_name){
2930 string data_dim2_name = sdmap->getDataDimension();
2931 int dim2_size = sd->obtain_dimsize_with_dimname(data_dim2_name);
2932 if((this->swaths).size() > 1)
2933 data_dim2_name = data_dim2_name+
"_"+sd->name;
2934 if(sd->GeoDim_in_vars ==
false && dim2_map_count == 0) {
2935 (fd->dims)[1]->name = data_dim2_name;
2936 (fd->dims)[1]->dimsize = dim2_size;
2937 geo_var_dim2.push_back((fd->dims)[1]);
2940 auto lon_dim =
new Dimension(data_dim2_name,dim2_size);
2941 geo_var_dim2.push_back(lon_dim);
2948 for(
int i = 0; i<lat_names.size();i++) {
2958void File::create_geo_vars(SwathDataset* sd,Field *orig_lat,Field*orig_lon,
2959 const vector<string>& lat_names,
const vector<string>& lon_names,
2960 vector<Dimension*>&geo_var_dim1,vector<Dimension*>&geo_var_dim2)
throw(Exception){
2971 for (vector<Field *>::iterator i = sd->geofields.begin();
2972 i!=sd->geofields.end();++i) {
2973 if((*i)->name ==
"Latitude")
2975 else if((*i)->name ==
"Longitude")
2981 string ll_ori_dim0_name = (orig_lon->dims)[0]->name;
2982 string ll_ori_dim1_name = (orig_lon->dims)[1]->name;
2983 int dmap_offset = 0;
2985 if(sd->GeoDim_in_vars ==
false) {
2989 (orig_lat->dims)[0]->name = geo_var_dim1[0]->name;
2990 (orig_lat->dims)[0]->dimsize = geo_var_dim1[0]->dimsize;
2991 (orig_lat->dims)[1]->name = geo_var_dim2[0]->name;
2992 (orig_lat->dims)[1]->dimsize = geo_var_dim2[0]->dimsize;
2996 (orig_lon->dims)[0]->name = geo_var_dim1[0]->name;
2997 (orig_lon->dims)[0]->dimsize = geo_var_dim1[0]->dimsize;
2998 (orig_lon->dims)[1]->name = geo_var_dim2[0]->name;
2999 (orig_lon->dims)[1]->dimsize = geo_var_dim2[0]->dimsize;
3000 string ll_datadim0_name = geo_var_dim1[0]->name;
3001 string ll_datadim1_name = geo_var_dim2[0]->name;
3002 if(this->swaths.size() >1) {
3003 string prefix_remove =
"_"+sd->name;
3004 ll_datadim0_name = ll_datadim0_name.substr(0,ll_datadim0_name.size()-prefix_remove.size());
3005 ll_datadim1_name = ll_datadim1_name.substr(0,ll_datadim1_name.size()-prefix_remove.size());
3009 if(
false == sd->obtain_dmap_offset_inc(ll_ori_dim0_name,ll_datadim0_name,dmap_offset,dmap_inc)){
3010 throw5(
"Cannot retrieve dimension map offset and inc ",sd->name,
3011 orig_lon->name,ll_ori_dim0_name,ll_datadim0_name);
3013 orig_lon->ll_dim0_inc = dmap_inc;
3014 orig_lon->ll_dim0_offset = dmap_offset;
3015 orig_lat->ll_dim0_inc = dmap_inc;
3016 orig_lat->ll_dim0_offset = dmap_offset;
3018 if(
false == sd->obtain_dmap_offset_inc(ll_ori_dim1_name,ll_datadim1_name,dmap_offset,dmap_inc)){
3019 throw5(
"Cannot retrieve dimension map offset and inc ",sd->name,
3020 orig_lon->name,ll_ori_dim1_name,ll_datadim1_name);
3022 orig_lon->ll_dim1_inc = dmap_inc;
3023 orig_lon->ll_dim1_offset = dmap_offset;
3024 orig_lat->ll_dim1_inc = dmap_inc;
3025 orig_lat->ll_dim1_offset = dmap_offset;
3027cerr<<
"orig_lon "<<orig_lon->name <<endl;
3028cerr<<
"orig_lon dim1 inc "<<orig_lon->ll_dim1_inc<<endl;
3029cerr<<
"orig_lon dim1 offset "<<orig_lon->ll_dim1_offset<<endl;
3030cerr<<
"orig_lon dim0 inc "<<orig_lon->ll_dim0_inc<<endl;
3031cerr<<
"orig_lon dim0 offset "<<orig_lon->ll_dim0_offset<<endl;
3039 if((this->swaths).size() >1) {
3040 (orig_lon->dims)[0]->name = (orig_lon->dims)[0]->name +
"_" + sd->name;
3041 (orig_lon->dims)[1]->name = (orig_lon->dims)[1]->name +
"_" + sd->name;
3046 if((this->swaths).size()>1) {
3047 orig_lat->name = lat_names[0];
3048 orig_lon->name = lon_names[0];
3052 orig_lat->fieldtype = 1;
3053 orig_lon->fieldtype = 2;
3056 for (
int i = 1; i <lat_names.size();i++) {
3057 auto newlat =
new Field();
3058 newlat->name = lat_names[i];
3059 (newlat->dims).push_back(geo_var_dim1[i]);
3060 (newlat->dims).push_back(geo_var_dim2[i]);
3061 newlat->fieldtype = 1;
3063 newlat->type = orig_lat->type;
3064 auto newlon =
new Field();
3065 newlon->name = lon_names[i];
3069 new Dimension(geo_var_dim1[i]->name,geo_var_dim1[i]->dimsize);
3071 new Dimension(geo_var_dim2[i]->name,geo_var_dim2[i]->dimsize);
3072 (newlon->dims).push_back(lon_dim1);
3073 (newlon->dims).push_back(lon_dim2);
3074 newlon->fieldtype = 2;
3076 newlon->type = orig_lon->type;
3078 string ll_datadim0_name = geo_var_dim1[i]->name;
3079 string ll_datadim1_name = geo_var_dim2[i]->name;
3080 if(this->swaths.size() >1) {
3081 string prefix_remove =
"_"+sd->name;
3082 ll_datadim0_name = ll_datadim0_name.substr(0,ll_datadim0_name.size()-prefix_remove.size());
3083 ll_datadim1_name = ll_datadim1_name.substr(0,ll_datadim1_name.size()-prefix_remove.size());
3087 if(
false == sd->obtain_dmap_offset_inc(ll_ori_dim0_name,ll_datadim0_name,dmap_offset,dmap_inc)){
3088 throw5(
"Cannot retrieve dimension map offset and inc ",sd->name,
3089 newlon->name,ll_ori_dim0_name,ll_datadim0_name);
3091 newlon->ll_dim0_inc = dmap_inc;
3092 newlon->ll_dim0_offset = dmap_offset;
3093 newlat->ll_dim0_inc = dmap_inc;
3094 newlat->ll_dim0_offset = dmap_offset;
3095 if(
false == sd->obtain_dmap_offset_inc(ll_ori_dim1_name,ll_datadim1_name,dmap_offset,dmap_inc)){
3096 throw5(
"Cannot retrieve dimension map offset and inc ",sd->name,
3097 newlon->name,ll_ori_dim0_name,ll_datadim1_name);
3099 newlon->ll_dim1_inc = dmap_inc;
3100 newlon->ll_dim1_offset = dmap_offset;
3101 newlat->ll_dim1_inc = dmap_inc;
3102 newlat->ll_dim1_offset = dmap_offset;
3104cerr<<
"newlon "<<newlon->name <<endl;
3105cerr<<
"newlon dim1 inc "<<newlon->ll_dim1_inc<<endl;
3106cerr<<
"newlon dim1 offset "<<newlon->ll_dim1_offset<<endl;
3107cerr<<
"newlon dim0 inc "<<newlon->ll_dim0_inc<<endl;
3108cerr<<
"newlon dim0 offset "<<newlon->ll_dim0_offset<<endl;
3110 sd->geofields.push_back(newlat);
3111 sd->geofields.push_back(newlon);
3115for (vector<Field *>::const_iterator j =
3116 sd->getGeoFields().begin();
3117 j != sd->getGeoFields().end(); ++j) {
3118cerr<<
"Field name: "<<(*j)->name <<endl;
3119cerr<<
"Dimension name and size: "<<endl;
3120 for(vector<Dimension *>::const_iterator k=
3121 (*j)->getDimensions().begin();k!=(*j)->getDimensions().end();++k)
3122cerr<<(*k)->getName() <<
": "<<(*k)->getSize() <<endl;
3131void File::update_swath_dims_for_dimmap(SwathDataset* sd,
const std::vector<Dimension*>&geo_var_dim1,
3132 const std::vector<Dimension*>&geo_var_dim2) {
3137 for (
const auto &gfield:sd->getGeoFields()) {
3139 if(gfield->fieldtype == 1 || gfield->fieldtype == 2)
3141 for (
const auto &dim:gfield->getDimensions()) {
3142 string new_dim_name = dim->name +
"_"+sd->name;
3143 if (find_dim_in_dims(geo_var_dim1,new_dim_name) ||
3144 find_dim_in_dims(geo_var_dim2,new_dim_name))
3145 dim->name = new_dim_name;
3149 for (
const auto &dfield:sd->getDataFields()){
3150 for (
const auto &dim:dfield->getDimensions()) {
3151 string new_dim_name = dim->name +
"_"+sd->name;
3152 if(find_dim_in_dims(geo_var_dim1,new_dim_name) ||
3153 find_dim_in_dims(geo_var_dim2,new_dim_name))
3154 dim->name = new_dim_name;
3159 for (
const auto &dim:sd->getDimensions()) {
3160 string new_dim_name = dim->name +
"_"+sd->name;
3161 if(find_dim_in_dims(geo_var_dim1,new_dim_name) ||
3162 find_dim_in_dims(geo_var_dim2,new_dim_name))
3163 dim->name = new_dim_name;
3173void File::create_swath_latlon_dim_cvar_map_for_dimmap(SwathDataset* sd, Field* ori_lat, Field* ori_lon)
throw(Exception) {
3175 bool one_swath = ((this->swaths).size() == 1);
3177 int num_extra_lat_lon_pairs = 0;
3180if(sd->GeoDim_in_vars ==
true)
3181 cerr<<
" swath name is "<<sd->name <<endl;
3185 if(sd->GeoDim_in_vars ==
false)
3186 num_extra_lat_lon_pairs--;
3188 num_extra_lat_lon_pairs += (sd->num_map)/2;
3190 vector<string> lat_names;
3191 create_geo_varnames_list(lat_names,sd->name,ori_lat->name,num_extra_lat_lon_pairs,one_swath);
3192 vector<string>lon_names;
3193 create_geo_varnames_list(lon_names,sd->name,ori_lon->name,num_extra_lat_lon_pairs,one_swath);
3194 vector<Dimension*> geo_var_dim1;
3195 vector<Dimension*> geo_var_dim2;
3198 create_geo_dim_var_maps(sd, ori_lat,lat_names,lon_names,geo_var_dim1,geo_var_dim2);
3199 create_geo_vars(sd,ori_lat,ori_lon,lat_names,lon_names,geo_var_dim1,geo_var_dim2);
3203 if((this->swaths).size() >1)
3204 update_swath_dims_for_dimmap(sd,geo_var_dim1,geo_var_dim2);
3211void File::Prepare(
const char *eosfile_path)
throw(Exception)
3219 auto numgrid = (
int)(this->grids.size());
3220 auto numswath = (
int)(this->swaths.size());
3223 throw2(
"the number of grid is less than 0", eosfile_path);
3229 string DIMXNAME = this->get_geodim_x_name();
3231 string DIMYNAME = this->get_geodim_y_name();
3233 string LATFIELDNAME = this->get_latfield_name();
3235 string LONFIELDNAME = this->get_lonfield_name();
3238 string GEOGRIDNAME =
"location";
3243 check_onelatlon_grids();
3246 for (
const auto &grid:this->grids)
3247 handle_one_grid_zdim(grid);
3251 if (
true == this->onelatlon)
3252 handle_onelatlon_grids();
3254 for (
const auto &grid:this->grids) {
3258 grid->setDimxName(DIMXNAME);
3259 grid->setDimyName(DIMYNAME);
3262 handle_one_grid_latlon(grid);
3267 handle_grid_dim_cvar_maps();
3270 handle_grid_coards();
3273 update_grid_field_corrected_dims();
3276 handle_grid_cf_attrs();
3279 handle_grid_SOM_projection();
3285 for (
const auto& grid:this->grids)
3286 grid->SetScaleType(grid->name);
3294 check_swath_dimmap(numswath);
3297 check_dm_geo_dims_in_vars();
3300 check_swath_dimmap_bk_compat(numswath);
3303 create_swath_latlon_dim_cvar_map();
3306 create_swath_nonll_dim_cvar_map();
3309 handle_swath_dim_cvar_maps();
3313 handle_swath_cf_attrs();
3316 for (
const auto &swath:this->swaths)
3317 swath->SetScaleType(swath->name);
3327void correct_unlimited_missing_zdim(GridDataset* gdset)
throw(Exception) {
3329 for (vector<Field *>::const_iterator j =
3330 gdset->getDataFields().begin();
3331 j != gdset->getDataFields().end(); ++j) {
3334 if ((*j)->getRank()==1 && (*j)->){
3344bool File::check_special_1d_grid() throw(Exception) {
3346 auto numgrid = (
int)(this->grids.size());
3347 auto numswath = (
int)(this->swaths.size());
3349 if (numgrid != 1 || numswath != 0)
3353 string DIMXNAME = this->get_geodim_x_name();
3354 string DIMYNAME = this->get_geodim_y_name();
3356 if(DIMXNAME !=
"XDim" || DIMYNAME !=
"YDim")
3359 int var_dimx_size = 0;
3360 int var_dimy_size = 0;
3362 GridDataset *mygrid = (this->grids)[0];
3364 int field_xydim_flag = 0;
3365 for (
const auto &dfield:mygrid->getDataFields()) {
3366 if(1==dfield->rank) {
3367 if(dfield->name ==
"XDim"){
3369 var_dimx_size = (dfield->getDimensions())[0]->getSize();
3371 if(dfield->name ==
"YDim"){
3373 var_dimy_size = (dfield->getDimensions())[0]->getSize();
3376 if(2==field_xydim_flag)
3380 if(field_xydim_flag !=2)
3384 int xdimsize = mygrid->getInfo().getX();
3385 int ydimsize = mygrid->getInfo().getY();
3387 if(var_dimx_size != xdimsize || var_dimy_size != ydimsize)
3395bool File::check_ll_in_coords(
const string& vname)
throw(Exception) {
3397 bool ret_val =
false;
3398 for (
const auto &swath:this->swaths) {
3399 for (
const auto &gfield:swath->getGeoFields()) {
3401 if(gfield->fieldtype == 1 || gfield->fieldtype == 2) {
3402 if(gfield->getNewName() == vname) {
3410 for (
const auto &dfield:swath->getDataFields()) {
3413 if(dfield->fieldtype == 1 || dfield->fieldtype == 2) {
3414 if(dfield->getNewName() == vname) {
3435void Dataset::SetScaleType(
const string & EOS2ObjName)
throw(Exception) {
3443 vector<string> modis_multi_scale_type;
3444 modis_multi_scale_type.push_back(
"L1B");
3445 modis_multi_scale_type.push_back(
"GEO");
3446 modis_multi_scale_type.push_back(
"BRDF");
3447 modis_multi_scale_type.push_back(
"0.05Deg");
3448 modis_multi_scale_type.push_back(
"Reflectance");
3449 modis_multi_scale_type.push_back(
"MOD17A2");
3450 modis_multi_scale_type.push_back(
"North");
3451 modis_multi_scale_type.push_back(
"South");
3452 modis_multi_scale_type.push_back(
"MOD_Swath_Sea_Ice");
3453 modis_multi_scale_type.push_back(
"MOD_Grid_MOD15A2");
3454 modis_multi_scale_type.push_back(
"MOD_Grid_MOD16A2");
3455 modis_multi_scale_type.push_back(
"MOD_Grid_MOD16A3");
3456 modis_multi_scale_type.push_back(
"MODIS_NACP_LAI");
3458 vector<string> modis_div_scale_type;
3461 modis_div_scale_type.push_back(
"VI");
3462 modis_div_scale_type.push_back(
"1km_2D");
3463 modis_div_scale_type.push_back(
"L2g_2d");
3464 modis_div_scale_type.push_back(
"CMG");
3465 modis_div_scale_type.push_back(
"MODIS SWATH TYPE L2");
3470 string modis_eq_scale_type =
"LST";
3471 string modis_equ_scale_lst_group1=
"MODIS_Grid_8Day_1km_LST21";
3472 string modis_equ_scale_lst_group2=
"MODIS_Grid_Daily_1km_LST21";
3474 string modis_divequ_scale_group =
"MODIS_Grid";
3475 string modis_div_scale_group =
"MOD_Grid";
3479 string modis_equ_scale_group =
"MODIS_Grid_1km_2D";
3481 if(EOS2ObjName==
"mod05" || EOS2ObjName==
"mod06" || EOS2ObjName==
"mod07"
3482 || EOS2ObjName==
"mod08" || EOS2ObjName==
"atml2")
3484 scaletype = MODIS_MUL_SCALE;
3498 if(EOS2ObjName.find(
"MOD")==0 || EOS2ObjName.find(
"mod")==0)
3500 size_t pos = EOS2ObjName.rfind(modis_eq_scale_type);
3501 if(pos != string::npos &&
3502 (pos== (EOS2ObjName.length()-modis_eq_scale_type.length())))
3504 scaletype = MODIS_EQ_SCALE;
3508 for(
unsigned int k=0; k<modis_multi_scale_type.size(); k++)
3510 pos = EOS2ObjName.rfind(modis_multi_scale_type[k]);
3511 if (pos !=string::npos &&
3512 (pos== (EOS2ObjName.length()-modis_multi_scale_type[k].length())))
3514 scaletype = MODIS_MUL_SCALE;
3519 for(
unsigned int k=0; k<modis_div_scale_type.size(); k++)
3521 pos = EOS2ObjName.rfind(modis_div_scale_type[k]);
3522 if (pos != string::npos &&
3523 (pos==(EOS2ObjName.length()-modis_div_scale_type[k].length()))){
3524 scaletype = MODIS_DIV_SCALE;
3529 if (EOS2ObjName !=
"MODIS_Grid_1km_2D")
3536 pos = EOS2ObjName.find(modis_divequ_scale_group);
3542 size_t eq_scale_pos = EOS2ObjName.find(modis_equ_scale_group)
3543 *EOS2ObjName.find(modis_equ_scale_lst_group1)
3544 *EOS2ObjName.find(modis_equ_scale_lst_group2);
3545 if (0 == eq_scale_pos)
3546 scaletype = MODIS_EQ_SCALE;
3548 scaletype = MODIS_DIV_SCALE;
3551 size_t div_scale_pos = EOS2ObjName.find(modis_div_scale_group);
3554 if ( 0 == div_scale_pos)
3555 scaletype = MODIS_DIV_SCALE;
3561 if (EOS2ObjName ==
"VIP_CMG_GRID")
3562 scaletype = MODIS_DIV_SCALE;
3565int Dataset::obtain_dimsize_with_dimname(
const string & dimname) {
3567 for(vector<Dimension *>::const_iterator k=
3568 this->getDimensions().begin();k!=this->getDimensions().end();++k){
3569 if((*k)->name == dimname){
3570 ret_value = (*k)->dimsize;
3580 for_each(this->dims.begin(), this->dims.end(), delete_elem());
3581 for_each(this->correcteddims.begin(), this->correcteddims.end(), delete_elem());
3587 for_each(this->dims.begin(), this->dims.end(), delete_elem());
3588 for_each(this->datafields.begin(), this->datafields.end(),
3590 for_each(this->attrs.begin(), this->attrs.end(), delete_elem());
3594void Dataset::ReadDimensions(int32 (*entries)(int32, int32, int32 *),
3595 int32 (*inq)(int32,
char *, int32 *),
3596 vector<Dimension *> &d_dims)
throw(Exception)
3606 if ((numdims = entries(this->datasetid, HDFE_NENTDIM, &bufsize)) == -1)
3607 throw2(
"dimension entry", this->name);
3611 vector<char> namelist;
3612 vector<int32> dimsize;
3614 namelist.resize(bufsize + 1);
3615 dimsize.resize(numdims);
3618 if (inq(this->datasetid, namelist.data(), dimsize.data()) == -1)
3619 throw2(
"inquire dimension", this->name);
3621 vector<string> dimnames;
3627 for (vector<string>::const_iterator i = dimnames.begin();
3628 i != dimnames.end(); ++i) {
3629 Dimension *dim =
new Dimension(*i, dimsize[count]);
3630 d_dims.push_back(dim);
3637void Dataset::ReadFields(int32 (*entries)(int32, int32, int32 *),
3638 int32 (*inq)(int32,
char *, int32 *, int32 *),
3640 (int32,
char *, int32 *, int32 *, int32 *,
char *),
3642 (int32,
char *, int32 *, int32 *, int32 *, VOIDP),
3643 intn (*getfill)(int32,
char *, VOIDP),
3645 vector<Field *> &fields)
throw(Exception)
3649 int32 numfields = 0;
3656 if ((numfields = entries(this->datasetid, geofield ?
3657 HDFE_NENTGFLD : HDFE_NENTDFLD, &bufsize)) == -1)
3658 throw2(
"field entry", this->name);
3662 if (numfields > 0) {
3663 vector<char> namelist;
3665 namelist.resize(bufsize + 1);
3668 if (inq(this->datasetid, namelist.data(),
nullptr,
nullptr) == -1)
3669 throw2(
"inquire field", this->name);
3671 vector<string> fieldnames;
3676 for (vector<string>::const_iterator i = fieldnames.begin();
3677 i != fieldnames.end(); ++i) {
3679 Field *field =
new Field();
3680 if(field ==
nullptr)
3681 throw1(
"The field is nullptr");
3684 bool throw_error =
false;
3693 if ((fldinfo(this->datasetid,
3694 const_cast<char *
>(field->name.c_str()),
3695 &field->rank, dimsize, &field->type, dimlist)) == -1){
3696 string fieldname_for_eh = field->name;
3698 err_msg =
"Obtain field info error for field name "+fieldname_for_eh;
3701 if(
false == throw_error) {
3703 vector<string> dimnames;
3707 if ((
int)dimnames.size() != field->rank) {
3709 err_msg =
"Dimension names size is not consistent with field rank. ";
3710 err_msg +=
"Field name is "+field->name;
3713 for (
int k = 0; k < field->rank; ++k) {
3714 Dimension *dim =
new Dimension(dimnames[k], dimsize[k]);
3715 field->dims.push_back(dim);
3719 field->filler.resize(DFKNTsize(field->type));
3720 if (getfill(this->datasetid,
3721 const_cast<char *
>(field->name.c_str()),
3722 &field->filler[0]) == -1)
3723 field->filler.clear();
3726 fields.push_back(field);
3730 if(
true == throw_error) {
3740void Dataset::ReadAttributes(int32 (*inq)(int32,
char *, int32 *),
3741 intn (*attrinfo)(int32,
char *, int32 *, int32 *),
3742 intn (*readattr)(int32,
char *, VOIDP),
3743 vector<Attribute *> &obj_attrs)
throw(Exception)
3752 if ((numattrs = inq(this->datasetid,
nullptr, &bufsize)) == -1)
3753 throw2(
"inquire attribute", this->name);
3757 vector<char> namelist;
3759 namelist.resize(bufsize + 1);
3762 if (inq(this->datasetid, namelist.data(), &bufsize) == -1)
3763 throw2(
"inquire attribute", this->name);
3765 vector<string> attrnames;
3770 for (vector<string>::const_iterator i = attrnames.begin();
3771 i != attrnames.end(); ++i) {
3772 Attribute *attr =
new Attribute();
3779 if (attrinfo(this->datasetid,
const_cast<char *
>
3780 (attr->name.c_str()), &attr->type, &count) == -1) {
3782 throw3(
"attribute info", this->name, attr->name);
3785 attr->count = count/DFKNTsize(attr->type);
3786 attr->value.resize(count);
3793 if (readattr(this->datasetid,
3794 const_cast<char *
>(attr->name.c_str()),
3795 &attr->value[0]) == -1) {
3797 throw3(
"read attribute", this->name, attr->name);
3801 obj_attrs.push_back(attr);
3807GridDataset::~GridDataset()
3809 if (this->datasetid != -1){
3810 GDdetach(this->datasetid);
3815GridDataset * GridDataset::Read(int32 fd,
const string &gridname)
3819 bool GD_fun_err =
false;
3820 GridDataset *grid =
new GridDataset(gridname);
3823 if ((grid->datasetid = GDattach(fd,
const_cast<char *
>(gridname.c_str())))
3825 err_msg =
"attach grid";
3833 Info &info = grid->info;
3834 if (GDgridinfo(grid->datasetid, &info.xdim, &info.ydim, info.upleft,
3835 info.lowright) == -1) {
3836 err_msg =
"grid info";
3844 Projection &proj = grid->proj;
3845 if (GDprojinfo(grid->datasetid, &proj.code, &proj.zone, &proj.sphere,
3846 proj.param) == -1) {
3847 err_msg =
"projection info";
3851 if (GDpixreginfo(grid->datasetid, &proj.pix) == -1) {
3852 err_msg =
"pixreg info";
3856 if (GDorigininfo(grid->datasetid, &proj.origin) == -1){
3857 err_msg =
"origin info";
3863 if(
true == GD_fun_err){
3865 throw2(err_msg,gridname);
3870 grid->ReadDimensions(GDnentries, GDinqdims, grid->dims);
3873 grid->ReadFields(GDnentries, GDinqfields, GDfieldinfo, GDreadfield,
3874 GDgetfillvalue,
false, grid->datafields);
3877 grid->ReadAttributes(GDinqattrs, GDattrinfo, GDreadattr, grid->attrs);
3887GridDataset::Calculated & GridDataset::getCalculated()
const
3889 if (this->calculated.grid !=
this)
3890 this->calculated.grid =
this;
3891 return this->calculated;
3894bool GridDataset::Calculated::isYDimMajor() throw(Exception)
3896 this->DetectMajorDimension();
3897 return this->ydimmajor;
3901bool GridDataset::Calculated::isOrthogonal() throw(Exception)
3904 this->ReadLongitudeLatitude();
3905 return this->orthogonal;
3909int GridDataset::Calculated::DetectFieldMajorDimension() throw(Exception)
3914 for (vector<Field *>::const_iterator i =
3915 this->grid->getDataFields().begin();
3916 i != this->grid->getDataFields().end(); ++i) {
3918 int xdimindex = -1, ydimindex = -1, index = 0;
3921 for (vector<Dimension *>::const_iterator j =
3922 (*i)->getDimensions().begin();
3923 j != (*i)->getDimensions().end(); ++j) {
3924 if ((*j)->getName() == this->grid->dimxname)
3926 else if ((*j)->getName() == this->grid->dimyname)
3930 if (xdimindex == -1 || ydimindex == -1)
3933 int major = ydimindex < xdimindex ? 1 : 0;
3944 else if (ym != major)
3945 throw2(
"inconsistent XDim, YDim order", this->grid->getName());
3951 throw2(
"lack of data fields", this->grid->getName());
3956void GridDataset::Calculated::DetectMajorDimension() throw(Exception)
3964 for (vector<Field *>::const_iterator i =
3965 this->grid->getDataFields().begin();
3966 i != this->grid->getDataFields().end(); ++i) {
3968 int xdimindex = -1, ydimindex = -1, index = 0;
3971 for (vector<Dimension *>::const_iterator j =
3972 (*i)->getDimensions().begin();
3973 j != (*i)->getDimensions().end(); ++j) {
3974 if ((*j)->getName() == this->grid->dimxname)
3976 else if ((*j)->getName() == this->grid->dimyname)
3980 if (xdimindex == -1 || ydimindex == -1)
3985 if(this->grid->getProjection().getCode() == GCTP_LAMAZ)
3988 major = ydimindex < xdimindex ? 1 : 0;
3999 else if (ym != major)
4000 throw2(
"inconsistent XDim, YDim order", this->grid->getName());
4005 throw2(
"lack of data fields", this->grid->getName());
4006 this->ydimmajor = ym != 0;
4015static bool IsDisjoint(
const vector<Field *> &l,
4016 const vector<Field *> &r)
4018 for (vector<Field *>::const_iterator i = l.begin(); i != l.end(); ++i)
4020 if (find(r.begin(), r.end(), *i) != r.end())
4028static bool IsDisjoint(vector<pair<Field *, string> > &l,
const vector<Field *> &r)
4030 for (vector<pair<Field *, string> >::const_iterator i =
4031 l.begin(); i != l.end(); ++i) {
4032 if (find(r.begin(), r.end(), i->first) != r.end())
4039static bool IsSubset(
const vector<Field *> &s,
const vector<Field *> &b)
4041 for (vector<Field *>::const_iterator i = s.begin(); i != s.end(); ++i)
4043 if (find(b.begin(), b.end(), *i) == b.end())
4050static bool IsSubset(vector<pair<Field *, string> > &s,
const vector<Field *> &b)
4052 for (vector<pair<Field *, string> >::const_iterator i
4053 = s.begin(); i != s.end(); ++i) {
4054 if (find(b.begin(), b.end(), i->first) == b.end())
4062SwathDataset::~SwathDataset()
4064 if (this->datasetid != -1) {
4065 SWdetach(this->datasetid);
4068 for_each(this->dimmaps.begin(), this->dimmaps.end(), delete_elem());
4069 for_each(this->indexmaps.begin(), this->indexmaps.end(),
4072 for_each(this->geofields.begin(), this->geofields.end(),
4079SwathDataset * SwathDataset::Read(int32 fd,
const string &swathname)
4082 SwathDataset *swath =
new SwathDataset(swathname);
4083 if(swath ==
nullptr)
4084 throw1(
"Cannot allocate HDF5 Swath object");
4087 if ((swath->datasetid = SWattach(fd,
4088 const_cast<char *
>(swathname.c_str())))
4091 throw2(
"attach swath", swathname);
4099 swath->ReadDimensions(SWnentries, SWinqdims, swath->dims);
4102 swath->ReadFields(SWnentries, SWinqdatafields, SWfieldinfo, SWreadfield,
4103 SWgetfillvalue,
false, swath->datafields);
4106 swath->ReadFields(SWnentries, SWinqgeofields, SWfieldinfo, SWreadfield,
4107 SWgetfillvalue,
true, swath->geofields);
4110 swath->ReadAttributes(SWinqattrs, SWattrinfo, SWreadattr, swath->attrs);
4113 swath->set_num_map(swath->ReadDimensionMaps(swath->dimmaps));
4116 swath->ReadIndexMaps(swath->indexmaps);
4128int SwathDataset::ReadDimensionMaps(vector<DimensionMap *> &swath_dimmaps)
4131 int32 nummaps, bufsize;
4134 if ((nummaps = SWnentries(this->datasetid, HDFE_NENTMAP, &bufsize)) == -1)
4135 throw2(
"dimmap entry", this->name);
4146 vector<char> namelist;
4147 vector<int32> offset, increment;
4149 namelist.resize(bufsize + 1);
4150 offset.resize(nummaps);
4151 increment.resize(nummaps);
4152 if (SWinqmaps(this->datasetid, namelist.data(), offset.data(), increment.data())
4154 throw2(
"inquire dimmap", this->name);
4156 vector<string> mapnames;
4159 for (vector<string>::const_iterator i = mapnames.begin();
4160 i != mapnames.end(); ++i) {
4161 vector<string> parts;
4163 if (parts.size() != 2)
4164 throw3(
"dimmap part", parts.size(),
4167 DimensionMap *dimmap =
new DimensionMap(parts[0], parts[1],
4170 swath_dimmaps.push_back(dimmap);
4178void SwathDataset::ReadIndexMaps(vector<IndexMap *> &swath_indexmaps)
4181 int32 numindices, bufsize;
4183 if ((numindices = SWnentries(this->datasetid, HDFE_NENTIMAP, &bufsize)) ==
4185 throw2(
"indexmap entry", this->name);
4186 if (numindices > 0) {
4188 vector<char> namelist;
4190 namelist.resize(bufsize + 1);
4191 if (SWinqidxmaps(this->datasetid, namelist.data(),
nullptr) == -1)
4192 throw2(
"inquire indexmap", this->name);
4194 vector<string> mapnames;
4196 for (vector<string>::const_iterator i = mapnames.begin();
4197 i != mapnames.end(); ++i) {
4198 IndexMap *indexmap =
new IndexMap();
4199 vector<string> parts;
4201 if (parts.size() != 2)
4202 throw3(
"indexmap part", parts.size(),
4204 indexmap->geo = parts[0];
4205 indexmap->data = parts[1];
4210 SWdiminfo(this->datasetid,
4211 const_cast<char *
>(indexmap->geo.c_str())))
4213 throw3(
"dimension info", this->name, indexmap->geo);
4214 indexmap->indices.resize(dimsize);
4215 if (SWidxmapinfo(this->datasetid,
4216 const_cast<char *
>(indexmap->geo.c_str()),
4217 const_cast<char *
>(indexmap->data.c_str()),
4218 &indexmap->indices[0]) == -1)
4219 throw4(
"indexmap info", this->name, indexmap->geo,
4223 swath_indexmaps.push_back(indexmap);
4229PointDataset::~PointDataset()
4233PointDataset * PointDataset::Read(int32 ,
const string &pointname)
4236 PointDataset *point =
new PointDataset(pointname);
4242bool Utility::ReadNamelist(
const char *path,
4243 int32 (*inq)(
char *,
char *, int32 *),
4244 vector<string> &names)
4246 char *fname =
const_cast<char *
>(path);
4250 if ((numobjs = inq(fname,
nullptr, &bufsize)) == -1)
return false;
4252 vector<char> buffer;
4253 buffer.resize(bufsize + 1);
4254 if (inq(fname, buffer.data(), &bufsize) == -1)
return false;
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)