bes Updated for version 3.20.13
HDFSPArrayGeoField.cc
1
2// This file is part of the hdf4 data handler for the OPeNDAP data server.
3// It retrieves the latitude and longitude fields for some special HDF4 data products.
4// The products include TRMML2_V6,TRMML3B_V6,CER_AVG,CER_ES4,CER_CDAY,CER_CGEO,CER_SRB,CER_SYN,CER_ZAVG,OBPGL2,OBPGL3
5// To know more information about these products,check HDFSP.h.
6// Each product stores lat/lon in different way, so we have to retrieve them differently.
7// Authors: MuQun Yang <myang6@hdfgroup.org>
8// Copyright (c) 2010-2012 The HDF Group
10
11#include "HDFSPArrayGeoField.h"
12#include <iostream>
13#include <sstream>
14#include <cassert>
15#include <libdap/debug.h>
16#include "hdf.h"
17#include "mfhdf.h"
18#include <libdap/InternalErr.h>
19#include <BESDebug.h>
20#include "HDFCFUtil.h"
21#include "HDF4RequestHandler.h"
22
23using namespace std;
24using namespace libdap;
25#define SIGNED_BYTE_TO_INT32 1
26
27// The following macros provide names of latitude and longitude for specific HDF4 products.
28#define NUM_PIXEL_NAME "Pixels per Scan Line"
29#define NUM_POINTS_LINE_NAME "Number of Pixel Control Points"
30#define NUM_SCAN_LINE_NAME "Number of Scan Lines"
31#define NUM_LAT_NAME "Number of Lines"
32#define NUM_LON_NAME "Number of Columns"
33#define LAT_STEP_NAME "Latitude Step"
34#define LON_STEP_NAME "Longitude Step"
35#define SWP_LAT_NAME "SW Point Latitude"
36#define SWP_LON_NAME "SW Point Longitude"
37
38
39bool HDFSPArrayGeoField::read ()
40{
41
42 BESDEBUG("h4","Coming to HDFSPArrayGeoField read "<<endl);
43
44 if(length() == 0)
45 return true;
46
47 // Declare offset, count and step
48 vector<int> offset;
49 offset.resize(rank);
50 vector<int> count;
51 count.resize(rank);
52 vector<int> step;
53 step.resize(rank);
54
55 // Obtain offset, step and count from the client expression constraint
56 int nelms = -1;
57 nelms = format_constraint (offset.data(), step.data(), count.data());
58
59 // Just declare offset,count and step in the int32 type.
60 vector<int32>offset32;
61 offset32.resize(rank);
62 vector<int32>count32;
63 count32.resize(rank);
64 vector<int32>step32;
65 step32.resize(rank);
66
67
68 // Just obtain the offset,count and step in the int32 type.
69 for (int i = 0; i < rank; i++) {
70 offset32[i] = (int32) offset[i];
71 count32[i] = (int32) count[i];
72 step32[i] = (int32) step[i];
73 }
74
75 // Loop through the functions to obtain lat/lon for the specific HDF4 products the handler
76 // supports.
77 switch (sptype) {
78
79 // TRMM swath
80 case TRMML2_V6:
81 {
82 readtrmml2_v6 (offset32.data(), count32.data(), step32.data(), nelms);
83 break;
84 }
85
86 // TRMM grid
87 case TRMML3A_V6:
88 {
89 readtrmml3a_v6 (offset32.data(), count32.data(), step32.data(), nelms);
90 break;
91 }
92
93 case TRMML3B_V6:
94 {
95 readtrmml3b_v6 (offset32.data(), count32.data(), step32.data(), nelms);
96 break;
97 }
98
99 case TRMML3C_V6:
100 {
101 readtrmml3c_v6 (offset32.data(), count32.data(), step32.data(), nelms);
102 break;
103 }
104
105
106 // TRMM V7 grid
107 case TRMML3S_V7:
108 case TRMML3M_V7:
109 {
110 readtrmml3_v7 (offset32.data(), step32.data(), nelms);
111 break;
112 }
113 // CERES CER_AVG_Aqua-FM3-MODIS,CER_AVG_Terra-FM1-MODIS products
114 case CER_AVG:
115 {
116 readceravgsyn (offset32.data(), count32.data(), step32.data(), nelms);
117 break;
118 }
119
120 // CERES CER_ES4_??? products
121 case CER_ES4:
122 {
123 readceres4ig (offset32.data(), count32.data(), step32.data(), nelms);
124 break;
125 }
126
127 // CERES CER_ISCCP-D2like-Day product
128 case CER_CDAY:
129 {
130 readcersavgid2 (offset.data(), count.data(), step.data(), nelms);
131 break;
132 }
133
134 // CERES CER_ISCCP-D2like-GEO product
135 case CER_CGEO:
136 {
137 readceres4ig (offset32.data(), count32.data(), step32.data(), nelms);
138 break;
139 }
140
141 // CERES CER_SRBAVG3_Aqua product
142 case CER_SRB:
143 {
144 // When longitude is fixed
145 if (rank == 1) {
146 readcersavgid1 (offset.data(), count.data(), step.data(), nelms);
147 }
148 // When longitude is not fixed
149 else if (rank == 2) {
150 readcersavgid2 (offset.data(), count.data(), step.data(), nelms);
151 }
152 break;
153 }
154
155 // CERES SYN Aqua products
156 case CER_SYN:
157 {
158 readceravgsyn (offset32.data(), count32.data(), step32.data(), nelms);
159 break;
160 }
161
162 // CERES Zonal Average products
163 case CER_ZAVG:
164 {
165 readcerzavg (offset32.data(), count32.data(), step32.data(), nelms);
166 break;
167 }
168
169 // OBPG level 2 products
170 case OBPGL2:
171 {
172 readobpgl2 (offset32.data(), count32.data(), step32.data(), nelms);
173 break;
174 }
175
176 // OBPG level 3 products
177 case OBPGL3:
178 {
179 readobpgl3 (offset.data(), step.data(), nelms);
180 break;
181 }
182
183 // We don't handle any OtherHDF products
184 case OTHERHDF:
185 {
186 throw InternalErr (__FILE__, __LINE__, "Unsupported HDF files");
187
188 }
189 default:
190 {
191
192 throw InternalErr (__FILE__, __LINE__, "Unsupported HDF files");
193
194 }
195 }
196
197 return true;
198
199}
200
201
202// Read TRMM level 2 lat/lon. We need to split geolocation field.
203// geolocation[YDIM][XDIM][0] is latitude
204// geolocation[YDIM][XDIM][1] is longitude
205
206void
207HDFSPArrayGeoField::readtrmml2_v6 (int32 * offset32, int32 * count32,
208 int32 * step32, int nelms)
209{
210
211#if 0
212 string check_pass_fileid_key_str="H4.EnablePassFileID";
213 bool check_pass_fileid_key = false;
214 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
215#endif
216
217 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
218
219 int32 sdid = -1;
220
221 if(false == check_pass_fileid_key) {
222 sdid = SDstart (const_cast < char *>(filename.c_str ()), DFACC_READ);
223 if (sdid < 0) {
224 ostringstream eherr;
225 eherr << "File " << filename.c_str () << " cannot be open.";
226 throw InternalErr (__FILE__, __LINE__, eherr.str ());
227 }
228 }
229 else
230 sdid = sdfd;
231
232 int32 sdsid = 0;
233
234 vector<int32>geooffset32;
235 geooffset32.resize(rank+1);
236
237 vector<int32>geocount32;
238 geocount32.resize(rank+1);
239
240 vector<int32>geostep32;
241 geostep32.resize(rank+1);
242
243 for (int i = 0; i < rank; i++) {
244 geooffset32[i] = offset32[i];
245 geocount32[i] = count32[i];
246 geostep32[i] = step32[i];
247 }
248
249 if (fieldtype == 1) {
250 geooffset32[rank] = 0;
251 geocount32[rank] = 1;
252 geostep32[rank] = 1;
253 }
254
255 if (fieldtype == 2) {
256 geooffset32[rank] = 1;
257 geocount32[rank] = 1;
258 geostep32[rank] = 1;
259 }
260
261 int32 sdsindex = SDreftoindex (sdid, (int32) fieldref);
262
263 if (sdsindex == -1) {
264 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
265 ostringstream eherr;
266 eherr << "SDS index " << sdsindex << " is not right.";
267 throw InternalErr (__FILE__, __LINE__, eherr.str ());
268 }
269
270 sdsid = SDselect (sdid, sdsindex);
271 if (sdsid < 0) {
272 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
273 ostringstream eherr;
274 eherr << "SDselect failed.";
275 throw InternalErr (__FILE__, __LINE__, eherr.str ());
276 }
277
278 int32 r = 0;
279
280 switch (dtype) {
281
282 case DFNT_INT8:
283 {
284 vector <int8> val;
285 val.resize(nelms);
286 r = SDreaddata (sdsid, geooffset32.data(), geostep32.data(), geocount32.data(), (void*)(val.data()));
287 if (r != 0) {
288 SDendaccess (sdsid);
289 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
290 ostringstream eherr;
291 eherr << "SDreaddata failed.";
292 throw InternalErr (__FILE__, __LINE__, eherr.str ());
293 }
294
295#ifndef SIGNED_BYTE_TO_INT32
296 set_value ((dods_byte *) val.data(), nelms);
297#else
298 vector<int32>newval;
299 newval.resize(nelms);
300
301 for (int counter = 0; counter < nelms; counter++)
302 newval[counter] = (int32) (val[counter]);
303 set_value ((dods_int32 *) newval.data(), nelms);
304#endif
305 }
306
307 break;
308 case DFNT_UINT8:
309 case DFNT_UCHAR8:
310 {
311 vector<uint8> val;
312 val.resize(nelms);
313 r = SDreaddata (sdsid, geooffset32.data(), geostep32.data(), geocount32.data(), (void*)(val.data()));
314 if (r != 0) {
315 SDendaccess (sdsid);
316 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
317 ostringstream eherr;
318 eherr << "SDreaddata failed";
319 throw InternalErr (__FILE__, __LINE__, eherr.str ());
320 }
321
322 set_value ((dods_byte *) val.data(), nelms);
323 }
324 break;
325
326 case DFNT_INT16:
327 {
328 vector<int16> val;
329 val.resize(nelms);
330 r = SDreaddata (sdsid, geooffset32.data(), geostep32.data(), geocount32.data(), (void*)(val.data()));
331 if (r != 0) {
332 SDendaccess (sdsid);
333 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
334 ostringstream eherr;
335 eherr << "SDreaddata failed";
336 throw InternalErr (__FILE__, __LINE__, eherr.str ());
337 }
338
339 set_value ((dods_int16 *) val.data(), nelms);
340 }
341 break;
342
343 case DFNT_UINT16:
344 {
345 vector<uint16>val;
346 val.resize(nelms);
347 r = SDreaddata (sdsid, geooffset32.data(), geostep32.data(), geocount32.data(), (void*)(val.data()));
348 if (r != 0) {
349 SDendaccess (sdsid);
350 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
351 ostringstream eherr;
352 eherr << "SDreaddata failed";
353 throw InternalErr (__FILE__, __LINE__, eherr.str ());
354 }
355
356 set_value ((dods_uint16 *) val.data(), nelms);
357 }
358 break;
359 case DFNT_INT32:
360 {
361 vector<int32>val;
362 val.resize(nelms);
363 r = SDreaddata (sdsid, geooffset32.data(), geostep32.data(), geocount32.data(), (void*)(val.data()));
364 if (r != 0) {
365 SDendaccess (sdsid);
366 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
367 ostringstream eherr;
368 eherr << "SDreaddata failed";
369 throw InternalErr (__FILE__, __LINE__, eherr.str ());
370 }
371
372 set_value ((dods_int32 *) val.data(), nelms);
373 }
374 break;
375 case DFNT_UINT32:
376 {
377 vector<uint32>val;
378 val.resize(nelms);
379 r = SDreaddata (sdsid, geooffset32.data(), geostep32.data(), geocount32.data(), (void*)(val.data()));
380 if (r != 0) {
381 SDendaccess (sdsid);
382 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
383 ostringstream eherr;
384 eherr << "SDreaddata failed";
385 throw InternalErr (__FILE__, __LINE__, eherr.str ());
386 }
387 set_value ((dods_uint32 *) val.data(), nelms);
388 }
389 break;
390 case DFNT_FLOAT32:
391 {
392 vector<float32>val;
393 val.resize(nelms);
394 r = SDreaddata (sdsid, geooffset32.data(), geostep32.data(), geocount32.data(), (void*)(val.data()));
395 if (r != 0) {
396 SDendaccess (sdsid);
397 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
398 ostringstream eherr;
399 eherr << "SDreaddata failed";
400 throw InternalErr (__FILE__, __LINE__, eherr.str ());
401 }
402
403 set_value ((dods_float32 *) val.data(), nelms);
404 }
405 break;
406 case DFNT_FLOAT64:
407 {
408 vector<float64>val;
409 val.resize(nelms);
410 r = SDreaddata (sdsid, geooffset32.data(), geostep32.data(), geocount32.data(), (void*)(val.data()));
411 if (r != 0) {
412 SDendaccess (sdsid);
413 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
414 ostringstream eherr;
415 eherr << "SDreaddata failed";
416 throw InternalErr (__FILE__, __LINE__, eherr.str ());
417 }
418
419 set_value ((dods_float64 *) val.data(), nelms);
420 }
421 break;
422 default:
423 SDendaccess (sdsid);
424 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
425 throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
426 }
427
428 r = SDendaccess (sdsid);
429 if (r != 0) {
430 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
431 ostringstream eherr;
432 eherr << "SDendaccess failed.";
433 throw InternalErr (__FILE__, __LINE__, eherr.str ());
434 }
435
436 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
437
438}
439
440// Read TRMM level 3a version 6 lat/lon.
441// We follow TRMM document to retrieve the lat and lon.
442void
443HDFSPArrayGeoField::readtrmml3a_v6 (int32 * offset32, int32 * count32,
444 int32 * step32, int nelms)
445{
446
447 const float slat = 89.5;
448 const float slon = 0.5;
449 vector<float> val;
450 val.resize(nelms);
451
452 if (fieldtype == 1) {//latitude
453
454 int icount = 0;
455 float sval = slat - offset32[0];
456
457 while (icount < (int) (count32[0])) {
458 val[icount] = sval - step32[0] * icount;
459 icount++;
460 }
461 }
462
463 if (fieldtype == 2) { //longitude
464 int icount = 0;
465 float sval = slon + offset32[0];
466
467 while (icount < (int) (count32[0])) {
468 val[icount] = sval + step32[0] * icount;
469 icount++;
470 }
471 }
472 set_value ((dods_float32 *) val.data(), nelms);
473}
474
475
476// TRMM level 3 case. Have to follow http://disc.sci.gsfc.nasa.gov/additional/faq/precipitation_faq.shtml#lat_lon
477// to calculate lat/lon.
478void
479HDFSPArrayGeoField::readtrmml3b_v6 (int32 * offset32, int32 * count32,
480 int32 * step32, int nelms)
481{
482
483 const float slat = -49.875;
484 const float slon = -179.875;
485 vector<float> val;
486 val.resize(nelms);
487
488 if (fieldtype == 1) {//latitude
489
490 int icount = 0;
491 float sval = slat + 0.25 * (int) (offset32[0]);
492
493 while (icount < (int) (count32[0])) {
494 val[icount] = sval + 0.25 * (int) (step32[0]) * icount;
495 icount++;
496 }
497 }
498
499 if (fieldtype == 2) { //longitude
500 int icount = 0;
501 float sval = slon + 0.25 * (int) (offset32[0]);
502
503 while (icount < (int) (count32[0])) {
504 val[icount] = sval + 0.25 * (int) (step32[0]) * icount;
505 icount++;
506 }
507 }
508 set_value ((dods_float32 *) val.data(), nelms);
509}
510
511// Read TRMM level 3c(CSH) version 6 lat/lon.
512// We follow TRMM document to retrieve the lat and lon.
513void
514HDFSPArrayGeoField::readtrmml3c_v6 (int32 * offset32, int32 * count32,
515 int32 * step32, int nelms)
516{
517
518 const float slat = -36.75;
519 const float slon = -179.75;
520 vector<float> val;
521 val.resize(nelms);
522
523 if (fieldtype == 1) {//latitude
524
525 int icount = 0;
526 float sval = slat + 0.5 * (int) (offset32[0]);
527
528 while (icount < (int) (count32[0])) {
529 val[icount] = sval + 0.5 * (int) (step32[0]) * icount;
530 icount++;
531 }
532 }
533
534 if (fieldtype == 2) { //longitude
535 int icount = 0;
536 float sval = slon + 0.5 * (int) (offset32[0]);
537
538 while (icount < (int) (count32[0])) {
539 val[icount] = sval + 0.5 * (int) (step32[0]) * icount;
540 icount++;
541 }
542 }
543 set_value ((dods_float32 *) val.data(), nelms);
544}
545// Read latlon of TRMM version 7 products.
546// Lat/lon parameters can be retrieved from attribute GridHeader.
547void
548HDFSPArrayGeoField::readtrmml3_v7 (int32 * offset32,
549 int32 * step32, int nelms)
550{
551
552#if 0
553 string check_pass_fileid_key_str="H4.EnablePassFileID";
554 bool check_pass_fileid_key = false;
555 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
556#endif
557
558 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
559
560 int32 sdid = -1;
561
562 if(false == check_pass_fileid_key) {
563 sdid = SDstart (const_cast < char *>(filename.c_str ()), DFACC_READ);
564 if (sdid < 0) {
565 ostringstream eherr;
566 eherr << "File " << filename.c_str () << " cannot be open.";
567 throw InternalErr (__FILE__, __LINE__, eherr.str ());
568 }
569 }
570 else
571 sdid = sdfd;
572
573
574 string gridinfo_name = "GridHeader";
575 intn status = 0;
576
577 if(fieldref != -1) {
578
579 if (fieldref >9) {
580 throw InternalErr (__FILE__,__LINE__,
581 "The maximum number of grids to be supported in the current implementation is 9.");
582 }
583
584 else {
585 ostringstream fieldref_str;
586 fieldref_str << fieldref;
587 gridinfo_name = gridinfo_name + fieldref_str.str();
588 }
589 }
590
591 int32 attr_index = 0;
592 attr_index = SDfindattr (sdid, gridinfo_name.c_str());
593 if (attr_index == FAIL) {
594 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
595 string err_mesg = "SDfindattr failed,should find attribute "+gridinfo_name+" .";
596 throw InternalErr (__FILE__, __LINE__, err_mesg);
597 }
598
599 int32 attr_dtype = 0;
600 int32 n_values = 0;
601
602 char attr_name[H4_MAX_NC_NAME];
603 status =
604 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
605 if (status == FAIL) {
606 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
607 throw InternalErr (__FILE__, __LINE__, "SDattrinfo failed ");
608 }
609
610 vector<char> attr_value;
611 attr_value.resize(n_values * DFKNTsize(attr_dtype));
612
613 status = SDreadattr (sdid, attr_index, attr_value.data());
614 if (status == FAIL) {
615 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
616 throw InternalErr (__FILE__, __LINE__, "SDreadattr failed ");
617 }
618
619 float lat_start = 0;;
620 float lon_start = 0.;
621 float lat_res = 0.;
622 float lon_res = 0.;
623
624 int latsize = 0;
625 int lonsize = 0;
626
627 HDFCFUtil::parser_trmm_v7_gridheader(attr_value,latsize,lonsize,
628 lat_start,lon_start,lat_res,lon_res,false);
629//cerr<<"lat_start is "<<lat_start <<endl;
630//cerr<<"lon_start is "<<lon_start<<endl;
631//cerr <<"lat_res is "<<lat_res <<endl;
632//cerr<<"lon_res is "<<lon_res <<endl;
633
634 if(0 == latsize || 0 == lonsize) {
635 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
636 throw InternalErr (__FILE__, __LINE__, "Either latitude or longitude size is 0. ");
637 }
638
639
640 vector<float>val;
641 val.resize(nelms);
642
643 if(fieldtype == 1) {
644 for (int i = 0; i < nelms; ++i)
645 val[i] = lat_start+offset32[0]*lat_res+lat_res/2 + i*lat_res*step32[0];
646 }
647 else if(fieldtype == 2) {
648 for (int i = 0; i < nelms; ++i)
649 val[i] = lon_start+offset32[0]*lon_res+lon_res/2 + i*lon_res*step32[0];
650 }
651
652 set_value ((dods_float32 *) val.data(), nelms);
653
654 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
655}
656
657
658
659// OBPG Level 2 lat/lon including CZCS, MODISA, MODIST, OCTS and SewWIFS.
660// We need to use information retrieved from the attribute to interpoloate
661// the latitude/longitude. This is similar to the Swath dimension map case.
662// "Pixels per Scan Line" and "Number of Pixel Control Points"
663// should be used to interpolate.
664// "Pixels per Scan Line" is the final number of elements for lat/lon along the 2nd dimension
665// "Number of Pixel Control Points" includes the original number of elements for lat/lon along
666// the 2nd dimension.
667void
668HDFSPArrayGeoField::readobpgl2 (int32 * offset32, int32 * count32,
669 int32 * step32, int nelms)
670{
671#if 0
672 string check_pass_fileid_key_str="H4.EnablePassFileID";
673 bool check_pass_fileid_key = false;
674 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
675#endif
676
677 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
678
679 int32 sdid = -1;
680
681 if(false == check_pass_fileid_key) {
682 sdid = SDstart (const_cast < char *>(filename.c_str ()), DFACC_READ);
683 if (sdid < 0) {
684 ostringstream eherr;
685 eherr << "File " << filename.c_str () << " cannot be open.";
686 throw InternalErr (__FILE__, __LINE__, eherr.str ());
687 }
688 }
689 else
690 sdid = sdfd;
691
692 int32 sdsid = -1;
693 intn status = 0;
694
695 // Read File attributes to otain the segment
696 int32 attr_index = 0;
697 int32 attr_dtype = 0;
698 int32 n_values = 0;
699 int32 num_pixel_data = 0;
700 int32 num_point_data = 0;
701 int32 num_scan_data = 0;
702
703 attr_index = SDfindattr (sdid, NUM_PIXEL_NAME);
704 if (attr_index == FAIL) {
705 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
706 string attr_name(NUM_PIXEL_NAME);
707 string err_mesg = "SDfindattr failed,should find attribute "+attr_name+" .";
708 throw InternalErr (__FILE__, __LINE__, err_mesg);
709 }
710
711 char attr_name[H4_MAX_NC_NAME];
712 status =
713 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
714 if (status == FAIL) {
715 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
716 throw InternalErr (__FILE__, __LINE__, "SDattrinfo failed ");
717 }
718
719 if (n_values != 1) {
720 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
721 throw InternalErr (__FILE__, __LINE__,
722 "Only one value of number of scan line ");
723 }
724
725 status = SDreadattr (sdid, attr_index, &num_pixel_data);
726 if (status == FAIL) {
727 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
728 throw InternalErr (__FILE__, __LINE__, "SDreadattr failed ");
729 }
730
731 attr_index = SDfindattr (sdid, NUM_POINTS_LINE_NAME);
732 if (attr_index == FAIL) {
733 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
734 string attr_name(NUM_POINTS_LINE_NAME);
735 string err_mesg = "SDfindattr failed,should find attribute "+attr_name+" .";
736 throw InternalErr (__FILE__, __LINE__, err_mesg);
737
738 }
739
740 status =
741 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
742 if (status == FAIL) {
743 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
744 throw InternalErr (__FILE__, __LINE__, "SDattrinfo failed ");
745 }
746
747 if (n_values != 1) {
748 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
749 throw InternalErr (__FILE__, __LINE__,
750 "Only one value of number of point ");
751 }
752
753 status = SDreadattr (sdid, attr_index, &num_point_data);
754 if (status == FAIL) {
755 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
756 throw InternalErr (__FILE__, __LINE__, "SDreadattr failed ");
757 }
758
759 attr_index = SDfindattr (sdid, NUM_SCAN_LINE_NAME);
760 if (attr_index == FAIL) {
761 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
762 string attr_name(NUM_SCAN_LINE_NAME);
763 string err_mesg = "SDfindattr failed,should find attribute "+attr_name+" .";
764 throw InternalErr (__FILE__, __LINE__, err_mesg);
765
766 }
767
768 status =
769 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
770 if (status == FAIL) {
771 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
772 throw InternalErr (__FILE__, __LINE__, "SDattrinfo failed ");
773 }
774
775 if (n_values != 1) {
776 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
777 throw InternalErr (__FILE__, __LINE__,"Only one value of number of point ");
778 }
779
780 status = SDreadattr (sdid, attr_index, &num_scan_data);
781 if (status == FAIL) {
782 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
783 throw InternalErr (__FILE__, __LINE__, "SDreadattr failed ");
784 }
785
786
787 if ( 0 == num_scan_data || 0 == num_point_data || 0 == num_pixel_data) {
788 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
789 throw InternalErr (__FILE__, __LINE__, "num_scan or num_point or num_pixel should not be zero. ");
790 }
791
792 if ( 1 == num_point_data && num_pixel_data != 1) {
793 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
794 throw InternalErr (__FILE__, __LINE__,
795 "num_point is 1 and num_pixel is not 1, interpolation cannot be done ");
796 }
797
798 bool compmapflag = false;
799 if (num_pixel_data == num_point_data)
800 compmapflag = true;
801
802 int32 sdsindex = SDreftoindex (sdid, (int32) fieldref);
803 if (sdsindex == -1) {
804 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
805 ostringstream eherr;
806 eherr << "SDS index " << sdsindex << " is not right.";
807 throw InternalErr (__FILE__, __LINE__, eherr.str ());
808 }
809
810 sdsid = SDselect (sdid, sdsindex);
811 if (sdsid < 0) {
812 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
813 ostringstream eherr;
814 eherr << "SDselect failed.";
815 throw InternalErr (__FILE__, __LINE__, eherr.str ());
816 }
817
818 int32 r = 0;
819
820 switch (dtype) {
821 case DFNT_INT8:
822 case DFNT_UINT8:
823 case DFNT_UCHAR8:
824 case DFNT_CHAR8:
825 case DFNT_INT16:
826 case DFNT_UINT16:
827 case DFNT_INT32:
828 case DFNT_UINT32:
829 case DFNT_FLOAT64:
830 SDendaccess (sdsid);
831 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
832 throw InternalErr (__FILE__, __LINE__,"datatype is not float, unsupported.");
833 case DFNT_FLOAT32:
834 {
835 vector<float32> val;
836 val.resize(nelms);
837 if (compmapflag) {
838 r = SDreaddata (sdsid, offset32, step32, count32, val.data());
839 if (r != 0) {
840 SDendaccess (sdsid);
841 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
842 ostringstream eherr;
843 eherr << "SDreaddata failed";
844 throw InternalErr (__FILE__, __LINE__, eherr.str ());
845 }
846 }
847 else {
848
849 int total_elm = num_scan_data * num_point_data;
850 vector<float32>orival;
851 orival.resize(total_elm);
852 int32 all_start[2],all_edges[2];
853
854 all_start[0] = 0;
855 all_start[1] = 0;
856 all_edges[0] = num_scan_data;
857 all_edges[1] = num_point_data;
858
859 r = SDreaddata (sdsid, all_start, NULL, all_edges,
860 orival.data());
861 if (r != 0) {
862 SDendaccess (sdsid);
863 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
864 ostringstream eherr;
865 eherr << "SDreaddata failed";
866 throw InternalErr (__FILE__, __LINE__, eherr.str ());
867 }
868 int interpolate_elm = num_scan_data *num_pixel_data;
869
870 vector<float32> interp_val;
871 interp_val.resize(interpolate_elm);
872
873 // Number of scan line doesn't change, so just interpolate according to the fast changing dimension
874 int tempseg = 0;
875 float tempdiff = 0;
876
877 if (num_pixel_data % num_point_data == 0)
878 tempseg = num_pixel_data / num_point_data;
879 else
880 tempseg = num_pixel_data / num_point_data + 1;
881
882 int last_tempseg =
883 (num_pixel_data%num_point_data)?(num_pixel_data-1-(tempseg*(num_point_data-2))):tempseg;
884
885 if ( 0 == last_tempseg || 0 == tempseg) {
886 SDendaccess(sdsid);
887 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
888 throw InternalErr(__FILE__,__LINE__,"Segments cannot be zero");
889 }
890
891 int interp_val_index = 0;
892
893 for (int i = 0; i <num_scan_data; i++) {
894
895 // All the segements except the last one
896 for( int j =0; j <num_point_data -2; j ++) {
897 tempdiff = orival[i*num_point_data+j+1] - orival[i*num_point_data+j];
898 for (int k = 0; k <tempseg; k++) {
899 interp_val[interp_val_index] = orival[i*num_point_data+j] +
900 tempdiff/tempseg *k;
901 interp_val_index++;
902 }
903 }
904
905 // The last segment
906 tempdiff = orival[i*num_point_data+num_point_data-1]
907 -orival[i*num_point_data+num_point_data-2];
908 for (int k = 0; k <last_tempseg; k++) {
909 interp_val[interp_val_index] = orival[i*num_point_data+num_point_data-2] +
910 tempdiff/last_tempseg *k;
911 interp_val_index++;
912 }
913
914 interp_val[interp_val_index] = orival[i*num_point_data+num_point_data-1];
915 interp_val_index++;
916
917 }
918
919 LatLon2DSubset(val.data(),num_scan_data,num_pixel_data,interp_val.data(),
920 offset32,count32,step32);
921
922 }
923 // Leave the following comments
924#if 0
925 // WE SHOULD INTERPOLATE ACCORDING TO THE FAST CHANGING DIMENSION
926 // THis method will save some memory, but it will cause greater error
927 // if the step is very large. However, we should still take advantage
928 // of the small memory approach in the future implementation. KY 2012-09-04
929 float tempdiff;
930 int i, j, k, k2;
931 int32 realcount2 = oricount32[1] - 1;
932
933 for (i = 0; i < (int) count32[0]; i++) {
934 for (j = 0; j < (int) realcount2 - 1; j++) {
935 tempdiff =
936 orival[i * (int) oricount32[1] + j + 1] -
937 orival[i * (int) oricount32[1] + j];
938 for (k = 0; k < tempnewseg; k++) {
939 val[i * (int) count32[1] + j * tempnewseg + k] =
940 orival[i * (int) oricount32[1] + j] +
941 tempdiff / tempnewseg * k;
942 }
943 }
944 tempdiff =
945 orival[i * (int) oricount32[1] + j + 1] -
946 orival[i * (int) oricount32[1] + j];
947 // There are three cases:
948 // 1. when the oricount is 1
949 // int lastseg = num_pixel_data - tempnewseg*((int)oricount32[1]-1)-1;
950 int lastseg =
951 (int) (count32[1] -
952 tempnewseg * (int) (realcount2 - 1));
953 for (k2 = 0; k2 <= lastseg; k2++)
954 val[i * (int) count32[1] + j * tempnewseg + k2] =
955 orival[i * (int) oricount32[1] + j] +
956 tempdiff / lastseg * k2;
957 }
958
959 delete[](float32 *) orival;
960 }
961#endif
962
963 set_value ((dods_float32 *) val.data(), nelms);
964 }
965 break;
966 default:
967 SDendaccess (sdsid);
968 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
969 throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
970 }
971
972 r = SDendaccess (sdsid);
973 if (r != 0) {
974 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
975 ostringstream eherr;
976 eherr << "SDendaccess failed.";
977 throw InternalErr (__FILE__, __LINE__, eherr.str ());
978 }
979
980 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
981}
982
983// Obtain lat/lon for OBPG CZCS, MODISA, MODIST,OCTS and SeaWIFS products
984// lat/lon should be calculated based on the file attribute.
985// Geographic projection: lat,lon starting point and also lat/lon steps.
986void
987HDFSPArrayGeoField::readobpgl3 (int *offset, int *step, int nelms)
988{
989
990#if 0
991 string check_pass_fileid_key_str="H4.EnablePassFileID";
992 bool check_pass_fileid_key = false;
993 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
994#endif
995
996 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
997
998 int32 sdid = -1;
999
1000 if(false == check_pass_fileid_key) {
1001 sdid = SDstart (const_cast < char *>(filename.c_str ()), DFACC_READ);
1002 if (sdid < 0) {
1003 ostringstream eherr;
1004 eherr << "File " << filename.c_str () << " cannot be open.";
1005 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1006 }
1007 }
1008 else
1009 sdid = sdfd;
1010
1011 intn status = 0;
1012
1013 // Read File attributes to otain the segment
1014 int32 attr_index = 0;
1015 int32 attr_dtype = 0;
1016 int32 n_values = 0;
1017 int32 num_lat_data = 0;
1018 int32 num_lon_data = 0;
1019 float32 lat_step = 0.;
1020 float32 lon_step = 0.;
1021 float32 swp_lat = 0.;
1022 float32 swp_lon = 0.;
1023
1024 // Obtain number of latitude
1025 attr_index = SDfindattr (sdid, NUM_LAT_NAME);
1026 if (attr_index == FAIL) {
1027 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1028 string attr_name(NUM_LAT_NAME);
1029 string err_mesg = "SDfindattr failed,should find attribute "+attr_name+" .";
1030 throw InternalErr (__FILE__, __LINE__, err_mesg);
1031 //throw InternalErr (__FILE__, __LINE__, "SDfindattr failed,should find this attribute. ");
1032 }
1033
1034 char attr_name[H4_MAX_NC_NAME];
1035 status =
1036 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
1037 if (status == FAIL) {
1038 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1039 throw InternalErr (__FILE__, __LINE__, "SDattrinfo failed ");
1040 }
1041
1042 if (n_values != 1) {
1043 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1044 throw InternalErr (__FILE__, __LINE__, "Only should have one value ");
1045 }
1046
1047 status = SDreadattr (sdid, attr_index, &num_lat_data);
1048 if (status == FAIL) {
1049 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1050 throw InternalErr (__FILE__, __LINE__, "SDreadattr failed ");
1051 }
1052
1053 // Obtain number of longitude
1054 attr_index = SDfindattr (sdid, NUM_LON_NAME);
1055 if (attr_index == FAIL) {
1056 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1057 string attr_name2(NUM_LON_NAME);
1058 string err_mesg = "SDfindattr failed,should find attribute "+attr_name2+" .";
1059 throw InternalErr (__FILE__, __LINE__, err_mesg);
1060// throw InternalErr (__FILE__, __LINE__, "SDfindattr failed ");
1061 }
1062
1063 status =
1064 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
1065 if (status == FAIL) {
1066 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1067 throw InternalErr (__FILE__, __LINE__, "SDattrinfo failed ");
1068 }
1069
1070 if (n_values != 1) {
1071 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1072 throw InternalErr (__FILE__, __LINE__, "Only should have one value ");
1073 }
1074
1075 status = SDreadattr (sdid, attr_index, &num_lon_data);
1076 if (status == FAIL) {
1077 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1078 throw InternalErr (__FILE__, __LINE__, "SDreadattr failed ");
1079 }
1080
1081 // obtain latitude step
1082 attr_index = SDfindattr (sdid, LAT_STEP_NAME);
1083 if (attr_index == FAIL) {
1084 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1085 string attr_name2(LAT_STEP_NAME);
1086 string err_mesg = "SDfindattr failed,should find attribute "+attr_name2+" .";
1087 throw InternalErr (__FILE__, __LINE__, err_mesg);
1088// throw InternalErr (__FILE__, __LINE__, "SDfindattr failed ");
1089 }
1090
1091 status =
1092 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
1093 if (status == FAIL) {
1094 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1095 throw InternalErr (__FILE__, __LINE__, "SDattrinfo failed ");
1096 }
1097
1098 if (n_values != 1) {
1099 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1100 throw InternalErr (__FILE__, __LINE__, "Only should have one value ");
1101 }
1102
1103 status = SDreadattr (sdid, attr_index, &lat_step);
1104 if (status == FAIL) {
1105 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1106 throw InternalErr (__FILE__, __LINE__, "SDreadattr failed ");
1107 }
1108
1109 // Obtain longitude step
1110 attr_index = SDfindattr (sdid, LON_STEP_NAME);
1111 if (attr_index == FAIL) {
1112 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1113 string attr_name2(LON_STEP_NAME);
1114 string err_mesg = "SDfindattr failed,should find attribute "+attr_name2+" .";
1115 throw InternalErr (__FILE__, __LINE__, err_mesg);
1116// throw InternalErr (__FILE__, __LINE__, "SDfindattr failed ");
1117 }
1118
1119 status =
1120 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
1121 if (status == FAIL) {
1122 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1123 throw InternalErr (__FILE__, __LINE__, "SDattrinfo failed ");
1124 }
1125
1126 if (n_values != 1) {
1127 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1128 throw InternalErr (__FILE__, __LINE__, "Only should have one value ");
1129 }
1130
1131 status = SDreadattr (sdid, attr_index, &lon_step);
1132 if (status == FAIL) {
1133 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1134 throw InternalErr (__FILE__, __LINE__, "SDreadattr failed ");
1135 }
1136
1137 // obtain south west corner latitude
1138 attr_index = SDfindattr (sdid, SWP_LAT_NAME);
1139 if (attr_index == FAIL) {
1140 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1141 string attr_name2(SWP_LAT_NAME);
1142 string err_mesg = "SDfindattr failed,should find attribute "+attr_name2+" .";
1143 throw InternalErr (__FILE__, __LINE__, err_mesg);
1144//throw InternalErr (__FILE__, __LINE__, "SDfindattr failed ");
1145 }
1146
1147 status =
1148 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
1149 if (status == FAIL) {
1150 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1151 throw InternalErr (__FILE__, __LINE__, "SDattrinfo failed ");
1152 }
1153
1154 if (n_values != 1) {
1155 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1156 throw InternalErr (__FILE__, __LINE__, "Only should have one value ");
1157 }
1158
1159 status = SDreadattr (sdid, attr_index, &swp_lat);
1160 if (status == FAIL) {
1161 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1162 throw InternalErr (__FILE__, __LINE__, "SDreadattr failed ");
1163 }
1164
1165 // obtain south west corner longitude
1166 attr_index = SDfindattr (sdid, SWP_LON_NAME);
1167 if (attr_index == FAIL) {
1168 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1169 string attr_name2(SWP_LON_NAME);
1170 string err_mesg = "SDfindattr failed,should find attribute "+attr_name2+" .";
1171 throw InternalErr (__FILE__, __LINE__, err_mesg);
1172//throw InternalErr (__FILE__, __LINE__, "SDfindattr failed,should find this attribute. ");
1173 }
1174
1175 status =
1176 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
1177 if (status == FAIL) {
1178 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1179 throw InternalErr (__FILE__, __LINE__, "SDattrinfo failed ");
1180 }
1181
1182 if (n_values != 1) {
1183 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1184 throw InternalErr (__FILE__, __LINE__, "Only should have one value ");
1185 }
1186
1187 status = SDreadattr (sdid, attr_index, &swp_lon);
1188 if (status == FAIL) {
1189 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1190 throw InternalErr (__FILE__, __LINE__, "SDreadattr failed ");
1191 }
1192
1193
1194 if (fieldtype == 1) {
1195
1196 vector<float32> allval;
1197 allval.resize(num_lat_data);
1198
1199 // The first index is from the north, so we need to reverse the calculation of the index
1200
1201 for (int j = 0; j < num_lat_data; j++)
1202 allval[j] = (num_lat_data - j - 1) * lat_step + swp_lat;
1203
1204 vector<float32>val;
1205 val.resize(nelms);
1206
1207 for (int k = 0; k < nelms; k++)
1208 val[k] = allval[(int) (offset[0] + step[0] * k)];
1209
1210 set_value ((dods_float32 *) val.data(), nelms);
1211 }
1212
1213 if (fieldtype == 2) {
1214
1215 vector<float32>allval;
1216 allval.resize(num_lon_data);
1217 for (int j = 0; j < num_lon_data; j++)
1218 allval[j] = swp_lon + j * lon_step;
1219
1220 vector<float32>val;
1221 val.resize(nelms);
1222 for (int k = 0; k < nelms; k++)
1223 val[k] = allval[(int) (offset[0] + step[0] * k)];
1224
1225 set_value ((dods_float32 *) val.data(), nelms);
1226
1227 }
1228
1229 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1230
1231}
1232
1233
1234// CERES SAVG and ISSP DAY case(CER_SAVG and CER_ISCCP_.._DAY )
1235// We need to calculate latitude/longitude based on CERES nested grid formula
1236// http://eosweb.larc.nasa.gov/PRODOCS/ceres/SRBAVG/Quality_Summaries/srbavg_ed2d/nestedgrid.html
1237// or https://eosweb.larc.nasa.gov/sites/default/files/project/ceres/quality_summaries/srbavg_ed2d/nestedgrid.pdf
1238void
1239HDFSPArrayGeoField::readcersavgid2 (int *offset, int *count, int *step,
1240 int nelms)
1241{
1242
1243 const int dimsize0 = 180;
1244 const int dimsize1 = 360;
1245 float32 val[count[0]][count[1]];
1246 float32 orival[dimsize0][dimsize1];
1247
1248 // Following CERES Nested grid
1249 // URL http://eosweb.larc.nasa.gov/PRODOCS/ceres/SRBAVG/Quality_Summaries/srbavg_ed2d/nestedgrid.html
1250 // or https://eosweb.larc.nasa.gov/sites/default/files/project/ceres/quality_summaries/srbavg_ed2d/nestedgrid.pdf
1251 if (fieldtype == 1) {// Calculate the latitude
1252 for (int i = 0; i < dimsize0; i++)
1253 for (int j = 0; j < dimsize1; j++)
1254 orival[i][j] = 89.5 - i;
1255
1256 for (int i = 0; i < count[0]; i++) {
1257 for (int j = 0; j < count[1]; j++) {
1258 val[i][j] =
1259 orival[offset[0] + step[0] * i][offset[1] + step[1] * j];
1260 }
1261 }
1262 }
1263
1264 if (fieldtype == 2) {// Calculate the longitude
1265
1266 int i = 0;
1267 int j = 0;
1268 int k = 0;
1269
1270 // Longitude extent
1271 int lonextent;
1272
1273 // Latitude index
1274 int latindex_south;
1275
1276 // Latitude index
1277 int latindex_north;
1278
1279 // Latitude range
1280 int latrange;
1281
1282
1283 //latitude 89-90 (both south and north) 1 value each part
1284 for (j = 0; j < dimsize1; j++) {
1285 orival[0][j] = -179.5;
1286 orival[dimsize0 - 1][j] = -179.5;
1287 }
1288
1289 //latitude 80-89 (both south and north) 45 values each part
1290 // longitude extent is 8
1291 lonextent = 8;
1292 // From 89 N to 80 N
1293 latindex_north = 1;
1294 latrange = 9;
1295 for (i = 0; i < latrange; i++)
1296 for (j = 0; j < (dimsize1 / lonextent); j++)
1297 for (k = 0; k < lonextent; k++)
1298 orival[i + latindex_north][j * lonextent + k] =
1299 -179.5 + lonextent * j;
1300
1301 // From 89 S to 80 S
1302 latindex_south = dimsize0 - 1 - latrange;
1303 for (i = 0; i < latrange; i++)
1304 for (j = 0; j < (dimsize1 / lonextent); j++)
1305 for (k = 0; k < lonextent; k++)
1306 orival[i + latindex_south][j * lonextent + k] =
1307 -179.5 + lonextent * j;
1308
1309 // From 80 N to 70 N
1310 latindex_north = latindex_north + latrange;
1311 latrange = 10;
1312 lonextent = 4;
1313
1314 for (i = 0; i < latrange; i++)
1315 for (j = 0; j < (dimsize1 / lonextent); j++)
1316 for (k = 0; k < lonextent; k++)
1317 orival[i + latindex_north][j * lonextent + k] =
1318 -179.5 + lonextent * j;
1319
1320 // From 80 S to 70 S
1321 latindex_south = latindex_south - latrange;
1322 for (i = 0; i < latrange; i++)
1323 for (j = 0; j < (dimsize1 / lonextent); j++)
1324 for (k = 0; k < lonextent; k++)
1325 orival[i + latindex_south][j * lonextent + k] =
1326 -179.5 + lonextent * j;
1327
1328 // From 70 N to 45 N
1329 latindex_north = latindex_north + latrange;
1330 latrange = 25;
1331 lonextent = 2;
1332
1333 for (i = 0; i < latrange; i++)
1334 for (j = 0; j < (dimsize1 / lonextent); j++)
1335 for (k = 0; k < lonextent; k++)
1336 orival[i + latindex_north][j * lonextent + k] =
1337 -179.5 + lonextent * j;
1338
1339 // From 70 S to 45 S
1340 latindex_south = latindex_south - latrange;
1341
1342 for (i = 0; i < latrange; i++)
1343 for (j = 0; j < (dimsize1 / lonextent); j++)
1344 for (k = 0; k < lonextent; k++)
1345 orival[i + latindex_south][j * lonextent + k] =
1346 -179.5 + lonextent * j;
1347
1348 // From 45 N to 0 N
1349 latindex_north = latindex_north + latrange;
1350 latrange = 45;
1351 lonextent = 1;
1352
1353 for (i = 0; i < latrange; i++)
1354 for (j = 0; j < (dimsize1 / lonextent); j++)
1355 for (k = 0; k < lonextent; k++)
1356 orival[i + latindex_north][j * lonextent + k] =
1357 -179.5 + lonextent * j;
1358
1359 // From 45 S to 0 S
1360 latindex_south = latindex_south - latrange;
1361 for (i = 0; i < latrange; i++)
1362 for (j = 0; j < (dimsize1 / lonextent); j++)
1363 for (k = 0; k < lonextent; k++)
1364 orival[i + latindex_south][j * lonextent + k] =
1365 -179.5 + lonextent * j;
1366
1367 for (i = 0; i < count[0]; i++) {
1368 for (j = 0; j < count[1]; j++) {
1369 val[i][j] =
1370 orival[offset[0] + step[0] * i][offset[1] + step[1] * j];
1371 }
1372 }
1373 }
1374 set_value ((dods_float32 *) (&val[0][0]), nelms);
1375
1376}
1377
1378// This function calculate zonal average(longitude is fixed) of CERES SAVG and CERES_ISCCP_DAY_LIKE.
1379void
1380HDFSPArrayGeoField::readcersavgid1 (int *offset, int *count, int *step,
1381 int nelms)
1382{
1383
1384 // Following CERES Nested grid
1385 // URL http://eosweb.larc.nasa.gov/PRODOCS/ceres/SRBAVG/Quality_Summaries/srbavg_ed2d/nestedgrid.html
1386 // or https://eosweb.larc.nasa.gov/sites/default/files/project/ceres/quality_summaries/srbavg_ed2d/nestedgrid.pdf
1387 if (fieldtype == 1) { // Calculate the latitude
1388 const int dimsize0 = 180;
1389 float32 val[count[0]];
1390 float32 orival[dimsize0];
1391
1392 for (int i = 0; i < dimsize0; i++)
1393 orival[i] = 89.5 - i;
1394
1395 for (int i = 0; i < count[0]; i++) {
1396 val[i] = orival[offset[0] + step[0] * i];
1397 }
1398 set_value ((dods_float32 *) (&val[0]), nelms);
1399
1400 }
1401
1402 if (fieldtype == 2) { // Calculate the longitude
1403
1404 // Assume the longitude is 0 in average
1405 float32 val = 0;
1406 if (nelms > 1)
1407 throw InternalErr (__FILE__, __LINE__,
1408 "the number of element must be 1");
1409 set_value ((dods_float32 *) (&val), nelms);
1410
1411 }
1412
1413}
1414
1415// Calculate CERES AVG and SYN lat and lon.
1416// We just need to retrieve lat/lon from the field.
1417void
1418HDFSPArrayGeoField::readceravgsyn (int32 * offset32, int32 * count32,
1419 int32 * step32, int nelms)
1420{
1421
1422#if 0
1423 string check_pass_fileid_key_str="H4.EnablePassFileID";
1424 bool check_pass_fileid_key = false;
1425 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
1426#endif
1427 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
1428
1429 int32 sdid = -1;
1430
1431 if(false == check_pass_fileid_key) {
1432 sdid = SDstart (const_cast < char *>(filename.c_str ()), DFACC_READ);
1433 if (sdid < 0) {
1434 ostringstream eherr;
1435 eherr << "File " << filename.c_str () << " cannot be open.";
1436 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1437 }
1438 }
1439 else
1440 sdid = sdfd;
1441
1442 int i = 0;
1443 int32 sdsid = -1;
1444
1445 int32 sdsindex = SDreftoindex (sdid, fieldref);
1446
1447 if (sdsindex == -1) {
1448 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1449 ostringstream eherr;
1450 eherr << "SDS index " << sdsindex << " is not right.";
1451 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1452 }
1453
1454 sdsid = SDselect (sdid, sdsindex);
1455 if (sdsid < 0) {
1456 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1457 ostringstream eherr;
1458 eherr << "SDselect failed.";
1459 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1460 }
1461
1462 int32 r;
1463
1464 switch (dtype) {
1465 case DFNT_INT8:
1466 case DFNT_UINT8:
1467 case DFNT_UCHAR8:
1468 case DFNT_CHAR8:
1469 case DFNT_INT16:
1470 case DFNT_UINT16:
1471 case DFNT_INT32:
1472 case DFNT_UINT32:
1473 SDendaccess (sdsid);
1474 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1475 throw InternalErr (__FILE__, __LINE__,
1476 "datatype is not float, unsupported.");
1477 case DFNT_FLOAT32:
1478 {
1479 vector<float32>val;
1480 val.resize(nelms);
1481 r = SDreaddata (sdsid, offset32, step32, count32, val.data());
1482 if (r != 0) {
1483 SDendaccess (sdsid);
1484 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1485 ostringstream eherr;
1486 eherr << "SDreaddata failed";
1487 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1488 }
1489
1490 if (fieldtype == 1) {
1491 for (i = 0; i < nelms; i++)
1492 val[i] = 90 - val[i];
1493 }
1494 if (fieldtype == 2) {
1495 for (i = 0; i < nelms; i++)
1496 if (val[i] > 180.0)
1497 val[i] = val[i] - 360.0;
1498 }
1499
1500 set_value ((dods_float32 *) val.data(), nelms);
1501 break;
1502 }
1503 case DFNT_FLOAT64:
1504 {
1505 vector<float64>val;
1506 val.resize(nelms);
1507
1508 r = SDreaddata (sdsid, offset32, step32, count32, val.data());
1509 if (r != 0) {
1510 SDendaccess (sdsid);
1511 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1512 ostringstream eherr;
1513 eherr << "SDreaddata failed";
1514 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1515 }
1516
1517 if (fieldtype == 1) {
1518 for (i = 0; i < nelms; i++)
1519 val[i] = 90 - val[i];
1520 }
1521 if (fieldtype == 2) {
1522 for (i = 0; i < nelms; i++)
1523 if (val[i] > 180.0)
1524 val[i] = val[i] - 360.0;
1525 }
1526 set_value ((dods_float64 *) val.data(), nelms);
1527 break;
1528 }
1529 default:
1530 SDendaccess (sdsid);
1531 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1532 throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
1533 }
1534
1535 r = SDendaccess (sdsid);
1536 if (r != 0) {
1537 ostringstream eherr;
1538 eherr << "SDendaccess failed.";
1539 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1540 }
1541
1542 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1543}
1544
1545// Calculate CERES ES4 and GEO lat/lon.
1546// We have to retrieve the original lat/lon and condense it from >1-D to 1-D.
1547void
1548HDFSPArrayGeoField::readceres4ig (int32 * offset32, int32 * count32,
1549 int32 * step32, int nelms)
1550{
1551#if 0
1552 string check_pass_fileid_key_str="H4.EnablePassFileID";
1553 bool check_pass_fileid_key = false;
1554 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
1555#endif
1556
1557 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
1558
1559 int32 sdid = -1;
1560
1561 if(false == check_pass_fileid_key) {
1562 sdid = SDstart (const_cast < char *>(filename.c_str ()), DFACC_READ);
1563 if (sdid < 0) {
1564 ostringstream eherr;
1565 eherr << "File " << filename.c_str () << " cannot be open.";
1566 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1567 }
1568 }
1569 else
1570 sdid = sdfd;
1571
1572
1573 int32 sdsid = -1;
1574 intn status = 0;
1575
1576 int32 sdsindex = SDreftoindex (sdid, (int32) fieldref);
1577 if (sdsindex == -1) {
1578 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1579 ostringstream eherr;
1580 eherr << "SDS index " << sdsindex << " is not right.";
1581 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1582 }
1583
1584 sdsid = SDselect (sdid, sdsindex);
1585 if (sdsid < 0) {
1586 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1587 ostringstream eherr;
1588 eherr << "SDselect failed.";
1589 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1590 }
1591
1592 int32 sdsrank = 0;
1593 int32 sds_dtype = 0;
1594 int32 n_attrs = 0;
1595 char sdsname[H4_MAX_NC_NAME];
1596 int32 dim_sizes[H4_MAX_VAR_DIMS];
1597
1598 status =
1599 SDgetinfo (sdsid, sdsname, &sdsrank, dim_sizes, &sds_dtype, &n_attrs);
1600 if (status < 0) {
1601 SDendaccess (sdsid);
1602 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1603 ostringstream eherr;
1604 eherr << "SDgetinfo failed.";
1605 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1606 }
1607
1608 vector<int32> orioffset32;
1609 vector<int32> oricount32;
1610 vector<int32> oristep32;
1611 orioffset32.resize(sdsrank);
1612 oricount32.resize(sdsrank);
1613 oristep32.resize(sdsrank);
1614
1615 int32 r;
1616
1617 switch (sds_dtype) {
1618 case DFNT_INT8:
1619 case DFNT_UINT8:
1620 case DFNT_UCHAR8:
1621 case DFNT_CHAR8:
1622 case DFNT_INT16:
1623 case DFNT_UINT16:
1624 case DFNT_INT32:
1625 case DFNT_UINT32:
1626 case DFNT_FLOAT64:
1627 SDendaccess (sdsid);
1628 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1629 throw InternalErr (__FILE__, __LINE__,
1630 "datatype is not float, unsupported.");
1631 case DFNT_FLOAT32:
1632 {
1633 vector<float32> val;
1634 val.resize(nelms);
1635 if (fieldtype == 1) {
1636 if (sptype == CER_CGEO) {
1637 if (sdsrank != 3) {
1638 SDendaccess (sdsid);
1639 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1640 throw InternalErr (__FILE__, __LINE__,
1641 "For CER_ISCCP-D2like-GEO case, lat/lon must be 3-D");
1642 }
1643 orioffset32[0] = 0;
1644
1645 // The second dimension of the original latitude should be condensed
1646 orioffset32[1] = offset32[0];
1647 orioffset32[2] = 0;
1648 oricount32[0] = 1;
1649 oricount32[1] = count32[0];
1650 oricount32[2] = 1;
1651 oristep32[0] = 1;
1652 oristep32[1] = step32[0];
1653 oristep32[2] = 1;
1654 }
1655 if (sptype == CER_ES4) {
1656 if (sdsrank != 2) {
1657 SDendaccess (sdsid);
1658 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1659 throw InternalErr (__FILE__, __LINE__,
1660 "For CER_ES4 case, lat/lon must be 2-D");
1661 }
1662 // The first dimension of the original latitude should be condensed
1663 orioffset32[0] = offset32[0];
1664 orioffset32[1] = 0;
1665 oricount32[0] = count32[0];
1666 oricount32[1] = 1;
1667 oristep32[0] = step32[0];
1668 oristep32[1] = 1;
1669 }
1670 }
1671
1672 if (fieldtype == 2) {
1673 if (sptype == CER_CGEO) {
1674 if (sdsrank != 3) {
1675 SDendaccess (sdsid);
1676 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1677 throw InternalErr (__FILE__, __LINE__,
1678 "For CER_ISCCP-D2like-GEO case, lat/lon must be 3-D");
1679 }
1680 orioffset32[0] = 0;
1681 orioffset32[2] = offset32[0];// The third dimension of the original latitude should be condensed
1682 orioffset32[1] = 0;
1683 oricount32[0] = 1;
1684 oricount32[2] = count32[0];
1685 oricount32[1] = 1;
1686 oristep32[0] = 1;
1687 oristep32[2] = step32[0];
1688 oristep32[1] = 1;
1689 }
1690 if (sptype == CER_ES4) {
1691 if (sdsrank != 2) {
1692 SDendaccess (sdsid);
1693 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1694 throw InternalErr (__FILE__, __LINE__,
1695 "For CER_ES4 case, lat/lon must be 2-D");
1696 }
1697 orioffset32[1] = offset32[0]; // The second dimension of the original latitude should be condensed
1698 orioffset32[0] = 0;
1699 oricount32[1] = count32[0];
1700 oricount32[0] = 1;
1701 oristep32[1] = step32[0];
1702 oristep32[0] = 1;
1703 }
1704 }
1705
1706 r = SDreaddata (sdsid, orioffset32.data(), oristep32.data(), oricount32.data(), val.data());
1707 if (r != 0) {
1708 SDendaccess (sdsid);
1709 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1710 ostringstream eherr;
1711 eherr << "SDreaddata failed";
1712 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1713 }
1714
1715 if (fieldtype == 1)
1716 for (int i = 0; i < nelms; i++)
1717 val[i] = 90 - val[i];
1718 if (fieldtype == 2) {
1719 // Since Panoply cannot handle the case when the longitude is jumped from 180 to -180
1720 // So turn it off and see if it works with other clients,
1721 // change my mind, should contact Panoply developer to solve this
1722 // Just check Panoply(3.2.1) with the latest release(1.9). This is no longer an issue.
1723 for (int i = 0; i < nelms; i++)
1724 if (val[i] > 180.0)
1725 val[i] = val[i] - 360.0;
1726 }
1727
1728 set_value ((dods_float32 *) val.data(), nelms);
1729 break;
1730 }
1731 default:
1732 SDendaccess (sdsid);
1733 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1734 throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
1735 }
1736
1737 r = SDendaccess (sdsid);
1738 if (r != 0) {
1739 ostringstream eherr;
1740 eherr << "SDendaccess failed.";
1741 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1742 }
1743
1744 HDFCFUtil::close_fileid(sdid,-1,-1,-1,check_pass_fileid_key);
1745}
1746
1747// Read CERES Zonal average latitude field
1748void
1749HDFSPArrayGeoField::readcerzavg (int32 * offset32, int32 * count32,
1750 int32 * step32, int nelms)
1751{
1752 if (fieldtype == 1) {
1753
1754 vector<float32>val;
1755 val.resize(nelms);
1756
1757 float32 latstep = 1.0;
1758
1759 for (int i = 0; i < nelms; i++)
1760 val[i] =
1761 89.5 - ((int) (offset32[0]) +
1762 ((int) (step32[0])) * i) * latstep;
1763 set_value ((dods_float32 *) val.data(), nelms);
1764 }
1765
1766 if (fieldtype == 2) {
1767 if (count32[0] != 1 || nelms != 1)
1768 throw InternalErr (__FILE__, __LINE__,
1769 "Longitude should only have one value for zonal mean");
1770
1771 // We don't need to specify the longitude value.
1772 float32 val = 0.;// our convention
1773 set_value ((dods_float32 *) & val, nelms);
1774 }
1775}
1776
1777
1778// Standard way of DAP handlers to pass the coordinates of the subsetted region to the handlers
1779// Return the number of elements to read.
1780int
1781HDFSPArrayGeoField::format_constraint (int *offset, int *step, int *count)
1782{
1783 long nels = 1;
1784 int id = 0;
1785
1786 Dim_iter p = dim_begin ();
1787 while (p != dim_end ()) {
1788
1789 int start = dimension_start (p, true);
1790 int stride = dimension_stride (p, true);
1791 int stop = dimension_stop (p, true);
1792
1793 // Check for illegal constraint
1794 if (start > stop) {
1795 ostringstream oss;
1796 oss << "Array/Grid hyperslab start point "<< start <<
1797 " is greater than stop point " << stop <<".";
1798 throw Error(malformed_expr, oss.str());
1799 }
1800
1801 offset[id] = start;
1802 step[id] = stride;
1803 count[id] = ((stop - start) / stride) + 1; // count of elements
1804 nels *= count[id]; // total number of values for variable
1805
1806 BESDEBUG ("h4",
1807 "=format_constraint():"
1808 << "id=" << id << " offset=" << offset[id]
1809 << " step=" << step[id]
1810 << " count=" << count[id]
1811 << endl);
1812
1813 id++;
1814 p++;
1815 }// while (p != dim_end ())
1816
1817 return nels;
1818}
1819
1820
1821
1822// Subset of latitude and longitude to follow the parameters from the DAP expression constraint
1823template < typename T >
1824void HDFSPArrayGeoField::LatLon2DSubset (T * outlatlon,
1825 int /*majordim //unused SBL 2/7/20*/,
1826 int minordim,
1827 T * latlon,
1828 int32 * offset,
1829 int32 * count,
1830 int32 * step)
1831{
1832#if 0
1833 // 'typeof' is supported only by gcc and not on OSX with --std=c++11.
1834 // jhrg 3/28/19
1835 T templatlonptr[majordim][minordim] (typeof templatlonptr) latlon;
1836#endif
1837
1838 int i = 0;
1839 int j = 0;
1840 int k = 0;
1841
1842 // do subsetting
1843 // Find the correct index
1844 const int dim0count = count[0];
1845 const int dim1count = count[1];
1846 int dim0index[dim0count];
1847 int dim1index[dim1count];
1848
1849 for (i = 0; i < count[0]; i++) // count[0] is the least changing dimension
1850 dim0index[i] = offset[0] + i * step[0];
1851
1852 for (j = 0; j < count[1]; j++)
1853 dim1index[j] = offset[1] + j * step[1];
1854
1855 // Now assign the subsetting data
1856 k = 0;
1857
1858 for (i = 0; i < count[0]; i++) {
1859 for (j = 0; j < count[1]; j++) {
1860#if 0
1861 outlatlon[k] = templatlonptr[dim0index[i]][dim1index[j]];
1862#endif
1863 outlatlon[k] = *(latlon + (dim0index[i] * minordim) + dim1index[j]); //[dim0index[i]][dim1index[j]];
1864 k++;
1865 }
1866 }
1867}
1868
static void close_fileid(int32 sdfd, int32 file_id, int32 gridfd, int32 swathfd, bool pass_fileid_key)
Definition: HDFCFUtil.cc:3650