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 25112 2011-12-29 21:44:54Z 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 }
430 
439 void GeoConstraint::reorder_longitude_map(int longitude_index_left)
440 {
441  double *tmp_lon = new double[d_lon_length];
442 
443  swap_vector_ends((char *) tmp_lon, (char *) d_lon, d_lon_length,
444  longitude_index_left, sizeof(double));
445 
446  memcpy(d_lon, tmp_lon, d_lon_length * sizeof(double));
447 
448  delete[]tmp_lon;
449 }
450 
451 static int
452 count_dimensions_except_longitude(Array &a)
453 {
454  int size = 1;
455  for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i)
456  size *= a.dimension_size(i, true);
457 
458  return size;
459 }
460 
478 void GeoConstraint::reorder_data_longitude_axis(Array &a, Array::Dim_iter lon_dim)
479 {
480 
481  if (!is_longitude_rightmost())
482  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.");
483 
484  DBG(cerr << "Constraint for the left half: " << get_longitude_index_left()
485  << ", " << get_lon_length() - 1 << endl);
486 
487  // Build a constraint for the left part and get those values
488  a.add_constraint(lon_dim, get_longitude_index_left(), 1,
489  get_lon_length() - 1);
490  a.set_read_p(false);
491  a.read();
492  DBG2(a.print_val(stderr));
493 
494  // Save the left-hand data to local storage
495  int left_size = a.width(); // width() == length() * element size
496  char *left_data = (char*)a.value(); // value() allocates and copies
497 
498  // Build a constraint for the 'right' part, which goes from the left edge
499  // of the array to the right index and read those data.
500  // (Don't call a.clear_constraint() since that will clear the constraint on
501  // all the dimensions while add_constraint() will replace a constraint on
502  // the given dimension).
503 
504  DBG(cerr << "Constraint for the right half: " << 0
505  << ", " << get_longitude_index_right() << endl);
506 
507  a.add_constraint(lon_dim, 0, 1, get_longitude_index_right());
508  a.set_read_p(false);
509  a.read();
510  DBG2(a.print_val(stderr));
511 
512  char *right_data = (char*)a.value();
513  int right_size = a.width();
514 
515  // Make one big lump O'data
516  d_array_data_size = left_size + right_size;
517  d_array_data = new char[d_array_data_size];
518 
519  // Assume COARDS conventions are being followed: lon varies fastest.
520  // These *_elements variables are actually elements * bytes/element since
521  // memcpy() uses bytes.
522  int elem_size = a.var()->width();
523  int left_row_size = (get_lon_length() - get_longitude_index_left()) * elem_size;
524  int right_row_size = (get_longitude_index_right() + 1) * elem_size;
525  int total_bytes_per_row = left_row_size + right_row_size;
526 
527  DBG2(cerr << "elem_size: " << elem_size << "; left & right size: "
528  << left_row_size << ", " << right_row_size << endl);
529 
530  // This will work for any number of dimension so long as longitude is the
531  // right-most array dimension.
532  int rows_to_copy = count_dimensions_except_longitude(a);
533  for (int i = 0; i < rows_to_copy; ++i) {
534  DBG(cerr << "Copying " << i << "th row" << endl);
535  DBG(cerr << "left memcpy: " << *(float *)(left_data + (left_row_size * i)) << endl);
536 
537  memcpy(d_array_data + (total_bytes_per_row * i),
538  left_data + (left_row_size * i),
539  left_row_size);
540 
541  DBG(cerr << "right memcpy: " << *(float *)(right_data + (right_row_size * i)) << endl);
542 
543  memcpy(d_array_data + (total_bytes_per_row * i) + left_row_size,
544  right_data + (right_row_size * i),
545  right_row_size);
546  }
547 
548  delete[]left_data;
549  delete[]right_data;
550 }
551 
554 GeoConstraint::GeoConstraint()
555  : d_array_data(0), d_array_data_size(0),
556  d_lat(0), d_lon(0),
557  d_bounding_box_set(false),
558  d_longitude_rightmost(false),
559  d_longitude_notation(unknown_notation),
560  d_latitude_sense(unknown_sense)
561 {
562  // Build sets of attribute values for easy searching. Maybe overkill???
563  d_coards_lat_units.insert("degrees_north");
564  d_coards_lat_units.insert("degree_north");
565  d_coards_lat_units.insert("degree_N");
566  d_coards_lat_units.insert("degrees_N");
567 
568  d_coards_lon_units.insert("degrees_east");
569  d_coards_lon_units.insert("degree_east");
570  d_coards_lon_units.insert("degrees_E");
571  d_coards_lon_units.insert("degree_E");
572 
573  d_lat_names.insert("COADSY");
574  d_lat_names.insert("lat");
575  d_lat_names.insert("Lat");
576  d_lat_names.insert("LAT");
577 
578  d_lon_names.insert("COADSX");
579  d_lon_names.insert("lon");
580  d_lon_names.insert("Lon");
581  d_lon_names.insert("LON");
582 }
583 
594 void GeoConstraint::set_bounding_box(double top, double left,
595  double bottom, double right)
596 {
597  // Ensure this method is called only once. What about pthreads?
598  // The method Array::reset_constraint() might make this so it could be
599  // called more than once. jhrg 8/30/06
600  if (d_bounding_box_set)
601  throw Error("It is not possible to register more than one geographical constraint on a variable.");
602 
603  // Record the 'sense' of the latitude for use here and later on.
604  d_latitude_sense = categorize_latitude();
605 #if 0
606  if (d_latitude_sense == inverted)
607  throw Error("geogrid() does not currently work with inverted data (data where the north pole is at a negative latitude value).");
608 #endif
609 
610  // Categorize the notation used by the bounding box (0/359 or -180/179).
611  d_longitude_notation = categorize_notation(left, right);
612 
613  // If the notation uses -180/179, transform the request to 0/359 notation.
614  if (d_longitude_notation == neg_pos)
616 
617  // If the grid uses -180/179, transform it to 0/359 as well. This will make
618  // subsequent logic easier and adds only a few extra operations, even with
619  // large maps.
620  Notation longitude_notation =
621  categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
622 
623  if (longitude_notation == neg_pos)
625 
626  if (!is_bounding_box_valid(left, top, right, bottom))
627  throw Error("The bounding box does not intersect any data within this Grid or Array. The\ngeographical extent of these data are from latitude "
628  + double_to_string(d_lat[0]) + " to "
629  + double_to_string(d_lat[d_lat_length-1])
630  + "\nand longitude " + double_to_string(d_lon[0])
631  + " to " + double_to_string(d_lon[d_lon_length-1])
632  + " while the bounding box provided was latitude "
633  + double_to_string(top) + " to "
634  + double_to_string(bottom) + "\nand longitude "
635  + double_to_string(left) + " to "
636  + double_to_string(right));
637 
638  // This is simpler than the longitude case because there's no need to
639  // test for several notations, no need to accommodate them in the return,
640  // no modulo arithmetic for the axis and no need to account for a
641  // constraint with two disconnected parts to be joined.
642  find_latitude_indeces(top, bottom, d_latitude_sense,
643  d_latitude_index_top, d_latitude_index_bottom);
644 
645 
646  // Find the longitude map indexes that correspond to the bounding box.
647  find_longitude_indeces(left, right, d_longitude_index_left,
648  d_longitude_index_right);
649 
650  DBG(cerr << "Bounding box (tlbr): " << d_latitude_index_top << ", "
651  << d_longitude_index_left << ", "
652  << d_latitude_index_bottom << ", "
653  << d_longitude_index_right << endl);
654 
655  d_bounding_box_set = true;
656 }
657 
658 } // namespace libdap