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