OPeNDAP Hyrax Back End Server (BES)  Updated for version 3.8.3
GeoConstraint.cc
Go to the documentation of this file.
1 
2 // -*- mode: c++; c-basic-offset:4 -*-
3 
4 // This file is part of libdap, A C++ implementation of the OPeNDAP Data
5 // Access Protocol.
6 
7 // Copyright (c) 2006 OPeNDAP, Inc.
8 // Author: James Gallagher <jgallagher@opendap.org>
9 //
10 // This library is free software; you can redistribute it and/or
11 // modify it under the terms of the GNU Lesser General Public
12 // License as published by the Free Software Foundation; either
13 // version 2.1 of the License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 //
24 // You can contact OPeNDAP, Inc. at PO Box 112, Saunderstown, RI. 02874-0112.
25 
26 // The Grid Selection Expression Clause class.
27 
28 
29 #include "config.h"
30 
31 #include <cstring>
32 #include <cmath>
33 #include <iostream>
34 #include <sstream>
35 #include <algorithm> // for find_if
36 
37 //#define DODS_DEBUG
38 //#define DODS_DEBUG2
39 
40 #include <Float64.h>
41 #include <Error.h>
42 #include <InternalErr.h>
43 #include <dods-datatypes.h>
44 #include <util.h>
45 #include <debug.h>
46 
47 #include "GeoConstraint.h"
48 
49 using namespace std;
50 
51 namespace libdap {
52 
57 class is_prefix
58 {
59 private:
60  string s;
61 public:
62  is_prefix(const string & in): s(in)
63  {}
64 
65  bool operator()(const string & prefix)
66  {
67  return s.find(prefix) == 0;
68  }
69 };
70 
81 bool
82 unit_or_name_match(set < string > units, set < string > names,
83  const string & var_units, const string & var_name)
84 {
85  return (units.find(var_units) != units.end()
86  || find_if(names.begin(), names.end(),
87  is_prefix(var_name)) != names.end());
88 }
89 
104 GeoConstraint::Notation
105 GeoConstraint::categorize_notation(const double left,
106  const double right) const
107 {
108  return (left < 0 || right < 0) ? neg_pos : pos;
109 }
110 
111 /* A private method to translate the longitude constraint from -180/179
112  notation to 0/359 notation.
113 
114  About the two notations:
115  0 180 360
116  |---------------------------|-------------------------|
117  0 180,-180 0
118 
119  so in the neg-pos notation (using the name I give it in this class) both
120  180 and -180 are the same longitude. And in the pos notation 0 and 360 are
121  the same.
122 
123  @param left Value-result parameter; the left side of the bounding box
124  @parm right Value-result parameter; the right side of the bounding box */
125 void
126 GeoConstraint::transform_constraint_to_pos_notation(double &left,
127  double &right) const
128 {
129  if (left < 0)
130  left += 360 ;
131 
132  if (right < 0)
133  right += 360;
134 }
135 
144 void GeoConstraint::transform_longitude_to_pos_notation()
145 {
146  // Assume earlier logic is correct (since the test is expensive)
147  // for each value, add 180
148  // Longitude could be represented using any of the numeric types
149  for (int i = 0; i < d_lon_length; ++i)
150  if (d_lon[i] < 0)
151  d_lon[i] += 360;
152 }
153 
163 void GeoConstraint::transform_longitude_to_neg_pos_notation()
164 {
165  for (int i = 0; i < d_lon_length; ++i)
166  if (d_lon[i] > 180)
167  d_lon[i] -= 360;
168 }
169 
170 bool GeoConstraint::is_bounding_box_valid(const double left, const double top,
171  const double right, const double bottom) const
172 {
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]))
175  return false;
176 
177  if (d_latitude_sense == normal) {
178  // When sense is normal, the largest values are at the origin.
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]))
181  return false;
182  }
183  else {
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]))
186  return false;
187  }
188 
189  return true;
190 }
191 
202 void GeoConstraint::find_longitude_indeces(double left, double right,
203  int &longitude_index_left, int &longitude_index_right) const
204 {
205  // Use all values 'modulo 360' to take into account the cases where the
206  // constraint and/or the longitude vector are given using values greater
207  // than 360 (i.e., when 380 is used to mean 20).
208  double t_left = fmod(left, 360.0);
209  double t_right = fmod(right, 360.0);
210 
211  // Find the place where 'longitude starts.' That is, what value of the
212  // index 'i' corresponds to the smallest value of d_lon. Why we do this:
213  // Some data sources use offset longitude axes so that the 'seam' is
214  // shifted to a place other than the date line.
215  int i = 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;
223  }
224  ++i;
225  }
226 
227  DBG2(cerr << "lon_origin_index: " << lon_origin_index << endl);
228 
229  // Scan from the index of the smallest value looking for the place where
230  // the value is greater than or equal to the left most point of the bounding
231  // box.
232  i = lon_origin_index;
233  while (fmod(d_lon[i], 360.0) < t_left) {
234  ++i;
235  i = i % d_lon_length;
236 
237  // If we cycle completely through all the values/indices, throw
238  if (i == lon_origin_index)
239  throw Error("geogrid: Could not find an index for the longitude value '" + double_to_string(left) + "'");
240  }
241 
242  if (fmod(d_lon[i], 360.0) == t_left)
243  longitude_index_left = i;
244  else
245  longitude_index_left = (i - 1) > 0 ? i - 1 : 0;
246 
247  DBG2(cerr << "longitude_index_left: " << longitude_index_left << endl);
248 
249  // Assume the vector is circular --> the largest value is next to the
250  // smallest.
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) {
254  // This is like modulus but for 'counting down'
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) + "'");
258  }
259 
260  if (fmod(d_lon[i], 360.0) == t_right)
261  longitude_index_right = i;
262  else
263  longitude_index_right = (i + 1) < d_lon_length - 1 ? i + 1 : d_lon_length - 1;
264 
265  DBG2(cerr << "longitude_index_right: " << longitude_index_right << endl);
266 }
267 
280 void GeoConstraint::find_latitude_indeces(double top, double bottom,
281  LatitudeSense sense,
282  int &latitude_index_top,
283  int &latitude_index_bottom) const
284 {
285  int i, j;
286 
287  if (sense == normal) {
288  i = 0;
289  while (i < d_lat_length - 1 && top < d_lat[i])
290  ++i;
291 
292  j = d_lat_length - 1;
293  while (j > 0 && bottom > d_lat[j])
294  --j;
295 
296  if (d_lat[i] == top)
297  latitude_index_top = i;
298  else
299  latitude_index_top = (i - 1) > 0 ? i - 1 : 0;
300 
301  if (d_lat[j] == bottom)
302  latitude_index_bottom = j;
303  else
304  latitude_index_bottom =
305  (j + 1) < d_lat_length - 1 ? j + 1 : d_lat_length - 1;
306  }
307  else {
308  i = d_lat_length - 1;
309  while (i > 0 && d_lat[i] > top)
310  --i;
311 
312  j = 0;
313  while (j < d_lat_length - 1 && d_lat[j] < bottom)
314  ++j;
315 
316  if (d_lat[i] == top)
317  latitude_index_top = i;
318  else
319  latitude_index_top = (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1;
320 
321  if (d_lat[j] == bottom)
322  latitude_index_bottom = j;
323  else
324  latitude_index_bottom = (j - 1) > 0 ? j - 1 : 0;
325  }
326 }
327 
331 GeoConstraint::LatitudeSense GeoConstraint::categorize_latitude() const
332 {
333  return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted;
334 }
335 
336 // Use 'index' as the pivot point. Move the points behind index to the front of
337 // the vector and those points in front of and at index to the rear.
338 static void
339 swap_vector_ends(char *dest, char *src, int len, int index, int elem_sz)
340 {
341  memcpy(dest, src + index * elem_sz, (len - index) * elem_sz);
342 
343  memcpy(dest + (len - index) * elem_sz, src, index * elem_sz);
344 }
345 
346 template<class T>
347 static void transpose(std::vector<std::vector<T> > a,
348  std::vector<std::vector<T> > b, int width, int height)
349 {
350  for (int i = 0; i < width; i++) {
351  for (int j = 0; j < height; j++) {
352  b[j][i] = a[i][j];
353  }
354  }
355 }
356 
364 void GeoConstraint::transpose_vector(double *src, const int length)
365 {
366  double *tmp = new double[length];
367 
368  int i = 0, j = length-1;
369  while (i < length)
370  tmp[j--] = src[i++];
371 
372  memcpy(src, tmp,length * sizeof(double));
373 
374  delete[] tmp;
375 }
376 
377 static int
378 count_size_except_latitude_and_longitude(Array &a)
379 {
380  if (a.dim_end() - a.dim_begin() <= 2) // < 2 is really an error...
381  return 1;
382 
383  int size = 1;
384  for (Array::Dim_iter i = a.dim_begin(); (i + 2) != a.dim_end(); ++i)
385  size *= a.dimension_size(i, true);
386 
387  return size;
388 }
389 
390 void GeoConstraint::flip_latitude_within_array(Array &a, int lat_length,
391  int lon_length)
392 {
393  if (!d_array_data) {
394  a.read();
395  d_array_data = static_cast<char*>(a.value());
396  d_array_data_size = a.width(true); // Bytes not elements
397  }
398 
399  int size = count_size_except_latitude_and_longitude(a);
400  // char *tmp_data = new char[d_array_data_size];
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);
404 
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);
410 
411  for (int i = 0; i < size; ++i) {
412  int lat = 0;
413  int s_lat = lat_length - 1;
414  // lon_length is the element size; memcpy() needs the number of bytes
415  int lon_size = array_elem_size * lon_length;
416  int offset = i * lat_lon_size;
417  while (s_lat > -1)
418  memcpy(&tmp_data[0] + offset + (lat++ * lon_size),
419  d_array_data + offset + (s_lat-- * lon_size),
420  lon_size);
421  }
422 
423  memcpy(d_array_data, &tmp_data[0], d_array_data_size);
424 }
425 
434 void GeoConstraint::reorder_longitude_map(int longitude_index_left)
435 {
436  double *tmp_lon = new double[d_lon_length];
437 
438  swap_vector_ends((char *) tmp_lon, (char *) d_lon, d_lon_length,
439  longitude_index_left, sizeof(double));
440 
441  memcpy(d_lon, tmp_lon, d_lon_length * sizeof(double));
442 
443  delete[]tmp_lon;
444 }
445 
446 static int
447 count_dimensions_except_longitude(Array &a)
448 {
449  int size = 1;
450  for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i)
451  size *= a.dimension_size(i, true);
452 
453  return size;
454 }
455 
473 void GeoConstraint::reorder_data_longitude_axis(Array &a, Array::Dim_iter lon_dim)
474 {
475 
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.");
478 
479  DBG(cerr << "Constraint for the left half: " << get_longitude_index_left()
480  << ", " << get_lon_length() - 1 << endl);
481 
482  // Build a constraint for the left part and get those values
483  a.add_constraint(lon_dim, get_longitude_index_left(), 1,
484  get_lon_length() - 1);
485  a.set_read_p(false);
486  a.read();
487  DBG2(a.print_val(stderr));
488 
489  // Save the left-hand data to local storage
490  int left_size = a.width(true); // width() == length() * element size
491  char *left_data = (char*)a.value(); // value() allocates and copies
492 
493  // Build a constraint for the 'right' part, which goes from the left edge
494  // of the array to the right index and read those data.
495  // (Don't call a.clear_constraint() since that will clear the constraint on
496  // all the dimensions while add_constraint() will replace a constraint on
497  // the given dimension).
498 
499  DBG(cerr << "Constraint for the right half: " << 0
500  << ", " << get_longitude_index_right() << endl);
501 
502  a.add_constraint(lon_dim, 0, 1, get_longitude_index_right());
503  a.set_read_p(false);
504  a.read();
505  DBG2(a.print_val(stderr));
506 
507  char *right_data = (char*)a.value();
508  int right_size = a.width(true);
509 
510  // Make one big lump O'data
511  d_array_data_size = left_size + right_size;
512  d_array_data = new char[d_array_data_size];
513 
514  // Assume COARDS conventions are being followed: lon varies fastest.
515  // These *_elements variables are actually elements * bytes/element since
516  // memcpy() uses bytes.
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;
521 
522  DBG2(cerr << "elem_size: " << elem_size << "; left & right size: "
523  << left_row_size << ", " << right_row_size << endl);
524 
525  // This will work for any number of dimension so long as longitude is the
526  // right-most array dimension.
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);
531 
532  memcpy(d_array_data + (total_bytes_per_row * i),
533  left_data + (left_row_size * i),
534  left_row_size);
535 
536  DBG(cerr << "right memcpy: " << *(float *)(right_data + (right_row_size * i)) << endl);
537 
538  memcpy(d_array_data + (total_bytes_per_row * i) + left_row_size,
539  right_data + (right_row_size * i),
540  right_row_size);
541  }
542 
543  delete[]left_data;
544  delete[]right_data;
545 }
546 
549 GeoConstraint::GeoConstraint()
550  : d_array_data(0), d_array_data_size(0),
551  d_lat(0), d_lon(0),
552  d_bounding_box_set(false),
553  d_longitude_rightmost(false),
554  d_longitude_notation(unknown_notation),
555  d_latitude_sense(unknown_sense)
556 {
557  // Build sets of attribute values for easy searching. Maybe overkill???
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");
562 
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");
567 
568  d_lat_names.insert("COADSY");
569  d_lat_names.insert("lat");
570  d_lat_names.insert("Lat");
571  d_lat_names.insert("LAT");
572 
573  d_lon_names.insert("COADSX");
574  d_lon_names.insert("lon");
575  d_lon_names.insert("Lon");
576  d_lon_names.insert("LON");
577 }
578 
589 void GeoConstraint::set_bounding_box(double top, double left,
590  double bottom, double right)
591 {
592  // Ensure this method is called only once. What about pthreads?
593  // The method Array::reset_constraint() might make this so it could be
594  // called more than once. jhrg 8/30/06
595  if (d_bounding_box_set)
596  throw Error("It is not possible to register more than one geographical constraint on a variable.");
597 
598  // Record the 'sense' of the latitude for use here and later on.
599  d_latitude_sense = categorize_latitude();
600 #if 0
601  if (d_latitude_sense == inverted)
602  throw Error("geogrid() does not currently work with inverted data (data where the north pole is at a negative latitude value).");
603 #endif
604 
605  // Categorize the notation used by the bounding box (0/359 or -180/179).
606  d_longitude_notation = categorize_notation(left, right);
607 
608  // If the notation uses -180/179, transform the request to 0/359 notation.
609  if (d_longitude_notation == neg_pos)
611 
612  // If the grid uses -180/179, transform it to 0/359 as well. This will make
613  // subsequent logic easier and adds only a few extra operations, even with
614  // large maps.
615  Notation longitude_notation =
616  categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
617 
618  if (longitude_notation == neg_pos)
620 
621  if (!is_bounding_box_valid(left, top, right, bottom))
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));
632 
633  // This is simpler than the longitude case because there's no need to
634  // test for several notations, no need to accommodate them in the return,
635  // no modulo arithmetic for the axis and no need to account for a
636  // constraint with two disconnected parts to be joined.
637  find_latitude_indeces(top, bottom, d_latitude_sense,
638  d_latitude_index_top, d_latitude_index_bottom);
639 
640 
641  // Find the longitude map indexes that correspond to the bounding box.
642  find_longitude_indeces(left, right, d_longitude_index_left,
643  d_longitude_index_right);
644 
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);
649 
650  d_bounding_box_set = true;
651 }
652 
653 } // namespace libdap
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...