bes  Updated for version 3.20.10
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 #include <cassert>
28 
29 #include <STARE.h>
30 #include <SpatialRange.h>
31 
32 //#include <hdf5.h>
33 #include <netcdf.h>
34 
35 #include <libdap/BaseType.h>
36 #include <libdap/Float64.h>
37 #include <libdap/Str.h>
38 #include <libdap/Array.h>
39 #include <libdap/Grid.h>
40 #include <libdap/Structure.h>
41 
42 #include <libdap/DMR.h>
43 #include <libdap/D4RValue.h>
44 
45 #include <libdap/Byte.h>
46 #include <libdap/Int16.h>
47 #include <libdap/Int32.h>
48 #include <libdap/UInt16.h>
49 #include <libdap/UInt32.h>
50 #include <libdap/Float32.h>
51 #include <libdap/UInt64.h>
52 #include <libdap/Int64.h>
53 
54 #include "BESDebug.h"
55 #include "BESUtil.h"
56 #include "BESInternalError.h"
57 #include "BESSyntaxUserError.h"
58 
59 #include "StareFunctions.h"
60 #include "GeoFile.h"
61 
62 // Used with BESDEBUG
63 #define STARE_FUNC "stare"
64 
65 using namespace libdap;
66 using namespace std;
67 
68 namespace functions {
69 
70 // These default values can be overridden using BES keys.
71 // See StareFunctions.h jhrg 5/21/20
72 // If stare_storage_path is empty, expect the sidecar file in the same
73 // place as the data file. jhrg 5/26/20
74 string stare_storage_path = "";
75 
76 // TODO Remove this once change-over is complete. jhrg 6/17/21
77 string stare_sidecar_suffix = "_stare.nc";
78 
90 ostream & operator << (ostream &out, const stare_matches &m)
91 {
92  assert(m.stare_indices.size() == m.row_indices.size()
93  && m.row_indices.size() == m.col_indices.size());
94 
95  auto ti = m.target_indices.begin();
96  auto si = m.stare_indices.begin();
97  auto xi = m.row_indices.begin();
98  auto yi = m.col_indices.begin();
99 
100  while (si != m.stare_indices.end()) {
101  out << "Target: " << *ti++ << ", Dataset Index: " << *si++ << ", coord: row: " << *xi++ << ", col: " << *yi++ << endl;
102  }
103 
104  return out;
105 }
106 
107 static void
108 extract_stare_index_array(Array *var, vector<STARE_ArrayIndexSpatialValue> &values)
109 {
110  if (var->var()->type() != dods_uint64_c)
111  throw BESSyntaxUserError("STARE server function passed an invalid Index array (" + var->name()
112  + " is type: " + var->var()->type_name() + ").", __FILE__, __LINE__);
113 
114  values.resize(var->length());
115  var->value((dods_uint64*)&values[0]); // Extract the values of 'var' to 'values'
116 }
117 
130 //int cmpSpatial(STARE_ArrayIndexSpatialValue a_, STARE_ArrayIndexSpatialValue b_);
131 
139 bool
140 target_in_dataset(const vector<STARE_ArrayIndexSpatialValue> &target_indices,
141  const vector<STARE_ArrayIndexSpatialValue> &data_stare_indices) {
142 #if 1
143  // this took 0.23s and worked.
144  // Changes to the range-for loop, fixed the type (was unsigned long long
145  // which works on OSX but not CentOS7). jhrg 11/5/19
146  for (const STARE_ArrayIndexSpatialValue &i : target_indices) {
147  for (const STARE_ArrayIndexSpatialValue &j :data_stare_indices ) {
148  // Check to see if the index 'i' overlaps the index 'j'. The cmpSpatial()
149  // function returns -1, 0, 1 depending on i in j, no overlap or, j in i.
150  // testing for !0 covers the general overlap case.
151  int result = cmpSpatial(i, j);
152  if (result != 0)
153  return true;
154  }
155  }
156 #else
157  // For one sample MOD05 file, this took 5.6s and failed to work.
158  // Problem: seg fault.
159 
160  // Initialize the skiplists for the search
161  auto r = SpatialRange(data_stare_indices);
162  cerr << "built spatial range" << endl;
163  for (const dods_uint64 &sid : target_indices) {
164  if( r.contains(sid) ) { return true; }
165  }
166 #endif
167  return false;
168 }
169 
186 unsigned int
187 count(const vector<STARE_ArrayIndexSpatialValue> &target_indices,
188  const vector<STARE_ArrayIndexSpatialValue> &dataset_indices,
189  bool all_dataset_matches /*= false*/) {
190  unsigned int counter = 0;
191  for (const auto &i : target_indices) {
192  for (const auto &j : dataset_indices)
193  // Here we are counting the number of target indices that overlap the
194  // dataset indices.
195  if (cmpSpatial(i, j) != 0) {
196  counter++;
197  BESDEBUG(STARE_FUNC, "Matching (dataset, target) indices: " << i << ", " << j << endl);
198  if (!all_dataset_matches)
199  break; // exit the inner loop
200  }
201  }
202 
203  return counter;
204 }
205 
214 unique_ptr<stare_matches>
215 stare_subset_helper(const vector<STARE_ArrayIndexSpatialValue> &target_indices,
216  const vector<STARE_ArrayIndexSpatialValue> &dataset_indices,
217  size_t dataset_rows, size_t dataset_cols)
218 {
219  //auto subset = new stare_matches;
220  unique_ptr<stare_matches> subset(new stare_matches());
221 
222  auto sid_iter = dataset_indices.begin();
223  for (size_t i = 0; i < dataset_rows; ++i) {
224  for (size_t j = 0; j < dataset_cols; ++j) {
225  auto sid = *sid_iter++;
226  for (const auto &target : target_indices) {
227  if (cmpSpatial(sid, target) != 0) { // != 0 --> sid is in target OR target is in sid
228  subset->add(i, j, sid, target);
229  }
230  }
231  }
232  }
233 
234  return subset;
235 }
236 
251 template <class T>
252 void stare_subset_array_helper(vector<T> &result_data, const vector<T> &src_data,
253  const vector<STARE_ArrayIndexSpatialValue> &target_indices,
254  const vector<STARE_ArrayIndexSpatialValue> &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 auto &i : dataset_indices) {
262  for (const auto &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 
291 template <class T>
292 void StareSubsetArrayFunction::build_masked_data(Array *dependent_var,
293  const vector<STARE_ArrayIndexSpatialValue> &dep_var_stare_indices,
294  const vector<STARE_ArrayIndexSpatialValue> &target_s_indices, T mask_value,
295  unique_ptr<Array> &result) {
296  vector<T> src_data(dependent_var->length());
297  dependent_var->read(); // TODO Do we need to call read() here? jhrg 6/16/20
298  dependent_var->value(&src_data[0]);
299 
300  //T mask_value = 0; // TODO This should use the value in mask_val_var. jhrg 6/16/20
301  vector<T> result_data(dependent_var->length(), mask_value);
302 
303  stare_subset_array_helper(result_data, src_data, target_s_indices, dep_var_stare_indices);
304 
305  result->set_value(result_data, result_data.size());
306 }
307 
308 void
309 read_stare_indices_from_function_argument(BaseType *raw_stare_indices,
310  vector<STARE_ArrayIndexSpatialValue> &s_indices) {
311 
312  Array *stare_indices = dynamic_cast<Array *>(raw_stare_indices);
313  if (stare_indices == nullptr)
314  throw BESSyntaxUserError(
315  "Expected an Array but found a " + raw_stare_indices->type_name(), __FILE__, __LINE__);
316 
317  if (stare_indices->var()->type() != dods_uint64_c)
318  throw BESSyntaxUserError(
319  "Expected an Array of UInt64 values but found an Array of " + stare_indices->var()->type_name(),
320  __FILE__, __LINE__);
321 
322  stare_indices->read();
323 
324  extract_stare_index_array(stare_indices, s_indices);
325 }
326 
337 BaseType *
338 StareIntersectionFunction::stare_intersection_dap4_function(D4RValueList *args, DMR &dmr)
339 {
340  if (args->size() != 2) {
341  ostringstream oss;
342  oss << "stare_intersection(): Expected two arguments, but got " << args->size();
343  throw BESSyntaxUserError(oss.str(), __FILE__, __LINE__);
344  }
345 
346  BaseType *dependent_var = args->get_rvalue(0)->value(dmr);
347  BaseType *raw_stare_indices = args->get_rvalue(1)->value(dmr);
348 
349  unique_ptr<GeoFile> gf(new GeoFile(dmr.filename()));
350 
351  // Read the stare indices for the dependent var from the sidecar file.
352  vector<STARE_ArrayIndexSpatialValue> dep_var_stare_indices;
353  gf->get_stare_indices(dependent_var->name(), dep_var_stare_indices);
354 
355  // Put the stare indices passed into the function into a vector<>
356  vector<STARE_ArrayIndexSpatialValue> target_s_indices;
357  read_stare_indices_from_function_argument(raw_stare_indices, target_s_indices);
358 
359  // Are any of the target indices covered by this variable
360  bool status = target_in_dataset(target_s_indices, dep_var_stare_indices);
361 
362  unique_ptr<Int32> result(new Int32("result"));
363  if (status) {
364  result->set_value(1);
365  }
366  else {
367  result->set_value(0);
368  }
369 
370  return result.release();
371 }
372 
391 BaseType *
392 StareCountFunction::stare_count_dap4_function(D4RValueList *args, DMR &dmr)
393 {
394  if (args->size() != 2) {
395  ostringstream oss;
396  oss << "stare_intersection(): Expected two arguments, but got " << args->size();
397  throw BESSyntaxUserError(oss.str(), __FILE__, __LINE__);
398  }
399 
400  BaseType *dependent_var = args->get_rvalue(0)->value(dmr);
401  BaseType *raw_stare_indices = args->get_rvalue(1)->value(dmr);
402 
403  unique_ptr<GeoFile> gf(new GeoFile(dmr.filename()));
404 
405  // Read the stare indices for the dependent var from the sidecar file.
406  vector<STARE_ArrayIndexSpatialValue> dep_var_stare_indices;
407  gf->get_stare_indices(dependent_var->name(), dep_var_stare_indices);
408 
409  // Put the stare indices passed into the function into a vector<>
410  vector<STARE_ArrayIndexSpatialValue> target_s_indices;
411  read_stare_indices_from_function_argument(raw_stare_indices, target_s_indices);
412 
413  int num = count(target_s_indices, dep_var_stare_indices, false);
414 
415  unique_ptr<Int32> result(new Int32("result"));
416  result->set_value(num);
417  return result.release();
418 }
419 
434 BaseType *
435 StareSubsetFunction::stare_subset_dap4_function(D4RValueList *args, DMR &dmr)
436 {
437  if (args->size() != 2) {
438  ostringstream oss;
439  oss << "stare_subset(): Expected two arguments, but got " << args->size();
440  throw BESSyntaxUserError(oss.str(), __FILE__, __LINE__);
441  }
442 
443  BaseType *dependent_var = args->get_rvalue(0)->value(dmr);
444  BaseType *raw_stare_indices = args->get_rvalue(1)->value(dmr);
445 
446  unique_ptr<GeoFile> gf(new GeoFile(dmr.filename()));
447 
448  // Read the stare indices for the dependent var from the sidecar file.
449  vector<STARE_ArrayIndexSpatialValue> dep_var_stare_indices;
450  gf->get_stare_indices(dependent_var->name(), dep_var_stare_indices);
451 
452  // Put the stare indices passed into the function into a vector<>
453  vector<STARE_ArrayIndexSpatialValue> target_s_indices;
454  read_stare_indices_from_function_argument(raw_stare_indices, target_s_indices);
455 
456  unique_ptr <stare_matches> subset = stare_subset_helper(target_s_indices, dep_var_stare_indices,
457  gf->get_variable_rows(dependent_var->name()),
458  gf->get_variable_cols(dependent_var->name()));
459 
460  // When no subset is found (none of the target indices match those in the dataset)
461  if (subset->stare_indices.size() == 0) {
462  subset->stare_indices.push_back(0);
463  subset->target_indices.push_back(0);
464  subset->row_indices.push_back(-1);
465  subset->col_indices.push_back(-1);
466  }
467  // Transfer values to a Structure
468  unique_ptr<Structure> result(new Structure("result"));
469 
470  unique_ptr<Array> stare(new Array("stare", new UInt64("stare")));
471  stare->set_value((dods_uint64*)&(subset->stare_indices[0]), subset->stare_indices.size());
472  stare->append_dim(subset->stare_indices.size());
473  result->add_var_nocopy(stare.release());
474 
475  unique_ptr<Array> target(new Array("target", new UInt64("target")));
476  target->set_value((dods_uint64*)&(subset->target_indices[0]), subset->target_indices.size());
477  target->append_dim(subset->target_indices.size());
478  result->add_var_nocopy(target.release());
479 
480  unique_ptr<Array> x(new Array("row", new Int32("row")));
481  x->set_value(subset->row_indices, subset->row_indices.size());
482  x->append_dim(subset->row_indices.size());
483  result->add_var_nocopy(x.release());
484 
485  unique_ptr<Array> y(new Array("col", new Int32("col")));
486  y->set_value(subset->col_indices, subset->col_indices.size());
487  y->append_dim(subset->col_indices.size());
488  result->add_var_nocopy(y.release());
489 
490  return result.release();
491 }
492 
501 double get_double_value(BaseType *btp)
502 {
503  switch (btp->type()) {
504  case libdap::dods_float64_c:
505  return dynamic_cast<Float64*>(btp)->value();
506  case libdap::dods_int64_c:
507  return dynamic_cast<Int64*>(btp)->value();
508  case libdap::dods_uint64_c:
509  return dynamic_cast<UInt64*>(btp)->value();
510  default: {
511  throw BESSyntaxUserError(string("Expected a constant value, but got ").append(btp->type_name())
512  .append(" instead."), __FILE__, __LINE__);
513  }
514  }
515 }
516 
517 BaseType *
518 StareSubsetArrayFunction::stare_subset_array_dap4_function(D4RValueList *args, DMR &dmr)
519 {
520  if (args->size() != 3) {
521  ostringstream oss;
522  oss << "stare_subset_array(): Expected three arguments, but got " << args->size();
523  throw BESSyntaxUserError(oss.str(), __FILE__, __LINE__);
524  }
525 
526  Array *dependent_var = dynamic_cast<Array*>(args->get_rvalue(0)->value(dmr));
527  if (!dependent_var)
528  throw BESSyntaxUserError("Expected an Array as teh first argument to stare_subset_array()", __FILE__, __LINE__);
529 
530  BaseType *mask_val_var = args->get_rvalue(1)->value(dmr);
531  double mask_value = get_double_value(mask_val_var);
532 
533  BaseType *raw_stare_indices = args->get_rvalue(2)->value(dmr);
534 
535  unique_ptr<GeoFile> gf(new GeoFile(dmr.filename()));
536 
537  // Read the stare indices for the dependent var from the sidecar file.
538  vector<STARE_ArrayIndexSpatialValue> dep_var_stare_indices;
539  gf->get_stare_indices(dependent_var->name(), dep_var_stare_indices);
540 
541  // Put the stare indices passed into the function into a vector<>
542  vector<STARE_ArrayIndexSpatialValue> target_s_indices;
543  read_stare_indices_from_function_argument(raw_stare_indices, target_s_indices);
544 
545  // ptr_duplicate() does not copy data values
546  unique_ptr<Array> result(dynamic_cast<Array*>(dependent_var->ptr_duplicate()));
547 
548  // FIXME Add more types. jhrg 6/17/20
549  switch(dependent_var->var()->type()) {
550  case dods_int16_c: {
551  build_masked_data<dods_int16>(dependent_var, dep_var_stare_indices, target_s_indices, mask_value, result);
552  break;
553  }
554  case dods_float32_c: {
555  build_masked_data<dods_float32>(dependent_var, dep_var_stare_indices, target_s_indices, mask_value, result);
556  break;
557  }
558 
559  default:
560  throw BESInternalError(string("stare_subset_array() failed: Unsupported array element type (")
561  + dependent_var->var()->type_name() + ").", __FILE__, __LINE__);
562  }
563 
564  return result.release();
565 }
566 
574 STARE_SpatialIntervals
575 stare_box_helper(const vector<point> &points, int resolution) {
576  LatLonDegrees64ValueVector latlonbox;
577  for (auto &p: points) {
578  latlonbox.push_back(LatLonDegrees64(p.lat, p.lon));
579  }
580 
581  STARE index(27, 6);
582  return index.CoverBoundingBoxFromLatLonDegrees(latlonbox, resolution);
583 }
584 
596 STARE_SpatialIntervals
597 stare_box_helper(const point &top_left, const point &bottom_right, int resolution) {
598  LatLonDegrees64ValueVector latlonbox;
599  latlonbox.push_back(LatLonDegrees64(top_left.lat, top_left.lon));
600  latlonbox.push_back(LatLonDegrees64(top_left.lat, bottom_right.lon));
601  latlonbox.push_back(LatLonDegrees64(bottom_right.lat, bottom_right.lon));
602  latlonbox.push_back(LatLonDegrees64(bottom_right.lat, top_left.lon));
603 
604  STARE index(27, 6); // Hack values.
605  return index.CoverBoundingBoxFromLatLonDegrees(latlonbox, resolution);
606 }
607 
608 BaseType *
609 StareBoxFunction::stare_box_dap4_function(libdap::D4RValueList *args, libdap::DMR &dmr)
610 {
611  STARE_SpatialIntervals sivs;
612 
613  if (args->size() == 4) {
614  // build cover from simple box
615  double tl_lat = get_double_value(args->get_rvalue(0)->value(dmr));
616  double tl_lon = get_double_value(args->get_rvalue(1)->value(dmr));
617  double br_lat = get_double_value(args->get_rvalue(2)->value(dmr));
618  double br_lon = get_double_value(args->get_rvalue(3)->value(dmr));
619 
620  sivs = stare_box_helper(point(tl_lat, tl_lon), point(br_lat, br_lon));
621  }
622 #if 1
623  else if (args->size() >= 6 && (args->size() % 2) == 0) {
624  // build cover from list of three or more points (lat, lon)
625  bool flip_flop = false;
626  double lat_value = -90;
627  vector<point> points;
628  for (auto &arg: *args) {
629  if (flip_flop) {
630  point pt(lat_value, get_double_value(arg->value(dmr)));
631  points.push_back(pt);
632  flip_flop = false;
633  }
634  else {
635  lat_value = get_double_value(arg->value(dmr));
636  flip_flop = true;
637  }
638  }
639 
640  sivs = stare_box_helper(points);
641  }
642 #endif
643  else {
644  ostringstream oss;
645  oss << "stare_box(): Expected four corner lat/lon values or a list of three or more points, but got "
646  << args->size() << " values.";
647  throw BESSyntaxUserError(oss.str(), __FILE__, __LINE__);
648  }
649 
650  unique_ptr<Array> cover(new Array("cover", new UInt64("cover")));
651  cover->set_value((dods_uint64*)(&sivs[0]), sivs.size());
652  cover->append_dim(sivs.size());
653 
654  return cover.release();
655 }
656 
657 
658 
659 
660 } // 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
Hold the result from the subset helper function as a collection of vectors.