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