bes Updated for version 3.21.1
The Backend Server (BES) is the lower two tiers of the Hyrax data server
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: Kent Yang <myang6@hdfgroup.org> Eunsoo Seo
5// Copyright (c) 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!=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 }
394
395 // For swath, we don't see any MODIS 1-D fields that have scale,offset
396 // and fill value attributes that need to be changed.
397 // So now we don't need to handle the vdata handling.
398 // Another reason is the possible change of the implementation
399 // of the SDS attribute handlings. That may be too costly.
400 // KY 2012-07-31
401
402 if ( false == field_is_vdata) {
403
404 char attrname[H4_MAX_NC_NAME + 1];
405 vector<char> attrbuf;
406
407 // Obtain attribute values.
408 if (false == isgeofile || false == check_pass_fileid_key)
409 sdfileid = sdfd;
410 else {
411 sdfileid = SDstart(filename.c_str (), DFACC_READ);
412 if (FAIL == sdfileid) {
413 ostringstream eherr;
414 eherr << "Cannot Start the SD interface for the file "
415 << filename;
416 throw InternalErr (__FILE__, __LINE__, eherr.str ());
417 }
418 }
419
420 int32 sdsindex = -1;
421 int32 sdsid;
422 sdsindex = SDnametoindex(sdfileid, fieldname.c_str());
423 if (FAIL == sdsindex) {
424 if (true == isgeofile || false == check_pass_fileid_key)
425 SDend(sdfileid);
426 ostringstream eherr;
427 eherr << "Cannot obtain the index of " << fieldname;
428 throw InternalErr (__FILE__, __LINE__, eherr.str ());
429 }
430
431 sdsid = SDselect(sdfileid, sdsindex);
432 if (FAIL == sdsid) {
433 if (true == isgeofile || false == check_pass_fileid_key)
434 SDend(sdfileid);
435 ostringstream eherr;
436 eherr << "Cannot obtain the SDS ID of " << fieldname;
437 throw InternalErr (__FILE__, __LINE__, eherr.str ());
438 }
439
440#if 0
441 char attrname[H4_MAX_NC_NAME + 1];
442 vector<char> attrbuf, attrbuf2;
443
444 // Here we cannot check if SDfindattr fails or not since even SDfindattr fails it doesn't mean
445 // errors happen. If no such attribute can be found, SDfindattr still returns FAIL.
446 // The correct way is to use SDgetinfo and SDattrinfo to check if attributes "radiance_scales" etc exist.
447 // For the time being, I won't do this, due to the performance reason and code simplity and also the
448 // very small chance of real FAIL for SDfindattr.
449 cf_general_attrindex = SDfindattr(sdsid, "radiance_scales");
450 cf_general_attrindex2 = SDfindattr(sdsid, "radiance_offsets");
451
452 // Obtain the values of radiance_scales and radiance_offsets if they are available.
453 if (cf_general_attrindex!=FAIL && cf_general_attrindex2!=FAIL)
454 {
455 intn ret = -1;
456 ret = SDattrinfo(sdsid, cf_general_attrindex, attrname, &attr_dtype, &temp_attrcount);
457 if (ret==FAIL)
458 {
459 SDendaccess(sdsid);
460 if (true == isgeofile)
461 SDend(sdfileid);
462 ostringstream eherr;
463 eherr << "Attribute 'radiance_scales' in " << fieldname.c_str () << " cannot be obtained.";
464 throw InternalErr (__FILE__, __LINE__, eherr.str ());
465 }
466 attrbuf.clear();
467 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
468 ret = SDreadattr(sdsid, cf_general_attrindex, (VOIDP)attrbuf.data());
469 if (ret==FAIL)
470 {
471 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
472 SDendaccess(sdsid);
473 if (true == isgeofile)
474 SDend(sdfileid);
475 ostringstream eherr;
476 eherr << "Attribute 'radiance_scales' in " << fieldname.c_str () << " cannot be obtained.";
477 throw InternalErr (__FILE__, __LINE__, eherr.str ());
478 }
479 ret = SDattrinfo(sdsid, cf_general_attrindex2, attrname, &attr_dtype, &temp_attrcount);
480 if (ret==FAIL)
481 {
482 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
483 SDendaccess(sdsid);
484 if (true == isgeofile)
485 SDend(sdfileid);
486 ostringstream eherr;
487 eherr << "Attribute 'radiance_offsets' in " << fieldname.c_str () << " cannot be obtained.";
488 throw InternalErr (__FILE__, __LINE__, eherr.str ());
489 }
490 attrbuf2.clear();
491 attrbuf2.resize(DFKNTsize(attr_dtype)*temp_attrcount);
492 ret = SDreadattr(sdsid, cf_general_attrindex2, (VOIDP)attrbuf2.data());
493 if (ret==FAIL)
494 {
495 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
496 SDendaccess(sdsid);
497 if (true == isgeofile)
498 SDend(sdfileid);
499 ostringstream eherr;
500 eherr << "Attribute 'radiance_offsets' in " << fieldname.c_str () << " cannot be obtained.";
501 throw InternalErr (__FILE__, __LINE__, eherr.str ());
502 }
503
504 // The following macro will obtain radiance_scales and radiance_offsets.
505 // Although the code is compact, it may not be easy to follow. The similar macro can also be found later.
506 switch(attr_dtype)
507 {
508#define GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(TYPE, CAST) \
509 case DFNT_##TYPE: \
510 { \
511 CAST *ptr = (CAST*)attrbuf.data(); \
512 CAST *ptr2 = (CAST*)attrbuf2.data(); \
513 radiance_scales = new float[temp_attrcount]; \
514 radiance_offsets = new float[temp_attrcount]; \
515 for(int l=0; l<temp_attrcount; l++) \
516 { \
517 radiance_scales[l] = ptr[l]; \
518 radiance_offsets[l] = ptr2[l]; \
519 } \
520 } \
521 break;
522 GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT32, float);
523 GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT64, double);
524 }
525#undef GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES
526 num_eles_of_an_attr = temp_attrcount; // Store the count of attributes.
527 }
528
529 // Obtain attribute values of reflectance_scales and reflectance_offsets if they are available.
530 cf_general_attrindex = SDfindattr(sdsid, "reflectance_scales");
531 cf_general_attrindex2 = SDfindattr(sdsid, "reflectance_offsets");
532 if (cf_general_attrindex!=FAIL && cf_general_attrindex2!=FAIL)
533 {
534 intn ret = -1;
535 ret = SDattrinfo(sdsid, cf_general_attrindex, attrname, &attr_dtype, &temp_attrcount);
536 if (ret==FAIL)
537 {
538 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
539 SDendaccess(sdsid);
540 if (true == isgeofile)
541 SDend(sdfileid);
542 ostringstream eherr;
543 eherr << "Attribute 'reflectance_scales' in " << fieldname.c_str () << " cannot be obtained.";
544 throw InternalErr (__FILE__, __LINE__, eherr.str ());
545 }
546 attrbuf.clear();
547 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
548 ret = SDreadattr(sdsid, cf_general_attrindex, (VOIDP)attrbuf.data());
549 if (ret==FAIL)
550 {
551 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
552 SDendaccess(sdsid);
553 if (true == isgeofile)
554 SDend(sdfileid);
555 ostringstream eherr;
556 eherr << "Attribute 'reflectance_scales' in " << fieldname.c_str () << " cannot be obtained.";
557 throw InternalErr (__FILE__, __LINE__, eherr.str ());
558 }
559
560 ret = SDattrinfo(sdsid, cf_general_attrindex2, attrname, &attr_dtype, &temp_attrcount);
561 if (ret==FAIL)
562 {
563 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
564 SDendaccess(sdsid);
565 if (true == isgeofile)
566 SDend(sdfileid);
567 ostringstream eherr;
568 eherr << "Attribute 'reflectance_offsets' in " << fieldname.c_str () << " cannot be obtained.";
569 throw InternalErr (__FILE__, __LINE__, eherr.str ());
570 }
571 attrbuf2.clear();
572 attrbuf2.resize(DFKNTsize(attr_dtype)*temp_attrcount);
573 ret = SDreadattr(sdsid, cf_general_attrindex2, (VOIDP)attrbuf2.data());
574 if (ret==FAIL)
575 {
576 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
577 SDendaccess(sdsid);
578 if (true == isgeofile)
579 SDend(sdfileid);
580 ostringstream eherr;
581 eherr << "Attribute 'reflectance_offsets' in " << fieldname.c_str () << " cannot be obtained.";
582 throw InternalErr (__FILE__, __LINE__, eherr.str ());
583 }
584 switch(attr_dtype)
585 {
586#define GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(TYPE, CAST) \
587 case DFNT_##TYPE: \
588 { \
589 CAST *ptr = (CAST*)attrbuf.data(); \
590 CAST *ptr2 = (CAST*)attrbuf2.data(); \
591 reflectance_scales = new float[temp_attrcount]; \
592 reflectance_offsets = new float[temp_attrcount]; \
593 for(int l=0; l<temp_attrcount; l++) \
594 { \
595 reflectance_scales[l] = ptr[l]; \
596 reflectance_offsets[l] = ptr2[l]; \
597 } \
598 } \
599 break;
600 GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT32, float);
601 GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT64, double);
602 }
603#undef GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES
604 num_eles_of_an_attr = temp_attrcount; // Store the count of attributes.
605 }
606
607#endif
608 // Obtain the value of attribute scale_factor.
609 scale_factor_attr_index = SDfindattr(sdsid, "scale_factor");
610 if (scale_factor_attr_index!=FAIL)
611 {
612 intn ret = -1;
613 ret = SDattrinfo(sdsid, scale_factor_attr_index, attrname,
614 &attr_dtype, &temp_attrcount);
615 if (ret==FAIL)
616 {
617 SDendaccess(sdsid);
618 if (true == isgeofile || false == check_pass_fileid_key)
619 SDend(sdfileid);
620 ostringstream eherr;
621 eherr << "Attribute 'scale_factor' in "
622 << fieldname.c_str () << " cannot be obtained.";
623 throw InternalErr (__FILE__, __LINE__, eherr.str ());
624 }
625 attrbuf.clear();
626 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
627 ret = SDreadattr(sdsid, scale_factor_attr_index, (VOIDP)attrbuf.data());
628 if (ret==FAIL)
629 {
630 SDendaccess(sdsid);
631 if (true == isgeofile || false == check_pass_fileid_key)
632 SDend(sdfileid);
633
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
640 switch(attr_dtype)
641 {
642#define GET_SCALE_FACTOR_ATTR_VALUE(TYPE, CAST) \
643 case DFNT_##TYPE: \
644 { \
645 CAST tmpvalue = *(CAST*)attrbuf.data(); \
646 scale = (float)tmpvalue; \
647 } \
648 break;
649 GET_SCALE_FACTOR_ATTR_VALUE(INT8, int8)
650 GET_SCALE_FACTOR_ATTR_VALUE(CHAR,int8)
651 GET_SCALE_FACTOR_ATTR_VALUE(UINT8, uint8)
652 GET_SCALE_FACTOR_ATTR_VALUE(UCHAR,uint8)
653 GET_SCALE_FACTOR_ATTR_VALUE(INT16, int16)
654 GET_SCALE_FACTOR_ATTR_VALUE(UINT16, uint16)
655 GET_SCALE_FACTOR_ATTR_VALUE(INT32, int32)
656 GET_SCALE_FACTOR_ATTR_VALUE(UINT32, uint32)
657 GET_SCALE_FACTOR_ATTR_VALUE(FLOAT32, float)
658 GET_SCALE_FACTOR_ATTR_VALUE(FLOAT64, double)
659 default:
660 throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
661
662
663 }
664#undef GET_SCALE_FACTOR_ATTR_VALUE
665 }
666
667 // Obtain the value of attribute add_offset
668 add_offset_attr_index = SDfindattr(sdsid, "add_offset");
669 if (add_offset_attr_index!=FAIL)
670 {
671 intn ret;
672 ret = SDattrinfo(sdsid, add_offset_attr_index, attrname,
673 &attr_dtype, &temp_attrcount);
674 if (ret==FAIL)
675 {
676 SDendaccess(sdsid);
677 if (true == isgeofile || false == check_pass_fileid_key)
678 SDend(sdfileid);
679
680 ostringstream eherr;
681 eherr << "Attribute 'add_offset' in " << fieldname.c_str ()
682 << " cannot be obtained.";
683 throw InternalErr (__FILE__, __LINE__, eherr.str ());
684 }
685 attrbuf.clear();
686 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
687 ret = SDreadattr(sdsid, add_offset_attr_index, (VOIDP)attrbuf.data());
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
700 switch(attr_dtype)
701 {
702#define GET_ADD_OFFSET_ATTR_VALUE(TYPE, CAST) \
703 case DFNT_##TYPE: \
704 { \
705 CAST tmpvalue = *(CAST*)attrbuf.data(); \
706 field_offset = (float)tmpvalue; \
707 } \
708 break;
709 GET_ADD_OFFSET_ATTR_VALUE(INT8, int8)
710 GET_ADD_OFFSET_ATTR_VALUE(CHAR,int8)
711 GET_ADD_OFFSET_ATTR_VALUE(UINT8, uint8)
712 GET_ADD_OFFSET_ATTR_VALUE(UCHAR,uint8)
713 GET_ADD_OFFSET_ATTR_VALUE(INT16, int16)
714 GET_ADD_OFFSET_ATTR_VALUE(UINT16, uint16)
715 GET_ADD_OFFSET_ATTR_VALUE(INT32, int32)
716 GET_ADD_OFFSET_ATTR_VALUE(UINT32, uint32)
717 GET_ADD_OFFSET_ATTR_VALUE(FLOAT32, float)
718 GET_ADD_OFFSET_ATTR_VALUE(FLOAT64, double)
719 default:
720 throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
721
722 }
723#undef GET_ADD_OFFSET_ATTR_VALUE
724 }
725
726 // Obtain the value of the attribute _FillValue
727 cf_fv_attrindex = SDfindattr(sdsid, "_FillValue");
728 if (cf_fv_attrindex!=FAIL)
729 {
730 intn ret;
731 ret = SDattrinfo(sdsid, cf_fv_attrindex, attrname, &attr_dtype, &temp_attrcount);
732 if (ret==FAIL)
733 {
734 SDendaccess(sdsid);
735 if (true == isgeofile || false == check_pass_fileid_key)
736 SDend(sdfileid);
737
738 ostringstream eherr;
739 eherr << "Attribute '_FillValue' in " << fieldname.c_str ()
740 << " cannot be obtained.";
741 throw InternalErr (__FILE__, __LINE__, eherr.str ());
742 }
743 attrbuf.clear();
744 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
745 ret = SDreadattr(sdsid, cf_fv_attrindex, (VOIDP)attrbuf.data());
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
758 switch(attr_dtype)
759 {
760#define GET_FILLVALUE_ATTR_VALUE(TYPE, CAST) \
761 case DFNT_##TYPE: \
762 { \
763 CAST tmpvalue = *(CAST*)attrbuf.data(); \
764 fillvalue = (float)tmpvalue; \
765 } \
766 break;
767 GET_FILLVALUE_ATTR_VALUE(INT8, int8)
768 GET_FILLVALUE_ATTR_VALUE(CHAR, int8)
769 GET_FILLVALUE_ATTR_VALUE(INT16, int16)
770 GET_FILLVALUE_ATTR_VALUE(INT32, int32)
771 GET_FILLVALUE_ATTR_VALUE(UINT8, uint8)
772 GET_FILLVALUE_ATTR_VALUE(UCHAR, uint8)
773 GET_FILLVALUE_ATTR_VALUE(UINT16, uint16)
774 GET_FILLVALUE_ATTR_VALUE(UINT32, uint32)
775 default:
776 throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
777
778 }
779#undef GET_FILLVALUE_ATTR_VALUE
780 }
781
782 // Retrieve valid_range,valid_range is normally represented as (valid_min,valid_max)
783 // for non-CF scale and offset rules, the data is always float. So we only
784 // need to change the data type to float.
785 cf_vr_attrindex = SDfindattr(sdsid, "valid_range");
786 if (cf_vr_attrindex!=FAIL)
787 {
788 intn ret;
789 ret = SDattrinfo(sdsid, cf_vr_attrindex, attrname, &attr_dtype, &temp_attrcount);
790 if (ret==FAIL)
791 {
792 SDendaccess(sdsid);
793 if (true == isgeofile || false == check_pass_fileid_key)
794 SDend(sdfileid);
795
796 ostringstream eherr;
797 eherr << "Attribute '_FillValue' in " << fieldname.c_str ()
798 << " cannot be obtained.";
799 throw InternalErr (__FILE__, __LINE__, eherr.str ());
800 }
801 attrbuf.clear();
802 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
803 ret = SDreadattr(sdsid, cf_vr_attrindex, (VOIDP)attrbuf.data());
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
816
817 string attrbuf_str(attrbuf.begin(),attrbuf.end());
818
819 switch(attr_dtype) {
820
821 case DFNT_CHAR:
822 {
823 // We need to treat the attribute data as characters or string.
824 // So find the separator.
825 size_t found = attrbuf_str.find_first_of(",");
826 size_t found_from_end = attrbuf_str.find_last_of(",");
827
828 if (string::npos == found){
829 SDendaccess(sdsid);
830 if (true == isgeofile || false == check_pass_fileid_key)
831 SDend(sdfileid);
832 throw InternalErr(__FILE__,__LINE__,"should find the separator ,");
833 }
834 if (found != found_from_end){
835 SDendaccess(sdsid);
836 if (true == isgeofile || false == check_pass_fileid_key)
837 SDend(sdfileid);
838 throw InternalErr(__FILE__,__LINE__,
839 "Only one separator , should be available.");
840 }
841
842 orig_valid_min = (float)(atof((attrbuf_str.substr(0,found)).c_str()));
843 orig_valid_max = (float)(atof((attrbuf_str.substr(found+1)).c_str()));
844
845 }
846 break;
847 case DFNT_INT8:
848 {
849 // We find a special case that even valid_range is logically
850 // interpreted as two elements,
851 // but the count of attribute elements is more than 2. The count
852 // actually is the number of
853 // characters stored as the attribute value. So we need to find
854 // the separator "," and then
855 // change the string before and after the separator into float numbers.
856 //
857 if (temp_attrcount >2) {
858
859 size_t found = attrbuf_str.find_first_of(",");
860 size_t found_from_end = attrbuf_str.find_last_of(",");
861
862 if (string::npos == found){
863 SDendaccess(sdsid);
864 if (true == isgeofile || false == check_pass_fileid_key)
865 SDend(sdfileid);
866 throw InternalErr(__FILE__,__LINE__,"should find the separator ,");
867 }
868 if (found != found_from_end){
869 SDendaccess(sdsid);
870 if (true == isgeofile || false == check_pass_fileid_key)
871 SDend(sdfileid);
872 throw InternalErr(__FILE__,__LINE__,
873 "Only one separator , should be available.");
874 }
875
876 orig_valid_min = (float)(atof((attrbuf_str.substr(0,found)).c_str()));
877 orig_valid_max = (float)(atof((attrbuf_str.substr(found+1)).c_str()));
878
879 }
880 else if (2 == temp_attrcount) {
881 orig_valid_min = (float)attrbuf[0];
882 orig_valid_max = (float)attrbuf[1];
883 }
884 else {
885 SDendaccess(sdsid);
886 if (true == isgeofile || false == check_pass_fileid_key)
887 SDend(sdfileid);
888 throw InternalErr(__FILE__,__LINE__,
889 "The number of attribute count should be greater than 1.");
890 }
891
892 }
893 break;
894
895 case DFNT_UINT8:
896 case DFNT_UCHAR:
897 {
898 if (temp_attrcount != 2) {
899 SDendaccess(sdsid);
900 if (true == isgeofile || false == check_pass_fileid_key)
901 SDend(sdfileid);
902
903 throw InternalErr(__FILE__,__LINE__,
904 "The number of attribute count should be 2 for the DFNT_UINT8 type.");
905 }
906
907 auto temp_valid_range = (unsigned char *)attrbuf.data();
908 orig_valid_min = (float)(temp_valid_range[0]);
909 orig_valid_max = (float)(temp_valid_range[1]);
910 }
911 break;
912
913 case DFNT_INT16:
914 {
915 if (temp_attrcount != 2) {
916 SDendaccess(sdsid);
917 if (true == isgeofile || false == check_pass_fileid_key)
918 SDend(sdfileid);
919
920 throw InternalErr(__FILE__,__LINE__,
921 "The number of attribute count should be 2 for the DFNT_INT16 type.");
922 }
923
924 auto temp_valid_range = (short *)attrbuf.data();
925 orig_valid_min = (float)(temp_valid_range[0]);
926 orig_valid_max = (float)(temp_valid_range[1]);
927 }
928 break;
929
930 case DFNT_UINT16:
931 {
932 if (temp_attrcount != 2) {
933 SDendaccess(sdsid);
934 if (true == isgeofile || false == check_pass_fileid_key)
935 SDend(sdfileid);
936
937 throw InternalErr(__FILE__,__LINE__,
938 "The number of attribute count should be 2 for the DFNT_UINT16 type.");
939 }
940
941 auto temp_valid_range = (unsigned short *)attrbuf.data();
942 orig_valid_min = (float)(temp_valid_range[0]);
943 orig_valid_max = (float)(temp_valid_range[1]);
944 }
945 break;
946
947 case DFNT_INT32:
948 {
949 if (temp_attrcount != 2) {
950 SDendaccess(sdsid);
951 if (true == isgeofile || false == check_pass_fileid_key)
952 SDend(sdfileid);
953
954 throw InternalErr(__FILE__,__LINE__,
955 "The number of attribute count should be 2 for the DFNT_INT32 type.");
956 }
957
958 auto temp_valid_range = (int *)attrbuf.data();
959 orig_valid_min = (float)(temp_valid_range[0]);
960 orig_valid_max = (float)(temp_valid_range[1]);
961 }
962 break;
963
964 case DFNT_UINT32:
965 {
966 if (temp_attrcount != 2) {
967 SDendaccess(sdsid);
968 if (true == isgeofile || false == check_pass_fileid_key)
969 SDend(sdfileid);
970
971 throw InternalErr(__FILE__,__LINE__,
972 "The number of attribute count should be 2 for the DFNT_UINT32 type.");
973 }
974
975 auto temp_valid_range = (unsigned int *)attrbuf.data();
976 orig_valid_min = (float)(temp_valid_range[0]);
977 orig_valid_max = (float)(temp_valid_range[1]);
978 }
979 break;
980
981 case DFNT_FLOAT32:
982 {
983 if (temp_attrcount != 2) {
984 SDendaccess(sdsid);
985 if (true == isgeofile || false == check_pass_fileid_key)
986 SDend(sdfileid);
987
988 throw InternalErr(__FILE__,__LINE__,
989 "The number of attribute count should be 2 for the DFNT_FLOAT32 type.");
990 }
991
992 auto temp_valid_range = (float *)attrbuf.data();
993 orig_valid_min = temp_valid_range[0];
994 orig_valid_max = temp_valid_range[1];
995 }
996 break;
997
998 case DFNT_FLOAT64:
999 {
1000 if (temp_attrcount != 2){
1001 SDendaccess(sdsid);
1002 if (true == isgeofile || false == check_pass_fileid_key)
1003 SDend(sdfileid);
1004
1005 throw InternalErr(__FILE__,__LINE__,
1006 "The number of attribute count should be 2 for the DFNT_FLOAT64 type.");
1007 }
1008 auto temp_valid_range = (double *)attrbuf.data();
1009
1010 // Notice: this approach will lose precision and possibly overflow.
1011 // Fortunately it is not a problem for MODIS data.
1012 // This part of code may not be called.
1013 // If it is called, mostly the value is within the floating-point range.
1014 // KY 2013-01-29
1015 orig_valid_min = (float)temp_valid_range[0];
1016 orig_valid_max = (float)temp_valid_range[1];
1017 }
1018 break;
1019 default: {
1020 SDendaccess(sdsid);
1021 if (true == isgeofile || false == check_pass_fileid_key)
1022 SDend(sdfileid);
1023 throw InternalErr(__FILE__,__LINE__,"Unsupported data type.");
1024 }
1025 }
1026 }
1027
1028 // Check if the data has the "Key" attribute.
1029 // We found that some NSIDC MODIS data(MOD29) used "Key" to identify some special values.
1030 // To get the values that are within the range identified by the "Key",
1031 // scale offset rules also need to be applied to those values
1032 // outside the original valid range. KY 2013-02-25
1033 int32 cf_mod_key_attrindex = SUCCEED;
1034 cf_mod_key_attrindex = SDfindattr(sdsid, "Key");
1035 if (cf_mod_key_attrindex !=FAIL) {
1036 has_Key_attr = true;
1037 }
1038
1039 attrbuf.clear();
1040 vector<char> attrbuf2;
1041
1042 // Here we cannot check if SDfindattr fails since even SDfindattr fails it doesn't mean
1043 // errors happen. If no such attribute can be found, SDfindattr still returns FAIL.
1044 // The correct way is to use SDgetinfo and SDattrinfo to check if attributes
1045 // "radiance_scales" etc exist.
1046 // For the time being, I won't do this, due to the performance reason and code simplity
1047 // and also the very small chance of real FAIL for SDfindattr.
1048 cf_modl1b_rr_attrindex = SDfindattr(sdsid, "radiance_scales");
1049 cf_modl1b_rr_attrindex2 = SDfindattr(sdsid, "radiance_offsets");
1050
1051 // Obtain the values of radiance_scales and radiance_offsets if they are available.
1052 if (cf_modl1b_rr_attrindex!=FAIL && cf_modl1b_rr_attrindex2!=FAIL)
1053 {
1054 intn ret = -1;
1055 ret = SDattrinfo(sdsid, cf_modl1b_rr_attrindex, attrname,
1056 &attr_dtype, &temp_attrcount);
1057 if (ret==FAIL)
1058 {
1059 SDendaccess(sdsid);
1060 if (true == isgeofile || false == check_pass_fileid_key)
1061 SDend(sdfileid);
1062 ostringstream eherr;
1063 eherr << "Attribute 'radiance_scales' in " << fieldname.c_str ()
1064 << " cannot be obtained.";
1065 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1066 }
1067 attrbuf.clear();
1068 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
1069 ret = SDreadattr(sdsid, cf_modl1b_rr_attrindex, (VOIDP)attrbuf.data());
1070 if (ret==FAIL)
1071 {
1072 SDendaccess(sdsid);
1073 if (true == isgeofile || false == check_pass_fileid_key)
1074 SDend(sdfileid);
1075 ostringstream eherr;
1076 eherr << "Attribute 'radiance_scales' in " << fieldname.c_str ()
1077 << " cannot be obtained.";
1078 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1079 }
1080 ret = SDattrinfo(sdsid, cf_modl1b_rr_attrindex2, attrname,
1081 &attr_dtype, &temp_attrcount);
1082 if (ret==FAIL)
1083 {
1084 SDendaccess(sdsid);
1085 if (true == isgeofile || false == check_pass_fileid_key)
1086 SDend(sdfileid);
1087 ostringstream eherr;
1088 eherr << "Attribute 'radiance_offsets' in "
1089 << fieldname.c_str () << " cannot be obtained.";
1090 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1091 }
1092 attrbuf2.clear();
1093 attrbuf2.resize(DFKNTsize(attr_dtype)*temp_attrcount);
1094 ret = SDreadattr(sdsid, cf_modl1b_rr_attrindex2, (VOIDP)attrbuf2.data());
1095 if (ret==FAIL)
1096 {
1097 SDendaccess(sdsid);
1098 if (true == isgeofile || false == check_pass_fileid_key)
1099 SDend(sdfileid);
1100 ostringstream eherr;
1101 eherr << "Attribute 'radiance_offsets' in "
1102 << fieldname.c_str () << " cannot be obtained.";
1103 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1104 }
1105
1106 // The following macro will obtain radiance_scales and radiance_offsets.
1107 // Although the code is compact, it may not be easy to follow.
1108 // The similar macro can also be found later.
1109 switch(attr_dtype)
1110 {
1111#define GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(TYPE, CAST) \
1112 case DFNT_##TYPE: \
1113 { \
1114 CAST *ptr = (CAST*)attrbuf.data(); \
1115 CAST *ptr2 = (CAST*)attrbuf2.data(); \
1116 radiance_scales = new float[temp_attrcount]; \
1117 radiance_offsets = new float[temp_attrcount]; \
1118 for(int l=0; l<temp_attrcount; l++) \
1119 { \
1120 radiance_scales[l] = ptr[l]; \
1121 radiance_offsets[l] = ptr2[l]; \
1122 } \
1123 } \
1124 break;
1125 GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT32, float)
1126 GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT64, double)
1127 default:
1128 throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
1129
1130 }
1131#undef GET_RADIANCE_SCALES_OFFSETS_ATTR_VALUES
1132 // Store the count of attributes.
1133 num_eles_of_an_attr = temp_attrcount;
1134 }
1135
1136 // Obtain attribute values of reflectance_scales
1137 // and reflectance_offsets if they are available.
1138 cf_modl1b_rr_attrindex = SDfindattr(sdsid, "reflectance_scales");
1139 cf_modl1b_rr_attrindex2 = SDfindattr(sdsid, "reflectance_offsets");
1140 if (cf_modl1b_rr_attrindex!=FAIL && cf_modl1b_rr_attrindex2!=FAIL)
1141 {
1142 intn ret = -1;
1143 ret = SDattrinfo(sdsid, cf_modl1b_rr_attrindex, attrname,
1144 &attr_dtype, &temp_attrcount);
1145 if (ret==FAIL)
1146 {
1147 release_mod1b_res(reflectance_scales,reflectance_offsets,
1148 radiance_scales,radiance_offsets);
1149 SDendaccess(sdsid);
1150 if (true == isgeofile || false == check_pass_fileid_key)
1151 SDend(sdfileid);
1152 ostringstream eherr;
1153 eherr << "Attribute 'reflectance_scales' in "
1154 << fieldname.c_str () << " cannot be obtained.";
1155 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1156 }
1157 attrbuf.clear();
1158 attrbuf.resize(DFKNTsize(attr_dtype)*temp_attrcount);
1159 ret = SDreadattr(sdsid, cf_modl1b_rr_attrindex, (VOIDP)attrbuf.data());
1160 if (ret==FAIL)
1161 {
1162 release_mod1b_res(reflectance_scales,reflectance_offsets,
1163 radiance_scales,radiance_offsets);
1164 SDendaccess(sdsid);
1165 if (true == isgeofile || false == check_pass_fileid_key)
1166 SDend(sdfileid);
1167 ostringstream eherr;
1168 eherr << "Attribute 'reflectance_scales' in "
1169 << fieldname.c_str () << " cannot be obtained.";
1170 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1171 }
1172
1173 ret = SDattrinfo(sdsid, cf_modl1b_rr_attrindex2, attrname,
1174 &attr_dtype, &temp_attrcount);
1175 if (ret==FAIL)
1176 {
1177 release_mod1b_res(reflectance_scales,reflectance_offsets,
1178 radiance_scales,radiance_offsets);
1179 SDendaccess(sdsid);
1180 if (true == isgeofile || false == check_pass_fileid_key)
1181 SDend(sdfileid);
1182 ostringstream eherr;
1183 eherr << "Attribute 'reflectance_offsets' in "
1184 << fieldname.c_str () << " cannot be obtained.";
1185 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1186 }
1187 attrbuf2.clear();
1188 attrbuf2.resize(DFKNTsize(attr_dtype)*temp_attrcount);
1189 ret = SDreadattr(sdsid, cf_modl1b_rr_attrindex2, (VOIDP)attrbuf2.data());
1190 if (ret==FAIL)
1191 {
1192 release_mod1b_res(reflectance_scales,reflectance_offsets,
1193 radiance_scales,radiance_offsets);
1194 SDendaccess(sdsid);
1195 if (true == isgeofile || false == check_pass_fileid_key)
1196 SDend(sdfileid);
1197 ostringstream eherr;
1198 eherr << "Attribute 'reflectance_offsets' in "
1199 << fieldname.c_str () << " cannot be obtained.";
1200 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1201 }
1202 switch(attr_dtype)
1203 {
1204#define GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(TYPE, CAST) \
1205 case DFNT_##TYPE: \
1206 { \
1207 CAST *ptr = (CAST*)attrbuf.data(); \
1208 CAST *ptr2 = (CAST*)attrbuf2.data(); \
1209 reflectance_scales = new float[temp_attrcount]; \
1210 reflectance_offsets = new float[temp_attrcount]; \
1211 for(int l=0; l<temp_attrcount; l++) \
1212 { \
1213 reflectance_scales[l] = ptr[l]; \
1214 reflectance_offsets[l] = ptr2[l]; \
1215 } \
1216 } \
1217 break;
1218 GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT32, float)
1219 GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES(FLOAT64, double)
1220 default:
1221 throw InternalErr(__FILE__,__LINE__,"unsupported data type.");
1222
1223 }
1224#undef GET_REFLECTANCE_SCALES_OFFSETS_ATTR_VALUES
1225 num_eles_of_an_attr = temp_attrcount; // Store the count of attributes.
1226 }
1227
1228 SDendaccess(sdsid);
1229
1230 BESDEBUG("h4","scale is "<<scale <<endl);
1231 BESDEBUG("h4","offset is "<<field_offset <<endl);
1232 BESDEBUG("h4","fillvalue is "<<fillvalue <<endl);
1233 }
1234 }
1235
1236 // According to our observations, it seems that MODIS products ALWAYS
1237 // use the "scale" factor to make bigger values smaller.
1238 // So for MODIS_MUL_SCALE products, if the scale of some variable is greater than 1,
1239 // it means that for this variable, the MODIS type for this variable may be MODIS_DIV_SCALE.
1240 // For the similar logic, we may need to change MODIS_DIV_SCALE to MODIS_MUL_SCALE
1241 // and MODIS_EQ_SCALE to MODIS_DIV_SCALE.
1242 // We indeed find such a case. HDF-EOS2 Grid MODIS_Grid_1km_2D of MOD(or MYD)09GA is
1243 // a MODIS_EQ_SCALE.
1244 // However,
1245 // the scale_factor of the variable Range_1 in the MOD09GA product is 25.
1246 // According to our observation,
1247 // this variable should be MODIS_DIV_SCALE.We verify this is true according to
1248 // MODIS 09 product document
1249 // http://modis-sr.ltdri.org/products/MOD09_UserGuide_v1_3.pdf.
1250 // Since this conclusion is based on our observation, we would like to add a BESlog to detect
1251 // if we find
1252 // the similar cases so that we can verify with the corresponding product documents to see if
1253 // this is true.
1254 // More updated information,
1255 // We just verified with the MOD09 data producer, the scale_factor for Range_1 is 25
1256 // but the equation is still multiplication instead of division.
1257 // So we have to make this as a special case and don't change the scale and offset settings
1258 // for Range_1 of MOD09 products.
1259 // KY 2014-01-13
1260
1261 if (SOType::MODIS_EQ_SCALE == sotype || SOType::MODIS_MUL_SCALE == sotype) {
1262 if (scale > 1) {
1263 bool need_change_scale = true;
1264 if (gridname!="") {
1265
1266 string temp_filename;
1267 if (filename.find("#") != string::npos)
1268 temp_filename =filename.substr(filename.find_last_of("#") + 1);
1269 else
1270 temp_filename = filename.substr(filename.find_last_of("/") +1);
1271
1272 if ((temp_filename.size() >5) && ((temp_filename.compare(0,5,"MOD09") == 0)
1273 ||(temp_filename.compare(0,5,"MYD09") == 0))) {
1274 if ((fieldname.size() >5) && fieldname.compare(0,5,"Range") == 0)
1275 need_change_scale = false;
1276 }
1277 // MOD16A2
1278 else if ((temp_filename.size() >7)&&
1279 ((temp_filename.compare(0,7,"MOD16A2") == 0)|| (temp_filename.compare(0,7,"MYD16A2")==0)||
1280 (temp_filename.compare(0,7,"MOD16A3") == 0)|| (temp_filename.compare(0,7,"MYD16A3")==0)))
1281 need_change_scale = false;
1282
1283
1284 }
1285 if (true == need_change_scale) {
1286 sotype = SOType::MODIS_DIV_SCALE;
1287 INFO_LOG( "The field " + fieldname + " scale factor is "
1288 + std::to_string(scale) + " . But the original scale factor type is MODIS_MUL_SCALE"
1289 + " or MODIS_EQ_SCALE. Now change it to MODIS_DIV_SCALE.");
1290 }
1291 }
1292 }
1293
1294 if (SOType::MODIS_DIV_SCALE == sotype) {
1295 if (scale < 1) {
1296 sotype = SOType::MODIS_MUL_SCALE;
1297 INFO_LOG("The field " + fieldname + " scale factor is "
1298 + std::to_string(scale) + ". But the original scale factor type is MODIS_DIV_SCALE. "
1299 + " Now change it to MODIS_MUL_SCALE.");
1300 }
1301 }
1302
1303// We need to loop through all datatpes to allocate the memory buffer for the data.
1304// It is hard to add comments to the macro. We may need to change them to general
1305// routines in the future.
1306// Some MODIS products use both valid_range(valid_min, valid_max) and fillvalues for data fields.
1307// When do recalculating,
1308// I check fillvalue first, then check valid_min and valid_max if they are available.
1309// The middle check is_special_value addresses the MODIS L1B special value.
1310// Updated: just find that the RECALCULATE will be done only when the valid_range
1311// attribute is present(if cf_vr_attrindex!=FAIL).
1312// This restriction is in theory not necessary, but for more MODIS data products,
1313// this restriction may be valid since valid_range pairs with scale/offset to identify
1314// the valid data values. KY 2014-02-19
1315//
1316#if 0
1317/* if ((float)tmptr[l] != fillvalue ) \
1318// { \
1319// f(false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[l]))\
1320// { \
1321// if (orig_valid_min<tmpval[l] && orig_valid_max>tmpval[l] \
1322// if (sotype==MODIS_MUL_SCALE) \
1323// tmpval[l] = (tmptr[l]-field_offset)*scale; \
1324// else if (sotype==MODIS_EQ_SCALE) \
1325// tmpval[l] = tmptr[l]*scale + field_offset; \
1326// else if (sotype==MODIS_DIV_SCALE) \
1327// tmpval[l] = (tmptr[l]-field_offset)/scale; \
1328// } \
1329*/
1330#endif
1331#define RECALCULATE(CAST, DODS_CAST, VAL) \
1332{ \
1333 bool change_data_value = false; \
1334 if (sotype!=SOType::DEFAULT_CF_EQU) \
1335 { \
1336 vector<float>tmpval; \
1337 tmpval.resize(nelms); \
1338 CAST tmptr = (CAST)VAL; \
1339 for(int l=0; l<nelms; l++) \
1340 tmpval[l] = (float)tmptr[l]; \
1341 bool special_case = false; \
1342 if (scale_factor_attr_index==FAIL) \
1343 if (num_eles_of_an_attr==1) \
1344 if ((radiance_scales!=nullptr) && (radiance_offsets!=nullptr)) \
1345 { \
1346 scale = radiance_scales[0]; \
1347 field_offset = radiance_offsets[0];\
1348 special_case = true; \
1349 } \
1350 if (((scale_factor_attr_index!=FAIL) && !((scale==1) && (field_offset==0))) || special_case) \
1351 { \
1352 float temp_scale = scale; \
1353 float temp_offset = field_offset; \
1354 if (sotype==SOType::MODIS_MUL_SCALE) \
1355 temp_offset = -1. *field_offset*temp_scale;\
1356 else if (sotype==SOType::MODIS_DIV_SCALE) \
1357 {\
1358 temp_scale = 1/scale; \
1359 temp_offset = -1. *field_offset*temp_scale;\
1360 }\
1361 for(int l=0; l<nelms; l++) \
1362 { \
1363 if (cf_vr_attrindex!=FAIL) \
1364 { \
1365 if ((float)tmptr[l] != fillvalue ) \
1366 { \
1367 if (false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[l]))\
1368 { \
1369 if (((orig_valid_min<=tmpval[l]) && (orig_valid_max>=tmpval[l])) || (true==has_Key_attr))\
1370 { \
1371 tmpval[l] = tmptr[l]*temp_scale + temp_offset; \
1372 } \
1373 } \
1374 } \
1375 } \
1376 } \
1377 change_data_value = true; \
1378 set_value((dods_float32 *)tmpval.data(), nelms); \
1379 } else if ((num_eles_of_an_attr>1) && (((radiance_scales!=nullptr) && (radiance_offsets!=nullptr)) || ((reflectance_scales!=nullptr) && (reflectance_offsets!=nullptr)))) \
1380 { \
1381 size_t dimindex=0; \
1382 if ( num_eles_of_an_attr!=tmp_dims[dimindex]) \
1383 { \
1384 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets); \
1385 ostringstream eherr; \
1386 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."; \
1387 throw InternalErr (__FILE__, __LINE__, eherr.str ()); \
1388 } \
1389 size_t start_index, end_index; \
1390 size_t nr_elems = nelms/count32[dimindex]; \
1391 start_index = offset32[dimindex]; \
1392 end_index = start_index+step32[dimindex]*(count32[dimindex]-1); \
1393 size_t index = 0;\
1394 for(size_t k=start_index; k<=end_index; k+=step32[dimindex]) \
1395 { \
1396 float tmpscale = (fieldname.find("Emissive")!=string::npos)? radiance_scales[k]: reflectance_scales[k]; \
1397 float tmpoffset = (fieldname.find("Emissive")!=string::npos)? radiance_offsets[k]: reflectance_offsets[k]; \
1398 for(size_t l=0; l<nr_elems; l++) \
1399 { \
1400 if (cf_vr_attrindex!=FAIL) \
1401 { \
1402 if (((float)tmptr[index])!=fillvalue) \
1403 { \
1404 if (false == HDFCFUtil::is_special_value(field_dtype,fillvalue,tmptr[index]))\
1405 { \
1406 if (sotype==SOType::MODIS_MUL_SCALE) \
1407 tmpval[index] = (tmptr[index]-tmpoffset)*tmpscale; \
1408 else if (sotype==SOType::MODIS_EQ_SCALE) \
1409 tmpval[index] = tmptr[index]*tmpscale+tmpoffset; \
1410 else if (sotype==SOType::MODIS_DIV_SCALE) \
1411 tmpval[index] = (tmptr[index]-tmpoffset)/tmpscale; \
1412 } \
1413 } \
1414 } \
1415 index++; \
1416 } \
1417 } \
1418 change_data_value = true; \
1419 set_value((dods_float32 *)tmpval.data(), nelms); \
1420 } \
1421 } \
1422 if (!change_data_value) \
1423 { \
1424 set_value ((DODS_CAST)VAL, nelms); \
1425 } \
1426}
1427 switch (field_dtype) {
1428 case DFNT_INT8:
1429 {
1430
1431 vector<int8>val;
1432 val.resize(nelms);
1433 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1434 offset32.data(), step32.data(), count32.data(), val.data());
1435 if (r != 0) {
1436 release_mod1b_res(reflectance_scales,reflectance_offsets,
1437 radiance_scales,radiance_offsets);
1438 ostringstream eherr;
1439 eherr << "field " << fieldname.c_str () << "cannot be read.";
1440 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1441 }
1442
1443#ifndef SIGNED_BYTE_TO_INT32
1444 RECALCULATE(int8*, dods_byte*, val.data());
1445#else
1446
1447 vector<int32>newval;
1448 newval.resize(nelms);
1449
1450 for (int counter = 0; counter < nelms; counter++)
1451 newval[counter] = (int32) (val[counter]);
1452
1453 RECALCULATE(int32*, dods_int32*, newval.data())
1454#endif
1455 }
1456 break;
1457 case DFNT_UINT8:
1458 case DFNT_UCHAR8:
1459 {
1460
1461 vector<uint8>val;
1462 val.resize(nelms);
1463 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1464 offset32.data(), step32.data(), count32.data(), val.data());
1465 if (r != 0) {
1466 release_mod1b_res(reflectance_scales,reflectance_offsets,
1467 radiance_scales,radiance_offsets);
1468 ostringstream eherr;
1469
1470 eherr << "field " << fieldname.c_str () << "cannot be read.";
1471 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1472 }
1473
1474 RECALCULATE(uint8*, dods_byte*, val.data())
1475 }
1476 break;
1477
1478 case DFNT_INT16:
1479 {
1480 vector<int16>val;
1481 val.resize(nelms);
1482 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1483 offset32.data(), step32.data(), count32.data(), val.data());
1484
1485 if (r != 0) {
1486
1487 release_mod1b_res(reflectance_scales,reflectance_offsets,
1488 radiance_scales,radiance_offsets);
1489 ostringstream eherr;
1490
1491 eherr << "field " << fieldname.c_str () << "cannot be read.";
1492 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1493 }
1494 RECALCULATE(int16*, dods_int16*, val.data())
1495 }
1496 break;
1497 case DFNT_UINT16:
1498 {
1499 vector<uint16>val;
1500 val.resize(nelms);
1501
1502 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1503 offset32.data(), step32.data(), count32.data(), val.data());
1504 if (r != 0) {
1505 release_mod1b_res(reflectance_scales,reflectance_offsets,
1506 radiance_scales,radiance_offsets);
1507 ostringstream eherr;
1508
1509 eherr << "field " << fieldname.c_str () << "cannot be read.";
1510 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1511 }
1512
1513 RECALCULATE(uint16*, dods_uint16*, val.data())
1514 }
1515 break;
1516 case DFNT_INT32:
1517 {
1518 vector<int32>val;
1519 val.resize(nelms);
1520 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1521 offset32.data(), step32.data(), count32.data(), val.data());
1522 if (r != 0) {
1523
1524 release_mod1b_res(reflectance_scales,reflectance_offsets,
1525 radiance_scales,radiance_offsets);
1526 ostringstream eherr;
1527
1528 eherr << "field " << fieldname.c_str () << "cannot be read.";
1529 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1530 }
1531
1532 RECALCULATE(int32*, dods_int32*, val.data())
1533 }
1534 break;
1535 case DFNT_UINT32:
1536 {
1537 vector<uint32>val;
1538 val.resize(nelms);
1539 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1540 offset32.data(), step32.data(), count32.data(), val.data());
1541 if (r != 0) {
1542
1543 release_mod1b_res(reflectance_scales,reflectance_offsets,
1544 radiance_scales,radiance_offsets);
1545 ostringstream eherr;
1546
1547 eherr << "field " << fieldname.c_str () << "cannot be read.";
1548 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1549 }
1550
1551 RECALCULATE(uint32*, dods_uint32*, val.data())
1552 }
1553 break;
1554 case DFNT_FLOAT32:
1555 {
1556 vector<float32>val;
1557 val.resize(nelms);
1558 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1559 offset32.data(), step32.data(), count32.data(), val.data());
1560 if (r != 0) {
1561
1562 release_mod1b_res(reflectance_scales,reflectance_offsets,
1563 radiance_scales,radiance_offsets);
1564 ostringstream eherr;
1565
1566 eherr << "field " << fieldname.c_str () << "cannot be read.";
1567 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1568 }
1569
1570 // Recalculate seems not necessary.
1571 RECALCULATE(float32*, dods_float32*, val.data())
1572 //set_value((dods_float32*)val.data(),nelms);
1573 }
1574 break;
1575 case DFNT_FLOAT64:
1576 {
1577 vector<float64>val;
1578 val.resize(nelms);
1579 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1580 offset32.data(), step32.data(), count32.data(), val.data());
1581 if (r != 0) {
1582
1583 release_mod1b_res(reflectance_scales,reflectance_offsets,
1584 radiance_scales,radiance_offsets);
1585 ostringstream eherr;
1586
1587 eherr << "field " << fieldname.c_str () << "cannot be read.";
1588 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1589 }
1590 set_value ((dods_float64 *) val.data(), nelms);
1591 }
1592 break;
1593 default:
1594 release_mod1b_res(reflectance_scales,reflectance_offsets,
1595 radiance_scales,radiance_offsets);
1596 throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
1597 }
1598
1599 release_mod1b_res(reflectance_scales,reflectance_offsets,radiance_scales,radiance_offsets);
1600
1601 // Somehow the macro RECALCULATE causes the interaction between gridid and sdfileid. SO
1602 // If I close the sdfileid earlier, gridid becomes invalid. So close the sdfileid now. KY 2014-10-24
1603 if (true == isgeofile || false == check_pass_fileid_key)
1604 SDend(sdfileid);
1605 //
1606 return false;
1607
1608}
1609
1610
1611int
1612HDFEOS2Array_RealField::write_dap_data_disable_scale_comp(int32 gridid,
1613 int nelms,
1614 int32 *offset32,
1615 int32*count32,
1616 int32*step32) {
1617
1618
1619 BESDEBUG("h4",
1620 "Coming to HDFEOS2_Array_RealField: write_dap_data_disable_scale_comp"
1621 <<endl);
1622
1623 // Define function pointers to handle both grid and swath
1624 intn (*fieldinfofunc) (int32, char *, int32 *, int32 *, int32 *, char *);
1625 intn (*readfieldfunc) (int32, char *, int32 *, int32 *, int32 *, void *);
1626
1627
1628 if (swathname == "") {
1629 fieldinfofunc = GDfieldinfo;
1630 readfieldfunc = GDreadfield;
1631
1632 }
1633 else if (gridname == "") {
1634 fieldinfofunc = SWfieldinfo;
1635 readfieldfunc = SWreadfield;
1636
1637 }
1638 else
1639 throw InternalErr (__FILE__, __LINE__, "It should be either grid or swath.");
1640
1641
1642 // tmp_rank and tmp_dimlist are two dummy variables
1643 // that are only used when calling fieldinfo.
1644 int32 tmp_rank = 0;
1645 char tmp_dimlist[1024];
1646
1647 // field dimension sizes
1648 int32 tmp_dims[rank];
1649
1650 // field data type
1651 int32 field_dtype = 0;
1652
1653 // returned value of HDF4 and HDF-EOS2 APIs
1654 intn r = 0;
1655
1656 // Obtain the field info. We mainly need the datatype information
1657 // to allocate the buffer to store the data
1658 r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
1659 &tmp_rank, tmp_dims, &field_dtype, tmp_dimlist);
1660 if (r != 0) {
1661 ostringstream eherr;
1662 eherr << "Field " << fieldname.c_str ()
1663 << " information cannot be obtained.";
1664 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1665 }
1666
1667
1668 switch (field_dtype) {
1669 case DFNT_INT8:
1670 {
1671 vector<int8>val;
1672 val.resize(nelms);
1673 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1674 offset32, step32, count32, val.data());
1675 if (r != 0) {
1676 ostringstream eherr;
1677 eherr << "field " << fieldname.c_str () << "cannot be read.";
1678 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1679 }
1680
1681#ifndef SIGNED_BYTE_TO_INT32
1682 set_value((dods_byte*)val.data(),nelms);
1683#else
1684
1685 vector<int32>newval;
1686 newval.resize(nelms);
1687
1688 for (int counter = 0; counter < nelms; counter++)
1689 newval[counter] = (int32) (val[counter]);
1690
1691 set_value((dods_int32*)newval.data(),nelms);
1692#endif
1693 }
1694 break;
1695 case DFNT_UINT8:
1696 case DFNT_UCHAR8:
1697 {
1698
1699 vector<uint8>val;
1700 val.resize(nelms);
1701 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1702 offset32, step32, count32, val.data());
1703 if (r != 0) {
1704
1705 ostringstream eherr;
1706 eherr << "field " << fieldname.c_str () << "cannot be read.";
1707 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1708 }
1709
1710 set_value((dods_byte*)val.data(),nelms);
1711 }
1712 break;
1713
1714 case DFNT_INT16:
1715 {
1716 vector<int16>val;
1717 val.resize(nelms);
1718 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1719 offset32, step32, count32, val.data());
1720
1721 if (r != 0) {
1722 ostringstream eherr;
1723 eherr << "field " << fieldname.c_str () << "cannot be read.";
1724 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1725 }
1726 set_value((dods_int16*)val.data(),nelms);
1727 }
1728 break;
1729 case DFNT_UINT16:
1730 {
1731 vector<uint16>val;
1732 val.resize(nelms);
1733 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1734 offset32, step32, count32, val.data());
1735 if (r != 0) {
1736 ostringstream eherr;
1737 eherr << "field " << fieldname.c_str () << "cannot be read.";
1738 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1739 }
1740
1741 set_value((dods_uint16*)val.data(),nelms);
1742 }
1743 break;
1744 case DFNT_INT32:
1745 {
1746 vector<int32>val;
1747 val.resize(nelms);
1748 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1749 offset32, step32, count32, val.data());
1750 if (r != 0) {
1751 ostringstream eherr;
1752 eherr << "field " << fieldname.c_str () << "cannot be read.";
1753 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1754 }
1755
1756 set_value((dods_int32*)val.data(),nelms);
1757 }
1758 break;
1759 case DFNT_UINT32:
1760 {
1761 vector<uint32>val;
1762 val.resize(nelms);
1763 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1764 offset32, step32, count32, val.data());
1765 if (r != 0) {
1766 ostringstream eherr;
1767 eherr << "field " << fieldname.c_str () << "cannot be read.";
1768 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1769 }
1770
1771 set_value((dods_uint32*)val.data(),nelms);
1772 }
1773 break;
1774 case DFNT_FLOAT32:
1775 {
1776 vector<float32>val;
1777 val.resize(nelms);
1778 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1779 offset32, step32, count32, val.data());
1780 if (r != 0) {
1781 ostringstream eherr;
1782 eherr << "field " << fieldname.c_str () << "cannot be read.";
1783 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1784 }
1785
1786 // Recalculate seems not necessary.
1787 set_value((dods_float32*)val.data(),nelms);
1788 }
1789 break;
1790 case DFNT_FLOAT64:
1791 {
1792 vector<float64>val;
1793 val.resize(nelms);
1794 r = readfieldfunc (gridid, const_cast < char *>(fieldname.c_str ()),
1795 offset32, step32, count32, val.data());
1796 if (r != 0) {
1797 ostringstream eherr;
1798 eherr << "field " << fieldname.c_str () << "cannot be read.";
1799 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1800 }
1801 set_value ((dods_float64 *) val.data(), nelms);
1802 }
1803 break;
1804 default:
1805 throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
1806 }
1807 return 0;
1808}
1809
1810// Standard way to pass the coordinates of the subsetted region from the client to the handlers
1811// Return the number of elements to read.
1812int
1813HDFEOS2Array_RealField::format_constraint (int *offset, int *step, int *count)
1814{
1815 int nels = 1;
1816 int id = 0;
1817
1818 Dim_iter p = dim_begin ();
1819 while (p != dim_end ()) {
1820
1821 int start = dimension_start (p, true);
1822 int stride = dimension_stride (p, true);
1823 int stop = dimension_stop (p, true);
1824
1825 // Check for illegal constraint
1826 if (start > stop) {
1827 ostringstream oss;
1828 oss << "Array/Grid hyperslab start point "<< start <<
1829 " is greater than stop point " << stop <<".";
1830 throw Error(malformed_expr, oss.str());
1831 }
1832
1833 offset[id] = start;
1834 step[id] = stride;
1835 count[id] = ((stop - start) / stride) + 1; // count of elements
1836 nels *= count[id]; // total number of values for variable
1837
1838 BESDEBUG ("h4",
1839 "=format_constraint():"
1840 << "id=" << id << " offset=" << offset[id]
1841 << " step=" << step[id]
1842 << " count=" << count[id]
1843 << endl);
1844
1845 id++;
1846 p++;
1847 }
1848
1849 return nels;
1850}
1851
1852
1853void HDFEOS2Array_RealField::close_fileid(const int gsfileid, const int sdfileid) const{
1854
1855 if (true == isgeofile || false == HDF4RequestHandler::get_pass_fileid()) {
1856
1857 if (sdfileid != -1)
1858 SDend(sdfileid);
1859
1860 if (gsfileid != -1){
1861 if (""==gridname)
1862 SWclose(gsfileid);
1863 if (""==swathname)
1864 GDclose(gsfileid);
1865 }
1866
1867 }
1868
1869}
1870
1871void HDFEOS2Array_RealField::release_mod1b_res(float*ref_scale,
1872 float*ref_offset,
1873 float*rad_scale,
1874 float*rad_offset) {
1875
1876 if (ref_scale != nullptr)
1877 delete[] ref_scale;
1878 if (ref_offset != nullptr)
1879 delete[] ref_offset;
1880 if (rad_scale != nullptr)
1881 delete[] rad_scale;
1882 if (rad_offset != nullptr)
1883 delete[] rad_offset;
1884
1885}
1886
1887
1888#endif
ostringstream
STL class.
vector
STL class.