bes Updated for version 3.20.13
GeoConstraint.cc
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#include "config.h"
29
30#include <cstring>
31#include <cmath>
32#include <iostream>
33#include <sstream>
34#include <algorithm> // for find_if
35
36//#define DODS_DEBUG
37//#define DODS_DEBUG2
38
39#include <libdap/Float64.h>
40#include <libdap/Array.h>
41#include <libdap/Error.h>
42#include <libdap/InternalErr.h>
43#include <libdap/dods-datatypes.h>
44#include <libdap/util.h>
45#include <libdap/debug.h>
46
47#include "GeoConstraint.h"
48
49using namespace std;
50using namespace libdap;
51
52namespace functions {
53
58class is_prefix
59{
60private:
61 string s;
62public:
63 is_prefix(const string & in): s(in)
64 {}
65
66 bool operator()(const string & prefix)
67 {
68 return s.find(prefix) == 0;
69 }
70};
71
82bool
83unit_or_name_match(set < string > units, set < string > names,
84 const string & var_units, const string & var_name)
85{
86 return (units.find(var_units) != units.end()
87 || find_if(names.begin(), names.end(),
88 is_prefix(var_name)) != names.end());
89}
90
107 const double right) const
108{
109 return (left < 0 || right < 0) ? neg_pos : pos;
110}
111
112/* A private method to translate the longitude constraint from -180/179
113 notation to 0/359 notation.
114
115 About the two notations:
116 0 180 360
117 |---------------------------|-------------------------|
118 0 180,-180 0
119
120 so in the neg-pos notation (using the name I give it in this class) both
121 180 and -180 are the same longitude. And in the pos notation 0 and 360 are
122 the same.
123
124 @param left Value-result parameter; the left side of the bounding box
125 @parm right Value-result parameter; the right side of the bounding box */
126void
127GeoConstraint::transform_constraint_to_pos_notation(double &left,
128 double &right) const
129{
130 if (left < 0)
131 left += 360 ;
132
133 if (right < 0)
134 right += 360;
135}
136
146{
147 // Assume earlier logic is correct (since the test is expensive)
148 // for each value, add 180
149 // Longitude could be represented using any of the numeric types
150 for (int i = 0; i < d_lon_length; ++i)
151 if (d_lon[i] < 0)
152 d_lon[i] += 360;
153}
154
165{
166 for (int i = 0; i < d_lon_length; ++i)
167 if (d_lon[i] > 180)
168 d_lon[i] -= 360;
169}
170
171bool GeoConstraint::is_bounding_box_valid(const double left, const double top,
172 const double right, const double bottom) const
173{
174 if ((left < d_lon[0] && right < d_lon[0])
175 || (left > d_lon[d_lon_length-1] && right > d_lon[d_lon_length-1]))
176 return false;
177
178 if (d_latitude_sense == normal) {
179 // When sense is normal, the largest values are at the origin.
180 if ((top > d_lat[0] && bottom > d_lat[0])
181 || (top < d_lat[d_lat_length-1] && bottom < d_lat[d_lat_length-1]))
182 return false;
183 }
184 else {
185 if ((top < d_lat[0] && bottom < d_lat[0])
186 || (top > d_lat[d_lat_length-1] && bottom > d_lat[d_lat_length-1]))
187 return false;
188 }
189
190 return true;
191}
192
203void GeoConstraint::find_longitude_indeces(double left, double right,
204 int &longitude_index_left, int &longitude_index_right) const
205{
206 // Use all values 'modulo 360' to take into account the cases where the
207 // constraint and/or the longitude vector are given using values greater
208 // than 360 (i.e., when 380 is used to mean 20).
209 double t_left = fmod(left, 360.0);
210 double t_right = fmod(right, 360.0);
211
212 // Find the place where 'longitude starts.' That is, what value of the
213 // index 'i' corresponds to the smallest value of d_lon. Why we do this:
214 // Some data sources use offset longitude axes so that the 'seam' is
215 // shifted to a place other than the date line.
216 int i = 0;
217 int lon_origin_index = 0;
218 double smallest_lon = fmod(d_lon[0], 360.0);
219 while (i < d_lon_length) {
220 double curent_lon_value = fmod(d_lon[i], 360.0);
221 if (smallest_lon > curent_lon_value) {
222 smallest_lon = curent_lon_value;
223 lon_origin_index = i;
224 }
225 ++i;
226 }
227
228 DBG2(cerr << "lon_origin_index: " << lon_origin_index << endl);
229
230 // Scan from the index of the smallest value looking for the place where
231 // the value is greater than or equal to the left most point of the bounding
232 // box.
233 i = lon_origin_index;
234 while (fmod(d_lon[i], 360.0) < t_left) {
235 ++i;
236 i = i % d_lon_length;
237
238 // If we cycle completely through all the values/indices, throw
239 if (i == lon_origin_index)
240 throw Error("geogrid: Could not find an index for the longitude value '" + double_to_string(left) + "'");
241 }
242
243 if (fmod(d_lon[i], 360.0) == t_left)
244 longitude_index_left = i;
245 else
246 longitude_index_left = (i - 1) > 0 ? i - 1 : 0;
247
248 DBG2(cerr << "longitude_index_left: " << longitude_index_left << endl);
249
250 // Assume the vector is circular --> the largest value is next to the
251 // smallest.
252 int largest_lon_index = (lon_origin_index - 1 + d_lon_length) % d_lon_length;
253 i = largest_lon_index;
254 while (fmod(d_lon[i], 360.0) > t_right) {
255 // This is like modulus but for 'counting down'
256 i = (i == 0) ? d_lon_length - 1 : i - 1;
257 if (i == largest_lon_index)
258 throw Error("geogrid: Could not find an index for the longitude value '" + double_to_string(right) + "'");
259 }
260
261 if (fmod(d_lon[i], 360.0) == t_right)
262 longitude_index_right = i;
263 else
264 longitude_index_right = (i + 1) < d_lon_length - 1 ? i + 1 : d_lon_length - 1;
265
266 DBG2(cerr << "longitude_index_right: " << longitude_index_right << endl);
267}
268
281void GeoConstraint::find_latitude_indeces(double top, double bottom,
282 LatitudeSense sense,
283 int &latitude_index_top,
284 int &latitude_index_bottom) const
285{
286 int i, j;
287
288 if (sense == normal) {
289 i = 0;
290 while (i < d_lat_length - 1 && top < d_lat[i])
291 ++i;
292
293 j = d_lat_length - 1;
294 while (j > 0 && bottom > d_lat[j])
295 --j;
296
297 if (d_lat[i] == top)
298 latitude_index_top = i;
299 else
300 latitude_index_top = (i - 1) > 0 ? i - 1 : 0;
301
302 if (d_lat[j] == bottom)
303 latitude_index_bottom = j;
304 else
305 latitude_index_bottom =
306 (j + 1) < d_lat_length - 1 ? j + 1 : d_lat_length - 1;
307 }
308 else {
309 i = d_lat_length - 1;
310 while (i > 0 && d_lat[i] > top)
311 --i;
312
313 j = 0;
314 while (j < d_lat_length - 1 && d_lat[j] < bottom)
315 ++j;
316
317 if (d_lat[i] == top)
318 latitude_index_top = i;
319 else
320 latitude_index_top = (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1;
321
322 if (d_lat[j] == bottom)
323 latitude_index_bottom = j;
324 else
325 latitude_index_bottom = (j - 1) > 0 ? j - 1 : 0;
326 }
327}
328
333{
334 return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted;
335}
336
337// Use 'index' as the pivot point. Move the points behind index to the front of
338// the vector and those points in front of and at index to the rear.
339static void
340swap_vector_ends(char *dest, char *src, int len, int index, int elem_sz)
341{
342 memcpy(dest, src + index * elem_sz, (len - index) * elem_sz);
343
344 memcpy(dest + (len - index) * elem_sz, src, index * elem_sz);
345}
346
347template<class T>
348static void transpose(std::vector<std::vector<T> > a,
349 std::vector<std::vector<T> > b, int width, int height)
350{
351 for (int i = 0; i < width; i++) {
352 for (int j = 0; j < height; j++) {
353 b[j][i] = a[i][j];
354 }
355 }
356}
357
365void GeoConstraint::transpose_vector(double *src, const int length)
366{
367 double *tmp = new double[length];
368
369 int i = 0, j = length-1;
370 while (i < length)
371 tmp[j--] = src[i++];
372
373 memcpy(src, tmp,length * sizeof(double));
374
375 delete[] tmp;
376}
377
378static int
379count_size_except_latitude_and_longitude(Array &a)
380{
381 if (a.dim_end() - a.dim_begin() <= 2) // < 2 is really an error...
382 return 1;
383
384 int size = 1;
385 for (Array::Dim_iter i = a.dim_begin(); (i + 2) != a.dim_end(); ++i)
386 size *= a.dimension_size(i, true);
387
388 return size;
389}
390
391void GeoConstraint::flip_latitude_within_array(Array &a, int lat_length,
392 int lon_length)
393{
394 if (!d_array_data) {
395 a.read();
396 d_array_data = static_cast<char*>(a.value());
397 d_array_data_size = a.width(true); // Bytes not elements
398 }
399
400 int size = count_size_except_latitude_and_longitude(a);
401 // char *tmp_data = new char[d_array_data_size];
402 vector<char> tmp_data(d_array_data_size);
403 int array_elem_size = a.var()->width(true);
404 int lat_lon_size = (d_array_data_size / size);
405
406 DBG(cerr << "lat, lon_length: " << lat_length << ", " << lon_length << endl);
407 DBG(cerr << "size: " << size << endl);
408 DBG(cerr << "d_array_data_size: " << d_array_data_size << endl);
409 DBG(cerr << "array_elem_size: " << array_elem_size<< endl);
410 DBG(cerr << "lat_lon_size: " << lat_lon_size<< endl);
411
412 for (int i = 0; i < size; ++i) {
413 int lat = 0;
414 int s_lat = lat_length - 1;
415 // lon_length is the element size; memcpy() needs the number of bytes
416 int lon_size = array_elem_size * lon_length;
417 int offset = i * lat_lon_size;
418 while (s_lat > -1)
419 memcpy(tmp_data.data() + offset + (lat++ * lon_size),
420 d_array_data + offset + (s_lat-- * lon_size),
421 lon_size);
422 }
423
424 memcpy(d_array_data, tmp_data.data(), d_array_data_size);
425}
426
435void GeoConstraint::reorder_longitude_map(int longitude_index_left)
436{
437 double *tmp_lon = new double[d_lon_length];
438
439 swap_vector_ends((char *) tmp_lon, (char *) d_lon, d_lon_length,
440 longitude_index_left, sizeof(double));
441
442 memcpy(d_lon, tmp_lon, d_lon_length * sizeof(double));
443
444 delete[]tmp_lon;
445}
446
447static int
448count_dimensions_except_longitude(Array &a)
449{
450 int size = 1;
451 for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i)
452 size *= a.dimension_size(i, true);
453
454 return size;
455}
456
474void GeoConstraint::reorder_data_longitude_axis(Array &a, Array::Dim_iter lon_dim)
475{
476
477 if (!is_longitude_rightmost())
478 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
480 DBG(cerr << "Constraint for the left half: " << get_longitude_index_left()
481 << ", " << get_lon_length() - 1 << endl);
482
483 // Build a constraint for the left part and get those values
484 a.add_constraint(lon_dim, get_longitude_index_left(), 1,
485 get_lon_length() - 1);
486 a.set_read_p(false);
487 a.read();
488 DBG2(a.print_val(stderr));
489
490 // Save the left-hand data to local storage
491 int left_size = a.width(true); // width() == length() * element size
492 char *left_data = (char*)a.value(); // value() allocates and copies
493
494 // Build a constraint for the 'right' part, which goes from the left edge
495 // of the array to the right index and read those data.
496 // (Don't call a.clear_constraint() since that will clear the constraint on
497 // all the dimensions while add_constraint() will replace a constraint on
498 // the given dimension).
499
500 DBG(cerr << "Constraint for the right half: " << 0
501 << ", " << get_longitude_index_right() << endl);
502
503 a.add_constraint(lon_dim, 0, 1, get_longitude_index_right());
504 a.set_read_p(false);
505 a.read();
506 DBG2(a.print_val(stderr));
507
508 char *right_data = (char*)a.value();
509 int right_size = a.width(true);
510
511 // Make one big lump O'data
512 d_array_data_size = left_size + right_size;
513 d_array_data = new char[d_array_data_size];
514
515 // Assume COARDS conventions are being followed: lon varies fastest.
516 // These *_elements variables are actually elements * bytes/element since
517 // memcpy() uses bytes.
518 int elem_size = a.var()->width(true);
519 int left_row_size = (get_lon_length() - get_longitude_index_left()) * elem_size;
520 int right_row_size = (get_longitude_index_right() + 1) * elem_size;
521 int total_bytes_per_row = left_row_size + right_row_size;
522
523 DBG2(cerr << "elem_size: " << elem_size << "; left & right size: "
524 << left_row_size << ", " << right_row_size << endl);
525
526 // This will work for any number of dimension so long as longitude is the
527 // right-most array dimension.
528 int rows_to_copy = count_dimensions_except_longitude(a);
529 for (int i = 0; i < rows_to_copy; ++i) {
530 DBG(cerr << "Copying " << i << "th row" << endl);
531 DBG(cerr << "left memcpy: " << *(float *)(left_data + (left_row_size * i)) << endl);
532
533 memcpy(d_array_data + (total_bytes_per_row * i),
534 left_data + (left_row_size * i),
535 left_row_size);
536
537 DBG(cerr << "right memcpy: " << *(float *)(right_data + (right_row_size * i)) << endl);
538
539 memcpy(d_array_data + (total_bytes_per_row * i) + left_row_size,
540 right_data + (right_row_size * i),
541 right_row_size);
542 }
543
544 delete[]left_data;
545 delete[]right_data;
546}
547
551 : d_array_data(0), d_array_data_size(0),
552 d_lat(0), d_lon(0), d_lat_length(0), d_lon_length(0),
553 d_latitude_index_top(0),
554 d_latitude_index_bottom(0),
555 d_longitude_index_left(0),
556 d_longitude_index_right(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
594void 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)
615 transform_constraint_to_pos_notation(left, right);
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
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)