OPeNDAP Hyrax Back End Server (BES)  Updated for version 3.8.3
HDFEOS2ArrayGridGeoField.cc
Go to the documentation of this file.
1 // retrieves the latitude and longitude of the HDF-EOS2 Grid
3 // Authors: MuQun Yang <myang6@hdfgroup.org>
4 // Copyright (c) 2009-2012 The HDF Group
6 #ifdef USE_HDFEOS2_LIB
7 
8 #include "HDFEOS2ArrayGridGeoField.h"
9 #include <stdio.h>
10 #include<stdlib.h>
11 #include<sys/stat.h>
12 #include <iostream>
13 #include <sstream>
14 #include <cassert>
15 #include <debug.h>
16 #include "HDFEOS2.h"
17 #include "InternalErr.h"
18 #include <BESDebug.h>
19 #include "HDFCFUtil.h"
20 
21 #include "misrproj.h"
22 #include "errormacros.h"
23 #include <proj.h>
24 #include<sys/types.h>
25 #include<fcntl.h>
26 #include<unistd.h>
27 
28 #include "BESH4MCache.h"
29 
30 using namespace std;
31 
32 
33 #define SIGNED_BYTE_TO_INT32 1
34 
35 
36 // These two functions are used to handle MISR products with the SOM projections.
37 extern "C" {
38  int inv_init(int insys, int inzone, double *inparm, int indatum, char *fn27, char *fn83, int *iflg, int (*inv_trans[])(double, double, double*, double*));
39  int sominv(double y, double x, double *lon, double *lat);
40 }
41 
42 bool
43 HDFEOS2ArrayGridGeoField::read ()
44 {
45 
46  BESDEBUG("h4","Coming to HDFEOS2ArrayGridGeoField read "<<endl);
47 
48  string check_pass_fileid_key_str="H4.EnablePassFileID";
49  bool check_pass_fileid_key = false;
50  check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
51 
52 
53  // Currently The latitude and longitude rank from HDF-EOS2 grid must be either 1-D or 2-D.
54  // However, For SOM projection the final rank will become 3.
55  if (rank < 1 || rank > 2) {
56  throw InternalErr (__FILE__, __LINE__, "The rank of geo field is greater than 2, currently we don't support 3-D lat/lon cases.");
57  }
58 
59  // MISR SOM file's final rank is 3. So declare a new variable.
60  int final_rank = -1;
61 
62  if (true == condenseddim)
63  final_rank = 1;
64  else if(4 == specialformat)// For the SOM projection, the final output of latitude/longitude rank should be 3.
65  final_rank = 3;
66  else
67  final_rank = rank;
68 
69  vector<int> offset;
70  offset.resize(final_rank);
71  vector<int> count;
72  count.resize(final_rank);
73  vector<int> step;
74  step.resize(final_rank);
75 
76  int nelms = -1;
77 
78  // Obtain the number of the subsetted elements
79  nelms = format_constraint (&offset[0], &step[0], &count[0]);
80 
81  // Define function pointers to handle both grid and swath Note: in
82  // this code, we only handle grid, implementing this way is to
83  // keep the same style as the read functions in other files.
84  int32 (*openfunc) (char *, intn);
85  int32 (*attachfunc) (int32, char *);
86  intn (*detachfunc) (int32);
87  intn (*fieldinfofunc) (int32, char *, int32 *, int32 *, int32 *, char *);
88  intn (*readfieldfunc) (int32, char *, int32 *, int32 *, int32 *, void *);
89 
90  string datasetname;
91  openfunc = GDopen;
92  attachfunc = GDattach;
93  detachfunc = GDdetach;
94  fieldinfofunc = GDfieldinfo;
95  readfieldfunc = GDreadfield;
96  datasetname = gridname;
97 
98  int32 gfid = -1;
99  int32 gridid = -1;
100 
101  /* Declare projection code, zone, etc grid parameters. */
102  int32 projcode = -1;
103  int32 zone = -1;
104  int32 sphere = -1;
105  float64 params[16];
106 
107  int32 xdim = 0;
108  int32 ydim = 0;
109 
110  float64 upleft[2];
111  float64 lowright[2];
112 
113  string cache_fpath="";
114  bool use_cache = false;
115 
116  // Check if passing file IDs to data
117  if(true == check_pass_fileid_key)
118  gfid = gridfd;
119  else {
120  gfid = openfunc (const_cast < char *>(filename.c_str ()), DFACC_READ);
121  if (gfid < 0) {
122  ostringstream eherr;
123  eherr << "File " << filename.c_str () << " cannot be open.";
124  throw InternalErr (__FILE__, __LINE__, eherr.str ());
125  }
126  }
127 
128  // Attach the grid id; make the grid valid.
129  gridid = attachfunc (gfid, const_cast < char *>(datasetname.c_str ()));
130  if (gridid < 0) {
131  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
132  ostringstream eherr;
133  eherr << "Grid " << datasetname.c_str () << " cannot be attached.";
134  throw InternalErr (__FILE__, __LINE__, eherr.str ());
135  }
136 
137  if(false == llflag) {
138 
139  // Cache
140  // Check if a BES key H4.EnableEOSGeoCacheFile is true, if yes, we will check
141  // if a lat/lon cache file exists and if we can read lat/lon from this file.
142  string check_eos_geo_cache_key = "H4.EnableEOSGeoCacheFile";
143  bool enable_eos_geo_cache_key = false;
144  enable_eos_geo_cache_key = HDFCFUtil::check_beskeys(check_eos_geo_cache_key);
145 
146  if(true == enable_eos_geo_cache_key) {
147 
148  use_cache = true;
149  BESH4Cache *llcache = BESH4Cache::get_instance();
150 
151  // Here we have a sanity check for the cached parameters:Cached directory,file prefix and cached directory size.
152  // Supposedly Hyrax BES cache feature should check this and the code exists. However, the
153  // current hyrax 1.9.7 doesn't provide this feature. KY 2014-10-24
154  string bescachedir= llcache->getCacheDirFromConfig();
155  string bescacheprefix = llcache->getCachePrefixFromConfig();
156  unsigned int cachesize = llcache->getCacheSizeFromConfig();
157 
158  if(("" == bescachedir)||(""==bescacheprefix)||(cachesize <=0)){
159  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
160  throw InternalErr (__FILE__, __LINE__, "Either the cached dir is empty or the prefix is NULL or the cache size is not set.");
161  }
162  else {
163  struct stat sb;
164  if(stat(bescachedir.c_str(),&sb) !=0) {
165  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
166  string err_mesg="The cached directory " + bescachedir;
167  err_mesg = err_mesg + " doesn't exist. ";
168  throw InternalErr(__FILE__,__LINE__,err_mesg);
169 
170  }
171  else {
172  if(true == S_ISDIR(sb.st_mode)) {
173  if(access(bescachedir.c_str(),R_OK|W_OK|X_OK) == -1) {
174  //if(access(bescachedir.c_str(),R_OK) == -1)
175  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
176  string err_mesg="The cached directory " + bescachedir;
177  err_mesg = err_mesg + " can NOT be read,written or executable.";
178  throw InternalErr(__FILE__,__LINE__,err_mesg);
179  }
180 
181  }
182  else {
183  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
184  string err_mesg="The cached directory " + bescachedir;
185  err_mesg = err_mesg + " is not a directory.";
186  throw InternalErr(__FILE__,__LINE__,err_mesg);
187 
188  }
189  }
190  }
191 
192  string cache_fname=llcache->getCachePrefixFromConfig();
193 
194  intn r = -1;
195  r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
196  if (r!=0) {
197  detachfunc(gridid);
198  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
199  throw InternalErr (__FILE__, __LINE__, "GDprojinfo failed");
200  }
201 
202  // Retrieve dimensions and X-Y coordinates of corners
203  if (GDgridinfo(gridid, &xdim, &ydim, upleft,
204  lowright) == -1) {
205  detachfunc(gridid);
206  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
207  throw InternalErr (__FILE__, __LINE__, "GDgridinfo failed");
208  }
209 
210  // Retrieve pixel registration information
211  int32 pixreg = 0;
212  r = GDpixreginfo (gridid, &pixreg);
213  if (r != 0) {
214  detachfunc(gridid);
215  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
216  ostringstream eherr;
217  eherr << "cannot obtain grid pixel registration info.";
218  throw InternalErr (__FILE__, __LINE__, eherr.str ());
219  }
220 
221  //Retrieve grid pixel origin
222  int32 origin = 0;
223  r = GDorigininfo (gridid, &origin);
224  if (r != 0) {
225  detachfunc(gridid);
226  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
227  ostringstream eherr;
228  eherr << "cannot obtain grid origin info.";
229  throw InternalErr (__FILE__, __LINE__, eherr.str ());
230  }
231 
232 
233  // Projection code,zone,sphere,pix,origin
234  cache_fname +=HDFCFUtil::get_int_str(projcode);
235  cache_fname +=HDFCFUtil::get_int_str(zone);
236  cache_fname +=HDFCFUtil::get_int_str(sphere);
237  cache_fname +=HDFCFUtil::get_int_str(pixreg);
238  cache_fname +=HDFCFUtil::get_int_str(origin);
239 
240 
241  // Xdimsize and ydimsize. Although it is rare, need to consider dim major.
242  // Whether latlon is ydim,xdim or xdim,ydim.
243  if(ydimmajor) {
244  cache_fname +=HDFCFUtil::get_int_str(ydim);
245  cache_fname +=HDFCFUtil::get_int_str(xdim);
246 
247  }
248  else {
249  cache_fname +=HDFCFUtil::get_int_str(xdim);
250  cache_fname +=HDFCFUtil::get_int_str(ydim);
251  }
252 
253  // upleft,lowright
254  // HDF-EOS upleft,lowright,params use DDDDMMMSSS.6 digits. So choose %17.6f.
255  cache_fname +=HDFCFUtil::get_double_str(upleft[0],17,6);
256  cache_fname +=HDFCFUtil::get_double_str(upleft[1],17,6);
257  cache_fname +=HDFCFUtil::get_double_str(lowright[0],17,6);
258  cache_fname +=HDFCFUtil::get_double_str(lowright[1],17,6);
259 
260  // According to HDF-EOS2 document, only 13 parameters are used.
261  for(int ipar = 0; ipar<13;ipar++) {
262 //cerr<<"params["<<ipar<<"] is "<<params[ipar]<<endl;
263  cache_fname+=HDFCFUtil::get_double_str(params[ipar],17,6);
264  }
265 
266  cache_fpath = bescachedir + "/"+ cache_fname;
267 #if 0
268 cerr<<"cache file path is "<<cache_fpath <<endl;
269 cerr<<"obtain file path from BESMCache "<<endl;
270 cerr<<"Name is "<<llcache->get_cache_file_name_h4(cache_fpath,false) <<endl;
271 int fd;
272 llcache->get_read_lock(cache_fpath,fd);
273 cerr<<"after testing get_read_lock"<<endl;
274 #endif
275 
276  try {
277  do {
278  int expected_file_size = 0;
279  if(GCTP_CEA == projcode || GCTP_GEO == projcode)
280  expected_file_size = (xdim+ydim)*sizeof(double);
281  else if(GCTP_SOM == projcode)
282  expected_file_size = xdim*ydim*NBLOCK*2*sizeof(double);
283  else
284  expected_file_size = xdim*ydim*2*sizeof(double);
285 
286  int fd = 0;
287  bool latlon_from_cache = llcache->get_data_from_cache(cache_fpath, expected_file_size,fd);
288 #if 0
289 if(true == latlon_from_cache)
290  cerr<<"the cached file exists: "<<endl;
291 else
292  cerr<<"the cached file doesn't exist "<< endl;
293 #endif
294  if(false == latlon_from_cache)
295  break;
296 
297  // Get the offset of lat/lon in the cached file. Since lat is stored first and lon is stored second,
298  // so offset_1d for lat/lon is different.
299  // We still need to consider different projections. 1D,2D,3D reading.Need also to consider dim major and special format.
300  size_t offset_1d = 0;
301 
302  // Get the count of the lat/lon from the cached file.
303  // Notice the data is read continuously. So starting from the offset point, we have to read all elements until the
304  // last points. The total count for the data point is bigger than the production of count and step.
305  int count_1d = 1;
306 
307  if(GCTP_CEA == projcode|| GCTP_GEO== projcode) {
308 
309  // It seems that for 1-D lat/lon, regardless of xdimmajor or ydimmajor. It is always Lat[YDim],Lon[XDim}, check getCorrectSubset
310  // So we don't need to consider the dimension major case.
311  offset_1d = (fieldtype == 1) ?offset[0] :(ydim+offset[0]);
312 #if 0
313  if(true == ydimmajor) {
314  offset_1d = (fieldtype == 1) ?offset[0] :(ydim*sizeof(double)+offset[0]);
315  }
316  else {
317  offset_1d = (fieldtype == 1) ?offset[0] :(xdim*sizeof(double)+offset[0]);
318  }
319 #endif
320  count_1d = 1+(count[0]-1)*step[0];
321  }
322  else if (GCTP_SOM == projcode) {
323 
324  if(true == ydimmajor) {
325  offset_1d = (fieldtype == 1)?(offset[0]*xdim*ydim+offset[1]*xdim+offset[2])
326  :(offset[0]*xdim*ydim+offset[1]*xdim+offset[2]+expected_file_size/2/sizeof(double));
327  }
328  else {
329  offset_1d = (fieldtype == 1)?(offset[0]*xdim*ydim+offset[1]*ydim+offset[2])
330  :(offset[0]*xdim*ydim+offset[1]*ydim+offset[2]+expected_file_size/2/sizeof(double));
331  }
332 
333  int total_count_dim0 = (count[0]-1)*step[0];
334  int total_count_dim1 = (count[1]-1)*step[1];
335  int total_count_dim2 = (count[2]-1)*step[2];
336  int total_dim1 = (true ==ydimmajor)?ydim:xdim;
337  int total_dim2 = (true ==ydimmajor)?xdim:ydim;
338 
339  // Flatten the 3-D index to 1-D
340  // This calculation can be generalized from nD to 1D
341  // but since we only use it here. Just keep it this way.
342  count_1d = 1 + total_count_dim0*total_dim1*total_dim2 + total_count_dim1*total_dim2 + total_count_dim2;
343 
344  }
345  else {// 2-D lat/lon case
346  if (true == ydimmajor)
347  offset_1d = (fieldtype == 1) ?(offset[0] * xdim + offset[1]):(expected_file_size/2/sizeof(double)+offset[0]*xdim+offset[1]);
348  else
349  offset_1d = (fieldtype == 1) ?(offset[0] * ydim + offset[1]):(expected_file_size/2/sizeof(double)+offset[0]*ydim+offset[1]);
350 
351  // Flatten the 2-D index to 1-D
352  int total_count_dim0 = (count[0]-1)*step[0];
353  int total_count_dim1 = (count[1]-1)*step[1];
354  int total_dim1 = (true ==ydimmajor)?xdim:ydim;
355 
356  count_1d = 1 + total_count_dim0*total_dim1 + total_count_dim1;
357  }
358 
359  // Assign a vector to store lat/lon
360  vector<double> latlon_1d;
361  latlon_1d.resize(count_1d);
362 
363  // Read lat/lon from the file.
364  //int fd;
365  //fd = open(cache_fpath.c_str(),O_RDONLY,0666);
366  // TODO: Use BESLog to report that the cached file cannot be read.
367  off_t fpos = lseek(fd,sizeof(double)*offset_1d,SEEK_SET);
368  if (-1 == fpos) {
369  llcache->unlock_and_close(cache_fpath);
370  llcache->purge_file(cache_fpath);
371  break;
372  }
373  ssize_t read_size = HDFCFUtil::read_vector_from_file(fd,latlon_1d,sizeof(double));
374  llcache->unlock_and_close(cache_fpath);
375  if((-1 == read_size) || ((size_t)read_size != count_1d*sizeof(double))) {
376  llcache->purge_file(cache_fpath);
377  break;
378  }
379 
380 // Leave the debugging comments for the time being.
381 #if 0
382  // ONLY READ the subset
383  FILE *pFile;
384  pFile = fopen(cache_fpath.c_str(),"rb");
385  if(NULL == pFile)
386  break;
387 
388  int ret_value = fseek(pFile,sizeof(double)*offset_1d,SEEK_SET);
389  if(ret_value != 0) {
390  // fall back to the original calculation.
391  fclose(pFile);
392  break;
393  }
394 cerr<<"field name is "<<fieldname <<endl;
395 cerr<<"fread is checked "<<endl;
396 cerr<<"offset_1d is "<<offset_1d <<endl;
397 cerr<<"count_1d is "<<count_1d <<endl;
398 
399 
400  ret_value = fread(&latlon_1d[0],sizeof(double),count_1d,pFile);
401  if(0 == ret_value) {
402  // fall back to the original calculation
403  cerr<<"fread fails "<<endl;
404  fclose(pFile);
405  break;
406  }
407 #endif
408 #if 0
409 for(int i =0;i<count_1d;i++)
410 cerr<<"latlon_1d["<<i<<"]"<<latlon_1d[i]<<endl;
411 #endif
412 
413  int total_count = 1;
414  for (int i_rank = 0; i_rank<final_rank;i_rank++)
415  total_count = total_count*count[i_rank];
416 
417  // We will see if there is a shortcut that the lat/lon is accessed with
418  // one-big-block. Actually this is the most common case. If we find
419  // such a case, we simply read the whole data into the latlon buffer and
420  // send it to BES.
421  if(total_count == count_1d) {
422  set_value((dods_float64*)&latlon_1d[0],nelms);
423  detachfunc(gridid);
424  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
425  return false;
426  }
427 
428  vector<double>latlon;
429  latlon.resize(total_count);
430 
431  // Retrieve latlon according to the projection
432  if(GCTP_CEA == projcode|| GCTP_GEO== projcode) {
433  for (int i = 0; i <total_count;i++)
434  latlon[i] = latlon_1d[i*step[0]];
435 
436  }
437  else if (GCTP_SOM == projcode) {
438 
439  for (int i =0; i<count[0];i++)
440  for(int j =0;j<count[1];j++)
441  for(int k=0;k<count[2];k++)
442  latlon[i*count[1]*count[2]+j*count[2]+k]=(true == ydimmajor)
443  ?latlon_1d[i*ydim*xdim*step[0]+j*xdim*step[1]+k*step[2]]
444  :latlon_1d[i*ydim*xdim*step[0]+j*ydim*step[1]+k*step[2]];
445  }
446  else {
447  for (int i =0; i<count[0];i++)
448  for(int j =0;j<count[1];j++)
449  latlon[i*count[1]+j]=(true == ydimmajor)
450  ?latlon_1d[i*xdim*step[0]+j*step[1]]
451  :latlon_1d[i*ydim*step[0]+j*step[1]];
452 
453  }
454 
455  set_value((dods_float64*)&latlon[0],nelms);
456  detachfunc(gridid);
457  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
458  return false;
459 
460  } while (0);
461 
462  }
463  catch(...) {
464  detachfunc(gridid);
465  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
466  throw;
467  }
468  //detachfunc(gridid);
469 
470  }
471 
472  }
473 
474 
475  // In this step, if use_cache is true, we always need to write the lat/lon into the cache.
476  // SOM projection should be calculated differently. If turning on the lat/lon cache feature, it also needs to be handled differently.
477  if(specialformat == 4) {// SOM projection
478  try {
479  CalculateSOMLatLon(gridid, &offset[0], &count[0], &step[0], nelms,cache_fpath,use_cache);
480  }
481  catch(...) {
482  detachfunc(gridid);
483  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
484  throw;
485  }
486  detachfunc(gridid);
487  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
488  return false;
489  }
490 
491  // We define offset,count and step in int32 datatype.
492  vector<int32>offset32;
493  offset32.resize(rank);
494 
495  vector<int32>count32;
496  count32.resize(rank);
497 
498  vector<int32>step32;
499  step32.resize(rank);
500 
501 
502  // Obtain offset32 with the correct rank, the rank of lat/lon of
503  // GEO and CEA projections in the file may be 2 instead of 1.
504  try {
505  getCorrectSubset (&offset[0], &count[0], &step[0], &offset32[0], &count32[0], &step32[0], condenseddim, ydimmajor, fieldtype, rank);
506  }
507  catch(...) {
508  detachfunc(gridid);
509  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
510  throw;
511  }
512 
513  // The following case handles when the lat/lon is not provided.
514  if (llflag == false) { // We have to calculate the lat/lon
515 
516  vector<float64>latlon;
517  latlon.resize(nelms);
518 
519  // If projection code etc. is not retrieved, retrieve them.
520  // When specialformat is 3, the file is a file of which the project code is set to -1, we need to skip it. KY 2014-09-11
521  if(projcode == -1 && specialformat !=3) {
522 
523 
524  intn r = 0;
525  r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
526  if (r!=0) {
527  detachfunc(gridid);
528  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
529  throw InternalErr (__FILE__, __LINE__, "GDprojinfo failed");
530  }
531 
532  // Retrieve dimensions and X-Y coordinates of corners
533  if (GDgridinfo(gridid, &xdim, &ydim, upleft,
534  lowright) == -1) {
535  detachfunc(gridid);
536  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
537  throw InternalErr (__FILE__, __LINE__, "GDgridinfo failed");
538  }
539  }
540 
541 
542  // Handle LAMAZ projection first.
543  if (GCTP_LAMAZ == projcode) {
544  try {
545  vector<double>latlon_all;
546  latlon_all.resize(xdim*ydim*2);
547 
548  CalculateLAMAZLatLon(gridid, fieldtype, &latlon[0], &latlon_all[0],&offset[0], &count[0], &step[0], nelms,use_cache);
549  if(true == use_cache) {
550 
551  BESH4Cache *llcache = BESH4Cache::get_instance();
552  llcache->write_cached_data(cache_fpath,xdim*ydim*2*sizeof(double),latlon_all);
553 
554  }
555  }
556  catch(...) {
557  detachfunc(gridid);
558  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
559  throw;
560  }
561  set_value ((dods_float64 *) &latlon[0], nelms);
562  detachfunc(gridid);
563  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
564  return false;
565  }
566 
567  // Aim to handle large MCD Grid such as 21600*43200 lat,lon
568  if (specialformat == 1) {
569 
570  try {
571  vector<double>latlon_all;
572  latlon_all.resize(xdim+ydim);
573 
574  CalculateLargeGeoLatLon(gridid, fieldtype,&latlon[0], &latlon_all[0],&offset[0], &count[0], &step[0], nelms,use_cache);
575  if(true == use_cache) {
576 
577  BESH4Cache *llcache = BESH4Cache::get_instance();
578  llcache->write_cached_data(cache_fpath,(xdim+ydim)*sizeof(double),latlon_all);
579 
580 // if(HDFCFUtil::write_vector_to_file(cache_fpath,latlon_all,sizeof(double)) != ((xdim+ydim))) {
581 #if 0
582  if(HDFCFUtil::write_vector_to_file2(cache_fpath,latlon_all,sizeof(double)) != ((xdim+ydim)*sizeof(double))) {
583  if(remove(cache_fpath.c_str()) !=0) {
584  throw InternalErr(__FILE__,__LINE__,"Cannot remove the cached file.");
585  }
586  }
587 #endif
588  }
589 
590  }
591  catch(...) {
592  detachfunc(gridid);
593  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
594  throw;
595  }
596  set_value((dods_float64 *)&latlon[0],nelms);
597  detachfunc(gridid);
598  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
599 
600  return false;
601  }
602 
603  // Now handle other cases,note the values will be written after the if-block
604  else if (specialformat == 3) {// Have to provide latitude and longitude by ourselves
605  try {
606  CalculateSpeLatLon (gridid, fieldtype, &latlon[0], &offset32[0], &count32[0], &step32[0], nelms);
607  }
608  catch(...) {
609  detachfunc(gridid);
610  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
611  throw;
612 
613  }
614  detachfunc(gridid);
615  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
616  }
617  else {// This is mostly general case, it will calculate lat/lon with GDij2ll.
618 
619  // Cache: check the flag and decide whether to calculate the lat/lon.
620  vector<double>latlon_all;
621 
622  if(GCTP_GEO == projcode || GCTP_CEA == projcode)
623  latlon_all.resize(xdim+ydim);
624  else
625  latlon_all.resize(xdim*ydim*2);
626 
627  CalculateLatLon (gridid, fieldtype, specialformat, &latlon[0],&latlon_all[0],
628  &offset32[0], &count32[0], &step32[0], nelms,use_cache);
629  if(true == use_cache) {
630 // size_t num_item_to_write = HDFCFUtil::write_vector_to_file2(cache_fpath,latlon_all,sizeof(double));
631 
632  size_t num_item_expected = 0;
633  if(GCTP_GEO == projcode || GCTP_CEA == projcode)
634  num_item_expected = xdim + ydim;
635  else
636  num_item_expected = xdim*ydim*2;
637 
638  BESH4Cache *llcache = BESH4Cache::get_instance();
639  llcache->write_cached_data(cache_fpath,num_item_expected*sizeof(double),latlon_all);
640 
641  }
642  detachfunc(gridid);
643  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
644  }
645 
646  // Some longitude values need to be corrected.
647  if (speciallon && fieldtype == 2) {
648  CorSpeLon (&latlon[0], nelms);
649  }
650 
651  set_value ((dods_float64 *) &latlon[0], nelms);
652 
653  return false;
654  }
655 
656 
657  // Now lat and lon are stored as HDF-EOS2 fields. We need to read the lat and lon values from the fields.
658  int32 tmp_rank = -1;
659  int32 tmp_dims[rank];
660  char tmp_dimlist[1024];
661  int32 type = -1;
662  intn r = -1;
663 
664  // Obtain field info.
665  r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
666  &tmp_rank, tmp_dims, &type, tmp_dimlist);
667 
668  if (r != 0) {
669  detachfunc(gridid);
670  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
671  ostringstream eherr;
672  eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
673  throw InternalErr (__FILE__, __LINE__, eherr.str ());
674  }
675 
676  // Retrieve dimensions and X-Y coordinates of corners
677  //int32 xdim = 0;
678  //int32 ydim = 0;
679 
680  //float64 upleft[2];
681  //float64 lowright[2];
682 
683  r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
684  if (r != 0) {
685  detachfunc(gridid);
686  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
687  ostringstream eherr;
688  eherr << "Grid " << datasetname.c_str () << " information cannot be obtained.";
689  throw InternalErr (__FILE__, __LINE__, eherr.str ());
690  }
691 
692  // Retrieve all GCTP projection information
693  //int32 projcode = -1;
694  //int32 zone = -1;
695  //int32 sphere = -1;
696  //float64 params[16];
697 
698  r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
699  if (r != 0) {
700  detachfunc(gridid);
701  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
702  ostringstream eherr;
703  eherr << "Grid " << datasetname.c_str () << " projection info. cannot be obtained.";
704  throw InternalErr (__FILE__, __LINE__, eherr.str ());
705  }
706 
707  if (projcode != GCTP_GEO) { // Just retrieve the data like other fields
708  // We have to loop through all datatype and read the lat/lon out.
709  switch (type) {
710  case DFNT_INT8:
711  {
712  vector<int8> val;
713  val.resize(nelms);
714  r = readfieldfunc (gridid,
715  const_cast < char *>(fieldname.c_str ()),
716  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
717  if (r != 0) {
718  detachfunc(gridid);
719  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
720  ostringstream eherr;
721  eherr << "field " << fieldname.c_str () << "cannot be read.";
722  throw InternalErr (__FILE__, __LINE__, eherr.str ());
723  }
724 
725  // DAP2 requires the map of SIGNED_BYTE to INT32 if
726  // SIGNED_BYTE_TO_INT32 is defined.
727 #ifndef SIGNED_BYTE_TO_INT32
728  set_value ((dods_byte *) &val[0], nelms);
729 #else
730  vector<int32>newval;
731  newval.resize(nelms);
732 
733  for (int counter = 0; counter < nelms; counter++)
734  newval[counter] = (int32) (val[counter]);
735 
736  set_value ((dods_int32 *) &newval[0], nelms);
737 #endif
738 
739  }
740  break;
741  case DFNT_UINT8:
742  case DFNT_UCHAR8:
743 
744  {
745  vector<uint8> val;
746  val.resize(nelms);
747  r = readfieldfunc (gridid,
748  const_cast < char *>(fieldname.c_str ()),
749  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
750  if (r != 0) {
751  detachfunc(gridid);
752  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
753  ostringstream eherr;
754  eherr << "field " << fieldname.c_str () << "cannot be read.";
755  throw InternalErr (__FILE__, __LINE__, eherr.str ());
756  }
757  set_value ((dods_byte *) &val[0], nelms);
758 
759  }
760  break;
761 
762  case DFNT_INT16:
763 
764  {
765  vector<int16> val;
766  val.resize(nelms);
767 
768  r = readfieldfunc (gridid,
769  const_cast < char *>(fieldname.c_str ()),
770  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
771  if (r != 0) {
772  detachfunc(gridid);
773  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
774  ostringstream eherr;
775  eherr << "field " << fieldname.c_str () << "cannot be read.";
776  throw InternalErr (__FILE__, __LINE__, eherr.str ());
777  }
778 
779  set_value ((dods_int16 *) &val[0], nelms);
780 
781  }
782  break;
783  case DFNT_UINT16:
784 
785  {
786  vector<uint16> val;
787  val.resize(nelms);
788 
789  r = readfieldfunc (gridid,
790  const_cast < char *>(fieldname.c_str ()),
791  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
792  if (r != 0) {
793  detachfunc(gridid);
794  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
795  ostringstream eherr;
796  eherr << "field " << fieldname.c_str () << "cannot be read.";
797  throw InternalErr (__FILE__, __LINE__, eherr.str ());
798  }
799 
800  set_value ((dods_uint16 *) &val[0], nelms);
801  }
802  break;
803  case DFNT_INT32:
804 
805  {
806  vector<int32> val;
807  val.resize(nelms);
808 
809  r = readfieldfunc (gridid,
810  const_cast < char *>(fieldname.c_str ()),
811  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
812  if (r != 0) {
813  detachfunc(gridid);
814  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
815  ostringstream eherr;
816  eherr << "field " << fieldname.c_str () << "cannot be read.";
817  throw InternalErr (__FILE__, __LINE__, eherr.str ());
818  }
819 
820  set_value ((dods_int32 *) &val[0], nelms);
821  }
822  break;
823  case DFNT_UINT32:
824 
825  {
826  vector<uint32> val;
827  val.resize(nelms);
828 
829  r = readfieldfunc (gridid,
830  const_cast < char *>(fieldname.c_str ()),
831  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
832  if (r != 0) {
833  detachfunc(gridid);
834  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
835  ostringstream eherr;
836  eherr << "field " << fieldname.c_str () << "cannot be read.";
837  throw InternalErr (__FILE__, __LINE__, eherr.str ());
838  }
839  set_value ((dods_uint32 *) &val[0], nelms);
840  }
841  break;
842  case DFNT_FLOAT32:
843 
844  {
845  vector<float32> val;
846  val.resize(nelms);
847 
848  r = readfieldfunc (gridid,
849  const_cast < char *>(fieldname.c_str ()),
850  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
851  if (r != 0) {
852  detachfunc(gridid);
853  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
854  ostringstream eherr;
855  eherr << "field " << fieldname.c_str () << "cannot be read.";
856  throw InternalErr (__FILE__, __LINE__, eherr.str ());
857  }
858 
859  set_value ((dods_float32 *) &val[0], nelms);
860  }
861  break;
862  case DFNT_FLOAT64:
863 
864  {
865  vector<float64> val;
866  val.resize(nelms);
867 
868  r = readfieldfunc (gridid,
869  const_cast < char *>(fieldname.c_str ()),
870  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
871  if (r != 0) {
872  detachfunc(gridid);
873  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
874  ostringstream eherr;
875  eherr << "field " << fieldname.c_str () << "cannot be read.";
876  throw InternalErr (__FILE__, __LINE__, eherr.str ());
877  }
878 
879  set_value ((dods_float64 *) &val[0], nelms);
880  }
881  break;
882  default:
883  {
884  detachfunc(gridid);
885  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
886  InternalErr (__FILE__, __LINE__, "unsupported data type.");
887  }
888 
889  }
890  }
891  else {// Only handle special cases for the Geographic Projection
892  // We find that lat/lon of the geographic projection in some
893  // files include fill values. So we recalculate lat/lon based
894  // on starting value,step values and number of steps.
895  // GDgetfillvalue will return 0 if having fill values.
896  // The other returned value indicates no fillvalue is found inside the lat or lon.
897  switch (type) {
898  case DFNT_INT8:
899  {
900  vector<int8> val;
901  val.resize(nelms);
902 
903  int8 fillvalue = 0;
904 
905  r = GDgetfillvalue (gridid,
906  const_cast < char *>(fieldname.c_str ()),
907  &fillvalue);
908  if (r == 0) {
909  int ifillvalue = fillvalue;
910 
911  vector <int8> temp_total_val;
912  //The previous size doesn't make sense since num_elems = xdim*ydim
913  temp_total_val.resize(xdim*ydim);
914  //temp_total_val.resize(xdim*ydim*4);
915 
916  r = readfieldfunc(gridid,
917  const_cast < char *>(fieldname.c_str ()),
918  NULL, NULL, NULL, (void *)(&temp_total_val[0]));
919 
920  if (r != 0) {
921  detachfunc(gridid);
922  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
923  ostringstream eherr;
924  eherr << "field " << fieldname.c_str () << "cannot be read.";
925  throw InternalErr (__FILE__, __LINE__, eherr.str ());
926  }
927 
928  try {
929  // Recalculate lat/lon for the geographic projection lat/lon that has fill values
930  HandleFillLatLon(temp_total_val, (int8*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
931  }
932  catch(...) {
933  detachfunc(gridid);
934  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
935  throw;
936  }
937 
938  }
939 
940  else {
941 
942  r = readfieldfunc (gridid,
943  const_cast < char *>(fieldname.c_str ()),
944  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
945  if (r != 0) {
946  detachfunc(gridid);
947  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
948  ostringstream eherr;
949  eherr << "field " << fieldname.c_str () << "cannot be read.";
950  throw InternalErr (__FILE__, __LINE__, eherr.str ());
951  }
952  }
953 
954  if (speciallon && fieldtype == 2)
955  CorSpeLon ((int8 *) &val[0], nelms);
956 
957 
958 #ifndef SIGNED_BYTE_TO_INT32
959  set_value ((dods_byte *) &val[0], nelms);
960 #else
961  vector<int32>newval;
962  newval.resize(nelms);
963 
964  for (int counter = 0; counter < nelms; counter++)
965  newval[counter] = (int32) (val[counter]);
966 
967  set_value ((dods_int32 *) &newval[0], nelms);
968 
969 #endif
970 
971  }
972  break;
973 
974  case DFNT_UINT8:
975  case DFNT_UCHAR8:
976  {
977  vector<uint8> val;
978  val.resize(nelms);
979 
980  uint8 fillvalue = 0;
981 
982  r = GDgetfillvalue (gridid,
983  const_cast < char *>(fieldname.c_str ()),
984  &fillvalue);
985 
986  if (r == 0) {
987 
988  int ifillvalue = fillvalue;
989  vector <uint8> temp_total_val;
990  temp_total_val.resize(xdim*ydim);
991 
992  r = readfieldfunc(gridid,
993  const_cast < char *>(fieldname.c_str ()),
994  NULL, NULL, NULL, (void *)(&temp_total_val[0]));
995 
996  if (r != 0) {
997  detachfunc(gridid);
998  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
999  ostringstream eherr;
1000  eherr << "field " << fieldname.c_str () << "cannot be read.";
1001  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1002  }
1003 
1004  try {
1005  HandleFillLatLon(temp_total_val, (uint8*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
1006  }
1007  catch(...) {
1008  detachfunc(gridid);
1009  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1010  throw;
1011  }
1012 
1013  }
1014 
1015  else {
1016 
1017  r = readfieldfunc (gridid,
1018  const_cast < char *>(fieldname.c_str ()),
1019  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
1020  if (r != 0) {
1021  detachfunc(gridid);
1022  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1023  ostringstream eherr;
1024  eherr << "field " << fieldname.c_str () << "cannot be read.";
1025  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1026  }
1027  }
1028 
1029  if (speciallon && fieldtype == 2)
1030  CorSpeLon ((uint8 *) &val[0], nelms);
1031  set_value ((dods_byte *) &val[0], nelms);
1032 
1033  }
1034  break;
1035 
1036  case DFNT_INT16:
1037  {
1038  vector<int16> val;
1039  val.resize(nelms);
1040 
1041  int16 fillvalue = 0;
1042 
1043  r = GDgetfillvalue (gridid,
1044  const_cast < char *>(fieldname.c_str ()),
1045  &fillvalue);
1046  if (r == 0) {
1047 
1048  int ifillvalue = fillvalue;
1049  vector <int16> temp_total_val;
1050  temp_total_val.resize(xdim*ydim);
1051 
1052  r = readfieldfunc(gridid,
1053  const_cast < char *>(fieldname.c_str ()),
1054  NULL, NULL, NULL, (void *)(&temp_total_val[0]));
1055 
1056  if (r != 0) {
1057  detachfunc(gridid);
1058  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1059  ostringstream eherr;
1060  eherr << "field " << fieldname.c_str () << "cannot be read.";
1061  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1062  }
1063 
1064  try {
1065  HandleFillLatLon(temp_total_val, (int16*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
1066  }
1067  catch(...) {
1068  detachfunc(gridid);
1069  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1070  throw;
1071  }
1072 
1073  }
1074 
1075  else {
1076 
1077  r = readfieldfunc (gridid,
1078  const_cast < char *>(fieldname.c_str ()),
1079  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
1080  if (r != 0) {
1081  detachfunc(gridid);
1082  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1083  ostringstream eherr;
1084  eherr << "field " << fieldname.c_str () << "cannot be read.";
1085  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1086  }
1087  }
1088 
1089 
1090  if (speciallon && fieldtype == 2)
1091  CorSpeLon ((int16 *) &val[0], nelms);
1092 
1093  set_value ((dods_int16 *) &val[0], nelms);
1094  }
1095  break;
1096  case DFNT_UINT16:
1097  {
1098  uint16 fillvalue = 0;
1099  vector<uint16> val;
1100  val.resize(nelms);
1101 
1102  r = GDgetfillvalue (gridid,
1103  const_cast < char *>(fieldname.c_str ()),
1104  &fillvalue);
1105 
1106  if (r == 0) {
1107 
1108  int ifillvalue = fillvalue;
1109 
1110  vector <uint16> temp_total_val;
1111  temp_total_val.resize(xdim*ydim);
1112 
1113  r = readfieldfunc(gridid,
1114  const_cast < char *>(fieldname.c_str ()),
1115  NULL, NULL, NULL, (void *)(&temp_total_val[0]));
1116 
1117  if (r != 0) {
1118  detachfunc(gridid);
1119  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1120  ostringstream eherr;
1121  eherr << "field " << fieldname.c_str () << "cannot be read.";
1122  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1123  }
1124 
1125  try {
1126  HandleFillLatLon(temp_total_val, (uint16*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
1127  }
1128  catch(...) {
1129  detachfunc(gridid);
1130  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1131  throw;
1132  }
1133  }
1134 
1135  else {
1136 
1137  r = readfieldfunc (gridid,
1138  const_cast < char *>(fieldname.c_str ()),
1139  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
1140  if (r != 0) {
1141  detachfunc(gridid);
1142  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1143  ostringstream eherr;
1144  eherr << "field " << fieldname.c_str () << "cannot be read.";
1145  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1146  }
1147  }
1148 
1149  if (speciallon && fieldtype == 2)
1150  CorSpeLon ((uint16 *) &val[0], nelms);
1151 
1152  set_value ((dods_uint16 *) &val[0], nelms);
1153 
1154  }
1155  break;
1156 
1157  case DFNT_INT32:
1158  {
1159  vector<int32> val;
1160  val.resize(nelms);
1161 
1162  int32 fillvalue = 0;
1163 
1164  r = GDgetfillvalue (gridid,
1165  const_cast < char *>(fieldname.c_str ()),
1166  &fillvalue);
1167  if (r == 0) {
1168 
1169  int ifillvalue = fillvalue;
1170 
1171  vector <int32> temp_total_val;
1172  temp_total_val.resize(xdim*ydim);
1173 
1174  r = readfieldfunc(gridid,
1175  const_cast < char *>(fieldname.c_str ()),
1176  NULL, NULL, NULL, (void *)(&temp_total_val[0]));
1177 
1178  if (r != 0) {
1179  ostringstream eherr;
1180  detachfunc(gridid);
1181  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1182  eherr << "field " << fieldname.c_str () << "cannot be read.";
1183  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1184  }
1185 
1186  try {
1187  HandleFillLatLon(temp_total_val, (int32*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
1188  }
1189  catch(...) {
1190  detachfunc(gridid);
1191  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1192  throw;
1193  }
1194 
1195  }
1196 
1197  else {
1198 
1199  r = readfieldfunc (gridid,
1200  const_cast < char *>(fieldname.c_str ()),
1201  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
1202  if (r != 0) {
1203  detachfunc(gridid);
1204  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1205  ostringstream eherr;
1206  eherr << "field " << fieldname.c_str () << "cannot be read.";
1207  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1208  }
1209 
1210  }
1211  if (speciallon && fieldtype == 2)
1212  CorSpeLon ((int32 *) &val[0], nelms);
1213 
1214  set_value ((dods_int32 *) &val[0], nelms);
1215 
1216  }
1217  break;
1218  case DFNT_UINT32:
1219 
1220  {
1221  vector<uint32> val;
1222  val.resize(nelms);
1223 
1224  uint32 fillvalue = 0;
1225 
1226  r = GDgetfillvalue (gridid,
1227  const_cast < char *>(fieldname.c_str ()),
1228  &fillvalue);
1229  if (r == 0) {
1230 
1231  // this may cause overflow. Although we don't find the overflow in the NASA HDF products, may still fix it later. KY 2012-8-20
1232  int ifillvalue = (int)fillvalue;
1233  vector <uint32> temp_total_val;
1234  temp_total_val.resize(xdim*ydim);
1235  r = readfieldfunc(gridid,
1236  const_cast < char *>(fieldname.c_str ()),
1237  NULL, NULL, NULL, (void *)(&temp_total_val[0]));
1238 
1239  if (r != 0) {
1240  detachfunc(gridid);
1241  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1242  ostringstream eherr;
1243  eherr << "field " << fieldname.c_str () << "cannot be read.";
1244  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1245  }
1246 
1247  try {
1248  HandleFillLatLon(temp_total_val, (uint32*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
1249 
1250  }
1251  catch(...) {
1252  detachfunc(gridid);
1253  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1254  throw;
1255  }
1256  }
1257 
1258  else {
1259 
1260  r = readfieldfunc (gridid,
1261  const_cast < char *>(fieldname.c_str ()),
1262  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
1263  if (r != 0) {
1264  detachfunc(gridid);
1265  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1266  ostringstream eherr;
1267  eherr << "field " << fieldname.c_str () << "cannot be read.";
1268  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1269  }
1270 
1271  }
1272  if (speciallon && fieldtype == 2)
1273  CorSpeLon ((uint32 *) &val[0], nelms);
1274 
1275  set_value ((dods_uint32 *) &val[0], nelms);
1276 
1277  }
1278  break;
1279  case DFNT_FLOAT32:
1280 
1281  {
1282  vector<float32> val;
1283  val.resize(nelms);
1284 
1285  float32 fillvalue =0;
1286  r = GDgetfillvalue (gridid,
1287  const_cast < char *>(fieldname.c_str ()),
1288  &fillvalue);
1289 
1290 
1291  if (r == 0) {
1292  // May cause overflow,not find this happen in NASA HDF files, may still need to handle later.
1293  // KY 2012-08-20
1294  int ifillvalue =(int)fillvalue;
1295 
1296  vector <float32> temp_total_val;
1297  temp_total_val.resize(xdim*ydim);
1298  //temp_total_val.resize(xdim*ydim*4);
1299 
1300  r = readfieldfunc(gridid,
1301  const_cast < char *>(fieldname.c_str ()),
1302  NULL, NULL, NULL, (void *)(&temp_total_val[0]));
1303 
1304  if (r != 0) {
1305  detachfunc(gridid);
1306  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1307  ostringstream eherr;
1308  eherr << "field " << fieldname.c_str () << "cannot be read.";
1309  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1310  }
1311 
1312  try {
1313  HandleFillLatLon(temp_total_val, (float32*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
1314  }
1315  catch(...) {
1316  detachfunc(gridid);
1317  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1318  throw;
1319  }
1320 
1321  }
1322  else {
1323 
1324  r = readfieldfunc (gridid,
1325  const_cast < char *>(fieldname.c_str ()),
1326  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
1327  if (r != 0) {
1328  detachfunc(gridid);
1329  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1330  ostringstream eherr;
1331  eherr << "field " << fieldname.c_str () << "cannot be read.";
1332  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1333  }
1334 
1335  }
1336  if (speciallon && fieldtype == 2)
1337  CorSpeLon ((float32 *) &val[0], nelms);
1338 
1339  set_value ((dods_float32 *) &val[0], nelms);
1340 
1341  }
1342  break;
1343  case DFNT_FLOAT64:
1344 
1345  {
1346  vector<float64> val;
1347  val.resize(nelms);
1348 
1349  float64 fillvalue = 0;
1350  r = GDgetfillvalue (gridid,
1351  const_cast < char *>(fieldname.c_str ()),
1352  &fillvalue);
1353  if (r == 0) {
1354 
1355  // May cause overflow,not find this happen in NASA HDF files, may still need to handle later.
1356  // KY 2012-08-20
1357 
1358  int ifillvalue = (int)fillvalue;
1359  vector <float64> temp_total_val;
1360  temp_total_val.resize(xdim*ydim);
1361  r = readfieldfunc(gridid,
1362  const_cast < char *>(fieldname.c_str ()),
1363  NULL, NULL, NULL, (void *)(&temp_total_val[0]));
1364 
1365  if (r != 0) {
1366  detachfunc(gridid);
1367  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1368  ostringstream eherr;
1369  eherr << "field " << fieldname.c_str () << "cannot be read.";
1370  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1371  }
1372 
1373  try {
1374  HandleFillLatLon(temp_total_val, (float64*)&val[0],ydimmajor,fieldtype,xdim,ydim,&offset32[0],&count32[0],&step32[0],ifillvalue);
1375  }
1376  catch(...) {
1377  detachfunc(gridid);
1378  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1379  throw;
1380  }
1381 
1382  }
1383 
1384  else {
1385 
1386  r = readfieldfunc (gridid,
1387  const_cast < char *>(fieldname.c_str ()),
1388  &offset32[0], &step32[0], &count32[0], (void*)(&val[0]));
1389  if (r != 0) {
1390  detachfunc(gridid);
1391  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1392  ostringstream eherr;
1393  eherr << "field " << fieldname.c_str () << "cannot be read.";
1394  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1395  }
1396 
1397  }
1398  if (speciallon && fieldtype == 2)
1399  CorSpeLon ((float64 *) &val[0], nelms);
1400 
1401  set_value ((dods_float64 *) &val[0], nelms);
1402 
1403  }
1404  break;
1405  default:
1406  detachfunc(gridid);
1407  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1408  InternalErr (__FILE__, __LINE__, "unsupported data type.");
1409  }
1410 
1411  }
1412 
1413  r = detachfunc (gridid);
1414  if (r != 0) {
1415  ostringstream eherr;
1416  eherr << "Grid " << datasetname.c_str () << " cannot be detached.";
1417  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1418  }
1419 
1420 
1421  HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1422 #if 0
1423  r = closefunc (gfid);
1424  if (r != 0) {
1425  ostringstream eherr;
1426 
1427  eherr << "Grid " << filename.c_str () << " cannot be closed.";
1428  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1429  }
1430 #endif
1431 
1432 
1433  return false;
1434 }
1435 
1436 // Standard way of DAP handlers to pass the coordinates of the subsetted region to the handlers
1437 // Return the number of elements to read.
1438 int
1439 HDFEOS2ArrayGridGeoField::format_constraint (int *offset, int *step,
1440  int *count)
1441 {
1442  long nels = 1;
1443  int id = 0;
1444 
1445  Dim_iter p = dim_begin ();
1446 
1447  while (p != dim_end ()) {
1448 
1449  int start = dimension_start (p, true);
1450  int stride = dimension_stride (p, true);
1451  int stop = dimension_stop (p, true);
1452 
1453 
1454  // Check for illegical constraint
1455  if (stride < 0 || start < 0 || stop < 0 || start > stop) {
1456  ostringstream oss;
1457 
1458  oss << "Array/Grid hyperslab indices are bad: [" << start <<
1459  ":" << stride << ":" << stop << "]";
1460  throw Error (malformed_expr, oss.str ());
1461  }
1462 
1463  // Check for an empty constraint and use the whole dimension if so.
1464  if (start == 0 && stop == 0 && stride == 0) {
1465  start = dimension_start (p, false);
1466  stride = dimension_stride (p, false);
1467  stop = dimension_stop (p, false);
1468  }
1469 
1470  offset[id] = start;
1471  step[id] = stride;
1472  count[id] = ((stop - start) / stride) + 1;// count of elements
1473  nels *= count[id]; // total number of values for variable
1474 
1475  BESDEBUG ("h4",
1476  "=format_constraint():"
1477  << "id=" << id << " offset=" << offset[id]
1478  << " step=" << step[id]
1479  << " count=" << count[id]
1480  << endl);
1481 
1482  id++;
1483  p++;
1484  }
1485 
1486  return nels;
1487 }
1488 
1489 // Calculate lat/lon based on HDF-EOS2 APIs.
1490 void
1491 HDFEOS2ArrayGridGeoField::CalculateLatLon (int32 gridid, int fieldtype,
1492  int specialformat,
1493  float64 * outlatlon,float64* latlon_all,
1494  int32 * offset, int32 * count,
1495  int32 * step, int nelms,bool write_latlon_cache)
1496 {
1497 
1498  // Retrieve dimensions and X-Y coordinates of corners
1499  int32 xdim = 0;
1500  int32 ydim = 0;
1501  int r = -1;
1502  float64 upleft[2];
1503  float64 lowright[2];
1504 
1505  r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
1506  if (r != 0) {
1507  ostringstream eherr;
1508  eherr << "cannot obtain grid information.";
1509  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1510  }
1511 
1512  // The coordinate values(MCD products) are set to -180.0, -90.0, etc.
1513  // We have to change them to DDDMMMSSS.SS format, so
1514  // we have to multiply them by 1000000.
1515  if (specialformat == 1) {
1516  upleft[0] = upleft[0] * 1000000;
1517  upleft[1] = upleft[1] * 1000000;
1518  lowright[0] = lowright[0] * 1000000;
1519  lowright[1] = lowright[1] * 1000000;
1520  }
1521 
1522  // The coordinate values(CERES TRMM) are set to default,which are zeros.
1523  // Based on the grid names and size, we find it covers the whole global.
1524  // So we set the corner coordinates to (-180000000.00,90000000.00) and
1525  // (180000000.00,-90000000.00).
1526  if (specialformat == 2) {
1527  upleft[0] = 0.0;
1528  upleft[1] = 90000000.0;
1529  lowright[0] = 360000000.0;
1530  lowright[1] = -90000000.0;
1531  }
1532 
1533  // Retrieve all GCTP projection information
1534  int32 projcode = 0;
1535  int32 zone = 0;
1536  int32 sphere = 0;
1537  float64 params[16];
1538 
1539  r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
1540  if (r != 0) {
1541  ostringstream eherr;
1542  eherr << "cannot obtain grid projection information";
1543  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1544  }
1545 
1546  // Retrieve pixel registration information
1547  int32 pixreg = 0;
1548 
1549  r = GDpixreginfo (gridid, &pixreg);
1550  if (r != 0) {
1551  ostringstream eherr;
1552  eherr << "cannot obtain grid pixel registration info.";
1553  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1554  }
1555 
1556 
1557  //Retrieve grid pixel origin
1558  int32 origin = 0;
1559 
1560  r = GDorigininfo (gridid, &origin);
1561  if (r != 0) {
1562  ostringstream eherr;
1563  eherr << "cannot obtain grid origin info.";
1564  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1565  }
1566 
1567  vector<int32>rows;
1568  vector<int32>cols;
1569  vector<float64>lon;
1570  vector<float64>lat;
1571  rows.resize(xdim*ydim);
1572  cols.resize(xdim*ydim);
1573  lon.resize(xdim*ydim);
1574  lat.resize(xdim*ydim);
1575 
1576 
1577  int i = 0, j = 0, k = 0;
1578 
1579  if (ydimmajor) {
1580  /* Fill two arguments, rows and columns */
1581  // rows cols
1582  // /- xdim -/ /- xdim -/
1583  // 0 0 0 ... 0 0 1 2 ... x
1584  // 1 1 1 ... 1 0 1 2 ... x
1585  // ... ...
1586  // y y y ... y 0 1 2 ... x
1587 
1588  for (k = j = 0; j < ydim; ++j) {
1589  for (i = 0; i < xdim; ++i) {
1590  rows[k] = j;
1591  cols[k] = i;
1592  ++k;
1593  }
1594  }
1595  }
1596  else {
1597  // rows cols
1598  // /- ydim -/ /- ydim -/
1599  // 0 1 2 ... y 0 0 0 ... y
1600  // 0 1 2 ... y 1 1 1 ... y
1601  // ... ...
1602  // 0 1 2 ... y 2 2 2 ... y
1603 
1604  for (k = j = 0; j < xdim; ++j) {
1605  for (i = 0; i < ydim; ++i) {
1606  rows[k] = i;
1607  cols[k] = j;
1608  ++k;
1609  }
1610  }
1611  }
1612 
1613 
1614  r = GDij2ll (projcode, zone, params, sphere, xdim, ydim, upleft, lowright,
1615  xdim * ydim, &rows[0], &cols[0], &lon[0], &lat[0], pixreg, origin);
1616 
1617  if (r != 0) {
1618  ostringstream eherr;
1619  eherr << "cannot calculate grid latitude and longitude";
1620  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1621  }
1622 
1623  // ADDING CACHE file routine,save lon and lat to a cached file. lat first, lon second.
1624  //vector<double>latlon_all;
1625  if(true == write_latlon_cache) {
1626  if(GCTP_CEA == projcode || GCTP_GEO == projcode) {
1627  //latlon_all.resize(xdim+ydim);
1628  vector<double>temp_lat;
1629  vector<double>temp_lon;
1630  int32 temp_offset[2];
1631  int32 temp_count[2];
1632  int32 temp_step[2];
1633  temp_offset[0] = 0;
1634  temp_offset[1] = 0;
1635  temp_step[0] = 1;
1636  temp_step[1] = 1;
1637  if(ydimmajor) {
1638  // Latitude
1639  temp_count[0] = ydim;
1640  temp_count[1] = 1;
1641  temp_lat.resize(ydim);
1642  LatLon2DSubset(&temp_lat[0],ydim,xdim,&lat[0],temp_offset,temp_count,temp_step);
1643 
1644  // Longitude
1645  temp_count[0] = 1;
1646  temp_count[1] = xdim;
1647  temp_lon.resize(xdim);
1648  LatLon2DSubset(&temp_lon[0],ydim,xdim,&lon[0],temp_offset,temp_count,temp_step);
1649 
1650  for(int i = 0; i<ydim;i++)
1651  latlon_all[i] = temp_lat[i];
1652 
1653  for(int i = 0; i<xdim;i++)
1654  latlon_all[i+ydim] = temp_lon[i];
1655 
1656  }
1657  else {
1658  // Latitude
1659  temp_count[1] = ydim;
1660  temp_count[0] = 1;
1661  temp_lat.resize(ydim);
1662  LatLon2DSubset(&temp_lat[0],xdim,ydim,&lat[0],temp_offset,temp_count,temp_step);
1663 
1664  // Longitude
1665  temp_count[1] = 1;
1666  temp_count[0] = xdim;
1667  temp_lon.resize(xdim);
1668  LatLon2DSubset(&temp_lon[0],xdim,ydim,&lon[0],temp_offset,temp_count,temp_step);
1669 
1670  for(int i = 0; i<ydim;i++)
1671  latlon_all[i] = temp_lat[i];
1672 
1673  for(int i = 0; i<xdim;i++)
1674  latlon_all[i+ydim] = temp_lon[i];
1675 
1676 
1677  }
1678  }
1679  else {
1680  memcpy((char*)(&latlon_all[0]),&lat[0],xdim*ydim*sizeof(double));
1681  memcpy((char*)(&latlon_all[0])+xdim*ydim*sizeof(double),&lon[0],xdim*ydim*sizeof(double));
1682  // memcpy(&latlon_all[0]+xdim*ydim*sizeof(double),&lon[0],xdim*ydim*sizeof(double));
1683 
1684  }
1685  }
1686 
1687  // 2-D Lat/Lon, need to decompose the data for subsetting.
1688  if (nelms == (xdim * ydim)) { // no subsetting return all, for the performance reason.
1689  if (fieldtype == 1)
1690  memcpy (outlatlon, &lat[0], xdim * ydim * sizeof (double));
1691  else
1692  memcpy (outlatlon, &lon[0], xdim * ydim * sizeof (double));
1693  }
1694  else { // Messy subsetting case, needs to know the major dimension
1695  if (ydimmajor) {
1696  if (fieldtype == 1) // Lat
1697  LatLon2DSubset (outlatlon, ydim, xdim, &lat[0], offset, count,
1698  step);
1699  else // Lon
1700  LatLon2DSubset (outlatlon, ydim, xdim, &lon[0], offset, count,
1701  step);
1702  }
1703  else {
1704  if (fieldtype == 1) // Lat
1705  LatLon2DSubset (outlatlon, xdim, ydim, &lat[0], offset, count,
1706  step);
1707  else // Lon
1708  LatLon2DSubset (outlatlon, xdim, ydim, &lon[0], offset, count,
1709  step);
1710  }
1711  }
1712 }
1713 
1714 
1715 // Map the subset of the lat/lon buffer to the corresponding 2D array.
1716 template<class T> void
1717 HDFEOS2ArrayGridGeoField::LatLon2DSubset (T * outlatlon, int majordim,
1718  int minordim, T * latlon,
1719  int32 * offset, int32 * count,
1720  int32 * step)
1721 {
1722  // float64 templatlon[majordim][minordim];
1723  T (*templatlonptr)[majordim][minordim] =
1724  (typeof templatlonptr) latlon;
1725  int i, j, k;
1726 
1727  // do subsetting
1728  // Find the correct index
1729  int dim0count = count[0];
1730  int dim1count = count[1];
1731  int dim0index[dim0count], dim1index[dim1count];
1732 
1733  for (i = 0; i < count[0]; i++) // count[0] is the least changing dimension
1734  dim0index[i] = offset[0] + i * step[0];
1735 
1736 
1737  for (j = 0; j < count[1]; j++)
1738  dim1index[j] = offset[1] + j * step[1];
1739 
1740  // Now assign the subsetting data
1741  k = 0;
1742 
1743  for (i = 0; i < count[0]; i++) {
1744  for (j = 0; j < count[1]; j++) {
1745  outlatlon[k] = (*templatlonptr)[dim0index[i]][dim1index[j]];
1746  k++;
1747 
1748  }
1749  }
1750 }
1751 
1752 // Some HDF-EOS2 geographic projection lat/lon fields have fill values.
1753 // This routine is used to replace those fill values by using the formula to calculate
1754 // the lat/lon of the geographic projection.
1755 template < class T > bool HDFEOS2ArrayGridGeoField::CorLatLon (T * latlon,
1756  int fieldtype,
1757  int elms,
1758  int fv)
1759 {
1760 
1761  // Since we only find the contiguous fill value of lat/lon from some position to the end
1762  // So to speed up the performance, the following algorithm is limited to that case.
1763 
1764  // The first two values cannot be fill value.
1765  // We find a special case :the first latitude(index 0) is a special value.
1766  // So we need to have three non-fill values to calculate the increment.
1767 
1768  if (elms < 3) {
1769  for (int i = 0; i < elms; i++)
1770  if ((int) (latlon[i]) == fv)
1771  return false;
1772  return true;
1773  }
1774 
1775  // Number of elements is greater than 3.
1776 
1777  for (int i = 0; i < 3; i++) // The first three elements should not include fill value.
1778  if ((int) (latlon[i]) == fv)
1779  return false;
1780 
1781  if ((int) (latlon[elms - 1]) != fv)
1782  return true;
1783 
1784  T increment = latlon[2] - latlon[1];
1785 
1786  int index = 0;
1787 
1788  // Find the first fill value
1789  index = findfirstfv (latlon, 0, elms - 1, fv);
1790  if (index < 2) {
1791  ostringstream eherr;
1792  eherr << "cannot calculate the fill value. ";
1793  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1794  }
1795 
1796  for (int i = index; i < elms; i++) {
1797 
1798  latlon[i] = latlon[i - 1] + increment;
1799 
1800  // The latitude must be within (-90,90)
1801  if (i != (elms - 1) && (fieldtype == 1) &&
1802  ((float) (latlon[i]) < -90.0 || (float) (latlon[i]) > 90.0))
1803  return false;
1804 
1805  // For longitude, since some files use (0,360)
1806  // some files use (-180,180), for simple check
1807  // we just choose (-180,360).
1808  // I haven't found longitude has missing values.
1809  if (i != (elms - 1) && (fieldtype == 2) &&
1810  ((float) (latlon[i]) < -180.0 || (float) (latlon[i]) > 360.0))
1811  return false;
1812  }
1813  if (fieldtype == 1) {
1814  if ((float) (latlon[elms - 1]) < -90.0)
1815  latlon[elms - 1] = (T)-90;
1816  if ((float) (latlon[elms - 1]) > 90.0)
1817  latlon[elms - 1] = (T)90;
1818  }
1819 
1820  if (fieldtype == 2) {
1821  if ((float) (latlon[elms - 1]) < -180.0)
1822  latlon[elms - 1] = (T)-180.0;
1823  if ((float) (latlon[elms - 1]) > 360.0)
1824  latlon[elms - 1] = (T)360.0;
1825  }
1826  return true;
1827 }
1828 
1829 // Make longitude (0-360) to (-180 - 180)
1830 template < class T > void
1831 HDFEOS2ArrayGridGeoField::CorSpeLon (T * lon, int xdim)
1832 {
1833  int i;
1834  float64 accuracy = 1e-3; // in case there is a lon value = 180.0 in the middle, make the error to be less than 1e-3.
1835  float64 temp = 0;
1836 
1837  // Check if this lon. field falls to the (0-360) case.
1838  int speindex = -1;
1839 
1840  for (i = 0; i < xdim; i++) {
1841  if ((double) lon[i] < 180.0)
1842  temp = 180.0 - (double) lon[i];
1843  if ((double) lon[i] > 180.0)
1844  temp = (double) lon[i] - 180.0;
1845 
1846  if (temp < accuracy) {
1847  speindex = i;
1848  break;
1849  }
1850  else if ((static_cast < double >(lon[i]) < 180.0)
1851  &&(static_cast<double>(lon[i + 1]) > 180.0)) {
1852  speindex = i;
1853  break;
1854  }
1855  else
1856  continue;
1857  }
1858 
1859  if (speindex != -1) {
1860  for (i = speindex + 1; i < xdim; i++) {
1861  lon[i] =
1862  static_cast < T > (static_cast < double >(lon[i]) - 360.0);
1863  }
1864  }
1865  return;
1866 }
1867 
1868 // Get correct subsetting indexes. This is especially useful when 2D lat/lon can be condensed to 1D.
1869 void
1870 HDFEOS2ArrayGridGeoField::getCorrectSubset (int *offset, int *count,
1871  int *step, int32 * offset32,
1872  int32 * count32, int32 * step32,
1873  bool condenseddim, bool ydimmajor,
1874  int fieldtype, int rank)
1875 {
1876 
1877  if (rank == 1) {
1878  offset32[0] = (int32) offset[0];
1879  count32[0] = (int32) count[0];
1880  step32[0] = (int32) step[0];
1881  }
1882  else if (condenseddim) {
1883 
1884  // Since offset,count and step for some dimensions will always
1885  // be 1, so first assign offset32,count32,step32 to 1.
1886  for (int i = 0; i < rank; i++) {
1887  offset32[i] = 0;
1888  count32[i] = 1;
1889  step32[i] = 1;
1890  }
1891 
1892  if (ydimmajor && fieldtype == 1) {// YDim major, User: Lat[YDim], File: Lat[YDim][XDim]
1893  offset32[0] = (int32) offset[0];
1894  count32[0] = (int32) count[0];
1895  step32[0] = (int32) step[0];
1896  }
1897  else if (ydimmajor && fieldtype == 2) { // YDim major, User: Lon[XDim],File: Lon[YDim][XDim]
1898  offset32[1] = (int32) offset[0];
1899  count32[1] = (int32) count[0];
1900  step32[1] = (int32) step[0];
1901  }
1902  else if (!ydimmajor && fieldtype == 1) {// XDim major, User: Lat[YDim], File: Lat[XDim][YDim]
1903  offset32[1] = (int32) offset[0];
1904  count32[1] = (int32) count[0];
1905  step32[1] = (int32) step[0];
1906  }
1907  else if (!ydimmajor && fieldtype == 2) {// XDim major, User: Lon[XDim], File: Lon[XDim][YDim]
1908  offset32[0] = (int32) offset[0];
1909  count32[0] = (int32) count[0];
1910  step32[0] = (int32) step[0];
1911  }
1912 
1913  else {// errors
1914  InternalErr (__FILE__, __LINE__,
1915  "Lat/lon subset is wrong for condensed lat/lon");
1916  }
1917  }
1918  else {
1919  for (int i = 0; i < rank; i++) {
1920  offset32[i] = (int32) offset[i];
1921  count32[i] = (int32) count[i];
1922  step32[i] = (int32) step[i];
1923  }
1924  }
1925 }
1926 
1927 // Correct latitude and longitude that have fill values. Although I only found this
1928 // happens for AIRS CO2 grids, I still implemented this as general as I can.
1929 
1930 template <class T> void
1931 HDFEOS2ArrayGridGeoField::HandleFillLatLon(vector<T> total_latlon, T* latlon,bool ydimmajor, int fieldtype, int32 xdim , int32 ydim, int32* offset, int32* count, int32* step, int fv) {
1932 
1933  class vector <T> temp_lat;
1934  class vector <T> temp_lon;
1935 
1936  if (true == ydimmajor) {
1937 
1938  if (1 == fieldtype) {
1939  temp_lat.resize(ydim);
1940  for (int i = 0; i <(int)ydim; i++)
1941  temp_lat[i] = total_latlon[i*xdim];
1942 
1943  if (false == CorLatLon(&temp_lat[0],fieldtype,ydim,fv))
1944  InternalErr(__FILE__,__LINE__,"Cannot handle the fill values in lat/lon correctly");
1945 
1946  for (int i = 0; i <(int)(count[0]); i++)
1947  latlon[i] = temp_lat[offset[0] + i* step[0]];
1948  }
1949  else {
1950 
1951  temp_lon.resize(xdim);
1952  for (int i = 0; i <(int)xdim; i++)
1953  temp_lon[i] = total_latlon[i];
1954 
1955 
1956  if (false == CorLatLon(&temp_lon[0],fieldtype,xdim,fv))
1957  InternalErr(__FILE__,__LINE__,"Cannot handle the fill values in lat/lon correctly");
1958 
1959  for (int i = 0; i <(int)(count[1]); i++)
1960  latlon[i] = temp_lon[offset[1] + i* step[1]];
1961 
1962  }
1963  }
1964  else {
1965 
1966  if (1 == fieldtype) {
1967  temp_lat.resize(xdim);
1968  for (int i = 0; i <(int)xdim; i++)
1969  temp_lat[i] = total_latlon[i];
1970 
1971  if (false == CorLatLon(&temp_lat[0],fieldtype,ydim,fv))
1972  InternalErr(__FILE__,__LINE__,"Cannot handle the fill values in lat/lon correctly");
1973 
1974  for (int i = 0; i <(int)(count[1]); i++)
1975  latlon[i] = temp_lat[offset[1] + i* step[1]];
1976  }
1977  else {
1978 
1979  temp_lon.resize(ydim);
1980  for (int i = 0; i <(int)ydim; i++)
1981  temp_lon[i] = total_latlon[i*xdim];
1982 
1983 
1984  if (false == CorLatLon(&temp_lon[0],fieldtype,xdim,fv))
1985  InternalErr(__FILE__,__LINE__,"Cannot handle the fill values in lat/lon correctly");
1986 
1987  for (int i = 0; i <(int)(count[0]); i++)
1988  latlon[i] = temp_lon[offset[0] + i* step[0]];
1989  }
1990 
1991  }
1992 }
1993 
1994 // A helper recursive function to find the first filled value index.
1995 template < class T > int
1996 HDFEOS2ArrayGridGeoField::findfirstfv (T * array, int start, int end,
1997  int fillvalue)
1998 {
1999 
2000  if (start == end || start == (end - 1)) {
2001  if (static_cast < int >(array[start]) == fillvalue)
2002  return start;
2003  else
2004  return end;
2005  }
2006  else {
2007  int current = (start + end) / 2;
2008 
2009  if (static_cast < int >(array[current]) == fillvalue)
2010  return findfirstfv (array, start, current, fillvalue);
2011  else
2012  return findfirstfv (array, current, end, fillvalue);
2013  }
2014 }
2015 
2016 // Calculate Special Latitude and Longitude.
2017 //One MOD13C2 file doesn't provide projection code
2018 // The upperleft and lowerright coordinates are all -1
2019 // We have to calculate lat/lon by ourselves.
2020 // Since it doesn't provide the project code, we double check their information
2021 // and find that it covers the whole globe with 0.05 degree resolution.
2022 // Lat. is from 90 to -90 and Lon is from -180 to 180.
2023 void
2024 HDFEOS2ArrayGridGeoField::CalculateSpeLatLon (int32 gridid, int fieldtype,
2025  float64 * outlatlon,
2026  int32 * offset32,
2027  int32 * count32, int32 * step32,
2028  int nelms)
2029 {
2030 
2031  // Retrieve dimensions and X-Y coordinates of corners
2032  int32 xdim = 0;
2033  int32 ydim = 0;
2034  int r = -1;
2035  float64 upleft[2];
2036  float64 lowright[2];
2037 
2038  r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
2039  if (r != 0) {
2040  ostringstream eherr;
2041  eherr << "cannot obtain grid information.";
2042  throw InternalErr (__FILE__, __LINE__, eherr.str ());
2043  }
2044  //Since this is a special calcuation out of using the GDij2ll function,
2045  // the rank is always assumed to be 2 and we condense to 1. So the
2046  // count for longitude should be count[1] instead of count[0]. See function GetCorSubset
2047 
2048  // Since the project parameters in StructMetadata are all set to be default, I will use
2049  // the default HDF-EOS2 cell center as the origin of the coordinate. See the HDF-EOS2 user's guide
2050  // for details. KY 2012-09-10
2051 
2052  if(0 == xdim || 0 == ydim)
2053  throw InternalErr(__FILE__,__LINE__,"xdim or ydim cannot be zero");
2054 
2055  if (fieldtype == 1) {
2056  double latstep = 180.0 / ydim;
2057 
2058  for (int i = 0; i < (int) (count32[0]); i++)
2059  outlatlon[i] = 90.0 -latstep/2 - latstep * (offset32[0] + i * step32[0]);
2060  }
2061  else {// Longitude should use count32[1] etc.
2062  double lonstep = 360.0 / xdim;
2063 
2064  for (int i = 0; i < (int) (count32[1]); i++)
2065  outlatlon[i] = -180.0 + lonstep/2 + lonstep * (offset32[1] + i * step32[1]);
2066  }
2067 }
2068 
2069 // Calculate latitude and longitude for the MISR SOM projection HDF-EOS2 product.
2070 // since the latitude and longitude of the SOM projection are 3-D, so we need to handle this projection in a special way.
2071 // Based on our current understanding, the third dimension size is always 180.
2072 // This is according to the MISR Lat/lon calculation document
2073 // at http://eosweb.larc.nasa.gov/PRODOCS/misr/DPS/DPS_v50_RevS.pdf
2074 void
2075 HDFEOS2ArrayGridGeoField::CalculateSOMLatLon(int32 gridid, int *start, int *count, int *step, int nelms,const string & cache_fpath,bool write_latlon_cache)
2076 {
2077  int32 projcode = -1;
2078  int32 zone = -1;
2079  int32 sphere = -1;
2080  float64 params[NPROJ];
2081  intn r = -1;
2082 
2083  r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
2084  if (r!=0)
2085  throw InternalErr (__FILE__, __LINE__, "GDprojinfo doesn't return the correct values");
2086 
2087  int MAXNDIM = 10;
2088  int32 dim[MAXNDIM];
2089  char dimlist[STRLEN];
2090  r = GDinqdims(gridid, dimlist, dim);
2091  // r is the number of dims. or 0.
2092  // So the valid returned value can be greater than 0. Only throw error when r is less than 0.
2093  if (r<0)
2094  throw InternalErr (__FILE__, __LINE__, "GDinqdims doesn't return the correct values");
2095 
2096  bool is_block_180 = false;
2097  for(int i=0; i<MAXNDIM; i++)
2098  {
2099  if(dim[i]==NBLOCK)
2100  {
2101  is_block_180 = true;
2102  break;
2103  }
2104  }
2105  if(false == is_block_180) {
2106  ostringstream eherr;
2107  eherr <<"Number of Block is not " << NBLOCK ;
2108  throw InternalErr(__FILE__,__LINE__,eherr.str());
2109  }
2110 
2111  int32 xdim = 0;
2112  int32 ydim = 0;
2113  float64 ulc[2];
2114  float64 lrc[2];
2115 
2116  r = GDgridinfo (gridid, &xdim, &ydim, ulc, lrc);
2117  if (r!=0)
2118  throw InternalErr(__FILE__,__LINE__,"GDgridinfo doesn't return the correct values");
2119 
2120 
2121  float32 offset[NOFFSET];
2122  char som_rw_code[]="r";
2123  r = GDblkSOMoffset(gridid, offset, NOFFSET, som_rw_code);
2124  if(r!=0)
2125  throw InternalErr(__FILE__,__LINE__,"GDblkSOMoffset doesn't return the correct values");
2126 
2127  int status = misr_init(NBLOCK, xdim, ydim, offset, ulc, lrc);
2128  if(status!=0)
2129  throw InternalErr(__FILE__,__LINE__,"misr_init doesn't return the correct values");
2130 
2131  long iflg = 0;
2132  int (*inv_trans[MAXPROJ+1])(double, double, double*, double*);
2133  inv_init((long)projcode, (long)zone, (double*)params, (long)sphere, NULL, NULL, (int*)&iflg, inv_trans);
2134  if(iflg)
2135  throw InternalErr(__FILE__,__LINE__,"inv_init doesn't return correct values");
2136 
2137  // Change to vector in the future. KY 2012-09-20
2138  //double *latlon = NULL;
2139  double somx = 0.;
2140  double somy = 0.;
2141  double lat_r = 0.;
2142  double lon_r = 0.;
2143  int i = 0;
2144  int j = 0;
2145  int k = 0;
2146  int b = 0;
2147  int npts=0;
2148  float l = 0;
2149  float s = 0;
2150 
2151  // Seems setting blockdim = 0 always, need to understand this more. KY 2012-09-20
2152  int blockdim=0; //20; //84.2115,84.2018, 84.192, ... //0 for all
2153  if(blockdim==0) //66.2263, 66.224, ....
2154  {
2155 
2156  if(true == write_latlon_cache) {
2157  vector<double>latlon_all;
2158  latlon_all.resize(xdim*ydim*NBLOCK*2);
2159  for(int i =1; i <NBLOCK+1;i++)
2160  for(int j=0;j<xdim;j++)
2161  for(int k=0;k<ydim;k++)
2162  {
2163  b = i;
2164  l = j;
2165  s = k;
2166  misrinv(b, l, s, &somx, &somy); /* (b,l.l,s.s) -> (X,Y) */
2167  sominv(somx, somy, &lon_r, &lat_r); /* (X,Y) -> (lat,lon) */
2168  latlon_all[npts] = lat_r*R2D;
2169  latlon_all[xdim*ydim*NBLOCK+npts] = lon_r*R2D;
2170  npts++;
2171 
2172  }
2173 #if 0
2174  // Not necessary here, it will be handled by the cached class.
2175  // Need to remove the file if the file size is not the size of the latlon array.
2176  //if(HDFCFUtil::write_vector_to_file(cache_fpath,latlon_all,sizeof(double)) != (xdim*ydim*NBLOCK*2)) {
2177  if(HDFCFUtil::write_vector_to_file2(cache_fpath,latlon_all,sizeof(double)) != (xdim*ydim*NBLOCK*2)*sizeof(double)) {
2178  if(remove(cache_fpath.c_str()) !=0) {
2179  throw InternalErr(__FILE__,__LINE__,"Cannot remove the cached file.");
2180  }
2181  }
2182 #endif
2183  BESH4Cache *llcache = BESH4Cache::get_instance();
2184  llcache->write_cached_data(cache_fpath,xdim*ydim*NBLOCK*2*sizeof(double),latlon_all);
2185 
2186  // Send the subset of latlon to DAP.
2187  vector<double>latlon;
2188  latlon.resize(nelms); //double[180*xdim*ydim];
2189  //int s1=start[0]+1, e1=s1+count[0]*step[0];
2190  //int s2=start[1], e2=s2+count[1]*step[1];
2191  //int s3=start[2], e3=s3+count[2]*step[2];
2192  //int s1=start[0]+1;
2193  //int s2=start[1];
2194  //int s3=start[2];
2195 
2196 
2197  npts =0;
2198  for(int i=0; i<count[0]; i++) //i = 1; i<180+1; i++)
2199  for(int j=0; j<count[1]; j++)//j=0; j<xdim; j++)
2200  for(int k=0; k<count[2]; k++)//k=0; k<ydim; k++)
2201  {
2202  if(fieldtype == 1) {
2203  latlon[npts] = latlon_all[start[0]*ydim*xdim+start[1]*ydim+start[2]+
2204  i*ydim*xdim*step[0]+j*ydim*step[1]+k*step[2]];
2205  }
2206  else {
2207  latlon[npts] = latlon_all[xdim*ydim*NBLOCK+start[0]*ydim*xdim+start[1]*ydim+start[2]+
2208  i*ydim*xdim*step[0]+j*ydim*step[1]+k*step[2]];
2209 
2210  }
2211  npts++;
2212  }
2213 
2214  set_value ((dods_float64 *) &latlon[0], nelms); //(180*xdim*ydim)); //nelms);
2215  }
2216  else {
2217  vector<double>latlon;
2218  latlon.resize(nelms); //double[180*xdim*ydim];
2219  int s1=start[0]+1, e1=s1+count[0]*step[0];
2220  int s2=start[1], e2=s2+count[1]*step[1];
2221  int s3=start[2], e3=s3+count[2]*step[2];
2222  for(i=s1; i<e1; i+=step[0]) //i = 1; i<180+1; i++)
2223  for(j=s2; j<e2; j+=step[1])//j=0; j<xdim; j++)
2224  for(k=s3; k<e3; k+=step[2])//k=0; k<ydim; k++)
2225  {
2226  b = i;
2227  l = j;
2228  s = k;
2229  misrinv(b, l, s, &somx, &somy); /* (b,l.l,s.s) -> (X,Y) */
2230  sominv(somx, somy, &lon_r, &lat_r); /* (X,Y) -> (lat,lon) */
2231  if(fieldtype==1)
2232  latlon[npts] = lat_r*R2D;
2233  else
2234  latlon[npts] = lon_r*R2D;
2235  npts++;
2236  }
2237  set_value ((dods_float64 *) &latlon[0], nelms); //(180*xdim*ydim)); //nelms);
2238  }
2239  }
2240  //if (latlon != NULL)
2241  // delete [] latlon;
2242 }
2243 
2244 // The following code aims to handle large MCD Grid(GCTP_GEO projection) such as 21600*43200 lat and lon.
2245 // These MODIS MCD files don't follow standard format for lat/lon (DDDMMMSSS);
2246 // they simply represent lat/lon as -180.0000000 or -90.000000.
2247 // HDF-EOS2 library won't give the correct value based on these value.
2248 // We need to calculate the latitude and longitude values.
2249 void
2250 HDFEOS2ArrayGridGeoField::CalculateLargeGeoLatLon(int32 gridid, int fieldtype, float64* latlon, float64* latlon_all,int *start, int *count, int *step, int nelms,bool write_latlon_cache)
2251 {
2252 
2253  int32 xdim = 0;
2254  int32 ydim = 0;
2255  float64 upleft[2];
2256  float64 lowright[2];
2257  int r = 0;
2258  r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
2259  if (r!=0) {
2260  throw InternalErr(__FILE__,__LINE__, "GDgridinfo failed");
2261  }
2262 
2263  if (0 == xdim || 0 == ydim) {
2264  throw InternalErr(__FILE__,__LINE__, "xdim or ydim should not be zero. ");
2265  }
2266 
2267  if (upleft[0]>180.0 || upleft[0] <-180.0 ||
2268  upleft[1]>90.0 || upleft[1] <-90.0 ||
2269  lowright[0] >180.0 || lowright[0] <-180.0 ||
2270  lowright[1] >90.0 || lowright[1] <-90.0) {
2271 
2272  throw InternalErr(__FILE__,__LINE__, "lat/lon corner points are out of range. ");
2273  }
2274 
2275  if (count[0] != nelms) {
2276  throw InternalErr(__FILE__,__LINE__, "rank is not 1 ");
2277  }
2278  float lat_step = (lowright[1] - upleft[1])/ydim;
2279  float lon_step = (lowright[0] - upleft[0])/xdim;
2280 
2281  if(true == write_latlon_cache) {
2282 
2283  for(int i = 0;i<ydim;i++)
2284  latlon_all[i] = upleft[1] + i*lat_step + lat_step/2;
2285 
2286  for(int i = 0;i<xdim;i++)
2287  latlon_all[i+ydim] = upleft[0] + i*lon_step + lon_step/2;
2288 
2289  }
2290 
2291  // Treat the origin of the coordinate as the center of the cell.
2292  // This has been the setting of MCD43 data. KY 2012-09-10
2293  if (1 == fieldtype) { //Latitude
2294  float start_lat = upleft[1] + start[0] *lat_step + lat_step/2;
2295  float step_lat = lat_step *step[0];
2296  for (int i = 0; i < count[0]; i++)
2297  latlon[i] = start_lat +i *step_lat;
2298  }
2299  else { // Longitude
2300  float start_lon = upleft[0] + start[0] *lon_step + lon_step/2;
2301  float step_lon = lon_step *step[0];
2302  for (int i = 0; i < count[0]; i++)
2303  latlon[i] = start_lon +i *step_lon;
2304  }
2305 
2306 }
2307 
2308 
2309 // Calculate latitude and longitude for LAMAZ projection lat/lon products.
2310 // GDij2ll returns infinite numbers over the north pole or the south pole.
2311 void
2312 HDFEOS2ArrayGridGeoField::CalculateLAMAZLatLon(int32 gridid, int fieldtype, float64* latlon, float64* latlon_all, int *start, int *count, int *step, int nelms,bool write_latlon_cache)
2313 {
2314  int32 xdim = 0;
2315  int32 ydim = 0;
2316  intn r = 0;
2317  float64 upleft[2];
2318  float64 lowright[2];
2319 
2320  r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
2321  if (r != 0)
2322  throw InternalErr(__FILE__,__LINE__,"GDgridinfo failed");
2323 
2324  vector<float64> tmp1;
2325  tmp1.resize(xdim*ydim);
2326  int32 tmp2[] = {0, 0};
2327  int32 tmp3[] = {xdim, ydim};
2328  int32 tmp4[] = {1, 1};
2329 
2330  CalculateLatLon (gridid, fieldtype, specialformat, &tmp1[0], latlon_all, tmp2, tmp3, tmp4, xdim*ydim,write_latlon_cache);
2331 
2332  if(write_latlon_cache == true) {
2333 
2334  vector<float64> temp_lat_all,lat_all;
2335  temp_lat_all.resize(xdim*ydim);
2336  lat_all.resize(xdim*ydim);
2337 
2338  vector<float64> temp_lon_all,lon_all;
2339  temp_lon_all.resize(xdim*ydim);
2340  lon_all.resize(xdim*ydim);
2341 
2342  for(int w=0; w < xdim*ydim; w++){
2343  temp_lat_all[w] = latlon_all[w];
2344  lat_all[w] = latlon_all[w];
2345  temp_lon_all[w] = latlon_all[w+xdim*ydim];
2346  lon_all[w] = latlon_all[w+xdim*ydim];
2347  }
2348 
2349  // If we find infinite number among lat or lon values, we use the nearest neighbor method to calculate lat or lon.
2350  if(ydimmajor) {
2351  for(int i=0; i<ydim; i++)//Lat
2352  for(int j=0; j<xdim; j++)
2353  if(isundef_lat(lat_all[i*xdim+j]))
2354  lat_all[i*xdim+j]=nearestNeighborLatVal(&temp_lat_all[0], i, j, ydim, xdim);
2355  for(int i=0; i<ydim; i++)
2356  for(int j=0; j<xdim; j++)
2357  if(isundef_lon(lon_all[i*xdim+j]))
2358  lon_all[i*xdim+j]=nearestNeighborLonVal(&temp_lon_all[0], i, j, ydim, xdim);
2359  }
2360  else { // end if(ydimmajor)
2361  for(int i=0; i<xdim; i++)
2362  for(int j=0; j<ydim; j++)
2363  if(isundef_lat(lat_all[i*ydim+j]))
2364  lat_all[i*ydim+j]=nearestNeighborLatVal(&temp_lat_all[0], i, j, xdim, ydim);
2365 
2366  for(int i=0; i<xdim; i++)
2367  for(int j=0; j<ydim; j++)
2368  if(isundef_lon(lon_all[i*ydim+j]))
2369  lon_all[i*ydim+j]=nearestNeighborLonVal(&temp_lon_all[0], i, j, xdim, ydim);
2370 
2371  }
2372 
2373  for(int i = 0; i<xdim*ydim;i++) {
2374  latlon_all[i] = lat_all[i];
2375  latlon_all[i+xdim*ydim] = lon_all[i];
2376  }
2377 
2378  }
2379 
2380  // Need to optimize the access of LAMAZ subset
2381  vector<float64> tmp5;
2382  tmp5.resize(xdim*ydim);
2383 
2384  for(int w=0; w < xdim*ydim; w++)
2385  tmp5[w] = tmp1[w];
2386 
2387  // If we find infinite number among lat or lon values, we use the nearest neighbor method to calculate lat or lon.
2388  if(ydimmajor) {
2389  if(fieldtype==1) {// Lat.
2390  for(int i=0; i<ydim; i++)
2391  for(int j=0; j<xdim; j++)
2392  if(isundef_lat(tmp1[i*xdim+j]))
2393  tmp1[i*xdim+j]=nearestNeighborLatVal(&tmp5[0], i, j, ydim, xdim);
2394  } else if(fieldtype==2){ // Lon.
2395  for(int i=0; i<ydim; i++)
2396  for(int j=0; j<xdim; j++)
2397  if(isundef_lon(tmp1[i*xdim+j]))
2398  tmp1[i*xdim+j]=nearestNeighborLonVal(&tmp5[0], i, j, ydim, xdim);
2399  }
2400  } else { // end if(ydimmajor)
2401  if(fieldtype==1) {
2402  for(int i=0; i<xdim; i++)
2403  for(int j=0; j<ydim; j++)
2404  if(isundef_lat(tmp1[i*ydim+j]))
2405  tmp1[i*ydim+j]=nearestNeighborLatVal(&tmp5[0], i, j, xdim, ydim);
2406  } else if(fieldtype==2) {
2407  for(int i=0; i<xdim; i++)
2408  for(int j=0; j<ydim; j++)
2409  if(isundef_lon(tmp1[i*ydim+j]))
2410  tmp1[i*ydim+j]=nearestNeighborLonVal(&tmp5[0], i, j, xdim, ydim);
2411  }
2412  }
2413 
2414  for(int i=start[0], k=0; i<start[0]+count[0]*step[0]; i+=step[0])
2415  for(int j=start[1]; j<start[1]+count[1]*step[1]; j+= step[1])
2416  latlon[k++] = tmp1[i*ydim+j];
2417 
2418 }
2419 #endif
2420 
STRLEN
#define STRLEN
Definition: misrproj.h:4
BESH4Cache::get_data_from_cache
bool get_data_from_cache(const string &cache_file_name, const int expected_file_size, int &fd)
Definition: BESH4MCache.cc:199
std
STL namespace.
HDFCFUtil::close_fileid
static void close_fileid(int32 sdfd, int32 file_id, int32 gridfd, int32 swathfd, bool pass_fileid_key)
Close HDF4 and HDF-EOS2 file IDs. For performance reasons, we want to keep HDF-EOS2/HDF4 IDs open for...
Definition: HDFCFUtil.cc:2881
BESH4Cache::getCachePrefixFromConfig
static string getCachePrefixFromConfig()
Definition: BESH4MCache.cc:46
NBLOCK
#define NBLOCK
Definition: misrproj.h:5
NOFFSET
#define NOFFSET
Definition: misrproj.h:6
BESH4Cache::getCacheSizeFromConfig
static unsigned long getCacheSizeFromConfig()
Definition: BESH4MCache.cc:26
HDFCFUtil::get_int_str
static std::string get_int_str(int)
Definition: HDFCFUtil.cc:2803
misrinv
int misrinv(const int block, const float line, const float sample, double *x, double *y)
BESH4Cache::write_cached_data
bool write_cached_data(const string &cache_file_name, const int expected_file_size, const std::vector< double > &val)
Definition: BESH4MCache.cc:222
HDFCFUtil::get_double_str
static std::string get_double_str(double, int, int)
Definition: HDFCFUtil.cc:2770
NULL
#define NULL
Definition: wcsUtil.h:65
R2D
#define R2D
Definition: misrproj.h:7
BESFileLockingCache::purge_file
virtual void purge_file(const string &file)
Purge a single file from the cache.
Definition: BESFileLockingCache.cc:842
BESH4Cache::get_instance
static BESH4Cache * get_instance()
Get the default instance of the BESH4Cache object.
Definition: BESH4MCache.cc:158
BESFileLockingCache::get_read_lock
virtual bool get_read_lock(const string &target, int &fd)
Get a read-only lock on the file if it exists.
Definition: BESFileLockingCache.cc:448
misr_init
int misr_init(const int nblock, const int nline, const int nsample, const float relOff[NOFFSET], const double ulc_coord[], const double lrc_coord[])
NPROJ
#define NPROJ
Definition: misrproj.h:9
HDFCFUtil::read_vector_from_file
static ssize_t read_vector_from_file(int fd, vector< double > &, size_t)
Definition: HDFCFUtil.cc:2863
HDFCFUtil::check_beskeys
static bool check_beskeys(const std::string &key)
Check the BES key. This function will check a BES key specified at the file h4.conf.in. If the key's value is either true or yes. The handler claims to find a key and will do some operations. Otherwise, will do different operations. For example, One may find a line H4.EnableCF=true at h4.conf.in. That means, the HDF4 handler will handle the HDF4 files by following CF conventions.
Definition: HDFCFUtil.cc:17
BESDEBUG
#define BESDEBUG(x, y)
macro used to send debug information to the debug stream
Definition: BESDebug.h:64
BESH4Cache::getCacheDirFromConfig
static string getCacheDirFromConfig()
Definition: BESH4MCache.cc:66
BESFileLockingCache::unlock_and_close
virtual void unlock_and_close(const string &target)
Unlock the named file.
Definition: BESFileLockingCache.cc:578