bes Updated for version 3.20.10
HDFEOS2Array_RealField.cc
1
2// This file is part of the hdf4 data handler for the OPeNDAP data server.
3// It retrieves the real field values.
4// Authors: MuQun Yang <myang6@hdfgroup.org> Eunsoo Seo
5// Copyright (c) 2010-2012 The HDF Group
7#ifdef USE_HDFEOS2_LIB
8#include "config.h"
9#include "config_hdf.h"
10
11#include <iostream>
12#include <sstream>
13#include <cassert>
14#include <libdap/debug.h>
15#include <libdap/InternalErr.h>
16#include <BESDebug.h>
17#include <BESLog.h>
18
19#include "HDFCFUtil.h"
20#include "HDFEOS2Array_RealField.h"
21#include "dodsutil.h"
22#include "HDF4RequestHandler.h"
23
24using namespace std;
25using namespace libdap;
26
27#define SIGNED_BYTE_TO_INT32 1
28
29bool
30HDFEOS2Array_RealField::read ()
31{
32
33 BESDEBUG("h4","Coming to HDFEOS2_Array_RealField read "<<endl);
34 if(length() == 0)
35 return true;
36
37#if 0
38 string check_pass_fileid_key_str="H4.EnablePassFileID";
39 bool check_pass_fileid_key = false;
40 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
41#endif
42
43 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
44
45 // Declare offset, count and step
46 vector<int>offset;
47 offset.resize(rank);
48 vector<int>count;
49 count.resize(rank);
50 vector<int>step;
51 step.resize(rank);
52
53 // Obtain offset,step and count from the client expression constraint
54 int nelms = 0;
55 nelms = format_constraint (&offset[0], &step[0], &count[0]);
56
57 // Just declare offset,count and step in the int32 type.
58 vector<int32>offset32;
59 offset32.resize(rank);
60 vector<int32>count32;
61 count32.resize(rank);
62 vector<int32>step32;
63 step32.resize(rank);
64
65 // Just obtain the offset,count and step in the datatype of int32.
66 for (int i = 0; i < rank; i++) {
67 offset32[i] = (int32) offset[i];
68 count32[i] = (int32) count[i];
69 step32[i] = (int32) step[i];
70 }
71
72 // Define function pointers to handle both grid and swath
73 int32 (*openfunc) (char *, intn);
74 intn (*closefunc) (int32);
75 int32 (*attachfunc) (int32, char *);
76 intn (*detachfunc) (int32);
77 intn (*fieldinfofunc) (int32, char *, int32 *, int32 *, int32 *, char *);
78
79 string datasetname;
80 if (swathname == "") {
81 openfunc = GDopen;
82 closefunc = GDclose;
83 attachfunc = GDattach;
84 detachfunc = GDdetach;
85 fieldinfofunc = GDfieldinfo;
86 datasetname = gridname;
87 }
88 else if (gridname == "") {
89 openfunc = SWopen;
90 closefunc = SWclose;
91 attachfunc = SWattach;
92 detachfunc = SWdetach;
93 fieldinfofunc = SWfieldinfo;
94 datasetname = swathname;
95 }
96 else
97 throw InternalErr (__FILE__, __LINE__, "It should be either grid or swath.");
98
99 // Note gfid and gridid represent either swath or grid.
100 int32 gfid = 0;
101 int32 gridid = 0;
102
103 if (true == isgeofile || false == check_pass_fileid_key) {
104
105 // Obtain the EOS object ID(either grid or swath)
106 gfid = openfunc (const_cast < char *>(filename.c_str ()), DFACC_READ);
107 if (gfid < 0) {
108 ostringstream eherr;
109 eherr << "File " << filename.c_str () << " cannot be open.";
110 throw InternalErr (__FILE__, __LINE__, eherr.str ());
111 }
112 }
113 else
114 gfid = gsfd;
115
116 // Attach the EOS object ID
117 gridid = attachfunc (gfid, const_cast < char *>(datasetname.c_str ()));
118 if (gridid < 0) {
119 close_fileid(gfid,-1);
120 ostringstream eherr;
121 eherr << "Grid/Swath " << datasetname.c_str () << " cannot be attached.";
122 throw InternalErr (__FILE__, __LINE__, eherr.str ());
123 }
124
125#if 0
126 string check_disable_scale_comp_key = "H4.DisableScaleOffsetComp";
127 bool turn_on_disable_scale_comp_key= false;
128 turn_on_disable_scale_comp_key = HDFCFUtil::check_beskeys(check_disable_scale_comp_key);
129#endif
130
131
132 bool is_modis_l1b = false;
133 if("MODIS_SWATH_Type_L1B" == swathname)
134 is_modis_l1b = true;
135
136 bool is_modis_vip = false;
137 if ("VIP_CMG_GRID" == gridname)
138 is_modis_vip = true;
139
140 bool field_is_vdata = false;
141
142 // HDF-EOS2 swath maps 1-D field as vdata. So we need to check if this field is vdata or SDS.
143 // Essentially we only call SDS attribute routines to retrieve MODIS scale,offset and
144 // fillvalue attributes since we don't
145 // find 1-D MODIS field has scale,offset and fillvalue attributes. We may need to visit
146 // this again in the future to see if we also need to support the handling of
147 // scale,offset,fillvalue via vdata routines. KY 2013-07-15
148 if (""==gridname) {
149
150 int32 tmp_rank = 0;
151 char tmp_dimlist[1024];
152 int32 tmp_dims[rank];
153 int32 field_dtype = 0;
154 intn r = 0;
155
156 r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
157 &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
158 if (r != 0) {
159 detachfunc(gridid);
160 close_fileid(gfid,-1);
161 ostringstream eherr;
162
163 eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
164 throw InternalErr (__FILE__, __LINE__, eherr.str ());
165 }
166
167 if (1 == tmp_rank)
168 field_is_vdata = true;
169 }
170
171
172 bool has_Key_attr = false;
173
174 if (false == field_is_vdata) {
175
176 // Obtain attribute values.
177 int32 sdfileid = -1;
178
179 if (true == isgeofile || false == check_pass_fileid_key) {
180
181 sdfileid = SDstart(const_cast < char *>(filename.c_str ()), DFACC_READ);
182
183 if (FAIL == sdfileid) {
184 detachfunc(gridid);
185 close_fileid(gfid,-1);
186 ostringstream eherr;
187 eherr << "Cannot Start the SD interface for the file " << filename <<endl;
188 }
189 }
190 else
191 sdfileid = sdfd;
192
193 int32 sdsindex = -1;
194 int32 sdsid = -1;
195 sdsindex = SDnametoindex(sdfileid, fieldname.c_str());
196 if (FAIL == sdsindex) {
197 detachfunc(gridid);
198 close_fileid(gfid,sdfileid);
199 ostringstream eherr;
200 eherr << "Cannot obtain the index of " << fieldname;
201 throw InternalErr (__FILE__, __LINE__, eherr.str ());
202 }
203
204 sdsid = SDselect(sdfileid, sdsindex);
205 if (FAIL == sdsid) {
206 detachfunc(gridid);
207 close_fileid(gfid,sdfileid);
208 ostringstream eherr;
209 eherr << "Cannot obtain the SDS ID of " << fieldname;
210 throw InternalErr (__FILE__, __LINE__, eherr.str ());
211 }
212
213 // Here we cannot check if SDfindattr fails since even SDfindattr fails it doesn't mean
214 // errors happen. If no such attribute can be found, SDfindattr still returns FAIL.
215 // The correct way is to use SDgetinfo and SDattrinfo to check if attributes
216 // "radiance_scales" etc exist.
217 // For the time being, I won't do this, due to the performance reason and code simplicity and also the
218 // very small chance of real FAIL for SDfindattr.
219 if(SDfindattr(sdsid, "Key")!=FAIL)
220 has_Key_attr = true;
221
222 // Close the interfaces
223 SDendaccess(sdsid);
224 if (true == isgeofile || false == check_pass_fileid_key)
225 SDend(sdfileid);
226 }
227
228 // USE a try-catch block to release the resources.
229 try {
230 if((false == is_modis_l1b) && (false == is_modis_vip)
231 &&(false == has_Key_attr) && (true == HDF4RequestHandler::get_disable_scaleoffset_comp()))
232 write_dap_data_disable_scale_comp(gridid,nelms,&offset32[0],&count32[0],&step32[0]);
233 else
234 write_dap_data_scale_comp(gridid,nelms,offset32,count32,step32);
235 }
236 catch(...) {
237 detachfunc(gridid);
238 close_fileid(gfid,-1);
239 throw;
240 }
241
242 int32 r = -1;
243 r = detachfunc (gridid);
244 if (r != 0) {
245 close_fileid(gfid,-1);
246 ostringstream eherr;
247 eherr << "Grid/Swath " << datasetname.c_str () << " cannot be detached.";
248 throw InternalErr (__FILE__, __LINE__, eherr.str ());
249 }
250
251
252 if(true == isgeofile || false == check_pass_fileid_key) {
253 r = closefunc (gfid);
254 if (r != 0) {
255 ostringstream eherr;
256 eherr << "Grid/Swath " << filename.c_str () << " cannot be closed.";
257 throw InternalErr (__FILE__, __LINE__, eherr.str ());
258 }
259 }
260
261 return false;
262}
263
264int
265HDFEOS2Array_RealField::write_dap_data_scale_comp(int32 gridid,
266 int nelms,
267 vector<int32>& offset32,
268 vector<int32>& count32,
269 vector<int32>& step32) {
270
271
272 BESDEBUG("h4",
273 "coming to HDFEOS2Array_RealField write_dap_data_scale_comp "
274 <<endl);
275
276#if 0
277 string check_pass_fileid_key_str="H4.EnablePassFileID";
278 bool check_pass_fileid_key = false;
279 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
280#endif
281 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
282
283 // Define function pointers to handle both grid and swath
284 intn (*fieldinfofunc) (int32, char *, int32 *, int32 *, int32 *, char *);
285 intn (*readfieldfunc) (int32, char *, int32 *, int32 *, int32 *, void *);
286
287
288 if (swathname == "") {
289 fieldinfofunc = GDfieldinfo;
290 readfieldfunc = GDreadfield;
291 }
292 else if (gridname == "") {
293 fieldinfofunc = SWfieldinfo;
294 readfieldfunc = SWreadfield;
295 }
296 else
297 throw InternalErr (__FILE__, __LINE__, "It should be either grid or swath.");
298
299 // tmp_rank and tmp_dimlist are two dummy variables that are only used when calling fieldinfo.
300 int32 tmp_rank = 0;
301 char tmp_dimlist[1024];
302
303 // field dimension sizes
304 int32 tmp_dims[rank];
305
306 // field data type
307 int32 field_dtype = 0;
308
309 // returned value of HDF4 and HDF-EOS2 APIs
310 intn r = 0;
311
312 // Obtain the field info. We mainly need the datatype information
313 // to allocate the buffer to store the data
314 r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
315 &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
316 if (r != 0) {
317 ostringstream eherr;
318
319 eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
320 throw InternalErr (__FILE__, __LINE__, eherr.str ());
321 }
322
323
324 // The following chunk of code until switch(field_dtype) handles MODIS level 1B,
325 // MOD29E1D Key and VIP products. The reason to keep the code this way is due to
326 // use of RECALCULATE macro. It is too much work to change it now. KY 2013-12-17
327 // MODIS level 1B reflectance and radiance fields have scale/offset arrays rather
328 // than one scale/offset value.
329 // So we need to handle these fields specially.
330 float *reflectance_offsets =NULL;
331 float *reflectance_scales =NULL;
332 float *radiance_offsets =NULL;
333 float *radiance_scales =NULL;
334
335 // Attribute datatype, reused for several attributes
336 int32 attr_dtype = 0;
337
338 // Number of elements for an attribute, reused
339 int32 temp_attrcount = 0;
340
341 // Number of elements in an attribute
342 int32 num_eles_of_an_attr = 0;
343
344 // Attribute(radiance_scales, reflectance_scales) index
345 int32 cf_modl1b_rr_attrindex = 0;
346
347 // Attribute (radiance_offsets) index
348 int32 cf_modl1b_rr_attrindex2 = 0;
349
350 // Attribute valid_range index
351 int32 cf_vr_attrindex = 0;
352
353 // Attribute fill value index
354 int32 cf_fv_attrindex = 0;
355
356 // Scale factor attribute index
357 int32 scale_factor_attr_index = 0;
358
359 // Add offset attribute index
360 int32 add_offset_attr_index = 0;
361
362 // Initialize scale
363 float scale = 1;
364
365 // Intialize field_offset
366 float field_offset = 0;
367
368 // Initialize fillvalue
369 float fillvalue = 0;
370
371 // Initialize the original valid_min
372 float orig_valid_min = 0;
373
374 // Initialize the original valid_max
375 float orig_valid_max = 0;
376
377 // Some NSIDC products use the "Key" attribute to identify
378 // the discrete valid values(land, cloud etc).
379 // Since the valid_range attribute in these products may treat values
380 // identified by the Key attribute as invalid,
381 // we need to handle them in a special way.So set a flag here.
382 bool has_Key_attr = false;
383
384 int32 sdfileid = -1;
385 if(sotype!=DEFAULT_CF_EQU) {
386
387 bool field_is_vdata = false;
388
389 // HDF-EOS2 swath maps 1-D field as vdata. So we need to check
390 // if this field is vdata or SDS.
391 // Essentially we only call SDS attribute routines to retrieve MODIS scale,
392 // offset and fillvalue
393 // attributes since we don't find 1-D MODIS field has scale,offset and
394 // fillvalue attributes.
395 // We may need to visit this again in the future to see
396 // if we also need to support the handling of scale,offset,fillvalue via
397 // vdata routines. KY 2013-07-15
398 if (""==gridname) {
399
400 r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
401 &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
402 if (r != 0) {
403 ostringstream eherr;
404 eherr << "Field " << fieldname.c_str ()
405 << " information cannot be obtained.";
406 throw InternalErr (__FILE__, __LINE__, eherr.str ());
407 }
408
409 if (1 == tmp_rank)
410 field_is_vdata = true;
411 }
413#if 0
414r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
415 &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
416if (r != 0) {
417 ostringstream eherr;
418
419 eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
420 throw InternalErr (__FILE__, __LINE__, eherr.str ());
421}
422
423cerr<<"tmp_rank is "<<tmp_rank <<endl;
424#endif
425
426
427 // For swath, we don't see any MODIS 1-D fields that have scale,offset
428 // and fill value attributes that need to be changed.
429 // So now we don't need to handle the vdata handling.
430 // Another reason is the possible change of the implementation
431 // of the SDS attribute handlings. That may be too costly.
432 // KY 2012-07-31
433
434 if( false == field_is_vdata) {
435
436 char attrname[H4_MAX_NC_NAME + 1];
437 vector<char> attrbuf;
438
439 // Obtain attribute values.
440 if(false == isgeofile || false == check_pass_fileid_key)
441 sdfileid = sdfd;
442 else {
443 sdfileid = SDstart(const_cast < char *>(filename.c_str ()), DFACC_READ);
444 if (FAIL == sdfileid) {
445 ostringstream eherr;
446 eherr << "Cannot Start the SD interface for the file "
447 << filename;
448 throw InternalErr (__FILE__, __LINE__, eherr.str ());
449 }
450 }
451
452 int32 sdsindex = -1;
453 int32 sdsid;
454 sdsindex = SDnametoindex(sdfileid, fieldname.c_str());
455 if (FAIL == sdsindex) {
456 if(true == isgeofile || false == check_pass_fileid_key)
457 SDend(sdfileid);
458 ostringstream eherr;
459 eherr << "Cannot obtain the index of " << fieldname;
460 throw InternalErr (__FILE__, __LINE__, eherr.str ());
461 }
462
463 sdsid = SDselect(sdfileid, sdsindex);
464 if (FAIL == sdsid) {
465 if (true == isgeofile || false == check_pass_fileid_key)
466 SDend(sdfileid);
467 ostringstream eherr;
468 eherr << "Cannot obtain the SDS ID of " << fieldname;
469 throw InternalErr (__FILE__, __LINE__, eherr.str ());
470 }
471
472#if 0
473 char attrname[H4_MAX_NC_NAME + 1];
474 vector<char> attrbuf, attrbuf2;
475
476 // Here we cannot check if SDfindattr fails or not since even SDfindattr fails it doesn't mean
477 // errors happen. If no such attribute can be found, SDfindattr still returns FAIL.
478 // The correct way is to use SDgetinfo and SDattrinfo to check if attributes "radiance_scales" etc exist.
479 // For the time being, I won't do this, due to the performance reason and code simplity and also the
480 // very small chance of real FAIL for SDfindattr.
481 cf_general_attrindex = SDfindattr(sdsid, "radiance_scales");
482 cf_general_attrindex2 = SDfindattr(sdsid, "radiance_offsets");
483
484 // Obtain the values of radiance_scales and radiance_offsets if they are available.
485 if(cf_general_attrindex!=FAIL && cf_general_attrindex2!=FAIL)
486 {
487 intn ret = -1;
488 ret = SDattrinfo(sdsid, cf_general_attrindex, attrname, &attr_dtype, &temp_attrcount);
489 if (ret==FAIL)
490 {
491 SDendaccess(sdsid);
492 if(true == isgeofile)
493 SDend(sdfileid);
494 ostringstream eherr;
495 eherr << "Attribute 'radiance_scales' in " << fieldname.c_str () << " cannot be obtained.";
496 throw InternalErr (__FILE__, __LINE__, eherr.str ());
497 }
498 attrbuf.clear();
499 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
500 ret = SDreadattr(sdsid, cf_general_attrindex, (VOIDP)&attrbuf[0]);
501 if (ret==FAIL)
502 {
503 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
504 SDendaccess(sdsid);
505 if (true == isgeofile)
506 SDend(sdfileid);
507 ostringstream eherr;
508 eherr << "Attribute 'radiance_scales' in " << fieldname.c_str () << " cannot be obtained.";
509 throw InternalErr (__FILE__, __LINE__, eherr.str ());
510 }
511 ret = SDattrinfo(sdsid, cf_general_attrindex2, attrname, &attr_dtype, &temp_attrcount);
512 if (ret==FAIL)
513 {
514 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
515 SDendaccess(sdsid);
516 if(true == isgeofile)
517 SDend(sdfileid);
518 ostringstream eherr;
519 eherr << "Attribute 'radiance_offsets' in " << fieldname.c_str () << " cannot be obtained.";
520 throw InternalErr (__FILE__, __LINE__, eherr.str ());
521 }
522 attrbuf2.clear();
523 attrbuf2.resize(DFKNTsize(attr_dtype)*temp_attrcount);
524 ret = SDreadattr(sdsid, cf_general_attrindex2, (VOIDP)&attrbuf2[0]);
525 if (ret==FAIL)
526 {
527 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
528 SDendaccess(sdsid);
529 if(true == isgeofile)
530 SDend(sdfileid);
531 ostringstream eherr;
532 eherr << "Attribute 'radiance_offsets' in " << fieldname.c_str () << " cannot be obtained.";
533 throw InternalErr (__FILE__, __LINE__, eherr.str ());
534 }
535
536 // The following macro will obtain radiance_scales and radiance_offsets.
537 // Although the code is compact, it may not be easy to follow. The similar macro can also be found later.
538 switch(attr_dtype)
539 {
540#define GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(TYPE, CAST) \
541 case DFNT_##TYPE: \
542 { \
543 CAST *ptr = (CAST*)&attrbuf[0]; \
544 CAST *ptr2 = (CAST*)&attrbuf2[0]; \
545 radiance_scales = new float[temp_attrcount]; \
546 radiance_offsets = new float[temp_attrcount]; \
547 for(int l=0; l<temp_attrcount; l++) \
548 { \
549 radiance_scales[l] = ptr[l]; \
550 radiance_offsets[l] = ptr2[l]; \
551 } \
552 } \
553 break;
554 GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT32, float);
555 GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT64, double);
556 }
557#undef GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES
558 num_eles_of_an_attr = temp_attrcount; // Store the count of attributes.
559 }
560
561 // Obtain attribute values of reflectance_scales and reflectance_offsets if they are available.
562 cf_general_attrindex = SDfindattr(sdsid, "reflectance_scales");
563 cf_general_attrindex2 = SDfindattr(sdsid, "reflectance_offsets");
564 if(cf_general_attrindex!=FAIL && cf_general_attrindex2!=FAIL)
565 {
566 intn ret = -1;
567 ret = SDattrinfo(sdsid, cf_general_attrindex, attrname, &attr_dtype, &temp_attrcount);
568 if (ret==FAIL)
569 {
570 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
571 SDendaccess(sdsid);
572 if(true == isgeofile)
573 SDend(sdfileid);
574 ostringstream eherr;
575 eherr << "Attribute 'reflectance_scales' in " << fieldname.c_str () << " cannot be obtained.";
576 throw InternalErr (__FILE__, __LINE__, eherr.str ());
577 }
578 attrbuf.clear();
579 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
580 ret = SDreadattr(sdsid, cf_general_attrindex, (VOIDP)&attrbuf[0]);
581 if (ret==FAIL)
582 {
583 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
584 SDendaccess(sdsid);
585 if(true == isgeofile)
586 SDend(sdfileid);
587 ostringstream eherr;
588 eherr << "Attribute 'reflectance_scales' in " << fieldname.c_str () << " cannot be obtained.";
589 throw InternalErr (__FILE__, __LINE__, eherr.str ());
590 }
591
592 ret = SDattrinfo(sdsid, cf_general_attrindex2, attrname, &attr_dtype, &temp_attrcount);
593 if (ret==FAIL)
594 {
595 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
596 SDendaccess(sdsid);
597 if(true == isgeofile)
598 SDend(sdfileid);
599 ostringstream eherr;
600 eherr << "Attribute 'reflectance_offsets' in " << fieldname.c_str () << " cannot be obtained.";
601 throw InternalErr (__FILE__, __LINE__, eherr.str ());
602 }
603 attrbuf2.clear();
604 attrbuf2.resize(DFKNTsize(attr_dtype)*temp_attrcount);
605 ret = SDreadattr(sdsid, cf_general_attrindex2, (VOIDP)&attrbuf2[0]);
606 if (ret==FAIL)
607 {
608 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
609 SDendaccess(sdsid);
610 if(true == isgeofile)
611 SDend(sdfileid);
612 ostringstream eherr;
613 eherr << "Attribute 'reflectance_offsets' in " << fieldname.c_str () << " cannot be obtained.";
614 throw InternalErr (__FILE__, __LINE__, eherr.str ());
615 }
616 switch(attr_dtype)
617 {
618#define GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(TYPE, CAST) \
619 case DFNT_##TYPE: \
620 { \
621 CAST *ptr = (CAST*)&attrbuf[0]; \
622 CAST *ptr2 = (CAST*)&attrbuf2[0]; \
623 reflectance_scales = new float[temp_attrcount]; \
624 reflectance_offsets = new float[temp_attrcount]; \
625 for(int l=0; l<temp_attrcount; l++) \
626 { \
627 reflectance_scales[l] = ptr[l]; \
628 reflectance_offsets[l] = ptr2[l]; \
629 } \
630 } \
631 break;
632 GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT32, float);
633 GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT64, double);
634 }
635#undef GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES
636 num_eles_of_an_attr = temp_attrcount; // Store the count of attributes.
637 }
638
639#endif
640 // Obtain the value of attribute scale_factor.
641 scale_factor_attr_index = SDfindattr(sdsid, "scale_factor");
642 if(scale_factor_attr_index!=FAIL)
643 {
644 intn ret = -1;
645 ret = SDattrinfo(sdsid, scale_factor_attr_index, attrname,
646 &attr_dtype, &temp_attrcount);
647 if (ret==FAIL)
648 {
649 SDendaccess(sdsid);
650 if(true == isgeofile || false == check_pass_fileid_key)
651 SDend(sdfileid);
652 ostringstream eherr;
653 eherr << "Attribute 'scale_factor' in "
654 << fieldname.c_str () << " cannot be obtained.";
655 throw InternalErr (__FILE__, __LINE__, eherr.str ());
656 }
657 attrbuf.clear();
658 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
659 ret = SDreadattr(sdsid, scale_factor_attr_index, (VOIDP)&attrbuf[0]);
660 if (ret==FAIL)
661 {
662 SDendaccess(sdsid);
663 if(true == isgeofile || false == check_pass_fileid_key)
664 SDend(sdfileid);
665
666 ostringstream eherr;
667 eherr << "Attribute 'scale_factor' in "
668 << fieldname.c_str () << " cannot be obtained.";
669 throw InternalErr (__FILE__, __LINE__, eherr.str ());
670 }
671
672 switch(attr_dtype)
673 {
674#define GET_SCALE_FACTOR_ATTR_VALUE(TYPE, CAST) \
675 case DFNT_##TYPE: \
676 { \
677 CAST tmpvalue = *(CAST*)&attrbuf[0]; \
678 scale = (float)tmpvalue; \
679 } \
680 break;
681 GET_SCALE_FACTOR_ATTR_VALUE(INT8, int8);
682 GET_SCALE_FACTOR_ATTR_VALUE(CHAR,int8);
683 GET_SCALE_FACTOR_ATTR_VALUE(UINT8, uint8);
684 GET_SCALE_FACTOR_ATTR_VALUE(UCHAR,uint8);
685 GET_SCALE_FACTOR_ATTR_VALUE(INT16, int16);
686 GET_SCALE_FACTOR_ATTR_VALUE(UINT16, uint16);
687 GET_SCALE_FACTOR_ATTR_VALUE(INT32, int32);
688 GET_SCALE_FACTOR_ATTR_VALUE(UINT32, uint32);
689 GET_SCALE_FACTOR_ATTR_VALUE(FLOAT32, float);
690 GET_SCALE_FACTOR_ATTR_VALUE(FLOAT64, double);
691 default:
692 throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
693
694
695 };
696#undef GET_SCALE_FACTOR_ATTR_VALUE
697 }
698
699 // Obtain the value of attribute add_offset
700 add_offset_attr_index = SDfindattr(sdsid, "add_offset");
701 if(add_offset_attr_index!=FAIL)
702 {
703 intn ret;
704 ret = SDattrinfo(sdsid, add_offset_attr_index, attrname,
705 &attr_dtype, &temp_attrcount);
706 if (ret==FAIL)
707 {
708 SDendaccess(sdsid);
709 if(true == isgeofile || false == check_pass_fileid_key)
710 SDend(sdfileid);
711
712 ostringstream eherr;
713 eherr << "Attribute 'add_offset' in " << fieldname.c_str ()
714 << " cannot be obtained.";
715 throw InternalErr (__FILE__, __LINE__, eherr.str ());
716 }
717 attrbuf.clear();
718 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
719 ret = SDreadattr(sdsid, add_offset_attr_index, (VOIDP)&attrbuf[0]);
720 if (ret==FAIL)
721 {
722 SDendaccess(sdsid);
723 if(true == isgeofile || false == check_pass_fileid_key)
724 SDend(sdfileid);
725
726 ostringstream eherr;
727 eherr << "Attribute 'add_offset' in " << fieldname.c_str ()
728 << " cannot be obtained.";
729 throw InternalErr (__FILE__, __LINE__, eherr.str ());
730 }
731
732 switch(attr_dtype)
733 {
734#define GET_ADD_OFFSET_ATTR_VALUE(TYPE, CAST) \
735 case DFNT_##TYPE: \
736 { \
737 CAST tmpvalue = *(CAST*)&attrbuf[0]; \
738 field_offset = (float)tmpvalue; \
739 } \
740 break;
741 GET_ADD_OFFSET_ATTR_VALUE(INT8, int8);
742 GET_ADD_OFFSET_ATTR_VALUE(CHAR,int8);
743 GET_ADD_OFFSET_ATTR_VALUE(UINT8, uint8);
744 GET_ADD_OFFSET_ATTR_VALUE(UCHAR,uint8);
745 GET_ADD_OFFSET_ATTR_VALUE(INT16, int16);
746 GET_ADD_OFFSET_ATTR_VALUE(UINT16, uint16);
747 GET_ADD_OFFSET_ATTR_VALUE(INT32, int32);
748 GET_ADD_OFFSET_ATTR_VALUE(UINT32, uint32);
749 GET_ADD_OFFSET_ATTR_VALUE(FLOAT32, float);
750 GET_ADD_OFFSET_ATTR_VALUE(FLOAT64, double);
751 default:
752 throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
753
754 };
755#undef GET_ADD_OFFSET_ATTR_VALUE
756 }
757
758 // Obtain the value of the attribute _FillValue
759 cf_fv_attrindex = SDfindattr(sdsid, "_FillValue");
760 if(cf_fv_attrindex!=FAIL)
761 {
762 intn ret;
763 ret = SDattrinfo(sdsid, cf_fv_attrindex, attrname, &attr_dtype, &temp_attrcount);
764 if (ret==FAIL)
765 {
766 SDendaccess(sdsid);
767 if(true == isgeofile || false == check_pass_fileid_key)
768 SDend(sdfileid);
769
770 ostringstream eherr;
771 eherr << "Attribute '_FillValue' in " << fieldname.c_str ()
772 << " cannot be obtained.";
773 throw InternalErr (__FILE__, __LINE__, eherr.str ());
774 }
775 attrbuf.clear();
776 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
777 ret = SDreadattr(sdsid, cf_fv_attrindex, (VOIDP)&attrbuf[0]);
778 if (ret==FAIL)
779 {
780 SDendaccess(sdsid);
781 if(true == isgeofile || false == check_pass_fileid_key)
782 SDend(sdfileid);
783
784 ostringstream eherr;
785 eherr << "Attribute '_FillValue' in " << fieldname.c_str ()
786 << " cannot be obtained.";
787 throw InternalErr (__FILE__, __LINE__, eherr.str ());
788 }
789
790 switch(attr_dtype)
791 {
792#define GET_FILLVALUE_ATTR_VALUE(TYPE, CAST) \
793 case DFNT_##TYPE: \
794 { \
795 CAST tmpvalue = *(CAST*)&attrbuf[0]; \
796 fillvalue = (float)tmpvalue; \
797 } \
798 break;
799 GET_FILLVALUE_ATTR_VALUE(INT8, int8);
800 GET_FILLVALUE_ATTR_VALUE(CHAR, int8);
801 GET_FILLVALUE_ATTR_VALUE(INT16, int16);
802 GET_FILLVALUE_ATTR_VALUE(INT32, int32);
803 GET_FILLVALUE_ATTR_VALUE(UINT8, uint8);
804 GET_FILLVALUE_ATTR_VALUE(UCHAR, uint8);
805 GET_FILLVALUE_ATTR_VALUE(UINT16, uint16);
806 GET_FILLVALUE_ATTR_VALUE(UINT32, uint32);
807 default:
808 throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
809
810 };
811#undef GET_FILLVALUE_ATTR_VALUE
812 }
813
814 // Retrieve valid_range,valid_range is normally represented as (valid_min,valid_max)
815 // for non-CF scale and offset rules, the data is always float. So we only
816 // need to change the data type to float.
817 cf_vr_attrindex = SDfindattr(sdsid, "valid_range");
818 if(cf_vr_attrindex!=FAIL)
819 {
820 intn ret;
821 ret = SDattrinfo(sdsid, cf_vr_attrindex, attrname, &attr_dtype, &temp_attrcount);
822 if (ret==FAIL)
823 {
824 SDendaccess(sdsid);
825 if(true == isgeofile || false == check_pass_fileid_key)
826 SDend(sdfileid);
827
828 ostringstream eherr;
829 eherr << "Attribute '_FillValue' in " << fieldname.c_str ()
830 << " cannot be obtained.";
831 throw InternalErr (__FILE__, __LINE__, eherr.str ());
832 }
833 attrbuf.clear();
834 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
835 ret = SDreadattr(sdsid, cf_vr_attrindex, (VOIDP)&attrbuf[0]);
836 if (ret==FAIL)
837 {
838 SDendaccess(sdsid);
839 if(true == isgeofile || false == check_pass_fileid_key)
840 SDend(sdfileid);
841
842 ostringstream eherr;
843 eherr << "Attribute '_FillValue' in " << fieldname.c_str ()
844 << " cannot be obtained.";
845 throw InternalErr (__FILE__, __LINE__, eherr.str ());
846 }
847
848
849 string attrbuf_str(attrbuf.begin(),attrbuf.end());
850
851 switch(attr_dtype) {
852
853 case DFNT_CHAR:
854 {
855 // We need to treat the attribute data as characters or string.
856 // So find the separator.
857 size_t found = attrbuf_str.find_first_of(",");
858 size_t found_from_end = attrbuf_str.find_last_of(",");
859
860 if (string::npos == found){
861 SDendaccess(sdsid);
862 if(true == isgeofile || false == check_pass_fileid_key)
863 SDend(sdfileid);
864 throw InternalErr(__FILE__,__LINE__,"should find the separator ,");
865 }
866 if (found != found_from_end){
867 SDendaccess(sdsid);
868 if(true == isgeofile || false == check_pass_fileid_key)
869 SDend(sdfileid);
870 throw InternalErr(__FILE__,__LINE__,
871 "Only one separator , should be available.");
872 }
873
874#if 0
875 //istringstream(attrbuf_str.substr(0,found))>> orig_valid_min;
876 //istringstream(attrbuf_str.substr(found+1))>> orig_valid_max;
877#endif
878
879 orig_valid_min = atof((attrbuf_str.substr(0,found)).c_str());
880 orig_valid_max = atof((attrbuf_str.substr(found+1)).c_str());
881
882 }
883 break;
884 case DFNT_INT8:
885 {
886 // We find a special case that even valid_range is logically
887 // interpreted as two elements,
888 // but the count of attribute elements is more than 2. The count
889 // actually is the number of
890 // characters stored as the attribute value. So we need to find
891 // the separator "," and then
892 // change the string before and after the separator into float numbers.
893 //
894 if (temp_attrcount >2) {
895
896 size_t found = attrbuf_str.find_first_of(",");
897 size_t found_from_end = attrbuf_str.find_last_of(",");
898
899 if (string::npos == found){
900 SDendaccess(sdsid);
901 if(true == isgeofile || false == check_pass_fileid_key)
902 SDend(sdfileid);
903 throw InternalErr(__FILE__,__LINE__,"should find the separator ,");
904 }
905 if (found != found_from_end){
906 SDendaccess(sdsid);
907 if(true == isgeofile || false == check_pass_fileid_key)
908 SDend(sdfileid);
909 throw InternalErr(__FILE__,__LINE__,
910 "Only one separator , should be available.");
911 }
912
913 //istringstream(attrbuf_str.substr(0,found))>> orig_valid_min;
914 //istringstream(attrbuf_str.substr(found+1))>> orig_valid_max;
915
916 orig_valid_min = atof((attrbuf_str.substr(0,found)).c_str());
917 orig_valid_max = atof((attrbuf_str.substr(found+1)).c_str());
918
919 }
920 else if (2 == temp_attrcount) {
921 orig_valid_min = (float)attrbuf[0];
922 orig_valid_max = (float)attrbuf[1];
923 }
924 else {
925 SDendaccess(sdsid);
926 if(true == isgeofile || false == check_pass_fileid_key)
927 SDend(sdfileid);
928 throw InternalErr(__FILE__,__LINE__,
929 "The number of attribute count should be greater than 1.");
930 }
931
932 }
933 break;
934
935 case DFNT_UINT8:
936 case DFNT_UCHAR:
937 {
938 if (temp_attrcount != 2) {
939 SDendaccess(sdsid);
940 if(true == isgeofile || false == check_pass_fileid_key)
941 SDend(sdfileid);
942
943 throw InternalErr(__FILE__,__LINE__,
944 "The number of attribute count should be 2 for the DFNT_UINT8 type.");
945 }
946
947 unsigned char* temp_valid_range = (unsigned char *)&attrbuf[0];
948 orig_valid_min = (float)(temp_valid_range[0]);
949 orig_valid_max = (float)(temp_valid_range[1]);
950 }
951 break;
952
953 case DFNT_INT16:
954 {
955 if (temp_attrcount != 2) {
956 SDendaccess(sdsid);
957 if(true == isgeofile || false == check_pass_fileid_key)
958 SDend(sdfileid);
959
960 throw InternalErr(__FILE__,__LINE__,
961 "The number of attribute count should be 2 for the DFNT_INT16 type.");
962 }
963
964 short* temp_valid_range = (short *)&attrbuf[0];
965 orig_valid_min = (float)(temp_valid_range[0]);
966 orig_valid_max = (float)(temp_valid_range[1]);
967 }
968 break;
969
970 case DFNT_UINT16:
971 {
972 if (temp_attrcount != 2) {
973 SDendaccess(sdsid);
974 if(true == isgeofile || false == check_pass_fileid_key)
975 SDend(sdfileid);
976
977 throw InternalErr(__FILE__,__LINE__,
978 "The number of attribute count should be 2 for the DFNT_UINT16 type.");
979 }
980
981 unsigned short* temp_valid_range = (unsigned short *)&attrbuf[0];
982 orig_valid_min = (float)(temp_valid_range[0]);
983 orig_valid_max = (float)(temp_valid_range[1]);
984 }
985 break;
986
987 case DFNT_INT32:
988 {
989 if (temp_attrcount != 2) {
990 SDendaccess(sdsid);
991 if(true == isgeofile || false == check_pass_fileid_key)
992 SDend(sdfileid);
993
994 throw InternalErr(__FILE__,__LINE__,
995 "The number of attribute count should be 2 for the DFNT_INT32 type.");
996 }
997
998 int* temp_valid_range = (int *)&attrbuf[0];
999 orig_valid_min = (float)(temp_valid_range[0]);
1000 orig_valid_max = (float)(temp_valid_range[1]);
1001 }
1002 break;
1003
1004 case DFNT_UINT32:
1005 {
1006 if (temp_attrcount != 2) {
1007 SDendaccess(sdsid);
1008 if(true == isgeofile || false == check_pass_fileid_key)
1009 SDend(sdfileid);
1010
1011 throw InternalErr(__FILE__,__LINE__,
1012 "The number of attribute count should be 2 for the DFNT_UINT32 type.");
1013 }
1014
1015 unsigned int* temp_valid_range = (unsigned int *)&attrbuf[0];
1016 orig_valid_min = (float)(temp_valid_range[0]);
1017 orig_valid_max = (float)(temp_valid_range[1]);
1018 }
1019 break;
1020
1021 case DFNT_FLOAT32:
1022 {
1023 if (temp_attrcount != 2) {
1024 SDendaccess(sdsid);
1025 if(true == isgeofile || false == check_pass_fileid_key)
1026 SDend(sdfileid);
1027
1028 throw InternalErr(__FILE__,__LINE__,
1029 "The number of attribute count should be 2 for the DFNT_FLOAT32 type.");
1030 }
1031
1032 float* temp_valid_range = (float *)&attrbuf[0];
1033 orig_valid_min = temp_valid_range[0];
1034 orig_valid_max = temp_valid_range[1];
1035 }
1036 break;
1037
1038 case DFNT_FLOAT64:
1039 {
1040 if (temp_attrcount != 2){
1041 SDendaccess(sdsid);
1042 if(true == isgeofile || false == check_pass_fileid_key)
1043 SDend(sdfileid);
1044
1045 throw InternalErr(__FILE__,__LINE__,
1046 "The number of attribute count should be 2 for the DFNT_FLOAT64 type.");
1047 }
1048 double* temp_valid_range = (double *)&attrbuf[0];
1049
1050 // Notice: this approach will lose precision and possibly overflow.
1051 // Fortunately it is not a problem for MODIS data.
1052 // This part of code may not be called.
1053 // If it is called, mostly the value is within the floating-point range.
1054 // KY 2013-01-29
1055 orig_valid_min = temp_valid_range[0];
1056 orig_valid_max = temp_valid_range[1];
1057 }
1058 break;
1059 default: {
1060 SDendaccess(sdsid);
1061 if(true == isgeofile || false == check_pass_fileid_key)
1062 SDend(sdfileid);
1063 throw InternalErr(__FILE__,__LINE__,"Unsupported data type.");
1064 }
1065 }
1066 }
1067
1068 // Check if the data has the "Key" attribute.
1069 // We found that some NSIDC MODIS data(MOD29) used "Key" to identify some special values.
1070 // To get the values that are within the range identified by the "Key",
1071 // scale offset rules also need to be applied to those values
1072 // outside the original valid range. KY 2013-02-25
1073 int32 cf_mod_key_attrindex = SUCCEED;
1074 cf_mod_key_attrindex = SDfindattr(sdsid, "Key");
1075 if(cf_mod_key_attrindex !=FAIL) {
1076 has_Key_attr = true;
1077 }
1078
1079 attrbuf.clear();
1080 vector<char> attrbuf2;
1081
1082 // Here we cannot check if SDfindattr fails since even SDfindattr fails it doesn't mean
1083 // errors happen. If no such attribute can be found, SDfindattr still returns FAIL.
1084 // The correct way is to use SDgetinfo and SDattrinfo to check if attributes
1085 // "radiance_scales" etc exist.
1086 // For the time being, I won't do this, due to the performance reason and code simplity
1087 // and also the very small chance of real FAIL for SDfindattr.
1088 cf_modl1b_rr_attrindex = SDfindattr(sdsid, "radiance_scales");
1089 cf_modl1b_rr_attrindex2 = SDfindattr(sdsid, "radiance_offsets");
1090
1091 // Obtain the values of radiance_scales and radiance_offsets if they are available.
1092 if(cf_modl1b_rr_attrindex!=FAIL && cf_modl1b_rr_attrindex2!=FAIL)
1093 {
1094 intn ret = -1;
1095 ret = SDattrinfo(sdsid, cf_modl1b_rr_attrindex, attrname,
1096 &attr_dtype, &temp_attrcount);
1097 if (ret==FAIL)
1098 {
1099 SDendaccess(sdsid);
1100 if(true == isgeofile || false == check_pass_fileid_key)
1101 SDend(sdfileid);
1102 ostringstream eherr;
1103 eherr << "Attribute 'radiance_scales' in " << fieldname.c_str ()
1104 << " cannot be obtained.";
1105 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1106 }
1107 attrbuf.clear();
1108 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
1109 ret = SDreadattr(sdsid, cf_modl1b_rr_attrindex, (VOIDP)&attrbuf[0]);
1110 if (ret==FAIL)
1111 {
1112 SDendaccess(sdsid);
1113 if (true == isgeofile || false == check_pass_fileid_key)
1114 SDend(sdfileid);
1115 ostringstream eherr;
1116 eherr << "Attribute 'radiance_scales' in " << fieldname.c_str ()
1117 << " cannot be obtained.";
1118 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1119 }
1120 ret = SDattrinfo(sdsid, cf_modl1b_rr_attrindex2, attrname,
1121 &attr_dtype, &temp_attrcount);
1122 if (ret==FAIL)
1123 {
1124 SDendaccess(sdsid);
1125 if(true == isgeofile || false == check_pass_fileid_key)
1126 SDend(sdfileid);
1127 ostringstream eherr;
1128 eherr << "Attribute 'radiance_offsets' in "
1129 << fieldname.c_str () << " cannot be obtained.";
1130 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1131 }
1132 attrbuf2.clear();
1133 attrbuf2.resize(DFKNTsize(attr_dtype)*temp_attrcount);
1134 ret = SDreadattr(sdsid, cf_modl1b_rr_attrindex2, (VOIDP)&attrbuf2[0]);
1135 if (ret==FAIL)
1136 {
1137 SDendaccess(sdsid);
1138 if(true == isgeofile || false == check_pass_fileid_key)
1139 SDend(sdfileid);
1140 ostringstream eherr;
1141 eherr << "Attribute 'radiance_offsets' in "
1142 << fieldname.c_str () << " cannot be obtained.";
1143 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1144 }
1145
1146 // The following macro will obtain radiance_scales and radiance_offsets.
1147 // Although the code is compact, it may not be easy to follow.
1148 // The similar macro can also be found later.
1149 switch(attr_dtype)
1150 {
1151#define GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(TYPE, CAST) \
1152 case DFNT_##TYPE: \
1153 { \
1154 CAST *ptr = (CAST*)&attrbuf[0]; \
1155 CAST *ptr2 = (CAST*)&attrbuf2[0]; \
1156 radiance_scales = new float[temp_attrcount]; \
1157 radiance_offsets = new float[temp_attrcount]; \
1158 for(int l=0; l<temp_attrcount; l++) \
1159 { \
1160 radiance_scales[l] = ptr[l]; \
1161 radiance_offsets[l] = ptr2[l]; \
1162 } \
1163 } \
1164 break;
1165 GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT32, float);
1166 GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT64, double);
1167 default:
1168 throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
1169
1170 }
1171#undef GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES
1172 // Store the count of attributes.
1173 num_eles_of_an_attr = temp_attrcount;
1174 }
1175
1176 // Obtain attribute values of reflectance_scales
1177 // and reflectance_offsets if they are available.
1178 cf_modl1b_rr_attrindex = SDfindattr(sdsid, "reflectance_scales");
1179 cf_modl1b_rr_attrindex2 = SDfindattr(sdsid, "reflectance_offsets");
1180 if(cf_modl1b_rr_attrindex!=FAIL && cf_modl1b_rr_attrindex2!=FAIL)
1181 {
1182 intn ret = -1;
1183 ret = SDattrinfo(sdsid, cf_modl1b_rr_attrindex, attrname,
1184 &attr_dtype, &temp_attrcount);
1185 if (ret==FAIL)
1186 {
1187 release_mod1b_res(reflectance_scales,reflectance_offsets,
1188 radiance_scales,radiance_offsets);
1189 SDendaccess(sdsid);
1190 if(true == isgeofile || false == check_pass_fileid_key)
1191 SDend(sdfileid);
1192 ostringstream eherr;
1193 eherr << "Attribute 'reflectance_scales' in "
1194 << fieldname.c_str () << " cannot be obtained.";
1195 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1196 }
1197 attrbuf.clear();
1198 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
1199 ret = SDreadattr(sdsid, cf_modl1b_rr_attrindex, (VOIDP)&attrbuf[0]);
1200 if (ret==FAIL)
1201 {
1202 release_mod1b_res(reflectance_scales,reflectance_offsets,
1203 radiance_scales,radiance_offsets);
1204 SDendaccess(sdsid);
1205 if(true == isgeofile || false == check_pass_fileid_key)
1206 SDend(sdfileid);
1207 ostringstream eherr;
1208 eherr << "Attribute 'reflectance_scales' in "
1209 << fieldname.c_str () << " cannot be obtained.";
1210 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1211 }
1212
1213 ret = SDattrinfo(sdsid, cf_modl1b_rr_attrindex2, attrname,
1214 &attr_dtype, &temp_attrcount);
1215 if (ret==FAIL)
1216 {
1217 release_mod1b_res(reflectance_scales,reflectance_offsets,
1218 radiance_scales,radiance_offsets);
1219 SDendaccess(sdsid);
1220 if(true == isgeofile || false == check_pass_fileid_key)
1221 SDend(sdfileid);
1222 ostringstream eherr;
1223 eherr << "Attribute 'reflectance_offsets' in "
1224 << fieldname.c_str () << " cannot be obtained.";
1225 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1226 }
1227 attrbuf2.clear();
1228 attrbuf2.resize(DFKNTsize(attr_dtype)*temp_attrcount);
1229 ret = SDreadattr(sdsid, cf_modl1b_rr_attrindex2, (VOIDP)&attrbuf2[0]);
1230 if (ret==FAIL)
1231 {
1232 release_mod1b_res(reflectance_scales,reflectance_offsets,
1233 radiance_scales,radiance_offsets);
1234 SDendaccess(sdsid);
1235 if(true == isgeofile || false == check_pass_fileid_key)
1236 SDend(sdfileid);
1237 ostringstream eherr;
1238 eherr << "Attribute 'reflectance_offsets' in "
1239 << fieldname.c_str () << " cannot be obtained.";
1240 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1241 }
1242 switch(attr_dtype)
1243 {
1244#define GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(TYPE, CAST) \
1245 case DFNT_##TYPE: \
1246 { \
1247 CAST *ptr = (CAST*)&attrbuf[0]; \
1248 CAST *ptr2 = (CAST*)&attrbuf2[0]; \
1249 reflectance_scales = new float[temp_attrcount]; \
1250 reflectance_offsets = new float[temp_attrcount]; \
1251 for(int l=0; l<temp_attrcount; l++) \
1252 { \
1253 reflectance_scales[l] = ptr[l]; \
1254 reflectance_offsets[l] = ptr2[l]; \
1255 } \
1256 } \
1257 break;
1258 GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT32, float);
1259 GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT64, double);
1260 default:
1261 throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
1262
1263 }
1264#undef GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES
1265 num_eles_of_an_attr = temp_attrcount; // Store the count of attributes.
1266 }
1267
1268 SDendaccess(sdsid);
1270#if 0
1271r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
1272 &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1273if (r != 0) {
1274 ostringstream eherr;
1275
1276 eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
1277 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1278}
1279
1280cerr<<"tmp_rank 3 is "<<tmp_rank <<endl;
1281#endif
1282 //if (true == isgeofile || false == check_pass_fileid_key)
1283 // SDend(sdfileid);
1285#if 0
1286r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
1287 &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1288if (r != 0) {
1289 ostringstream eherr;
1290
1291 eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
1292 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1293}
1294
1295cerr<<"tmp_rank 4 is "<<tmp_rank <<endl;
1296#endif
1297
1298 BESDEBUG("h4","scale is "<<scale <<endl);
1299 BESDEBUG("h4","offset is "<<field_offset <<endl);
1300 BESDEBUG("h4","fillvalue is "<<fillvalue <<endl);
1301 }
1302 }
1303
1304 // According to our observations, it seems that MODIS products ALWAYS
1305 // use the "scale" factor to make bigger values smaller.
1306 // So for MODIS_MUL_SCALE products, if the scale of some variable is greater than 1,
1307 // it means that for this variable, the MODIS type for this variable may be MODIS_DIV_SCALE.
1308 // For the similar logic, we may need to change MODIS_DIV_SCALE to MODIS_MUL_SCALE
1309 // and MODIS_EQ_SCALE to MODIS_DIV_SCALE.
1310 // We indeed find such a case. HDF-EOS2 Grid MODIS_Grid_1km_2D of MOD(or MYD)09GA is
1311 // a MODIS_EQ_SCALE.
1312 // However,
1313 // the scale_factor of the variable Range_1 in the MOD09GA product is 25.
1314 // According to our observation,
1315 // this variable should be MODIS_DIV_SCALE.We verify this is true according to
1316 // MODIS 09 product document
1317 // http://modis-sr.ltdri.org/products/MOD09_UserGuide_v1_3.pdf.
1318 // Since this conclusion is based on our observation, we would like to add a BESlog to detect
1319 // if we find
1320 // the similar cases so that we can verify with the corresponding product documents to see if
1321 // this is true.
1322 // More updated information,
1323 // We just verified with the MOD09 data producer, the scale_factor for Range_1 is 25
1324 // but the equation is still multiplication instead of division.
1325 // So we have to make this as a special case and don't change the scale and offset settings
1326 // for Range_1 of MOD09 products.
1327 // KY 2014-01-13
1328
1329 if (MODIS_EQ_SCALE == sotype || MODIS_MUL_SCALE == sotype) {
1330 if (scale > 1) {
1331 bool need_change_scale = true;
1332 if(gridname!="") {
1333
1334 string temp_filename;
1335 if (filename.find("#") != string::npos)
1336 temp_filename =filename.substr(filename.find_last_of("#") + 1);
1337 else
1338 temp_filename = filename.substr(filename.find_last_of("/") +1);
1339
1340 if ((temp_filename.size() >5) && ((temp_filename.compare(0,5,"MOD09") == 0)
1341 ||(temp_filename.compare(0,5,"MYD09") == 0))) {
1342 if ((fieldname.size() >5) && fieldname.compare(0,5,"Range") == 0)
1343 need_change_scale = false;
1344 }
1345 // MOD16A2
1346 else if((temp_filename.size() >7)&&
1347 ((temp_filename.compare(0,7,"MOD16A2") == 0)|| (temp_filename.compare(0,7,"MYD16A2")==0)||
1348 (temp_filename.compare(0,7,"MOD16A3") == 0)|| (temp_filename.compare(0,7,"MYD16A3")==0)))
1349 need_change_scale = false;
1350
1351
1352 }
1353 if(true == need_change_scale) {
1354 sotype = MODIS_DIV_SCALE;
1355 (*BESLog::TheLog())
1356 << "The field " << fieldname << " scale factor is "
1357 << scale << " ."<<endl
1358 << " But the original scale factor type is MODIS_MUL_SCALE or MODIS_EQ_SCALE. "
1359 << endl
1360 << " Now change it to MODIS_DIV_SCALE. "<<endl;
1361 }
1362 }
1363 }
1364
1365 if (MODIS_DIV_SCALE == sotype) {
1366 if (scale < 1) {
1367 sotype = MODIS_MUL_SCALE;
1368 (*BESLog::TheLog())<< "The field " << fieldname << " scale factor is "
1369 << scale << " ."<<endl
1370 << " But the original scale factor type is MODIS_DIV_SCALE. "
1371 << endl
1372 << " Now change it to MODIS_MUL_SCALE. "<<endl;
1373 }
1374 }
1375#if 0
1377r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
1378 &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1379if (r != 0) {
1380 ostringstream eherr;
1381
1382 eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
1383 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1384}
1385
1386cerr<<"tmp_rank 2 is "<<tmp_rank <<endl;
1387#endif
1388
1389#if 0
1390// We need to loop through all datatpes to allocate the memory buffer for the data.
1391// It is hard to add comments to the macro. We may need to change them to general routines in the future.
1392// Some MODIS products use both valid_range(valid_min, valid_max) and fillvalues for data fields. When do recalculating,
1393// I check fillvalue first, then check valid_min and valid_max if they are available.
1394// The middle check is_special_value addresses the MODIS L1B special value.
1395// ////////////////////////////////////////////////////////////////////////////////////
1396/* if((float)tmptr[l] != fillvalue ) \
1397// { \
1398// if(false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[l]))\
1399// { \
1400// if (orig_valid_min<tmpval[l] && orig_valid_max>tmpval[l] \
1401// if(sotype==MODIS_MUL_SCALE) \
1402// tmpval[l] = (tmptr[l]-field_offset)*scale; \
1403// else if(sotype==MODIS_EQ_SCALE) \
1404// tmpval[l] = tmptr[l]*scale + field_offset; \
1405// else if(sotype==MODIS_DIV_SCALE) \
1406// tmpval[l] = (tmptr[l]-field_offset)/scale; \
1407// } \
1408*/
1409#define RECALCULATE(CAST, DODS_CAST, VAL) \
1410{ \
1411 bool change_data_value = false; \
1412 if(sotype!=DEFAULT_CF_EQU) \
1413 { \
1414 vector<float>tmpval; \
1415 tmpval.resize(nelms); \
1416 CAST tmptr = (CAST)VAL; \
1417 for(int l=0; l<nelms; l++) \
1418 tmpval[l] = (float)tmptr[l]; \
1419 bool special_case = false; \
1420 if(scale_factor_attr_index==FAIL) \
1421 if(num_eles_of_an_attr==1) \
1422 if(radiance_scales!=NULL && radiance_offsets!=NULL) \
1423 { \
1424 scale = radiance_scales[0]; \
1425 field_offset = radiance_offsets[0];\
1426 special_case = true; \
1427 } \
1428 if((scale_factor_attr_index!=FAIL && !(scale==1 && field_offset==0)) || special_case) \
1429 { \
1430 for(int l=0; l<nelms; l++) \
1431 { \
1432 if(cf_vr_attrindex!=FAIL) \
1433 { \
1434 if((float)tmptr[l] != fillvalue ) \
1435 { \
1436 if(false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[l]))\
1437 { \
1438 if ((orig_valid_min<=tmpval[l] && orig_valid_max>=tmpval[l]) || (true==has_Key_attr))\
1439 { \
1440 if(sotype==MODIS_MUL_SCALE) \
1441 tmpval[l] = (tmptr[l]-field_offset)*scale; \
1442 else if(sotype==MODIS_EQ_SCALE) \
1443 tmpval[l] = tmptr[l]*scale + field_offset; \
1444 else if(sotype==MODIS_DIV_SCALE) \
1445 tmpval[l] = (tmptr[l]-field_offset)/scale;\
1446 } \
1447 } \
1448 } \
1449 } \
1450 } \
1451 change_data_value = true; \
1452 set_value((dods_float32 *)&tmpval[0], nelms); \
1453 } else if(num_eles_of_an_attr>1 && (radiance_scales!=NULL && radiance_offsets!=NULL) || (reflectance_scales!=NULL && reflectance_offsets!=NULL)) \
1454 { \
1455 size_t dimindex=0; \
1456 if( num_eles_of_an_attr!=tmp_dims[dimindex]) \
1457 { \
1458 ostringstream eherr; \
1459 eherr << "The number of Z-Dimension scale attribute is not equal to the size of the first dimension in " << fieldname.c_str() << ". These two values must be equal."; \
1460 throw InternalErr (__FILE__, __LINE__, eherr.str ()); \
1461 } \
1462 size_t start_index, end_index; \
1463 size_t nr_elems = nelms/count32[dimindex]; \
1464 start_index = offset32[dimindex]; \
1465 end_index = start_index+step32[dimindex]*(count32[dimindex]-1); \
1466 size_t index = 0;\
1467 for(size_t k=start_index; k<=end_index; k+=step32[dimindex]) \
1468 { \
1469 float tmpscale = (fieldname.find("Emissive")!=string::npos)? radiance_scales[k]: reflectance_scales[k]; \
1470 float tmpoffset = (fieldname.find("Emissive")!=string::npos)? radiance_offsets[k]: reflectance_offsets[k]; \
1471 for(size_t l=0; l<nr_elems; l++) \
1472 { \
1473 if(cf_vr_attrindex!=FAIL) \
1474 { \
1475 if(((float)tmptr[index])!=fillvalue) \
1476 { \
1477 if(false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[index]))\
1478 { \
1479 if(sotype==MODIS_MUL_SCALE) \
1480 tmpval[index] = (tmptr[index]-tmpoffset)*tmpscale; \
1481 else if(sotype==MODIS_EQ_SCALE) \
1482 tmpval[index] = tmptr[index]*tmpscale+tmpoffset; \
1483 else if(sotype==MODIS_DIV_SCALE) \
1484 tmpval[index] = (tmptr[index]-tmpoffset)/tmpscale; \
1485 } \
1486 } \
1487 } \
1488 index++; \
1489 } \
1490 } \
1491 change_data_value = true; \
1492 set_value((dods_float32 *)&tmpval[0], nelms); \
1493 } \
1494 } \
1495 if(!change_data_value) \
1496 { \
1497 set_value ((DODS_CAST)VAL, nelms); \
1498 } \
1499}
1500#endif
1501
1502// We need to loop through all datatpes to allocate the memory buffer for the data.
1503// It is hard to add comments to the macro. We may need to change them to general
1504// routines in the future.
1505// Some MODIS products use both valid_range(valid_min, valid_max) and fillvalues for data fields.
1506// When do recalculating,
1507// I check fillvalue first, then check valid_min and valid_max if they are available.
1508// The middle check is_special_value addresses the MODIS L1B special value.
1509// Updated: just find that the RECALCULATE will be done only when the valid_range
1510// attribute is present(if cf_vr_attrindex!=FAIL).
1511// This restriction is in theory not necessary, but for more MODIS data products,
1512// this restriction may be valid since valid_range pairs with scale/offset to identify
1513// the valid data values. KY 2014-02-19
1514//
1515/* if((float)tmptr[l] != fillvalue ) \
1516// { \
1517// f(false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[l]))\
1518// { \
1519// if (orig_valid_min<tmpval[l] && orig_valid_max>tmpval[l] \
1520// if(sotype==MODIS_MUL_SCALE) \
1521// tmpval[l] = (tmptr[l]-field_offset)*scale; \
1522// else if(sotype==MODIS_EQ_SCALE) \
1523// tmpval[l] = tmptr[l]*scale + field_offset; \
1524// else if(sotype==MODIS_DIV_SCALE) \
1525// tmpval[l] = (tmptr[l]-field_offset)/scale; \
1526// } \
1527*/
1528#define RECALCULATE(CAST, DODS_CAST, VAL) \
1529{ \
1530 bool change_data_value = false; \
1531 if(sotype!=DEFAULT_CF_EQU) \
1532 { \
1533 vector<float>tmpval; \
1534 tmpval.resize(nelms); \
1535 CAST tmptr = (CAST)VAL; \
1536 for(int l=0; l<nelms; l++) \
1537 tmpval[l] = (float)tmptr[l]; \
1538 bool special_case = false; \
1539 if(scale_factor_attr_index==FAIL) \
1540 if(num_eles_of_an_attr==1) \
1541 if((radiance_scales!=NULL) && (radiance_offsets!=NULL)) \
1542 { \
1543 scale = radiance_scales[0]; \
1544 field_offset = radiance_offsets[0];\
1545 special_case = true; \
1546 } \
1547 if(((scale_factor_attr_index!=FAIL) && !((scale==1) && (field_offset==0))) || special_case) \
1548 { \
1549 float temp_scale = scale; \
1550 float temp_offset = field_offset; \
1551 if(sotype==MODIS_MUL_SCALE) \
1552 temp_offset = -1. *field_offset*temp_scale;\
1553 else if (sotype==MODIS_DIV_SCALE) \
1554 {\
1555 temp_scale = 1/scale; \
1556 temp_offset = -1. *field_offset*temp_scale;\
1557 }\
1558 for(int l=0; l<nelms; l++) \
1559 { \
1560 if(cf_vr_attrindex!=FAIL) \
1561 { \
1562 if((float)tmptr[l] != fillvalue ) \
1563 { \
1564 if(false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[l]))\
1565 { \
1566 if (((orig_valid_min<=tmpval[l]) && (orig_valid_max>=tmpval[l])) || (true==has_Key_attr))\
1567 { \
1568 tmpval[l] = tmptr[l]*temp_scale + temp_offset; \
1569 } \
1570 } \
1571 } \
1572 } \
1573 } \
1574 change_data_value = true; \
1575 set_value((dods_float32 *)&tmpval[0], nelms); \
1576 } else if((num_eles_of_an_attr>1) && (((radiance_scales!=NULL) && (radiance_offsets!=NULL)) || ((reflectance_scales!=NULL) && (reflectance_offsets!=NULL)))) \
1577 { \
1578 size_t dimindex=0; \
1579 if( num_eles_of_an_attr!=tmp_dims[dimindex]) \
1580 { \
1581 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets); \
1582 ostringstream eherr; \
1583 eherr << "The number of Z-Dimension scale attribute is not equal to the size of the first dimension in " << fieldname.c_str() << ". These two values must be equal."; \
1584 throw InternalErr (__FILE__, __LINE__, eherr.str ()); \
1585 } \
1586 size_t start_index, end_index; \
1587 size_t nr_elems = nelms/count32[dimindex]; \
1588 start_index = offset32[dimindex]; \
1589 end_index = start_index+step32[dimindex]*(count32[dimindex]-1); \
1590 size_t index = 0;\
1591 for(size_t k=start_index; k<=end_index; k+=step32[dimindex]) \
1592 { \
1593 float tmpscale = (fieldname.find("Emissive")!=string::npos)? radiance_scales[k]: reflectance_scales[k]; \
1594 float tmpoffset = (fieldname.find("Emissive")!=string::npos)? radiance_offsets[k]: reflectance_offsets[k]; \
1595 for(size_t l=0; l<nr_elems; l++) \
1596 { \
1597 if(cf_vr_attrindex!=FAIL) \
1598 { \
1599 if(((float)tmptr[index])!=fillvalue) \
1600 { \
1601 if(false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[index]))\
1602 { \
1603 if(sotype==MODIS_MUL_SCALE) \
1604 tmpval[index] = (tmptr[index]-tmpoffset)*tmpscale; \
1605 else if(sotype==MODIS_EQ_SCALE) \
1606 tmpval[index] = tmptr[index]*tmpscale+tmpoffset; \
1607 else if(sotype==MODIS_DIV_SCALE) \
1608 tmpval[index] = (tmptr[index]-tmpoffset)/tmpscale; \
1609 } \
1610 } \
1611 } \
1612 index++; \
1613 } \
1614 } \
1615 change_data_value = true; \
1616 set_value((dods_float32 *)&tmpval[0], nelms); \
1617 } \
1618 } \
1619 if(!change_data_value) \
1620 { \
1621 set_value ((DODS_CAST)VAL, nelms); \
1622 } \
1623}
1624#if 0
1626r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
1627 &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1628if (r != 0) {
1629 ostringstream eherr;
1630
1631 eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
1632 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1633}
1634
1635cerr<<"tmp_rank again is "<<tmp_rank <<endl;
1636#endif
1637 switch (field_dtype) {
1638 case DFNT_INT8:
1639 {
1640
1641 vector<int8>val;
1642 val.resize(nelms);
1643 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1644 &offset32[0], &step32[0], &count32[0], &val[0]);
1645 if (r != 0) {
1646 release_mod1b_res(reflectance_scales,reflectance_offsets,
1647 radiance_scales,radiance_offsets);
1648 ostringstream eherr;
1649 eherr << "field " << fieldname.c_str () << "cannot be read.";
1650 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1651 }
1652
1653#ifndef SIGNED_BYTE_TO_INT32
1654 RECALCULATE(int8*, dods_byte*, &val[0]);
1655#else
1656
1657 vector<int32>newval;
1658 newval.resize(nelms);
1659
1660 for (int counter = 0; counter < nelms; counter++)
1661 newval[counter] = (int32) (val[counter]);
1662
1663 RECALCULATE(int32*, dods_int32*, &newval[0]);
1664#endif
1665 }
1666 break;
1667 case DFNT_UINT8:
1668 case DFNT_UCHAR8:
1669 {
1670
1671 vector<uint8>val;
1672 val.resize(nelms);
1673 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1674 &offset32[0], &step32[0], &count32[0], &val[0]);
1675 if (r != 0) {
1676 release_mod1b_res(reflectance_scales,reflectance_offsets,
1677 radiance_scales,radiance_offsets);
1678 ostringstream eherr;
1679
1680 eherr << "field " << fieldname.c_str () << "cannot be read.";
1681 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1682 }
1683
1684 RECALCULATE(uint8*, dods_byte*, &val[0]);
1685 }
1686 break;
1687
1688 case DFNT_INT16:
1689 {
1690 vector<int16>val;
1691 val.resize(nelms);
1692 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1693 &offset32[0], &step32[0], &count32[0], &val[0]);
1694
1695 if (r != 0) {
1696
1697 release_mod1b_res(reflectance_scales,reflectance_offsets,
1698 radiance_scales,radiance_offsets);
1699 ostringstream eherr;
1700
1701 eherr << "field " << fieldname.c_str () << "cannot be read.";
1702 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1703 }
1704 RECALCULATE(int16*, dods_int16*, &val[0]);
1705 }
1706 break;
1707 case DFNT_UINT16:
1708 {
1709 vector<uint16>val;
1710 val.resize(nelms);
1711#if 0
1712cerr<<"gridid is "<<gridid <<endl;
1713int32 tmp_rank = 0;
1714char tmp_dimlist[1024];
1715int32 tmp_dims[rank];
1716int32 field_dtype = 0;
1717intn r = 0;
1718
1719r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
1720 &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1721if (r != 0) {
1722 ostringstream eherr;
1723
1724 eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
1725 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1726}
1727#endif
1728
1729 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1730 &offset32[0], &step32[0], &count32[0], &val[0]);
1731 if (r != 0) {
1732 release_mod1b_res(reflectance_scales,reflectance_offsets,
1733 radiance_scales,radiance_offsets);
1734 ostringstream eherr;
1735
1736 eherr << "field " << fieldname.c_str () << "cannot be read.";
1737 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1738 }
1739
1740 RECALCULATE(uint16*, dods_uint16*, &val[0]);
1741 }
1742 break;
1743 case DFNT_INT32:
1744 {
1745 vector<int32>val;
1746 val.resize(nelms);
1747 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1748 &offset32[0], &step32[0], &count32[0], &val[0]);
1749 if (r != 0) {
1750
1751 release_mod1b_res(reflectance_scales,reflectance_offsets,
1752 radiance_scales,radiance_offsets);
1753 ostringstream eherr;
1754
1755 eherr << "field " << fieldname.c_str () << "cannot be read.";
1756 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1757 }
1758
1759 RECALCULATE(int32*, dods_int32*, &val[0]);
1760 }
1761 break;
1762 case DFNT_UINT32:
1763 {
1764 vector<uint32>val;
1765 val.resize(nelms);
1766 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1767 &offset32[0], &step32[0], &count32[0], &val[0]);
1768 if (r != 0) {
1769
1770 release_mod1b_res(reflectance_scales,reflectance_offsets,
1771 radiance_scales,radiance_offsets);
1772 ostringstream eherr;
1773
1774 eherr << "field " << fieldname.c_str () << "cannot be read.";
1775 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1776 }
1777
1778 RECALCULATE(uint32*, dods_uint32*, &val[0]);
1779 }
1780 break;
1781 case DFNT_FLOAT32:
1782 {
1783 vector<float32>val;
1784 val.resize(nelms);
1785 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1786 &offset32[0], &step32[0], &count32[0], &val[0]);
1787 if (r != 0) {
1788
1789 release_mod1b_res(reflectance_scales,reflectance_offsets,
1790 radiance_scales,radiance_offsets);
1791 ostringstream eherr;
1792
1793 eherr << "field " << fieldname.c_str () << "cannot be read.";
1794 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1795 }
1796
1797 // Recalculate seems not necessary.
1798 RECALCULATE(float32*, dods_float32*, &val[0]);
1799 //set_value((dods_float32*)&val[0],nelms);
1800 }
1801 break;
1802 case DFNT_FLOAT64:
1803 {
1804 vector<float64>val;
1805 val.resize(nelms);
1806 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1807 &offset32[0], &step32[0], &count32[0], &val[0]);
1808 if (r != 0) {
1809
1810 release_mod1b_res(reflectance_scales,reflectance_offsets,
1811 radiance_scales,radiance_offsets);
1812 ostringstream eherr;
1813
1814 eherr << "field " << fieldname.c_str () << "cannot be read.";
1815 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1816 }
1817 set_value ((dods_float64 *) &val[0], nelms);
1818 }
1819 break;
1820 default:
1821 release_mod1b_res(reflectance_scales,reflectance_offsets,
1822 radiance_scales,radiance_offsets);
1823 throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
1824 }
1825
1826 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
1827#if 0
1828 if(reflectance_scales!=NULL)
1829 {
1830 delete[] reflectance_offsets;
1831 delete[] reflectance_scales;
1832 }
1833
1834 if(radiance_scales!=NULL)
1835 {
1836 delete[] radiance_offsets;
1837 delete[] radiance_scales;
1838 }
1839#endif
1840
1841 // Somehow the macro RECALCULATE causes the interaction between gridid and sdfileid. SO
1842 // If I close the sdfileid earlier, gridid becomes invalid. So close the sdfileid now. KY 2014-10-24
1843 if (true == isgeofile || false == check_pass_fileid_key)
1844 SDend(sdfileid);
1845 //
1846 return false;
1847
1848}
1849
1850
1851int
1852HDFEOS2Array_RealField::write_dap_data_disable_scale_comp(int32 gridid,
1853 int nelms,
1854 int32 *offset32,
1855 int32*count32,
1856 int32*step32) {
1857
1858
1859 BESDEBUG("h4",
1860 "Coming to HDFEOS2_Array_RealField: write_dap_data_disable_scale_comp"
1861 <<endl);
1862
1863 // Define function pointers to handle both grid and swath
1864 intn (*fieldinfofunc) (int32, char *, int32 *, int32 *, int32 *, char *);
1865 intn (*readfieldfunc) (int32, char *, int32 *, int32 *, int32 *, void *);
1866
1867
1868 if (swathname == "") {
1869 fieldinfofunc = GDfieldinfo;
1870 readfieldfunc = GDreadfield;
1871
1872 }
1873 else if (gridname == "") {
1874 fieldinfofunc = SWfieldinfo;
1875 readfieldfunc = SWreadfield;
1876
1877 }
1878 else
1879 throw InternalErr (__FILE__, __LINE__, "It should be either grid or swath.");
1880
1881
1882 // tmp_rank and tmp_dimlist are two dummy variables
1883 // that are only used when calling fieldinfo.
1884 int32 tmp_rank = 0;
1885 char tmp_dimlist[1024];
1886
1887 // field dimension sizes
1888 int32 tmp_dims[rank];
1889
1890 // field data type
1891 int32 field_dtype = 0;
1892
1893 // returned value of HDF4 and HDF-EOS2 APIs
1894 intn r = 0;
1895
1896 // Obtain the field info. We mainly need the datatype information
1897 // to allocate the buffer to store the data
1898 r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
1899 &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1900 if (r != 0) {
1901 ostringstream eherr;
1902 eherr << "Field " << fieldname.c_str ()
1903 << " information cannot be obtained.";
1904 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1905 }
1906
1907
1908 switch (field_dtype) {
1909 case DFNT_INT8:
1910 {
1911 vector<int8>val;
1912 val.resize(nelms);
1913 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1914 offset32, step32, count32, &val[0]);
1915 if (r != 0) {
1916 ostringstream eherr;
1917 eherr << "field " << fieldname.c_str () << "cannot be read.";
1918 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1919 }
1920
1921#ifndef SIGNED_BYTE_TO_INT32
1922 set_value((dods_byte*)&val[0],nelms);
1923#else
1924
1925 vector<int32>newval;
1926 newval.resize(nelms);
1927
1928 for (int counter = 0; counter < nelms; counter++)
1929 newval[counter] = (int32) (val[counter]);
1930
1931 set_value((dods_int32*)&newval[0],nelms);
1932#endif
1933 }
1934 break;
1935 case DFNT_UINT8:
1936 case DFNT_UCHAR8:
1937 {
1938
1939 vector<uint8>val;
1940 val.resize(nelms);
1941 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1942 offset32, step32, count32, &val[0]);
1943 if (r != 0) {
1944
1945 ostringstream eherr;
1946 eherr << "field " << fieldname.c_str () << "cannot be read.";
1947 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1948 }
1949
1950 set_value((dods_byte*)&val[0],nelms);
1951 }
1952 break;
1953
1954 case DFNT_INT16:
1955 {
1956 vector<int16>val;
1957 val.resize(nelms);
1958 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1959 offset32, step32, count32, &val[0]);
1960
1961 if (r != 0) {
1962 ostringstream eherr;
1963 eherr << "field " << fieldname.c_str () << "cannot be read.";
1964 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1965 }
1966 set_value((dods_int16*)&val[0],nelms);
1967 }
1968 break;
1969 case DFNT_UINT16:
1970 {
1971 vector<uint16>val;
1972 val.resize(nelms);
1973 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1974 offset32, step32, count32, &val[0]);
1975 if (r != 0) {
1976 ostringstream eherr;
1977 eherr << "field " << fieldname.c_str () << "cannot be read.";
1978 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1979 }
1980
1981 set_value((dods_uint16*)&val[0],nelms);
1982 }
1983 break;
1984 case DFNT_INT32:
1985 {
1986 vector<int32>val;
1987 val.resize(nelms);
1988 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1989 offset32, step32, count32, &val[0]);
1990 if (r != 0) {
1991 ostringstream eherr;
1992 eherr << "field " << fieldname.c_str () << "cannot be read.";
1993 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1994 }
1995
1996 set_value((dods_int32*)&val[0],nelms);
1997 }
1998 break;
1999 case DFNT_UINT32:
2000 {
2001 vector<uint32>val;
2002 val.resize(nelms);
2003 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
2004 offset32, step32, count32, &val[0]);
2005 if (r != 0) {
2006 ostringstream eherr;
2007 eherr << "field " << fieldname.c_str () << "cannot be read.";
2008 throw InternalErr (__FILE__, __LINE__, eherr.str ());
2009 }
2010
2011 set_value((dods_uint32*)&val[0],nelms);
2012 }
2013 break;
2014 case DFNT_FLOAT32:
2015 {
2016 vector<float32>val;
2017 val.resize(nelms);
2018 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
2019 offset32, step32, count32, &val[0]);
2020 if (r != 0) {
2021 ostringstream eherr;
2022 eherr << "field " << fieldname.c_str () << "cannot be read.";
2023 throw InternalErr (__FILE__, __LINE__, eherr.str ());
2024 }
2025
2026 // Recalculate seems not necessary.
2027 set_value((dods_float32*)&val[0],nelms);
2028 }
2029 break;
2030 case DFNT_FLOAT64:
2031 {
2032 vector<float64>val;
2033 val.resize(nelms);
2034 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
2035 offset32, step32, count32, &val[0]);
2036 if (r != 0) {
2037 ostringstream eherr;
2038 eherr << "field " << fieldname.c_str () << "cannot be read.";
2039 throw InternalErr (__FILE__, __LINE__, eherr.str ());
2040 }
2041 set_value ((dods_float64 *) &val[0], nelms);
2042 }
2043 break;
2044 default:
2045 throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
2046 }
2047 return 0;
2048}
2049#if 0
2050 r = detachfunc (gridid);
2051 if (r != 0) {
2052 closefunc(gfid);
2053 ostringstream eherr;
2054
2055 eherr << "Grid/Swath " << datasetname.c_str () << " cannot be detached.";
2056 throw InternalErr (__FILE__, __LINE__, eherr.str ());
2057 }
2058
2059
2060 r = closefunc (gfid);
2061 if (r != 0) {
2062 ostringstream eherr;
2063
2064 eherr << "Grid/Swath " << filename.c_str () << " cannot be closed.";
2065 throw InternalErr (__FILE__, __LINE__, eherr.str ());
2066 }
2067
2068 return false;
2069}
2070#endif
2071
2072// Standard way to pass the coordinates of the subsetted region from the client to the handlers
2073// Return the number of elements to read.
2074int
2075HDFEOS2Array_RealField::format_constraint (int *offset, int *step, int *count)
2076{
2077 long nels = 1;
2078 int id = 0;
2079
2080 Dim_iter p = dim_begin ();
2081 while (p != dim_end ()) {
2082
2083 int start = dimension_start (p, true);
2084 int stride = dimension_stride (p, true);
2085 int stop = dimension_stop (p, true);
2086
2087 // Check for illegal constraint
2088 if (start > stop) {
2089 ostringstream oss;
2090 oss << "Array/Grid hyperslab start point "<< start <<
2091 " is greater than stop point " << stop <<".";
2092 throw Error(malformed_expr, oss.str());
2093 }
2094
2095 offset[id] = start;
2096 step[id] = stride;
2097 count[id] = ((stop - start) / stride) + 1; // count of elements
2098 nels *= count[id]; // total number of values for variable
2099
2100 BESDEBUG ("h4",
2101 "=format_constraint():"
2102 << "id=" << id << " offset=" << offset[id]
2103 << " step=" << step[id]
2104 << " count=" << count[id]
2105 << endl);
2106
2107 id++;
2108 p++;
2109 }// while (p != dim_end ())
2110
2111 return nels;
2112}
2113
2114
2115void HDFEOS2Array_RealField::close_fileid(const int gsfileid, const int sdfileid) {
2116
2117#if 0
2118 string check_pass_fileid_key_str="H4.EnablePassFileID";
2119 bool check_pass_fileid_key = false;
2120 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
2121#endif
2122
2123
2124 //if(true == isgeofile || false == check_pass_fileid_key) {
2125 if(true == isgeofile || false == HDF4RequestHandler::get_pass_fileid()) {
2126
2127 if(sdfileid != -1)
2128 SDend(sdfileid);
2129
2130 if(gsfileid != -1){
2131 if(""==gridname)
2132 SWclose(gsfileid);
2133 if (""==swathname)
2134 GDclose(gsfileid);
2135 }
2136
2137 }
2138
2139}
2140
2141void HDFEOS2Array_RealField::release_mod1b_res(float*ref_scale,
2142 float*ref_offset,
2143 float*rad_scale,
2144 float*rad_offset) {
2145
2146 if(ref_scale != NULL)
2147 delete[] ref_scale;
2148 if(ref_offset != NULL)
2149 delete[] ref_offset;
2150 if(rad_scale != NULL)
2151 delete[] rad_scale;
2152 if(rad_offset != NULL)
2153 delete[] rad_offset;
2154
2155}
2156
2157
2158#endif
void close_fileid(hid_t fid)
Definition: h5get.cc:433