30 #include <SpatialRange.h>
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>
42 #include <libdap/DMR.h>
43 #include <libdap/D4RValue.h>
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>
56 #include "BESInternalError.h"
57 #include "BESSyntaxUserError.h"
59 #include "StareFunctions.h"
63 #define STARE_FUNC "stare"
74 string stare_storage_path =
"";
77 string stare_sidecar_suffix =
"_stare.nc";
92 assert(m.stare_indices.size() == m.row_indices.size()
93 && m.row_indices.size() == m.col_indices.size());
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();
100 while (si != m.stare_indices.end()) {
101 out <<
"Target: " << *ti++ <<
", Dataset Index: " << *si++ <<
", coord: row: " << *xi++ <<
", col: " << *yi++ << endl;
108 extract_stare_index_array(Array *var, vector<STARE_ArrayIndexSpatialValue> &values)
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__);
114 values.resize(var->length());
115 var->value((dods_uint64*)&values[0]);
140 target_in_dataset(
const vector<STARE_ArrayIndexSpatialValue> &target_indices,
141 const vector<STARE_ArrayIndexSpatialValue> &data_stare_indices) {
146 for (
const STARE_ArrayIndexSpatialValue &i : target_indices) {
147 for (
const STARE_ArrayIndexSpatialValue &j :data_stare_indices ) {
151 int result = cmpSpatial(i, j);
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; }
187 count(
const vector<STARE_ArrayIndexSpatialValue> &target_indices,
188 const vector<STARE_ArrayIndexSpatialValue> &dataset_indices,
189 bool all_dataset_matches ) {
190 unsigned int counter = 0;
191 for (
const auto &i : target_indices) {
192 for (
const auto &j : dataset_indices)
195 if (cmpSpatial(i, j) != 0) {
197 BESDEBUG(STARE_FUNC,
"Matching (dataset, target) indices: " << i <<
", " << j << endl);
198 if (!all_dataset_matches)
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)
220 unique_ptr<stare_matches> subset(
new stare_matches());
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) {
228 subset->add(i, j, sid, target);
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)
256 assert(dataset_indices.size() == src_data.size());
257 assert(dataset_indices.size() == result_data.size());
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) {
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();
298 dependent_var->value(&src_data[0]);
301 vector<T> result_data(dependent_var->length(), mask_value);
303 stare_subset_array_helper(result_data, src_data, target_s_indices, dep_var_stare_indices);
305 result->set_value(result_data, result_data.size());
309 read_stare_indices_from_function_argument(BaseType *raw_stare_indices,
310 vector<STARE_ArrayIndexSpatialValue> &s_indices) {
312 Array *stare_indices =
dynamic_cast<Array *
>(raw_stare_indices);
313 if (stare_indices ==
nullptr)
315 "Expected an Array but found a " + raw_stare_indices->type_name(), __FILE__, __LINE__);
317 if (stare_indices->var()->type() != dods_uint64_c)
319 "Expected an Array of UInt64 values but found an Array of " + stare_indices->var()->type_name(),
322 stare_indices->read();
324 extract_stare_index_array(stare_indices, s_indices);
338 StareIntersectionFunction::stare_intersection_dap4_function(D4RValueList *args, DMR &dmr)
340 if (args->size() != 2) {
342 oss <<
"stare_intersection(): Expected two arguments, but got " << args->size();
346 BaseType *dependent_var = args->get_rvalue(0)->value(dmr);
347 BaseType *raw_stare_indices = args->get_rvalue(1)->value(dmr);
349 unique_ptr<GeoFile> gf(
new GeoFile(dmr.filename()));
352 vector<STARE_ArrayIndexSpatialValue> dep_var_stare_indices;
353 gf->get_stare_indices(dependent_var->name(), dep_var_stare_indices);
356 vector<STARE_ArrayIndexSpatialValue> target_s_indices;
357 read_stare_indices_from_function_argument(raw_stare_indices, target_s_indices);
360 bool status = target_in_dataset(target_s_indices, dep_var_stare_indices);
362 unique_ptr<Int32> result(
new Int32(
"result"));
364 result->set_value(1);
367 result->set_value(0);
370 return result.release();
392 StareCountFunction::stare_count_dap4_function(D4RValueList *args, DMR &dmr)
394 if (args->size() != 2) {
396 oss <<
"stare_intersection(): Expected two arguments, but got " << args->size();
400 BaseType *dependent_var = args->get_rvalue(0)->value(dmr);
401 BaseType *raw_stare_indices = args->get_rvalue(1)->value(dmr);
403 unique_ptr<GeoFile> gf(
new GeoFile(dmr.filename()));
406 vector<STARE_ArrayIndexSpatialValue> dep_var_stare_indices;
407 gf->get_stare_indices(dependent_var->name(), dep_var_stare_indices);
410 vector<STARE_ArrayIndexSpatialValue> target_s_indices;
411 read_stare_indices_from_function_argument(raw_stare_indices, target_s_indices);
413 int num = count(target_s_indices, dep_var_stare_indices,
false);
415 unique_ptr<Int32> result(
new Int32(
"result"));
416 result->set_value(num);
417 return result.release();
435 StareSubsetFunction::stare_subset_dap4_function(D4RValueList *args, DMR &dmr)
437 if (args->size() != 2) {
439 oss <<
"stare_subset(): Expected two arguments, but got " << args->size();
443 BaseType *dependent_var = args->get_rvalue(0)->value(dmr);
444 BaseType *raw_stare_indices = args->get_rvalue(1)->value(dmr);
446 unique_ptr<GeoFile> gf(
new GeoFile(dmr.filename()));
449 vector<STARE_ArrayIndexSpatialValue> dep_var_stare_indices;
450 gf->get_stare_indices(dependent_var->name(), dep_var_stare_indices);
453 vector<STARE_ArrayIndexSpatialValue> target_s_indices;
454 read_stare_indices_from_function_argument(raw_stare_indices, target_s_indices);
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()));
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);
468 unique_ptr<Structure> result(
new Structure(
"result"));
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());
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());
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());
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());
490 return result.release();
501 double get_double_value(BaseType *btp)
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();
511 throw BESSyntaxUserError(
string(
"Expected a constant value, but got ").append(btp->type_name())
512 .append(
" instead."), __FILE__, __LINE__);
518 StareSubsetArrayFunction::stare_subset_array_dap4_function(D4RValueList *args, DMR &dmr)
520 if (args->size() != 3) {
522 oss <<
"stare_subset_array(): Expected three arguments, but got " << args->size();
526 Array *dependent_var =
dynamic_cast<Array*
>(args->get_rvalue(0)->value(dmr));
528 throw BESSyntaxUserError(
"Expected an Array as teh first argument to stare_subset_array()", __FILE__, __LINE__);
530 BaseType *mask_val_var = args->get_rvalue(1)->value(dmr);
531 double mask_value = get_double_value(mask_val_var);
533 BaseType *raw_stare_indices = args->get_rvalue(2)->value(dmr);
535 unique_ptr<GeoFile> gf(
new GeoFile(dmr.filename()));
538 vector<STARE_ArrayIndexSpatialValue> dep_var_stare_indices;
539 gf->get_stare_indices(dependent_var->name(), dep_var_stare_indices);
542 vector<STARE_ArrayIndexSpatialValue> target_s_indices;
543 read_stare_indices_from_function_argument(raw_stare_indices, target_s_indices);
546 unique_ptr<Array> result(
dynamic_cast<Array*
>(dependent_var->ptr_duplicate()));
549 switch(dependent_var->var()->type()) {
551 build_masked_data<dods_int16>(dependent_var, dep_var_stare_indices, target_s_indices, mask_value, result);
554 case dods_float32_c: {
555 build_masked_data<dods_float32>(dependent_var, dep_var_stare_indices, target_s_indices, mask_value, result);
560 throw BESInternalError(
string(
"stare_subset_array() failed: Unsupported array element type (")
561 + dependent_var->var()->type_name() +
").", __FILE__, __LINE__);
564 return result.release();
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));
582 return index.CoverBoundingBoxFromLatLonDegrees(latlonbox, resolution);
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));
605 return index.CoverBoundingBoxFromLatLonDegrees(latlonbox, resolution);
609 StareBoxFunction::stare_box_dap4_function(libdap::D4RValueList *args, libdap::DMR &dmr)
611 STARE_SpatialIntervals sivs;
613 if (args->size() == 4) {
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));
620 sivs = stare_box_helper(point(tl_lat, tl_lon), point(br_lat, br_lon));
623 else if (args->size() >= 6 && (args->size() % 2) == 0) {
625 bool flip_flop =
false;
626 double lat_value = -90;
627 vector<point> points;
628 for (
auto &arg: *args) {
630 point pt(lat_value, get_double_value(arg->value(dmr)));
631 points.push_back(pt);
635 lat_value = get_double_value(arg->value(dmr));
640 sivs = stare_box_helper(points);
645 oss <<
"stare_box(): Expected four corner lat/lon values or a list of three or more points, but got "
646 << args->size() <<
" values.";
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());
654 return cover.release();
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.