libdap++  Updated for version 3.8.2
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 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 static char id[] not_used =
32  { "$Id: GeoConstraint.cc 24370 2011-03-28 16:21:32Z jimg $"
33  };
34 
35 #include <cstring>
36 #include <cmath>
37 #include <iostream>
38 #include <sstream>
39 #include <algorithm> // for find_if
40 
41 //#define DODS_DEBUG
42 //#define DODS_DEBUG2
43 
44 #include "debug.h"
45 #include "dods-datatypes.h"
46 #include "GeoConstraint.h"
47 #include "Float64.h"
48 
49 #include "Error.h"
50 #include "InternalErr.h"
51 #include "ce_functions.h"
52 #include "util.h"
53 
54 using namespace std;
55 
56 namespace libdap {
57 
62 class is_prefix
63 {
64 private:
65  string s;
66 public:
67  is_prefix(const string & in): s(in)
68  {}
69 
70  bool operator()(const string & prefix)
71  {
72  return s.find(prefix) == 0;
73  }
74 };
75 
86 bool
87 unit_or_name_match(set < string > units, set < string > names,
88  const string & var_units, const string & var_name)
89 {
90  return (units.find(var_units) != units.end()
91  || find_if(names.begin(), names.end(),
92  is_prefix(var_name)) != names.end());
93 }
94 
109 GeoConstraint::Notation
110 GeoConstraint::categorize_notation(const double left,
111  const double right) const
112 {
113  return (left < 0 || right < 0) ? neg_pos : pos;
114 }
115 
116 /* A private method to translate the longitude constraint from -180/179
117  notation to 0/359 notation.
118 
119  About the two notations:
120  0 180 360
121  |---------------------------|-------------------------|
122  0 180,-180 0
123 
124  so in the neg-pos notation (using the name I give it in this class) both
125  180 and -180 are the same longitude. And in the pos notation 0 and 360 are
126  the same.
127 
128  @param left Value-result parameter; the left side of the bounding box
129  @parm right Value-result parameter; the right side of the bounding box */
130 void
131 GeoConstraint::transform_constraint_to_pos_notation(double &left,
132  double &right) const
133 {
134  if (left < 0)
135  left += 360 ;
136 
137  if (right < 0)
138  right += 360;
139 }
140 
149 void GeoConstraint::transform_longitude_to_pos_notation()
150 {
151  // Assume earlier logic is correct (since the test is expensive)
152  // for each value, add 180
153  // Longitude could be represented using any of the numeric types
154  for (int i = 0; i < d_lon_length; ++i)
155  if (d_lon[i] < 0)
156  d_lon[i] += 360;
157 }
158 
168 void GeoConstraint::transform_longitude_to_neg_pos_notation()
169 {
170  for (int i = 0; i < d_lon_length; ++i)
171  if (d_lon[i] > 180)
172  d_lon[i] -= 360;
173 }
174 
175 bool GeoConstraint::is_bounding_box_valid(const double left, const double top,
176  const double right, const double bottom) const
177 {
178  if ((left < d_lon[0] && right < d_lon[0])
179  || (left > d_lon[d_lon_length-1] && right > d_lon[d_lon_length-1]))
180  return false;
181 
182  if (d_latitude_sense == normal) {
183  // When sense is normal, the largest values are at the origin.
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  else {
189  if ((top < d_lat[0] && bottom < d_lat[0])
190  || (top > d_lat[d_lat_length-1] && bottom > d_lat[d_lat_length-1]))
191  return false;
192  }
193 
194  return true;
195 }
196 
207 void GeoConstraint::find_longitude_indeces(double left, double right,
208  int &longitude_index_left, int &longitude_index_right) const
209 {
210  // Use all values 'modulo 360' to take into account the cases where the
211  // constraint and/or the longitude vector are given using values greater
212  // than 360 (i.e., when 380 is used to mean 20).
213  double t_left = fmod(left, 360.0);
214  double t_right = fmod(right, 360.0);
215 
216  // Find the place where 'longitude starts.' That is, what value of the
217  // index 'i' corresponds to the smallest value of d_lon. Why we do this:
218  // Some data sources use offset longitude axes so that the 'seam' is
219  // shifted to a place other than the date line.
220  int i = 0;
221  int lon_origin_index = 0;
222  double smallest_lon = fmod(d_lon[0], 360.0);
223  while (i < d_lon_length) {
224  double curent_lon_value = fmod(d_lon[i], 360.0);
225  if (smallest_lon > curent_lon_value) {
226  smallest_lon = curent_lon_value;
227  lon_origin_index = i;
228  }
229  ++i;
230  }
231 
232  DBG2(cerr << "lon_origin_index: " << lon_origin_index << endl);
233 
234  // Scan from the index of the smallest value looking for the place where
235  // the value is greater than or equal to the left most point of the bounding
236  // box.
237  i = lon_origin_index;
238  while (fmod(d_lon[i], 360.0) < t_left) {
239  ++i;
240  i = i % d_lon_length;
241 
242  // If we cycle completely through all the values/indices, throw
243  if (i == lon_origin_index)
244  throw Error("geogrid: Could not find an index for the longitude value '" + double_to_string(left) + "'");
245  }
246 
247  if (fmod(d_lon[i], 360.0) == t_left)
248  longitude_index_left = i;
249  else
250  longitude_index_left = (i - 1) > 0 ? i - 1 : 0;
251 
252  DBG2(cerr << "longitude_index_left: " << longitude_index_left << endl);
253 
254  // Assume the vector is circular --> the largest value is next to the
255  // smallest.
256  int largest_lon_index = (lon_origin_index - 1 + d_lon_length) % d_lon_length;
257  i = largest_lon_index;
258  while (fmod(d_lon[i], 360.0) > t_right) {
259  // This is like modulus but for 'counting down'
260  i = (i == 0) ? d_lon_length - 1 : i - 1;
261  if (i == largest_lon_index)
262  throw Error("geogrid: Could not find an index for the longitude value '" + double_to_string(right) + "'");
263  }
264 
265  if (fmod(d_lon[i], 360.0) == t_right)
266  longitude_index_right = i;
267  else
268  longitude_index_right = (i + 1) < d_lon_length - 1 ? i + 1 : d_lon_length - 1;
269 
270  DBG2(cerr << "longitude_index_right: " << longitude_index_right << endl);
271 }
272 
285 void GeoConstraint::find_latitude_indeces(double top, double bottom,
286  LatitudeSense sense,
287  int &latitude_index_top,
288  int &latitude_index_bottom) const
289 {
290  int i, j;
291 
292  if (sense == normal) {
293  i = 0;
294  while (i < d_lat_length - 1 && top < d_lat[i])
295  ++i;
296 
297  j = d_lat_length - 1;
298  while (j > 0 && bottom > d_lat[j])
299  --j;
300 
301  if (d_lat[i] == top)
302  latitude_index_top = i;
303  else
304  latitude_index_top = (i - 1) > 0 ? i - 1 : 0;
305 
306  if (d_lat[j] == bottom)
307  latitude_index_bottom = j;
308  else
309  latitude_index_bottom =
310  (j + 1) < d_lat_length - 1 ? j + 1 : d_lat_length - 1;
311  }
312  else {
313  i = d_lat_length - 1;
314  while (i > 0 && d_lat[i] > top)
315  --i;
316 
317  j = 0;
318  while (j < d_lat_length - 1 && d_lat[j] < bottom)
319  ++j;
320 
321  if (d_lat[i] == top)
322  latitude_index_top = i;
323  else
324  latitude_index_top = (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1;
325 
326  if (d_lat[j] == bottom)
327  latitude_index_bottom = j;
328  else
329  latitude_index_bottom = (j - 1) > 0 ? j - 1 : 0;
330  }
331 }
332 
336 GeoConstraint::LatitudeSense GeoConstraint::categorize_latitude() const
337 {
338  return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted;
339 }
340 
341 // Use 'index' as the pivot point. Move the points behind index to the front of
342 // the vector and those points in front of and at index to the rear.
343 static void
344 swap_vector_ends(char *dest, char *src, int len, int index, int elem_sz)
345 {
346  memcpy(dest, src + index * elem_sz, (len - index) * elem_sz);
347 
348  memcpy(dest + (len - index) * elem_sz, src, index * elem_sz);
349 }
350 
351 template<class T>
352 static void transpose(std::vector<std::vector<T> > a,
353  std::vector<std::vector<T> > b, int width, int height)
354 {
355  for (int i = 0; i < width; i++) {
356  for (int j = 0; j < height; j++) {
357  b[j][i] = a[i][j];
358  }
359  }
360 }
361 
369 void GeoConstraint::transpose_vector(double *src, const int length)
370 {
371  double *tmp = new double[length];
372 
373  int i = 0, j = length-1;
374  while (i < length)
375  tmp[j--] = src[i++];
376 
377  memcpy(src, tmp,length * sizeof(double));
378 
379  delete[] tmp;
380 }
381 
382 static int
383 count_size_except_latitude_and_longitude(Array &a)
384 {
385  if (a.dim_end() - a.dim_begin() <= 2) // < 2 is really an error...
386  return 1;
387 
388  int size = 1;
389  for (Array::Dim_iter i = a.dim_begin(); (i + 2) != a.dim_end(); ++i)
390  size *= a.dimension_size(i, true);
391 
392  return size;
393 }
394 
395 void GeoConstraint::flip_latitude_within_array(Array &a, int lat_length,
396  int lon_length)
397 {
398  if (!d_array_data) {
399  a.read();
400  d_array_data = static_cast<char*>(a.value());
401  d_array_data_size = a.width(); // Bytes not elements
402  }
403 
404  int size = count_size_except_latitude_and_longitude(a);
405  // char *tmp_data = new char[d_array_data_size];
406  vector<char> tmp_data(d_array_data_size);
407  int array_elem_size = a.var()->width();
408  int lat_lon_size = (d_array_data_size / size);
409 
410  DBG(cerr << "lat, lon_length: " << lat_length << ", " << lon_length << endl);
411  DBG(cerr << "size: " << size << endl);
412  DBG(cerr << "d_array_data_size: " << d_array_data_size << endl);
413  DBG(cerr << "array_elem_size: " << array_elem_size<< endl);
414  DBG(cerr << "lat_lon_size: " << lat_lon_size<< endl);
415 
416  for (int i = 0; i < size; ++i) {
417  int lat = 0;
418  int s_lat = lat_length - 1;
419  // lon_length is the element size; memcpy() needs the number of bytes
420  int lon_size = array_elem_size * lon_length;
421  int offset = i * lat_lon_size;
422  while (s_lat > -1)
423  memcpy(&tmp_data[0] + offset + (lat++ * lon_size),
424  d_array_data + offset + (s_lat-- * lon_size),
425  lon_size);
426  }
427 
428  memcpy(d_array_data, &tmp_data[0], d_array_data_size);
429  // delete [] tmp_data;
430 }
431 
440 void GeoConstraint::reorder_longitude_map(int longitude_index_left)
441 {
442  double *tmp_lon = new double[d_lon_length];
443 
444  swap_vector_ends((char *) tmp_lon, (char *) d_lon, d_lon_length,
445  longitude_index_left, sizeof(double));
446 
447  memcpy(d_lon, tmp_lon, d_lon_length * sizeof(double));
448 
449  delete[]tmp_lon;
450 }
451 
452 static int
453 count_dimensions_except_longitude(Array &a)
454 {
455  int size = 1;
456  for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i)
457  size *= a.dimension_size(i, true);
458 
459  return size;
460 }
461 
479 void GeoConstraint::reorder_data_longitude_axis(Array &a, Array::Dim_iter lon_dim)
480 {
481 
482  if (!is_longitude_rightmost())
483  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.");
484 
485  DBG(cerr << "Constraint for the left half: " << get_longitude_index_left()
486  << ", " << get_lon_length() - 1 << endl);
487 
488  // Build a constraint for the left part and get those values
489  a.add_constraint(lon_dim, get_longitude_index_left(), 1,
490  get_lon_length() - 1);
491  a.set_read_p(false);
492  a.read();
493  DBG2(a.print_val(stderr));
494 
495  // Save the left-hand data to local storage
496  int left_size = a.width(); // width() == length() * element size
497  char *left_data = (char*)a.value(); // value() allocates and copies
498 
499  // Build a constraint for the 'right' part, which goes from the left edge
500  // of the array to the right index and read those data.
501  // (Don't call a.clear_constraint() since that will clear the constraint on
502  // all the dimensions while add_constraint() will replace a constraint on
503  // the given dimension).
504 
505  DBG(cerr << "Constraint for the right half: " << 0
506  << ", " << get_longitude_index_right() << endl);
507 
508  a.add_constraint(lon_dim, 0, 1, get_longitude_index_right());
509  a.set_read_p(false);
510  a.read();
511  DBG2(a.print_val(stderr));
512 
513  char *right_data = (char*)a.value();
514  int right_size = a.width();
515 
516  // Make one big lump O'data
517  d_array_data_size = left_size + right_size;
518  d_array_data = new char[d_array_data_size];
519 
520  // Assume COARDS conventions are being followed: lon varies fastest.
521  // These *_elements variables are actually elements * bytes/element since
522  // memcpy() uses bytes.
523  int elem_size = a.var()->width();
524  int left_row_size = (get_lon_length() - get_longitude_index_left()) * elem_size;
525  int right_row_size = (get_longitude_index_right() + 1) * elem_size;
526  int total_bytes_per_row = left_row_size + right_row_size;
527 
528  DBG2(cerr << "elem_size: " << elem_size << "; left & right size: "
529  << left_row_size << ", " << right_row_size << endl);
530 
531  // This will work for any number of dimension so long as longitude is the
532  // right-most array dimension.
533  int rows_to_copy = count_dimensions_except_longitude(a);
534  for (int i = 0; i < rows_to_copy; ++i) {
535  DBG(cerr << "Copying " << i << "th row" << endl);
536  DBG(cerr << "left memcpy: " << *(float *)(left_data + (left_row_size * i)) << endl);
537 
538  memcpy(d_array_data + (total_bytes_per_row * i),
539  left_data + (left_row_size * i),
540  left_row_size);
541 
542  DBG(cerr << "right memcpy: " << *(float *)(right_data + (right_row_size * i)) << endl);
543 
544  memcpy(d_array_data + (total_bytes_per_row * i) + left_row_size,
545  right_data + (right_row_size * i),
546  right_row_size);
547  }
548 
549  delete[]left_data;
550  delete[]right_data;
551 }
552 
555 GeoConstraint::GeoConstraint()
556  : d_array_data(0), d_array_data_size(0),
557  d_lat(0), d_lon(0),
558  d_bounding_box_set(false),
559  d_longitude_rightmost(false),
560  d_longitude_notation(unknown_notation),
561  d_latitude_sense(unknown_sense)
562 {
563  // Build sets of attribute values for easy searching. Maybe overkill???
564  d_coards_lat_units.insert("degrees_north");
565  d_coards_lat_units.insert("degree_north");
566  d_coards_lat_units.insert("degree_N");
567  d_coards_lat_units.insert("degrees_N");
568 
569  d_coards_lon_units.insert("degrees_east");
570  d_coards_lon_units.insert("degree_east");
571  d_coards_lon_units.insert("degrees_E");
572  d_coards_lon_units.insert("degree_E");
573 
574  d_lat_names.insert("COADSY");
575  d_lat_names.insert("lat");
576  d_lat_names.insert("Lat");
577  d_lat_names.insert("LAT");
578 
579  d_lon_names.insert("COADSX");
580  d_lon_names.insert("lon");
581  d_lon_names.insert("Lon");
582  d_lon_names.insert("LON");
583 }
584 
595 void GeoConstraint::set_bounding_box(double top, double left,
596  double bottom, double right)
597 {
598  // Ensure this method is called only once. What about pthreads?
599  // The method Array::reset_constraint() might make this so it could be
600  // called more than once. jhrg 8/30/06
601  if (d_bounding_box_set)
602  throw Error("It is not possible to register more than one geographical constraint on a variable.");
603 
604  // Record the 'sense' of the latitude for use here and later on.
605  d_latitude_sense = categorize_latitude();
606 #if 0
607  if (d_latitude_sense == inverted)
608  throw Error("geogrid() does not currently work with inverted data (data where the north pole is at a negative latitude value).");
609 #endif
610 
611  // Categorize the notation used by the bounding box (0/359 or -180/179).
612  d_longitude_notation = categorize_notation(left, right);
613 
614  // If the notation uses -180/179, transform the request to 0/359 notation.
615  if (d_longitude_notation == neg_pos)
617 
618  // If the grid uses -180/179, transform it to 0/359 as well. This will make
619  // subsequent logic easier and adds only a few extra operations, even with
620  // large maps.
621  Notation longitude_notation =
622  categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
623 
624  if (longitude_notation == neg_pos)
626 
627  if (!is_bounding_box_valid(left, top, right, bottom))
628  throw Error("The bounding box does not intersect any data within this Grid or Array. The\ngeographical extent of these data are from latitude "
629  + double_to_string(d_lat[0]) + " to "
630  + double_to_string(d_lat[d_lat_length-1])
631  + "\nand longitude " + double_to_string(d_lon[0])
632  + " to " + double_to_string(d_lon[d_lon_length-1])
633  + " while the bounding box provided was latitude "
634  + double_to_string(top) + " to "
635  + double_to_string(bottom) + "\nand longitude "
636  + double_to_string(left) + " to "
637  + double_to_string(right));
638 
639  // This is simpler than the longitude case because there's no need to
640  // test for several notations, no need to accommodate them in the return,
641  // no modulo arithmetic for the axis and no need to account for a
642  // constraint with two disconnected parts to be joined.
643  find_latitude_indeces(top, bottom, d_latitude_sense,
644  d_latitude_index_top, d_latitude_index_bottom);
645 
646 
647  // Find the longitude map indexes that correspond to the bounding box.
648  find_longitude_indeces(left, right, d_longitude_index_left,
649  d_longitude_index_right);
650 
651  DBG(cerr << "Bounding box (tlbr): " << d_latitude_index_top << ", "
652  << d_longitude_index_left << ", "
653  << d_latitude_index_bottom << ", "
654  << d_longitude_index_right << endl);
655 
656  d_bounding_box_set = true;
657 }
658 
659 } // namespace libdap