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
65using namespace libdap;
66using namespace std;
67
68namespace 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
74string stare_storage_path = "";
75
76// TODO Remove this once change-over is complete. jhrg 6/17/21
77string stare_sidecar_suffix = "_stare.nc";
78
90ostream & 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
107static void
108extract_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
139bool
140target_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
186unsigned int
187count(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
214unique_ptr<stare_matches>
215stare_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
251template <class T>
252void 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
291template <class T>
292void 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
308void
309read_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
337BaseType *
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
391BaseType *
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
434BaseType *
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
501double 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
517BaseType *
518StareSubsetArrayFunction::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
574STARE_SpatialIntervals
575stare_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
596STARE_SpatialIntervals
597stare_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
608BaseType *
609StareBoxFunction::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
static libdap::BaseType * stare_count_dap4_function(libdap::D4RValueList *args, libdap::DMR &dmr)
Count the number of STARE indices in the arg that overlap the indices of this dataset.
static libdap::BaseType * stare_intersection_dap4_function(libdap::D4RValueList *args, libdap::DMR &dmr)
Return true/false indicating that the given stare indices intersect the variables.
static libdap::BaseType * stare_subset_dap4_function(libdap::D4RValueList *args, libdap::DMR &dmr)
For the given target STARE indices, return the overlapping dataset X, Y, and STARE indices.
Hold the result from the subset helper function as a collection of vectors.