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