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