8 #include "HDFEOS2ArrayGridGeoField.h"
18 #include "InternalErr.h"
20 #include "HDFCFUtil.h"
23 #include "errormacros.h"
25 #include <sys/types.h>
29 #include "BESH4MCache.h"
30 #include "HDF4RequestHandler.h"
35 #define SIGNED_BYTE_TO_INT32 1
39 int inv_init(
int insys,
int inzone,
double *inparm,
int indatum,
char *fn27,
char *fn83,
int *iflg,
int (*inv_trans[])(
double,
double,
double*,
double*));
40 int sominv(
double y,
double x,
double *lon,
double *lat);
44 HDFEOS2ArrayGridGeoField::read ()
46 BESDEBUG(
"h4",
"Coming to HDFEOS2ArrayGridGeoField read "<<endl);
51 string check_pass_fileid_key_str=
"H4.EnablePassFileID";
52 bool check_pass_fileid_key =
false;
53 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
55 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
59 if (rank < 1 || rank > 2) {
60 throw InternalErr (__FILE__, __LINE__,
"The rank of geo field is greater than 2, currently we don't support 3-D lat/lon cases.");
66 if (
true == condenseddim)
68 else if(4 == specialformat)
74 offset.resize(final_rank);
76 count.resize(final_rank);
78 step.resize(final_rank);
83 nelms = format_constraint (&offset[0], &step[0], &count[0]);
88 int32 (*openfunc) (
char *, intn);
89 int32 (*attachfunc) (int32,
char *);
90 intn (*detachfunc) (int32);
91 intn (*fieldinfofunc) (int32,
char *, int32 *, int32 *, int32 *,
char *);
92 intn (*readfieldfunc) (int32,
char *, int32 *, int32 *, int32 *,
void *);
96 attachfunc = GDattach;
97 detachfunc = GDdetach;
98 fieldinfofunc = GDfieldinfo;
99 readfieldfunc = GDreadfield;
100 datasetname = gridname;
117 string cache_fpath=
"";
118 bool use_cache =
false;
121 if(
true == check_pass_fileid_key)
124 gfid = openfunc (
const_cast < char *
>(filename.c_str ()), DFACC_READ);
127 eherr <<
"File " << filename.c_str () <<
" cannot be open.";
128 throw InternalErr (__FILE__, __LINE__, eherr.str ());
133 gridid = attachfunc (gfid,
const_cast < char *
>(datasetname.c_str ()));
137 eherr <<
"Grid " << datasetname.c_str () <<
" cannot be attached.";
138 throw InternalErr (__FILE__, __LINE__, eherr.str ());
141 if(
false == llflag) {
148 string check_eos_geo_cache_key =
"H4.EnableEOSGeoCacheFile";
149 bool enable_eos_geo_cache_key =
false;
150 enable_eos_geo_cache_key = HDFCFUtil::check_beskeys(check_eos_geo_cache_key);
153 if(
true == HDF4RequestHandler::get_enable_eosgeo_cachefile()) {
162 string bescachedir= llcache->getCacheDirFromConfig();
163 string bescacheprefix = llcache->getCachePrefixFromConfig();
164 unsigned long cachesize = llcache->getCacheSizeFromConfig();
167 string bescachedir = HDF4RequestHandler::get_cache_latlon_path();
168 string bescacheprefix = HDF4RequestHandler::get_cache_latlon_prefix();
169 long cachesize = HDF4RequestHandler::get_cache_latlon_size();
172 if((
"" == bescachedir)||(
""==bescacheprefix)||(cachesize <=0)){
174 throw InternalErr (__FILE__, __LINE__,
"Either the cached dir is empty or the prefix is NULL or the cache size is not set.");
178 if(stat(bescachedir.c_str(),&sb) !=0) {
180 string err_mesg=
"The cached directory " + bescachedir;
181 err_mesg = err_mesg +
" doesn't exist. ";
182 throw InternalErr(__FILE__,__LINE__,err_mesg);
186 if(
true == S_ISDIR(sb.st_mode)) {
187 if(access(bescachedir.c_str(),R_OK|W_OK|X_OK) == -1) {
190 string err_mesg=
"The cached directory " + bescachedir;
191 err_mesg = err_mesg +
" can NOT be read,written or executable.";
192 throw InternalErr(__FILE__,__LINE__,err_mesg);
198 string err_mesg=
"The cached directory " + bescachedir;
199 err_mesg = err_mesg +
" is not a directory.";
200 throw InternalErr(__FILE__,__LINE__,err_mesg);
206 string cache_fname=HDF4RequestHandler::get_cache_latlon_prefix();
209 r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
213 throw InternalErr (__FILE__, __LINE__,
"GDprojinfo failed");
217 if (GDgridinfo(gridid, &xdim, &ydim, upleft,
221 throw InternalErr (__FILE__, __LINE__,
"GDgridinfo failed");
226 r = GDpixreginfo (gridid, &pixreg);
231 eherr <<
"cannot obtain grid pixel registration info.";
232 throw InternalErr (__FILE__, __LINE__, eherr.str ());
237 r = GDorigininfo (gridid, &origin);
242 eherr <<
"cannot obtain grid origin info.";
243 throw InternalErr (__FILE__, __LINE__, eherr.str ());
248 cache_fname +=HDFCFUtil::get_int_str(projcode);
249 cache_fname +=HDFCFUtil::get_int_str(zone);
250 cache_fname +=HDFCFUtil::get_int_str(sphere);
251 cache_fname +=HDFCFUtil::get_int_str(pixreg);
252 cache_fname +=HDFCFUtil::get_int_str(origin);
258 cache_fname +=HDFCFUtil::get_int_str(ydim);
259 cache_fname +=HDFCFUtil::get_int_str(xdim);
263 cache_fname +=HDFCFUtil::get_int_str(xdim);
264 cache_fname +=HDFCFUtil::get_int_str(ydim);
269 cache_fname +=HDFCFUtil::get_double_str(upleft[0],17,6);
270 cache_fname +=HDFCFUtil::get_double_str(upleft[1],17,6);
271 cache_fname +=HDFCFUtil::get_double_str(lowright[0],17,6);
272 cache_fname +=HDFCFUtil::get_double_str(lowright[1],17,6);
275 for(
int ipar = 0; ipar<13;ipar++)
276 cache_fname+=HDFCFUtil::get_double_str(params[ipar],17,6);
279 cache_fpath = bescachedir +
"/"+ cache_fname;
281 cerr<<
"cache file path is "<<cache_fpath <<endl;
282 cerr<<
"obtain file path from BESMCache "<<endl;
283 cerr<<
"Name is "<<llcache->get_cache_file_name_h4(cache_fpath,
false) <<endl;
286 cerr<<
"after testing get_read_lock"<<endl;
291 int expected_file_size = 0;
292 if(GCTP_CEA == projcode || GCTP_GEO == projcode)
293 expected_file_size = (xdim+ydim)*
sizeof(
double);
294 else if(GCTP_SOM == projcode)
295 expected_file_size = xdim*ydim*NBLOCK*2*
sizeof(double);
297 expected_file_size = xdim*ydim*2*
sizeof(double);
300 bool latlon_from_cache = llcache->get_data_from_cache(cache_fpath, expected_file_size,fd);
302 if(
true == latlon_from_cache)
303 cerr<<
"the cached file exists: "<<endl;
305 cerr<<
"the cached file doesn't exist "<< endl;
307 if(
false == latlon_from_cache)
313 size_t offset_1d = 0;
320 if(GCTP_CEA == projcode|| GCTP_GEO== projcode) {
324 offset_1d = (fieldtype == 1) ?offset[0] :(ydim+offset[0]);
326 if(
true == ydimmajor) {
327 offset_1d = (fieldtype == 1) ?offset[0] :(ydim*
sizeof(
double)+offset[0]);
330 offset_1d = (fieldtype == 1) ?offset[0] :(xdim*
sizeof(
double)+offset[0]);
333 count_1d = 1+(count[0]-1)*step[0];
335 else if (GCTP_SOM == projcode) {
337 if(
true == ydimmajor) {
338 offset_1d = (fieldtype == 1)?(offset[0]*xdim*ydim+offset[1]*xdim+offset[2])
339 :(offset[0]*xdim*ydim+offset[1]*xdim+offset[2]+expected_file_size/2/
sizeof(double));
342 offset_1d = (fieldtype == 1)?(offset[0]*xdim*ydim+offset[1]*ydim+offset[2])
343 :(offset[0]*xdim*ydim+offset[1]*ydim+offset[2]+expected_file_size/2/
sizeof(double));
346 int total_count_dim0 = (count[0]-1)*step[0];
347 int total_count_dim1 = (count[1]-1)*step[1];
348 int total_count_dim2 = (count[2]-1)*step[2];
349 int total_dim1 = (
true ==ydimmajor)?ydim:xdim;
350 int total_dim2 = (
true ==ydimmajor)?xdim:ydim;
355 count_1d = 1 + total_count_dim0*total_dim1*total_dim2 + total_count_dim1*total_dim2 + total_count_dim2;
359 if (
true == ydimmajor)
360 offset_1d = (fieldtype == 1) ?(offset[0] * xdim + offset[1]):(expected_file_size/2/
sizeof(double)+offset[0]*xdim+offset[1]);
362 offset_1d = (fieldtype == 1) ?(offset[0] * ydim + offset[1]):(expected_file_size/2/
sizeof(double)+offset[0]*ydim+offset[1]);
365 int total_count_dim0 = (count[0]-1)*step[0];
366 int total_count_dim1 = (count[1]-1)*step[1];
367 int total_dim1 = (
true ==ydimmajor)?xdim:ydim;
369 count_1d = 1 + total_count_dim0*total_dim1 + total_count_dim1;
373 vector<double> latlon_1d;
374 latlon_1d.resize(count_1d);
380 off_t fpos = lseek(fd,
sizeof(
double)*offset_1d,SEEK_SET);
386 ssize_t read_size = HDFCFUtil::read_vector_from_file(fd,latlon_1d,
sizeof(
double));
388 if((-1 == read_size) || ((
size_t)read_size != count_1d*
sizeof(
double))) {
397 pFile = fopen(cache_fpath.c_str(),
"rb");
401 int ret_value = fseek(pFile,
sizeof(
double)*offset_1d,SEEK_SET);
407 cerr<<
"field name is "<<fieldname <<endl;
408 cerr<<
"fread is checked "<<endl;
409 cerr<<
"offset_1d is "<<offset_1d <<endl;
410 cerr<<
"count_1d is "<<count_1d <<endl;
413 ret_value = fread(&latlon_1d[0],
sizeof(
double),count_1d,pFile);
416 cerr<<
"fread fails "<<endl;
422 for(
int i =0;i<count_1d;i++)
423 cerr<<
"latlon_1d["<<i<<
"]"<<latlon_1d[i]<<endl;
427 for (
int i_rank = 0; i_rank<final_rank;i_rank++)
428 total_count = total_count*count[i_rank];
434 if(total_count == count_1d) {
435 set_value((dods_float64*)&latlon_1d[0],nelms);
441 vector<double>latlon;
442 latlon.resize(total_count);
445 if(GCTP_CEA == projcode|| GCTP_GEO== projcode) {
446 for (
int i = 0; i <total_count;i++)
447 latlon[i] = latlon_1d[i*step[0]];
450 else if (GCTP_SOM == projcode) {
452 for (
int i =0; i<count[0];i++)
453 for(
int j =0;j<count[1];j++)
454 for(
int k=0;k<count[2];k++)
455 latlon[i*count[1]*count[2]+j*count[2]+k]=(
true == ydimmajor)
456 ?latlon_1d[i*ydim*xdim*step[0]+j*xdim*step[1]+k*step[2]]
457 :latlon_1d[i*ydim*xdim*step[0]+j*ydim*step[1]+k*step[2]];
460 for (
int i =0; i<count[0];i++)
461 for(
int j =0;j<count[1];j++)
462 latlon[i*count[1]+j]=(
true == ydimmajor)
463 ?latlon_1d[i*xdim*step[0]+j*step[1]]
464 :latlon_1d[i*ydim*step[0]+j*step[1]];
468 set_value((dods_float64*)&latlon[0],nelms);
489 if(specialformat == 4) {
491 CalculateSOMLatLon(gridid, &offset[0], &count[0], &step[0], nelms,cache_fpath,use_cache);
504 vector<int32>offset32;
505 offset32.resize(rank);
507 vector<int32>count32;
508 count32.resize(rank);
517 getCorrectSubset (&offset[0], &count[0], &step[0], &offset32[0], &count32[0], &step32[0], condenseddim, ydimmajor, fieldtype, rank);
526 if (llflag ==
false) {
528 vector<float64>latlon;
529 latlon.resize(nelms);
533 if(projcode == -1 && specialformat !=3) {
537 r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
541 throw InternalErr (__FILE__, __LINE__,
"GDprojinfo failed");
545 if (GDgridinfo(gridid, &xdim, &ydim, upleft,
549 throw InternalErr (__FILE__, __LINE__,
"GDgridinfo failed");
555 if (GCTP_LAMAZ == projcode) {
557 vector<double>latlon_all;
558 latlon_all.resize(xdim*ydim*2);
560 CalculateLAMAZLatLon(gridid, fieldtype, &latlon[0], &latlon_all[0],&offset[0], &count[0], &step[0], use_cache);
561 if(
true == use_cache) {
564 llcache->write_cached_data(cache_fpath,xdim*ydim*2*
sizeof(
double),latlon_all);
573 set_value ((dods_float64 *) &latlon[0], nelms);
580 if (specialformat == 1) {
583 vector<double>latlon_all;
584 latlon_all.resize(xdim+ydim);
586 CalculateLargeGeoLatLon(gridid, fieldtype,&latlon[0], &latlon_all[0],&offset[0], &count[0], &step[0], nelms,use_cache);
587 if(
true == use_cache) {
590 llcache->write_cached_data(cache_fpath,(xdim+ydim)*
sizeof(
double),latlon_all);
594 if(HDFCFUtil::write_vector_to_file2(cache_fpath,latlon_all,
sizeof(
double)) != ((xdim+ydim)*
sizeof(
double))) {
595 if(remove(cache_fpath.c_str()) !=0) {
596 throw InternalErr(__FILE__,__LINE__,
"Cannot remove the cached file.");
608 set_value((dods_float64 *)&latlon[0],nelms);
616 else if (specialformat == 3) {
618 CalculateSpeLatLon (gridid, fieldtype, &latlon[0], &offset32[0], &count32[0], &step32[0]);
632 vector<double>latlon_all;
634 if(GCTP_GEO == projcode || GCTP_CEA == projcode)
635 latlon_all.resize(xdim+ydim);
637 latlon_all.resize(xdim*ydim*2);
639 CalculateLatLon (gridid, fieldtype, specialformat, &latlon[0],&latlon_all[0],
640 &offset32[0], &count32[0], &step32[0], nelms,use_cache);
641 if(
true == use_cache) {
644 size_t num_item_expected = 0;
645 if(GCTP_GEO == projcode || GCTP_CEA == projcode)
646 num_item_expected = xdim + ydim;
648 num_item_expected = xdim*ydim*2;
651 llcache->write_cached_data(cache_fpath,num_item_expected*
sizeof(
double),latlon_all);
657 if (speciallon && fieldtype == 2)
658 CorSpeLon(&latlon[0], nelms);
663 set_value ((dods_float64 *) &latlon[0], nelms);
671 vector<int32> tmp_dims;
672 tmp_dims.resize(rank);
674 char tmp_dimlist[1024];
679 r = fieldinfofunc (gridid,
const_cast < char *
>(fieldname.c_str ()),
680 &tmp_rank, &tmp_dims[0], &type, tmp_dimlist);
686 eherr <<
"Field " << fieldname.c_str () <<
" information cannot be obtained.";
687 throw InternalErr (__FILE__, __LINE__, eherr.str ());
698 r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
703 eherr <<
"Grid " << datasetname.c_str () <<
" information cannot be obtained.";
704 throw InternalErr (__FILE__, __LINE__, eherr.str ());
713 r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
718 eherr <<
"Grid " << datasetname.c_str () <<
" projection info. cannot be obtained.";
719 throw InternalErr (__FILE__, __LINE__, eherr.str ());
722 if (projcode != GCTP_GEO) {
729 r = readfieldfunc (gridid,
730 const_cast < char *
>(fieldname.c_str ()),
731 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
736 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
737 throw InternalErr (__FILE__, __LINE__, eherr.str ());
742 #ifndef SIGNED_BYTE_TO_INT32
743 set_value ((dods_byte *) &val[0], nelms);
746 newval.resize(nelms);
748 for (
int counter = 0; counter < nelms; counter++)
749 newval[counter] = (int32) (val[counter]);
751 set_value ((dods_int32 *) &newval[0], nelms);
762 r = readfieldfunc (gridid,
763 const_cast < char *
>(fieldname.c_str ()),
764 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
769 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
770 throw InternalErr (__FILE__, __LINE__, eherr.str ());
772 set_value ((dods_byte *) &val[0], nelms);
783 r = readfieldfunc (gridid,
784 const_cast < char *
>(fieldname.c_str ()),
785 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
790 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
791 throw InternalErr (__FILE__, __LINE__, eherr.str ());
794 set_value ((dods_int16 *) &val[0], nelms);
804 r = readfieldfunc (gridid,
805 const_cast < char *
>(fieldname.c_str ()),
806 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
811 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
812 throw InternalErr (__FILE__, __LINE__, eherr.str ());
815 set_value ((dods_uint16 *) &val[0], nelms);
824 r = readfieldfunc (gridid,
825 const_cast < char *
>(fieldname.c_str ()),
826 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
831 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
832 throw InternalErr (__FILE__, __LINE__, eherr.str ());
835 set_value ((dods_int32 *) &val[0], nelms);
844 r = readfieldfunc (gridid,
845 const_cast < char *
>(fieldname.c_str ()),
846 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
851 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
852 throw InternalErr (__FILE__, __LINE__, eherr.str ());
854 set_value ((dods_uint32 *) &val[0], nelms);
863 r = readfieldfunc (gridid,
864 const_cast < char *
>(fieldname.c_str ()),
865 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
870 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
871 throw InternalErr (__FILE__, __LINE__, eherr.str ());
874 set_value ((dods_float32 *) &val[0], nelms);
883 r = readfieldfunc (gridid,
884 const_cast < char *
>(fieldname.c_str ()),
885 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
890 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
891 throw InternalErr (__FILE__, __LINE__, eherr.str ());
894 set_value ((dods_float64 *) &val[0], nelms);
901 throw InternalErr (__FILE__, __LINE__,
"unsupported data type.");
920 r = GDgetfillvalue (gridid,
921 const_cast < char *
>(fieldname.c_str ()),
924 int ifillvalue = fillvalue;
926 vector <int8> temp_total_val;
928 temp_total_val.resize(xdim*ydim);
931 r = readfieldfunc(gridid,
932 const_cast < char *
>(fieldname.c_str ()),
933 NULL, NULL, NULL, (
void *)(&temp_total_val[0]));
939 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
940 throw InternalErr (__FILE__, __LINE__, eherr.str ());
945 HandleFillLatLon(temp_total_val, (int8*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
957 r = readfieldfunc (gridid,
958 const_cast < char *
>(fieldname.c_str ()),
959 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
964 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
965 throw InternalErr (__FILE__, __LINE__, eherr.str ());
969 if (speciallon && fieldtype == 2)
970 CorSpeLon ((int8 *) &val[0], nelms);
973 #ifndef SIGNED_BYTE_TO_INT32
974 set_value ((dods_byte *) &val[0], nelms);
977 newval.resize(nelms);
979 for (
int counter = 0; counter < nelms; counter++)
980 newval[counter] = (int32) (val[counter]);
982 set_value ((dods_int32 *) &newval[0], nelms);
997 r = GDgetfillvalue (gridid,
998 const_cast < char *
>(fieldname.c_str ()),
1003 int ifillvalue = fillvalue;
1004 vector <uint8> temp_total_val;
1005 temp_total_val.resize(xdim*ydim);
1007 r = readfieldfunc(gridid,
1008 const_cast < char *
>(fieldname.c_str ()),
1009 NULL, NULL, NULL, (
void *)(&temp_total_val[0]));
1014 ostringstream eherr;
1015 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
1016 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1020 HandleFillLatLon(temp_total_val, (uint8*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
1032 r = readfieldfunc (gridid,
1033 const_cast < char *
>(fieldname.c_str ()),
1034 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
1038 ostringstream eherr;
1039 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
1040 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1044 if (speciallon && fieldtype == 2)
1045 CorSpeLon ((uint8 *) &val[0], nelms);
1046 set_value ((dods_byte *) &val[0], nelms);
1056 int16 fillvalue = 0;
1058 r = GDgetfillvalue (gridid,
1059 const_cast < char *
>(fieldname.c_str ()),
1063 int ifillvalue = fillvalue;
1064 vector <int16> temp_total_val;
1065 temp_total_val.resize(xdim*ydim);
1067 r = readfieldfunc(gridid,
1068 const_cast < char *
>(fieldname.c_str ()),
1069 NULL, NULL, NULL, (
void *)(&temp_total_val[0]));
1074 ostringstream eherr;
1075 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
1076 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1080 HandleFillLatLon(temp_total_val, (int16*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
1092 r = readfieldfunc (gridid,
1093 const_cast < char *
>(fieldname.c_str ()),
1094 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
1098 ostringstream eherr;
1099 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
1100 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1105 if (speciallon && fieldtype == 2)
1106 CorSpeLon ((int16 *) &val[0], nelms);
1108 set_value ((dods_int16 *) &val[0], nelms);
1113 uint16 fillvalue = 0;
1117 r = GDgetfillvalue (gridid,
1118 const_cast < char *
>(fieldname.c_str ()),
1123 int ifillvalue = fillvalue;
1125 vector <uint16> temp_total_val;
1126 temp_total_val.resize(xdim*ydim);
1128 r = readfieldfunc(gridid,
1129 const_cast < char *
>(fieldname.c_str ()),
1130 NULL, NULL, NULL, (
void *)(&temp_total_val[0]));
1135 ostringstream eherr;
1136 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
1137 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1141 HandleFillLatLon(temp_total_val, (uint16*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
1152 r = readfieldfunc (gridid,
1153 const_cast < char *
>(fieldname.c_str ()),
1154 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
1158 ostringstream eherr;
1159 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
1160 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1164 if (speciallon && fieldtype == 2)
1165 CorSpeLon ((uint16 *) &val[0], nelms);
1167 set_value ((dods_uint16 *) &val[0], nelms);
1177 int32 fillvalue = 0;
1179 r = GDgetfillvalue (gridid,
1180 const_cast < char *
>(fieldname.c_str ()),
1184 int ifillvalue = fillvalue;
1186 vector <int32> temp_total_val;
1187 temp_total_val.resize(xdim*ydim);
1189 r = readfieldfunc(gridid,
1190 const_cast < char *
>(fieldname.c_str ()),
1191 NULL, NULL, NULL, (
void *)(&temp_total_val[0]));
1194 ostringstream eherr;
1197 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
1198 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1202 HandleFillLatLon(temp_total_val, (int32*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
1214 r = readfieldfunc (gridid,
1215 const_cast < char *
>(fieldname.c_str ()),
1216 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
1220 ostringstream eherr;
1221 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
1222 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1226 if (speciallon && fieldtype == 2)
1227 CorSpeLon ((int32 *) &val[0], nelms);
1229 set_value ((dods_int32 *) &val[0], nelms);
1239 uint32 fillvalue = 0;
1241 r = GDgetfillvalue (gridid,
1242 const_cast < char *
>(fieldname.c_str ()),
1247 int ifillvalue = (int)fillvalue;
1248 vector <uint32> temp_total_val;
1249 temp_total_val.resize(xdim*ydim);
1250 r = readfieldfunc(gridid,
1251 const_cast < char *
>(fieldname.c_str ()),
1252 NULL, NULL, NULL, (
void *)(&temp_total_val[0]));
1257 ostringstream eherr;
1258 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
1259 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1263 HandleFillLatLon(temp_total_val, (uint32*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
1275 r = readfieldfunc (gridid,
1276 const_cast < char *
>(fieldname.c_str ()),
1277 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
1281 ostringstream eherr;
1282 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
1283 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1287 if (speciallon && fieldtype == 2)
1288 CorSpeLon ((uint32 *) &val[0], nelms);
1290 set_value ((dods_uint32 *) &val[0], nelms);
1297 vector<float32> val;
1300 float32 fillvalue =0;
1301 r = GDgetfillvalue (gridid,
1302 const_cast < char *
>(fieldname.c_str ()),
1309 int ifillvalue =(int)fillvalue;
1311 vector <float32> temp_total_val;
1312 temp_total_val.resize(xdim*ydim);
1315 r = readfieldfunc(gridid,
1316 const_cast < char *
>(fieldname.c_str ()),
1317 NULL, NULL, NULL, (
void *)(&temp_total_val[0]));
1322 ostringstream eherr;
1323 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
1324 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1328 HandleFillLatLon(temp_total_val, (float32*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
1339 r = readfieldfunc (gridid,
1340 const_cast < char *
>(fieldname.c_str ()),
1341 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
1345 ostringstream eherr;
1346 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
1347 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1351 if (speciallon && fieldtype == 2)
1352 CorSpeLon ((float32 *) &val[0], nelms);
1354 set_value ((dods_float32 *) &val[0], nelms);
1361 vector<float64> val;
1364 float64 fillvalue = 0;
1365 r = GDgetfillvalue (gridid,
1366 const_cast < char *
>(fieldname.c_str ()),
1373 int ifillvalue = (int)fillvalue;
1374 vector <float64> temp_total_val;
1375 temp_total_val.resize(xdim*ydim);
1376 r = readfieldfunc(gridid,
1377 const_cast < char *
>(fieldname.c_str ()),
1378 NULL, NULL, NULL, (
void *)(&temp_total_val[0]));
1383 ostringstream eherr;
1384 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
1385 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1389 HandleFillLatLon(temp_total_val, (float64*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
1401 r = readfieldfunc (gridid,
1402 const_cast < char *
>(fieldname.c_str ()),
1403 &offset32[0], &step32[0], &count32[0], (
void*)(&val[0]));
1407 ostringstream eherr;
1408 eherr <<
"field " << fieldname.c_str () <<
"cannot be read.";
1409 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1413 if (speciallon && fieldtype == 2)
1414 CorSpeLon ((float64 *) &val[0], nelms);
1416 set_value ((dods_float64 *) &val[0], nelms);
1423 throw InternalErr (__FILE__, __LINE__,
"unsupported data type.");
1428 r = detachfunc (gridid);
1430 ostringstream eherr;
1431 eherr <<
"Grid " << datasetname.c_str () <<
" cannot be detached.";
1432 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1438 r = closefunc (gfid);
1440 ostringstream eherr;
1442 eherr <<
"Grid " << filename.c_str () <<
" cannot be closed.";
1443 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1454 HDFEOS2ArrayGridGeoField::format_constraint (
int *offset,
int *step,
1461 Dim_iter p = dim_begin ();
1462 while (p != dim_end ()) {
1464 int start = dimension_start (p,
true);
1465 int stride = dimension_stride (p,
true);
1466 int stop = dimension_stop (p,
true);
1471 oss <<
"Array/Grid hyperslab start point "<< start <<
1472 " is greater than stop point " << stop <<
".";
1473 throw Error(malformed_expr, oss.str());
1478 count[id] = ((stop - start) / stride) + 1;
1482 "=format_constraint():"
1483 <<
"id=" <<
id <<
" offset=" << offset[
id]
1484 <<
" step=" << step[
id]
1485 <<
" count=" << count[
id]
1498 HDFEOS2ArrayGridGeoField::CalculateLatLon (int32 gridid,
int g_fieldtype,
1499 int g_specialformat,
1500 float64 * outlatlon,float64* latlon_all,
1501 int32 * offset, int32 * count,
1502 int32 * step,
int nelms,
bool write_latlon_cache)
1510 float64 lowright[2];
1512 r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
1514 ostringstream eherr;
1515 eherr <<
"cannot obtain grid information.";
1516 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1522 if (g_specialformat == 1) {
1523 upleft[0] = upleft[0] * 1000000;
1524 upleft[1] = upleft[1] * 1000000;
1525 lowright[0] = lowright[0] * 1000000;
1526 lowright[1] = lowright[1] * 1000000;
1533 if (g_specialformat == 2) {
1535 upleft[1] = 90000000.0;
1536 lowright[0] = 360000000.0;
1537 lowright[1] = -90000000.0;
1546 r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
1548 ostringstream eherr;
1549 eherr <<
"cannot obtain grid projection information";
1550 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1556 r = GDpixreginfo (gridid, &pixreg);
1558 ostringstream eherr;
1559 eherr <<
"cannot obtain grid pixel registration info.";
1560 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1567 r = GDorigininfo (gridid, &origin);
1569 ostringstream eherr;
1570 eherr <<
"cannot obtain grid origin info.";
1571 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1578 rows.resize(xdim*ydim);
1579 cols.resize(xdim*ydim);
1580 lon.resize(xdim*ydim);
1581 lat.resize(xdim*ydim);
1597 for (k = j = 0; j < ydim; ++j) {
1598 for (i = 0; i < xdim; ++i) {
1613 for (k = j = 0; j < xdim; ++j) {
1614 for (i = 0; i < ydim; ++i) {
1623 r = GDij2ll (projcode, zone, params, sphere, xdim, ydim, upleft, lowright,
1624 xdim * ydim, &rows[0], &cols[0], &lon[0], &lat[0], pixreg, origin);
1627 ostringstream eherr;
1628 eherr <<
"cannot calculate grid latitude and longitude";
1629 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1633 if(
true == write_latlon_cache) {
1634 if(GCTP_CEA == projcode || GCTP_GEO == projcode) {
1635 vector<double>temp_lat;
1636 vector<double>temp_lon;
1637 int32 temp_offset[2];
1638 int32 temp_count[2];
1646 temp_count[0] = ydim;
1648 temp_lat.resize(ydim);
1649 LatLon2DSubset(&temp_lat[0],ydim,xdim,&lat[0],temp_offset,temp_count,temp_step);
1653 temp_count[1] = xdim;
1654 temp_lon.resize(xdim);
1655 LatLon2DSubset(&temp_lon[0],ydim,xdim,&lon[0],temp_offset,temp_count,temp_step);
1657 for(i = 0; i<ydim;i++)
1658 latlon_all[i] = temp_lat[i];
1663 if(speciallon ==
true) {
1664 CorSpeLon(&temp_lon[0],xdim);
1667 for(i = 0; i<xdim;i++)
1668 latlon_all[i+ydim] = temp_lon[i];
1673 temp_count[1] = ydim;
1675 temp_lat.resize(ydim);
1676 LatLon2DSubset(&temp_lat[0],xdim,ydim,&lat[0],temp_offset,temp_count,temp_step);
1680 temp_count[0] = xdim;
1681 temp_lon.resize(xdim);
1682 LatLon2DSubset(&temp_lon[0],xdim,ydim,&lon[0],temp_offset,temp_count,temp_step);
1684 for(i = 0; i<ydim;i++)
1685 latlon_all[i] = temp_lat[i];
1690 if(speciallon ==
true)
1691 CorSpeLon(&temp_lon[0],xdim);
1693 for(i = 0; i<xdim;i++)
1694 latlon_all[i+ydim] = temp_lon[i];
1699 memcpy((
char*)(&latlon_all[0]),&lat[0],xdim*ydim*
sizeof(
double));
1700 memcpy((
char*)(&latlon_all[0])+xdim*ydim*
sizeof(
double),&lon[0],xdim*ydim*
sizeof(
double));
1707 if (nelms == (xdim * ydim)) {
1708 if (g_fieldtype == 1)
1709 memcpy (outlatlon, &lat[0], xdim * ydim *
sizeof (
double));
1711 memcpy (outlatlon, &lon[0], xdim * ydim *
sizeof (
double));
1715 if (g_fieldtype == 1)
1716 LatLon2DSubset (outlatlon, ydim, xdim, &lat[0], offset, count,
1719 LatLon2DSubset (outlatlon, ydim, xdim, &lon[0], offset, count,
1723 if (g_fieldtype == 1)
1724 LatLon2DSubset (outlatlon, xdim, ydim, &lat[0], offset, count,
1727 LatLon2DSubset (outlatlon, xdim, ydim, &lon[0], offset, count,
1735 template<
class T>
void
1736 HDFEOS2ArrayGridGeoField::LatLon2DSubset (T * outlatlon,
int majordim,
1737 int minordim, T * latlon,
1738 int32 * offset, int32 * count,
1742 T (*templatlonptr)[majordim][minordim] =
reinterpret_cast<T *[majordim][minordim]
>(latlon);
1749 int dim0count = count[0];
1750 int dim1count = count[1];
1751 int dim0index[dim0count], dim1index[dim1count];
1753 for (i = 0; i < count[0]; i++)
1754 dim0index[i] = offset[0] + i * step[0];
1757 for (j = 0; j < count[1]; j++)
1758 dim1index[j] = offset[1] + j * step[1];
1763 for (i = 0; i < count[0]; i++) {
1764 for (j = 0; j < count[1]; j++) {
1767 outlatlon[k] = (*templatlonptr)[dim0index[i]][dim1index[j]];
1769 outlatlon[k] = *(latlon + (dim0index[i] * minordim) + dim1index[j]);
1779 template <
class T >
bool HDFEOS2ArrayGridGeoField::CorLatLon (T * latlon,
1793 for (
int i = 0; i < elms; i++)
1794 if ((
int) (latlon[i]) == fv)
1801 for (
int i = 0; i < 3; i++)
1802 if ((
int) (latlon[i]) == fv)
1805 if ((
int) (latlon[elms - 1]) != fv)
1808 T increment = latlon[2] - latlon[1];
1813 index = findfirstfv (latlon, 0, elms - 1, fv);
1815 ostringstream eherr;
1816 eherr <<
"cannot calculate the fill value. ";
1817 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1820 for (
int i = index; i < elms; i++) {
1822 latlon[i] = latlon[i - 1] + increment;
1825 if (i != (elms - 1) && (g_fieldtype == 1) &&
1826 ((
float) (latlon[i]) < -90.0 || (
float) (latlon[i]) > 90.0))
1833 if (i != (elms - 1) && (g_fieldtype == 2) &&
1834 ((
float) (latlon[i]) < -180.0 || (
float) (latlon[i]) > 360.0))
1837 if (g_fieldtype == 1) {
1838 if ((
float) (latlon[elms - 1]) < -90.0)
1839 latlon[elms - 1] = (T)-90;
1840 if ((
float) (latlon[elms - 1]) > 90.0)
1841 latlon[elms - 1] = (T)90;
1844 if (g_fieldtype == 2) {
1845 if ((
float) (latlon[elms - 1]) < -180.0)
1846 latlon[elms - 1] = (T)-180.0;
1847 if ((
float) (latlon[elms - 1]) > 360.0)
1848 latlon[elms - 1] = (T)360.0;
1854 template <
class T >
void
1855 HDFEOS2ArrayGridGeoField::CorSpeLon (T * lon,
int xdim)
1858 float64 accuracy = 1e-3;
1864 for (i = 0; i < xdim; i++) {
1865 if ((
double) lon[i] < 180.0)
1866 temp = 180.0 - (double) lon[i];
1867 if ((
double) lon[i] > 180.0)
1868 temp = (double) lon[i] - 180.0;
1870 if (temp < accuracy) {
1874 else if ((
static_cast < double >(lon[i]) < 180.0)
1875 &&(
static_cast<double>(lon[i + 1]) > 180.0)) {
1883 if (speindex != -1) {
1884 for (i = speindex + 1; i < xdim; i++) {
1886 static_cast < T
> (
static_cast < double >(lon[i]) - 360.0);
1894 HDFEOS2ArrayGridGeoField::getCorrectSubset (
int *offset,
int *count,
1895 int *step, int32 * offset32,
1896 int32 * count32, int32 * step32,
1897 bool gf_condenseddim,
bool gf_ydimmajor,
1898 int gf_fieldtype,
int gf_rank)
1902 offset32[0] = (int32) offset[0];
1903 count32[0] = (int32) count[0];
1904 step32[0] = (int32) step[0];
1906 else if (gf_condenseddim) {
1910 for (
int i = 0; i < gf_rank; i++) {
1916 if (gf_ydimmajor && gf_fieldtype == 1) {
1917 offset32[0] = (int32) offset[0];
1918 count32[0] = (int32) count[0];
1919 step32[0] = (int32) step[0];
1921 else if (gf_ydimmajor && gf_fieldtype == 2) {
1922 offset32[1] = (int32) offset[0];
1923 count32[1] = (int32) count[0];
1924 step32[1] = (int32) step[0];
1926 else if (!gf_ydimmajor && gf_fieldtype == 1) {
1927 offset32[1] = (int32) offset[0];
1928 count32[1] = (int32) count[0];
1929 step32[1] = (int32) step[0];
1931 else if (!gf_ydimmajor && gf_fieldtype == 2) {
1932 offset32[0] = (int32) offset[0];
1933 count32[0] = (int32) count[0];
1934 step32[0] = (int32) step[0];
1938 throw InternalErr (__FILE__, __LINE__,
1939 "Lat/lon subset is wrong for condensed lat/lon");
1943 for (
int i = 0; i < gf_rank; i++) {
1944 offset32[i] = (int32) offset[i];
1945 count32[i] = (int32) count[i];
1946 step32[i] = (int32) step[i];
1954 template <
class T>
void
1955 HDFEOS2ArrayGridGeoField::HandleFillLatLon(vector<T> total_latlon, T* latlon,
bool gf_ydimmajor,
int gf_fieldtype, int32 xdim , int32 ydim, int32* offset, int32* count, int32* step,
int fv) {
1957 class vector <T> temp_lat;
1958 class vector <T> temp_lon;
1960 if (
true == gf_ydimmajor) {
1962 if (1 == gf_fieldtype) {
1963 temp_lat.resize(ydim);
1964 for (
int i = 0; i <(int)ydim; i++)
1965 temp_lat[i] = total_latlon[i*xdim];
1967 if (
false == CorLatLon(&temp_lat[0],gf_fieldtype,ydim,fv))
1968 throw InternalErr(__FILE__,__LINE__,
"Cannot handle the fill values in lat/lon correctly");
1970 for (
int i = 0; i <(int)(count[0]); i++)
1971 latlon[i] = temp_lat[offset[0] + i* step[0]];
1975 temp_lon.resize(xdim);
1976 for (
int i = 0; i <(int)xdim; i++)
1977 temp_lon[i] = total_latlon[i];
1980 if (
false == CorLatLon(&temp_lon[0],gf_fieldtype,xdim,fv))
1981 throw InternalErr(__FILE__,__LINE__,
"Cannot handle the fill values in lat/lon correctly");
1983 for (
int i = 0; i <(int)(count[1]); i++)
1984 latlon[i] = temp_lon[offset[1] + i* step[1]];
1990 if (1 == gf_fieldtype) {
1991 temp_lat.resize(xdim);
1992 for (
int i = 0; i <(int)xdim; i++)
1993 temp_lat[i] = total_latlon[i];
1995 if (
false == CorLatLon(&temp_lat[0],gf_fieldtype,ydim,fv))
1996 throw InternalErr(__FILE__,__LINE__,
"Cannot handle the fill values in lat/lon correctly");
1998 for (
int i = 0; i <(int)(count[1]); i++)
1999 latlon[i] = temp_lat[offset[1] + i* step[1]];
2003 temp_lon.resize(ydim);
2004 for (
int i = 0; i <(int)ydim; i++)
2005 temp_lon[i] = total_latlon[i*xdim];
2008 if (
false == CorLatLon(&temp_lon[0],gf_fieldtype,xdim,fv))
2009 throw InternalErr(__FILE__,__LINE__,
"Cannot handle the fill values in lat/lon correctly");
2011 for (
int i = 0; i <(int)(count[0]); i++)
2012 latlon[i] = temp_lon[offset[0] + i* step[0]];
2019 template <
class T >
int
2020 HDFEOS2ArrayGridGeoField::findfirstfv (T * array,
int start,
int end,
2024 if (start == end || start == (end - 1)) {
2025 if (
static_cast < int >(array[start]) == fillvalue)
2031 int current = (start + end) / 2;
2033 if (
static_cast < int >(array[current]) == fillvalue)
2034 return findfirstfv (array, start, current, fillvalue);
2036 return findfirstfv (array, current, end, fillvalue);
2048 HDFEOS2ArrayGridGeoField::CalculateSpeLatLon (int32 gridid,
int gf_fieldtype,
2049 float64 * outlatlon,
2051 int32 * count32, int32 * step32)
2059 float64 lowright[2];
2061 r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
2063 ostringstream eherr;
2064 eherr <<
"cannot obtain grid information.";
2065 throw InternalErr (__FILE__, __LINE__, eherr.str ());
2075 if(0 == xdim || 0 == ydim)
2076 throw InternalErr(__FILE__,__LINE__,
"xdim or ydim cannot be zero");
2078 if (gf_fieldtype == 1) {
2079 double latstep = 180.0 / ydim;
2081 for (
int i = 0; i < (int) (count32[0]); i++)
2082 outlatlon[i] = 90.0 -latstep/2 - latstep * (offset32[0] + i * step32[0]);
2085 double lonstep = 360.0 / xdim;
2087 for (
int i = 0; i < (int) (count32[1]); i++)
2088 outlatlon[i] = -180.0 + lonstep/2 + lonstep * (offset32[1] + i * step32[1]);
2098 HDFEOS2ArrayGridGeoField::CalculateSOMLatLon(int32 gridid,
int *start,
int *count,
int *step,
int nelms,
const string & cache_fpath,
bool write_latlon_cache)
2100 int32 projcode = -1;
2103 float64 params[NPROJ];
2106 r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
2108 throw InternalErr (__FILE__, __LINE__,
"GDprojinfo doesn't return the correct values");
2112 char dimlist[STRLEN];
2113 r = GDinqdims(gridid, dimlist, dim);
2117 throw InternalErr (__FILE__, __LINE__,
"GDinqdims doesn't return the correct values");
2119 bool is_block_180 =
false;
2120 for(
int i=0; i<MAXNDIM; i++)
2124 is_block_180 =
true;
2128 if(
false == is_block_180) {
2129 ostringstream eherr;
2130 eherr <<
"Number of Block is not " << NBLOCK ;
2131 throw InternalErr(__FILE__,__LINE__,eherr.str());
2139 r = GDgridinfo (gridid, &xdim, &ydim, ulc, lrc);
2141 throw InternalErr(__FILE__,__LINE__,
"GDgridinfo doesn't return the correct values");
2144 float32 offset[NOFFSET];
2145 char som_rw_code[]=
"r";
2146 r = GDblkSOMoffset(gridid, offset, NOFFSET, som_rw_code);
2148 throw InternalErr(__FILE__,__LINE__,
"GDblkSOMoffset doesn't return the correct values");
2150 int status = misr_init(NBLOCK, xdim, ydim, offset, ulc, lrc);
2152 throw InternalErr(__FILE__,__LINE__,
"misr_init doesn't return the correct values");
2155 int (*inv_trans[MAXPROJ+1])(double, double,
double*,
double*);
2156 inv_init((
long)projcode, (
long)zone, (
double*)params, (
long)sphere, NULL, NULL, (
int*)&iflg, inv_trans);
2158 throw InternalErr(__FILE__,__LINE__,
"inv_init doesn't return correct values");
2178 if(
true == write_latlon_cache) {
2179 vector<double>latlon_all;
2180 latlon_all.resize(xdim*ydim*NBLOCK*2);
2181 for(i =1; i <NBLOCK+1;i++)
2188 misrinv(b, l, s, &somx, &somy);
2189 sominv(somx, somy, &lon_r, &lat_r);
2190 latlon_all[npts] = lat_r*R2D;
2191 latlon_all[xdim*ydim*NBLOCK+npts] = lon_r*R2D;
2199 if(HDFCFUtil::write_vector_to_file2(cache_fpath,latlon_all,
sizeof(
double)) != (xdim*ydim*NBLOCK*2)*
sizeof(
double)) {
2200 if(remove(cache_fpath.c_str()) !=0) {
2201 throw InternalErr(__FILE__,__LINE__,
"Cannot remove the cached file.");
2206 llcache->write_cached_data(cache_fpath,xdim*ydim*NBLOCK*2*
sizeof(
double),latlon_all);
2209 vector<double>latlon;
2210 latlon.resize(nelms);
2220 for(i=0; i<count[0]; i++)
2221 for(j=0; j<count[1]; j++)
2222 for(k=0; k<count[2]; k++)
2224 if(fieldtype == 1) {
2225 latlon[npts] = latlon_all[start[0]*ydim*xdim+start[1]*ydim+start[2]+
2226 i*ydim*xdim*step[0]+j*ydim*step[1]+k*step[2]];
2229 latlon[npts] = latlon_all[xdim*ydim*NBLOCK+start[0]*ydim*xdim+start[1]*ydim+start[2]+
2230 i*ydim*xdim*step[0]+j*ydim*step[1]+k*step[2]];
2236 set_value ((dods_float64 *) &latlon[0], nelms);
2239 vector<double>latlon;
2240 latlon.resize(nelms);
2242 int e1=s1+count[0]*step[0];
2244 int e2=s2+count[1]*step[1];
2246 int e3=s3+count[2]*step[2];
2247 for(i=s1; i<e1; i+=step[0])
2248 for(j=s2; j<e2; j+=step[1])
2249 for(k=s3; k<e3; k+=step[2])
2254 misrinv(b, l, s, &somx, &somy);
2255 sominv(somx, somy, &lon_r, &lat_r);
2257 latlon[npts] = lat_r*R2D;
2259 latlon[npts] = lon_r*R2D;
2262 set_value ((dods_float64 *) &latlon[0], nelms);
2275 HDFEOS2ArrayGridGeoField::CalculateLargeGeoLatLon(int32 gridid,
int gf_fieldtype, float64* latlon, float64* latlon_all,
int *start,
int *count,
int *step,
int nelms,
bool write_latlon_cache)
2281 float64 lowright[2];
2283 r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
2285 throw InternalErr(__FILE__,__LINE__,
"GDgridinfo failed");
2288 if (0 == xdim || 0 == ydim) {
2289 throw InternalErr(__FILE__,__LINE__,
"xdim or ydim should not be zero. ");
2292 if (upleft[0]>180.0 || upleft[0] <-180.0 ||
2293 upleft[1]>90.0 || upleft[1] <-90.0 ||
2294 lowright[0] >180.0 || lowright[0] <-180.0 ||
2295 lowright[1] >90.0 || lowright[1] <-90.0) {
2297 throw InternalErr(__FILE__,__LINE__,
"lat/lon corner points are out of range. ");
2300 if (count[0] != nelms) {
2301 throw InternalErr(__FILE__,__LINE__,
"rank is not 1 ");
2303 float lat_step = (lowright[1] - upleft[1])/ydim;
2304 float lon_step = (lowright[0] - upleft[0])/xdim;
2306 if(
true == write_latlon_cache) {
2308 for(
int i = 0;i<ydim;i++)
2309 latlon_all[i] = upleft[1] + i*lat_step + lat_step/2;
2311 for(
int i = 0;i<xdim;i++)
2312 latlon_all[i+ydim] = upleft[0] + i*lon_step + lon_step/2;
2318 if (1 == gf_fieldtype) {
2319 float start_lat = upleft[1] + start[0] *lat_step + lat_step/2;
2320 float step_lat = lat_step *step[0];
2321 for (
int i = 0; i < count[0]; i++)
2322 latlon[i] = start_lat +i *step_lat;
2325 float start_lon = upleft[0] + start[0] *lon_step + lon_step/2;
2326 float step_lon = lon_step *step[0];
2327 for (
int i = 0; i < count[0]; i++)
2328 latlon[i] = start_lon +i *step_lon;
2337 HDFEOS2ArrayGridGeoField::CalculateLAMAZLatLon(int32 gridid,
int gf_fieldtype, float64* latlon, float64* latlon_all,
int *start,
int *count,
int *step,
bool write_latlon_cache)
2343 float64 lowright[2];
2345 r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
2347 throw InternalErr(__FILE__,__LINE__,
"GDgridinfo failed");
2349 vector<float64> tmp1;
2350 tmp1.resize(xdim*ydim);
2351 int32 tmp2[] = {0, 0};
2352 int32 tmp3[] = {xdim, ydim};
2353 int32 tmp4[] = {1, 1};
2355 CalculateLatLon (gridid, gf_fieldtype, specialformat, &tmp1[0], latlon_all, tmp2, tmp3, tmp4, xdim*ydim,write_latlon_cache);
2357 if(write_latlon_cache ==
true) {
2359 vector<float64> temp_lat_all;
2360 vector<float64> lat_all;
2361 temp_lat_all.resize(xdim*ydim);
2362 lat_all.resize(xdim*ydim);
2364 vector<float64> temp_lon_all;
2365 vector<float64> lon_all;
2366 temp_lon_all.resize(xdim*ydim);
2367 lon_all.resize(xdim*ydim);
2369 for(
int w=0; w < xdim*ydim; w++){
2370 temp_lat_all[w] = latlon_all[w];
2371 lat_all[w] = latlon_all[w];
2372 temp_lon_all[w] = latlon_all[w+xdim*ydim];
2373 lon_all[w] = latlon_all[w+xdim*ydim];
2378 for(
int i=0; i<ydim; i++)
2379 for(
int j=0; j<xdim; j++)
2380 if(isundef_lat(lat_all[i*xdim+j]))
2381 lat_all[i*xdim+j]=nearestNeighborLatVal(&temp_lat_all[0], i, j, ydim, xdim);
2382 for(
int i=0; i<ydim; i++)
2383 for(
int j=0; j<xdim; j++)
2384 if(isundef_lon(lon_all[i*xdim+j]))
2385 lon_all[i*xdim+j]=nearestNeighborLonVal(&temp_lon_all[0], i, j, ydim, xdim);
2388 for(
int i=0; i<xdim; i++)
2389 for(
int j=0; j<ydim; j++)
2390 if(isundef_lat(lat_all[i*ydim+j]))
2391 lat_all[i*ydim+j]=nearestNeighborLatVal(&temp_lat_all[0], i, j, xdim, ydim);
2393 for(
int i=0; i<xdim; i++)
2394 for(
int j=0; j<ydim; j++)
2395 if(isundef_lon(lon_all[i*ydim+j]))
2396 lon_all[i*ydim+j]=nearestNeighborLonVal(&temp_lon_all[0], i, j, xdim, ydim);
2400 for(
int i = 0; i<xdim*ydim;i++) {
2401 latlon_all[i] = lat_all[i];
2402 latlon_all[i+xdim*ydim] = lon_all[i];
2408 vector<float64> tmp5;
2409 tmp5.resize(xdim*ydim);
2411 for(
int w=0; w < xdim*ydim; w++)
2416 if(gf_fieldtype==1) {
2417 for(
int i=0; i<ydim; i++)
2418 for(
int j=0; j<xdim; j++)
2419 if(isundef_lat(tmp1[i*xdim+j]))
2420 tmp1[i*xdim+j]=nearestNeighborLatVal(&tmp5[0], i, j, ydim, xdim);
2421 }
else if(gf_fieldtype==2){
2422 for(
int i=0; i<ydim; i++)
2423 for(
int j=0; j<xdim; j++)
2424 if(isundef_lon(tmp1[i*xdim+j]))
2425 tmp1[i*xdim+j]=nearestNeighborLonVal(&tmp5[0], i, j, ydim, xdim);
2428 if(gf_fieldtype==1) {
2429 for(
int i=0; i<xdim; i++)
2430 for(
int j=0; j<ydim; j++)
2431 if(isundef_lat(tmp1[i*ydim+j]))
2432 tmp1[i*ydim+j]=nearestNeighborLatVal(&tmp5[0], i, j, xdim, ydim);
2433 }
else if(gf_fieldtype==2) {
2434 for(
int i=0; i<xdim; i++)
2435 for(
int j=0; j<ydim; j++)
2436 if(isundef_lon(tmp1[i*ydim+j]))
2437 tmp1[i*ydim+j]=nearestNeighborLonVal(&tmp5[0], i, j, xdim, ydim);
2441 for(
int i=start[0], k=0; i<start[0]+count[0]*step[0]; i+=step[0])
2442 for(
int j=start[1]; j<start[1]+count[1]*step[1]; j+= step[1])
2443 latlon[k++] = tmp1[i*ydim+j];