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