bes Updated for version 3.20.10
HDFEOS2ArraySwathDimMapField.cc
1
2
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 <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;
32bool
33HDFEOS2ArraySwathDimMapField::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
193else {
194for(int i = 0; i <dimmaps.size();i++) {
195cerr<<"dimmap geodim is: "<<dimmaps[i].geodim <<endl;
196cerr<<"dimmap datadim is: "<<dimmaps[i].datadim <<endl;
197cerr<<"offset is: "<<dimmaps[i].offset <<endl;
198cerr<<"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.
284int
285HDFEOS2ArraySwathDimMapField::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.
326template < class T > int
327HDFEOS2ArraySwathDimMapField::
328GetFieldValue (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.
393template < class T > int
394HDFEOS2ArraySwathDimMapField::_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
492template < class T >
493bool 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
515template < class T >
516bool 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
534template < class T >
535bool 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
580template < class T >
581bool 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
633int
634HDFEOS2ArraySwathDimMapField::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
1360int
1361HDFEOS2ArraySwathDimMapField::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
1598void 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
1619bool 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:433
static void Split(const char *s, int len, char sep, std::vector< std::string > &names)
Definition: HDFCFUtil.cc:80
Definition: HDFCFUtil.h:52