bes  Updated for version 3.20.6
HDF5CFUtil.cc
Go to the documentation of this file.
1 // This file is part of the hdf5_handler implementing for the CF-compliant
2 // Copyright (c) 2011-2016 The HDF Group, Inc. and OPeNDAP, Inc.
3 //
4 // This is free software; you can redistribute it and/or modify it under the
5 // terms of the GNU Lesser General Public License as published by the Free
6 // Software Foundation; either version 2.1 of the License, or (at your
7 // option) any later version.
8 //
9 // This software is distributed in the hope that it will be useful, but
10 // WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
11 // or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
12 // License for more details.
13 //
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
17 //
18 // You can contact OPeNDAP, Inc. at PO Box 112, Saunderstown, RI. 02874-0112.
19 // You can contact The HDF Group, Inc. at 1800 South Oak Street,
20 // Suite 203, Champaign, IL 61820
21 
31 
32 #include "HDF5CFUtil.h"
33 //#include "HE5GridPara.h"
34 #include "HDF5RequestHandler.h"
35 #include <set>
36 #include <sstream>
37 #include <algorithm>
38 #include <stdlib.h>
39 #include <unistd.h>
40 #include <math.h>
41 #include<InternalErr.h>
42 
43 using namespace std;
44 using namespace libdap;
45 // For using GCTP to calculate the lat/lon
46 extern "C" {
47 int hinv_init(int insys, int inzone, double *inparm, int indatum, char *fn27, char *fn83, int *iflg, int (*hinv_trans[])(double, double, double*, double*));
48 
49 int hfor_init(int outsys, int outzone, double *outparm, int outdatum, char *fn27, char *fn83, int *iflg, int (*hfor_trans[])(double, double, double *, double *));
50 
51 }
52 
53 H5DataType
55 {
56  size_t size = 0;
57  int sign = -2;
58 
59  switch (H5Tget_class(h5_type_id)) {
60 
61  case H5T_INTEGER:
62  size = H5Tget_size(h5_type_id);
63  sign = H5Tget_sign(h5_type_id);
64 
65  if (size == 1) { // Either signed char or unsigned char
66  if (sign == H5T_SGN_2)
67  return H5CHAR;
68  else
69  return H5UCHAR;
70  }
71  else if (size == 2) {
72  if (sign == H5T_SGN_2)
73  return H5INT16;
74  else
75  return H5UINT16;
76  }
77  else if (size == 4) {
78  if (sign == H5T_SGN_2)
79  return H5INT32;
80  else
81  return H5UINT32;
82  }
83  else if (size == 8) {
84  if (sign == H5T_SGN_2)
85  return H5INT64;
86  else
87  return H5UINT64;
88  }
89  else return H5UNSUPTYPE;
90 
91  case H5T_FLOAT:
92  size = H5Tget_size(h5_type_id);
93 
94  if (size == 4) return H5FLOAT32;
95  else if (size == 8) return H5FLOAT64;
96  else return H5UNSUPTYPE;
97 
98  case H5T_STRING:
99  if (H5Tis_variable_str(h5_type_id))
100  return H5VSTRING;
101  else return H5FSTRING;
102 
103  case H5T_REFERENCE:
104  return H5REFERENCE;
105 
106 
107  case H5T_COMPOUND:
108  return H5COMPOUND;
109 
110  case H5T_ARRAY:
111  return H5ARRAY;
112 
113  default:
114  return H5UNSUPTYPE;
115 
116  }
117 }
118 
119 size_t HDF5CFUtil::H5_numeric_atomic_type_size(H5DataType h5type) {
120 
121  switch(h5type) {
122  case H5CHAR:
123  case H5UCHAR:
124  return 1;
125  case H5INT16:
126  case H5UINT16:
127  return 2;
128  case H5INT32:
129  case H5UINT32:
130  case H5FLOAT32:
131  return 4;
132  case H5FLOAT64:
133  case H5INT64:
134  case H5UINT64:
135  return 8;
136  default:
137  throw InternalErr(__FILE__,__LINE__,"This routine doesn't support to return the size of this datatype");
138  }
139 
140 }
141 
142 #if 0
143 bool HDF5CFUtil::use_lrdata_mem_cache(H5DataType h5type, CVType cvtype, bool islatlon ) {
144  if(h5type != H5CHAR && h5type !=H5UCHAR && h5type!=H5INT16 && h5type !=H5UINT16 &&
145  h5type != H5INT32 && h5type !=H5UINT32 && h5type !=H5FLOAT32 && h5type!=H5FLOAT64 &&
146  h5type != H5INT64 && h5type !=H5UINT64)
147  return false;
148  else {
149  if(cvtype != CV_UNSUPPORTED)
150  return true;
151  // TODO; if varpath is specified by the user, this should return true!
152  else if(varpath == "")
153  return false;
154  else
155  return false;
156 
157  }
158 
159 }
160 #endif
161 
162 // Check if we cna use data memory cache
163 // TODO: This functio is not used, we will check it in the next release.
164 bool HDF5CFUtil::use_data_mem_cache(H5DataType h5type, CVType cvtype, const string &varpath) {
165  if(h5type != H5CHAR && h5type !=H5UCHAR && h5type!=H5INT16 && h5type !=H5UINT16 &&
166  h5type != H5INT32 && h5type !=H5UINT32 && h5type !=H5FLOAT32 && h5type!=H5FLOAT64 &&
167  h5type != H5INT64 && h5type !=H5UINT64)
168  return false;
169  else {
170  if(cvtype != CV_UNSUPPORTED)
171  return true;
172  // TODO; if varpath is specified by the user, this should return true!
173  else if(varpath == "")
174  return false;
175  else
176  return false;
177 
178  }
179 
180 }
181 
182 bool HDF5CFUtil::cf_strict_support_type(H5DataType dtype) {
183  if ((H5UNSUPTYPE == dtype)||(H5REFERENCE == dtype)
184  || (H5COMPOUND == dtype) || (H5ARRAY == dtype))
185  // Important comments for the future work.
186  // Try to suport 64-bit integer for DAP4 CF, check the starting code at get_dmr_64bit_int()
187  //"|| (H5INT64 == dtype) ||(H5UINT64 == dtype))"
188  return false;
189  else
190  return true;
191 }
192 
193 bool HDF5CFUtil::cf_dap2_support_numeric_type(H5DataType dtype) {
194  if ((H5UNSUPTYPE == dtype)||(H5REFERENCE == dtype)
195  || (H5COMPOUND == dtype) || (H5ARRAY == dtype)
196  || (H5INT64 == dtype) ||(H5UINT64 == dtype)
197  || (H5FSTRING == dtype) ||(H5VSTRING == dtype))
198  return false;
199  else
200  return true;
201 }
202 
203 string HDF5CFUtil::obtain_string_after_lastslash(const string s) {
204 
205  string ret_str="";
206  size_t last_fslash_pos = s.find_last_of("/");
207  if (string::npos != last_fslash_pos &&
208  last_fslash_pos != (s.size()-1))
209  ret_str=s.substr(last_fslash_pos+1);
210  return ret_str;
211 }
212 
213 string HDF5CFUtil::obtain_string_before_lastslash(const string & s) {
214 
215  string ret_str="";
216  size_t last_fslash_pos = s.find_last_of("/");
217  if (string::npos != last_fslash_pos)
218  ret_str=s.substr(0,last_fslash_pos+1);
219  return ret_str;
220 
221 }
222 
223 string HDF5CFUtil::trim_string(hid_t ty_id,const string s, int num_sect, size_t sect_size, vector<size_t>& sect_newsize) {
224 
225  string temp_sect_str = "";
226  string temp_sect_newstr = "";
227  string final_str="";
228 
229  for (int i = 0; i < num_sect; i++) {
230 
231  if (i != (num_sect-1))
232  temp_sect_str = s.substr(i*sect_size,sect_size);
233  else
234  temp_sect_str = s.substr((num_sect-1)*sect_size,s.size()-(num_sect-1)*sect_size);
235 
236  size_t temp_pos = 0;
237 
238  if (H5T_STR_NULLTERM == H5Tget_strpad(ty_id))
239  temp_pos = temp_sect_str.find_first_of('\0');
240  else if (H5T_STR_SPACEPAD == H5Tget_strpad(ty_id))
241  temp_pos = temp_sect_str.find_last_not_of(' ')+1;
242  else temp_pos = temp_sect_str.find_last_not_of('0')+1;// NULL PAD
243 
244  if (temp_pos != string::npos) {
245 
246  // Here I only add a space at the end of each section for the SPACEPAD case,
247  // but not for NULL TERM
248  // or NULL PAD. Two reasons: C++ string doesn't need NULL TERM.
249  // We don't know and don't see any NULL PAD applications for NASA data.
250  if (H5T_STR_SPACEPAD == H5Tget_strpad(ty_id)) {
251 
252  if (temp_pos == temp_sect_str.size())
253  temp_sect_newstr = temp_sect_str +" ";
254  else
255  temp_sect_newstr = temp_sect_str.substr(0,temp_pos+1);
256 
257  sect_newsize[i] = temp_pos +1;
258  }
259  else {
260  temp_sect_newstr = temp_sect_str.substr(0,temp_pos);
261  sect_newsize[i] = temp_pos ;
262  }
263 
264  }
265 
266  else {// NULL is not found, adding a NULL at the end of this string.
267 
268  temp_sect_newstr = temp_sect_str;
269 
270  // Here I only add a space for the SPACEPAD, but not for NULL TERM
271  // or NULL PAD. Two reasons: C++ string doesn't need NULL TERM.
272  // We don't know and don't see any NULL PAD applications for NASA data.
273  if (H5T_STR_SPACEPAD == H5Tget_strpad(ty_id)) {
274  temp_sect_newstr.resize(temp_sect_str.size()+1);
275  temp_sect_newstr.append(1,' ');
276  sect_newsize[i] = sect_size + 1;
277  }
278  else
279  sect_newsize[i] = sect_size;
280  }
281  final_str+=temp_sect_newstr;
282  }
283 
284  return final_str;
285 }
286 
287 string HDF5CFUtil::remove_substrings(string str,const string &substr) {
288 
289  string::size_type i = str.find(substr);
290  while (i != std::string::npos) {
291  str.erase(i, substr.size());
292  i = str.find(substr, i);
293  }
294  return str;
295 }
296 // Obtain the unique name for the clashed names and save it to set namelist.
297 void HDF5CFUtil::gen_unique_name(string &str,set<string>& namelist, int&clash_index) {
298 
299  pair<set<string>::iterator,bool> ret;
300  string newstr = "";
301  stringstream sclash_index;
302  sclash_index << clash_index;
303  newstr = str + sclash_index.str();
304 
305  ret = namelist.insert(newstr);
306  if (false == ret.second) {
307  clash_index++;
308  gen_unique_name(str,namelist,clash_index);
309  }
310  else
311  str = newstr;
312 }
313 
314 void
315 HDF5CFUtil::Split_helper(vector<string> &tokens, const string &text, const char sep)
316 {
317  string::size_type start = 0, end = 0;
318  while ((end = text.find(sep, start)) != string::npos) {
319  tokens.push_back(text.substr(start, end - start));
320  start = end + 1;
321  }
322  tokens.push_back(text.substr(start));
323 }
324 
325 
326 // From a string separated by a separator to a list of string,
327 // for example, split "ab,c" to {"ab","c"}
328 void
329 HDF5CFUtil::Split(const char *s, int len, char sep, std::vector<std::string> &names)
330 {
331  names.clear();
332  for (int i = 0, j = 0; j <= len; ++j) {
333  if ((j == len && len) || s[j] == sep) {
334  string elem(s + i, j - i);
335  names.push_back(elem);
336  i = j + 1;
337  continue;
338  }
339  }
340 }
341 
342 
343 // Assume sz is Null terminated string.
344 void
345 HDF5CFUtil::Split(const char *sz, char sep, std::vector<std::string> &names)
346 {
347  Split(sz, (int)strlen(sz), sep, names);
348 }
349 
350 void HDF5CFUtil::parser_gpm_l3_gridheader(const vector<char>& value,
351  int& latsize, int&lonsize,
352  float& lat_start, float& lon_start,
353  float& lat_res, float& lon_res,
354  bool check_reg_orig ){
355 
356  float lat_north = 0.;
357  float lat_south = 0.;
358  float lon_east = 0.;
359  float lon_west = 0.;
360 
361  vector<string> ind_elems;
362  char sep='\n';
363  HDF5CFUtil::Split(&value[0],sep,ind_elems);
364 
365  // The number of elements in the GridHeader is 9. The string vector will add a leftover. So the size should be 10.
366  if(ind_elems.size()!=10)
367  throw InternalErr(__FILE__,__LINE__,"The number of elements in the TRMM level 3 GridHeader is not right.");
368 
369  if(false == check_reg_orig) {
370  if (0 != ind_elems[1].find("Registration=CENTER"))
371  throw InternalErr(__FILE__,__LINE__,"The TRMM grid registration is not center.");
372  }
373 
374  if (0 == ind_elems[2].find("LatitudeResolution")){
375 
376  size_t equal_pos = ind_elems[2].find_first_of('=');
377  if(string::npos == equal_pos)
378  throw InternalErr(__FILE__,__LINE__,"Cannot find latitude resolution for TRMM level 3 products");
379 
380  size_t scolon_pos = ind_elems[2].find_first_of(';');
381  if(string::npos == scolon_pos)
382  throw InternalErr(__FILE__,__LINE__,"Cannot find latitude resolution for TRMM level 3 products");
383  if (equal_pos < scolon_pos){
384 
385  string latres_str = ind_elems[2].substr(equal_pos+1,scolon_pos-equal_pos-1);
386  lat_res = strtof(latres_str.c_str(),NULL);
387  }
388  else
389  throw InternalErr(__FILE__,__LINE__,"latitude resolution is not right for TRMM level 3 products");
390  }
391  else
392  throw InternalErr(__FILE__,__LINE__,"The TRMM grid LatitudeResolution doesn't exist.");
393 
394  if (0 == ind_elems[3].find("LongitudeResolution")){
395 
396  size_t equal_pos = ind_elems[3].find_first_of('=');
397  if(string::npos == equal_pos)
398  throw InternalErr(__FILE__,__LINE__,"Cannot find longitude resolution for TRMM level 3 products");
399 
400  size_t scolon_pos = ind_elems[3].find_first_of(';');
401  if(string::npos == scolon_pos)
402  throw InternalErr(__FILE__,__LINE__,"Cannot find longitude resolution for TRMM level 3 products");
403  if (equal_pos < scolon_pos){
404  string lonres_str = ind_elems[3].substr(equal_pos+1,scolon_pos-equal_pos-1);
405  lon_res = strtof(lonres_str.c_str(),NULL);
406  }
407  else
408  throw InternalErr(__FILE__,__LINE__,"longitude resolution is not right for TRMM level 3 products");
409  }
410  else
411  throw InternalErr(__FILE__,__LINE__,"The TRMM grid LongitudeResolution doesn't exist.");
412 
413  if (0 == ind_elems[4].find("NorthBoundingCoordinate")){
414 
415  size_t equal_pos = ind_elems[4].find_first_of('=');
416  if(string::npos == equal_pos)
417  throw InternalErr(__FILE__,__LINE__,"Cannot find latitude resolution for TRMM level 3 products");
418 
419  size_t scolon_pos = ind_elems[4].find_first_of(';');
420  if(string::npos == scolon_pos)
421  throw InternalErr(__FILE__,__LINE__,"Cannot find latitude resolution for TRMM level 3 products");
422  if (equal_pos < scolon_pos){
423  string north_bounding_str = ind_elems[4].substr(equal_pos+1,scolon_pos-equal_pos-1);
424  lat_north = strtof(north_bounding_str.c_str(),NULL);
425  }
426  else
427  throw InternalErr(__FILE__,__LINE__,"NorthBoundingCoordinate is not right for TRMM level 3 products");
428 
429  }
430  else
431  throw InternalErr(__FILE__,__LINE__,"The TRMM grid NorthBoundingCoordinate doesn't exist.");
432 
433  if (0 == ind_elems[5].find("SouthBoundingCoordinate")){
434 
435  size_t equal_pos = ind_elems[5].find_first_of('=');
436  if(string::npos == equal_pos)
437  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for TRMM level 3 products");
438 
439  size_t scolon_pos = ind_elems[5].find_first_of(';');
440  if(string::npos == scolon_pos)
441  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for TRMM level 3 products");
442  if (equal_pos < scolon_pos){
443  string lat_south_str = ind_elems[5].substr(equal_pos+1,scolon_pos-equal_pos-1);
444  lat_south = strtof(lat_south_str.c_str(),NULL);
445  }
446  else
447  throw InternalErr(__FILE__,__LINE__,"south bound coordinate is not right for TRMM level 3 products");
448 
449  }
450  else
451  throw InternalErr(__FILE__,__LINE__,"The TRMM grid SouthBoundingCoordinate doesn't exist.");
452 
453  if (0 == ind_elems[6].find("EastBoundingCoordinate")){
454 
455  size_t equal_pos = ind_elems[6].find_first_of('=');
456  if(string::npos == equal_pos)
457  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for TRMM level 3 products");
458 
459  size_t scolon_pos = ind_elems[6].find_first_of(';');
460  if(string::npos == scolon_pos)
461  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for TRMM level 3 products");
462  if (equal_pos < scolon_pos){
463  string lon_east_str = ind_elems[6].substr(equal_pos+1,scolon_pos-equal_pos-1);
464  lon_east = strtof(lon_east_str.c_str(),NULL);
465  }
466  else
467  throw InternalErr(__FILE__,__LINE__,"south bound coordinate is not right for TRMM level 3 products");
468 
469  }
470  else
471  throw InternalErr(__FILE__,__LINE__,"The TRMM grid EastBoundingCoordinate doesn't exist.");
472 
473  if (0 == ind_elems[7].find("WestBoundingCoordinate")){
474 
475  size_t equal_pos = ind_elems[7].find_first_of('=');
476  if(string::npos == equal_pos)
477  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for TRMM level 3 products");
478 
479  size_t scolon_pos = ind_elems[7].find_first_of(';');
480  if(string::npos == scolon_pos)
481  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for TRMM level 3 products");
482  if (equal_pos < scolon_pos){
483  string lon_west_str = ind_elems[7].substr(equal_pos+1,scolon_pos-equal_pos-1);
484  lon_west = strtof(lon_west_str.c_str(),NULL);
485  }
486  else
487  throw InternalErr(__FILE__,__LINE__,"south bound coordinate is not right for TRMM level 3 products");
488 
489  }
490  else
491  throw InternalErr(__FILE__,__LINE__,"The TRMM grid WestBoundingCoordinate doesn't exist.");
492 
493  if (false == check_reg_orig) {
494  if (0 != ind_elems[8].find("Origin=SOUTHWEST"))
495  throw InternalErr(__FILE__,__LINE__,"The TRMM grid origin is not SOUTHWEST.");
496  }
497 
498  // Since we only treat the case when the Registration is center, so the size should be the
499  // regular number size - 1.
500  latsize =(int)((lat_north-lat_south)/lat_res);
501  lonsize =(int)((lon_east-lon_west)/lon_res);
502  lat_start = lat_south;
503  lon_start = lon_west;
504 }
505 
506 void HDF5CFUtil::close_fileid(hid_t file_id,bool pass_fileid) {
507  if(false == pass_fileid) {
508  if(file_id != -1)
509  H5Fclose(file_id);
510  }
511 
512 }
513 
514 // Somehow the conversion of double to c++ string with sprintf causes the memory error in
515 // the testing code. I used the following code to do the conversion. Most part of the code
516 // in reverse, int_to_str and dtoa are adapted from geeksforgeeks.org
517 
518 // reverses a string 'str' of length 'len'
519 void HDF5CFUtil::rev_str(char *str, int len)
520 {
521  int i=0;
522  int j=len-1;
523  int temp = 0;
524  while (i<j)
525  {
526  temp = str[i];
527  str[i] = str[j];
528  str[j] = temp;
529  i++;
530  j--;
531  }
532 }
533 
534 // Converts a given integer x to string str[]. d is the number
535 // of digits required in output. If d is more than the number
536 // of digits in x, then 0s are added at the beginning.
537 int HDF5CFUtil::int_to_str(int x, char str[], int d)
538 {
539  int i = 0;
540  while (x)
541  {
542  str[i++] = (x%10) + '0';
543  x = x/10;
544  }
545 
546  // If number of digits required is more, then
547  // add 0s at the beginning
548  while (i < d)
549  str[i++] = '0';
550 
551  rev_str(str, i);
552  str[i] = '\0';
553  return i;
554 }
555 
556 // Converts a double floating point number to string.
557 void HDF5CFUtil::dtoa(double n, char *res, int afterpoint)
558 {
559  // Extract integer part
560  int ipart = (int)n;
561 
562  // Extract the double part
563  double fpart = n - (double)ipart;
564 
565  // convert integer part to string
566  int i = int_to_str(ipart, res, 0);
567 
568  // check for display option after point
569  if (afterpoint != 0)
570  {
571  res[i] = '.'; // add dot
572 
573  // Get the value of fraction part upto given no.
574  // of points after dot. The third parameter is needed
575  // to handle cases like 233.007
576  fpart = fpart * pow(10, afterpoint);
577 
578  // A round-error of 1 is found when casting to the integer for some numbers.
579  // We need to correct it.
580  int final_fpart = (int)fpart;
581  if(fpart -(int)fpart >0.5)
582  final_fpart = (int)fpart +1;
583  int_to_str(final_fpart, res + i + 1, afterpoint);
584  }
585 }
586 
587 
588 // Used to generate EOS geolocation cache name
589 string HDF5CFUtil::get_int_str(int x) {
590 
591  string str;
592  if(x > 0 && x <10)
593  str.push_back(x+'0');
594 
595  else if (x >10 && x<100) {
596  str.push_back(x/10+'0');
597  str.push_back(x%10+'0');
598  }
599  else {
600  int num_digit = 0;
601  int abs_x = (x<0)?-x:x;
602  while(abs_x/=10)
603  num_digit++;
604  if(x<=0)
605  num_digit++;
606  vector<char> buf;
607  buf.resize(num_digit);
608  snprintf(&buf[0],num_digit,"%d",x);
609  str.assign(&buf[0]);
610 
611  }
612 
613  return str;
614 
615 }
616 
617 //Used to generate EOS5 lat/lon cache name
618 string HDF5CFUtil::get_double_str(double x,int total_digit,int after_point) {
619 
620  string str;
621  if(x!=0) {
622  //char res[total_digit];
623  vector<char> res;
624  res.resize(total_digit);
625  for(int i = 0; i<total_digit;i++)
626  res[i] = '\0';
627  if (x<0) {
628  str.push_back('-');
629  dtoa(-x,&res[0],after_point);
630  for(int i = 0; i<total_digit;i++) {
631  if(res[i] != '\0')
632  str.push_back(res[i]);
633  }
634  }
635  else {
636  dtoa(x, &res[0], after_point);
637  //dtoa(x, res, after_point);
638  for(int i = 0; i<total_digit;i++) {
639  if(res[i] != '\0')
640  str.push_back(res[i]);
641  }
642  }
643 
644  }
645  else
646  str.push_back('0');
647 
648  return str;
649 
650 }
651 
652 
653 // This function is adapted from the HDF-EOS library.
654 int GDij2ll(int projcode, int zonecode, double projparm[],
655  int spherecode, int xdimsize, int ydimsize,
656  double upleftpt[], double lowrightpt[],
657  int npnts, int row[], int col[],
658  double longitude[], double latitude[], EOS5GridPRType pixcen, EOS5GridOriginType pixcnr)
659 {
660 
661  int errorcode = 0;
662  // In the original GCTP library, the function pointer names should be inv_trans and for_trans.
663  // However, since Hyrax supports both GDAL(including the HDF-EOS2 driver) and HDF handlers,
664  // on some machines, the functions inside the HDF-EOS2 driver will be called in run-time and wrong lat/lon
665  // values may be generated. To avoid, we change the function pointer names inside the GCTP library.
666  int(*hinv_trans[100]) (double,double,double*,double*);
667  int(*hfor_trans[100]) (double,double,double*,double*); /* GCTP function pointer */
668  double arg1, arg2;
669  double pixadjX = 0.; /* Pixel adjustment (x) */
670  double pixadjY = 0.; /* Pixel adjustment (y) */
671  double lonrad0 = 0.; /* Longitude in radians of upleft point */
672  double latrad0 = 0.; /* Latitude in radians of upleft point */
673  double scaleX = 0.; /* X scale factor */
674  double scaleY = 0.; /* Y scale factor */
675  double lonrad = 0.; /* Longitude in radians of point */
676  double latrad = 0.; /* Latitude in radians of point */
677  double xMtr0, yMtr0, xMtr1, yMtr1;
678 
679 
680 
681  /* Compute adjustment of position within pixel */
682  /* ------------------------------------------- */
683  if (pixcen == HE5_HDFE_CENTER)
684  {
685  /* Pixel defined at center */
686  /* ----------------------- */
687  pixadjX = 0.5;
688  pixadjY = 0.5;
689  }
690  else
691  {
692  switch (pixcnr)
693  {
694 
695  case HE5_HDFE_GD_UL:
696  {
697  /* Pixel defined at upper left corner */
698  /* ---------------------------------- */
699  pixadjX = 0.0;
700  pixadjY = 0.0;
701  break;
702  }
703 
704  case HE5_HDFE_GD_UR:
705  {
706  /* Pixel defined at upper right corner */
707  /* ----------------------------------- */
708  pixadjX = 1.0;
709  pixadjY = 0.0;
710  break;
711  }
712 
713  case HE5_HDFE_GD_LL:
714  {
715  /* Pixel defined at lower left corner */
716  /* ---------------------------------- */
717  pixadjX = 0.0;
718  pixadjY = 1.0;
719  break;
720  }
721 
722  case HE5_HDFE_GD_LR:
723  {
724 
725  /* Pixel defined at lower right corner */
726  /* ----------------------------------- */
727  pixadjX = 1.0;
728  pixadjY = 1.0;
729  break;
730  }
731 
732  default:
733  {
734  throw InternalErr(__FILE__,__LINE__,"Unknown pixel corner to retrieve lat/lon from a grid.");
735  }
736  }
737  }
738 
739 
740 
741  // If projection not GEO or BCEA call GCTP initialization routine
742  if (projcode != HE5_GCTP_GEO && projcode != HE5_GCTP_BCEA)
743  {
744 
745  scaleX = (lowrightpt[0] - upleftpt[0]) / xdimsize;
746  scaleY = (lowrightpt[1] - upleftpt[1]) / ydimsize;
747  string eastFile = HDF5RequestHandler::get_stp_east_filename();
748  string northFile = HDF5RequestHandler::get_stp_north_filename();
749 
750  hinv_init(projcode, zonecode, projparm, spherecode, (char*)eastFile.c_str(), (char*)northFile.c_str(),
751  &errorcode, hinv_trans);
752 
753 
754  /* Report error if any */
755  /* ------------------- */
756  if (errorcode != 0)
757  {
758  throw InternalErr(__FILE__,__LINE__,"GCTP hinv_init Error to retrieve lat/lon from a grid.");
759 
760  }
761  else
762  {
763  /* For each point ... */
764  /* ------------------ */
765  for (int i = 0; i < npnts; i++)
766  {
767  /* Convert from meters to lon/lat (radians) using GCTP */
768  /* --------------------------------------------------- */
769 #if 0
770  /*errorcode = hinv_trans[projcode] ((col[i] + pixadjX) * scaleX + upleftpt[0], (row[i] + pixadjY) * scaleY + upleftpt[1], &lonrad, &latrad);*/
771 #endif
772 
773  /* modified previous line to the following for the linux64 with -fPIC in cmpilation.
774  Whithout the change col[] and row[] values are ridiclous numbers, resulting a strange
775  number (very big) for arg1 and arg2. But with (int) typecast they become normal integers,
776  resulting in a acceptable values for arg1 and arg2. The problem was discovered during the
777  lat/lon geolocating of an hdfeos5 file with 64-bit hadview plug-in, developped for linux64.
778  */
779  arg1 = (((int)col[i] + pixadjX) * scaleX + upleftpt[0]);
780  arg2 = (((int)row[i] + pixadjY) * scaleY + upleftpt[1]);
781  errorcode = hinv_trans[projcode] (arg1, arg2, &lonrad, &latrad);
782 
783  /* Report error if any */
784  /* ------------------- */
785  if (errorcode != 0)
786  {
787 
788  if(projcode == HE5_GCTP_LAMAZ) {
789  longitude[i] = 1.0e51;
790  latitude[i] = 1.0e51;
791  }
792  else
793  throw InternalErr(__FILE__,__LINE__,"GCTP hinv_trans Error to retrieve lat/lon from a grid.");
794 
795  }
796  else
797  {
798 
799  /* Convert from radians to decimal degrees */
800  /* --------------------------------------- */
801  longitude[i] = HE5_EHconvAng(lonrad, HE5_HDFE_RAD_DEG);
802  latitude[i] = HE5_EHconvAng(latrad, HE5_HDFE_RAD_DEG);
803 
804  }
805  }
806  }
807  }
808  else if (projcode == HE5_GCTP_BCEA)
809  {
810  /* BCEA projection */
811  /* -------------- */
812 
813  /* Note: upleftpt and lowrightpt are in packed degrees, so they
814  must be converted to meters for this projection */
815 
816  /* Initialize forward transformation */
817  /* --------------------------------- */
818  hfor_init(projcode, zonecode, projparm, spherecode, NULL, NULL,&errorcode, hfor_trans);
819 
820  /* Report error if any */
821  /* ------------------- */
822  if (errorcode != 0)
823  {
824  throw InternalErr(__FILE__,__LINE__,"GCTP hfor_init Error to retrieve lat/lon from a grid.");
825  }
826 
827  /* Convert upleft and lowright X coords from DMS to radians */
828  /* -------------------------------------------------------- */
829  lonrad0 =HE5_EHconvAng(upleftpt[0], HE5_HDFE_DMS_RAD);
830  lonrad = HE5_EHconvAng(lowrightpt[0], HE5_HDFE_DMS_RAD);
831 
832  /* Convert upleft and lowright Y coords from DMS to radians */
833  /* -------------------------------------------------------- */
834  latrad0 = HE5_EHconvAng(upleftpt[1], HE5_HDFE_DMS_RAD);
835  latrad = HE5_EHconvAng(lowrightpt[1], HE5_HDFE_DMS_RAD);
836 
837  /* Convert form lon/lat to meters(or whatever unit is, i.e unit
838  of r_major and r_minor) using GCTP */
839  /* ----------------------------------------- */
840  errorcode = hfor_trans[projcode] (lonrad0, latrad0, &xMtr0, &yMtr0);
841 
842  /* Report error if any */
843  if (errorcode != 0)
844  {
845  throw InternalErr(__FILE__,__LINE__,"GCTP hfor_trans Error to retrieve lat/lon from a grid.");
846 
847  }
848 
849  /* Convert from lon/lat to meters or whatever unit is, i.e unit
850  of r_major and r_minor) using GCTP */
851  /* ----------------------------------------- */
852  errorcode = hfor_trans[projcode] (lonrad, latrad, &xMtr1, &yMtr1);
853 
854  /* Report error if any */
855  if (errorcode != 0)
856  {
857  throw InternalErr(__FILE__,__LINE__,"GCTP hfor_trans Error to retrieve lat/lon from a grid.");
858  }
859 
860  /* Compute x scale factor */
861  /* ---------------------- */
862  scaleX = (xMtr1 - xMtr0) / xdimsize;
863 
864  /* Compute y scale factor */
865  /* ---------------------- */
866  scaleY = (yMtr1 - yMtr0) / ydimsize;
867 
868  /* Initialize inverse transformation */
869  /* --------------------------------- */
870  hinv_init(projcode, zonecode, projparm, spherecode, NULL, NULL, &errorcode, hinv_trans);
871  /* Report error if any */
872  /* ------------------- */
873  if (errorcode != 0)
874  {
875  throw InternalErr(__FILE__,__LINE__,"GCTP hinv_init Error to retrieve lat/lon from a grid.");
876  }
877  /* For each point ... */
878  /* ------------------ */
879  for (int i = 0; i < npnts; i++)
880  {
881  /* Convert from meters (or any units that r_major and
882  r_minor has) to lon/lat (radians) using GCTP */
883  /* --------------------------------------------------- */
884  errorcode = hinv_trans[projcode] (
885  (col[i] + pixadjX) * scaleX + xMtr0,
886  (row[i] + pixadjY) * scaleY + yMtr0,
887  &lonrad, &latrad);
888 
889  /* Report error if any */
890  /* ------------------- */
891  if (errorcode != 0)
892  {
893 #if 0
894  /* status = -1;
895  sprintf(errbuf, "GCTP Error: %li\n", errorcode);
896  H5Epush(__FILE__, "HE5_GDij2ll", __LINE__, H5E_ARGS, H5E_BADVALUE , errbuf);
897  HE5_EHprint(errbuf, __FILE__, __LINE__);
898  return (status); */
899 #endif
900  longitude[i] = 1.0e51; /* PGSd_GCT_IN_ERROR */
901  latitude[i] = 1.0e51; /* PGSd_GCT_IN_ERROR */
902  }
903 
904  /* Convert from radians to decimal degrees */
905  /* --------------------------------------- */
906  longitude[i] = HE5_EHconvAng(lonrad, HE5_HDFE_RAD_DEG);
907  latitude[i] = HE5_EHconvAng(latrad, HE5_HDFE_RAD_DEG);
908  }
909  }
910 
911  else if (projcode == HE5_GCTP_GEO)
912  {
913  /* GEO projection */
914  /* -------------- */
915 
916  /*
917  * Note: lonrad, lonrad0, latrad, latrad0 are actually in degrees for
918  * the GEO projection case.
919  */
920 
921 
922  /* Convert upleft and lowright X coords from DMS to degrees */
923  /* -------------------------------------------------------- */
924  lonrad0 = HE5_EHconvAng(upleftpt[0], HE5_HDFE_DMS_DEG);
925  lonrad = HE5_EHconvAng(lowrightpt[0], HE5_HDFE_DMS_DEG);
926 
927  /* Compute x scale factor */
928  /* ---------------------- */
929  scaleX = (lonrad - lonrad0) / xdimsize;
930 
931  /* Convert upleft and lowright Y coords from DMS to degrees */
932  /* -------------------------------------------------------- */
933  latrad0 = HE5_EHconvAng(upleftpt[1], HE5_HDFE_DMS_DEG);
934  latrad = HE5_EHconvAng(lowrightpt[1], HE5_HDFE_DMS_DEG);
935 
936  /* Compute y scale factor */
937  /* ---------------------- */
938  scaleY = (latrad - latrad0) / ydimsize;
939 
940  /* For each point ... */
941  /* ------------------ */
942  for (int i = 0; i < npnts; i++)
943  {
944  /* Convert to lon/lat (decimal degrees) */
945  /* ------------------------------------ */
946  longitude[i] = (col[i] + pixadjX) * scaleX + lonrad0;
947  latitude[i] = (row[i] + pixadjY) * scaleY + latrad0;
948  }
949  }
950 
951 
952 #if 0
953  hinv_init(projcode, zonecode, projparm, spherecode, eastFile, northFile,
954  (int *)&errorcode, hinv_trans);
955 
956  for (int i = 0; i < npnts; i++)
957  {
958  /* Convert from meters (or any units that r_major and
959  * r_minor has) to lon/lat (radians) using GCTP */
960  /* --------------------------------------------------- */
961  errorcode =
962  hinv_trans[projcode] (
963  upleftpt[0],
964  upleftpt[1],
965  &lonrad, &latrad);
966 
967  }
968 #endif
969 
970  return 0;
971 
972 }
973 
974 
975 // convert angle degree to radian.
976 double
977 HE5_EHconvAng(double inAngle, int code)
978 {
979  long min = 0; /* Truncated Minutes */
980  long deg = 0; /* Truncated Degrees */
981 
982  double sec = 0.; /* Seconds */
983  double outAngle = 0.; /* Angle in desired units */
984  double pi = 3.14159265358979324;/* Pi */
985  double r2d = 180 / pi; /* Rad-deg conversion */
986  double d2r = 1 / r2d; /* Deg-rad conversion */
987 
988  switch (code)
989  {
990 
991  /* Convert radians to degrees */
992  /* -------------------------- */
993  case HE5_HDFE_RAD_DEG:
994  outAngle = inAngle * r2d;
995  break;
996 
997  /* Convert degrees to radians */
998  /* -------------------------- */
999  case HE5_HDFE_DEG_RAD:
1000  outAngle = inAngle * d2r;
1001  break;
1002 
1003 
1004  /* Convert packed degrees to degrees */
1005  /* --------------------------------- */
1006  case HE5_HDFE_DMS_DEG:
1007  deg = (long)(inAngle / 1000000);
1008  min = (long)((inAngle - deg * 1000000) / 1000);
1009  sec = (inAngle - deg * 1000000 - min * 1000);
1010  outAngle = deg + min / 60.0 + sec / 3600.0;
1011  break;
1012 
1013 
1014  /* Convert degrees to packed degrees */
1015  /* --------------------------------- */
1016  case HE5_HDFE_DEG_DMS:
1017  {
1018  deg = (long)inAngle;
1019  min = (long)((inAngle - deg) * 60);
1020  sec = (inAngle - deg - min / 60.0) * 3600;
1021 #if 0
1022 /*
1023  if ((int)sec == 60)
1024  {
1025  sec = sec - 60;
1026  min = min + 1;
1027  }
1028 */
1029 #endif
1030  if ( fabs(sec - 0.0) < 1.e-7 )
1031  {
1032  sec = 0.0;
1033  }
1034 
1035  if ( (fabs(sec - 60) < 1.e-7 ) || ( sec > 60.0 ))
1036  {
1037  sec = sec - 60;
1038  min = min + 1;
1039  if(sec < 0.0)
1040  {
1041  sec = 0.0;
1042  }
1043  }
1044  if (min == 60)
1045  {
1046  min = min - 60;
1047  deg = deg + 1;
1048  }
1049  outAngle = deg * 1000000 + min * 1000 + sec;
1050  }
1051  break;
1052 
1053 
1054  /* Convert radians to packed degrees */
1055  /* --------------------------------- */
1056  case HE5_HDFE_RAD_DMS:
1057  {
1058  inAngle = inAngle * r2d;
1059  deg = (long)inAngle;
1060  min = (long)((inAngle - deg) * 60);
1061  sec = ((inAngle - deg - min / 60.0) * 3600);
1062 #if 0
1063 /*
1064  if ((int)sec == 60)
1065  {
1066  sec = sec - 60;
1067  min = min + 1;
1068  }
1069 */
1070 #endif
1071  if ( fabs(sec - 0.0) < 1.e-7 )
1072  {
1073  sec = 0.0;
1074  }
1075 
1076  if ( (fabs(sec - 60) < 1.e-7 ) || ( sec > 60.0 ))
1077  {
1078  sec = sec - 60;
1079  min = min + 1;
1080  if(sec < 0.0)
1081  {
1082  sec = 0.0;
1083  }
1084  }
1085  if (min == 60)
1086  {
1087  min = min - 60;
1088  deg = deg + 1;
1089  }
1090  outAngle = deg * 1000000 + min * 1000 + sec;
1091  }
1092  break;
1093 
1094 
1095  /* Convert packed degrees to radians */
1096  /* --------------------------------- */
1097  case HE5_HDFE_DMS_RAD:
1098  deg = (long)(inAngle / 1000000);
1099  min = (long)((inAngle - deg * 1000000) / 1000);
1100  sec = (inAngle - deg * 1000000 - min * 1000);
1101  outAngle = deg + min / 60.0 + sec / 3600.0;
1102  outAngle = outAngle * d2r;
1103  break;
1104  default:
1105  break;
1106  }
1107  return (outAngle);
1108 }
1109 
1110 
1111 
1112 
1113 
1114 #if 0
1115 inline size_t
1117 HDF5CFUtil::INDEX_nD_TO_1D (const std::vector <size_t > &dims,
1118  const std::vector < size_t > &pos)
1119 {
1120  //
1121  // int a[10][20][30]; // & a[1][2][3] == a + (20*30+1 + 30*2 + 1 *3);
1122  // int b[10][2]; // &b[1][2] == b + (20*1 + 2);
1123  //
1124  if(dims.size () != pos.size ())
1125  throw InternalErr(__FILE__,__LINE__,"dimension error in INDEX_nD_TO_1D routine.");
1126  size_t sum = 0;
1127  size_t start = 1;
1128 
1129  for (size_t p = 0; p < pos.size (); p++) {
1130  size_t m = 1;
1131 
1132  for (size_t j = start; j < dims.size (); j++)
1133  m *= dims[j];
1134  sum += m * pos[p];
1135  start++;
1136  }
1137  return sum;
1138 }
1139 #endif
1140 
1142 //
1143 // \param input Input variable
1144 // \param dim dimension info of the input
1145 // \param start start indexes of each dim
1146 // \param stride stride of each dim
1147 // \param edge count of each dim
1148 // \param poutput output variable
1149 // \parrm index dimension index
1150 // \return 0 if successful. -1 otherwise.
1151 //
1152 //
1153 #if 0
1154 template<typename T>
1155 int HDF5CFUtil::subset(
1156  const T input[],
1157  int rank,
1158  vector<int> & dim,
1159  int start[],
1160  int stride[],
1161  int edge[],
1162  std::vector<T> *poutput,
1163  vector<int>& pos,
1164  int index)
1165 {
1166  for(int k=0; k<edge[index]; k++)
1167  {
1168  pos[index] = start[index] + k*stride[index];
1169  if(index+1<rank)
1170  subset(input, rank, dim, start, stride, edge, poutput,pos,index+1);
1171  if(index==rank-1)
1172  {
1173  poutput->push_back(input[INDEX_nD_TO_1D( dim, pos)]);
1174  }
1175  } // end of for
1176  return 0;
1177 } // end of template<typename T> static int
1178 #endif
1179 
1180 // Need to wrap a 'read buffer' from a pure file call here since read() is also a DAP function to read DAP data.
1181 ssize_t HDF5CFUtil::read_buffer_from_file(int fd, void*buf, size_t total_read) {
1182 
1183  ssize_t ret_val = read(fd,buf,total_read);
1184  return ret_val;
1185 }
1186 
1187 // Obtain the cache name. The clashing is rare given that fname is unique.The "_" may cause clashing in theory.
1188 string HDF5CFUtil::obtain_cache_fname(const string & fprefix, const string &fname, const string &vname) {
1189 
1190  string cache_fname = fprefix;
1191 
1192  string correct_fname = fname;
1193  std::replace(correct_fname.begin(),correct_fname.end(),'/','_');
1194 
1195  string correct_vname = vname;
1196 
1197  // Replace the '/' to '_'
1198  std::replace(correct_vname.begin(),correct_vname.end(),'/','_');
1199 
1200  // Replace the ' ' to to '_" since space is not good for a file name
1201  std::replace(correct_vname.begin(),correct_vname.end(),' ','_');
1202 
1203 
1204  cache_fname = cache_fname +correct_fname +correct_vname;
1205 
1206  return cache_fname;
1207 }
1208 size_t INDEX_nD_TO_1D (const std::vector < size_t > &dims,
1209  const std::vector < size_t > &pos){
1210  //
1211  // "int a[10][20][30] // & a[1][2][3] == a + (20*30+1 + 30*2 + 1 *3)"
1212  // "int b[10][2] // &b[1][2] == b + (20*1 + 2)"
1213  //
1214  if(dims.size () != pos.size ())
1215  throw InternalErr(__FILE__,__LINE__,"dimension error in INDEX_nD_TO_1D routine.");
1216  size_t sum = 0;
1217  size_t start = 1;
1218 
1219  for (size_t p = 0; p < pos.size (); p++) {
1220  size_t m = 1;
1221 
1222  for (size_t j = start; j < dims.size (); j++)
1223  m *= dims[j];
1224  sum += m * pos[p];
1225  start++;
1226  }
1227  return sum;
1228 }
1229 
1230 
1231 
HDF5CFUtil::H5type_to_H5DAPtype
static H5DataType H5type_to_H5DAPtype(hid_t h5_type_id)
Map HDF5 Datatype to the intermediate H5DAPtype for the future use.
Definition: HDF5CFUtil.cc:54
HDF5CFUtil.h
This file includes several helper functions for translating HDF5 to CF-compliant.
libdap
Definition: BESDapFunctionResponseCache.h:35
HDF5RequestHandler.h
include the entry functions to execute the handlers
HDF5CFUtil::read_buffer_from_file
static ssize_t read_buffer_from_file(int fd, void *buf, size_t)
Getting a subset of a variable.
Definition: HDF5CFUtil.cc:1181
HDF5CFUtil::Split
static void Split(const char *s, int len, char sep, std::vector< std::string > &names)
Definition: HDF5CFUtil.cc:329
HDF5CFUtil::trim_string
static std::string trim_string(hid_t dtypeid, const std::string s, int num_sect, size_t section_size, std::vector< size_t > &sect_newsize)
Definition: HDF5CFUtil.cc:223