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 16088 2007-03-28 21:42:19Z jimg $"
00033     };
00034 
00035 #include <math.h>
00036 #include <string.h>
00037 
00038 #include <iostream>
00039 #include <sstream>
00040 #ifdef WIN32
00041 #include <algorithm>  //  for find_if
00042 #endif
00043 
00044 //#define DODS_DEBUG2
00045 
00046 #include "debug.h"
00047 #include "dods-datatypes.h"
00048 #include "GeoConstraint.h"
00049 #include "Float64.h"
00050 
00051 #include "Error.h"
00052 #include "InternalErr.h"
00053 #include "ce_functions.h"
00054 #include "util.h"
00055 
00056 using namespace std;
00057 using namespace libdap;
00058 
00065 void remove_quotes(string & value)
00066 {
00067     if (!value.empty() && value[0] == '"'
00068         && value[value.length() - 1] == '"') {
00069         value.erase(0, 1);
00070         value.erase(value.length() - 1);
00071     }
00072 
00073     return;
00074 }
00075 
00080 class is_prefix
00081 {
00082 private:
00083     string s;
00084 public:
00085     is_prefix(const string & in): s(in)
00086     {}
00087 
00088     bool operator()(const string & prefix)
00089     {
00090         return s.find(prefix) == 0;
00091     }
00092 };
00093 
00104 bool
00105 unit_or_name_match(set < string > units, set < string > names,
00106                        const string & var_units, const string & var_name)
00107 {
00108     return (units.find(var_units) != units.end()
00109             || find_if(names.begin(), names.end(),
00110                        is_prefix(var_name)) != names.end());
00111 }
00112 
00127 GeoConstraint::Notation
00128 GeoConstraint::categorize_notation(double left,
00129                                    double right) const
00130 {
00131     return (left < 0 || right < 0) ? neg_pos : pos;
00132 }
00133 
00134 /* A private method to translate the longitude constraint from -180/179
00135    notation to 0/359 notation.
00136    @param left Value-result parameter; the left side of the bounding box
00137    @parm right Value-result parameter; the right side of the bounding box */
00138 void
00139 GeoConstraint::transform_constraint_to_pos_notation(double &left,
00140         double &right) const
00141 {
00142     if (left < 0 || right < 0) {        // double check
00143         left += 180;
00144         right += 180;
00145     }
00146 }
00147 
00151 void GeoConstraint::transform_longitude_to_pos_notation()
00152 {
00153     // Assume earlier logic is correct (since the test is expensive)
00154     // for each value, add 180
00155     // Longitude could be represented using any of the numeric types
00156     for (int i = 0; i < d_lon_length; ++i)
00157         d_lon[i] += 180;
00158 }
00159 
00163 void GeoConstraint::transform_longitude_to_neg_pos_notation()
00164 {
00165     for (int i = 0; i < d_lon_length; ++i)
00166         d_lon[i] -= 180;
00167 }
00168 
00169 bool GeoConstraint::is_bounding_box_valid(double left, double top,
00170         double right, double bottom) const
00171 {
00172     if ((left < d_lon[0] && right < d_lon[0])
00173         || (left > d_lon[d_lon_length-1] && right > d_lon[d_lon_length-1]))
00174         return false;
00175 
00176     if (d_latitude_sense == normal) {
00177         // When sense is normal, the largest values are at the origin.
00178         if ((top > d_lat[0] && bottom > d_lat[0])
00179             || (top < d_lat[d_lat_length-1] && bottom < d_lat[d_lat_length-1]))
00180             return false;
00181     }
00182     else {
00183         if ((top < d_lat[0] && bottom < d_lat[0])
00184             || (top > d_lat[d_lat_length-1] && bottom > d_lat[d_lat_length-1]))
00185             return false;
00186     }
00187 
00188     return true;
00189 }
00190 
00201 void GeoConstraint::find_longitude_indeces(double left, double right,
00202         int &longitude_index_left,
00203         int &longitude_index_right) const
00204 {
00205     // Use all values 'modulo 360' to take into account the cases where the
00206     // constraint and/or the longitude vector are given using values greater
00207     // than 360 (i.e., when 380 is used to mean 20).
00208     double t_left = fmod(left, 360.0);
00209     double t_right = fmod(right, 360.0);
00210 
00211     // Find the place where 'longitude starts.' That is, what value of the
00212     // index 'i' corresponds to the smallest value of d_lon. Why we do this:
00213     // Some data sources use offset longitude axes so that the 'seam' is
00214     // shifted to a place other than the date line.
00215     int i = 0;
00216     int smallest_lon_index = 0;
00217     double smallest_lon = fmod(d_lon[smallest_lon_index], 360.0);
00218     while (i < d_lon_length - 1) {
00219         if (smallest_lon > fmod(d_lon[i], 360.0)) {
00220             smallest_lon_index = i;
00221             smallest_lon = fmod(d_lon[smallest_lon_index], 360.0);
00222         }
00223         ++i;
00224     }
00225     DBG2(cerr << "smallest_lon_index: " << smallest_lon_index << endl);
00226 
00227     // Scan from the index of the smallest value looking for the place where
00228     // the value is greater than or equal to the left most point of the bounding
00229     // box.
00230     i = smallest_lon_index;
00231     bool done = false;
00232 
00233     DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) < t_left: "
00234          << fmod(d_lon[i], 360.0) << " < " << t_left << endl);
00235 
00236     while (!done && fmod(d_lon[i], 360.0) < t_left) {
00237 
00238         DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) < t_left: "
00239              << fmod(d_lon[i], 360.0) << " < " << t_left << endl);
00240 
00241         ++i;
00242         i = i % d_lon_length;
00243         if (i == smallest_lon_index)
00244             done = true;
00245     }
00246     if (fmod(d_lon[i], 360.0) == t_left)
00247         longitude_index_left = i;
00248     else
00249         longitude_index_left = (i - 1) > 0 ? i - 1 : 0;
00250 
00251     // Assume the vector is cirular --> the largest value is next to the
00252     // smallest.
00253     int largest_lon_index =
00254         (smallest_lon_index - 1 + d_lon_length) % d_lon_length;
00255     DBG2(cerr << "largest_lon_index: " << largest_lon_index << endl);
00256     i = largest_lon_index;
00257     done = false;
00258 
00259     DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) > t_right: "
00260          << fmod(d_lon[i], 360.0) << " > " << t_right << endl);
00261 
00262     while (!done && fmod(d_lon[i], 360.0) > t_right) {
00263 
00264         DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) > t_right: "
00265              << fmod(d_lon[i], 360.0) << " > " << t_right << endl);
00266 
00267         // This is like modulus but for 'counting down'
00268         i = (i == 0) ? d_lon_length - 1 : i - 1;
00269         if (i == largest_lon_index)
00270             done = true;
00271     }
00272     if (fmod(d_lon[i], 360.0) == t_right)
00273         longitude_index_right = i;
00274     else
00275         longitude_index_right =
00276             (i + 1) < d_lon_length - 1 ? i + 1 : d_lon_length - 1;
00277 }
00278 
00289 void GeoConstraint::find_latitude_indeces(double top, double bottom,
00290         LatitudeSense sense,
00291         int &latitude_index_top,
00292         int &latitude_index_bottom) const
00293 {
00294     // Scan from the top down
00295     int i = 0;
00296     if (sense == normal) {
00297         while (i < d_lat_length - 1 && d_lat[i] > top)
00298             ++i;
00299         if (d_lat[i] == top)
00300             latitude_index_top = i;
00301         else
00302             latitude_index_top = (i - 1) > 0 ? i - 1 : 0;
00303     }
00304     else {
00305         while (i < d_lat_length - 1 && d_lat[i] < top)
00306             ++i;
00307         latitude_index_top = i;
00308     }
00309 
00310 
00311     // and from the bottom up
00312     i = d_lat_length - 1;
00313     if (sense == normal) {
00314         while (i > 0 && d_lat[i] < bottom)
00315             --i;
00316         if (d_lat[i] == bottom)
00317             latitude_index_bottom = i;
00318         else
00319             latitude_index_bottom =
00320                 (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1;
00321     }
00322     else {
00323         while (i > 0 && d_lat[i] > bottom)
00324             --i;
00325         latitude_index_bottom = i;
00326     }
00327 }
00328 
00332 GeoConstraint::LatitudeSense GeoConstraint::categorize_latitude() const
00333 {
00334     return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted;
00335 }
00336 
00337 #if 0
00338 
00345 void GeoConstraint::set_bounding_box_latitude(double top, double bottom)
00346 {
00347     // Record the 'sense' of the latitude for use here and later on.
00348     d_latitude_sense = categorize_latitude();
00349 
00350     // This is simpler than the longitude case because there's no need to test
00351     // for several notations, no need to accommodate them in the return, no
00352     // modulo arithmetic for the axis and no need to account for a constraint with
00353     // two disconnected parts to be joined.
00354     find_latitude_indeces(top, bottom, d_latitude_sense,
00355                           d_latitude_index_top, d_latitude_index_bottom);
00356 }
00357 #endif
00358 
00359 // Use 'index' as the pivot point. Move the behind index to the front of
00360 // the vector and those points in front of and at index to the rear.
00361 static void
00362 swap_vector_ends(char *dest, char *src, int len, int index, int elem_sz)
00363 {
00364     memcpy(dest, src + index * elem_sz, (len - index) * elem_sz);
00365 
00366     memcpy(dest + (len - index) * elem_sz, src, index * elem_sz);
00367 }
00368 
00377 void GeoConstraint::reorder_longitude_map(int longitude_index_left)
00378 {
00379     double *tmp_lon = new double[d_lon_length];
00380 
00381     swap_vector_ends((char *) tmp_lon, (char *) d_lon, d_lon_length,
00382                      longitude_index_left, sizeof(double));
00383 
00384     memcpy(d_lon, tmp_lon, d_lon_length * sizeof(double));
00385 
00386     delete[]tmp_lon;
00387 }
00388 
00389 static int
00390 count_dimensions_except_longitude(Array &a)
00391 {
00392     int size = 1;
00393     for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i)
00394         size *= a.dimension_size(i, true);
00395 
00396     return size;
00397 }
00398 
00417 void GeoConstraint::reorder_data_longitude_axis(Array &a)
00418 {
00419 
00420     if (!get_longitude_rightmost())
00421         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.");
00422 
00423     DBG(cerr << "Constraint for the left half: " << get_longitude_index_left()
00424         << ", " << get_lon_length() - 1 << endl);
00425 
00426     a.add_constraint(d_lon_dim, get_longitude_index_left(), 1,
00427                      get_lon_length() - 1);
00428     a.set_read_p(false);
00429     a.read(get_dataset());
00430     DBG2(a.print_val(stderr));
00431     char *left_data = 0;
00432     int left_size = a.buf2val((void **) & left_data);
00433 
00434     // Build a constraint for the 'right' part, which
00435     // goes from the left edge of the array to the right index and read those
00436     // data.
00437     //a.clear_constraint();
00438 
00439     DBG(cerr << "Constraint for the right half: " << 0
00440         << ", " << get_longitude_index_right() << endl);
00441 
00442     a.add_constraint(d_lon_dim, 0, 1, get_longitude_index_right());
00443     a.set_read_p(false);
00444     a.read(get_dataset());
00445     DBG2(a.print_val(stderr));
00446     char *right_data = 0;
00447     int right_size = a.buf2val((void **) & right_data);
00448 
00449     // Make one big lump O'data
00450     d_array_data_size = left_size + right_size;
00451     d_array_data = new char[d_array_data_size];
00452 
00453     // Assume COARDS convensions are being followed: lon varies fastest.
00454     // These *_elements variables are actually elements * bytes/element since
00455     // memcpy() uses bytes.
00456     int elem_width = a.var()->width();
00457     int left_elements = (get_lon_length() - get_longitude_index_left()) * elem_width;
00458     int right_elements = (get_longitude_index_right() + 1) * elem_width;
00459     int total_elements_per_row = left_elements + right_elements;
00460 
00461     // Interleve the left and right_data vectors. jhrg 8/31/06
00462     int rows_to_copy = count_dimensions_except_longitude(a);
00463     for (int i = 0; i < rows_to_copy; ++i) {
00464         memcpy(d_array_data + (total_elements_per_row * i),
00465                left_data + (left_elements * i),
00466                left_elements);
00467         memcpy(d_array_data + left_elements + (total_elements_per_row * i),
00468                right_data + (right_elements * i),
00469                right_elements);
00470     }
00471 
00472     delete[]left_data;
00473     delete[]right_data;
00474 }
00475 
00476 #if 0
00477 
00490 void GeoConstraint::set_bounding_box_longitude(double left, double right)
00491 {
00492     // Categorize the notation used by the bounding box (0/359 or -180/179).
00493     d_longitude_notation = categorize_notation(left, right);
00494 
00495     // If the notation uses -180/179, transform the request to 0/359 notation.
00496     if (d_longitude_notation == neg_pos)
00497         transform_constraint_to_pos_notation(left, right);
00498 
00499     // If the grid uses -180/179, transform it to 0/359 as well. This will make
00500     // subsequent logic easier and adds only a few extra operations, even with
00501     // large maps.
00502     Notation longitude_notation =
00503         categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
00504 
00505     if (longitude_notation == neg_pos)
00506         transform_longitude_to_pos_notation();
00507 
00508     // Find the longitude map indeces that correspond to the bounding box.
00509     find_longitude_indeces(left, right, d_longitude_index_left,
00510                            d_longitude_index_right);
00511 }
00512 #endif
00513 
00521 GeoConstraint::GeoConstraint(const string & ds_name)
00522         : d_dataset(ds_name), d_array_data(0), d_array_data_size(0),
00523         d_lat(0), d_lon(0),
00524         d_bounding_box_set(false),
00525         d_longitude_notation(unknown_notation),
00526         d_latitude_sense(unknown_sense)
00527 {
00528     // Build sets of attribute values for easy searching. Maybe overkill???
00529     d_coards_lat_units.insert("degrees_north");
00530     d_coards_lat_units.insert("degree_north");
00531     d_coards_lat_units.insert("degree_N");
00532     d_coards_lat_units.insert("degrees_N");
00533 
00534     d_coards_lon_units.insert("degrees_east");
00535     d_coards_lon_units.insert("degree_east");
00536     d_coards_lon_units.insert("degrees_E");
00537     d_coards_lon_units.insert("degree_E");
00538 
00539     d_lat_names.insert("COADSY");
00540     d_lat_names.insert("lat");
00541     d_lat_names.insert("Lat");
00542     d_lat_names.insert("LAT");
00543 
00544     d_lon_names.insert("COADSX");
00545     d_lon_names.insert("lon");
00546     d_lon_names.insert("Lon");
00547     d_lon_names.insert("LON");
00548 }
00549 
00560 void GeoConstraint::set_bounding_box(double left, double top,
00561                                      double right, double bottom)
00562 {
00563     // Ensure this method is called only once. What about pthreads?
00564     // The method Array::reset_constraint() might make this so it could be
00565     // called more than once. jhrg 8/30/06
00566     if (d_bounding_box_set)
00567         throw
00568         InternalErr
00569         ("It is not possible to register more than one geographical constraint on a variable.");
00570 
00571     // Record the 'sense' of the latitude for use here and later on.
00572     d_latitude_sense = categorize_latitude();
00573 
00574     // Categorize the notation used by the bounding box (0/359 or -180/179).
00575     d_longitude_notation = categorize_notation(left, right);
00576 
00577     // If the notation uses -180/179, transform the request to 0/359 notation.
00578     if (d_longitude_notation == neg_pos)
00579         transform_constraint_to_pos_notation(left, right);
00580 
00581     // If the grid uses -180/179, transform it to 0/359 as well. This will make
00582     // subsequent logic easier and adds only a few extra operations, even with
00583     // large maps.
00584     Notation longitude_notation =
00585         categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
00586 
00587     if (longitude_notation == neg_pos)
00588         transform_longitude_to_pos_notation();
00589 
00590     if (!is_bounding_box_valid(left, top, right, bottom))
00591         throw Error("The bounding box does not intersect any data within this Grid or Array. The\ngeographical extent of these data are from latitude "
00592                     + double_to_string(d_lat[0]) + " to " 
00593                     + double_to_string(d_lat[d_lat_length-1]) 
00594                     + "\nand longitude " + double_to_string(d_lon[0]) 
00595                     + " to " + double_to_string(d_lon[d_lon_length-1])
00596                     + " while the bounding box provided was latitude "
00597                     + double_to_string(top) + " to " 
00598                     + double_to_string(bottom) + "\nand longitude " 
00599                     + double_to_string(left) + " to " 
00600                     + double_to_string(right));
00601 
00602     // This is simpler than the longitude case because there's no need to
00603     // test for several notations, no need to accommodate them in the return,
00604     // no modulo arithmetic for the axis and no need to account for a
00605     // constraint with two disconnected parts to be joined.
00606     find_latitude_indeces(top, bottom, d_latitude_sense,
00607                           d_latitude_index_top, d_latitude_index_bottom);
00608 
00609 
00610     // Find the longitude map indeces that correspond to the bounding box.
00611     find_longitude_indeces(left, right, d_longitude_index_left,
00612                            d_longitude_index_right);
00613 
00614     DBG(cerr << "Bounding box (tlbr): " << d_latitude_index_top << ", "
00615         << d_longitude_index_left << ", "
00616         << d_latitude_index_bottom << ", "
00617         << d_longitude_index_right << endl);
00618 
00619     d_bounding_box_set = true;
00620 }

Generated on Wed Jun 27 12:56:39 2007 for libdap++ by  doxygen 1.4.7