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