42#include <InternalErr.h>
44#include <dods-datatypes.h>
47#include "GeoConstraint.h"
63 is_prefix(
const string &in) : s(in) {}
65 bool operator()(
const string &prefix) {
return s.find(prefix) == 0; }
78bool unit_or_name_match(set<string> units, set<string> names,
const string &var_units,
const string &var_name) {
79 return (units.find(var_units) != units.end() ||
80 find_if(names.begin(), names.end(), is_prefix(var_name)) != names.end());
98 return (left < 0 || right < 0) ? neg_pos : pos;
115void GeoConstraint::transform_constraint_to_pos_notation(
double &left,
double &right)
const {
135 for (
int i = 0; i < d_lon_length; ++i)
150 for (
int i = 0; i < d_lon_length; ++i)
155bool GeoConstraint::is_bounding_box_valid(
const double left,
const double top,
const double right,
156 const double bottom)
const {
157 if ((left < d_lon[0] && right < d_lon[0]) || (left > d_lon[d_lon_length - 1] && right > d_lon[d_lon_length - 1]))
160 if (d_latitude_sense == normal) {
162 if ((top > d_lat[0] && bottom > d_lat[0]) ||
163 (top < d_lat[d_lat_length - 1] && bottom < d_lat[d_lat_length - 1]))
166 if ((top < d_lat[0] && bottom < d_lat[0]) ||
167 (top > d_lat[d_lat_length - 1] && bottom > d_lat[d_lat_length - 1]))
185 int &longitude_index_right)
const {
189 double t_left = fmod(left, 360.0);
190 double t_right = fmod(right, 360.0);
197 int lon_origin_index = 0;
198 double smallest_lon = fmod(d_lon[0], 360.0);
199 while (i < d_lon_length) {
200 double curent_lon_value = fmod(d_lon[i], 360.0);
201 if (smallest_lon > curent_lon_value) {
202 smallest_lon = curent_lon_value;
203 lon_origin_index = i;
208 DBG2(cerr <<
"lon_origin_index: " << lon_origin_index << endl);
213 i = lon_origin_index;
214 while (fmod(d_lon[i], 360.0) < t_left) {
216 i = i % d_lon_length;
219 if (i == lon_origin_index)
220 throw Error(
"geogrid: Could not find an index for the longitude value '" + double_to_string(left) +
"'");
223 if (fmod(d_lon[i], 360.0) == t_left)
224 longitude_index_left = i;
226 longitude_index_left = (i - 1) > 0 ? i - 1 : 0;
228 DBG2(cerr <<
"longitude_index_left: " << longitude_index_left << endl);
232 int largest_lon_index = (lon_origin_index - 1 + d_lon_length) % d_lon_length;
233 i = largest_lon_index;
234 while (fmod(d_lon[i], 360.0) > t_right) {
236 i = (i == 0) ? d_lon_length - 1 : i - 1;
237 if (i == largest_lon_index)
238 throw Error(
"geogrid: Could not find an index for the longitude value '" + double_to_string(right) +
"'");
241 if (fmod(d_lon[i], 360.0) == t_right)
242 longitude_index_right = i;
244 longitude_index_right = (i + 1) < d_lon_length - 1 ? i + 1 : d_lon_length - 1;
246 DBG2(cerr <<
"longitude_index_right: " << longitude_index_right << endl);
262 int &latitude_index_bottom)
const {
265 if (sense == normal) {
267 while (i < d_lat_length - 1 && top < d_lat[i])
270 j = d_lat_length - 1;
271 while (j > 0 && bottom > d_lat[j])
275 latitude_index_top = i;
277 latitude_index_top = (i - 1) > 0 ? i - 1 : 0;
279 if (d_lat[j] == bottom)
280 latitude_index_bottom = j;
282 latitude_index_bottom = (j + 1) < d_lat_length - 1 ? j + 1 : d_lat_length - 1;
284 i = d_lat_length - 1;
285 while (i > 0 && d_lat[i] > top)
289 while (j < d_lat_length - 1 && d_lat[j] < bottom)
293 latitude_index_top = i;
295 latitude_index_top = (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1;
297 if (d_lat[j] == bottom)
298 latitude_index_bottom = j;
300 latitude_index_bottom = (j - 1) > 0 ? j - 1 : 0;
308 return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted;
313static void swap_vector_ends(
char *dest,
char *src,
int len,
int index,
int elem_sz) {
314 memcpy(dest, src + index * elem_sz, (len - index) * elem_sz);
316 memcpy(dest + (len - index) * elem_sz, src, index * elem_sz);
320static void transpose(std::vector<std::vector<T>> a, std::vector<std::vector<T>> b,
int width,
int height) {
321 for (
int i = 0; i < width; i++) {
322 for (
int j = 0; j < height; j++) {
336 double *tmp =
new double[length];
338 int i = 0, j = length - 1;
342 memcpy(src, tmp, length *
sizeof(
double));
347static int count_size_except_latitude_and_longitude(
Array &a) {
358void GeoConstraint::flip_latitude_within_array(
Array &a,
int lat_length,
int lon_length) {
361 d_array_data =
static_cast<char *
>(a.value());
362 d_array_data_size = a.
width(
true);
365 int size = count_size_except_latitude_and_longitude(a);
367 vector<char> tmp_data(d_array_data_size);
368 int array_elem_size = a.
var()->
width(
true);
369 int lat_lon_size = (d_array_data_size / size);
371 DBG(cerr <<
"lat, lon_length: " << lat_length <<
", " << lon_length << endl);
372 DBG(cerr <<
"size: " << size << endl);
373 DBG(cerr <<
"d_array_data_size: " << d_array_data_size << endl);
374 DBG(cerr <<
"array_elem_size: " << array_elem_size << endl);
375 DBG(cerr <<
"lat_lon_size: " << lat_lon_size << endl);
377 for (
int i = 0; i < size; ++i) {
379 int s_lat = lat_length - 1;
381 int lon_size = array_elem_size * lon_length;
382 int offset = i * lat_lon_size;
384 memcpy(tmp_data.data() + offset + (lat++ * lon_size), d_array_data + offset + (s_lat-- * lon_size),
388 memcpy(d_array_data, tmp_data.data(), d_array_data_size);
400 double *tmp_lon =
new double[d_lon_length];
402 swap_vector_ends((
char *)tmp_lon, (
char *)d_lon, d_lon_length, longitude_index_left,
sizeof(
double));
404 memcpy(d_lon, tmp_lon, d_lon_length *
sizeof(
double));
409static int count_dimensions_except_longitude(
Array &a) {
436 if (!is_longitude_rightmost())
437 throw Error(
"This grid does not have Longitude as its rightmost dimension, the geogrid()\ndoes not support "
438 "constraints that wrap around the edges of this type of grid.");
440 DBG(cerr <<
"Constraint for the left half: " << get_longitude_index_left() <<
", " << get_lon_length() - 1 << endl);
443 a.
add_constraint(lon_dim, get_longitude_index_left(), 1, get_lon_length() - 1);
449 int left_size = a.
width(
true);
450 char *left_data = (
char *)a.value();
458 DBG(cerr <<
"Constraint for the right half: " << 0 <<
", " << get_longitude_index_right() << endl);
465 char *right_data = (
char *)a.value();
466 int right_size = a.
width(
true);
469 d_array_data_size = left_size + right_size;
470 d_array_data =
new char[d_array_data_size];
475 int elem_size = a.
var()->
width(
true);
476 int left_row_size = (get_lon_length() - get_longitude_index_left()) * elem_size;
477 int right_row_size = (get_longitude_index_right() + 1) * elem_size;
478 int total_bytes_per_row = left_row_size + right_row_size;
480 DBG2(cerr <<
"elem_size: " << elem_size <<
"; left & right size: " << left_row_size <<
", " << right_row_size
485 int rows_to_copy = count_dimensions_except_longitude(a);
486 for (
int i = 0; i < rows_to_copy; ++i) {
487 DBG(cerr <<
"Copying " << i <<
"th row" << endl);
488 DBG(cerr <<
"left memcpy: " << *(
float *)(left_data + (left_row_size * i)) << endl);
490 memcpy(d_array_data + (total_bytes_per_row * i), left_data + (left_row_size * i), left_row_size);
492 DBG(cerr <<
"right memcpy: " << *(
float *)(right_data + (right_row_size * i)) << endl);
494 memcpy(d_array_data + (total_bytes_per_row * i) + left_row_size, right_data + (right_row_size * i),
505 : d_array_data(0), d_array_data_size(0), d_lat(0), d_lon(0), d_lat_length(0), d_lon_length(0),
506 d_latitude_index_top(0), d_latitude_index_bottom(0), d_longitude_index_left(0), d_longitude_index_right(0),
507 d_bounding_box_set(false), d_longitude_rightmost(false), d_longitude_notation(unknown_notation),
508 d_latitude_sense(unknown_sense) {
510 d_coards_lat_units.insert(
"degrees_north");
511 d_coards_lat_units.insert(
"degree_north");
512 d_coards_lat_units.insert(
"degree_N");
513 d_coards_lat_units.insert(
"degrees_N");
515 d_coards_lon_units.insert(
"degrees_east");
516 d_coards_lon_units.insert(
"degree_east");
517 d_coards_lon_units.insert(
"degrees_E");
518 d_coards_lon_units.insert(
"degree_E");
520 d_lat_names.insert(
"COADSY");
521 d_lat_names.insert(
"lat");
522 d_lat_names.insert(
"Lat");
523 d_lat_names.insert(
"LAT");
525 d_lon_names.insert(
"COADSX");
526 d_lon_names.insert(
"lon");
527 d_lon_names.insert(
"Lon");
528 d_lon_names.insert(
"LON");
545 if (d_bounding_box_set)
546 throw Error(
"It is not possible to register more than one geographical constraint on a variable.");
551 if (d_latitude_sense == inverted)
552 throw Error(
"geogrid() does not currently work with inverted data (data where the north pole is at a negative latitude value).");
559 if (d_longitude_notation == neg_pos)
560 transform_constraint_to_pos_notation(left, right);
567 if (longitude_notation == neg_pos)
570 if (!is_bounding_box_valid(left, top, right, bottom))
571 throw Error(
"The bounding box does not intersect any data within this Grid or Array. The\ngeographical extent "
572 "of these data are from latitude " +
573 double_to_string(d_lat[0]) +
" to " + double_to_string(d_lat[d_lat_length - 1]) +
574 "\nand longitude " + double_to_string(d_lon[0]) +
" to " +
575 double_to_string(d_lon[d_lon_length - 1]) +
" while the bounding box provided was latitude " +
576 double_to_string(top) +
" to " + double_to_string(bottom) +
"\nand longitude " +
577 double_to_string(left) +
" to " + double_to_string(right));
588 DBG(cerr <<
"Bounding box (tlbr): " << d_latitude_index_top <<
", " << d_longitude_index_left <<
", "
589 << d_latitude_index_bottom <<
", " << d_longitude_index_right << endl);
591 d_bounding_box_set =
true;
GeoConstraint()
Initialize GeoConstraint.
virtual void transform_longitude_to_pos_notation()
virtual void reorder_longitude_map(int longitude_index_left)
virtual void transform_longitude_to_neg_pos_notation()
virtual void reorder_data_longitude_axis(libdap::Array &a, libdap::Array::Dim_iter lon_dim)
void find_latitude_indeces(double top, double bottom, LatitudeSense sense, int &latitude_index_top, int &latitude_index_bottom) const
Notation categorize_notation(const double left, const double right) const
void find_longitude_indeces(double left, double right, int &longitude_index_left, int &longitude_index_right) const
virtual void transpose_vector(double *src, const int length)
virtual LatitudeSense categorize_latitude() const
void set_bounding_box(double top, double left, double bottom, double right)
A multidimensional array of identical data types.
virtual void add_constraint(Dim_iter i, int start, int stride, int stop)
Adds a constraint to an Array dimension.
void print_val(ostream &out, string space="", bool print_decl_p=true) override
Prints the value of the variable.
std::vector< dimension >::iterator Dim_iter
virtual int dimension_size(Dim_iter i, bool constrained=false)
Returns the size of the dimension.
virtual bool read()
Read data into a local buffer.
virtual unsigned int width(bool constrained=false) const
How many bytes does this variable use Return the number of bytes of storage this variable uses....
A class for error processing.
unsigned int width(bool constrained=false) const override
Returns the width of the data, in bytes.
void set_read_p(bool state) override
Indicates that the data is ready to send.
BaseType * var(const string &name="", bool exact_match=true, btp_stack *s=nullptr) override
top level DAP object to house generic methods