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