bes  Updated for version 3.20.6
HDFEOS2ArraySwathDimMapField.cc
1 
3 // Retrieves the latitude and longitude of the HDF-EOS2 Swath with dimension map
4 // Authors: MuQun Yang <myang6@hdfgroup.org>
5 // Copyright (c) 2010-2012 The HDF Group
7 
8 // Currently the handling of swath data fields with dimension maps is the same as
9 // other data fields(HDFEOS2Array_RealField.cc etc)
10 // The reason to keep it in separate is, in theory, that data fields with dimension map
11 // may need special handlings.
12 // So we will leave it this way for now. It may be removed in the future.
13 // HDFEOS2Array_RealField.cc may be used.
14 // KY 2014-02-19
15 
16 #ifdef USE_HDFEOS2_LIB
17 
18 #include <iostream>
19 #include <sstream>
20 #include <cassert>
21 #include <debug.h>
22 #include "InternalErr.h"
23 #include "BESDebug.h"
24 #include <BESLog.h>
25 #include "HDFEOS2ArraySwathDimMapField.h"
26 #include "HDF4RequestHandler.h"
27 #define SIGNED_BYTE_TO_INT32 1
28 
29 using namespace std;
30 bool
31 HDFEOS2ArraySwathDimMapField::read ()
32 {
33 
34  BESDEBUG("h4","Coming to HDFEOS2ArraySwathDimMapField read "<<endl);
35  if(length() == 0)
36  return true;
37 
38 #if 0
39  string check_pass_fileid_key_str="H4.EnablePassFileID";
40  bool check_pass_fileid_key = false;
41  check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
42 #endif
43 
44  bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
45 
46  // Declare offset, count and step
47  vector<int>offset;
48  offset.resize(rank);
49 
50  vector<int>count;
51  count.resize(rank);
52 
53  vector<int>step;
54  step.resize(rank);
55 
56  // Obtain offset,step and count from the client expression constraint
57  int nelms = format_constraint(&offset[0],&step[0],&count[0]);
58 
59  // Just declare offset,count and step in the int32 type.
60  vector<int32>offset32;
61  offset32.resize(rank);
62 
63  vector<int32>count32;
64  count32.resize(rank);
65 
66  vector<int32>step32;
67  step32.resize(rank);
68 
69  // Just obtain the offset,count and step in the datatype of int32.
70  for (int i = 0; i < rank; i++) {
71  offset32[i] = (int32) offset[i];
72  count32[i] = (int32) count[i];
73  step32[i] = (int32) step[i];
74  }
75 
76  // Define function pointers to handle both grid and swath
77  int32 (*openfunc) (char *, intn);
78  intn (*closefunc) (int32);
79  int32 (*attachfunc) (int32, char *);
80  intn (*detachfunc) (int32);
81 
82  string datasetname;
83 
84  if (swathname == "") {
85  throw InternalErr (__FILE__, __LINE__, "It should be either grid or swath.");
86  }
87  else if (gridname == "") {
88  openfunc = SWopen;
89  closefunc = SWclose;
90  attachfunc = SWattach;
91  detachfunc = SWdetach;
92  datasetname = swathname;
93  }
94  else {
95  throw InternalErr (__FILE__, __LINE__, "It should be either grid or swath.");
96  }
97 
98  // Swath ID, swathid is actually in this case only the id of latitude and longitude.
99  int32 sfid = -1;
100  int32 swathid = -1;
101 
102  if (true == isgeofile || false == check_pass_fileid_key) {
103 
104  // Open, attach and obtain datatype information based on HDF-EOS2 APIs.
105  sfid = openfunc (const_cast < char *>(filename.c_str ()), DFACC_READ);
106 
107  if (sfid < 0) {
108  ostringstream eherr;
109  eherr << "File " << filename.c_str () << " cannot be open.";
110  throw InternalErr (__FILE__, __LINE__, eherr.str ());
111  }
112  }
113  else
114  sfid = swfd;
115 
116  swathid = attachfunc (sfid, const_cast < char *>(datasetname.c_str ()));
117  if (swathid < 0) {
118  close_fileid (sfid,-1);
119  ostringstream eherr;
120  eherr << "Grid/Swath " << datasetname.c_str () << " cannot be attached.";
121  throw InternalErr (__FILE__, __LINE__, eherr.str ());
122  }
123 
124  // dimmaps was set to be empty in hdfdesc.cc if the extra geolocation file also
125  // uses the dimension map
126  // This is because the dimmaps may be different in the MODIS geolocation file.
127  // So we cannot just pass
128  // the dimmaps to this class.
129  // Here we then obtain the dimension map info. in the geolocation file.
130  if(true == dimmaps.empty()) {
131 
132  int32 nummaps = 0;
133  int32 bufsize = 0;
134 
135  // Obtain number of dimension maps and the buffer size.
136  if ((nummaps = SWnentries(swathid, HDFE_NENTMAP, &bufsize)) == -1){
137  detachfunc(swathid);
138  close_fileid(sfid,-1);
139  throw InternalErr (__FILE__, __LINE__, "cannot obtain the number of dimmaps");
140  }
141 
142  if (nummaps <= 0){
143  detachfunc(swathid);
144  close_fileid(sfid,-1);
145  throw InternalErr (__FILE__,__LINE__,
146  "Number of dimension maps should be greater than 0");
147  }
148 
149  vector<char> namelist;
150  vector<int32> map_offset;
151  vector<int32> increment;
152 
153  namelist.resize(bufsize + 1);
154  map_offset.resize(nummaps);
155  increment.resize(nummaps);
156  if (SWinqmaps(swathid, &namelist[0], &map_offset[0], &increment[0])
157  == -1) {
158  detachfunc(swathid);
159  close_fileid(sfid,-1);
160  throw InternalErr (__FILE__,__LINE__,"fail to inquiry dimension maps");
161  }
162 
163  vector<string> mapnames;
164  HDFCFUtil::Split(&namelist[0], bufsize, ',', mapnames);
165  int map_count = 0;
166  for (vector<string>::const_iterator i = mapnames.begin();
167  i != mapnames.end(); ++i) {
168  vector<string> parts;
169  HDFCFUtil::Split(i->c_str(), '/', parts);
170  if (parts.size() != 2){
171  detachfunc(swathid);
172  close_fileid(sfid,-1);
173  throw InternalErr (__FILE__,__LINE__,"the dimmaps should only include two parts");
174  }
175 
176  struct dimmap_entry tempdimmap;
177  tempdimmap.geodim = parts[0];
178  tempdimmap.datadim = parts[1];
179  tempdimmap.offset = map_offset[map_count];
180  tempdimmap.inc = increment[map_count];
181 
182  dimmaps.push_back(tempdimmap);
183  ++map_count;
184  }
185  }
186 
187  if (sotype!=DEFAULT_CF_EQU) {
188 
189  if("MODIS_SWATH_Type_L1B" == swathname) {
190 
191  string emissive_str = "Emissive";
192  string RefSB_str = "RefSB";
193  bool is_emissive_field = false;
194  bool is_refsb_field = false;
195 
196  if(fieldname.find(emissive_str)!=string::npos) {
197  if(0 == fieldname.compare(fieldname.size()-emissive_str.size(),
198  emissive_str.size(),emissive_str))
199  is_emissive_field = true;
200  }
201 
202  if(fieldname.find(RefSB_str)!=string::npos) {
203  if(0 == fieldname.compare(fieldname.size()-RefSB_str.size(),
204  RefSB_str.size(),RefSB_str))
205  is_refsb_field = true;
206  }
207 
208  if ((true == is_emissive_field) || (true == is_refsb_field)) {
209  detachfunc(swathid);
210  close_fileid(sfid,-1);
211  throw InternalErr (__FILE__, __LINE__,
212  "Currently don't support MODIS Level 1B swath dim. map for data ");
213  }
214  }
215  }
216 
217  bool is_modis1b = false;
218  if("MODIS_SWATH_Type_L1B" == swathname)
219  is_modis1b = true;
220 #if 0
221  string check_disable_scale_comp_key = "H4.DisableScaleOffsetComp";
222  bool turn_on_disable_scale_comp_key= false;
223  turn_on_disable_scale_comp_key = HDFCFUtil::check_beskeys(check_disable_scale_comp_key);
224 #endif
225 
226  try {
227  //if(true == turn_on_disable_scale_comp_key && false== is_modis1b)
228  if(true == HDF4RequestHandler::get_disable_scaleoffset_comp() && false== is_modis1b)
229  write_dap_data_disable_scale_comp(swathid,nelms,offset32,count32,step32);
230  else
231  write_dap_data_scale_comp(swathid,nelms,offset32,count32,step32);
232  }
233  catch(...) {
234  detachfunc(swathid);
235  close_fileid(sfid,-1);
236  throw;
237  }
238 
239  intn r = 0;
240  r = detachfunc (swathid);
241  if (r != 0) {
242  close_fileid(sfid,-1);
243  ostringstream eherr;
244 
245  eherr << "Grid/Swath " << datasetname.c_str () << " cannot be detached.";
246  throw InternalErr (__FILE__, __LINE__, eherr.str ());
247  }
248 
249 
250  if(true == isgeofile || false == check_pass_fileid_key) {
251  r = closefunc (sfid);
252  if (r != 0) {
253  ostringstream eherr;
254  eherr << "Grid/Swath " << filename.c_str () << " cannot be closed.";
255  throw InternalErr (__FILE__, __LINE__, eherr.str ());
256  }
257  }
258 
259 
260  return false;
261 }
262 
263 // Standard way of DAP handlers to pass the coordinates of the subsetted region to the handlers
264 // Return the number of elements to read.
265 int
266 HDFEOS2ArraySwathDimMapField::format_constraint (int *offset, int *step, int *count)
267 {
268  long nels = 1;
269  int id = 0;
270 
271  Dim_iter p = dim_begin ();
272  while (p != dim_end ()) {
273 
274  int start = dimension_start (p, true);
275  int stride = dimension_stride (p, true);
276  int stop = dimension_stop (p, true);
277 
278  // Check for illegal constraint
279  if (start > stop) {
280  ostringstream oss;
281  oss << "Array/Grid hyperslab start point "<< start <<
282  " is greater than stop point " << stop <<".";
283  throw Error(malformed_expr, oss.str());
284  }
285 
286  offset[id] = start;
287  step[id] = stride;
288  count[id] = ((stop - start) / stride) + 1; // count of elements
289  nels *= count[id]; // total number of values for variable
290 
291  BESDEBUG ("h4",
292  "=format_constraint():"
293  << "id=" << id << " offset=" << offset[id]
294  << " step=" << step[id]
295  << " count=" << count[id]
296  << endl);
297 
298  id++;
299  p++;
300  }// while (p != dim_end ())
301 
302  return nels;
303 }
304 
305 // Get latitude and longitude fields.
306 // It will call expand_dimmap_field to interpolate latitude and longitude.
307 template < class T > int
308 HDFEOS2ArraySwathDimMapField::
309 GetFieldValue (int32 swathid, const string & geofieldname,
310  vector < struct dimmap_entry >&sw_dimmaps,
311  vector < T > &vals, vector<int32>&newdims)
312 {
313 
314  int32 ret = -1;
315  int32 size = -1;
316  int32 sw_rank = -1;
317  int32 dims[130];
318  int32 type = -1;
319 
320  // Two dimensions for lat/lon; each dimension name is < 64 characters,
321  // The dimension names are separated by a comma.
322  char dimlist[130];
323  ret = SWfieldinfo (swathid, const_cast < char *>(geofieldname.c_str ()),
324  &sw_rank, dims, &type, dimlist);
325  if (ret != 0)
326  return -1;
327 
328  size = 1;
329  for (int i = 0; i <sw_rank; i++)
330  size *= dims[i];
331 
332  vals.resize (size);
333 
334  ret = SWreadfield (swathid, const_cast < char *>(geofieldname.c_str ()),
335  NULL, NULL, NULL, (void *) &vals[0]);
336  if (ret != 0)
337  return -1;
338 
339  vector < string > dimname;
340  HDFCFUtil::Split (dimlist, ',', dimname);
341 
342  for (int i = 0; i < sw_rank; i++) {
343  vector < struct dimmap_entry >::iterator it;
344 
345  for (it = sw_dimmaps.begin (); it != sw_dimmaps.end (); it++) {
346  if (it->geodim == dimname[i]) {
347  int32 ddimsize = SWdiminfo (swathid, (char *) it->datadim.c_str ());
348  if (ddimsize == -1)
349  return -1;
350  int r;
351 
352  r = _expand_dimmap_field (&vals, sw_rank, dims, i, ddimsize, it->offset, it->inc);
353  if (r != 0)
354  return -1;
355  }
356  }
357  }
358 
359  // dims[] are expanded already.
360  for (int i = 0; i < sw_rank; i++) {
361  //cerr<<"i "<< i << " "<< dims[i] <<endl;
362  if (dims[i] < 0)
363  return -1;
364  newdims[i] = dims[i];
365  }
366 
367  return 0;
368 }
369 
370 // expand the dimension map field.
371 template < class T > int
372 HDFEOS2ArraySwathDimMapField::_expand_dimmap_field (vector < T >
373  *pvals, int32 sw_rank,
374  int32 dimsa[],
375  int dimindex,
376  int32 ddimsize,
377  int32 offset,
378  int32 inc)
379 {
380  vector < T > orig = *pvals;
381  vector < int32 > pos;
382  vector < int32 > dims;
383  vector < int32 > newdims;
384  pos.resize (sw_rank);
385  dims.resize (sw_rank);
386 
387  for (int i = 0; i < sw_rank; i++) {
388  pos[i] = 0;
389  dims[i] = dimsa[i];
390  }
391  newdims = dims;
392  newdims[dimindex] = ddimsize;
393  dimsa[dimindex] = ddimsize;
394 
395  int newsize = 1;
396 
397  for (int i = 0; i < sw_rank; i++) {
398  newsize *= newdims[i];
399  }
400  pvals->clear ();
401  pvals->resize (newsize);
402 
403  for (;;) {
404  // if end
405  if (pos[0] == dims[0]) {
406  // we past then end
407  break;
408  }
409  else if (pos[dimindex] == 0) {
410  // extract 1D values
411  vector < T > v;
412  for (int i = 0; i < dims[dimindex]; i++) {
413  pos[dimindex] = i;
414  v.push_back (orig[INDEX_nD_TO_1D (dims, pos)]);
415  }
416  // expand them
417 
418  vector < T > w;
419  for (int32 j = 0; j < ddimsize; j++) {
420  int32 i = (j - offset) / inc;
421  T f;
422 
423  if (i * inc + offset == j) // perfect match
424  {
425  f = (v[i]);
426  }
427  else {
428  int32 i1 = 0;
429  int32 i2 = (i<=0)?1:0;
430  int32 j1 = 0;
431  int32 j2 = 0;
432 
433 #if 0
434  if (i <= 0) {
435  //i1 = 0;
436  i2 = 1;
437  }
438 #endif
439  if ((unsigned int) i + 1 >= v.size ()) {
440  i1 = v.size () - 2;
441  i2 = v.size () - 1;
442  }
443  else {
444  i1 = i;
445  i2 = i + 1;
446  }
447  j1 = i1 * inc + offset;
448  j2 = i2 * inc + offset;
449  f = (((j - j1) * v[i2] + (j2 - j) * v[i1]) / (j2 - j1));
450  }
451  w.push_back (f);
452  pos[dimindex] = j;
453  (*pvals)[INDEX_nD_TO_1D (newdims, pos)] = f;
454  }
455  pos[dimindex] = 0;
456  }
457  // next pos
458  pos[sw_rank - 1]++;
459  for (int i = sw_rank - 1; i > 0; i--) {
460  if (pos[i] == dims[i]) {
461  pos[i] = 0;
462  pos[i - 1]++;
463  }
464  }
465  }
466 
467  return 0;
468 }
469 
470 template < class T >
471 bool HDFEOS2ArraySwathDimMapField::FieldSubset (T * outlatlon,
472  const vector<int32>&newdims,
473  T * latlon,
474  int32 * offset,
475  int32 * count,
476  int32 * step)
477 {
478 
479  if (newdims.size() == 1)
480  Field1DSubset(outlatlon,newdims[0],latlon,offset,count,step);
481  else if (newdims.size() == 2)
482  Field2DSubset(outlatlon,newdims[0],newdims[1],latlon,offset,count,step);
483  else if (newdims.size() == 3)
484  Field3DSubset(outlatlon,newdims,latlon,offset,count,step);
485  else
486  throw InternalErr(__FILE__, __LINE__,
487  "Currently doesn't support rank >3 when interpolating with dimension map");
488 
489  return true;
490 }
491 
492 // Subset of 1-D field to follow the parameters from the DAP expression constraint
493 template < class T >
494 bool HDFEOS2ArraySwathDimMapField::Field1DSubset (T * outlatlon,
495  const int majordim,
496  T * latlon,
497  int32 * offset,
498  int32 * count,
499  int32 * step)
500 {
501  if (majordim < count[0])
502  throw InternalErr(__FILE__, __LINE__,
503  "The number of elements is greater than the total dimensional size");
504 
505  for (int i = 0; i < count[0]; i++)
506  outlatlon[i] = latlon[offset[0]+i*step[0]];
507  return true;
508 
509 }
510 // Subset of latitude and longitude to follow the parameters
511 // from the DAP expression constraint
512 template < class T >
513 bool HDFEOS2ArraySwathDimMapField::Field2DSubset (T * outlatlon,
514  const int majordim,
515  const int minordim,
516  T * latlon,
517  int32 * offset,
518  int32 * count,
519  int32 * step)
520 {
521 #if 0
522  T (*templatlonptr)[majordim][minordim] = (T *[][]) latlon;
523 #endif
524  int i = 0;
525  int j = 0;
526 
527  // do subsetting
528  // Find the correct index
529  int dim0count = count[0];
530  int dim1count = count[1];
531 
532  int dim0index[dim0count];
533  int dim1index[dim1count];
534 
535  for (i = 0; i < count[0]; i++) // count[0] is the least changing dimension
536  dim0index[i] = offset[0] + i * step[0];
537 
538 
539  for (j = 0; j < count[1]; j++)
540  dim1index[j] = offset[1] + j * step[1];
541 
542  // Now assign the subsetting data
543  int k = 0;
544 
545  for (i = 0; i < count[0]; i++) {
546  for (j = 0; j < count[1]; j++) {
547 #if 0
548  outlatlon[k] = (*templatlonptr)[dim0index[i]][dim1index[j]];
549 #endif
550  outlatlon[k] = *(latlon + (dim0index[i] * minordim) + dim1index[j]);
551  k++;
552  }
553  }
554  return true;
555 }
556 
557 // Subsetting the field to follow the parameters from the DAP expression constraint
558 template < class T >
559 bool HDFEOS2ArraySwathDimMapField::Field3DSubset (T * outlatlon,
560  const vector<int32>& newdims,
561  T * latlon,
562  int32 * offset,
563  int32 * count,
564  int32 * step)
565 {
566  if (newdims.size() !=3)
567  throw InternalErr(__FILE__, __LINE__,
568  "the rank must be 3 to call this function");
569 #if 0
570  T (*templatlonptr)[newdims[0]][newdims[1]][newdims[2]] = (T *[][][]) latlon;
571 #endif
572  int i = 0;
573  int j = 0;
574  int k = 0;
575 
576  // do subsetting
577  // Find the correct index
578  int dim0count = count[0];
579  int dim1count = count[1];
580  int dim2count = count[2];
581 
582  int dim0index[dim0count], dim1index[dim1count],dim2index[dim2count];
583 
584  for (i = 0; i < count[0]; i++) // count[0] is the least changing dimension
585  dim0index[i] = offset[0] + i * step[0];
586 
587 
588  for (j = 0; j < count[1]; j++)
589  dim1index[j] = offset[1] + j * step[1];
590 
591  for (k = 0; k < count[2]; k++)
592  dim2index[k] = offset[2] + k * step[2];
593 
594  // Now assign the subsetting data
595  int l = 0;
596 
597  for (i = 0; i < count[0]; i++) {
598  for (j = 0; j < count[1]; j++) {
599  for (k =0; k < count[2]; k++) {
600 #if 0
601  outlatlon[l] = (*templatlonptr)[dim0index[i]][dim1index[j]][dim2index[k]];
602 #endif
603  outlatlon[l] = *(latlon + (dim0index[i] * newdims[1] * newdims[2]) + (dim1index[j] * newdims[2])+ dim2index[k]);
604  l++;
605  }
606  }
607  }
608  return true;
609 }
610 
611 int
612 HDFEOS2ArraySwathDimMapField::write_dap_data_scale_comp(int32 swathid,
613  int nelms,
614  vector<int32>& offset32,
615  vector<int32>& count32,
616  vector<int32>& step32) {
617 
618 #if 0
619  string check_pass_fileid_key_str="H4.EnablePassFileID";
620  bool check_pass_fileid_key = false;
621  check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
622 #endif
623 
624  bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
625 
626  // Define function pointers to handle both grid and swath
627  intn (*fieldinfofunc) (int32, char *, int32 *, int32 *, int32 *, char *);
628 
629 
630  fieldinfofunc = SWfieldinfo;
631 
632  int32 attrtype = -1;
633  int32 attrcount = -1;
634  int32 attrindex = -1;
635 
636  int32 scale_factor_attr_index = -1;
637  int32 add_offset_attr_index =-1;
638 
639  float scale=1;
640  float offset2=0;
641  float fillvalue = 0.;
642 
643  if (sotype!=DEFAULT_CF_EQU) {
644 
645  // Obtain attribute values.
646  int32 sdfileid = -1;
647 
648  if (true == isgeofile || false == check_pass_fileid_key) {
649  sdfileid = SDstart(const_cast < char *>(filename.c_str ()), DFACC_READ);
650  if (FAIL == sdfileid) {
651  ostringstream eherr;
652  eherr << "Cannot Start the SD interface for the file " << filename <<endl;
653  throw InternalErr (__FILE__, __LINE__, eherr.str ());
654  }
655  }
656  else
657  sdfileid = sdfd;
658 
659  int32 sdsindex = -1;
660  int32 sdsid = -1;
661 
662  sdsindex = SDnametoindex(sdfileid, fieldname.c_str());
663  if (FAIL == sdsindex) {
664  if(true == isgeofile || false == check_pass_fileid_key)
665  SDend(sdfileid);
666  ostringstream eherr;
667  eherr << "Cannot obtain the index of " << fieldname;
668  throw InternalErr (__FILE__, __LINE__, eherr.str ());
669  }
670 
671  sdsid = SDselect(sdfileid, sdsindex);
672  if (FAIL == sdsid) {
673  if(true == isgeofile || false == check_pass_fileid_key)
674  SDend(sdfileid);
675  ostringstream eherr;
676  eherr << "Cannot obtain the SDS ID of " << fieldname;
677  throw InternalErr (__FILE__, __LINE__, eherr.str ());
678  }
679 
680  char attrname[H4_MAX_NC_NAME + 1];
681  vector<char> attrbuf;
682  vector<char> attrbuf2;
683 
684  scale_factor_attr_index = SDfindattr(sdsid, "scale_factor");
685  if(scale_factor_attr_index!=FAIL)
686  {
687  intn ret = 0;
688  ret = SDattrinfo(sdsid, scale_factor_attr_index, attrname, &attrtype, &attrcount);
689  if (ret==FAIL)
690  {
691  SDendaccess(sdsid);
692  if(true == isgeofile || false == check_pass_fileid_key)
693  SDend(sdfileid);
694  ostringstream eherr;
695  eherr << "Attribute 'scale_factor' in "
696  << fieldname.c_str () << " cannot be obtained.";
697  throw InternalErr (__FILE__, __LINE__, eherr.str ());
698  }
699 
700  attrbuf.clear();
701  attrbuf.resize(DFKNTsize(attrtype)*attrcount);
702  ret = SDreadattr(sdsid, scale_factor_attr_index, (VOIDP)&attrbuf[0]);
703  if (ret==FAIL)
704  {
705  SDendaccess(sdsid);
706  if(true == isgeofile || false == check_pass_fileid_key)
707  SDend(sdfileid);
708  ostringstream eherr;
709  eherr << "Attribute 'scale_factor' in "
710  << fieldname.c_str () << " cannot be obtained.";
711  throw InternalErr (__FILE__, __LINE__, eherr.str ());
712  }
713 
714  // Appears that the assumption for the datatype of scale_factor
715  // is either float or double
716  // for this type of MODIS files. So far we haven't found any problems.
717  // Maybe this is okay.
718  // KY 2013-12-19
719  switch(attrtype)
720  {
721 #define GET_SCALE_FACTOR_ATTR_VALUE(TYPE, CAST) \
722  case DFNT_##TYPE: \
723  { \
724  CAST tmpvalue = *(CAST*)&attrbuf[0]; \
725  scale = (float)tmpvalue; \
726  } \
727  break;
728  GET_SCALE_FACTOR_ATTR_VALUE(FLOAT32, float);
729  GET_SCALE_FACTOR_ATTR_VALUE(FLOAT64, double);
730  default:
731  throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
732 
733  }
734 
735 #undef GET_SCALE_FACTOR_ATTR_VALUE
736  }
737 
738  add_offset_attr_index = SDfindattr(sdsid, "add_offset");
739  if(add_offset_attr_index!=FAIL)
740  {
741  intn ret = 0;
742  ret = SDattrinfo(sdsid, add_offset_attr_index, attrname, &attrtype, &attrcount);
743  if (ret==FAIL)
744  {
745  SDendaccess(sdsid);
746  if(true == isgeofile || false == check_pass_fileid_key)
747  SDend(sdfileid);
748  ostringstream eherr;
749  eherr << "Attribute 'add_offset' in "
750  << fieldname.c_str () << " cannot be obtained.";
751  throw InternalErr (__FILE__, __LINE__, eherr.str ());
752  }
753  attrbuf.clear();
754  attrbuf.resize(DFKNTsize(attrtype)*attrcount);
755  ret = SDreadattr(sdsid, add_offset_attr_index, (VOIDP)&attrbuf[0]);
756  if (ret==FAIL)
757  {
758  SDendaccess(sdsid);
759  if(true == isgeofile || false == check_pass_fileid_key)
760  SDend(sdfileid);
761  ostringstream eherr;
762  eherr << "Attribute 'add_offset' in "
763  << fieldname.c_str () << " cannot be obtained.";
764  throw InternalErr (__FILE__, __LINE__, eherr.str ());
765  }
766  switch(attrtype)
767  {
768 #define GET_ADD_OFFSET_ATTR_VALUE(TYPE, CAST) \
769  case DFNT_##TYPE: \
770  { \
771  CAST tmpvalue = *(CAST*)&attrbuf[0]; \
772  offset2 = (float)tmpvalue; \
773  } \
774  break;
775  GET_ADD_OFFSET_ATTR_VALUE(FLOAT32, float);
776  GET_ADD_OFFSET_ATTR_VALUE(FLOAT64, double);
777  default:
778  throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
779  }
780 #undef GET_ADD_OFFSET_ATTR_VALUE
781  }
782 
783  attrindex = SDfindattr(sdsid, "_FillValue");
784  if(sotype!=DEFAULT_CF_EQU && attrindex!=FAIL)
785  {
786  intn ret = 0;
787  ret = SDattrinfo(sdsid, attrindex, attrname, &attrtype, &attrcount);
788  if (ret==FAIL)
789  {
790  SDendaccess(sdsid);
791  if(true == isgeofile || false == check_pass_fileid_key)
792  SDend(sdfileid);
793  ostringstream eherr;
794  eherr << "Attribute '_FillValue' in "
795  << fieldname.c_str () << " cannot be obtained.";
796  throw InternalErr (__FILE__, __LINE__, eherr.str ());
797  }
798  attrbuf.clear();
799  attrbuf.resize(DFKNTsize(attrtype)*attrcount);
800  ret = SDreadattr(sdsid, attrindex, (VOIDP)&attrbuf[0]);
801  if (ret==FAIL)
802  {
803  SDendaccess(sdsid);
804  if(true == isgeofile || false == check_pass_fileid_key)
805  SDend(sdfileid);
806  ostringstream eherr;
807  eherr << "Attribute '_FillValue' in "
808  << fieldname.c_str () << " cannot be obtained.";
809  throw InternalErr (__FILE__, __LINE__, eherr.str ());
810  }
811 
812  switch(attrtype)
813  {
814 #define GET_FILLVALUE_ATTR_VALUE(TYPE, CAST) \
815  case DFNT_##TYPE: \
816  { \
817  CAST tmpvalue = *(CAST*)&attrbuf[0]; \
818  fillvalue = (float)tmpvalue; \
819  } \
820  break;
821  GET_FILLVALUE_ATTR_VALUE(INT8, int8);
822  GET_FILLVALUE_ATTR_VALUE(INT16, int16);
823  GET_FILLVALUE_ATTR_VALUE(INT32, int32);
824  GET_FILLVALUE_ATTR_VALUE(UINT8, uint8);
825  GET_FILLVALUE_ATTR_VALUE(UINT16, uint16);
826  GET_FILLVALUE_ATTR_VALUE(UINT32, uint32);
827  // Float and double are not considered. Handle them later.
828  default:
829  ;
830  // throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
831 
832  }
833 #undef GET_FILLVALUE_ATTR_VALUE
834  }
835 
836 #if 0
837 
838  // There is a controversy if we need to apply the valid_range to the data, for the time being comment this out.
839  // KY 2013-12-19
840  float orig_valid_min = 0.;
841  float orig_valid_max = 0.;
842 
843 
844  // Retrieve valid_range,valid_range is normally represented as (valid_min,valid_max)
845  // for non-CF scale and offset rules, the data is always float. So we only
846  // need to change the data type to float.
847  attrindex = SDfindattr(sdsid, "valid_range");
848  if(attrindex!=FAIL)
849  {
850  intn ret;
851  ret = SDattrinfo(sdsid, attrindex, attrname, &attrtype, &attrcount);
852  if (ret==FAIL)
853  {
854  detachfunc(gridid);
855  closefunc(gfid);
856  SDendaccess(sdsid);
857  SDend(sdfileid);
858  ostringstream eherr;
859  eherr << "Attribute '_FillValue' in " << fieldname.c_str () << " cannot be obtained.";
860  throw InternalErr (__FILE__, __LINE__, eherr.str ());
861  }
862  attrbuf.clear();
863  attrbuf.resize(DFKNTsize(attrtype)*attrcount);
864  ret = SDreadattr(sdsid, attrindex, (VOIDP)&attrbuf[0]);
865  if (ret==FAIL)
866  {
867  detachfunc(gridid);
868  closefunc(gfid);
869  SDendaccess(sdsid);
870  SDend(sdfileid);
871  ostringstream eherr;
872  eherr << "Attribute '_FillValue' in " << fieldname.c_str () << " cannot be obtained.";
873  throw InternalErr (__FILE__, __LINE__, eherr.str ());
874  }
875 
876  string attrbuf_str(attrbuf.begin(),attrbuf.end());
877 
878  switch(attrtype) {
879 
880  case DFNT_CHAR:
881  {
882  // We need to treat the attribute data as characters or string.
883  // So find the separator.
884  size_t found = attrbuf_str.find_first_of(",");
885  size_t found_from_end = attrbuf_str.find_last_of(",");
886 
887  if (string::npos == found)
888  throw InternalErr(__FILE__,__LINE__,"should find the separator ,");
889  if (found != found_from_end)
890  throw InternalErr(__FILE__,__LINE__,"Only one separator , should be available.");
891 
892  //istringstream(attrbuf_str.substr(0,found))>> orig_valid_min;
893  //istringstream(attrbuf_str.substr(found+1))>> orig_valid_max;
894 
895  orig_valid_min = atof((attrbuf_str.substr(0,found)).c_str());
896  orig_valid_max = atof((attrbuf_str.substr(found+1)).c_str());
897 
898  }
899  break;
900 
901  case DFNT_INT8:
902  {
903  if (2 == temp_attrcount) {
904  orig_valid_min = (float)attrbuf[0];
905  orig_valid_max = (float)attrbuf[1];
906  }
907  else
908  throw InternalErr(__FILE__,__LINE__,"The number of attribute count should be greater than 1.");
909 
910  }
911  break;
912 
913  case DFNT_UINT8:
914  case DFNT_UCHAR:
915  {
916  if (temp_attrcount != 2)
917  throw InternalErr(__FILE__,__LINE__,"The number of attribute count should be 2 for the DFNT_UINT8 type.");
918 
919  unsigned char* temp_valid_range = (unsigned char *)&attrbuf[0];
920  orig_valid_min = (float)(temp_valid_range[0]);
921  orig_valid_max = (float)(temp_valid_range[1]);
922  }
923  break;
924 
925  case DFNT_INT16:
926  {
927  if (temp_attrcount != 2)
928  throw InternalErr(__FILE__,__LINE__,"The number of attribute count should be 2 for the DFNT_INT16 type.");
929 
930  short* temp_valid_range = (short *)&attrbuf[0];
931  orig_valid_min = (float)(temp_valid_range[0]);
932  orig_valid_max = (float)(temp_valid_range[1]);
933  }
934  break;
935 
936  case DFNT_UINT16:
937  {
938  if (temp_attrcount != 2)
939  throw InternalErr(__FILE__,__LINE__,"The number of attribute count should be 2 for the DFNT_UINT16 type.");
940 
941  unsigned short* temp_valid_range = (unsigned short *)&attrbuf[0];
942  orig_valid_min = (float)(temp_valid_range[0]);
943  orig_valid_max = (float)(temp_valid_range[1]);
944  }
945  break;
946 
947  case DFNT_INT32:
948  {
949  if (temp_attrcount != 2)
950  throw InternalErr(__FILE__,__LINE__,"The number of attribute count should be 2 for the DFNT_INT32 type.");
951 
952  int* temp_valid_range = (int *)&attrbuf[0];
953  orig_valid_min = (float)(temp_valid_range[0]);
954  orig_valid_max = (float)(temp_valid_range[1]);
955  }
956  break;
957 
958  case DFNT_UINT32:
959  {
960  if (temp_attrcount != 2)
961  throw InternalErr(__FILE__,__LINE__,"The number of attribute count should be 2 for the DFNT_UINT32 type.");
962 
963  unsigned int* temp_valid_range = (unsigned int *)&attrbuf[0];
964  orig_valid_min = (float)(temp_valid_range[0]);
965  orig_valid_max = (float)(temp_valid_range[1]);
966  }
967  break;
968 
969  case DFNT_FLOAT32:
970  {
971  if (temp_attrcount != 2)
972  throw InternalErr(__FILE__,__LINE__,"The number of attribute count should be 2 for the DFNT_FLOAT32 type.");
973 
974  float* temp_valid_range = (float *)&attrbuf[0];
975  orig_valid_min = temp_valid_range[0];
976  orig_valid_max = temp_valid_range[1];
977  }
978  break;
979 
980  case DFNT_FLOAT64:
981  {
982  if (temp_attrcount != 2)
983  throw InternalErr(__FILE__,__LINE__,"The number of attribute count should be 2 for the DFNT_FLOAT32 type.");
984  double* temp_valid_range = (double *)&attrbuf[0];
985 
986  // Notice: this approach will lose precision and possibly overflow. Fortunately it is not a problem for MODIS data.
987  // This part of code may not be called. If it is called, mostly the value is within the floating-point range.
988  // KY 2013-01-29
989  orig_valid_min = temp_valid_range[0];
990  orig_valid_max = temp_valid_range[1];
991  }
992  break;
993  default:
994  throw InternalErr(__FILE__,__LINE__,"Unsupported data type.");
995  }
996  }
997 
998 #endif
999 
1000 
1001 #if 0
1002  // For testing only.
1003  //cerr << "scale=" << scale << endl;
1004  //cerr << "offset=" << offset2 << endl;
1005  //cerr << "fillvalue=" << fillvalue << endl;
1006 #endif
1007 
1008  SDendaccess(sdsid);
1009  if(true == isgeofile || false == check_pass_fileid_key)
1010  SDend(sdfileid);
1011  }
1012 
1013  // According to our observations, it seems that MODIS products ALWAYS
1014  // use the "scale" factor to make bigger values smaller.
1015  // So for MODIS_MUL_SCALE products, if the scale of some variable is greater than 1,
1016  // it means that for this variable, the MODIS type for this variable may be MODIS_DIV_SCALE.
1017  // For the similar logic, we may need to change MODIS_DIV_SCALE to
1018  // MODIS_MUL_SCALE and MODIS_EQ_SCALE to MODIS_DIV_SCALE.
1019  // We indeed find such a case. HDF-EOS2 Grid MODIS_Grid_1km_2D of MOD(or MYD)09GA is
1020  // a MODIS_EQ_SCALE.
1021  // However,
1022  // the scale_factor of the variable Range_1 in the MOD09GA product is 25.
1023  // According to our observation, this variable should be MODIS_DIV_SCALE.
1024  // We verify this is true according to MODIS 09 product document
1025  // http://modis-sr.ltdri.org/products/MOD09_UserGuide_v1_3.pdf.
1026  // Since this conclusion is based on our observation, we would like to add
1027  // a BESlog to detect if we find
1028  // the similar cases so that we can verify with the corresponding product
1029  // documents to see if this is true.
1030  //
1031  // More information,
1032  // We just verified with the MOD09 data producer, the scale_factor for Range_1
1033  // and Range_c is 25 but the
1034  // equation is still multiplication instead of division. So we have to make this
1035  // as a special case that we don't want to change the scale and offset settings.
1036  // KY 2014-01-13
1037  // However, since this function only handles swath and we haven't found an outlier
1038  // for a swath case, we still keep the old ways.
1039 
1040 
1041  if (MODIS_EQ_SCALE == sotype || MODIS_MUL_SCALE == sotype) {
1042  if (scale > 1) {
1043  sotype = MODIS_DIV_SCALE;
1044  (*BESLog::TheLog())<< "The field " << fieldname << " scale factor is "
1045  << scale << endl
1046  << " But the original scale factor type is MODIS_MUL_SCALE or MODIS_EQ_SCALE. "
1047  << endl
1048  << " Now change it to MODIS_DIV_SCALE. "<<endl;
1049  }
1050  }
1051 
1052  if (MODIS_DIV_SCALE == sotype) {
1053  if (scale < 1) {
1054  sotype = MODIS_MUL_SCALE;
1055  (*BESLog::TheLog())<< "The field " << fieldname << " scale factor is "
1056  << scale << endl
1057  << " But the original scale factor type is MODIS_DIV_SCALE. "
1058  << endl
1059  << " Now change it to MODIS_MUL_SCALE. "<<endl;
1060  }
1061  }
1062 
1064 // RECALCULATE formula
1065 // if(sotype==MODIS_MUL_SCALE) \
1066 // tmpval[l] = (tmptr[l]-field_offset)*scale; \
1067 // else if(sotype==MODIS_EQ_SCALE) \
1068 // tmpval[l] = tmptr[l]*scale + field_offset; \
1069 // else if(sotype==MODIS_DIV_SCALE) \
1070 // tmpval[l] = (tmptr[l]-field_offset)/scale; \
1071 
1073 
1074 #define RECALCULATE(CAST, DODS_CAST, VAL) \
1075 { \
1076  bool change_data_value = false; \
1077  if(sotype!=DEFAULT_CF_EQU) \
1078  { \
1079  if(scale_factor_attr_index!=FAIL && !(scale==1 && offset2==0)) \
1080  { \
1081  vector<float>tmpval; \
1082  tmpval.resize(nelms); \
1083  CAST tmptr = (CAST)VAL; \
1084  for(int l=0; l<nelms; l++) \
1085  tmpval[l] = (float)tmptr[l]; \
1086  float temp_scale = scale; \
1087  float temp_offset = offset2; \
1088  if(sotype==MODIS_MUL_SCALE) \
1089  temp_offset = -1. *offset2*temp_scale;\
1090  else if (sotype==MODIS_DIV_SCALE) \
1091  {\
1092  temp_scale = 1/scale; \
1093  temp_offset = -1. *temp_scale *offset2;\
1094  }\
1095  for(int l=0; l<nelms; l++) \
1096  if(attrindex!=FAIL && ((float)tmptr[l])!=fillvalue) \
1097  tmpval[l] = tmptr[l]*temp_scale + temp_offset; \
1098  change_data_value = true; \
1099  set_value((dods_float32 *)&tmpval[0], nelms); \
1100  } \
1101  } \
1102  if(!change_data_value) \
1103  { \
1104  set_value ((DODS_CAST)VAL, nelms); \
1105  } \
1106 }
1107 
1108  // tmp_rank and tmp_dimlist are two dummy variables that are
1109  // only used when calling fieldinfo.
1110  int32 tmp_rank = 0;
1111  char tmp_dimlist[1024];
1112 
1113  // field dimension sizes
1114  int32 tmp_dims[rank];
1115 
1116  // field data type
1117  int32 field_dtype = 0;
1118 
1119  // returned value of HDF4 and HDF-EOS2 APIs
1120  intn r = 0;
1121 
1122  // Obtain the field info. We mainly need the datatype information
1123  // to allocate the buffer to store the data
1124  r = fieldinfofunc (swathid, const_cast < char *>(fieldname.c_str ()),
1125  &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1126  if (r != 0) {
1127  ostringstream eherr;
1128  eherr << "Field " << fieldname.c_str ()
1129  << " information cannot be obtained.";
1130  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1131  }
1132 
1133 
1134  // int32 majordimsize, minordimsize;
1135  vector<int32> newdims;
1136  newdims.resize(rank);
1137 
1138  // Loop through the data type.
1139  switch (field_dtype) {
1140 
1141  case DFNT_INT8:
1142  {
1143  // Obtaining the total value and interpolating the data
1144  // according to dimension map
1145  vector < int8 > total_val8;
1146  r = GetFieldValue (swathid, fieldname, dimmaps, total_val8, newdims);
1147  if (r != 0) {
1148  ostringstream eherr;
1149  eherr << "field " << fieldname.c_str ()
1150  << "cannot be read.";
1151  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1152  }
1153 
1154  check_num_elems_constraint(nelms,newdims);
1155 
1156  vector<int8>val8;
1157  val8.resize(nelms);
1158 
1159  FieldSubset (&val8[0], newdims, &total_val8[0],
1160  &offset32[0], &count32[0], &step32[0]);
1161 
1162 #ifndef SIGNED_BYTE_TO_INT32
1163  RECALCULATE(int8*, dods_byte*, &val8[0]);
1164 #else
1165  vector<int32>newval;
1166  newval.resize(nelms);
1167 
1168  for (int counter = 0; counter < nelms; counter++)
1169  newval[counter] = (int32) (val8[counter]);
1170 
1171  RECALCULATE(int32*, dods_int32*, &newval[0]);
1172 #endif
1173  }
1174  break;
1175  case DFNT_UINT8:
1176  case DFNT_UCHAR8:
1177  {
1178  // Obtaining the total value and interpolating the data
1179  // according to dimension map
1180  vector < uint8 > total_val_u8;
1181  r = GetFieldValue (swathid, fieldname, dimmaps, total_val_u8, newdims);
1182  if (r != 0) {
1183  ostringstream eherr;
1184  eherr << "field " << fieldname.c_str () << "cannot be read.";
1185  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1186  }
1187 
1188  check_num_elems_constraint(nelms,newdims);
1189  vector<uint8>val_u8;
1190  val_u8.resize(nelms);
1191 
1192  FieldSubset (&val_u8[0], newdims, &total_val_u8[0],
1193  &offset32[0], &count32[0], &step32[0]);
1194  RECALCULATE(uint8*, dods_byte*, &val_u8[0]);
1195  }
1196  break;
1197  case DFNT_INT16:
1198  {
1199  // Obtaining the total value and interpolating the data
1200  // according to dimension map
1201  vector < int16 > total_val16;
1202  r = GetFieldValue (swathid, fieldname, dimmaps, total_val16, newdims);
1203  if (r != 0) {
1204  ostringstream eherr;
1205  eherr << "field " << fieldname.c_str () << "cannot be read.";
1206  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1207  }
1208 
1209  check_num_elems_constraint(nelms,newdims);
1210  vector<int16>val16;
1211  val16.resize(nelms);
1212 
1213  FieldSubset (&val16[0], newdims, &total_val16[0],
1214  &offset32[0], &count32[0], &step32[0]);
1215 
1216  RECALCULATE(int16*, dods_int16*, &val16[0]);
1217  }
1218  break;
1219  case DFNT_UINT16:
1220  {
1221  // Obtaining the total value and interpolating the data
1222  // according to dimension map
1223  vector < uint16 > total_val_u16;
1224  r = GetFieldValue (swathid, fieldname, dimmaps, total_val_u16, newdims);
1225  if (r != 0) {
1226  ostringstream eherr;
1227 
1228  eherr << "field " << fieldname.c_str () << "cannot be read.";
1229  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1230  }
1231 
1232  check_num_elems_constraint(nelms,newdims);
1233  vector<uint16>val_u16;
1234  val_u16.resize(nelms);
1235 
1236  FieldSubset (&val_u16[0], newdims, &total_val_u16[0],
1237  &offset32[0], &count32[0], &step32[0]);
1238  RECALCULATE(uint16*, dods_uint16*, &val_u16[0]);
1239 
1240  }
1241  break;
1242  case DFNT_INT32:
1243  {
1244  // Obtaining the total value and interpolating the data
1245  // according to dimension map
1246  vector < int32 > total_val32;
1247  r = GetFieldValue (swathid, fieldname, dimmaps, total_val32, newdims);
1248  if (r != 0) {
1249  ostringstream eherr;
1250 
1251  eherr << "field " << fieldname.c_str () << "cannot be read.";
1252  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1253  }
1254 
1255  check_num_elems_constraint(nelms,newdims);
1256  vector<int32> val32;
1257  val32.resize(nelms);
1258 
1259  FieldSubset (&val32[0], newdims, &total_val32[0],
1260  &offset32[0], &count32[0], &step32[0]);
1261 
1262  RECALCULATE(int32*, dods_int32*, &val32[0]);
1263  }
1264  break;
1265  case DFNT_UINT32:
1266  {
1267  // Obtaining the total value and interpolating the data
1268  // according to dimension map
1269  // Notice the total_val_u32 is allocated inside the GetFieldValue.
1270  vector < uint32 > total_val_u32;
1271  r = GetFieldValue (swathid, fieldname, dimmaps, total_val_u32, newdims);
1272  if (r != 0) {
1273  ostringstream eherr;
1274  eherr << "field " << fieldname.c_str () << "cannot be read.";
1275  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1276  }
1277 
1278  check_num_elems_constraint(nelms,newdims);
1279  vector<uint32>val_u32;
1280  val_u32.resize(nelms);
1281 
1282  FieldSubset (&val_u32[0], newdims, &total_val_u32[0],
1283  &offset32[0], &count32[0], &step32[0]);
1284  RECALCULATE(uint32*, dods_uint32*, &val_u32[0]);
1285 
1286  }
1287  break;
1288  case DFNT_FLOAT32:
1289  {
1290  // Obtaining the total value and interpolating the data
1291  // according to dimension map
1292  vector < float32 > total_val_f32;
1293  r = GetFieldValue (swathid, fieldname, dimmaps, total_val_f32, newdims);
1294  if (r != 0) {
1295  ostringstream eherr;
1296  eherr << "field " << fieldname.c_str () << "cannot be read.";
1297  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1298  }
1299 
1300  check_num_elems_constraint(nelms,newdims);
1301  vector<float32>val_f32;
1302  val_f32.resize(nelms);
1303 
1304  FieldSubset (&val_f32[0], newdims, &total_val_f32[0],
1305  &offset32[0], &count32[0], &step32[0]);
1306  RECALCULATE(float32*, dods_float32*, &val_f32[0]);
1307  }
1308  break;
1309  case DFNT_FLOAT64:
1310  {
1311  // Obtaining the total value and interpolating the data according to dimension map
1312  vector < float64 > total_val_f64;
1313  r = GetFieldValue (swathid, fieldname, dimmaps, total_val_f64, newdims);
1314  if (r != 0) {
1315  ostringstream eherr;
1316  eherr << "field " << fieldname.c_str () << "cannot be read.";
1317  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1318  }
1319 
1320  check_num_elems_constraint(nelms,newdims);
1321  vector<float64>val_f64;
1322  val_f64.resize(nelms);
1323  FieldSubset (&val_f64[0], newdims, &total_val_f64[0],
1324  &offset32[0], &count32[0], &step32[0]);
1325  RECALCULATE(float64*, dods_float64*, &val_f64[0]);
1326 
1327  }
1328  break;
1329  default:
1330  {
1331  throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
1332  }
1333  }
1334 
1335  return 0;
1336 }
1337 
1338 int
1339 HDFEOS2ArraySwathDimMapField::write_dap_data_disable_scale_comp(int32 swathid,
1340  int nelms,
1341  vector<int32>& offset32,
1342  vector<int32>& count32,
1343  vector<int32>& step32) {
1344 
1345  // Define function pointers to handle both grid and swath
1346  intn (*fieldinfofunc) (int32, char *, int32 *, int32 *, int32 *, char *);
1347 
1348  fieldinfofunc = SWfieldinfo;
1349 
1350 
1351  // tmp_rank and tmp_dimlist are two dummy variables
1352  // that are only used when calling fieldinfo.
1353  int32 tmp_rank = 0;
1354  char tmp_dimlist[1024];
1355 
1356  // field dimension sizes
1357  int32 tmp_dims[rank];
1358 
1359  // field data type
1360  int32 field_dtype = 0;
1361 
1362  // returned value of HDF4 and HDF-EOS2 APIs
1363  intn r = 0;
1364 
1365  // Obtain the field info. We mainly need the datatype information
1366  // to allocate the buffer to store the data
1367  r = fieldinfofunc (swathid, const_cast < char *>(fieldname.c_str ()),
1368  &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1369  if (r != 0) {
1370  ostringstream eherr;
1371  eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
1372  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1373  }
1374 
1375 
1376  // int32 majordimsize, minordimsize;
1377  vector<int32> newdims;
1378  newdims.resize(rank);
1379 
1380  // Loop through the data type.
1381  switch (field_dtype) {
1382 
1383  case DFNT_INT8:
1384  {
1385  // Obtaining the total value and interpolating the data according to dimension map
1386  vector < int8 > total_val8;
1387  r = GetFieldValue (swathid, fieldname, dimmaps, total_val8, newdims);
1388  if (r != 0) {
1389  ostringstream eherr;
1390  eherr << "field " << fieldname.c_str () << "cannot be read.";
1391  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1392  }
1393 
1394  check_num_elems_constraint(nelms,newdims);
1395 
1396  vector<int8>val8;
1397  val8.resize(nelms);
1398  FieldSubset (&val8[0], newdims, &total_val8[0],
1399  &offset32[0], &count32[0], &step32[0]);
1400 
1401 
1402 #ifndef SIGNED_BYTE_TO_INT32
1403  set_value((dods_byte*)&val8[0],nelms);
1404 #else
1405  vector<int32>newval;
1406  newval.resize(nelms);
1407 
1408  for (int counter = 0; counter < nelms; counter++)
1409  newval[counter] = (int32) (val8[counter]);
1410 
1411  set_value ((dods_int32 *) &newval[0], nelms);
1412 #endif
1413  }
1414  break;
1415  case DFNT_UINT8:
1416  case DFNT_UCHAR8:
1417  {
1418  // Obtaining the total value and interpolating the data according to dimension map
1419  vector < uint8 > total_val_u8;
1420  r = GetFieldValue (swathid, fieldname, dimmaps, total_val_u8, newdims);
1421  if (r != 0) {
1422  ostringstream eherr;
1423  eherr << "field " << fieldname.c_str () << "cannot be read.";
1424  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1425  }
1426 
1427  check_num_elems_constraint(nelms,newdims);
1428  vector<uint8>val_u8;
1429  val_u8.resize(nelms);
1430 
1431  FieldSubset (&val_u8[0], newdims, &total_val_u8[0],
1432  &offset32[0], &count32[0], &step32[0]);
1433  set_value ((dods_byte *) &val_u8[0], nelms);
1434  }
1435  break;
1436  case DFNT_INT16:
1437  {
1438  // Obtaining the total value and interpolating the data according to dimension map
1439  vector < int16 > total_val16;
1440  r = GetFieldValue (swathid, fieldname, dimmaps, total_val16, newdims);
1441  if (r != 0) {
1442  ostringstream eherr;
1443  eherr << "field " << fieldname.c_str () << "cannot be read.";
1444  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1445  }
1446 
1447  check_num_elems_constraint(nelms,newdims);
1448  vector<int16>val16;
1449  val16.resize(nelms);
1450 
1451  FieldSubset (&val16[0], newdims, &total_val16[0],
1452  &offset32[0], &count32[0], &step32[0]);
1453 
1454  set_value ((dods_int16 *) &val16[0], nelms);
1455  }
1456  break;
1457  case DFNT_UINT16:
1458  {
1459  // Obtaining the total value and interpolating the data according to dimension map
1460  vector < uint16 > total_val_u16;
1461  r = GetFieldValue (swathid, fieldname, dimmaps, total_val_u16, newdims);
1462  if (r != 0) {
1463  ostringstream eherr;
1464  eherr << "field " << fieldname.c_str () << "cannot be read.";
1465  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1466  }
1467 
1468  check_num_elems_constraint(nelms,newdims);
1469  vector<uint16>val_u16;
1470  val_u16.resize(nelms);
1471 
1472  FieldSubset (&val_u16[0], newdims, &total_val_u16[0],
1473  &offset32[0], &count32[0], &step32[0]);
1474  set_value ((dods_uint16 *) &val_u16[0], nelms);
1475 
1476  }
1477  break;
1478  case DFNT_INT32:
1479  {
1480  // Obtaining the total value and interpolating the data according to dimension map
1481  vector < int32 > total_val32;
1482  r = GetFieldValue (swathid, fieldname, dimmaps, total_val32, newdims);
1483  if (r != 0) {
1484  ostringstream eherr;
1485 
1486  eherr << "field " << fieldname.c_str () << "cannot be read.";
1487  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1488  }
1489 
1490  check_num_elems_constraint(nelms,newdims);
1491  vector<int32> val32;
1492  val32.resize(nelms);
1493 
1494  FieldSubset (&val32[0], newdims, &total_val32[0],
1495  &offset32[0], &count32[0], &step32[0]);
1496  set_value ((dods_int32 *) &val32[0], nelms);
1497  }
1498  break;
1499  case DFNT_UINT32:
1500  {
1501  // Obtaining the total value and interpolating the data according to dimension map
1502  // Notice the total_val_u32 is allocated inside the GetFieldValue.
1503  vector < uint32 > total_val_u32;
1504  r = GetFieldValue (swathid, fieldname, dimmaps, total_val_u32, newdims);
1505  if (r != 0) {
1506  ostringstream eherr;
1507 
1508  eherr << "field " << fieldname.c_str () << "cannot be read.";
1509  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1510  }
1511 
1512  check_num_elems_constraint(nelms,newdims);
1513  vector<uint32>val_u32;
1514  val_u32.resize(nelms);
1515 
1516  FieldSubset (&val_u32[0], newdims, &total_val_u32[0],
1517  &offset32[0], &count32[0], &step32[0]);
1518  set_value ((dods_uint32 *) &val_u32[0], nelms);
1519 
1520  }
1521  break;
1522  case DFNT_FLOAT32:
1523  {
1524  // Obtaining the total value and interpolating the data according to dimension map
1525  vector < float32 > total_val_f32;
1526  r = GetFieldValue (swathid, fieldname, dimmaps, total_val_f32, newdims);
1527  if (r != 0) {
1528  ostringstream eherr;
1529 
1530  eherr << "field " << fieldname.c_str () << "cannot be read.";
1531  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1532  }
1533 
1534  check_num_elems_constraint(nelms,newdims);
1535  vector<float32>val_f32;
1536  val_f32.resize(nelms);
1537 
1538  FieldSubset (&val_f32[0], newdims, &total_val_f32[0],
1539  &offset32[0], &count32[0], &step32[0]);
1540 
1541  set_value ((dods_float32 *) &val_f32[0], nelms);
1542  }
1543  break;
1544  case DFNT_FLOAT64:
1545  {
1546  // Obtaining the total value and interpolating the data according to dimension map
1547  vector < float64 > total_val_f64;
1548  r = GetFieldValue (swathid, fieldname, dimmaps, total_val_f64, newdims);
1549  if (r != 0) {
1550  ostringstream eherr;
1551 
1552  eherr << "field " << fieldname.c_str () << "cannot be read.";
1553  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1554  }
1555 
1556  check_num_elems_constraint(nelms,newdims);
1557  vector<float64>val_f64;
1558  val_f64.resize(nelms);
1559  FieldSubset (&val_f64[0], newdims, &total_val_f64[0],
1560  &offset32[0], &count32[0], &step32[0]);
1561  set_value ((dods_float64 *) &val_f64[0], nelms);
1562 
1563  }
1564  break;
1565  default:
1566  {
1567  throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
1568  }
1569  }
1570 
1571  return 0;
1573 //#endif
1574 }
1575 
1576 void HDFEOS2ArraySwathDimMapField::close_fileid(const int32 swfileid, const int32 sdfileid) {
1577 
1578 #if 0
1579  string check_pass_fileid_key_str="H4.EnablePassFileID";
1580  bool check_pass_fileid_key = false;
1581  check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
1582 #endif
1583 
1584  //if(true == isgeofile || false == check_pass_fileid_key) {
1585  if(true == isgeofile || false == HDF4RequestHandler::get_pass_fileid()) {
1586 
1587  if(sdfileid != -1)
1588  SDend(sdfileid);
1589 
1590  if(swfileid != -1)
1591  SWclose(swfileid);
1592 
1593  }
1594 
1595 }
1596 
1597 bool HDFEOS2ArraySwathDimMapField::check_num_elems_constraint(const int num_elems,
1598  const vector<int32>&newdims) {
1599 
1600  int total_dim_size = 1;
1601  for (int i =0;i<rank;i++)
1602  total_dim_size*=newdims[i];
1603 
1604  if(total_dim_size < num_elems) {
1605  ostringstream eherr;
1606  eherr << "The total number of elements for the array " << total_dim_size
1607  << "is less than the user-requested number of elements " << num_elems;
1608  throw InternalErr (__FILE__, __LINE__, eherr.str ());
1609  }
1610 
1611  return false;
1612 
1613 }
1614 #endif
close_fileid
void close_fileid(hid_t fid)
Definition: h5get.cc:414
HDFCFUtil::Split
static void Split(const char *s, int len, char sep, std::vector< std::string > &names)
Definition: HDFCFUtil.cc:77
dimmap_entry
Definition: HDFCFUtil.h:51
Error