42 #include <InternalErr.h>
43 #include <dods-datatypes.h>
62 is_prefix(
const string & in): s(in)
65 bool operator()(
const string & prefix)
67 return s.find(prefix) == 0;
83 const string & var_units,
const string & var_name)
85 return (units.find(var_units) != units.end()
86 || find_if(names.begin(), names.end(),
87 is_prefix(var_name)) != names.end());
104 GeoConstraint::Notation
105 GeoConstraint::categorize_notation(
const double left,
106 const double right)
const
108 return (left < 0 || right < 0) ? neg_pos : pos;
126 GeoConstraint::transform_constraint_to_pos_notation(
double &left,
144 void GeoConstraint::transform_longitude_to_pos_notation()
149 for (
int i = 0; i < d_lon_length; ++i)
163 void GeoConstraint::transform_longitude_to_neg_pos_notation()
165 for (
int i = 0; i < d_lon_length; ++i)
170 bool GeoConstraint::is_bounding_box_valid(
const double left,
const double top,
171 const double right,
const double bottom)
const
173 if ((left < d_lon[0] && right < d_lon[0])
174 || (left > d_lon[d_lon_length-1] && right > d_lon[d_lon_length-1]))
177 if (d_latitude_sense == normal) {
179 if ((top > d_lat[0] && bottom > d_lat[0])
180 || (top < d_lat[d_lat_length-1] && bottom < d_lat[d_lat_length-1]))
184 if ((top < d_lat[0] && bottom < d_lat[0])
185 || (top > d_lat[d_lat_length-1] && bottom > d_lat[d_lat_length-1]))
202 void GeoConstraint::find_longitude_indeces(
double left,
double right,
203 int &longitude_index_left,
int &longitude_index_right)
const
208 double t_left = fmod(left, 360.0);
209 double t_right = fmod(right, 360.0);
216 int lon_origin_index = 0;
217 double smallest_lon = fmod(d_lon[0], 360.0);
218 while (i < d_lon_length) {
219 double curent_lon_value = fmod(d_lon[i], 360.0);
220 if (smallest_lon > curent_lon_value) {
221 smallest_lon = curent_lon_value;
222 lon_origin_index = i;
227 DBG2(cerr <<
"lon_origin_index: " << lon_origin_index << endl);
232 i = lon_origin_index;
233 while (fmod(d_lon[i], 360.0) < t_left) {
235 i = i % d_lon_length;
238 if (i == lon_origin_index)
239 throw Error(
"geogrid: Could not find an index for the longitude value '" + double_to_string(left) +
"'");
242 if (fmod(d_lon[i], 360.0) == t_left)
243 longitude_index_left = i;
245 longitude_index_left = (i - 1) > 0 ? i - 1 : 0;
247 DBG2(cerr <<
"longitude_index_left: " << longitude_index_left << endl);
251 int largest_lon_index = (lon_origin_index - 1 + d_lon_length) % d_lon_length;
252 i = largest_lon_index;
253 while (fmod(d_lon[i], 360.0) > t_right) {
255 i = (i == 0) ? d_lon_length - 1 : i - 1;
256 if (i == largest_lon_index)
257 throw Error(
"geogrid: Could not find an index for the longitude value '" + double_to_string(right) +
"'");
260 if (fmod(d_lon[i], 360.0) == t_right)
261 longitude_index_right = i;
263 longitude_index_right = (i + 1) < d_lon_length - 1 ? i + 1 : d_lon_length - 1;
265 DBG2(cerr <<
"longitude_index_right: " << longitude_index_right << endl);
280 void GeoConstraint::find_latitude_indeces(
double top,
double bottom,
282 int &latitude_index_top,
283 int &latitude_index_bottom)
const
287 if (sense == normal) {
289 while (i < d_lat_length - 1 && top < d_lat[i])
292 j = d_lat_length - 1;
293 while (j > 0 && bottom > d_lat[j])
297 latitude_index_top = i;
299 latitude_index_top = (i - 1) > 0 ? i - 1 : 0;
301 if (d_lat[j] == bottom)
302 latitude_index_bottom = j;
304 latitude_index_bottom =
305 (j + 1) < d_lat_length - 1 ? j + 1 : d_lat_length - 1;
308 i = d_lat_length - 1;
309 while (i > 0 && d_lat[i] > top)
313 while (j < d_lat_length - 1 && d_lat[j] < bottom)
317 latitude_index_top = i;
319 latitude_index_top = (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1;
321 if (d_lat[j] == bottom)
322 latitude_index_bottom = j;
324 latitude_index_bottom = (j - 1) > 0 ? j - 1 : 0;
333 return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted;
339 swap_vector_ends(
char *dest,
char *src,
int len,
int index,
int elem_sz)
341 memcpy(dest, src + index * elem_sz, (len - index) * elem_sz);
343 memcpy(dest + (len - index) * elem_sz, src, index * elem_sz);
347 static void transpose(std::vector<std::vector<T> > a,
348 std::vector<std::vector<T> > b,
int width,
int height)
350 for (
int i = 0; i < width; i++) {
351 for (
int j = 0; j < height; j++) {
364 void GeoConstraint::transpose_vector(
double *src,
const int length)
366 double *tmp =
new double[length];
368 int i = 0, j = length-1;
372 memcpy(src, tmp,length *
sizeof(
double));
378 count_size_except_latitude_and_longitude(Array &a)
380 if (a.dim_end() - a.dim_begin() <= 2)
384 for (Array::Dim_iter i = a.dim_begin(); (i + 2) != a.dim_end(); ++i)
385 size *= a.dimension_size(i,
true);
390 void GeoConstraint::flip_latitude_within_array(Array &a,
int lat_length,
395 d_array_data =
static_cast<char*
>(a.value());
396 d_array_data_size = a.width(
true);
399 int size = count_size_except_latitude_and_longitude(a);
401 vector<char> tmp_data(d_array_data_size);
402 int array_elem_size = a.var()->width(
true);
403 int lat_lon_size = (d_array_data_size / size);
405 DBG(cerr <<
"lat, lon_length: " << lat_length <<
", " << lon_length << endl);
406 DBG(cerr <<
"size: " << size << endl);
407 DBG(cerr <<
"d_array_data_size: " << d_array_data_size << endl);
408 DBG(cerr <<
"array_elem_size: " << array_elem_size<< endl);
409 DBG(cerr <<
"lat_lon_size: " << lat_lon_size<< endl);
411 for (
int i = 0; i < size; ++i) {
413 int s_lat = lat_length - 1;
415 int lon_size = array_elem_size * lon_length;
416 int offset = i * lat_lon_size;
418 memcpy(&tmp_data[0] + offset + (lat++ * lon_size),
419 d_array_data + offset + (s_lat-- * lon_size),
423 memcpy(d_array_data, &tmp_data[0], d_array_data_size);
434 void GeoConstraint::reorder_longitude_map(
int longitude_index_left)
436 double *tmp_lon =
new double[d_lon_length];
438 swap_vector_ends((
char *) tmp_lon, (
char *) d_lon, d_lon_length,
439 longitude_index_left,
sizeof(
double));
441 memcpy(d_lon, tmp_lon, d_lon_length *
sizeof(
double));
447 count_dimensions_except_longitude(Array &a)
450 for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i)
451 size *= a.dimension_size(i,
true);
473 void GeoConstraint::reorder_data_longitude_axis(Array &a, Array::Dim_iter lon_dim)
476 if (!is_longitude_rightmost())
477 throw Error(
"This grid does not have Longitude as its rightmost dimension, the geogrid()\ndoes not support constraints that wrap around the edges of this type of grid.");
479 DBG(cerr <<
"Constraint for the left half: " << get_longitude_index_left()
480 <<
", " << get_lon_length() - 1 << endl);
483 a.add_constraint(lon_dim, get_longitude_index_left(), 1,
484 get_lon_length() - 1);
487 DBG2(a.print_val(stderr));
490 int left_size = a.width(
true);
491 char *left_data = (
char*)a.value();
499 DBG(cerr <<
"Constraint for the right half: " << 0
500 <<
", " << get_longitude_index_right() << endl);
502 a.add_constraint(lon_dim, 0, 1, get_longitude_index_right());
505 DBG2(a.print_val(stderr));
507 char *right_data = (
char*)a.value();
508 int right_size = a.width(
true);
511 d_array_data_size = left_size + right_size;
512 d_array_data =
new char[d_array_data_size];
517 int elem_size = a.var()->width(
true);
518 int left_row_size = (get_lon_length() - get_longitude_index_left()) * elem_size;
519 int right_row_size = (get_longitude_index_right() + 1) * elem_size;
520 int total_bytes_per_row = left_row_size + right_row_size;
522 DBG2(cerr <<
"elem_size: " << elem_size <<
"; left & right size: "
523 << left_row_size <<
", " << right_row_size << endl);
527 int rows_to_copy = count_dimensions_except_longitude(a);
528 for (
int i = 0; i < rows_to_copy; ++i) {
529 DBG(cerr <<
"Copying " << i <<
"th row" << endl);
530 DBG(cerr <<
"left memcpy: " << *(
float *)(left_data + (left_row_size * i)) << endl);
532 memcpy(d_array_data + (total_bytes_per_row * i),
533 left_data + (left_row_size * i),
536 DBG(cerr <<
"right memcpy: " << *(
float *)(right_data + (right_row_size * i)) << endl);
538 memcpy(d_array_data + (total_bytes_per_row * i) + left_row_size,
539 right_data + (right_row_size * i),
549 GeoConstraint::GeoConstraint()
550 : d_array_data(0), d_array_data_size(0),
552 d_bounding_box_set(false),
553 d_longitude_rightmost(false),
554 d_longitude_notation(unknown_notation),
555 d_latitude_sense(unknown_sense)
558 d_coards_lat_units.insert(
"degrees_north");
559 d_coards_lat_units.insert(
"degree_north");
560 d_coards_lat_units.insert(
"degree_N");
561 d_coards_lat_units.insert(
"degrees_N");
563 d_coards_lon_units.insert(
"degrees_east");
564 d_coards_lon_units.insert(
"degree_east");
565 d_coards_lon_units.insert(
"degrees_E");
566 d_coards_lon_units.insert(
"degree_E");
568 d_lat_names.insert(
"COADSY");
569 d_lat_names.insert(
"lat");
570 d_lat_names.insert(
"Lat");
571 d_lat_names.insert(
"LAT");
573 d_lon_names.insert(
"COADSX");
574 d_lon_names.insert(
"lon");
575 d_lon_names.insert(
"Lon");
576 d_lon_names.insert(
"LON");
590 double bottom,
double right)
595 if (d_bounding_box_set)
596 throw Error(
"It is not possible to register more than one geographical constraint on a variable.");
602 throw Error(
"geogrid() does not currently work with inverted data (data where the north pole is at a negative latitude value).");
609 if (d_longitude_notation ==
neg_pos)
618 if (longitude_notation ==
neg_pos)
622 throw Error(
"The bounding box does not intersect any data within this Grid or Array. The\ngeographical extent of these data are from latitude "
623 + double_to_string(d_lat[0]) +
" to "
624 + double_to_string(d_lat[d_lat_length-1])
625 +
"\nand longitude " + double_to_string(d_lon[0])
626 +
" to " + double_to_string(d_lon[d_lon_length-1])
627 +
" while the bounding box provided was latitude "
628 + double_to_string(top) +
" to "
629 + double_to_string(bottom) +
"\nand longitude "
630 + double_to_string(left) +
" to "
631 + double_to_string(right));
638 d_latitude_index_top, d_latitude_index_bottom);
643 d_longitude_index_right);
645 DBG(cerr <<
"Bounding box (tlbr): " << d_latitude_index_top <<
", "
646 << d_longitude_index_left <<
", "
647 << d_latitude_index_bottom <<
", "
648 << d_longitude_index_right << endl);
650 d_bounding_box_set =
true;
void transform_constraint_to_pos_notation(double &left, double &right) const
void find_longitude_indeces(double left, double right, int &longitude_index_left, int &longitude_index_right) const
Scan from the left to the right, and the right to the left, looking for the left and right bounding b...
void set_bounding_box(double top, double left, double bottom, double right)
Set the bounding box for this constraint.
void find_latitude_indeces(double top, double bottom, LatitudeSense sense, int &latitude_index_top, int &latitude_index_bottom) const
Scan from the top to the bottom, and the bottom to the top, looking for the top and bottom bounding b...
bool unit_or_name_match(set< string > units, set< string > names, const string &var_units, const string &var_name)
Look in the containers which hold the units attributes and variable name prefixes which are considere...
virtual void transform_longitude_to_pos_notation()
Given that the Grid has a longitude map that uses the 'neg_pos' notation, transform it to the 'pos' n...
virtual LatitudeSense categorize_latitude() const
Take a look at the latitude vector values and record whether the world is normal or upside down...
LatitudeSense
Most of the time, latitude starts at the top of an array with positive values and ends up at the bott...
virtual bool is_bounding_box_valid(const double left, const double top, const double right, const double bottom) const
Notation
The longitude extents of the constraint bounding box can be expressed two ways: using a 0/359 notatio...
Notation categorize_notation(const double left, const double right) const
A private method that determines if the longitude part of the bounding box uses 0/359 or -180/179 not...