bes  Updated for version 3.20.8
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 if ((H5INT64 == dtype) || (H5UINT64 == dtype)) {
190  if (HDF5RequestHandler::get_dmr_long_int()==false)
191  return false;
192  else
193  return true;
194  }
195  else
196  return true;
197 }
198 
199 bool HDF5CFUtil::cf_dap2_support_numeric_type(H5DataType dtype) {
200  if ((H5UNSUPTYPE == dtype)||(H5REFERENCE == dtype)
201  || (H5COMPOUND == dtype) || (H5ARRAY == dtype)
202  || (H5INT64 == dtype) ||(H5UINT64 == dtype)
203  || (H5FSTRING == dtype) ||(H5VSTRING == dtype))
204  return false;
205  else
206  return true;
207 }
208 
209 string HDF5CFUtil::obtain_string_after_lastslash(const string s) {
210 
211  string ret_str="";
212  size_t last_fslash_pos = s.find_last_of("/");
213  if (string::npos != last_fslash_pos &&
214  last_fslash_pos != (s.size()-1))
215  ret_str=s.substr(last_fslash_pos+1);
216  return ret_str;
217 }
218 
219 string HDF5CFUtil::obtain_string_before_lastslash(const string & s) {
220 
221  string ret_str="";
222  size_t last_fslash_pos = s.find_last_of("/");
223  if (string::npos != last_fslash_pos)
224  ret_str=s.substr(0,last_fslash_pos+1);
225  return ret_str;
226 
227 }
228 
229 string HDF5CFUtil::trim_string(hid_t ty_id,const string s, int num_sect, size_t sect_size, vector<size_t>& sect_newsize) {
230 
231  string temp_sect_str = "";
232  string temp_sect_newstr = "";
233  string final_str="";
234 
235  for (int i = 0; i < num_sect; i++) {
236 
237  if (i != (num_sect-1))
238  temp_sect_str = s.substr(i*sect_size,sect_size);
239  else
240  temp_sect_str = s.substr((num_sect-1)*sect_size,s.size()-(num_sect-1)*sect_size);
241 
242  size_t temp_pos = 0;
243 
244  if (H5T_STR_NULLTERM == H5Tget_strpad(ty_id))
245  temp_pos = temp_sect_str.find_first_of('\0');
246  else if (H5T_STR_SPACEPAD == H5Tget_strpad(ty_id))
247  temp_pos = temp_sect_str.find_last_not_of(' ')+1;
248  else temp_pos = temp_sect_str.find_last_not_of('0')+1;// NULL PAD
249 
250  if (temp_pos != string::npos) {
251 
252  // Here I only add a space at the end of each section for the SPACEPAD case,
253  // but not for NULL TERM
254  // or NULL PAD. Two reasons: C++ string doesn't need NULL TERM.
255  // We don't know and don't see any NULL PAD applications for NASA data.
256  if (H5T_STR_SPACEPAD == H5Tget_strpad(ty_id)) {
257 
258  if (temp_pos == temp_sect_str.size())
259  temp_sect_newstr = temp_sect_str +" ";
260  else
261  temp_sect_newstr = temp_sect_str.substr(0,temp_pos+1);
262 
263  sect_newsize[i] = temp_pos +1;
264  }
265  else {
266  temp_sect_newstr = temp_sect_str.substr(0,temp_pos);
267  sect_newsize[i] = temp_pos ;
268  }
269 
270  }
271 
272  else {// NULL is not found, adding a NULL at the end of this string.
273 
274  temp_sect_newstr = temp_sect_str;
275 
276  // Here I only add a space for the SPACEPAD, but not for NULL TERM
277  // or NULL PAD. Two reasons: C++ string doesn't need NULL TERM.
278  // We don't know and don't see any NULL PAD applications for NASA data.
279  if (H5T_STR_SPACEPAD == H5Tget_strpad(ty_id)) {
280  temp_sect_newstr.resize(temp_sect_str.size()+1);
281  temp_sect_newstr.append(1,' ');
282  sect_newsize[i] = sect_size + 1;
283  }
284  else
285  sect_newsize[i] = sect_size;
286  }
287  final_str+=temp_sect_newstr;
288  }
289 
290  return final_str;
291 }
292 
293 string HDF5CFUtil::remove_substrings(string str,const string &substr) {
294 
295  string::size_type i = str.find(substr);
296  while (i != std::string::npos) {
297  str.erase(i, substr.size());
298  i = str.find(substr, i);
299  }
300  return str;
301 }
302 // Obtain the unique name for the clashed names and save it to set namelist.
303 void HDF5CFUtil::gen_unique_name(string &str,set<string>& namelist, int&clash_index) {
304 
305  pair<set<string>::iterator,bool> ret;
306  string newstr = "";
307  stringstream sclash_index;
308  sclash_index << clash_index;
309  newstr = str + sclash_index.str();
310 
311  ret = namelist.insert(newstr);
312  if (false == ret.second) {
313  clash_index++;
314  gen_unique_name(str,namelist,clash_index);
315  }
316  else
317  str = newstr;
318 }
319 
320 void
321 HDF5CFUtil::Split_helper(vector<string> &tokens, const string &text, const char sep)
322 {
323  string::size_type start = 0, end = 0;
324  while ((end = text.find(sep, start)) != string::npos) {
325  tokens.push_back(text.substr(start, end - start));
326  start = end + 1;
327  }
328  tokens.push_back(text.substr(start));
329 }
330 
331 
332 // From a string separated by a separator to a list of string,
333 // for example, split "ab,c" to {"ab","c"}
334 void
335 HDF5CFUtil::Split(const char *s, int len, char sep, std::vector<std::string> &names)
336 {
337  names.clear();
338  for (int i = 0, j = 0; j <= len; ++j) {
339  if ((j == len && len) || s[j] == sep) {
340  string elem(s + i, j - i);
341  names.push_back(elem);
342  i = j + 1;
343  continue;
344  }
345  }
346 }
347 
348 
349 // Assume sz is Null terminated string.
350 void
351 HDF5CFUtil::Split(const char *sz, char sep, std::vector<std::string> &names)
352 {
353  Split(sz, (int)strlen(sz), sep, names);
354 }
355 
356 void HDF5CFUtil::parser_gpm_l3_gridheader(const vector<char>& value,
357  int& latsize, int&lonsize,
358  float& lat_start, float& lon_start,
359  float& lat_res, float& lon_res,
360  bool check_reg_orig ){
361 
362  float lat_north = 0.;
363  float lat_south = 0.;
364  float lon_east = 0.;
365  float lon_west = 0.;
366 
367  vector<string> ind_elems;
368  char sep='\n';
369  HDF5CFUtil::Split(&value[0],sep,ind_elems);
370 
371  // The number of elements in the GridHeader is 9. The string vector will add a leftover. So the size should be 10.
372  if(ind_elems.size()!=10)
373  throw InternalErr(__FILE__,__LINE__,"The number of elements in the GPM level 3 GridHeader is not right.");
374 
375  if(false == check_reg_orig) {
376  if (0 != ind_elems[1].find("Registration=CENTER"))
377  throw InternalErr(__FILE__,__LINE__,"The GPM grid registration is not center.");
378  }
379 
380  if (0 == ind_elems[2].find("LatitudeResolution")){
381 
382  size_t equal_pos = ind_elems[2].find_first_of('=');
383  if(string::npos == equal_pos)
384  throw InternalErr(__FILE__,__LINE__,"Cannot find latitude resolution for GPM level 3 products");
385 
386  size_t scolon_pos = ind_elems[2].find_first_of(';');
387  if(string::npos == scolon_pos)
388  throw InternalErr(__FILE__,__LINE__,"Cannot find latitude resolution for GPM level 3 products");
389  if (equal_pos < scolon_pos){
390 
391  string latres_str = ind_elems[2].substr(equal_pos+1,scolon_pos-equal_pos-1);
392  lat_res = strtof(latres_str.c_str(),NULL);
393  }
394  else
395  throw InternalErr(__FILE__,__LINE__,"latitude resolution is not right for GPM level 3 products");
396  }
397  else
398  throw InternalErr(__FILE__,__LINE__,"The GPM grid LatitudeResolution doesn't exist.");
399 
400  if (0 == ind_elems[3].find("LongitudeResolution")){
401 
402  size_t equal_pos = ind_elems[3].find_first_of('=');
403  if(string::npos == equal_pos)
404  throw InternalErr(__FILE__,__LINE__,"Cannot find longitude resolution for GPM level 3 products");
405 
406  size_t scolon_pos = ind_elems[3].find_first_of(';');
407  if(string::npos == scolon_pos)
408  throw InternalErr(__FILE__,__LINE__,"Cannot find longitude resolution for GPM level 3 products");
409  if (equal_pos < scolon_pos){
410  string lonres_str = ind_elems[3].substr(equal_pos+1,scolon_pos-equal_pos-1);
411  lon_res = strtof(lonres_str.c_str(),NULL);
412  }
413  else
414  throw InternalErr(__FILE__,__LINE__,"longitude resolution is not right for GPM level 3 products");
415  }
416  else
417  throw InternalErr(__FILE__,__LINE__,"The GPM grid LongitudeResolution doesn't exist.");
418 
419  if (0 == ind_elems[4].find("NorthBoundingCoordinate")){
420 
421  size_t equal_pos = ind_elems[4].find_first_of('=');
422  if(string::npos == equal_pos)
423  throw InternalErr(__FILE__,__LINE__,"Cannot find latitude resolution for GPM level 3 products");
424 
425  size_t scolon_pos = ind_elems[4].find_first_of(';');
426  if(string::npos == scolon_pos)
427  throw InternalErr(__FILE__,__LINE__,"Cannot find latitude resolution for GPM level 3 products");
428  if (equal_pos < scolon_pos){
429  string north_bounding_str = ind_elems[4].substr(equal_pos+1,scolon_pos-equal_pos-1);
430  lat_north = strtof(north_bounding_str.c_str(),NULL);
431  }
432  else
433  throw InternalErr(__FILE__,__LINE__,"NorthBoundingCoordinate is not right for GPM level 3 products");
434 
435  }
436  else
437  throw InternalErr(__FILE__,__LINE__,"The GPM grid NorthBoundingCoordinate doesn't exist.");
438 
439  if (0 == ind_elems[5].find("SouthBoundingCoordinate")){
440 
441  size_t equal_pos = ind_elems[5].find_first_of('=');
442  if(string::npos == equal_pos)
443  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for GPM level 3 products");
444 
445  size_t scolon_pos = ind_elems[5].find_first_of(';');
446  if(string::npos == scolon_pos)
447  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for GPM level 3 products");
448  if (equal_pos < scolon_pos){
449  string lat_south_str = ind_elems[5].substr(equal_pos+1,scolon_pos-equal_pos-1);
450  lat_south = strtof(lat_south_str.c_str(),NULL);
451  }
452  else
453  throw InternalErr(__FILE__,__LINE__,"south bound coordinate is not right for GPM level 3 products");
454 
455  }
456  else
457  throw InternalErr(__FILE__,__LINE__,"The GPM grid SouthBoundingCoordinate doesn't exist.");
458 
459  if (0 == ind_elems[6].find("EastBoundingCoordinate")){
460 
461  size_t equal_pos = ind_elems[6].find_first_of('=');
462  if(string::npos == equal_pos)
463  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for GPM level 3 products");
464 
465  size_t scolon_pos = ind_elems[6].find_first_of(';');
466  if(string::npos == scolon_pos)
467  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for GPM level 3 products");
468  if (equal_pos < scolon_pos){
469  string lon_east_str = ind_elems[6].substr(equal_pos+1,scolon_pos-equal_pos-1);
470  lon_east = strtof(lon_east_str.c_str(),NULL);
471  }
472  else
473  throw InternalErr(__FILE__,__LINE__,"south bound coordinate is not right for GPM level 3 products");
474 
475  }
476  else
477  throw InternalErr(__FILE__,__LINE__,"The GPM grid EastBoundingCoordinate doesn't exist.");
478 
479  if (0 == ind_elems[7].find("WestBoundingCoordinate")){
480 
481  size_t equal_pos = ind_elems[7].find_first_of('=');
482  if(string::npos == equal_pos)
483  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for GPM level 3 products");
484 
485  size_t scolon_pos = ind_elems[7].find_first_of(';');
486  if(string::npos == scolon_pos)
487  throw InternalErr(__FILE__,__LINE__,"Cannot find south bound coordinate for GPM level 3 products");
488  if (equal_pos < scolon_pos){
489  string lon_west_str = ind_elems[7].substr(equal_pos+1,scolon_pos-equal_pos-1);
490  lon_west = strtof(lon_west_str.c_str(),NULL);
491  }
492  else
493  throw InternalErr(__FILE__,__LINE__,"south bound coordinate is not right for GPM level 3 products");
494 
495  }
496  else
497  throw InternalErr(__FILE__,__LINE__,"The GPM grid WestBoundingCoordinate doesn't exist.");
498 
499  if (false == check_reg_orig) {
500  if (0 != ind_elems[8].find("Origin=SOUTHWEST"))
501  throw InternalErr(__FILE__,__LINE__,"The GPM grid origin is not SOUTHWEST.");
502  }
503 
504  // Since we only treat the case when the Registration is center, so the size should be the
505  // regular number size - 1.
506  latsize =(int)((lat_north-lat_south)/lat_res);
507  lonsize =(int)((lon_east-lon_west)/lon_res);
508  lat_start = lat_south;
509  lon_start = lon_west;
510 }
511 
512 void HDF5CFUtil::close_fileid(hid_t file_id,bool pass_fileid) {
513  if(false == pass_fileid) {
514  if(file_id != -1)
515  H5Fclose(file_id);
516  }
517 
518 }
519 
520 // Somehow the conversion of double to c++ string with sprintf causes the memory error in
521 // the testing code. I used the following code to do the conversion. Most part of the code
522 // in reverse, int_to_str and dtoa are adapted from geeksforgeeks.org
523 
524 // reverses a string 'str' of length 'len'
525 void HDF5CFUtil::rev_str(char *str, int len)
526 {
527  int i=0;
528  int j=len-1;
529  int temp = 0;
530  while (i<j)
531  {
532  temp = str[i];
533  str[i] = str[j];
534  str[j] = temp;
535  i++;
536  j--;
537  }
538 }
539 
540 // Converts a given integer x to string str[]. d is the number
541 // of digits required in output. If d is more than the number
542 // of digits in x, then 0s are added at the beginning.
543 int HDF5CFUtil::int_to_str(int x, char str[], int d)
544 {
545  int i = 0;
546  while (x)
547  {
548  str[i++] = (x%10) + '0';
549  x = x/10;
550  }
551 
552  // If number of digits required is more, then
553  // add 0s at the beginning
554  while (i < d)
555  str[i++] = '0';
556 
557  rev_str(str, i);
558  str[i] = '\0';
559  return i;
560 }
561 
562 // Converts a double floating point number to string.
563 void HDF5CFUtil::dtoa(double n, char *res, int afterpoint)
564 {
565  // Extract integer part
566  int ipart = (int)n;
567 
568  // Extract the double part
569  double fpart = n - (double)ipart;
570 
571  // convert integer part to string
572  int i = int_to_str(ipart, res, 0);
573 
574  // check for display option after point
575  if (afterpoint != 0)
576  {
577  res[i] = '.'; // add dot
578 
579  // Get the value of fraction part upto given no.
580  // of points after dot. The third parameter is needed
581  // to handle cases like 233.007
582  fpart = fpart * pow(10, afterpoint);
583 
584  // A round-error of 1 is found when casting to the integer for some numbers.
585  // We need to correct it.
586  int final_fpart = (int)fpart;
587  if(fpart -(int)fpart >0.5)
588  final_fpart = (int)fpart +1;
589  int_to_str(final_fpart, res + i + 1, afterpoint);
590  }
591 }
592 
593 
594 // Used to generate EOS geolocation cache name
595 string HDF5CFUtil::get_int_str(int x) {
596 
597  string str;
598  if(x > 0 && x <10)
599  str.push_back(x+'0');
600 
601  else if (x >10 && x<100) {
602  str.push_back(x/10+'0');
603  str.push_back(x%10+'0');
604  }
605  else {
606  int num_digit = 0;
607  int abs_x = (x<0)?-x:x;
608  while(abs_x/=10)
609  num_digit++;
610  if(x<=0)
611  num_digit++;
612  vector<char> buf;
613  buf.resize(num_digit);
614  snprintf(&buf[0],num_digit,"%d",x);
615  str.assign(&buf[0]);
616 
617  }
618 
619  return str;
620 
621 }
622 
623 //Used to generate EOS5 lat/lon cache name
624 string HDF5CFUtil::get_double_str(double x,int total_digit,int after_point) {
625 
626  string str;
627  if(x!=0) {
628  //char res[total_digit];
629  vector<char> res;
630  res.resize(total_digit);
631  for(int i = 0; i<total_digit;i++)
632  res[i] = '\0';
633  if (x<0) {
634  str.push_back('-');
635  dtoa(-x,&res[0],after_point);
636  for(int i = 0; i<total_digit;i++) {
637  if(res[i] != '\0')
638  str.push_back(res[i]);
639  }
640  }
641  else {
642  dtoa(x, &res[0], after_point);
643  //dtoa(x, res, after_point);
644  for(int i = 0; i<total_digit;i++) {
645  if(res[i] != '\0')
646  str.push_back(res[i]);
647  }
648  }
649 
650  }
651  else
652  str.push_back('0');
653 
654  return str;
655 
656 }
657 
658 
659 // This function is adapted from the HDF-EOS library.
660 int GDij2ll(int projcode, int zonecode, double projparm[],
661  int spherecode, int xdimsize, int ydimsize,
662  double upleftpt[], double lowrightpt[],
663  int npnts, int row[], int col[],
664  double longitude[], double latitude[], EOS5GridPRType pixcen, EOS5GridOriginType pixcnr)
665 {
666 
667  int errorcode = 0;
668  // In the original GCTP library, the function pointer names should be inv_trans and for_trans.
669  // However, since Hyrax supports both GDAL(including the HDF-EOS2 driver) and HDF handlers,
670  // on some machines, the functions inside the HDF-EOS2 driver will be called in run-time and wrong lat/lon
671  // values may be generated. To avoid, we change the function pointer names inside the GCTP library.
672  int(*hinv_trans[100]) (double,double,double*,double*);
673  int(*hfor_trans[100]) (double,double,double*,double*); /* GCTP function pointer */
674  double arg1, arg2;
675  double pixadjX = 0.; /* Pixel adjustment (x) */
676  double pixadjY = 0.; /* Pixel adjustment (y) */
677  double lonrad0 = 0.; /* Longitude in radians of upleft point */
678  double latrad0 = 0.; /* Latitude in radians of upleft point */
679  double scaleX = 0.; /* X scale factor */
680  double scaleY = 0.; /* Y scale factor */
681  double lonrad = 0.; /* Longitude in radians of point */
682  double latrad = 0.; /* Latitude in radians of point */
683  double xMtr0, yMtr0, xMtr1, yMtr1;
684 
685 
686 
687  /* Compute adjustment of position within pixel */
688  /* ------------------------------------------- */
689  if (pixcen == HE5_HDFE_CENTER)
690  {
691  /* Pixel defined at center */
692  /* ----------------------- */
693  pixadjX = 0.5;
694  pixadjY = 0.5;
695  }
696  else
697  {
698  switch (pixcnr)
699  {
700 
701  case HE5_HDFE_GD_UL:
702  {
703  /* Pixel defined at upper left corner */
704  /* ---------------------------------- */
705  pixadjX = 0.0;
706  pixadjY = 0.0;
707  break;
708  }
709 
710  case HE5_HDFE_GD_UR:
711  {
712  /* Pixel defined at upper right corner */
713  /* ----------------------------------- */
714  pixadjX = 1.0;
715  pixadjY = 0.0;
716  break;
717  }
718 
719  case HE5_HDFE_GD_LL:
720  {
721  /* Pixel defined at lower left corner */
722  /* ---------------------------------- */
723  pixadjX = 0.0;
724  pixadjY = 1.0;
725  break;
726  }
727 
728  case HE5_HDFE_GD_LR:
729  {
730 
731  /* Pixel defined at lower right corner */
732  /* ----------------------------------- */
733  pixadjX = 1.0;
734  pixadjY = 1.0;
735  break;
736  }
737 
738  default:
739  {
740  throw InternalErr(__FILE__,__LINE__,"Unknown pixel corner to retrieve lat/lon from a grid.");
741  }
742  }
743  }
744 
745 
746 
747  // If projection not GEO or BCEA call GCTP initialization routine
748  if (projcode != HE5_GCTP_GEO && projcode != HE5_GCTP_BCEA)
749  {
750 
751  scaleX = (lowrightpt[0] - upleftpt[0]) / xdimsize;
752  scaleY = (lowrightpt[1] - upleftpt[1]) / ydimsize;
753  string eastFile = HDF5RequestHandler::get_stp_east_filename();
754  string northFile = HDF5RequestHandler::get_stp_north_filename();
755 
756  hinv_init(projcode, zonecode, projparm, spherecode, (char*)eastFile.c_str(), (char*)northFile.c_str(),
757  &errorcode, hinv_trans);
758 
759 
760  /* Report error if any */
761  /* ------------------- */
762  if (errorcode != 0)
763  {
764  throw InternalErr(__FILE__,__LINE__,"GCTP hinv_init Error to retrieve lat/lon from a grid.");
765 
766  }
767  else
768  {
769  /* For each point ... */
770  /* ------------------ */
771  for (int i = 0; i < npnts; i++)
772  {
773  /* Convert from meters to lon/lat (radians) using GCTP */
774  /* --------------------------------------------------- */
775 #if 0
776  /*errorcode = hinv_trans[projcode] ((col[i] + pixadjX) * scaleX + upleftpt[0], (row[i] + pixadjY) * scaleY + upleftpt[1], &lonrad, &latrad);*/
777 #endif
778 
779  /* modified previous line to the following for the linux64 with -fPIC in cmpilation.
780  Whithout the change col[] and row[] values are ridiclous numbers, resulting a strange
781  number (very big) for arg1 and arg2. But with (int) typecast they become normal integers,
782  resulting in a acceptable values for arg1 and arg2. The problem was discovered during the
783  lat/lon geolocating of an hdfeos5 file with 64-bit hadview plug-in, developped for linux64.
784  */
785  arg1 = (((int)col[i] + pixadjX) * scaleX + upleftpt[0]);
786  arg2 = (((int)row[i] + pixadjY) * scaleY + upleftpt[1]);
787  errorcode = hinv_trans[projcode] (arg1, arg2, &lonrad, &latrad);
788 
789  /* Report error if any */
790  /* ------------------- */
791  if (errorcode != 0)
792  {
793 
794  if(projcode == HE5_GCTP_LAMAZ) {
795  longitude[i] = 1.0e51;
796  latitude[i] = 1.0e51;
797  }
798  else
799  throw InternalErr(__FILE__,__LINE__,"GCTP hinv_trans Error to retrieve lat/lon from a grid.");
800 
801  }
802  else
803  {
804 
805  /* Convert from radians to decimal degrees */
806  /* --------------------------------------- */
807  longitude[i] = HE5_EHconvAng(lonrad, HE5_HDFE_RAD_DEG);
808  latitude[i] = HE5_EHconvAng(latrad, HE5_HDFE_RAD_DEG);
809 
810  }
811  }
812  }
813  }
814  else if (projcode == HE5_GCTP_BCEA)
815  {
816  /* BCEA projection */
817  /* -------------- */
818 
819  /* Note: upleftpt and lowrightpt are in packed degrees, so they
820  must be converted to meters for this projection */
821 
822  /* Initialize forward transformation */
823  /* --------------------------------- */
824  hfor_init(projcode, zonecode, projparm, spherecode, NULL, NULL,&errorcode, hfor_trans);
825 
826  /* Report error if any */
827  /* ------------------- */
828  if (errorcode != 0)
829  {
830  throw InternalErr(__FILE__,__LINE__,"GCTP hfor_init Error to retrieve lat/lon from a grid.");
831  }
832 
833  /* Convert upleft and lowright X coords from DMS to radians */
834  /* -------------------------------------------------------- */
835  lonrad0 =HE5_EHconvAng(upleftpt[0], HE5_HDFE_DMS_RAD);
836  lonrad = HE5_EHconvAng(lowrightpt[0], HE5_HDFE_DMS_RAD);
837 
838  /* Convert upleft and lowright Y coords from DMS to radians */
839  /* -------------------------------------------------------- */
840  latrad0 = HE5_EHconvAng(upleftpt[1], HE5_HDFE_DMS_RAD);
841  latrad = HE5_EHconvAng(lowrightpt[1], HE5_HDFE_DMS_RAD);
842 
843  /* Convert form lon/lat to meters(or whatever unit is, i.e unit
844  of r_major and r_minor) using GCTP */
845  /* ----------------------------------------- */
846  errorcode = hfor_trans[projcode] (lonrad0, latrad0, &xMtr0, &yMtr0);
847 
848  /* Report error if any */
849  if (errorcode != 0)
850  {
851  throw InternalErr(__FILE__,__LINE__,"GCTP hfor_trans Error to retrieve lat/lon from a grid.");
852 
853  }
854 
855  /* Convert from lon/lat to meters or whatever unit is, i.e unit
856  of r_major and r_minor) using GCTP */
857  /* ----------------------------------------- */
858  errorcode = hfor_trans[projcode] (lonrad, latrad, &xMtr1, &yMtr1);
859 
860  /* Report error if any */
861  if (errorcode != 0)
862  {
863  throw InternalErr(__FILE__,__LINE__,"GCTP hfor_trans Error to retrieve lat/lon from a grid.");
864  }
865 
866  /* Compute x scale factor */
867  /* ---------------------- */
868  scaleX = (xMtr1 - xMtr0) / xdimsize;
869 
870  /* Compute y scale factor */
871  /* ---------------------- */
872  scaleY = (yMtr1 - yMtr0) / ydimsize;
873 
874  /* Initialize inverse transformation */
875  /* --------------------------------- */
876  hinv_init(projcode, zonecode, projparm, spherecode, NULL, NULL, &errorcode, hinv_trans);
877  /* Report error if any */
878  /* ------------------- */
879  if (errorcode != 0)
880  {
881  throw InternalErr(__FILE__,__LINE__,"GCTP hinv_init Error to retrieve lat/lon from a grid.");
882  }
883  /* For each point ... */
884  /* ------------------ */
885  for (int i = 0; i < npnts; i++)
886  {
887  /* Convert from meters (or any units that r_major and
888  r_minor has) to lon/lat (radians) using GCTP */
889  /* --------------------------------------------------- */
890  errorcode = hinv_trans[projcode] (
891  (col[i] + pixadjX) * scaleX + xMtr0,
892  (row[i] + pixadjY) * scaleY + yMtr0,
893  &lonrad, &latrad);
894 
895  /* Report error if any */
896  /* ------------------- */
897  if (errorcode != 0)
898  {
899 #if 0
900  /* status = -1;
901  sprintf(errbuf, "GCTP Error: %li\n", errorcode);
902  H5Epush(__FILE__, "HE5_GDij2ll", __LINE__, H5E_ARGS, H5E_BADVALUE , errbuf);
903  HE5_EHprint(errbuf, __FILE__, __LINE__);
904  return (status); */
905 #endif
906  longitude[i] = 1.0e51; /* PGSd_GCT_IN_ERROR */
907  latitude[i] = 1.0e51; /* PGSd_GCT_IN_ERROR */
908  }
909 
910  /* Convert from radians to decimal degrees */
911  /* --------------------------------------- */
912  longitude[i] = HE5_EHconvAng(lonrad, HE5_HDFE_RAD_DEG);
913  latitude[i] = HE5_EHconvAng(latrad, HE5_HDFE_RAD_DEG);
914  }
915  }
916 
917  else if (projcode == HE5_GCTP_GEO)
918  {
919  /* GEO projection */
920  /* -------------- */
921 
922  /*
923  * Note: lonrad, lonrad0, latrad, latrad0 are actually in degrees for
924  * the GEO projection case.
925  */
926 
927 
928  /* Convert upleft and lowright X coords from DMS to degrees */
929  /* -------------------------------------------------------- */
930  lonrad0 = HE5_EHconvAng(upleftpt[0], HE5_HDFE_DMS_DEG);
931  lonrad = HE5_EHconvAng(lowrightpt[0], HE5_HDFE_DMS_DEG);
932 
933  /* Compute x scale factor */
934  /* ---------------------- */
935  scaleX = (lonrad - lonrad0) / xdimsize;
936 
937  /* Convert upleft and lowright Y coords from DMS to degrees */
938  /* -------------------------------------------------------- */
939  latrad0 = HE5_EHconvAng(upleftpt[1], HE5_HDFE_DMS_DEG);
940  latrad = HE5_EHconvAng(lowrightpt[1], HE5_HDFE_DMS_DEG);
941 
942  /* Compute y scale factor */
943  /* ---------------------- */
944  scaleY = (latrad - latrad0) / ydimsize;
945 
946  /* For each point ... */
947  /* ------------------ */
948  for (int i = 0; i < npnts; i++)
949  {
950  /* Convert to lon/lat (decimal degrees) */
951  /* ------------------------------------ */
952  longitude[i] = (col[i] + pixadjX) * scaleX + lonrad0;
953  latitude[i] = (row[i] + pixadjY) * scaleY + latrad0;
954  }
955  }
956 
957 
958 #if 0
959  hinv_init(projcode, zonecode, projparm, spherecode, eastFile, northFile,
960  (int *)&errorcode, hinv_trans);
961 
962  for (int i = 0; i < npnts; i++)
963  {
964  /* Convert from meters (or any units that r_major and
965  * r_minor has) to lon/lat (radians) using GCTP */
966  /* --------------------------------------------------- */
967  errorcode =
968  hinv_trans[projcode] (
969  upleftpt[0],
970  upleftpt[1],
971  &lonrad, &latrad);
972 
973  }
974 #endif
975 
976  return 0;
977 
978 }
979 
980 
981 // convert angle degree to radian.
982 double
983 HE5_EHconvAng(double inAngle, int code)
984 {
985  long min = 0; /* Truncated Minutes */
986  long deg = 0; /* Truncated Degrees */
987 
988  double sec = 0.; /* Seconds */
989  double outAngle = 0.; /* Angle in desired units */
990  double pi = 3.14159265358979324;/* Pi */
991  double r2d = 180 / pi; /* Rad-deg conversion */
992  double d2r = 1 / r2d; /* Deg-rad conversion */
993 
994  switch (code)
995  {
996 
997  /* Convert radians to degrees */
998  /* -------------------------- */
999  case HE5_HDFE_RAD_DEG:
1000  outAngle = inAngle * r2d;
1001  break;
1002 
1003  /* Convert degrees to radians */
1004  /* -------------------------- */
1005  case HE5_HDFE_DEG_RAD:
1006  outAngle = inAngle * d2r;
1007  break;
1008 
1009 
1010  /* Convert packed degrees to degrees */
1011  /* --------------------------------- */
1012  case HE5_HDFE_DMS_DEG:
1013  deg = (long)(inAngle / 1000000);
1014  min = (long)((inAngle - deg * 1000000) / 1000);
1015  sec = (inAngle - deg * 1000000 - min * 1000);
1016  outAngle = deg + min / 60.0 + sec / 3600.0;
1017  break;
1018 
1019 
1020  /* Convert degrees to packed degrees */
1021  /* --------------------------------- */
1022  case HE5_HDFE_DEG_DMS:
1023  {
1024  deg = (long)inAngle;
1025  min = (long)((inAngle - deg) * 60);
1026  sec = (inAngle - deg - min / 60.0) * 3600;
1027 #if 0
1028 /*
1029  if ((int)sec == 60)
1030  {
1031  sec = sec - 60;
1032  min = min + 1;
1033  }
1034 */
1035 #endif
1036  if ( fabs(sec - 0.0) < 1.e-7 )
1037  {
1038  sec = 0.0;
1039  }
1040 
1041  if ( (fabs(sec - 60) < 1.e-7 ) || ( sec > 60.0 ))
1042  {
1043  sec = sec - 60;
1044  min = min + 1;
1045  if(sec < 0.0)
1046  {
1047  sec = 0.0;
1048  }
1049  }
1050  if (min == 60)
1051  {
1052  min = min - 60;
1053  deg = deg + 1;
1054  }
1055  outAngle = deg * 1000000 + min * 1000 + sec;
1056  }
1057  break;
1058 
1059 
1060  /* Convert radians to packed degrees */
1061  /* --------------------------------- */
1062  case HE5_HDFE_RAD_DMS:
1063  {
1064  inAngle = inAngle * r2d;
1065  deg = (long)inAngle;
1066  min = (long)((inAngle - deg) * 60);
1067  sec = ((inAngle - deg - min / 60.0) * 3600);
1068 #if 0
1069 /*
1070  if ((int)sec == 60)
1071  {
1072  sec = sec - 60;
1073  min = min + 1;
1074  }
1075 */
1076 #endif
1077  if ( fabs(sec - 0.0) < 1.e-7 )
1078  {
1079  sec = 0.0;
1080  }
1081 
1082  if ( (fabs(sec - 60) < 1.e-7 ) || ( sec > 60.0 ))
1083  {
1084  sec = sec - 60;
1085  min = min + 1;
1086  if(sec < 0.0)
1087  {
1088  sec = 0.0;
1089  }
1090  }
1091  if (min == 60)
1092  {
1093  min = min - 60;
1094  deg = deg + 1;
1095  }
1096  outAngle = deg * 1000000 + min * 1000 + sec;
1097  }
1098  break;
1099 
1100 
1101  /* Convert packed degrees to radians */
1102  /* --------------------------------- */
1103  case HE5_HDFE_DMS_RAD:
1104  deg = (long)(inAngle / 1000000);
1105  min = (long)((inAngle - deg * 1000000) / 1000);
1106  sec = (inAngle - deg * 1000000 - min * 1000);
1107  outAngle = deg + min / 60.0 + sec / 3600.0;
1108  outAngle = outAngle * d2r;
1109  break;
1110  default:
1111  break;
1112  }
1113  return (outAngle);
1114 }
1115 
1116 
1117 
1118 
1119 
1120 #if 0
1122 inline size_t
1123 HDF5CFUtil::INDEX_nD_TO_1D (const std::vector <size_t > &dims,
1124  const std::vector < size_t > &pos)
1125 {
1126  //
1127  // int a[10][20][30]; // & a[1][2][3] == a + (20*30+1 + 30*2 + 1 *3);
1128  // int b[10][2]; // &b[1][2] == b + (20*1 + 2);
1129  //
1130  if(dims.size () != pos.size ())
1131  throw InternalErr(__FILE__,__LINE__,"dimension error in INDEX_nD_TO_1D routine.");
1132  size_t sum = 0;
1133  size_t start = 1;
1134 
1135  for (size_t p = 0; p < pos.size (); p++) {
1136  size_t m = 1;
1137 
1138  for (size_t j = start; j < dims.size (); j++)
1139  m *= dims[j];
1140  sum += m * pos[p];
1141  start++;
1142  }
1143  return sum;
1144 }
1145 #endif
1146 
1148 //
1149 // \param input Input variable
1150 // \param dim dimension info of the input
1151 // \param start start indexes of each dim
1152 // \param stride stride of each dim
1153 // \param edge count of each dim
1154 // \param poutput output variable
1155 // \parrm index dimension index
1156 // \return 0 if successful. -1 otherwise.
1157 //
1158 //
1159 #if 0
1160 template<typename T>
1161 int HDF5CFUtil::subset(
1162  const T input[],
1163  int rank,
1164  vector<int> & dim,
1165  int start[],
1166  int stride[],
1167  int edge[],
1168  std::vector<T> *poutput,
1169  vector<int>& pos,
1170  int index)
1171 {
1172  for(int k=0; k<edge[index]; k++)
1173  {
1174  pos[index] = start[index] + k*stride[index];
1175  if(index+1<rank)
1176  subset(input, rank, dim, start, stride, edge, poutput,pos,index+1);
1177  if(index==rank-1)
1178  {
1179  poutput->push_back(input[INDEX_nD_TO_1D( dim, pos)]);
1180  }
1181  } // end of for
1182  return 0;
1183 } // end of template<typename T> static int
1184 #endif
1185 
1186 // Need to wrap a 'read buffer' from a pure file call here since read() is also a DAP function to read DAP data.
1187 ssize_t HDF5CFUtil::read_buffer_from_file(int fd, void*buf, size_t total_read) {
1188 
1189  ssize_t ret_val = read(fd,buf,total_read);
1190  return ret_val;
1191 }
1192 
1193 // Obtain the cache name. The clashing is rare given that fname is unique.The "_" may cause clashing in theory.
1194 string HDF5CFUtil::obtain_cache_fname(const string & fprefix, const string &fname, const string &vname) {
1195 
1196  string cache_fname = fprefix;
1197 
1198  string correct_fname = fname;
1199  std::replace(correct_fname.begin(),correct_fname.end(),'/','_');
1200 
1201  string correct_vname = vname;
1202 
1203  // Replace the '/' to '_'
1204  std::replace(correct_vname.begin(),correct_vname.end(),'/','_');
1205 
1206  // Replace the ' ' to to '_" since space is not good for a file name
1207  std::replace(correct_vname.begin(),correct_vname.end(),' ','_');
1208 
1209 
1210  cache_fname = cache_fname +correct_fname +correct_vname;
1211 
1212  return cache_fname;
1213 }
1214 size_t INDEX_nD_TO_1D (const std::vector < size_t > &dims,
1215  const std::vector < size_t > &pos){
1216  //
1217  // "int a[10][20][30] // & a[1][2][3] == a + (20*30+1 + 30*2 + 1 *3)"
1218  // "int b[10][2] // &b[1][2] == b + (20*1 + 2)"
1219  //
1220  if(dims.size () != pos.size ())
1221  throw InternalErr(__FILE__,__LINE__,"dimension error in INDEX_nD_TO_1D routine.");
1222  size_t sum = 0;
1223  size_t start = 1;
1224 
1225  for (size_t p = 0; p < pos.size (); p++) {
1226  size_t m = 1;
1227 
1228  for (size_t j = start; j < dims.size (); j++)
1229  m *= dims[j];
1230  sum += m * pos[p];
1231  start++;
1232  }
1233  return sum;
1234 }
1235 
1236 void HDF5CFUtil::get_relpath_pos(const string& temp_str, const string& relpath, vector<size_t>&s_pos) {
1237 
1238 
1239  //vector<size_t> positions; // holds all the positions that sub occurs within str
1240 
1241  size_t pos = temp_str.find(relpath, 0);
1242  while(pos != string::npos)
1243  {
1244  s_pos.push_back(pos);
1245 //cout<<"pos is "<<pos <<endl;
1246  pos = temp_str.find(relpath,pos+1);
1247  }
1248 //cout<<"pos.size() is "<<s_pos.size() <<endl;
1249 
1250 
1251 }
1252 
1253 void HDF5CFUtil::cha_co(string &co,const string & vpath) {
1254 
1255  string sep="/";
1256  string rp_sep="../";
1257  if(vpath.find(sep,1)!=string::npos) {
1258  // if finding '/' in the co;
1259  if(co.find(sep)!=string::npos) {
1260  // if finding '../', reduce the path.
1261  if(co.find(rp_sep)!=string::npos) {
1262  vector<size_t>var_sep_pos;
1263  get_relpath_pos(vpath,sep,var_sep_pos);
1264  vector<size_t>co_rp_sep_pos;
1265  get_relpath_pos(co,rp_sep,co_rp_sep_pos);
1266  if(co_rp_sep_pos[0]==0) {
1267  // Obtain the '../' position at co
1268  if(co_rp_sep_pos.size() <var_sep_pos.size()) {
1269  size_t var_prefix_pos=var_sep_pos[var_sep_pos.size()-co_rp_sep_pos.size()-1];
1270 //cout<<"var_prefix_pos is "<<var_prefix_pos <<endl;
1271  string var_prefix=vpath.substr(1,var_prefix_pos);
1272 //cout<<"var_prefix is "<<var_prefix <<endl;
1273  string co_suffix = co.substr(co_rp_sep_pos[co_rp_sep_pos.size()-1]+rp_sep.size());
1274 //cout<<"co_suffix is "<<co_suffix <<endl;
1275  co = var_prefix+co_suffix;
1276  }
1277 //cout<<"co is "<<co<<endl;;
1278  }
1279  }
1280 
1281  }
1282  else {// co no path, add fullpath
1283  string var_prefix = obtain_string_before_lastslash(vpath).substr(1);
1284  co = var_prefix +co;
1285 //cout<<"co full is " <<co <<endl;
1286 
1287 
1288  }
1289  }
1290 
1291 }
1292 
This file includes several helper functions for translating HDF5 to CF-compliant.
include the entry functions to execute the handlers
static void Split(const char *s, int len, char sep, std::vector< std::string > &names)
Definition: HDF5CFUtil.cc:335
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
static ssize_t read_buffer_from_file(int fd, void *buf, size_t)
Getting a subset of a variable.
Definition: HDF5CFUtil.cc:1187
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:229