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