bes  Updated for version 3.20.8
StareFunctions.cc
1 // This file is part of libdap, A C++ implementation of the OPeNDAP Data
2 // Access Protocol.
3 
4 // Copyright (c) 2019 OPeNDAP, Inc.
5 // Authors: Kodi Neumiller <kneumiller@opendap.org>
6 //
7 // This library is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 //
12 // This library is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 // Lesser General Public License for more details.
16 //
17 // You should have received a copy of the GNU Lesser General Public
18 // License along with this library; if not, write to the Free Software
19 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20 //
21 // You can contact OPeNDAP, Inc. at PO Box 112, Saunderstown, RI. 02874-0112.
22 
23 #include "config.h"
24 
25 #include <sstream>
26 #include <memory>
27 
28 #include <STARE.h>
29 #include <hdf5.h>
30 
31 #include <BaseType.h>
32 #include <Float64.h>
33 #include <Str.h>
34 #include <Array.h>
35 #include <Grid.h>
36 #include <Structure.h>
37 
38 #include <DMR.h>
39 #include <D4RValue.h>
40 
41 #include <Byte.h>
42 #include <Int16.h>
43 #include <Int32.h>
44 #include <UInt16.h>
45 #include <UInt32.h>
46 #include <Float32.h>
47 #include <Int64.h>
48 #include <UInt64.h>
49 
50 #include "BESDebug.h"
51 #include "BESUtil.h"
52 #include "BESInternalError.h"
53 #include "BESSyntaxUserError.h"
54 
55 #include "StareFunctions.h"
56 
57 // Used with BESDEBUG
58 #define STARE "stare"
59 
60 using namespace libdap;
61 using namespace std;
62 
63 namespace functions {
64 
65 // These default values can be overridden using BES keys.
66 // See StareFunctions.h jhrg 5/21/20
67 // If stare_storage_path is empty, expect the sidecar file in the same
68 // place as the data file. jhrg 5/26/20
69 string stare_storage_path = "";
70 string stare_sidecar_suffix = "_sidecar";
71 
83 ostream & operator << (ostream &out, const stare_matches &m)
84 {
85  assert(m.stare_indices.size() == m.x_indices.size()
86  && m.x_indices.size() == m.y_indices.size());
87 
88  auto ti = m.target_indices.begin();
89  auto si = m.stare_indices.begin();
90  auto xi = m.x_indices.begin();
91  auto yi = m.y_indices.begin();
92 
93  while (si != m.stare_indices.end()) {
94  out << "Target: " << *ti++ << ", Dataset Index: " << *si++ << ", coord: x: " << *xi++ << ", y: " << *yi++ << endl;
95  }
96 
97  return out;
98 }
99 
100 #if 0
101 // May need to be moved to libdap/util
102 // This helper function assumes 'var' is the correct size.
103 // Made this static to limit its scope to this file. jhrg 11/7/19
104 
111 static vector<dods_uint64> *extract_uint64_array(Array *var) {
112  assert(var);
113 
114  int length = var->length();
115 
116  auto *newVar = new vector<dods_uint64>;
117  newVar->resize(length);
118  var->value(&(*newVar)[0]); // Extract the values of 'var' to 'newVar'
119 
120  return newVar;
121 }
122 #endif
123 
124 void
125 extract_uint64_array(Array *var, vector<dods_uint64> &values) {
126 
127  values.resize(var->length());
128  var->value(&values[0]); // Extract the values of 'var' to 'values'
129 }
130 
131 
144 //int cmpSpatial(STARE_ArrayIndexSpatialValue a_, STARE_ArrayIndexSpatialValue b_) {
145 
152 bool
153 target_in_dataset(const vector<dods_uint64> &targetIndices, const vector<dods_uint64> &dataStareIndices) {
154  // Changes to the range-for loop, fixed the type (was unsigned long long
155  // which works on OSX but not CentOS7). jhrg 11/5/19
156  for (const dods_uint64 &i : targetIndices) {
157  for (const dods_uint64 &j :dataStareIndices ) {
158  // Check to see if the index 'i' overlaps the index 'j'. The cmpSpatial()
159  // function returns -1, 0, 1 depending on i in j, no overlap or, j in i.
160  // testing for !0 covers the general overlap case.
161  int result = cmpSpatial(i, j);
162  if (result != 0)
163  return true;
164  }
165  }
166 
167  return false;
168 }
169 
184 unsigned int
185 count(const vector<dods_uint64> &target_indices, const vector<dods_uint64> &dataset_indices, bool all_target_matches /*= false*/) {
186  unsigned int counter = 0;
187  for (const dods_uint64 &i : dataset_indices) {
188  for (const dods_uint64 &j : target_indices)
189  // Here we are counting the number of target indices that overlap the
190  // dataset indices.
191  if (cmpSpatial(i, j) != 0) {
192  counter++;
193  BESDEBUG(STARE, "Matching (dataset, target) indices: " << i << ", " << j << endl);
194  if (!all_target_matches)
195  break; // exit the inner loop
196  }
197  }
198 
199  return counter;
200 }
201 
212 unique_ptr<stare_matches>
213 stare_subset_helper(const vector<dods_uint64> &target_indices, const vector<dods_uint64> &dataset_indices,
214  const vector<int> &dataset_x_coords, const vector<int> &dataset_y_coords)
215 {
216  assert(dataset_indices.size() == dataset_x_coords.size());
217  assert(dataset_indices.size() == dataset_y_coords.size());
218 
219  //auto subset = new stare_matches;
220  unique_ptr<stare_matches> subset(new stare_matches());
221 
222  auto x = dataset_x_coords.begin();
223  auto y = dataset_y_coords.begin();
224  for (const dods_uint64 &i : dataset_indices) {
225  for (const dods_uint64 &j : target_indices) {
226  if (cmpSpatial(i, j) != 0) { // != 0 --> i is in j OR j is in i
227  subset->add(*x, *y, i, j);
228  // TODO Add a break call here? jhrg 6/17/20
229  }
230  }
231  ++x;
232  ++y;
233  }
234 
235  return subset;
236 }
237 
252 template <class T>
253 void stare_subset_array_helper(vector<T> &result_data, const vector<T> &src_data,
254  const vector<dods_uint64> &target_indices, const vector<dods_uint64> &dataset_indices)
255 {
256  assert(dataset_indices.size() == src_data.size());
257  assert(dataset_indices.size() == result_data.size());
258 
259  auto r = result_data.begin();
260  auto s = src_data.begin();
261  for (const dods_uint64 &i : dataset_indices) {
262  for (const dods_uint64 &j : target_indices) {
263  if (cmpSpatial(i, j) != 0) { // != 0 --> i is in j OR j is in i
264  *r = *s;
265  break;
266  }
267  }
268  ++r; ++s;
269  }
270 }
271 
289 template <class T>
290 void StareSubsetArrayFunction::build_masked_data(Array *dependent_var, const vector<dods_uint64> &dep_var_stare_indices,
291  const vector<dods_uint64> &target_s_indices, unique_ptr<Array> &result) {
292  vector<T> src_data(dependent_var->length());
293  dependent_var->read(); // TODO Do we need to call read() here? jhrg 6/16/20
294  dependent_var->value(&src_data[0]);
295 
296  T mask_value = 0; // TODO This should use the value in mask_val_var. jhrg 6/16/20
297  vector<T> result_data(dependent_var->length(), mask_value);
298 
299  stare_subset_array_helper(result_data, src_data, target_s_indices, dep_var_stare_indices);
300 
301  result->set_value(result_data, result_data.size());
302 }
303 
314 string
315 get_sidecar_file_pathname(const string &pathName, const string &token)
316 {
317  string granuleName = pathName;
318  if (!stare_storage_path.empty()) {
319  granuleName = pathName.substr(pathName.find_last_of('/') + 1);
320  }
321 
322  size_t findDot = granuleName.find_last_of('.');
323  // Added extraction of the extension since the files won't always be *.h5
324  // also switched to .append() instead of '+' because the former is faster.
325  // jhrg 11/5/19
326  string extension = granuleName.substr(findDot); // ext includes the dot
327  granuleName = granuleName.substr(0, findDot).append(token).append(extension);
328 
329  if (!stare_storage_path.empty()) {
330  // Above the path has been removed
331  return BESUtil::pathConcat(stare_storage_path, granuleName);
332  }
333  else {
334  // stare_storage_path is empty, granuleName is the full path
335  return granuleName;
336  }
337 
338 }
339 
345 void
346 get_sidecar_int32_values(const string &filename, const string &variable, vector<dods_int32> &values)
347 {
348  //Read the file and store the datasets
349  hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
350  if (file < 0)
351  throw BESInternalError("Could not open file " + filename, __FILE__, __LINE__);
352 
353  hid_t dataset = H5Dopen(file, variable.c_str(), H5P_DEFAULT);
354  if (dataset < 0)
355  throw BESInternalError(string("Could not open dataset: ").append(variable), __FILE__, __LINE__);
356 
357  hid_t dspace = H5Dget_space(dataset);
358  const int ndims = H5Sget_simple_extent_ndims(dspace);
359  vector<hsize_t> dims(ndims);
360 
361  //Get the size of the dimension so that we know how big to make the memory space
362  //Each of the dataspaces should be the same size, if in the future they are different
363  // sizes then the size of each dataspace will need to be calculated.
364  H5Sget_simple_extent_dims(dspace, &dims[0], NULL);
365 
366  //We need to get the filespace and memspace before reading the values from each dataset
367  hid_t filespace = H5Dget_space(dataset);
368 
369  hid_t memspace = H5Screate_simple(ndims, &dims[0], NULL);
370 
371  //Get the number of elements in the dataspace and use that to appropriate the proper size of the vectors
372  values.resize(H5Sget_select_npoints(filespace));
373 
374  //Read the data file and store the values of each dataset into an array
375  H5Dread(dataset, H5T_NATIVE_INT, memspace, filespace, H5P_DEFAULT, &values[0]);
376 }
377 
384 void
385 get_sidecar_uint64_values(const string &filename, BaseType */*variable*/, vector<dods_uint64> &values)
386 {
387  //Read the file and store the datasets
388  hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
389  if (file < 0)
390  throw BESInternalError("Could not open file " + filename, __FILE__, __LINE__);
391 
392  // Here we look up the name of the stare index data for 'variable.' For now,
393  // use the name stored in 's_index_name'. jhrg 6/3/20
394 
395  hid_t dataset = H5Dopen(file, s_index_name.c_str(), H5P_DEFAULT);
396  if (dataset < 0)
397  throw BESInternalError(string("Could not open dataset: ").append(s_index_name), __FILE__, __LINE__);
398 
399  hid_t dspace = H5Dget_space(dataset);
400  const int ndims = H5Sget_simple_extent_ndims(dspace);
401  vector<hsize_t> dims(ndims);
402 
403  //Get the size of the dimension so that we know how big to make the memory space
404  //Each of the dataspaces should be the same size, if in the future they are different
405  // sizes then the size of each dataspace will need to be calculated.
406  H5Sget_simple_extent_dims(dspace, &dims[0], NULL);
407 
408  //We need to get the filespace and memspace before reading the values from each dataset
409  hid_t filespace = H5Dget_space(dataset);
410 
411  hid_t memspace = H5Screate_simple(ndims, &dims[0], NULL);
412 
413  //Get the number of elements in the dataspace and use that to appropriate the proper size of the vectors
414  values.resize(H5Sget_select_npoints(filespace));
415 
416  //Read the data file and store the values of each dataset into an array
417  H5Dread(dataset, H5T_NATIVE_ULLONG, memspace, filespace, H5P_DEFAULT, &values[0]);
418 }
419 
420 void
421 read_stare_indices_from_function_argument(BaseType *raw_stare_indices, vector<dods_uint64>&s_indices) {
422 
423  Array *stare_indices = dynamic_cast<Array *>(raw_stare_indices);
424  if (stare_indices == nullptr)
425  throw BESSyntaxUserError(
426  "Expected an Array but found a " + raw_stare_indices->type_name(), __FILE__, __LINE__);
427 
428  if (stare_indices->var()->type() != dods_uint64_c)
429  throw BESSyntaxUserError(
430  "Expected an Array of UInt64 values but found an Array of " + stare_indices->var()->type_name(),
431  __FILE__, __LINE__);
432 
433  stare_indices->read();
434 
435  extract_uint64_array(stare_indices, s_indices);
436 }
437 
448 BaseType *
449 StareIntersectionFunction::stare_intersection_dap4_function(D4RValueList *args, DMR &dmr)
450 {
451  if (args->size() != 2) {
452  ostringstream oss;
453  oss << "stare_intersection(): Expected two arguments, but got " << args->size();
454  throw BESSyntaxUserError(oss.str(), __FILE__, __LINE__);
455  }
456 
457  //Find the filename from the dmr
458  string fullPath = get_sidecar_file_pathname(dmr.filename(), stare_sidecar_suffix);
459 
460  BaseType *dependent_var = args->get_rvalue(0)->value(dmr);
461  BaseType *raw_stare_indices = args->get_rvalue(1)->value(dmr);
462 
463  //Read the data file and store the values of each dataset into an array
464  vector<dods_uint64> dep_var_stare_indices;
465  get_sidecar_uint64_values(fullPath, dependent_var, dep_var_stare_indices);
466 
467  // TODO: We can dump the values in 'stare_indices' here
468  vector<dods_uint64> target_s_indices;
469  read_stare_indices_from_function_argument(raw_stare_indices, target_s_indices);
470 
471  bool status = target_in_dataset(target_s_indices, dep_var_stare_indices);
472 
473 #if 0
474  Int32 *result = new Int32("result");
475 #endif
476  unique_ptr<Int32> result(new Int32("result"));
477  if (status) {
478  result->set_value(1);
479  }
480  else {
481  result->set_value(0);
482  }
483 
484  return result.release();
485 }
486 
505 BaseType *
506 StareCountFunction::stare_count_dap4_function(D4RValueList *args, DMR &dmr)
507 {
508  if (args->size() != 2) {
509  ostringstream oss;
510  oss << "stare_intersection(): Expected two arguments, but got " << args->size();
511  throw BESSyntaxUserError(oss.str(), __FILE__, __LINE__);
512  }
513 
514  //Find the filename from the dmr
515  string fullPath = get_sidecar_file_pathname(dmr.filename(), stare_sidecar_suffix);
516 
517  BaseType *dependent_var = args->get_rvalue(0)->value(dmr);
518  BaseType *raw_stare_indices = args->get_rvalue(1)->value(dmr);
519 
520  //Read the data file and store the values of each dataset into an array
521  vector<dods_uint64> dep_var_stare_indices;
522  get_sidecar_uint64_values(fullPath, dependent_var, dep_var_stare_indices);
523 
524  // TODO: We can dump the values in 'stare_indices' here
525  vector<dods_uint64> target_s_indices;
526  read_stare_indices_from_function_argument(raw_stare_indices, target_s_indices);
527 
528  int num = count(target_s_indices, dep_var_stare_indices);
529 
530 #if 0
531  Int32 *result = new Int32("result");
532 #endif
533  unique_ptr<Int32> result(new Int32("result"));
534  result->set_value(num);
535  return result.release();
536 }
537 
552 BaseType *
553 StareSubsetFunction::stare_subset_dap4_function(D4RValueList *args, DMR &dmr)
554 {
555  if (args->size() != 2) {
556  ostringstream oss;
557  oss << "stare_subset(): Expected two arguments, but got " << args->size();
558  throw BESSyntaxUserError(oss.str(), __FILE__, __LINE__);
559  }
560 
561  //Find the filename from the dmr
562  string fullPath = get_sidecar_file_pathname(dmr.filename(), stare_sidecar_suffix);
563 
564  BaseType *dependent_var = args->get_rvalue(0)->value(dmr);
565  BaseType *raw_stare_indices = args->get_rvalue(1)->value(dmr);
566 
567  //Read the data file and store the values of each dataset into an array
568  vector<dods_uint64> dep_var_stare_indices;
569  get_sidecar_uint64_values(fullPath, dependent_var, dep_var_stare_indices);
570 
571  // TODO: We can dump the values in 'stare_indices' here
572  vector<dods_uint64> target_s_indices;
573  read_stare_indices_from_function_argument(raw_stare_indices, target_s_indices);
574 
575  vector<dods_int32> dataset_x_coords;
576  get_sidecar_int32_values(fullPath, "X", dataset_x_coords);
577  vector<dods_int32> dataset_y_coords;
578  get_sidecar_int32_values(fullPath, "Y", dataset_y_coords);
579 
580  unique_ptr <stare_matches> subset = stare_subset_helper(target_s_indices, dep_var_stare_indices, dataset_x_coords, dataset_y_coords);
581 
582  // When no subset is found (none of the target indices match those in the dataset)
583  if (subset->stare_indices.size() == 0) {
584  subset->stare_indices.push_back(0);
585  subset->target_indices.push_back(0);
586  subset->x_indices.push_back(-1);
587  subset->y_indices.push_back(-1);
588  }
589  // Transfer values to a Structure
590  unique_ptr<Structure> result(new Structure("result"));
591 
592  unique_ptr<Array> stare(new Array("stare", new UInt64("stare")));
593  stare->set_value(&(subset->stare_indices[0]), subset->stare_indices.size());
594  stare->append_dim(subset->stare_indices.size());
595  result->add_var_nocopy(stare.release());
596 
597  unique_ptr<Array> target(new Array("target", new UInt64("target")));
598  target->set_value(&(subset->target_indices[0]), subset->target_indices.size());
599  target->append_dim(subset->target_indices.size());
600  result->add_var_nocopy(target.release());
601 
602  unique_ptr<Array> x(new Array("x", new Int32("x")));
603  x->set_value(subset->x_indices, subset->x_indices.size());
604  x->append_dim(subset->x_indices.size());
605  result->add_var_nocopy(x.release());
606 
607  unique_ptr<Array> y(new Array("y", new Int32("y")));
608  y->set_value(subset->y_indices, subset->y_indices.size());
609  y->append_dim(subset->y_indices.size());
610  result->add_var_nocopy(y.release());
611 
612  return result.release();
613 }
614 
615 BaseType *
616 StareSubsetArrayFunction::stare_subset_array_dap4_function(D4RValueList *args, DMR &dmr)
617 {
618  if (args->size() != 3) {
619  ostringstream oss;
620  oss << "stare_subset_array(): Expected three arguments, but got " << args->size();
621  throw BESSyntaxUserError(oss.str(), __FILE__, __LINE__);
622  }
623 
624  //Find the filename from the dmr
625  string fullPath = get_sidecar_file_pathname(dmr.filename(), stare_sidecar_suffix);
626 
627  Array *dependent_var = dynamic_cast<Array*>(args->get_rvalue(0)->value(dmr));
628  if (!dependent_var)
629  throw BESSyntaxUserError("stare_subset_array() expected an Array as the first argument.", __FILE__, __LINE__);
630 
631  // TODO Get/use this.
632  BaseType *mask_val_var = args->get_rvalue(1)->value(dmr);
633 
634  Array *raw_stare_indices = dynamic_cast<Array*>(args->get_rvalue(2)->value(dmr));
635  if (!raw_stare_indices)
636  throw BESSyntaxUserError("stare_subset_array() expected an Array as the third argument.", __FILE__, __LINE__);
637 
638  //Read the data file and store the values of each dataset into an array
639  vector<dods_uint64> dep_var_stare_indices;
640  get_sidecar_uint64_values(fullPath, dependent_var, dep_var_stare_indices);
641 
642  vector<dods_uint64> target_s_indices;
643  read_stare_indices_from_function_argument(raw_stare_indices, target_s_indices);
644 
645  // ptr_duplicate() does not copy data values
646  unique_ptr<Array> result(static_cast<Array*>(dependent_var->ptr_duplicate()));
647 
648  // TODO Add more types. jhrg 6/17/20
649  switch(dependent_var->var()->type()) {
650  case dods_int16_c: {
651  build_masked_data<dods_int16>(dependent_var, dep_var_stare_indices, target_s_indices, result);
652  break;
653  }
654  case dods_float32_c: {
655  build_masked_data<dods_float32>(dependent_var, dep_var_stare_indices, target_s_indices, result);
656  break;
657  }
658 
659  default:
660  throw BESInternalError(string("stare_subset_array() failed: Unsupported array element type (")
661  + dependent_var->var()->type_name() + ").", __FILE__, __LINE__);
662  }
663 
664  return result.release();
665 }
666 
667 } // namespace functions
exception thrown if internal error encountered
error thrown if there is a user syntax error in the request or any other user error
static std::string pathConcat(const std::string &firstPart, const std::string &secondPart, char separator='/')
Concatenate path fragments making sure that they are separated by a single '/' character.
Definition: BESUtil.cc:772
Hold the result from the subset helper function as a collection of vectors.