bes Updated for version 3.20.13
HDFEOS5CFMissLLArray.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 "config_hdf5.h"
33#include <sys/stat.h>
34#include <iostream>
35#include <sstream>
36#include <cassert>
37#include <BESDebug.h>
38#include <libdap/InternalErr.h>
39//#include <libdap/Array.h>
40
42#include "HDF5RequestHandler.h"
43
44using namespace std;
45using namespace libdap;
46
47BaseType *HDFEOS5CFMissLLArray::ptr_duplicate()
48{
49 return new HDFEOS5CFMissLLArray(*this);
50}
51
52bool HDFEOS5CFMissLLArray::read()
53{
54
55 BESDEBUG("h5","Coming to HDFEOS5CFMissLLArray read "<<endl);
56 if(nullptr == HDF5RequestHandler::get_lrdata_mem_cache())
57 read_data_NOT_from_mem_cache(false,nullptr);
58 else {
59 vector<string> cur_lrd_non_cache_dir_list;
60 HDF5RequestHandler::get_lrd_non_cache_dir_list(cur_lrd_non_cache_dir_list);
61
62 string cache_key;
63 // Check if this file is included in the non-cache directory
64 if( (cur_lrd_non_cache_dir_list.empty()) ||
65 ("" == check_str_sect_in_list(cur_lrd_non_cache_dir_list,filename,'/'))) {
66 short cache_flag = 2;
67 vector<string> cur_cache_dlist;
68 HDF5RequestHandler::get_lrd_cache_dir_list(cur_cache_dlist);
69 string cache_dir = check_str_sect_in_list(cur_cache_dlist,filename,'/');
70 if(cache_dir != ""){
71 cache_key = cache_dir + varname;
72 cache_flag = 3;
73 }
74 else
75 cache_key = filename + varname;
76
77 // Need to obtain the total number of elements.
78 // Currently only trivial geographic projection is supported.
79 // So the total number of elements for LAT is ydimsize,
80 // the total number of elements for LON is xdimsize.
81 if(cvartype == CV_LAT_MISS)
82 handle_data_with_mem_cache(H5FLOAT32,(size_t)ydimsize,cache_flag,cache_key,false);
83 else
84 handle_data_with_mem_cache(H5FLOAT32,(size_t)xdimsize,cache_flag,cache_key,false);
85 }
86 else
87 read_data_NOT_from_mem_cache(false,nullptr);
88 }
89 return true;
90}
91
92void HDFEOS5CFMissLLArray::read_data_NOT_from_mem_cache(bool add_cache,void*buf){
93
94 BESDEBUG("h5","Coming to read_data_NOT_from_mem_cache "<<endl);
95
96 // First handle geographic projection. No GCTP is needed.Since the calculation
97 // of lat/lon is really simple for this case. No need to provide the disk cache
98 // unless the calculation takes too long. We will see.
99 if(eos5_projcode == HE5_GCTP_GEO) {
100 read_data_NOT_from_mem_cache_geo(add_cache,buf);
101 return;
102 }
103
104 int nelms = -1;
105 vector<int>offset;
106 vector<int>count;
107 vector<int>step;
108
109 if (rank <= 0)
110 throw InternalErr (__FILE__, __LINE__,
111 "The number of dimension of this variable should be greater than 0");
112 else {
113 offset.resize(rank);
114 count.resize(rank);
115 step.resize(rank);
116 nelms = format_constraint (offset.data(), step.data(), count.data());
117 }
118
119 if (nelms <= 0)
120 throw InternalErr (__FILE__, __LINE__,
121 "The number of elments is negative.");
122
123 vector<size_t>pos(rank,0);
124 for (int i = 0; i< rank; i++)
125 pos[i] = offset[i];
126
127 vector<size_t>dimsizes;
128 dimsizes.push_back(ydimsize);
129 dimsizes.push_back(xdimsize);
130 int total_elms = xdimsize*ydimsize;
131
132 double upleft[2];
133 double lowright[2];
134 vector<int>rows;
135 vector<int>cols;
136 vector<double>lon;
137 vector<double>lat;
138 rows.resize(xdimsize*ydimsize);
139 cols.resize(xdimsize*ydimsize);
140 lon.resize(xdimsize*ydimsize);
141 lat.resize(xdimsize*ydimsize);
142
143 upleft[0] = point_left;
144 upleft[1] = point_upper;
145 lowright[0] = point_right;
146 lowright[1] = point_lower;
147
148
149 int j = 0;
150 int r = -1;
151
152 for (int k = j = 0; j < ydimsize; ++j) {
153 for (int i = 0; i < xdimsize; ++i) {
154 rows[k] = j;
155 cols[k] = i;
156 ++k;
157 }
158 }
159
160 BESDEBUG("h5", " Before calling GDij2ll, check all projection parameters. " << endl);
161 BESDEBUG("h5", " eos5_projcode is " << eos5_projcode <<endl);
162 BESDEBUG("h5", " eos5_zone is " << eos5_zone <<endl);
163 BESDEBUG("h5", " eos5_params[0] is " << eos5_params[0] <<endl);
164 BESDEBUG("h5", " eos5_params[1] is " << eos5_params[1] <<endl);
165 BESDEBUG("h5", " eos5_sphere is " << eos5_sphere <<endl);
166 BESDEBUG("h5", " xdimsize is " << xdimsize <<endl);
167 BESDEBUG("h5", " ydimsize is " << ydimsize <<endl);
168 BESDEBUG("h5", " eos5_pixelreg is " << eos5_pixelreg <<endl);
169 BESDEBUG("h5", " eos5_origin is " << eos5_origin <<endl);
170 BESDEBUG("h5", " upleft[0] is " << upleft[0] <<endl);
171 BESDEBUG("h5", " upleft[1] is " << upleft[1] <<endl);
172 BESDEBUG("h5", " lowright[0] is " << lowright[0] <<endl);
173 BESDEBUG("h5", " lowright[1] is " << lowright[1] <<endl);
174
175#if 0
176 cerr<< " eos5_params[0] is " << eos5_params[0] <<endl;
177 cerr<< " eos5_params[1] is " << eos5_params[1] <<endl;
178 cerr<< " eos5_sphere is " << eos5_sphere <<endl;
179 cerr<< " eos5_zone is " << eos5_zone <<endl;
180 cerr<< " Before calling GDij2ll, check all projection parameters. " << endl;
181 cerr<< " eos5_zone is " << eos5_zone <<endl;
182 cerr<< " eos5_params[0] is " << eos5_params[0] <<endl;
183 cerr<< " eos5_params[1] is " << eos5_params[1] <<endl;
184 BESDEBUG("h5", " xdimsize is " << xdimsize <<endl);
185 BESDEBUG("h5", " ydimsize is " << ydimsize <<endl);
186 BESDEBUG("h5", " eos5_pixelreg is " << eos5_pixelreg <<endl);
187 BESDEBUG("h5", " eos5_origin is " << eos5_origin <<endl);
188 BESDEBUG("h5", " upleft[0] is " << upleft[0] <<endl);
189 BESDEBUG("h5", " upleft[1] is " << upleft[1] <<endl);
190 BESDEBUG("h5", " lowright[0] is " << lowright[0] <<endl);
191 BESDEBUG("h5", " lowright[1] is " << lowright[1] <<endl);
192#endif
193
194 // boolean to mark if the value can be read from the cache.
195 // Not necessary from programming point of view. Use it to
196 // make the code flow clear.
197 bool ll_read_from_cache = true;
198
199 // Check if geo-cache is turned on.
200 bool use_latlon_cache = HDF5RequestHandler::get_use_eosgeo_cachefile();
201
202
203 // Adding more code to read latlon from cache here.
204 if(use_latlon_cache == true) {
205
206 // Cache name is made in a special way to make sure the same lat/lon
207 // fall into the same cache file.
208 string cache_fname = obtain_ll_cache_name();
209
210 // Obtain eos-geo disk cache size, dir and prefix.
211 long ll_disk_cache_size = HDF5RequestHandler::get_latlon_disk_cache_size();
212 string ll_disk_cache_dir = HDF5RequestHandler::get_latlon_disk_cache_dir();
213 string ll_disk_cache_prefix = HDF5RequestHandler::get_latlon_disk_cachefile_prefix();
214
215 // Expected cache file size
216 int expected_file_size = 2*xdimsize*ydimsize*sizeof(double);
217 HDF5DiskCache *ll_cache = HDF5DiskCache::get_instance(ll_disk_cache_size,ll_disk_cache_dir,ll_disk_cache_prefix);
218 int fd = 0;
219 ll_read_from_cache = ll_cache->get_data_from_cache(cache_fname, expected_file_size,fd);
220
221 if(ll_read_from_cache == true) {
222
223 BESDEBUG("h5", " Read latitude and longitude from a disk cache. " <<endl);
224 size_t var_offset = 0;
225 // Longitude is stored after latitude.
226 if(CV_LON_MISS == cvartype)
227 var_offset = xdimsize*ydimsize*sizeof(double);
228
229 vector<double> var_value;
230 var_value.resize(xdimsize*ydimsize);
231
232#if 0
233 //Latitude starts from 0, longitude starts from xdimsize*ydimsize*sizeof(double);
234#endif
235 off_t fpos = lseek(fd,var_offset,SEEK_SET);
236 if(fpos == -1) {
237 throw InternalErr (__FILE__, __LINE__,
238 "Cannot seek the cached file offset.");
239
240 }
241 ssize_t ret_val = HDF5CFUtil::read_buffer_from_file(fd,(void*)var_value.data(),var_value.size()*sizeof(double));
242 ll_cache->unlock_and_close(cache_fname);
243
244 // If the cached file is not read correctly, we should purge the file.
245 if((-1 == ret_val) || ((size_t)ret_val != (xdimsize*ydimsize*sizeof(double)))) {
246 ll_cache->purge_file(cache_fname);
247 ll_read_from_cache = false;
248 }
249 else {
250 // short-cut, no need to do subset.
251 if(total_elms == nelms)
252 set_value((dods_float64 *)var_value.data(),total_elms);
253 else {
254 vector<double>val;
255 subset<double>(
256 var_value.data(),
257 rank,
258 dimsizes,
259 offset.data(),
260 step.data(),
261 count.data(),
262 &val,
263 pos,
264 0);
265 set_value((dods_float64 *)val.data(),nelms);
266 }
267 return;
268 }
269 }
270 else
271 ll_read_from_cache = false;
272 }
273
274 // Calculate Lat/lon by using GCTP
275 r = GDij2ll (eos5_projcode, eos5_zone, eos5_params.data(), eos5_sphere, xdimsize, ydimsize, upleft, lowright,
276 xdimsize * ydimsize, rows.data(), cols.data(), lon.data(), lat.data(), eos5_pixelreg, eos5_origin);
277 if (r != 0) {
278 ostringstream eherr;
279 eherr << "cannot calculate grid latitude and longitude";
280 throw InternalErr (__FILE__, __LINE__, eherr.str ());
281 }
282
283
284 // ll_read_from_cache may be redundant. It is good to remind that we will generate the cache file
285 // when reading from the cache fails.
286 // We are using cache and need to write data to the cache.
287 if(use_latlon_cache == true && ll_read_from_cache == false) {
288 string cache_fname = obtain_ll_cache_name();
289 long ll_disk_cache_size = HDF5RequestHandler::get_latlon_disk_cache_size();
290 string ll_disk_cache_dir = HDF5RequestHandler::get_latlon_disk_cache_dir();
291 string ll_disk_cache_prefix = HDF5RequestHandler::get_latlon_disk_cachefile_prefix();
292
293 BESDEBUG("h5", " Write EOS5 grid latitude and longitude to a disk cache. " <<endl);
294
295 HDF5DiskCache *ll_cache = HDF5DiskCache::get_instance(ll_disk_cache_size,ll_disk_cache_dir,ll_disk_cache_prefix);
296
297 // Merge vector lat and lon. lat first.
298 vector <double>latlon;
299 latlon.reserve(xdimsize*ydimsize*2);
300 latlon.insert(latlon.end(),lat.begin(),lat.end());
301 latlon.insert(latlon.end(),lon.begin(),lon.end());
302 ll_cache->write_cached_data(cache_fname,2*xdimsize*ydimsize*sizeof(double),latlon);
303
304 }
305 BESDEBUG("h5", " The first value of lon is " << lon[0] <<endl);
306 BESDEBUG("h5", " The first value of lat is " << lat[0] <<endl);
307
308#if 0
309 vector<size_t>pos(rank,0);
310 for (int i = 0; i< rank; i++)
311 pos[i] = offset[i];
312
313 vector<size_t>dimsizes;
314 dimsizes.push_back(ydimsize);
315 dimsizes.push_back(xdimsize);
316 int total_elms = xdimsize*ydimsize;
317#endif
318
319
320 if(CV_LON_MISS == cvartype) {
321 if(total_elms == nelms)
322 set_value((dods_float64 *)lon.data(),total_elms);
323 else {
324 vector<double>val;
325 subset<double>(
326 lon.data(),
327 rank,
328 dimsizes,
329 offset.data(),
330 step.data(),
331 count.data(),
332 &val,
333 pos,
334 0);
335 set_value((dods_float64 *)val.data(),nelms);
336 }
337
338 }
339 else if(CV_LAT_MISS == cvartype) {
340
341 if(total_elms == nelms)
342 set_value((dods_float64 *)lat.data(),total_elms);
343 else {
344 vector<double>val;
345 subset<double>(
346 lat.data(),
347 rank,
348 dimsizes,
349 offset.data(),
350 step.data(),
351 count.data(),
352 &val,
353 pos,
354 0);
355 set_value((dods_float64 *)val.data(),nelms);
356 }
357 }
358 return;
359
360}
361
362string HDFEOS5CFMissLLArray::obtain_ll_cache_name() {
363
364 BESDEBUG("h5","Coming to obtain_ll_cache_name "<<endl);
365
366
367 // Here we have a sanity check for the cached parameters:Cached directory,file prefix and cached directory size.
368 // Supposedly Hyrax BES cache feature should check this and the code exists. However, the
369 string bescachedir = HDF5RequestHandler::get_latlon_disk_cache_dir();
370 string bescacheprefix = HDF5RequestHandler::get_latlon_disk_cachefile_prefix();
371 long cachesize = HDF5RequestHandler::get_latlon_disk_cache_size();
372
373 if(("" == bescachedir)||(""==bescacheprefix)||(cachesize <=0)){
374 throw InternalErr (__FILE__, __LINE__, "Either the cached dir is empty or the prefix is nullptr or the cache size is not set.");
375 }
376 else {
377 struct stat sb;
378 if(stat(bescachedir.c_str(),&sb) !=0) {
379 string err_mesg="The cached directory " + bescachedir;
380 err_mesg = err_mesg + " doesn't exist. ";
381 throw InternalErr(__FILE__,__LINE__,err_mesg);
382
383 }
384 else {
385 if(true == S_ISDIR(sb.st_mode)) {
386 if(access(bescachedir.c_str(),R_OK|W_OK|X_OK) == -1) {
387 string err_mesg="The cached directory " + bescachedir;
388 err_mesg = err_mesg + " can NOT be read,written or executable.";
389 throw InternalErr(__FILE__,__LINE__,err_mesg);
390 }
391 }
392 else {
393 string err_mesg="The cached directory " + bescachedir;
394 err_mesg = err_mesg + " is not a directory.";
395 throw InternalErr(__FILE__,__LINE__,err_mesg);
396 }
397 }
398 }
399
400 string cache_fname=HDF5RequestHandler::get_latlon_disk_cachefile_prefix();
401
402 // Projection code,zone,sphere,pix,origin
403 cache_fname +=HDF5CFUtil::get_int_str(eos5_projcode);
404 cache_fname +=HDF5CFUtil::get_int_str(eos5_zone);
405 cache_fname +=HDF5CFUtil::get_int_str(eos5_sphere);
406 cache_fname +=HDF5CFUtil::get_int_str(eos5_pixelreg);
407 cache_fname +=HDF5CFUtil::get_int_str(eos5_origin);
408
409 cache_fname +=HDF5CFUtil::get_int_str(ydimsize);
410 cache_fname +=HDF5CFUtil::get_int_str(xdimsize);
411
412 // upleft,lowright
413 // HDF-EOS upleft,lowright,params use DDDDMMMSSS.6 digits. So choose %17.6f.
414 cache_fname +=HDF5CFUtil::get_double_str(point_left,17,6);
415 cache_fname +=HDF5CFUtil::get_double_str(point_upper,17,6);
416 cache_fname +=HDF5CFUtil::get_double_str(point_right,17,6);
417 cache_fname +=HDF5CFUtil::get_double_str(point_lower,17,6);
418
419 // According to HDF-EOS2 document, only 13 parameters are used.
420 for(int ipar = 0; ipar<13;ipar++) {
421 cache_fname+=HDF5CFUtil::get_double_str(eos5_params[ipar],17,6);
422 }
423
424 string cache_fpath = bescachedir + "/"+ cache_fname;
425 return cache_fpath;
426
427}
428void HDFEOS5CFMissLLArray::read_data_NOT_from_mem_cache_geo(bool add_cache,void*buf){
429
430 BESDEBUG("h5","Coming to read_data_NOT_from_mem_cache_geo "<<endl);
431 int nelms = -1;
432 vector<int>offset;
433 vector<int>count;
434 vector<int>step;
435
436
437 if (rank <= 0)
438 throw InternalErr (__FILE__, __LINE__,
439 "The number of dimension of this variable should be greater than 0");
440 else {
441
442 offset.resize(rank);
443 count.resize(rank);
444 step.resize(rank);
445 nelms = format_constraint (offset.data(), step.data(), count.data());
446 }
447
448 if (nelms <= 0)
449 throw InternalErr (__FILE__, __LINE__,
450 "The number of elments is negative.");
451
452 float start = 0.0;
453 float end = 0.0;
454
455 vector<float>val;
456 val.resize(nelms);
457
458
459 if (CV_LAT_MISS == cvartype) {
460
461 if (HE5_HDFE_GD_UL == eos5_origin || HE5_HDFE_GD_UR == eos5_origin) {
462
463 start = point_upper;
464 end = point_lower;
465
466 }
467 else {// (gridorigin == HE5_HDFE_GD_LL || gridorigin == HE5_HDFE_GD_LR)
468
469 start = point_lower;
470 end = point_upper;
471 }
472
473 if(ydimsize <=0)
474 throw InternalErr (__FILE__, __LINE__,
475 "The number of elments should be greater than 0.");
476
477 float lat_step = (end - start) /ydimsize;
478
479 // Now offset,step and val will always be valid. line 74 and 85 assure this.
480 if ( HE5_HDFE_CENTER == eos5_pixelreg ) {
481 for (int i = 0; i < nelms; i++)
482 val[i] = ((offset[0]+i*step[0] + 0.5F) * lat_step + start) / 1000000.0F;
483
484 // If the memory cache is turned on, we have to save all values to the buf
485 if(add_cache == true) {
486 vector<float>total_val;
487 total_val.resize(ydimsize);
488 for (int total_i = 0; total_i < ydimsize; total_i++)
489 total_val[total_i] = ((total_i + 0.5F) * lat_step + start) / 1000000.0F;
490 // Note: the float is size 4
491 memcpy(buf,total_val.data(),4*ydimsize);
492 }
493
494 }
495 else { // HE5_HDFE_CORNER
496 for (int i = 0; i < nelms; i++)
497 val[i] = ((float)(offset[0]+i * step[0])*lat_step + start) / 1000000.0F;
498
499 // If the memory cache is turned on, we have to save all values to the buf
500 if(add_cache == true) {
501 vector<float>total_val;
502 total_val.resize(ydimsize);
503 for (int total_i = 0; total_i < ydimsize; total_i++)
504 total_val[total_i] = ((float)(total_i) * lat_step + start) / 1000000.0F;
505 // Note: the float is size 4
506 memcpy(buf,total_val.data(),4*ydimsize);
507 }
508
509 }
510 }
511
512 if (CV_LON_MISS == cvartype) {
513
514 if (HE5_HDFE_GD_UL == eos5_origin || HE5_HDFE_GD_LL == eos5_origin) {
515
516 start = point_left;
517 end = point_right;
518
519 }
520 else {// (gridorigin == HE5_HDFE_GD_UR || gridorigin == HE5_HDFE_GD_LR)
521
522 start = point_right;
523 end = point_left;
524 }
525
526 if(xdimsize <=0)
527 throw InternalErr (__FILE__, __LINE__,
528 "The number of elments should be greater than 0.");
529 float lon_step = (end - start) /xdimsize;
530
531 if (HE5_HDFE_CENTER == eos5_pixelreg) {
532 for (int i = 0; i < nelms; i++)
533 val[i] = ((offset[0] + i *step[0] + 0.5F) * lon_step + start ) / 1000000.0F;
534
535 // If the memory cache is turned on, we have to save all values to the buf
536 if(add_cache == true) {
537 vector<float>total_val;
538 total_val.resize(xdimsize);
539 for (int total_i = 0; total_i < xdimsize; total_i++)
540 total_val[total_i] = ((total_i+0.5F) * lon_step + start) / 1000000.0F;
541 // Note: the float is size 4
542 memcpy(buf,total_val.data(),4*xdimsize);
543 }
544
545 }
546 else { // HE5_HDFE_CORNER
547 for (int i = 0; i < nelms; i++)
548 val[i] = ((float)(offset[0]+i*step[0]) * lon_step + start) / 1000000.0F;
549
550 // If the memory cache is turned on, we have to save all values to the buf
551 if(add_cache == true) {
552 vector<float>total_val;
553 total_val.resize(xdimsize);
554 for (int total_i = 0; total_i < xdimsize; total_i++)
555 total_val[total_i] = ((float)(total_i) * lon_step + start) / 1000000.0F;
556 // Note: the float is size 4
557 memcpy(buf,total_val.data(),4*xdimsize);
558 }
559 }
560 }
561
562#if 0
563for (int i =0; i <nelms; i++)
564"h5","final data val "<< i <<" is " << val[i] <<endl;
565#endif
566
567 set_value ((dods_float32 *) val.data(), nelms);
568
569
570 return;
571}
572
573#if 0
574// parse constraint expr. and make hdf5 coordinate point location.
575// return number of elements to read.
576int
577HDFEOS5CFMissLLArray::format_constraint (int *offset, int *step, int *count)
578{
579 long nels = 1;
580 int id = 0;
581
582 Dim_iter p = dim_begin ();
583
584 while (p != dim_end ()) {
585
586 int start = dimension_start (p, true);
587 int stride = dimension_stride (p, true);
588 int stop = dimension_stop (p, true);
589
590 // Check for illegal constraint
591 if (start > stop) {
592 ostringstream oss;
593
594 oss << "Array/Grid hyperslab start point "<< start <<
595 " is greater than stop point " << stop <<".";
596 throw Error(malformed_expr, oss.str());
597 }
598
599
600
601 offset[id] = start;
602 step[id] = stride;
603 count[id] = ((stop - start) / stride) + 1; // count of elements
604 nels *= count[id]; // total number of values for variable
605
606 BESDEBUG ("h5",
607 "=format_constraint():"
608 << "id=" << id << " offset=" << offset[id]
609 << " step=" << step[id]
610 << " count=" << count[id]
611 << endl);
612
613 id++;
614 p++;
615 }
616
617 return nels;
618}
619
620#endif
include the entry functions to execute the handlers
This class specifies the retrieval of the missing lat/lon values for HDF-EOS5 products.
static HDF5DiskCache * get_instance(const long, const std::string &, const std::string &)
static ssize_t read_buffer_from_file(int fd, void *buf, size_t)
Getting a subset of a variable.
Definition: HDF5CFUtil.cc:1194