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