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