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