bes  Updated for version 3.20.6
HDFSPArrayGeoField.cc
1 // 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 <debug.h>
16 #include "hdf.h"
17 #include "mfhdf.h"
18 #include "InternalErr.h"
19 #include <BESDebug.h>
20 #include "HDFCFUtil.h"
21 #include "HDF4RequestHandler.h"
22 
23 using namespace std;
24 using 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 
39 bool 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[0], &step[0], &count[0]);
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[0], &count32[0], &step32[0], nelms);
83  break;
84  }
85 
86  // TRMM grid
87  case TRMML3A_V6:
88  {
89  readtrmml3a_v6 (&offset32[0], &count32[0], &step32[0], nelms);
90  break;
91  }
92 
93  case TRMML3B_V6:
94  {
95  readtrmml3b_v6 (&offset32[0], &count32[0], &step32[0], nelms);
96  break;
97  }
98 
99  case TRMML3C_V6:
100  {
101  readtrmml3c_v6 (&offset32[0], &count32[0], &step32[0], nelms);
102  break;
103  }
104 
105 
106  // TRMM V7 grid
107  case TRMML3S_V7:
108  case TRMML3M_V7:
109  {
110  readtrmml3_v7 (&offset32[0], &step32[0], 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[0], &count32[0], &step32[0], nelms);
117  break;
118  }
119 
120  // CERES CER_ES4_??? products
121  case CER_ES4:
122  {
123  readceres4ig (&offset32[0], &count32[0], &step32[0], nelms);
124  break;
125  }
126 
127  // CERES CER_ISCCP-D2like-Day product
128  case CER_CDAY:
129  {
130  readcersavgid2 (&offset[0], &count[0], &step[0], nelms);
131  break;
132  }
133 
134  // CERES CER_ISCCP-D2like-GEO product
135  case CER_CGEO:
136  {
137  readceres4ig (&offset32[0], &count32[0], &step32[0], 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[0], &count[0], &step[0], nelms);
147  }
148  // When longitude is not fixed
149  else if (rank == 2) {
150  readcersavgid2 (&offset[0], &count[0], &step[0], nelms);
151  }
152  break;
153  }
154 
155  // CERES SYN Aqua products
156  case CER_SYN:
157  {
158  readceravgsyn (&offset32[0], &count32[0], &step32[0], nelms);
159  break;
160  }
161 
162  // CERES Zonal Average products
163  case CER_ZAVG:
164  {
165  readcerzavg (&offset32[0], &count32[0], &step32[0], nelms);
166  break;
167  }
168 
169  // OBPG level 2 products
170  case OBPGL2:
171  {
172  readobpgl2 (&offset32[0], &count32[0], &step32[0], nelms);
173  break;
174  }
175 
176  // OBPG level 3 products
177  case OBPGL3:
178  {
179  readobpgl3 (&offset[0], &step[0], 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 
206 void
207 HDFSPArrayGeoField::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[0], &geostep32[0], &geocount32[0], (void*)(&val[0]));
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[0], 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[0], 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[0], &geostep32[0], &geocount32[0], (void*)(&val[0]));
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[0], nelms);
323  }
324  break;
325 
326  case DFNT_INT16:
327  {
328  vector<int16> val;
329  val.resize(nelms);
330  r = SDreaddata (sdsid, &geooffset32[0], &geostep32[0], &geocount32[0], (void*)(&val[0]));
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[0], nelms);
340  }
341  break;
342 
343  case DFNT_UINT16:
344  {
345  vector<uint16>val;
346  val.resize(nelms);
347  r = SDreaddata (sdsid, &geooffset32[0], &geostep32[0], &geocount32[0], (void*)(&val[0]));
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[0], nelms);
357  }
358  break;
359  case DFNT_INT32:
360  {
361  vector<int32>val;
362  val.resize(nelms);
363  r = SDreaddata (sdsid, &geooffset32[0], &geostep32[0], &geocount32[0], (void*)(&val[0]));
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[0], nelms);
373  }
374  break;
375  case DFNT_UINT32:
376  {
377  vector<uint32>val;
378  val.resize(nelms);
379  r = SDreaddata (sdsid, &geooffset32[0], &geostep32[0], &geocount32[0], (void*)(&val[0]));
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[0], nelms);
388  }
389  break;
390  case DFNT_FLOAT32:
391  {
392  vector<float32>val;
393  val.resize(nelms);
394  r = SDreaddata (sdsid, &geooffset32[0], &geostep32[0], &geocount32[0], (void*)(&val[0]));
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[0], nelms);
404  }
405  break;
406  case DFNT_FLOAT64:
407  {
408  vector<float64>val;
409  val.resize(nelms);
410  r = SDreaddata (sdsid, &geooffset32[0], &geostep32[0], &geocount32[0], (void*)(&val[0]));
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[0], 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.
442 void
443 HDFSPArrayGeoField::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[0], 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.
478 void
479 HDFSPArrayGeoField::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[0], nelms);
509 }
510 
511 // Read TRMM level 3c(CSH) version 6 lat/lon.
512 // We follow TRMM document to retrieve the lat and lon.
513 void
514 HDFSPArrayGeoField::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[0], nelms);
544 }
545 // Read latlon of TRMM version 7 products.
546 // Lat/lon parameters can be retrieved from attribute GridHeader.
547 void
548 HDFSPArrayGeoField::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[0]);
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[0], 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.
667 void
668 HDFSPArrayGeoField::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[0]);
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[0]);
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[0],num_scan_data,num_pixel_data,&interp_val[0],
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[0], 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.
986 void
987 HDFSPArrayGeoField::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[0], 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[0], 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
1238 void
1239 HDFSPArrayGeoField::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.
1379 void
1380 HDFSPArrayGeoField::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.
1417 void
1418 HDFSPArrayGeoField::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[0]);
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[0], 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[0]);
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[0], 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.
1547 void
1548 HDFSPArrayGeoField::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[0], &oristep32[0], &oricount32[0], &val[0]);
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[0], 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
1748 void
1749 HDFSPArrayGeoField::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[0], 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.
1780 int
1781 HDFSPArrayGeoField::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
1823 template < typename T >
1824 void HDFSPArrayGeoField::LatLon2DSubset (T * outlatlon,
1825  int majordim,
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 
libdap
Definition: BESDapFunctionResponseCache.h:35
Error
HDFCFUtil::close_fileid
static void close_fileid(int32 sdfd, int32 file_id, int32 gridfd, int32 swathfd, bool pass_fileid_key)
Definition: HDFCFUtil.cc:3685