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