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
43using namespace std;
44using namespace libdap;
45// For using GCTP to calculate the lat/lon
46extern "C" {
47int hinv_init(int insys, int inzone, double *inparm, int indatum, char *fn27, char *fn83, int *iflg, int (*hinv_trans[])(double, double, double*, double*));
48
49int 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
53H5DataType
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
119size_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
143bool 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.
164bool 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
182bool 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
199bool 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
215string 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
225string 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
235string 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
299string 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.
309void 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
326void
327HDF5CFUtil::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"}
340void
341HDF5CFUtil::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.
356void
357HDF5CFUtil::Split(const char *sz, char sep, std::vector<std::string> &names)
358{
359 Split(sz, (int)strlen(sz), sep, names);
360}
361
362void 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
520void 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'
533void 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.
551int 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.
571void 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
603string 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
632string 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.
668int 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.
990double
991HE5_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
1130inline size_t
1131HDF5CFUtil::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
1168template<typename T>
1169int 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.
1195ssize_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.
1202string 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}
1222size_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
1244void 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
1261void 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