libdap++ Updated for version 3.8.2

GeoConstraint.cc

Go to the documentation of this file.
00001 
00002 // -*- mode: c++; c-basic-offset:4 -*-
00003 
00004 // This file is part of libdap, A C++ implementation of the OPeNDAP Data
00005 // Access Protocol.
00006 
00007 // Copyright (c) 2006 OPeNDAP, Inc.
00008 // Author: James Gallagher <jgallagher@opendap.org>
00009 //
00010 // This library is free software; you can redistribute it and/or
00011 // modify it under the terms of the GNU Lesser General Public
00012 // License as published by the Free Software Foundation; either
00013 // version 2.1 of the License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful,
00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00023 //
00024 // You can contact OPeNDAP, Inc. at PO Box 112, Saunderstown, RI. 02874-0112.
00025 
00026 // The Grid Selection Expression Clause class.
00027 
00028 
00029 #include "config.h"
00030 
00031 static char id[] not_used =
00032     { "$Id: GeoConstraint.cc 19914 2008-11-25 17:26:22Z jimg $"
00033     };
00034 
00035 #include <cstring>
00036 #include <cmath>
00037 #include <iostream>
00038 #include <sstream>
00039 #include <algorithm>  //  for find_if
00040 
00041 //#define DODS_DEBUG2
00042 
00043 #include "debug.h"
00044 #include "dods-datatypes.h"
00045 #include "GeoConstraint.h"
00046 #include "Float64.h"
00047 
00048 #include "Error.h"
00049 #include "InternalErr.h"
00050 #include "ce_functions.h"
00051 #include "util.h"
00052 
00053 using namespace std;
00054 
00055 namespace libdap {
00056 
00061 class is_prefix
00062 {
00063 private:
00064     string s;
00065 public:
00066     is_prefix(const string & in): s(in)
00067     {}
00068 
00069     bool operator()(const string & prefix)
00070     {
00071         return s.find(prefix) == 0;
00072     }
00073 };
00074 
00085 bool
00086 unit_or_name_match(set < string > units, set < string > names,
00087                        const string & var_units, const string & var_name)
00088 {
00089     return (units.find(var_units) != units.end()
00090             || find_if(names.begin(), names.end(),
00091                        is_prefix(var_name)) != names.end());
00092 }
00093 
00108 GeoConstraint::Notation
00109 GeoConstraint::categorize_notation(double left,
00110                                    double right) const
00111 {
00112     return (left < 0 || right < 0) ? neg_pos : pos;
00113 }
00114 
00115 /* A private method to translate the longitude constraint from -180/179
00116    notation to 0/359 notation.
00117    @param left Value-result parameter; the left side of the bounding box
00118    @parm right Value-result parameter; the right side of the bounding box */
00119 void
00120 GeoConstraint::transform_constraint_to_pos_notation(double &left,
00121         double &right) const
00122 {
00123     if (left < 0 || right < 0) {        // double check
00124         left += 180;
00125         right += 180;
00126     }
00127 }
00128 
00132 void GeoConstraint::transform_longitude_to_pos_notation()
00133 {
00134     // Assume earlier logic is correct (since the test is expensive)
00135     // for each value, add 180
00136     // Longitude could be represented using any of the numeric types
00137     for (int i = 0; i < d_lon_length; ++i)
00138         d_lon[i] += 180;
00139 }
00140 
00144 void GeoConstraint::transform_longitude_to_neg_pos_notation()
00145 {
00146     for (int i = 0; i < d_lon_length; ++i)
00147         d_lon[i] -= 180;
00148 }
00149 
00150 bool GeoConstraint::is_bounding_box_valid(double left, double top,
00151         double right, double bottom) const
00152 {
00153     if ((left < d_lon[0] && right < d_lon[0])
00154         || (left > d_lon[d_lon_length-1] && right > d_lon[d_lon_length-1]))
00155         return false;
00156 
00157     if (d_latitude_sense == normal) {
00158         // When sense is normal, the largest values are at the origin.
00159         if ((top > d_lat[0] && bottom > d_lat[0])
00160             || (top < d_lat[d_lat_length-1] && bottom < d_lat[d_lat_length-1]))
00161             return false;
00162     }
00163     else {
00164         if ((top < d_lat[0] && bottom < d_lat[0])
00165             || (top > d_lat[d_lat_length-1] && bottom > d_lat[d_lat_length-1]))
00166             return false;
00167     }
00168 
00169     return true;
00170 }
00171 
00182 void GeoConstraint::find_longitude_indeces(double left, double right,
00183         int &longitude_index_left,
00184         int &longitude_index_right) const
00185 {
00186     // Use all values 'modulo 360' to take into account the cases where the
00187     // constraint and/or the longitude vector are given using values greater
00188     // than 360 (i.e., when 380 is used to mean 20).
00189     double t_left = fmod(left, 360.0);
00190     double t_right = fmod(right, 360.0);
00191 
00192     // Find the place where 'longitude starts.' That is, what value of the
00193     // index 'i' corresponds to the smallest value of d_lon. Why we do this:
00194     // Some data sources use offset longitude axes so that the 'seam' is
00195     // shifted to a place other than the date line.
00196     int i = 0;
00197     int smallest_lon_index = 0;
00198     double smallest_lon = fmod(d_lon[smallest_lon_index], 360.0);
00199     while (i < d_lon_length - 1) {
00200         if (smallest_lon > fmod(d_lon[i], 360.0)) {
00201             smallest_lon_index = i;
00202             smallest_lon = fmod(d_lon[smallest_lon_index], 360.0);
00203         }
00204         ++i;
00205     }
00206     DBG2(cerr << "smallest_lon_index: " << smallest_lon_index << endl);
00207 
00208     // Scan from the index of the smallest value looking for the place where
00209     // the value is greater than or equal to the left most point of the bounding
00210     // box.
00211     i = smallest_lon_index;
00212     bool done = false;
00213 
00214     DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) < t_left: "
00215          << fmod(d_lon[i], 360.0) << " < " << t_left << endl);
00216 
00217     while (!done && fmod(d_lon[i], 360.0) < t_left) {
00218 
00219         DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) < t_left: "
00220              << fmod(d_lon[i], 360.0) << " < " << t_left << endl);
00221 
00222         ++i;
00223         i = i % d_lon_length;
00224         if (i == smallest_lon_index)
00225             done = true;
00226     }
00227     if (fmod(d_lon[i], 360.0) == t_left)
00228         longitude_index_left = i;
00229     else
00230         longitude_index_left = (i - 1) > 0 ? i - 1 : 0;
00231 
00232     // Assume the vector is cirular --> the largest value is next to the
00233     // smallest.
00234     int largest_lon_index =
00235         (smallest_lon_index - 1 + d_lon_length) % d_lon_length;
00236     DBG2(cerr << "largest_lon_index: " << largest_lon_index << endl);
00237     i = largest_lon_index;
00238     done = false;
00239 
00240     DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) > t_right: "
00241          << fmod(d_lon[i], 360.0) << " > " << t_right << endl);
00242 
00243     while (!done && fmod(d_lon[i], 360.0) > t_right) {
00244 
00245         DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) > t_right: "
00246              << fmod(d_lon[i], 360.0) << " > " << t_right << endl);
00247 
00248         // This is like modulus but for 'counting down'
00249         i = (i == 0) ? d_lon_length - 1 : i - 1;
00250         if (i == largest_lon_index)
00251             done = true;
00252     }
00253     if (fmod(d_lon[i], 360.0) == t_right)
00254         longitude_index_right = i;
00255     else
00256         longitude_index_right =
00257             (i + 1) < d_lon_length - 1 ? i + 1 : d_lon_length - 1;
00258 }
00259 
00272 void GeoConstraint::find_latitude_indeces(double top, double bottom,
00273         LatitudeSense sense,
00274         int &latitude_index_top,
00275         int &latitude_index_bottom) const
00276 {
00277     // Scan from the top down
00278     int i = 0;
00279     if (sense == normal) {
00280         while (i < d_lat_length - 1 && d_lat[i] > top)
00281             ++i;
00282         if (d_lat[i] == top)
00283             latitude_index_top = i;
00284         else
00285             latitude_index_top = (i - 1) > 0 ? i - 1 : 0;
00286     }
00287     else {
00288         while (i < d_lat_length - 1 && d_lat[i] < top)
00289             ++i;
00290         latitude_index_top = i;
00291     }
00292 
00293 
00294     // and from the bottom up
00295     i = d_lat_length - 1;
00296     if (sense == normal) {
00297         while (i > 0 && d_lat[i] < bottom)
00298             --i;
00299         if (d_lat[i] == bottom)
00300             latitude_index_bottom = i;
00301         else
00302             latitude_index_bottom =
00303                 (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1;
00304     }
00305     else {
00306         while (i > 0 && d_lat[i] > bottom)
00307             --i;
00308         latitude_index_bottom = i;
00309     }
00310 }
00311 
00315 GeoConstraint::LatitudeSense GeoConstraint::categorize_latitude() const
00316 {
00317     return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted;
00318 }
00319 
00320 #if 0
00321 
00328 void GeoConstraint::set_bounding_box_latitude(double top, double bottom)
00329 {
00330     // Record the 'sense' of the latitude for use here and later on.
00331     d_latitude_sense = categorize_latitude();
00332 
00333     // This is simpler than the longitude case because there's no need to test
00334     // for several notations, no need to accommodate them in the return, no
00335     // modulo arithmetic for the axis and no need to account for a constraint with
00336     // two disconnected parts to be joined.
00337     find_latitude_indeces(top, bottom, d_latitude_sense,
00338                           d_latitude_index_top, d_latitude_index_bottom);
00339 }
00340 #endif
00341 
00342 // Use 'index' as the pivot point. Move the points behind index to the front of
00343 // the vector and those points in front of and at index to the rear.
00344 static void
00345 swap_vector_ends(char *dest, char *src, int len, int index, int elem_sz)
00346 {
00347     memcpy(dest, src + index * elem_sz, (len - index) * elem_sz);
00348 
00349     memcpy(dest + (len - index) * elem_sz, src, index * elem_sz);
00350 }
00351 
00360 void GeoConstraint::reorder_longitude_map(int longitude_index_left)
00361 {
00362     double *tmp_lon = new double[d_lon_length];
00363 
00364     swap_vector_ends((char *) tmp_lon, (char *) d_lon, d_lon_length,
00365                      longitude_index_left, sizeof(double));
00366 
00367     memcpy(d_lon, tmp_lon, d_lon_length * sizeof(double));
00368 
00369     delete[]tmp_lon;
00370 }
00371 
00372 static int
00373 count_dimensions_except_longitude(Array &a)
00374 {
00375     int size = 1;
00376     for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i)
00377         size *= a.dimension_size(i, true);
00378 
00379     return size;
00380 }
00381 
00400 void GeoConstraint::reorder_data_longitude_axis(Array &a)
00401 {
00402 
00403     if (!get_longitude_rightmost())
00404         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.");
00405 
00406     DBG(cerr << "Constraint for the left half: " << get_longitude_index_left()
00407         << ", " << get_lon_length() - 1 << endl);
00408 
00409     a.add_constraint(d_lon_dim, get_longitude_index_left(), 1,
00410                      get_lon_length() - 1);
00411     a.set_read_p(false);
00412     a.read();
00413     DBG2(a.print_val(stderr));
00414 #if 0
00415     char *left_data = 0;
00416     int left_size = a.buf2val((void **) & left_data);
00417 #endif
00418     char *left_data = (char*)a.value();
00419     int left_size = a.length();
00420 
00421     // Build a constraint for the 'right' part, which
00422     // goes from the left edge of the array to the right index and read those
00423     // data.
00424     //a.clear_constraint();
00425 
00426     DBG(cerr << "Constraint for the right half: " << 0
00427         << ", " << get_longitude_index_right() << endl);
00428 
00429     a.add_constraint(d_lon_dim, 0, 1, get_longitude_index_right());
00430     a.set_read_p(false);
00431     a.read();
00432     DBG2(a.print_val(stderr));
00433 #if 0
00434     char *right_data = 0;
00435     int right_size = a.buf2val((void **) & right_data);
00436 #endif
00437     char *right_data = (char*)a.value();
00438     int right_size = a.length();
00439 
00440     // Make one big lump O'data
00441     d_array_data_size = left_size + right_size;
00442     d_array_data = new char[d_array_data_size];
00443 
00444     // Assume COARDS conventions are being followed: lon varies fastest.
00445     // These *_elements variables are actually elements * bytes/element since
00446     // memcpy() uses bytes.
00447     int elem_width = a.var()->width();
00448     int left_elements = (get_lon_length() - get_longitude_index_left()) * elem_width;
00449     int right_elements = (get_longitude_index_right() + 1) * elem_width;
00450     int total_elements_per_row = left_elements + right_elements;
00451 
00452     // Interleave the left and right_data vectors. jhrg 8/31/06
00453     int rows_to_copy = count_dimensions_except_longitude(a);
00454     for (int i = 0; i < rows_to_copy; ++i) {
00455         memcpy(d_array_data + (total_elements_per_row * i),
00456                left_data + (left_elements * i),
00457                left_elements);
00458         memcpy(d_array_data + left_elements + (total_elements_per_row * i),
00459                right_data + (right_elements * i),
00460                right_elements);
00461     }
00462 
00463     delete[]left_data;
00464     delete[]right_data;
00465 }
00466 
00467 #if 0
00468 
00481 void GeoConstraint::set_bounding_box_longitude(double left, double right)
00482 {
00483     // Categorize the notation used by the bounding box (0/359 or -180/179).
00484     d_longitude_notation = categorize_notation(left, right);
00485 
00486     // If the notation uses -180/179, transform the request to 0/359 notation.
00487     if (d_longitude_notation == neg_pos)
00488         transform_constraint_to_pos_notation(left, right);
00489 
00490     // If the grid uses -180/179, transform it to 0/359 as well. This will make
00491     // subsequent logic easier and adds only a few extra operations, even with
00492     // large maps.
00493     Notation longitude_notation =
00494         categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
00495 
00496     if (longitude_notation == neg_pos)
00497         transform_longitude_to_pos_notation();
00498 
00499     // Find the longitude map indeces that correspond to the bounding box.
00500     find_longitude_indeces(left, right, d_longitude_index_left,
00501                            d_longitude_index_right);
00502 }
00503 #endif
00504 
00510 GeoConstraint::GeoConstraint()
00511         : d_array_data(0), d_array_data_size(0),
00512         d_lat(0), d_lon(0),
00513         d_bounding_box_set(false),
00514         d_longitude_notation(unknown_notation),
00515         d_latitude_sense(unknown_sense)
00516 {
00517     // Build sets of attribute values for easy searching. Maybe overkill???
00518     d_coards_lat_units.insert("degrees_north");
00519     d_coards_lat_units.insert("degree_north");
00520     d_coards_lat_units.insert("degree_N");
00521     d_coards_lat_units.insert("degrees_N");
00522 
00523     d_coards_lon_units.insert("degrees_east");
00524     d_coards_lon_units.insert("degree_east");
00525     d_coards_lon_units.insert("degrees_E");
00526     d_coards_lon_units.insert("degree_E");
00527 
00528     d_lat_names.insert("COADSY");
00529     d_lat_names.insert("lat");
00530     d_lat_names.insert("Lat");
00531     d_lat_names.insert("LAT");
00532 
00533     d_lon_names.insert("COADSX");
00534     d_lon_names.insert("lon");
00535     d_lon_names.insert("Lon");
00536     d_lon_names.insert("LON");
00537 }
00538 
00549 void GeoConstraint::set_bounding_box(double left, double top,
00550                                      double right, double bottom)
00551 {
00552     // Ensure this method is called only once. What about pthreads?
00553     // The method Array::reset_constraint() might make this so it could be
00554     // called more than once. jhrg 8/30/06
00555     if (d_bounding_box_set)
00556         throw
00557         InternalErr
00558         ("It is not possible to register more than one geographical constraint on a variable.");
00559 
00560     // Record the 'sense' of the latitude for use here and later on.
00561     d_latitude_sense = categorize_latitude();
00562 
00563     // Categorize the notation used by the bounding box (0/359 or -180/179).
00564     d_longitude_notation = categorize_notation(left, right);
00565 
00566     // If the notation uses -180/179, transform the request to 0/359 notation.
00567     if (d_longitude_notation == neg_pos)
00568         transform_constraint_to_pos_notation(left, right);
00569 
00570     // If the grid uses -180/179, transform it to 0/359 as well. This will make
00571     // subsequent logic easier and adds only a few extra operations, even with
00572     // large maps.
00573     Notation longitude_notation =
00574         categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
00575 
00576     if (longitude_notation == neg_pos)
00577         transform_longitude_to_pos_notation();
00578 
00579     if (!is_bounding_box_valid(left, top, right, bottom))
00580         throw Error("The bounding box does not intersect any data within this Grid or Array. The\ngeographical extent of these data are from latitude "
00581                     + double_to_string(d_lat[0]) + " to "
00582                     + double_to_string(d_lat[d_lat_length-1])
00583                     + "\nand longitude " + double_to_string(d_lon[0])
00584                     + " to " + double_to_string(d_lon[d_lon_length-1])
00585                     + " while the bounding box provided was latitude "
00586                     + double_to_string(top) + " to "
00587                     + double_to_string(bottom) + "\nand longitude "
00588                     + double_to_string(left) + " to "
00589                     + double_to_string(right));
00590 
00591     // This is simpler than the longitude case because there's no need to
00592     // test for several notations, no need to accommodate them in the return,
00593     // no modulo arithmetic for the axis and no need to account for a
00594     // constraint with two disconnected parts to be joined.
00595     find_latitude_indeces(top, bottom, d_latitude_sense,
00596                           d_latitude_index_top, d_latitude_index_bottom);
00597 
00598 
00599     // Find the longitude map indeces that correspond to the bounding box.
00600     find_longitude_indeces(left, right, d_longitude_index_left,
00601                            d_longitude_index_right);
00602 
00603     DBG(cerr << "Bounding box (tlbr): " << d_latitude_index_top << ", "
00604         << d_longitude_index_left << ", "
00605         << d_latitude_index_bottom << ", "
00606         << d_longitude_index_right << endl);
00607 
00608     d_bounding_box_set = true;
00609 }
00610 
00611 } // namespace libdap