bes  Updated for version 3.20.6
HDFEOS2Array_RealField.cc
1 // This file is part of the hdf4 data handler for the OPeNDAP data server.
3 // It retrieves the real field values.
4 // Authors: MuQun Yang <myang6@hdfgroup.org> Eunsoo Seo
5 // Copyright (c) 2010-2012 The HDF Group
7 #ifdef USE_HDFEOS2_LIB
8 
9 #include <iostream>
10 #include <sstream>
11 #include <cassert>
12 #include <debug.h>
13 #include "InternalErr.h"
14 #include <BESDebug.h>
15 #include <BESLog.h>
16 
17 #include "HDFCFUtil.h"
18 #include "HDFEOS2Array_RealField.h"
19 #include "dodsutil.h"
20 #include "HDF4RequestHandler.h"
21 
22 using namespace std;
23 using namespace libdap;
24 
25 #define SIGNED_BYTE_TO_INT32 1
26 
27 bool
28 HDFEOS2Array_RealField::read ()
29 {
30 
31  BESDEBUG("h4","Coming to HDFEOS2_Array_RealField read "<<endl);
32  if(length() == 0)
33  return true;
34 
35 #if 0
36  string check_pass_fileid_key_str="H4.EnablePassFileID";
37  bool check_pass_fileid_key = false;
38  check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
39 #endif
40 
41  bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
42 
43  // Declare offset, count and step
44  vector<int>offset;
45  offset.resize(rank);
46  vector<int>count;
47  count.resize(rank);
48  vector<int>step;
49  step.resize(rank);
50 
51  // Obtain offset,step and count from the client expression constraint
52  int nelms = 0;
53  nelms = format_constraint (&offset[0], &step[0], &count[0]);
54 
55  // Just declare offset,count and step in the int32 type.
56  vector<int32>offset32;
57  offset32.resize(rank);
58  vector<int32>count32;
59  count32.resize(rank);
60  vector<int32>step32;
61  step32.resize(rank);
62 
63  // Just obtain the offset,count and step in the datatype of int32.
64  for (int i = 0; i < rank; i++) {
65  offset32[i] = (int32) offset[i];
66  count32[i] = (int32) count[i];
67  step32[i] = (int32) step[i];
68  }
69 
70  // Define function pointers to handle both grid and swath
71  int32 (*openfunc) (char *, intn);
72  intn (*closefunc) (int32);
73  int32 (*attachfunc) (int32, char *);
74  intn (*detachfunc) (int32);
75  intn (*fieldinfofunc) (int32, char *, int32 *, int32 *, int32 *, char *);
76 
77  string datasetname;
78  if (swathname == "") {
79  openfunc = GDopen;
80  closefunc = GDclose;
81  attachfunc = GDattach;
82  detachfunc = GDdetach;
83  fieldinfofunc = GDfieldinfo;
84  datasetname = gridname;
85  }
86  else if (gridname == "") {
87  openfunc = SWopen;
88  closefunc = SWclose;
89  attachfunc = SWattach;
90  detachfunc = SWdetach;
91  fieldinfofunc = SWfieldinfo;
92  datasetname = swathname;
93  }
94  else
95  throw InternalErr (__FILE__, __LINE__, "It should be either grid or swath.");
96 
97  // Note gfid and gridid represent either swath or grid.
98  int32 gfid = 0;
99  int32 gridid = 0;
100 
101  if (true == isgeofile || false == check_pass_fileid_key) {
102 
103  // Obtain the EOS object ID(either grid or swath)
104  gfid = openfunc (const_cast < char *>(filename.c_str ()), DFACC_READ);
105  if (gfid < 0) {
106  ostringstream eherr;
107  eherr << "File " << filename.c_str () << " cannot be open.";
108  throw InternalErr (__FILE__, __LINE__, eherr.str ());
109  }
110  }
111  else
112  gfid = gsfd;
113 
114  // Attach the EOS object ID
115  gridid = attachfunc (gfid, const_cast < char *>(datasetname.c_str ()));
116  if (gridid < 0) {
117  close_fileid(gfid,-1);
118  ostringstream eherr;
119  eherr << "Grid/Swath " << datasetname.c_str () << " cannot be attached.";
120  throw InternalErr (__FILE__, __LINE__, eherr.str ());
121  }
122 
123 #if 0
124  string check_disable_scale_comp_key = "H4.DisableScaleOffsetComp";
125  bool turn_on_disable_scale_comp_key= false;
126  turn_on_disable_scale_comp_key = HDFCFUtil::check_beskeys(check_disable_scale_comp_key);
127 #endif
128 
129 
130  bool is_modis_l1b = false;
131  if("MODIS_SWATH_Type_L1B" == swathname)
132  is_modis_l1b = true;
133 
134  bool is_modis_vip = false;
135  if ("VIP_CMG_GRID" == gridname)
136  is_modis_vip = true;
137 
138  bool field_is_vdata = false;
139 
140  // HDF-EOS2 swath maps 1-D field as vdata. So we need to check if this field is vdata or SDS.
141  // Essentially we only call SDS attribute routines to retrieve MODIS scale,offset and
142  // fillvalue attributes since we don't
143  // find 1-D MODIS field has scale,offset and fillvalue attributes. We may need to visit
144  // this again in the future to see if we also need to support the handling of
145  // scale,offset,fillvalue via vdata routines. KY 2013-07-15
146  if (""==gridname) {
147 
148  int32 tmp_rank = 0;
149  char tmp_dimlist[1024];
150  int32 tmp_dims[rank];
151  int32 field_dtype = 0;
152  intn r = 0;
153 
154  r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
155  &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
156  if (r != 0) {
157  detachfunc(gridid);
158  close_fileid(gfid,-1);
159  ostringstream eherr;
160 
161  eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
162  throw InternalErr (__FILE__, __LINE__, eherr.str ());
163  }
164 
165  if (1 == tmp_rank)
166  field_is_vdata = true;
167  }
168 
169 
170  bool has_Key_attr = false;
171 
172  if (false == field_is_vdata) {
173 
174  // Obtain attribute values.
175  int32 sdfileid = -1;
176 
177  if (true == isgeofile || false == check_pass_fileid_key) {
178 
179  sdfileid = SDstart(const_cast < char *>(filename.c_str ()), DFACC_READ);
180 
181  if (FAIL == sdfileid) {
182  detachfunc(gridid);
183  close_fileid(gfid,-1);
184  ostringstream eherr;
185  eherr << "Cannot Start the SD interface for the file " << filename <<endl;
186  }
187  }
188  else
189  sdfileid = sdfd;
190 
191  int32 sdsindex = -1;
192  int32 sdsid = -1;
193  sdsindex = SDnametoindex(sdfileid, fieldname.c_str());
194  if (FAIL == sdsindex) {
195  detachfunc(gridid);
196  close_fileid(gfid,sdfileid);
197  ostringstream eherr;
198  eherr << "Cannot obtain the index of " << fieldname;
199  throw InternalErr (__FILE__, __LINE__, eherr.str ());
200  }
201 
202  sdsid = SDselect(sdfileid, sdsindex);
203  if (FAIL == sdsid) {
204  detachfunc(gridid);
205  close_fileid(gfid,sdfileid);
206  ostringstream eherr;
207  eherr << "Cannot obtain the SDS ID of " << fieldname;
208  throw InternalErr (__FILE__, __LINE__, eherr.str ());
209  }
210 
211  // Here we cannot check if SDfindattr fails since even SDfindattr fails it doesn't mean
212  // errors happen. If no such attribute can be found, SDfindattr still returns FAIL.
213  // The correct way is to use SDgetinfo and SDattrinfo to check if attributes
214  // "radiance_scales" etc exist.
215  // For the time being, I won't do this, due to the performance reason and code simplicity and also the
216  // very small chance of real FAIL for SDfindattr.
217  if(SDfindattr(sdsid, "Key")!=FAIL)
218  has_Key_attr = true;
219 
220  // Close the interfaces
221  SDendaccess(sdsid);
222  if (true == isgeofile || false == check_pass_fileid_key)
223  SDend(sdfileid);
224  }
225 
226  // USE a try-catch block to release the resources.
227  try {
228  if((false == is_modis_l1b) && (false == is_modis_vip)
229  &&(false == has_Key_attr) && (true == HDF4RequestHandler::get_disable_scaleoffset_comp()))
230  write_dap_data_disable_scale_comp(gridid,nelms,&offset32[0],&count32[0],&step32[0]);
231  else
232  write_dap_data_scale_comp(gridid,nelms,offset32,count32,step32);
233  }
234  catch(...) {
235  detachfunc(gridid);
236  close_fileid(gfid,-1);
237  throw;
238  }
239 
240  int32 r = -1;
241  r = detachfunc (gridid);
242  if (r != 0) {
243  close_fileid(gfid,-1);
244  ostringstream eherr;
245  eherr << "Grid/Swath " << datasetname.c_str () << " cannot be detached.";
246  throw InternalErr (__FILE__, __LINE__, eherr.str ());
247  }
248 
249 
250  if(true == isgeofile || false == check_pass_fileid_key) {
251  r = closefunc (gfid);
252  if (r != 0) {
253  ostringstream eherr;
254  eherr << "Grid/Swath " << filename.c_str () << " cannot be closed.";
255  throw InternalErr (__FILE__, __LINE__, eherr.str ());
256  }
257  }
258 
259  return false;
260 }
261 
262 int
263 HDFEOS2Array_RealField::write_dap_data_scale_comp(int32 gridid,
264  int nelms,
265  vector<int32>& offset32,
266  vector<int32>& count32,
267  vector<int32>& step32) {
268 
269 
270  BESDEBUG("h4",
271  "coming to HDFEOS2Array_RealField write_dap_data_scale_comp "
272  <<endl);
273 
274 #if 0
275  string check_pass_fileid_key_str="H4.EnablePassFileID";
276  bool check_pass_fileid_key = false;
277  check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
278 #endif
279  bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
280 
281  // Define function pointers to handle both grid and swath
282  intn (*fieldinfofunc) (int32, char *, int32 *, int32 *, int32 *, char *);
283  intn (*readfieldfunc) (int32, char *, int32 *, int32 *, int32 *, void *);
284 
285 
286  if (swathname == "") {
287  fieldinfofunc = GDfieldinfo;
288  readfieldfunc = GDreadfield;
289  }
290  else if (gridname == "") {
291  fieldinfofunc = SWfieldinfo;
292  readfieldfunc = SWreadfield;
293  }
294  else
295  throw InternalErr (__FILE__, __LINE__, "It should be either grid or swath.");
296 
297  // tmp_rank and tmp_dimlist are two dummy variables that are only used when calling fieldinfo.
298  int32 tmp_rank = 0;
299  char tmp_dimlist[1024];
300 
301  // field dimension sizes
302  int32 tmp_dims[rank];
303 
304  // field data type
305  int32 field_dtype = 0;
306 
307  // returned value of HDF4 and HDF-EOS2 APIs
308  intn r = 0;
309 
310  // Obtain the field info. We mainly need the datatype information
311  // to allocate the buffer to store the data
312  r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
313  &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
314  if (r != 0) {
315  ostringstream eherr;
316 
317  eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
318  throw InternalErr (__FILE__, __LINE__, eherr.str ());
319  }
320 
321 
322  // The following chunk of code until switch(field_dtype) handles MODIS level 1B,
323  // MOD29E1D Key and VIP products. The reason to keep the code this way is due to
324  // use of RECALCULATE macro. It is too much work to change it now. KY 2013-12-17
325  // MODIS level 1B reflectance and radiance fields have scale/offset arrays rather
326  // than one scale/offset value.
327  // So we need to handle these fields specially.
328  float *reflectance_offsets =NULL;
329  float *reflectance_scales =NULL;
330  float *radiance_offsets =NULL;
331  float *radiance_scales =NULL;
332 
333  // Attribute datatype, reused for several attributes
334  int32 attr_dtype = 0;
335 
336  // Number of elements for an attribute, reused
337  int32 temp_attrcount = 0;
338 
339  // Number of elements in an attribute
340  int32 num_eles_of_an_attr = 0;
341 
342  // Attribute(radiance_scales, reflectance_scales) index
343  int32 cf_modl1b_rr_attrindex = 0;
344 
345  // Attribute (radiance_offsets) index
346  int32 cf_modl1b_rr_attrindex2 = 0;
347 
348  // Attribute valid_range index
349  int32 cf_vr_attrindex = 0;
350 
351  // Attribute fill value index
352  int32 cf_fv_attrindex = 0;
353 
354  // Scale factor attribute index
355  int32 scale_factor_attr_index = 0;
356 
357  // Add offset attribute index
358  int32 add_offset_attr_index = 0;
359 
360  // Initialize scale
361  float scale = 1;
362 
363  // Intialize field_offset
364  float field_offset = 0;
365 
366  // Initialize fillvalue
367  float fillvalue = 0;
368 
369  // Initialize the original valid_min
370  float orig_valid_min = 0;
371 
372  // Initialize the original valid_max
373  float orig_valid_max = 0;
374 
375  // Some NSIDC products use the "Key" attribute to identify
376  // the discrete valid values(land, cloud etc).
377  // Since the valid_range attribute in these products may treat values
378  // identified by the Key attribute as invalid,
379  // we need to handle them in a special way.So set a flag here.
380  bool has_Key_attr = false;
381 
382  int32 sdfileid = -1;
383  if(sotype!=DEFAULT_CF_EQU) {
384 
385  bool field_is_vdata = false;
386 
387  // HDF-EOS2 swath maps 1-D field as vdata. So we need to check
388  // if this field is vdata or SDS.
389  // Essentially we only call SDS attribute routines to retrieve MODIS scale,
390  // offset and fillvalue
391  // attributes since we don't find 1-D MODIS field has scale,offset and
392  // fillvalue attributes.
393  // We may need to visit this again in the future to see
394  // if we also need to support the handling of scale,offset,fillvalue via
395  // vdata routines. KY 2013-07-15
396  if (""==gridname) {
397 
398  r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
399  &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
400  if (r != 0) {
401  ostringstream eherr;
402  eherr << "Field " << fieldname.c_str ()
403  << " information cannot be obtained.";
404  throw InternalErr (__FILE__, __LINE__, eherr.str ());
405  }
406 
407  if (1 == tmp_rank)
408  field_is_vdata = true;
409  }
411 #if 0
412 r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
413  &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
414 if (r != 0) {
415  ostringstream eherr;
416 
417  eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
418  throw InternalErr (__FILE__, __LINE__, eherr.str ());
419 }
420 
421 cerr<<"tmp_rank is "<<tmp_rank <<endl;
422 #endif
423 
424 
425  // For swath, we don't see any MODIS 1-D fields that have scale,offset
426  // and fill value attributes that need to be changed.
427  // So now we don't need to handle the vdata handling.
428  // Another reason is the possible change of the implementation
429  // of the SDS attribute handlings. That may be too costly.
430  // KY 2012-07-31
431 
432  if( false == field_is_vdata) {
433 
434  char attrname[H4_MAX_NC_NAME + 1];
435  vector<char> attrbuf;
436 
437  // Obtain attribute values.
438  if(false == isgeofile || false == check_pass_fileid_key)
439  sdfileid = sdfd;
440  else {
441  sdfileid = SDstart(const_cast < char *>(filename.c_str ()), DFACC_READ);
442  if (FAIL == sdfileid) {
443  ostringstream eherr;
444  eherr << "Cannot Start the SD interface for the file "
445  << filename;
446  throw InternalErr (__FILE__, __LINE__, eherr.str ());
447  }
448  }
449 
450  int32 sdsindex = -1;
451  int32 sdsid;
452  sdsindex = SDnametoindex(sdfileid, fieldname.c_str());
453  if (FAIL == sdsindex) {
454  if(true == isgeofile || false == check_pass_fileid_key)
455  SDend(sdfileid);
456  ostringstream eherr;
457  eherr << "Cannot obtain the index of " << fieldname;
458  throw InternalErr (__FILE__, __LINE__, eherr.str ());
459  }
460 
461  sdsid = SDselect(sdfileid, sdsindex);
462  if (FAIL == sdsid) {
463  if (true == isgeofile || false == check_pass_fileid_key)
464  SDend(sdfileid);
465  ostringstream eherr;
466  eherr << "Cannot obtain the SDS ID of " << fieldname;
467  throw InternalErr (__FILE__, __LINE__, eherr.str ());
468  }
469 
470 #if 0
471  char attrname[H4_MAX_NC_NAME + 1];
472  vector<char> attrbuf, attrbuf2;
473 
474  // Here we cannot check if SDfindattr fails or not since even SDfindattr fails it doesn't mean
475  // errors happen. If no such attribute can be found, SDfindattr still returns FAIL.
476  // The correct way is to use SDgetinfo and SDattrinfo to check if attributes "radiance_scales" etc exist.
477  // For the time being, I won't do this, due to the performance reason and code simplity and also the
478  // very small chance of real FAIL for SDfindattr.
479  cf_general_attrindex = SDfindattr(sdsid, "radiance_scales");
480  cf_general_attrindex2 = SDfindattr(sdsid, "radiance_offsets");
481 
482  // Obtain the values of radiance_scales and radiance_offsets if they are available.
483  if(cf_general_attrindex!=FAIL && cf_general_attrindex2!=FAIL)
484  {
485  intn ret = -1;
486  ret = SDattrinfo(sdsid, cf_general_attrindex, attrname, &attr_dtype, &temp_attrcount);
487  if (ret==FAIL)
488  {
489  SDendaccess(sdsid);
490  if(true == isgeofile)
491  SDend(sdfileid);
492  ostringstream eherr;
493  eherr << "Attribute 'radiance_scales' in " << fieldname.c_str () << " cannot be obtained.";
494  throw InternalErr (__FILE__, __LINE__, eherr.str ());
495  }
496  attrbuf.clear();
497  attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
498  ret = SDreadattr(sdsid, cf_general_attrindex, (VOIDP)&attrbuf[0]);
499  if (ret==FAIL)
500  {
501  release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
502  SDendaccess(sdsid);
503  if (true == isgeofile)
504  SDend(sdfileid);
505  ostringstream eherr;
506  eherr << "Attribute 'radiance_scales' in " << fieldname.c_str () << " cannot be obtained.";
507  throw InternalErr (__FILE__, __LINE__, eherr.str ());
508  }
509  ret = SDattrinfo(sdsid, cf_general_attrindex2, attrname, &attr_dtype, &temp_attrcount);
510  if (ret==FAIL)
511  {
512  release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
513  SDendaccess(sdsid);
514  if(true == isgeofile)
515  SDend(sdfileid);
516  ostringstream eherr;
517  eherr << "Attribute 'radiance_offsets' in " << fieldname.c_str () << " cannot be obtained.";
518  throw InternalErr (__FILE__, __LINE__, eherr.str ());
519  }
520  attrbuf2.clear();
521  attrbuf2.resize(DFKNTsize(attr_dtype)*temp_attrcount);
522  ret = SDreadattr(sdsid, cf_general_attrindex2, (VOIDP)&attrbuf2[0]);
523  if (ret==FAIL)
524  {
525  release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
526  SDendaccess(sdsid);
527  if(true == isgeofile)
528  SDend(sdfileid);
529  ostringstream eherr;
530  eherr << "Attribute 'radiance_offsets' in " << fieldname.c_str () << " cannot be obtained.";
531  throw InternalErr (__FILE__, __LINE__, eherr.str ());
532  }
533 
534  // The following macro will obtain radiance_scales and radiance_offsets.
535  // Although the code is compact, it may not be easy to follow. The similar macro can also be found later.
536  switch(attr_dtype)
537  {
538 #define GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(TYPE, CAST) \
539  case DFNT_##TYPE: \
540  { \
541  CAST *ptr = (CAST*)&attrbuf[0]; \
542  CAST *ptr2 = (CAST*)&attrbuf2[0]; \
543  radiance_scales = new float[temp_attrcount]; \
544  radiance_offsets = new float[temp_attrcount]; \
545  for(int l=0; l<temp_attrcount; l++) \
546  { \
547  radiance_scales[l] = ptr[l]; \
548  radiance_offsets[l] = ptr2[l]; \
549  } \
550  } \
551  break;
552  GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT32, float);
553  GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT64, double);
554  }
555 #undef GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES
556  num_eles_of_an_attr = temp_attrcount; // Store the count of attributes.
557  }
558 
559  // Obtain attribute values of reflectance_scales and reflectance_offsets if they are available.
560  cf_general_attrindex = SDfindattr(sdsid, "reflectance_scales");
561  cf_general_attrindex2 = SDfindattr(sdsid, "reflectance_offsets");
562  if(cf_general_attrindex!=FAIL && cf_general_attrindex2!=FAIL)
563  {
564  intn ret = -1;
565  ret = SDattrinfo(sdsid, cf_general_attrindex, attrname, &attr_dtype, &temp_attrcount);
566  if (ret==FAIL)
567  {
568  release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
569  SDendaccess(sdsid);
570  if(true == isgeofile)
571  SDend(sdfileid);
572  ostringstream eherr;
573  eherr << "Attribute 'reflectance_scales' in " << fieldname.c_str () << " cannot be obtained.";
574  throw InternalErr (__FILE__, __LINE__, eherr.str ());
575  }
576  attrbuf.clear();
577  attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
578  ret = SDreadattr(sdsid, cf_general_attrindex, (VOIDP)&attrbuf[0]);
579  if (ret==FAIL)
580  {
581  release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
582  SDendaccess(sdsid);
583  if(true == isgeofile)
584  SDend(sdfileid);
585  ostringstream eherr;
586  eherr << "Attribute 'reflectance_scales' in " << fieldname.c_str () << " cannot be obtained.";
587  throw InternalErr (__FILE__, __LINE__, eherr.str ());
588  }
589 
590  ret = SDattrinfo(sdsid, cf_general_attrindex2, attrname, &attr_dtype, &temp_attrcount);
591  if (ret==FAIL)
592  {
593  release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
594  SDendaccess(sdsid);
595  if(true == isgeofile)
596  SDend(sdfileid);
597  ostringstream eherr;
598  eherr << "Attribute 'reflectance_offsets' in " << fieldname.c_str () << " cannot be obtained.";
599  throw InternalErr (__FILE__, __LINE__, eherr.str ());
600  }
601  attrbuf2.clear();
602  attrbuf2.resize(DFKNTsize(attr_dtype)*temp_attrcount);
603  ret = SDreadattr(sdsid, cf_general_attrindex2, (VOIDP)&attrbuf2[0]);
604  if (ret==FAIL)
605  {
606  release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
607  SDendaccess(sdsid);
608  if(true == isgeofile)
609  SDend(sdfileid);
610  ostringstream eherr;
611  eherr << "Attribute 'reflectance_offsets' in " << fieldname.c_str () << " cannot be obtained.";
612  throw InternalErr (__FILE__, __LINE__, eherr.str ());
613  }
614  switch(attr_dtype)
615  {
616 #define GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(TYPE, CAST) \
617  case DFNT_##TYPE: \
618  { \
619  CAST *ptr = (CAST*)&attrbuf[0]; \
620  CAST *ptr2 = (CAST*)&attrbuf2[0]; \
621  reflectance_scales = new float[temp_attrcount]; \
622  reflectance_offsets = new float[temp_attrcount]; \
623  for(int l=0; l<temp_attrcount; l++) \
624  { \
625  reflectance_scales[l] = ptr[l]; \
626  reflectance_offsets[l] = ptr2[l]; \
627  } \
628  } \
629  break;
630  GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT32, float);
631  GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT64, double);
632  }
633 #undef GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES
634  num_eles_of_an_attr = temp_attrcount; // Store the count of attributes.
635  }
636 
637 #endif
638  // Obtain the value of attribute scale_factor.
639  scale_factor_attr_index = SDfindattr(sdsid, "scale_factor");
640  if(scale_factor_attr_index!=FAIL)
641  {
642  intn ret = -1;
643  ret = SDattrinfo(sdsid, scale_factor_attr_index, attrname,
644  &attr_dtype, &temp_attrcount);
645  if (ret==FAIL)
646  {
647  SDendaccess(sdsid);
648  if(true == isgeofile || false == check_pass_fileid_key)
649  SDend(sdfileid);
650  ostringstream eherr;
651  eherr << "Attribute 'scale_factor' in "
652  << fieldname.c_str () << " cannot be obtained.";
653  throw InternalErr (__FILE__, __LINE__, eherr.str ());
654  }
655  attrbuf.clear();
656  attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
657  ret = SDreadattr(sdsid, scale_factor_attr_index, (VOIDP)&attrbuf[0]);
658  if (ret==FAIL)
659  {
660  SDendaccess(sdsid);
661  if(true == isgeofile || false == check_pass_fileid_key)
662  SDend(sdfileid);
663 
664  ostringstream eherr;
665  eherr << "Attribute 'scale_factor' in "
666  << fieldname.c_str () << " cannot be obtained.";
667  throw InternalErr (__FILE__, __LINE__, eherr.str ());
668  }
669 
670  switch(attr_dtype)
671  {
672 #define GET_SCALE_FACTOR_ATTR_VALUE(TYPE, CAST) \
673  case DFNT_##TYPE: \
674  { \
675  CAST tmpvalue = *(CAST*)&attrbuf[0]; \
676  scale = (float)tmpvalue; \
677  } \
678  break;
679  GET_SCALE_FACTOR_ATTR_VALUE(INT8, int8);
680  GET_SCALE_FACTOR_ATTR_VALUE(CHAR,int8);
681  GET_SCALE_FACTOR_ATTR_VALUE(UINT8, uint8);
682  GET_SCALE_FACTOR_ATTR_VALUE(UCHAR,uint8);
683  GET_SCALE_FACTOR_ATTR_VALUE(INT16, int16);
684  GET_SCALE_FACTOR_ATTR_VALUE(UINT16, uint16);
685  GET_SCALE_FACTOR_ATTR_VALUE(INT32, int32);
686  GET_SCALE_FACTOR_ATTR_VALUE(UINT32, uint32);
687  GET_SCALE_FACTOR_ATTR_VALUE(FLOAT32, float);
688  GET_SCALE_FACTOR_ATTR_VALUE(FLOAT64, double);
689  default:
690  throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
691 
692 
693  };
694 #undef GET_SCALE_FACTOR_ATTR_VALUE
695  }
696 
697  // Obtain the value of attribute add_offset
698  add_offset_attr_index = SDfindattr(sdsid, "add_offset");
699  if(add_offset_attr_index!=FAIL)
700  {
701  intn ret;
702  ret = SDattrinfo(sdsid, add_offset_attr_index, attrname,
703  &attr_dtype, &temp_attrcount);
704  if (ret==FAIL)
705  {
706  SDendaccess(sdsid);
707  if(true == isgeofile || false == check_pass_fileid_key)
708  SDend(sdfileid);
709 
710  ostringstream eherr;
711  eherr << "Attribute 'add_offset' in " << fieldname.c_str ()
712  << " cannot be obtained.";
713  throw InternalErr (__FILE__, __LINE__, eherr.str ());
714  }
715  attrbuf.clear();
716  attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
717  ret = SDreadattr(sdsid, add_offset_attr_index, (VOIDP)&attrbuf[0]);
718  if (ret==FAIL)
719  {
720  SDendaccess(sdsid);
721  if(true == isgeofile || false == check_pass_fileid_key)
722  SDend(sdfileid);
723 
724  ostringstream eherr;
725  eherr << "Attribute 'add_offset' in " << fieldname.c_str ()
726  << " cannot be obtained.";
727  throw InternalErr (__FILE__, __LINE__, eherr.str ());
728  }
729 
730  switch(attr_dtype)
731  {
732 #define GET_ADD_OFFSET_ATTR_VALUE(TYPE, CAST) \
733  case DFNT_##TYPE: \
734  { \
735  CAST tmpvalue = *(CAST*)&attrbuf[0]; \
736  field_offset = (float)tmpvalue; \
737  } \
738  break;
739  GET_ADD_OFFSET_ATTR_VALUE(INT8, int8);
740  GET_ADD_OFFSET_ATTR_VALUE(CHAR,int8);
741  GET_ADD_OFFSET_ATTR_VALUE(UINT8, uint8);
742  GET_ADD_OFFSET_ATTR_VALUE(UCHAR,uint8);
743  GET_ADD_OFFSET_ATTR_VALUE(INT16, int16);
744  GET_ADD_OFFSET_ATTR_VALUE(UINT16, uint16);
745  GET_ADD_OFFSET_ATTR_VALUE(INT32, int32);
746  GET_ADD_OFFSET_ATTR_VALUE(UINT32, uint32);
747  GET_ADD_OFFSET_ATTR_VALUE(FLOAT32, float);
748  GET_ADD_OFFSET_ATTR_VALUE(FLOAT64, double);
749  default:
750  throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
751 
752  };
753 #undef GET_ADD_OFFSET_ATTR_VALUE
754  }
755 
756  // Obtain the value of the attribute _FillValue
757  cf_fv_attrindex = SDfindattr(sdsid, "_FillValue");
758  if(cf_fv_attrindex!=FAIL)
759  {
760  intn ret;
761  ret = SDattrinfo(sdsid, cf_fv_attrindex, attrname, &attr_dtype, &temp_attrcount);
762  if (ret==FAIL)
763  {
764  SDendaccess(sdsid);
765  if(true == isgeofile || false == check_pass_fileid_key)
766  SDend(sdfileid);
767 
768  ostringstream eherr;
769  eherr << "Attribute '_FillValue' in " << fieldname.c_str ()
770  << " cannot be obtained.";
771  throw InternalErr (__FILE__, __LINE__, eherr.str ());
772  }
773  attrbuf.clear();
774  attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
775  ret = SDreadattr(sdsid, cf_fv_attrindex, (VOIDP)&attrbuf[0]);
776  if (ret==FAIL)
777  {
778  SDendaccess(sdsid);
779  if(true == isgeofile || false == check_pass_fileid_key)
780  SDend(sdfileid);
781 
782  ostringstream eherr;
783  eherr << "Attribute '_FillValue' in " << fieldname.c_str ()
784  << " cannot be obtained.";
785  throw InternalErr (__FILE__, __LINE__, eherr.str ());
786  }
787 
788  switch(attr_dtype)
789  {
790 #define GET_FILLVALUE_ATTR_VALUE(TYPE, CAST) \
791  case DFNT_##TYPE: \
792  { \
793  CAST tmpvalue = *(CAST*)&attrbuf[0]; \
794  fillvalue = (float)tmpvalue; \
795  } \
796  break;
797  GET_FILLVALUE_ATTR_VALUE(INT8, int8);
798  GET_FILLVALUE_ATTR_VALUE(CHAR, int8);
799  GET_FILLVALUE_ATTR_VALUE(INT16, int16);
800  GET_FILLVALUE_ATTR_VALUE(INT32, int32);
801  GET_FILLVALUE_ATTR_VALUE(UINT8, uint8);
802  GET_FILLVALUE_ATTR_VALUE(UCHAR, uint8);
803  GET_FILLVALUE_ATTR_VALUE(UINT16, uint16);
804  GET_FILLVALUE_ATTR_VALUE(UINT32, uint32);
805  default:
806  throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
807 
808  };
809 #undef GET_FILLVALUE_ATTR_VALUE
810  }
811 
812  // Retrieve valid_range,valid_range is normally represented as (valid_min,valid_max)
813  // for non-CF scale and offset rules, the data is always float. So we only
814  // need to change the data type to float.
815  cf_vr_attrindex = SDfindattr(sdsid, "valid_range");
816  if(cf_vr_attrindex!=FAIL)
817  {
818  intn ret;
819  ret = SDattrinfo(sdsid, cf_vr_attrindex, attrname, &attr_dtype, &temp_attrcount);
820  if (ret==FAIL)
821  {
822  SDendaccess(sdsid);
823  if(true == isgeofile || false == check_pass_fileid_key)
824  SDend(sdfileid);
825 
826  ostringstream eherr;
827  eherr << "Attribute '_FillValue' in " << fieldname.c_str ()
828  << " cannot be obtained.";
829  throw InternalErr (__FILE__, __LINE__, eherr.str ());
830  }
831  attrbuf.clear();
832  attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
833  ret = SDreadattr(sdsid, cf_vr_attrindex, (VOIDP)&attrbuf[0]);
834  if (ret==FAIL)
835  {
836  SDendaccess(sdsid);
837  if(true == isgeofile || false == check_pass_fileid_key)
838  SDend(sdfileid);
839 
840  ostringstream eherr;
841  eherr << "Attribute '_FillValue' in " << fieldname.c_str ()
842  << " cannot be obtained.";
843  throw InternalErr (__FILE__, __LINE__, eherr.str ());
844  }
845 
846 
847  string attrbuf_str(attrbuf.begin(),attrbuf.end());
848 
849  switch(attr_dtype) {
850 
851  case DFNT_CHAR:
852  {
853  // We need to treat the attribute data as characters or string.
854  // So find the separator.
855  size_t found = attrbuf_str.find_first_of(",");
856  size_t found_from_end = attrbuf_str.find_last_of(",");
857 
858  if (string::npos == found){
859  SDendaccess(sdsid);
860  if(true == isgeofile || false == check_pass_fileid_key)
861  SDend(sdfileid);
862  throw InternalErr(__FILE__,__LINE__,"should find the separator ,");
863  }
864  if (found != found_from_end){
865  SDendaccess(sdsid);
866  if(true == isgeofile || false == check_pass_fileid_key)
867  SDend(sdfileid);
868  throw InternalErr(__FILE__,__LINE__,
869  "Only one separator , should be available.");
870  }
871 
872 #if 0
873  //istringstream(attrbuf_str.substr(0,found))>> orig_valid_min;
874  //istringstream(attrbuf_str.substr(found+1))>> orig_valid_max;
875 #endif
876 
877  orig_valid_min = atof((attrbuf_str.substr(0,found)).c_str());
878  orig_valid_max = atof((attrbuf_str.substr(found+1)).c_str());
879 
880  }
881  break;
882  case DFNT_INT8:
883  {
884  // We find a special case that even valid_range is logically
885  // interpreted as two elements,
886  // but the count of attribute elements is more than 2. The count
887  // actually is the number of
888  // characters stored as the attribute value. So we need to find
889  // the separator "," and then
890  // change the string before and after the separator into float numbers.
891  //
892  if (temp_attrcount >2) {
893 
894  size_t found = attrbuf_str.find_first_of(",");
895  size_t found_from_end = attrbuf_str.find_last_of(",");
896 
897  if (string::npos == found){
898  SDendaccess(sdsid);
899  if(true == isgeofile || false == check_pass_fileid_key)
900  SDend(sdfileid);
901  throw InternalErr(__FILE__,__LINE__,"should find the separator ,");
902  }
903  if (found != found_from_end){
904  SDendaccess(sdsid);
905  if(true == isgeofile || false == check_pass_fileid_key)
906  SDend(sdfileid);
907  throw InternalErr(__FILE__,__LINE__,
908  "Only one separator , should be available.");
909  }
910 
911  //istringstream(attrbuf_str.substr(0,found))>> orig_valid_min;
912  //istringstream(attrbuf_str.substr(found+1))>> orig_valid_max;
913 
914  orig_valid_min = atof((attrbuf_str.substr(0,found)).c_str());
915  orig_valid_max = atof((attrbuf_str.substr(found+1)).c_str());
916 
917  }
918  else if (2 == temp_attrcount) {
919  orig_valid_min = (float)attrbuf[0];
920  orig_valid_max = (float)attrbuf[1];
921  }
922  else {
923  SDendaccess(sdsid);
924  if(true == isgeofile || false == check_pass_fileid_key)
925  SDend(sdfileid);
926  throw InternalErr(__FILE__,__LINE__,
927  "The number of attribute count should be greater than 1.");
928  }
929 
930  }
931  break;
932 
933  case DFNT_UINT8:
934  case DFNT_UCHAR:
935  {
936  if (temp_attrcount != 2) {
937  SDendaccess(sdsid);
938  if(true == isgeofile || false == check_pass_fileid_key)
939  SDend(sdfileid);
940 
941  throw InternalErr(__FILE__,__LINE__,
942  "The number of attribute count should be 2 for the DFNT_UINT8 type.");
943  }
944 
945  unsigned char* temp_valid_range = (unsigned char *)&attrbuf[0];
946  orig_valid_min = (float)(temp_valid_range[0]);
947  orig_valid_max = (float)(temp_valid_range[1]);
948  }
949  break;
950 
951  case DFNT_INT16:
952  {
953  if (temp_attrcount != 2) {
954  SDendaccess(sdsid);
955  if(true == isgeofile || false == check_pass_fileid_key)
956  SDend(sdfileid);
957 
958  throw InternalErr(__FILE__,__LINE__,
959  "The number of attribute count should be 2 for the DFNT_INT16 type.");
960  }
961 
962  short* temp_valid_range = (short *)&attrbuf[0];
963  orig_valid_min = (float)(temp_valid_range[0]);
964  orig_valid_max = (float)(temp_valid_range[1]);
965  }
966  break;
967 
968  case DFNT_UINT16:
969  {
970  if (temp_attrcount != 2) {
971  SDendaccess(sdsid);
972  if(true == isgeofile || false == check_pass_fileid_key)
973  SDend(sdfileid);
974 
975  throw InternalErr(__FILE__,__LINE__,
976  "The number of attribute count should be 2 for the DFNT_UINT16 type.");
977  }
978 
979  unsigned short* temp_valid_range = (unsigned short *)&attrbuf[0];
980  orig_valid_min = (float)(temp_valid_range[0]);
981  orig_valid_max = (float)(temp_valid_range[1]);
982  }
983  break;
984 
985  case DFNT_INT32:
986  {
987  if (temp_attrcount != 2) {
988  SDendaccess(sdsid);
989  if(true == isgeofile || false == check_pass_fileid_key)
990  SDend(sdfileid);
991 
992  throw InternalErr(__FILE__,__LINE__,
993  "The number of attribute count should be 2 for the DFNT_INT32 type.");
994  }
995 
996  int* temp_valid_range = (int *)&attrbuf[0];
997  orig_valid_min = (float)(temp_valid_range[0]);
998  orig_valid_max = (float)(temp_valid_range[1]);
999  }
1000  break;
1001 
1002  case DFNT_UINT32:
1003  {
1004  if (temp_attrcount != 2) {
1005  SDendaccess(sdsid);
1006  if(true == isgeofile || false == check_pass_fileid_key)
1007  SDend(sdfileid);
1008 
1009  throw InternalErr(__FILE__,__LINE__,
1010  "The number of attribute count should be 2 for the DFNT_UINT32 type.");
1011  }
1012 
1013  unsigned int* temp_valid_range = (unsigned int *)&attrbuf[0];
1014  orig_valid_min = (float)(temp_valid_range[0]);
1015  orig_valid_max = (float)(temp_valid_range[1]);
1016  }
1017  break;
1018 
1019  case DFNT_FLOAT32:
1020  {
1021  if (temp_attrcount != 2) {
1022  SDendaccess(sdsid);
1023  if(true == isgeofile || false == check_pass_fileid_key)
1024  SDend(sdfileid);
1025 
1026  throw InternalErr(__FILE__,__LINE__,
1027  "The number of attribute count should be 2 for the DFNT_FLOAT32 type.");
1028  }
1029 
1030  float* temp_valid_range = (float *)&attrbuf[0];
1031  orig_valid_min = temp_valid_range[0];
1032  orig_valid_max = temp_valid_range[1];
1033  }
1034  break;
1035 
1036  case DFNT_FLOAT64:
1037  {
1038  if (temp_attrcount != 2){
1039  SDendaccess(sdsid);
1040  if(true == isgeofile || false == check_pass_fileid_key)
1041  SDend(sdfileid);
1042 
1043  throw InternalErr(__FILE__,__LINE__,
1044  "The number of attribute count should be 2 for the DFNT_FLOAT64 type.");
1045  }
1046  double* temp_valid_range = (double *)&attrbuf[0];
1047 
1048  // Notice: this approach will lose precision and possibly overflow.
1049  // Fortunately it is not a problem for MODIS data.
1050  // This part of code may not be called.
1051  // If it is called, mostly the value is within the floating-point range.
1052  // KY 2013-01-29
1053  orig_valid_min = temp_valid_range[0];
1054  orig_valid_max = temp_valid_range[1];
1055  }
1056  break;
1057  default: {
1058  SDendaccess(sdsid);
1059  if(true == isgeofile || false == check_pass_fileid_key)
1060  SDend(sdfileid);
1061  throw InternalErr(__FILE__,__LINE__,"Unsupported data type.");
1062  }
1063  }
1064  }
1065 
1066  // Check if the data has the "Key" attribute.
1067  // We found that some NSIDC MODIS data(MOD29) used "Key" to identify some special values.
1068  // To get the values that are within the range identified by the "Key",
1069  // scale offset rules also need to be applied to those values
1070  // outside the original valid range. KY 2013-02-25
1071  int32 cf_mod_key_attrindex = SUCCEED;
1072  cf_mod_key_attrindex = SDfindattr(sdsid, "Key");
1073  if(cf_mod_key_attrindex !=FAIL) {
1074  has_Key_attr = true;
1075  }
1076 
1077  attrbuf.clear();
1078  vector<char> attrbuf2;
1079 
1080  // Here we cannot check if SDfindattr fails since even SDfindattr fails it doesn't mean
1081  // errors happen. If no such attribute can be found, SDfindattr still returns FAIL.
1082  // The correct way is to use SDgetinfo and SDattrinfo to check if attributes
1083  // "radiance_scales" etc exist.
1084  // For the time being, I won't do this, due to the performance reason and code simplity
1085  // and also the very small chance of real FAIL for SDfindattr.
1086  cf_modl1b_rr_attrindex = SDfindattr(sdsid, "radiance_scales");
1087  cf_modl1b_rr_attrindex2 = SDfindattr(sdsid, "radiance_offsets");
1088 
1089  // Obtain the values of radiance_scales and radiance_offsets if they are available.
1090  if(cf_modl1b_rr_attrindex!=FAIL && cf_modl1b_rr_attrindex2!=FAIL)
1091  {
1092  intn ret = -1;
1093  ret = SDattrinfo(sdsid, cf_modl1b_rr_attrindex, attrname,
1094  &attr_dtype, &temp_attrcount);
1095  if (ret==FAIL)
1096  {
1097  SDendaccess(sdsid);
1098  if(true == isgeofile || false == check_pass_fileid_key)
1099  SDend(sdfileid);
1100  ostringstream eherr;
1101  eherr << "Attribute 'radiance_scales' in " << fieldname.c_str ()
1102  << " cannot be obtained.";
1103  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1104  }
1105  attrbuf.clear();
1106  attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
1107  ret = SDreadattr(sdsid, cf_modl1b_rr_attrindex, (VOIDP)&attrbuf[0]);
1108  if (ret==FAIL)
1109  {
1110  SDendaccess(sdsid);
1111  if (true == isgeofile || false == check_pass_fileid_key)
1112  SDend(sdfileid);
1113  ostringstream eherr;
1114  eherr << "Attribute 'radiance_scales' in " << fieldname.c_str ()
1115  << " cannot be obtained.";
1116  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1117  }
1118  ret = SDattrinfo(sdsid, cf_modl1b_rr_attrindex2, attrname,
1119  &attr_dtype, &temp_attrcount);
1120  if (ret==FAIL)
1121  {
1122  SDendaccess(sdsid);
1123  if(true == isgeofile || false == check_pass_fileid_key)
1124  SDend(sdfileid);
1125  ostringstream eherr;
1126  eherr << "Attribute 'radiance_offsets' in "
1127  << fieldname.c_str () << " cannot be obtained.";
1128  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1129  }
1130  attrbuf2.clear();
1131  attrbuf2.resize(DFKNTsize(attr_dtype)*temp_attrcount);
1132  ret = SDreadattr(sdsid, cf_modl1b_rr_attrindex2, (VOIDP)&attrbuf2[0]);
1133  if (ret==FAIL)
1134  {
1135  SDendaccess(sdsid);
1136  if(true == isgeofile || false == check_pass_fileid_key)
1137  SDend(sdfileid);
1138  ostringstream eherr;
1139  eherr << "Attribute 'radiance_offsets' in "
1140  << fieldname.c_str () << " cannot be obtained.";
1141  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1142  }
1143 
1144  // The following macro will obtain radiance_scales and radiance_offsets.
1145  // Although the code is compact, it may not be easy to follow.
1146  // The similar macro can also be found later.
1147  switch(attr_dtype)
1148  {
1149 #define GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(TYPE, CAST) \
1150  case DFNT_##TYPE: \
1151  { \
1152  CAST *ptr = (CAST*)&attrbuf[0]; \
1153  CAST *ptr2 = (CAST*)&attrbuf2[0]; \
1154  radiance_scales = new float[temp_attrcount]; \
1155  radiance_offsets = new float[temp_attrcount]; \
1156  for(int l=0; l<temp_attrcount; l++) \
1157  { \
1158  radiance_scales[l] = ptr[l]; \
1159  radiance_offsets[l] = ptr2[l]; \
1160  } \
1161  } \
1162  break;
1163  GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT32, float);
1164  GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT64, double);
1165  default:
1166  throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
1167 
1168  }
1169 #undef GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES
1170  // Store the count of attributes.
1171  num_eles_of_an_attr = temp_attrcount;
1172  }
1173 
1174  // Obtain attribute values of reflectance_scales
1175  // and reflectance_offsets if they are available.
1176  cf_modl1b_rr_attrindex = SDfindattr(sdsid, "reflectance_scales");
1177  cf_modl1b_rr_attrindex2 = SDfindattr(sdsid, "reflectance_offsets");
1178  if(cf_modl1b_rr_attrindex!=FAIL && cf_modl1b_rr_attrindex2!=FAIL)
1179  {
1180  intn ret = -1;
1181  ret = SDattrinfo(sdsid, cf_modl1b_rr_attrindex, attrname,
1182  &attr_dtype, &temp_attrcount);
1183  if (ret==FAIL)
1184  {
1185  release_mod1b_res(reflectance_scales,reflectance_offsets,
1186  radiance_scales,radiance_offsets);
1187  SDendaccess(sdsid);
1188  if(true == isgeofile || false == check_pass_fileid_key)
1189  SDend(sdfileid);
1190  ostringstream eherr;
1191  eherr << "Attribute 'reflectance_scales' in "
1192  << fieldname.c_str () << " cannot be obtained.";
1193  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1194  }
1195  attrbuf.clear();
1196  attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
1197  ret = SDreadattr(sdsid, cf_modl1b_rr_attrindex, (VOIDP)&attrbuf[0]);
1198  if (ret==FAIL)
1199  {
1200  release_mod1b_res(reflectance_scales,reflectance_offsets,
1201  radiance_scales,radiance_offsets);
1202  SDendaccess(sdsid);
1203  if(true == isgeofile || false == check_pass_fileid_key)
1204  SDend(sdfileid);
1205  ostringstream eherr;
1206  eherr << "Attribute 'reflectance_scales' in "
1207  << fieldname.c_str () << " cannot be obtained.";
1208  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1209  }
1210 
1211  ret = SDattrinfo(sdsid, cf_modl1b_rr_attrindex2, attrname,
1212  &attr_dtype, &temp_attrcount);
1213  if (ret==FAIL)
1214  {
1215  release_mod1b_res(reflectance_scales,reflectance_offsets,
1216  radiance_scales,radiance_offsets);
1217  SDendaccess(sdsid);
1218  if(true == isgeofile || false == check_pass_fileid_key)
1219  SDend(sdfileid);
1220  ostringstream eherr;
1221  eherr << "Attribute 'reflectance_offsets' in "
1222  << fieldname.c_str () << " cannot be obtained.";
1223  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1224  }
1225  attrbuf2.clear();
1226  attrbuf2.resize(DFKNTsize(attr_dtype)*temp_attrcount);
1227  ret = SDreadattr(sdsid, cf_modl1b_rr_attrindex2, (VOIDP)&attrbuf2[0]);
1228  if (ret==FAIL)
1229  {
1230  release_mod1b_res(reflectance_scales,reflectance_offsets,
1231  radiance_scales,radiance_offsets);
1232  SDendaccess(sdsid);
1233  if(true == isgeofile || false == check_pass_fileid_key)
1234  SDend(sdfileid);
1235  ostringstream eherr;
1236  eherr << "Attribute 'reflectance_offsets' in "
1237  << fieldname.c_str () << " cannot be obtained.";
1238  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1239  }
1240  switch(attr_dtype)
1241  {
1242 #define GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(TYPE, CAST) \
1243  case DFNT_##TYPE: \
1244  { \
1245  CAST *ptr = (CAST*)&attrbuf[0]; \
1246  CAST *ptr2 = (CAST*)&attrbuf2[0]; \
1247  reflectance_scales = new float[temp_attrcount]; \
1248  reflectance_offsets = new float[temp_attrcount]; \
1249  for(int l=0; l<temp_attrcount; l++) \
1250  { \
1251  reflectance_scales[l] = ptr[l]; \
1252  reflectance_offsets[l] = ptr2[l]; \
1253  } \
1254  } \
1255  break;
1256  GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT32, float);
1257  GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT64, double);
1258  default:
1259  throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
1260 
1261  }
1262 #undef GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES
1263  num_eles_of_an_attr = temp_attrcount; // Store the count of attributes.
1264  }
1265 
1266  SDendaccess(sdsid);
1268 #if 0
1269 r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
1270  &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1271 if (r != 0) {
1272  ostringstream eherr;
1273 
1274  eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
1275  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1276 }
1277 
1278 cerr<<"tmp_rank 3 is "<<tmp_rank <<endl;
1279 #endif
1280  //if (true == isgeofile || false == check_pass_fileid_key)
1281  // SDend(sdfileid);
1283 #if 0
1284 r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
1285  &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1286 if (r != 0) {
1287  ostringstream eherr;
1288 
1289  eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
1290  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1291 }
1292 
1293 cerr<<"tmp_rank 4 is "<<tmp_rank <<endl;
1294 #endif
1295 
1296  BESDEBUG("h4","scale is "<<scale <<endl);
1297  BESDEBUG("h4","offset is "<<field_offset <<endl);
1298  BESDEBUG("h4","fillvalue is "<<fillvalue <<endl);
1299  }
1300  }
1301 
1302  // According to our observations, it seems that MODIS products ALWAYS
1303  // use the "scale" factor to make bigger values smaller.
1304  // So for MODIS_MUL_SCALE products, if the scale of some variable is greater than 1,
1305  // it means that for this variable, the MODIS type for this variable may be MODIS_DIV_SCALE.
1306  // For the similar logic, we may need to change MODIS_DIV_SCALE to MODIS_MUL_SCALE
1307  // and MODIS_EQ_SCALE to MODIS_DIV_SCALE.
1308  // We indeed find such a case. HDF-EOS2 Grid MODIS_Grid_1km_2D of MOD(or MYD)09GA is
1309  // a MODIS_EQ_SCALE.
1310  // However,
1311  // the scale_factor of the variable Range_1 in the MOD09GA product is 25.
1312  // According to our observation,
1313  // this variable should be MODIS_DIV_SCALE.We verify this is true according to
1314  // MODIS 09 product document
1315  // http://modis-sr.ltdri.org/products/MOD09_UserGuide_v1_3.pdf.
1316  // Since this conclusion is based on our observation, we would like to add a BESlog to detect
1317  // if we find
1318  // the similar cases so that we can verify with the corresponding product documents to see if
1319  // this is true.
1320  // More updated information,
1321  // We just verified with the MOD09 data producer, the scale_factor for Range_1 is 25
1322  // but the equation is still multiplication instead of division.
1323  // So we have to make this as a special case and don't change the scale and offset settings
1324  // for Range_1 of MOD09 products.
1325  // KY 2014-01-13
1326 
1327  if (MODIS_EQ_SCALE == sotype || MODIS_MUL_SCALE == sotype) {
1328  if (scale > 1) {
1329  bool need_change_scale = true;
1330  if(gridname!="") {
1331 
1332  string temp_filename;
1333  if (filename.find("#") != string::npos)
1334  temp_filename =filename.substr(filename.find_last_of("#") + 1);
1335  else
1336  temp_filename = filename.substr(filename.find_last_of("/") +1);
1337 
1338  if ((temp_filename.size() >5) && ((temp_filename.compare(0,5,"MOD09") == 0)
1339  ||(temp_filename.compare(0,5,"MYD09") == 0))) {
1340  if ((fieldname.size() >5) && fieldname.compare(0,5,"Range") == 0)
1341  need_change_scale = false;
1342  }
1343  // MOD16A2
1344  else if((temp_filename.size() >7)&&
1345  ((temp_filename.compare(0,7,"MOD16A2") == 0)|| (temp_filename.compare(0,7,"MYD16A2")==0)||
1346  (temp_filename.compare(0,7,"MOD16A3") == 0)|| (temp_filename.compare(0,7,"MYD16A3")==0)))
1347  need_change_scale = false;
1348 
1349 
1350  }
1351  if(true == need_change_scale) {
1352  sotype = MODIS_DIV_SCALE;
1353  (*BESLog::TheLog())
1354  << "The field " << fieldname << " scale factor is "
1355  << scale << " ."<<endl
1356  << " But the original scale factor type is MODIS_MUL_SCALE or MODIS_EQ_SCALE. "
1357  << endl
1358  << " Now change it to MODIS_DIV_SCALE. "<<endl;
1359  }
1360  }
1361  }
1362 
1363  if (MODIS_DIV_SCALE == sotype) {
1364  if (scale < 1) {
1365  sotype = MODIS_MUL_SCALE;
1366  (*BESLog::TheLog())<< "The field " << fieldname << " scale factor is "
1367  << scale << " ."<<endl
1368  << " But the original scale factor type is MODIS_DIV_SCALE. "
1369  << endl
1370  << " Now change it to MODIS_MUL_SCALE. "<<endl;
1371  }
1372  }
1373 #if 0
1374 r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
1376  &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1377 if (r != 0) {
1378  ostringstream eherr;
1379 
1380  eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
1381  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1382 }
1383 
1384 cerr<<"tmp_rank 2 is "<<tmp_rank <<endl;
1385 #endif
1386 
1387 #if 0
1388  // We need to loop through all datatpes to allocate the memory buffer for the data.
1389 // It is hard to add comments to the macro. We may need to change them to general routines in the future.
1390 // Some MODIS products use both valid_range(valid_min, valid_max) and fillvalues for data fields. When do recalculating,
1391 // I check fillvalue first, then check valid_min and valid_max if they are available.
1392 // The middle check is_special_value addresses the MODIS L1B special value.
1393 // ////////////////////////////////////////////////////////////////////////////////////
1394 // if((float)tmptr[l] != fillvalue ) \
1395 // { \
1396 // if(false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[l]))\
1397  // { \
1398 // if (orig_valid_min<tmpval[l] && orig_valid_max>tmpval[l] \
1399 // if(sotype==MODIS_MUL_SCALE) \
1400 // tmpval[l] = (tmptr[l]-field_offset)*scale; \
1401 // else if(sotype==MODIS_EQ_SCALE) \
1402 // tmpval[l] = tmptr[l]*scale + field_offset; \
1403 // else if(sotype==MODIS_DIV_SCALE) \
1404 // tmpval[l] = (tmptr[l]-field_offset)/scale; \
1405 // } \
1406 #define RECALCULATE(CAST, DODS_CAST, VAL) \
1408 { \
1409  bool change_data_value = false; \
1410  if(sotype!=DEFAULT_CF_EQU) \
1411  { \
1412  vector<float>tmpval; \
1413  tmpval.resize(nelms); \
1414  CAST tmptr = (CAST)VAL; \
1415  for(int l=0; l<nelms; l++) \
1416  tmpval[l] = (float)tmptr[l]; \
1417  bool special_case = false; \
1418  if(scale_factor_attr_index==FAIL) \
1419  if(num_eles_of_an_attr==1) \
1420  if(radiance_scales!=NULL && radiance_offsets!=NULL) \
1421  { \
1422  scale = radiance_scales[0]; \
1423  field_offset = radiance_offsets[0];\
1424  special_case = true; \
1425  } \
1426  if((scale_factor_attr_index!=FAIL && !(scale==1 && field_offset==0)) || special_case) \
1427  { \
1428  for(int l=0; l<nelms; l++) \
1429  { \
1430  if(cf_vr_attrindex!=FAIL) \
1431  { \
1432  if((float)tmptr[l] != fillvalue ) \
1433  { \
1434  if(false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[l]))\
1435  { \
1436  if ((orig_valid_min<=tmpval[l] && orig_valid_max>=tmpval[l]) || (true==has_Key_attr))\
1437  { \
1438  if(sotype==MODIS_MUL_SCALE) \
1439  tmpval[l] = (tmptr[l]-field_offset)*scale; \
1440  else if(sotype==MODIS_EQ_SCALE) \
1441  tmpval[l] = tmptr[l]*scale + field_offset; \
1442  else if(sotype==MODIS_DIV_SCALE) \
1443  tmpval[l] = (tmptr[l]-field_offset)/scale;\
1444  } \
1445  } \
1446  } \
1447  } \
1448  } \
1449  change_data_value = true; \
1450  set_value((dods_float32 *)&tmpval[0], nelms); \
1451  } else if(num_eles_of_an_attr>1 && (radiance_scales!=NULL && radiance_offsets!=NULL) || (reflectance_scales!=NULL && reflectance_offsets!=NULL)) \
1452  { \
1453  size_t dimindex=0; \
1454  if( num_eles_of_an_attr!=tmp_dims[dimindex]) \
1455  { \
1456  ostringstream eherr; \
1457  eherr << "The number of Z-Dimension scale attribute is not equal to the size of the first dimension in " << fieldname.c_str() << ". These two values must be equal."; \
1458  throw InternalErr (__FILE__, __LINE__, eherr.str ()); \
1459  } \
1460  size_t start_index, end_index; \
1461  size_t nr_elems = nelms/count32[dimindex]; \
1462  start_index = offset32[dimindex]; \
1463  end_index = start_index+step32[dimindex]*(count32[dimindex]-1); \
1464  size_t index = 0;\
1465  for(size_t k=start_index; k<=end_index; k+=step32[dimindex]) \
1466  { \
1467  float tmpscale = (fieldname.find("Emissive")!=string::npos)? radiance_scales[k]: reflectance_scales[k]; \
1468  float tmpoffset = (fieldname.find("Emissive")!=string::npos)? radiance_offsets[k]: reflectance_offsets[k]; \
1469  for(size_t l=0; l<nr_elems; l++) \
1470  { \
1471  if(cf_vr_attrindex!=FAIL) \
1472  { \
1473  if(((float)tmptr[index])!=fillvalue) \
1474  { \
1475  if(false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[index]))\
1476  { \
1477  if(sotype==MODIS_MUL_SCALE) \
1478  tmpval[index] = (tmptr[index]-tmpoffset)*tmpscale; \
1479  else if(sotype==MODIS_EQ_SCALE) \
1480  tmpval[index] = tmptr[index]*tmpscale+tmpoffset; \
1481  else if(sotype==MODIS_DIV_SCALE) \
1482  tmpval[index] = (tmptr[index]-tmpoffset)/tmpscale; \
1483  } \
1484  } \
1485  } \
1486  index++; \
1487  } \
1488  } \
1489  change_data_value = true; \
1490  set_value((dods_float32 *)&tmpval[0], nelms); \
1491  } \
1492  } \
1493  if(!change_data_value) \
1494  { \
1495  set_value ((DODS_CAST)VAL, nelms); \
1496  } \
1497 }
1498 #endif
1499 
1500 // We need to loop through all datatpes to allocate the memory buffer for the data.
1501 // It is hard to add comments to the macro. We may need to change them to general
1502 // routines in the future.
1503 // Some MODIS products use both valid_range(valid_min, valid_max) and fillvalues for data fields.
1504 // When do recalculating,
1505 // I check fillvalue first, then check valid_min and valid_max if they are available.
1506 // The middle check is_special_value addresses the MODIS L1B special value.
1507 // Updated: just find that the RECALCULATE will be done only when the valid_range
1508 // attribute is present(if cf_vr_attrindex!=FAIL).
1509 // This restriction is in theory not necessary, but for more MODIS data products,
1510 // this restriction may be valid since valid_range pairs with scale/offset to identify
1511 // the valid data values. KY 2014-02-19
1513 // if((float)tmptr[l] != fillvalue ) \
1514 // { \
1515 // f(false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[l]))\
1516  // { \
1517 // if (orig_valid_min<tmpval[l] && orig_valid_max>tmpval[l] \
1518 // if(sotype==MODIS_MUL_SCALE) \
1519 // tmpval[l] = (tmptr[l]-field_offset)*scale; \
1520 // else if(sotype==MODIS_EQ_SCALE) \
1521 // tmpval[l] = tmptr[l]*scale + field_offset; \
1522 // else if(sotype==MODIS_DIV_SCALE) \
1523 // tmpval[l] = (tmptr[l]-field_offset)/scale; \
1524 // } \
1525 #define RECALCULATE(CAST, DODS_CAST, VAL) \
1527 { \
1528  bool change_data_value = false; \
1529  if(sotype!=DEFAULT_CF_EQU) \
1530  { \
1531  vector<float>tmpval; \
1532  tmpval.resize(nelms); \
1533  CAST tmptr = (CAST)VAL; \
1534  for(int l=0; l<nelms; l++) \
1535  tmpval[l] = (float)tmptr[l]; \
1536  bool special_case = false; \
1537  if(scale_factor_attr_index==FAIL) \
1538  if(num_eles_of_an_attr==1) \
1539  if((radiance_scales!=NULL) && (radiance_offsets!=NULL)) \
1540  { \
1541  scale = radiance_scales[0]; \
1542  field_offset = radiance_offsets[0];\
1543  special_case = true; \
1544  } \
1545  if(((scale_factor_attr_index!=FAIL) && !((scale==1) && (field_offset==0))) || special_case) \
1546  { \
1547  float temp_scale = scale; \
1548  float temp_offset = field_offset; \
1549  if(sotype==MODIS_MUL_SCALE) \
1550  temp_offset = -1. *field_offset*temp_scale;\
1551  else if (sotype==MODIS_DIV_SCALE) \
1552  {\
1553  temp_scale = 1/scale; \
1554  temp_offset = -1. *field_offset*temp_scale;\
1555  }\
1556  for(int l=0; l<nelms; l++) \
1557  { \
1558  if(cf_vr_attrindex!=FAIL) \
1559  { \
1560  if((float)tmptr[l] != fillvalue ) \
1561  { \
1562  if(false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[l]))\
1563  { \
1564  if (((orig_valid_min<=tmpval[l]) && (orig_valid_max>=tmpval[l])) || (true==has_Key_attr))\
1565  { \
1566  tmpval[l] = tmptr[l]*temp_scale + temp_offset; \
1567  } \
1568  } \
1569  } \
1570  } \
1571  } \
1572  change_data_value = true; \
1573  set_value((dods_float32 *)&tmpval[0], nelms); \
1574  } else if((num_eles_of_an_attr>1) && (((radiance_scales!=NULL) && (radiance_offsets!=NULL)) || ((reflectance_scales!=NULL) && (reflectance_offsets!=NULL)))) \
1575  { \
1576  size_t dimindex=0; \
1577  if( num_eles_of_an_attr!=tmp_dims[dimindex]) \
1578  { \
1579  release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets); \
1580  ostringstream eherr; \
1581  eherr << "The number of Z-Dimension scale attribute is not equal to the size of the first dimension in " << fieldname.c_str() << ". These two values must be equal."; \
1582  throw InternalErr (__FILE__, __LINE__, eherr.str ()); \
1583  } \
1584  size_t start_index, end_index; \
1585  size_t nr_elems = nelms/count32[dimindex]; \
1586  start_index = offset32[dimindex]; \
1587  end_index = start_index+step32[dimindex]*(count32[dimindex]-1); \
1588  size_t index = 0;\
1589  for(size_t k=start_index; k<=end_index; k+=step32[dimindex]) \
1590  { \
1591  float tmpscale = (fieldname.find("Emissive")!=string::npos)? radiance_scales[k]: reflectance_scales[k]; \
1592  float tmpoffset = (fieldname.find("Emissive")!=string::npos)? radiance_offsets[k]: reflectance_offsets[k]; \
1593  for(size_t l=0; l<nr_elems; l++) \
1594  { \
1595  if(cf_vr_attrindex!=FAIL) \
1596  { \
1597  if(((float)tmptr[index])!=fillvalue) \
1598  { \
1599  if(false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[index]))\
1600  { \
1601  if(sotype==MODIS_MUL_SCALE) \
1602  tmpval[index] = (tmptr[index]-tmpoffset)*tmpscale; \
1603  else if(sotype==MODIS_EQ_SCALE) \
1604  tmpval[index] = tmptr[index]*tmpscale+tmpoffset; \
1605  else if(sotype==MODIS_DIV_SCALE) \
1606  tmpval[index] = (tmptr[index]-tmpoffset)/tmpscale; \
1607  } \
1608  } \
1609  } \
1610  index++; \
1611  } \
1612  } \
1613  change_data_value = true; \
1614  set_value((dods_float32 *)&tmpval[0], nelms); \
1615  } \
1616  } \
1617  if(!change_data_value) \
1618  { \
1619  set_value ((DODS_CAST)VAL, nelms); \
1620  } \
1621 }
1622 #if 0
1623 r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
1625  &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1626 if (r != 0) {
1627  ostringstream eherr;
1628 
1629  eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
1630  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1631 }
1632 
1633 cerr<<"tmp_rank again is "<<tmp_rank <<endl;
1634 #endif
1635  switch (field_dtype) {
1636  case DFNT_INT8:
1637  {
1638 
1639  vector<int8>val;
1640  val.resize(nelms);
1641  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1642  &offset32[0], &step32[0], &count32[0], &val[0]);
1643  if (r != 0) {
1644  release_mod1b_res(reflectance_scales,reflectance_offsets,
1645  radiance_scales,radiance_offsets);
1646  ostringstream eherr;
1647  eherr << "field " << fieldname.c_str () << "cannot be read.";
1648  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1649  }
1650 
1651 #ifndef SIGNED_BYTE_TO_INT32
1652  RECALCULATE(int8*, dods_byte*, &val[0]);
1653 #else
1654 
1655  vector<int32>newval;
1656  newval.resize(nelms);
1657 
1658  for (int counter = 0; counter < nelms; counter++)
1659  newval[counter] = (int32) (val[counter]);
1660 
1661  RECALCULATE(int32*, dods_int32*, &newval[0]);
1662 #endif
1663  }
1664  break;
1665  case DFNT_UINT8:
1666  case DFNT_UCHAR8:
1667  {
1668 
1669  vector<uint8>val;
1670  val.resize(nelms);
1671  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1672  &offset32[0], &step32[0], &count32[0], &val[0]);
1673  if (r != 0) {
1674  release_mod1b_res(reflectance_scales,reflectance_offsets,
1675  radiance_scales,radiance_offsets);
1676  ostringstream eherr;
1677 
1678  eherr << "field " << fieldname.c_str () << "cannot be read.";
1679  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1680  }
1681 
1682  RECALCULATE(uint8*, dods_byte*, &val[0]);
1683  }
1684  break;
1685 
1686  case DFNT_INT16:
1687  {
1688  vector<int16>val;
1689  val.resize(nelms);
1690  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1691  &offset32[0], &step32[0], &count32[0], &val[0]);
1692 
1693  if (r != 0) {
1694 
1695  release_mod1b_res(reflectance_scales,reflectance_offsets,
1696  radiance_scales,radiance_offsets);
1697  ostringstream eherr;
1698 
1699  eherr << "field " << fieldname.c_str () << "cannot be read.";
1700  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1701  }
1702  RECALCULATE(int16*, dods_int16*, &val[0]);
1703  }
1704  break;
1705  case DFNT_UINT16:
1706  {
1707  vector<uint16>val;
1708  val.resize(nelms);
1709 #if 0
1710 cerr<<"gridid is "<<gridid <<endl;
1711 int32 tmp_rank = 0;
1712 char tmp_dimlist[1024];
1713 int32 tmp_dims[rank];
1714 int32 field_dtype = 0;
1715 intn r = 0;
1716 
1717 r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
1718  &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1719 if (r != 0) {
1720  ostringstream eherr;
1721 
1722  eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
1723  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1724 }
1725 #endif
1726 
1727  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1728  &offset32[0], &step32[0], &count32[0], &val[0]);
1729  if (r != 0) {
1730  release_mod1b_res(reflectance_scales,reflectance_offsets,
1731  radiance_scales,radiance_offsets);
1732  ostringstream eherr;
1733 
1734  eherr << "field " << fieldname.c_str () << "cannot be read.";
1735  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1736  }
1737 
1738  RECALCULATE(uint16*, dods_uint16*, &val[0]);
1739  }
1740  break;
1741  case DFNT_INT32:
1742  {
1743  vector<int32>val;
1744  val.resize(nelms);
1745  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1746  &offset32[0], &step32[0], &count32[0], &val[0]);
1747  if (r != 0) {
1748 
1749  release_mod1b_res(reflectance_scales,reflectance_offsets,
1750  radiance_scales,radiance_offsets);
1751  ostringstream eherr;
1752 
1753  eherr << "field " << fieldname.c_str () << "cannot be read.";
1754  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1755  }
1756 
1757  RECALCULATE(int32*, dods_int32*, &val[0]);
1758  }
1759  break;
1760  case DFNT_UINT32:
1761  {
1762  vector<uint32>val;
1763  val.resize(nelms);
1764  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1765  &offset32[0], &step32[0], &count32[0], &val[0]);
1766  if (r != 0) {
1767 
1768  release_mod1b_res(reflectance_scales,reflectance_offsets,
1769  radiance_scales,radiance_offsets);
1770  ostringstream eherr;
1771 
1772  eherr << "field " << fieldname.c_str () << "cannot be read.";
1773  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1774  }
1775 
1776  RECALCULATE(uint32*, dods_uint32*, &val[0]);
1777  }
1778  break;
1779  case DFNT_FLOAT32:
1780  {
1781  vector<float32>val;
1782  val.resize(nelms);
1783  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1784  &offset32[0], &step32[0], &count32[0], &val[0]);
1785  if (r != 0) {
1786 
1787  release_mod1b_res(reflectance_scales,reflectance_offsets,
1788  radiance_scales,radiance_offsets);
1789  ostringstream eherr;
1790 
1791  eherr << "field " << fieldname.c_str () << "cannot be read.";
1792  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1793  }
1794 
1795  // Recalculate seems not necessary.
1796  RECALCULATE(float32*, dods_float32*, &val[0]);
1797  //set_value((dods_float32*)&val[0],nelms);
1798  }
1799  break;
1800  case DFNT_FLOAT64:
1801  {
1802  vector<float64>val;
1803  val.resize(nelms);
1804  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1805  &offset32[0], &step32[0], &count32[0], &val[0]);
1806  if (r != 0) {
1807 
1808  release_mod1b_res(reflectance_scales,reflectance_offsets,
1809  radiance_scales,radiance_offsets);
1810  ostringstream eherr;
1811 
1812  eherr << "field " << fieldname.c_str () << "cannot be read.";
1813  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1814  }
1815  set_value ((dods_float64 *) &val[0], nelms);
1816  }
1817  break;
1818  default:
1819  release_mod1b_res(reflectance_scales,reflectance_offsets,
1820  radiance_scales,radiance_offsets);
1821  throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
1822  }
1823 
1824  release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
1825 #if 0
1826  if(reflectance_scales!=NULL)
1827  {
1828  delete[] reflectance_offsets;
1829  delete[] reflectance_scales;
1830  }
1831 
1832  if(radiance_scales!=NULL)
1833  {
1834  delete[] radiance_offsets;
1835  delete[] radiance_scales;
1836  }
1837 #endif
1838 
1839  // Somehow the macro RECALCULATE causes the interaction between gridid and sdfileid. SO
1840  // If I close the sdfileid earlier, gridid becomes invalid. So close the sdfileid now. KY 2014-10-24
1841  if (true == isgeofile || false == check_pass_fileid_key)
1842  SDend(sdfileid);
1843  //
1844  return false;
1845 
1846 }
1847 
1848 
1849 int
1850 HDFEOS2Array_RealField::write_dap_data_disable_scale_comp(int32 gridid,
1851  int nelms,
1852  int32 *offset32,
1853  int32*count32,
1854  int32*step32) {
1855 
1856 
1857  BESDEBUG("h4",
1858  "Coming to HDFEOS2_Array_RealField: write_dap_data_disable_scale_comp"
1859  <<endl);
1860 
1861  // Define function pointers to handle both grid and swath
1862  intn (*fieldinfofunc) (int32, char *, int32 *, int32 *, int32 *, char *);
1863  intn (*readfieldfunc) (int32, char *, int32 *, int32 *, int32 *, void *);
1864 
1865 
1866  if (swathname == "") {
1867  fieldinfofunc = GDfieldinfo;
1868  readfieldfunc = GDreadfield;
1869 
1870  }
1871  else if (gridname == "") {
1872  fieldinfofunc = SWfieldinfo;
1873  readfieldfunc = SWreadfield;
1874 
1875  }
1876  else
1877  throw InternalErr (__FILE__, __LINE__, "It should be either grid or swath.");
1878 
1879 
1880  // tmp_rank and tmp_dimlist are two dummy variables
1881  // that are only used when calling fieldinfo.
1882  int32 tmp_rank = 0;
1883  char tmp_dimlist[1024];
1884 
1885  // field dimension sizes
1886  int32 tmp_dims[rank];
1887 
1888  // field data type
1889  int32 field_dtype = 0;
1890 
1891  // returned value of HDF4 and HDF-EOS2 APIs
1892  intn r = 0;
1893 
1894  // Obtain the field info. We mainly need the datatype information
1895  // to allocate the buffer to store the data
1896  r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
1897  &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1898  if (r != 0) {
1899  ostringstream eherr;
1900  eherr << "Field " << fieldname.c_str ()
1901  << " information cannot be obtained.";
1902  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1903  }
1904 
1905 
1906  switch (field_dtype) {
1907  case DFNT_INT8:
1908  {
1909  vector<int8>val;
1910  val.resize(nelms);
1911  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1912  offset32, step32, count32, &val[0]);
1913  if (r != 0) {
1914  ostringstream eherr;
1915  eherr << "field " << fieldname.c_str () << "cannot be read.";
1916  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1917  }
1918 
1919 #ifndef SIGNED_BYTE_TO_INT32
1920  set_value((dods_byte*)&val[0],nelms);
1921 #else
1922 
1923  vector<int32>newval;
1924  newval.resize(nelms);
1925 
1926  for (int counter = 0; counter < nelms; counter++)
1927  newval[counter] = (int32) (val[counter]);
1928 
1929  set_value((dods_int32*)&newval[0],nelms);
1930 #endif
1931  }
1932  break;
1933  case DFNT_UINT8:
1934  case DFNT_UCHAR8:
1935  {
1936 
1937  vector<uint8>val;
1938  val.resize(nelms);
1939  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1940  offset32, step32, count32, &val[0]);
1941  if (r != 0) {
1942 
1943  ostringstream eherr;
1944  eherr << "field " << fieldname.c_str () << "cannot be read.";
1945  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1946  }
1947 
1948  set_value((dods_byte*)&val[0],nelms);
1949  }
1950  break;
1951 
1952  case DFNT_INT16:
1953  {
1954  vector<int16>val;
1955  val.resize(nelms);
1956  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1957  offset32, step32, count32, &val[0]);
1958 
1959  if (r != 0) {
1960  ostringstream eherr;
1961  eherr << "field " << fieldname.c_str () << "cannot be read.";
1962  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1963  }
1964  set_value((dods_int16*)&val[0],nelms);
1965  }
1966  break;
1967  case DFNT_UINT16:
1968  {
1969  vector<uint16>val;
1970  val.resize(nelms);
1971  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1972  offset32, step32, count32, &val[0]);
1973  if (r != 0) {
1974  ostringstream eherr;
1975  eherr << "field " << fieldname.c_str () << "cannot be read.";
1976  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1977  }
1978 
1979  set_value((dods_uint16*)&val[0],nelms);
1980  }
1981  break;
1982  case DFNT_INT32:
1983  {
1984  vector<int32>val;
1985  val.resize(nelms);
1986  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1987  offset32, step32, count32, &val[0]);
1988  if (r != 0) {
1989  ostringstream eherr;
1990  eherr << "field " << fieldname.c_str () << "cannot be read.";
1991  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1992  }
1993 
1994  set_value((dods_int32*)&val[0],nelms);
1995  }
1996  break;
1997  case DFNT_UINT32:
1998  {
1999  vector<uint32>val;
2000  val.resize(nelms);
2001  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
2002  offset32, step32, count32, &val[0]);
2003  if (r != 0) {
2004  ostringstream eherr;
2005  eherr << "field " << fieldname.c_str () << "cannot be read.";
2006  throw InternalErr (__FILE__, __LINE__, eherr.str ());
2007  }
2008 
2009  set_value((dods_uint32*)&val[0],nelms);
2010  }
2011  break;
2012  case DFNT_FLOAT32:
2013  {
2014  vector<float32>val;
2015  val.resize(nelms);
2016  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
2017  offset32, step32, count32, &val[0]);
2018  if (r != 0) {
2019  ostringstream eherr;
2020  eherr << "field " << fieldname.c_str () << "cannot be read.";
2021  throw InternalErr (__FILE__, __LINE__, eherr.str ());
2022  }
2023 
2024  // Recalculate seems not necessary.
2025  set_value((dods_float32*)&val[0],nelms);
2026  }
2027  break;
2028  case DFNT_FLOAT64:
2029  {
2030  vector<float64>val;
2031  val.resize(nelms);
2032  r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
2033  offset32, step32, count32, &val[0]);
2034  if (r != 0) {
2035  ostringstream eherr;
2036  eherr << "field " << fieldname.c_str () << "cannot be read.";
2037  throw InternalErr (__FILE__, __LINE__, eherr.str ());
2038  }
2039  set_value ((dods_float64 *) &val[0], nelms);
2040  }
2041  break;
2042  default:
2043  throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
2044  }
2045  return 0;
2046 }
2047 #if 0
2048  r = detachfunc (gridid);
2049  if (r != 0) {
2050  closefunc(gfid);
2051  ostringstream eherr;
2052 
2053  eherr << "Grid/Swath " << datasetname.c_str () << " cannot be detached.";
2054  throw InternalErr (__FILE__, __LINE__, eherr.str ());
2055  }
2056 
2057 
2058  r = closefunc (gfid);
2059  if (r != 0) {
2060  ostringstream eherr;
2061 
2062  eherr << "Grid/Swath " << filename.c_str () << " cannot be closed.";
2063  throw InternalErr (__FILE__, __LINE__, eherr.str ());
2064  }
2065 
2066  return false;
2067 }
2068 #endif
2069 
2070 // Standard way to pass the coordinates of the subsetted region from the client to the handlers
2071 // Return the number of elements to read.
2072 int
2073 HDFEOS2Array_RealField::format_constraint (int *offset, int *step, int *count)
2074 {
2075  long nels = 1;
2076  int id = 0;
2077 
2078  Dim_iter p = dim_begin ();
2079  while (p != dim_end ()) {
2080 
2081  int start = dimension_start (p, true);
2082  int stride = dimension_stride (p, true);
2083  int stop = dimension_stop (p, true);
2084 
2085  // Check for illegal constraint
2086  if (start > stop) {
2087  ostringstream oss;
2088  oss << "Array/Grid hyperslab start point "<< start <<
2089  " is greater than stop point " << stop <<".";
2090  throw Error(malformed_expr, oss.str());
2091  }
2092 
2093  offset[id] = start;
2094  step[id] = stride;
2095  count[id] = ((stop - start) / stride) + 1; // count of elements
2096  nels *= count[id]; // total number of values for variable
2097 
2098  BESDEBUG ("h4",
2099  "=format_constraint():"
2100  << "id=" << id << " offset=" << offset[id]
2101  << " step=" << step[id]
2102  << " count=" << count[id]
2103  << endl);
2104 
2105  id++;
2106  p++;
2107  }// while (p != dim_end ())
2108 
2109  return nels;
2110 }
2111 
2112 
2113 void HDFEOS2Array_RealField::close_fileid(const int gsfileid, const int sdfileid) {
2114 
2115 #if 0
2116  string check_pass_fileid_key_str="H4.EnablePassFileID";
2117  bool check_pass_fileid_key = false;
2118  check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
2119 #endif
2120 
2121 
2122  //if(true == isgeofile || false == check_pass_fileid_key) {
2123  if(true == isgeofile || false == HDF4RequestHandler::get_pass_fileid()) {
2124 
2125  if(sdfileid != -1)
2126  SDend(sdfileid);
2127 
2128  if(gsfileid != -1){
2129  if(""==gridname)
2130  SWclose(gsfileid);
2131  if (""==swathname)
2132  GDclose(gsfileid);
2133  }
2134 
2135  }
2136 
2137 }
2138 
2139 void HDFEOS2Array_RealField::release_mod1b_res(float*ref_scale,
2140  float*ref_offset,
2141  float*rad_scale,
2142  float*rad_offset) {
2143 
2144  if(ref_scale != NULL)
2145  delete[] ref_scale;
2146  if(ref_offset != NULL)
2147  delete[] ref_offset;
2148  if(rad_scale != NULL)
2149  delete[] rad_scale;
2150  if(rad_offset != NULL)
2151  delete[] rad_offset;
2152 
2153 }
2154 
2155 
2156 #endif
close_fileid
void close_fileid(hid_t fid)
Definition: h5get.cc:414
libdap
Definition: BESDapFunctionResponseCache.h:35
Error