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