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 16967 2007-08-20 23:39:01Z 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 
00291 void GeoConstraint::find_latitude_indeces(double top, double bottom,
00292         LatitudeSense sense,
00293         int &latitude_index_top,
00294         int &latitude_index_bottom) const
00295 {
00296     // Scan from the top down
00297     int i = 0;
00298     if (sense == normal) {
00299         while (i < d_lat_length - 1 && d_lat[i] > top)
00300             ++i;
00301         if (d_lat[i] == top)
00302             latitude_index_top = i;
00303         else
00304             latitude_index_top = (i - 1) > 0 ? i - 1 : 0;
00305     }
00306     else {
00307         while (i < d_lat_length - 1 && d_lat[i] < top)
00308             ++i;
00309         latitude_index_top = i;
00310     }
00311 
00312 
00313     // and from the bottom up
00314     i = d_lat_length - 1;
00315     if (sense == normal) {
00316         while (i > 0 && d_lat[i] < bottom)
00317             --i;
00318         if (d_lat[i] == bottom)
00319             latitude_index_bottom = i;
00320         else
00321             latitude_index_bottom =
00322                 (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1;
00323     }
00324     else {
00325         while (i > 0 && d_lat[i] > bottom)
00326             --i;
00327         latitude_index_bottom = i;
00328     }
00329 }
00330 
00334 GeoConstraint::LatitudeSense GeoConstraint::categorize_latitude() const
00335 {
00336     return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted;
00337 }
00338 
00339 #if 0
00340 
00347 void GeoConstraint::set_bounding_box_latitude(double top, double bottom)
00348 {
00349     // Record the 'sense' of the latitude for use here and later on.
00350     d_latitude_sense = categorize_latitude();
00351 
00352     // This is simpler than the longitude case because there's no need to test
00353     // for several notations, no need to accommodate them in the return, no
00354     // modulo arithmetic for the axis and no need to account for a constraint with
00355     // two disconnected parts to be joined.
00356     find_latitude_indeces(top, bottom, d_latitude_sense,
00357                           d_latitude_index_top, d_latitude_index_bottom);
00358 }
00359 #endif
00360 
00361 // Use 'index' as the pivot point. Move the points behind index to the front of
00362 // the vector and those points in front of and at index to the rear.
00363 static void
00364 swap_vector_ends(char *dest, char *src, int len, int index, int elem_sz)
00365 {
00366     memcpy(dest, src + index * elem_sz, (len - index) * elem_sz);
00367 
00368     memcpy(dest + (len - index) * elem_sz, src, index * elem_sz);
00369 }
00370 
00379 void GeoConstraint::reorder_longitude_map(int longitude_index_left)
00380 {
00381     double *tmp_lon = new double[d_lon_length];
00382 
00383     swap_vector_ends((char *) tmp_lon, (char *) d_lon, d_lon_length,
00384                      longitude_index_left, sizeof(double));
00385 
00386     memcpy(d_lon, tmp_lon, d_lon_length * sizeof(double));
00387 
00388     delete[]tmp_lon;
00389 }
00390 
00391 static int
00392 count_dimensions_except_longitude(Array &a)
00393 {
00394     int size = 1;
00395     for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i)
00396         size *= a.dimension_size(i, true);
00397 
00398     return size;
00399 }
00400 
00419 void GeoConstraint::reorder_data_longitude_axis(Array &a)
00420 {
00421 
00422     if (!get_longitude_rightmost())
00423         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.");
00424 
00425     DBG(cerr << "Constraint for the left half: " << get_longitude_index_left()
00426         << ", " << get_lon_length() - 1 << endl);
00427 
00428     a.add_constraint(d_lon_dim, get_longitude_index_left(), 1,
00429                      get_lon_length() - 1);
00430     a.set_read_p(false);
00431     a.read(get_dataset());
00432     DBG2(a.print_val(stderr));
00433 #if 0
00434     char *left_data = 0;
00435     int left_size = a.buf2val((void **) & left_data);
00436 #endif
00437     char *left_data = (char*)a.value();
00438     int left_size = a.length();
00439     
00440     // Build a constraint for the 'right' part, which
00441     // goes from the left edge of the array to the right index and read those
00442     // data.
00443     //a.clear_constraint();
00444 
00445     DBG(cerr << "Constraint for the right half: " << 0
00446         << ", " << get_longitude_index_right() << endl);
00447 
00448     a.add_constraint(d_lon_dim, 0, 1, get_longitude_index_right());
00449     a.set_read_p(false);
00450     a.read(get_dataset());
00451     DBG2(a.print_val(stderr));
00452 #if 0
00453     char *right_data = 0;
00454     int right_size = a.buf2val((void **) & right_data);
00455 #endif
00456     char *right_data = (char*)a.value();
00457     int right_size = a.length();
00458 
00459     // Make one big lump O'data
00460     d_array_data_size = left_size + right_size;
00461     d_array_data = new char[d_array_data_size];
00462 
00463     // Assume COARDS convensions are being followed: lon varies fastest.
00464     // These *_elements variables are actually elements * bytes/element since
00465     // memcpy() uses bytes.
00466     int elem_width = a.var()->width();
00467     int left_elements = (get_lon_length() - get_longitude_index_left()) * elem_width;
00468     int right_elements = (get_longitude_index_right() + 1) * elem_width;
00469     int total_elements_per_row = left_elements + right_elements;
00470 
00471     // Interleve the left and right_data vectors. jhrg 8/31/06
00472     int rows_to_copy = count_dimensions_except_longitude(a);
00473     for (int i = 0; i < rows_to_copy; ++i) {
00474         memcpy(d_array_data + (total_elements_per_row * i),
00475                left_data + (left_elements * i),
00476                left_elements);
00477         memcpy(d_array_data + left_elements + (total_elements_per_row * i),
00478                right_data + (right_elements * i),
00479                right_elements);
00480     }
00481 
00482     delete[]left_data;
00483     delete[]right_data;
00484 }
00485 
00486 #if 0
00487 
00500 void GeoConstraint::set_bounding_box_longitude(double left, double right)
00501 {
00502     // Categorize the notation used by the bounding box (0/359 or -180/179).
00503     d_longitude_notation = categorize_notation(left, right);
00504 
00505     // If the notation uses -180/179, transform the request to 0/359 notation.
00506     if (d_longitude_notation == neg_pos)
00507         transform_constraint_to_pos_notation(left, right);
00508 
00509     // If the grid uses -180/179, transform it to 0/359 as well. This will make
00510     // subsequent logic easier and adds only a few extra operations, even with
00511     // large maps.
00512     Notation longitude_notation =
00513         categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
00514 
00515     if (longitude_notation == neg_pos)
00516         transform_longitude_to_pos_notation();
00517 
00518     // Find the longitude map indeces that correspond to the bounding box.
00519     find_longitude_indeces(left, right, d_longitude_index_left,
00520                            d_longitude_index_right);
00521 }
00522 #endif
00523 
00528 GeoConstraint::GeoConstraint(const string & ds_name)
00529         : d_dataset(ds_name), d_array_data(0), d_array_data_size(0),
00530         d_lat(0), d_lon(0),
00531         d_bounding_box_set(false),
00532         d_longitude_notation(unknown_notation),
00533         d_latitude_sense(unknown_sense)
00534 {
00535     // Build sets of attribute values for easy searching. Maybe overkill???
00536     d_coards_lat_units.insert("degrees_north");
00537     d_coards_lat_units.insert("degree_north");
00538     d_coards_lat_units.insert("degree_N");
00539     d_coards_lat_units.insert("degrees_N");
00540 
00541     d_coards_lon_units.insert("degrees_east");
00542     d_coards_lon_units.insert("degree_east");
00543     d_coards_lon_units.insert("degrees_E");
00544     d_coards_lon_units.insert("degree_E");
00545 
00546     d_lat_names.insert("COADSY");
00547     d_lat_names.insert("lat");
00548     d_lat_names.insert("Lat");
00549     d_lat_names.insert("LAT");
00550 
00551     d_lon_names.insert("COADSX");
00552     d_lon_names.insert("lon");
00553     d_lon_names.insert("Lon");
00554     d_lon_names.insert("LON");
00555 }
00556 
00567 void GeoConstraint::set_bounding_box(double left, double top,
00568                                      double right, double bottom)
00569 {
00570     // Ensure this method is called only once. What about pthreads?
00571     // The method Array::reset_constraint() might make this so it could be
00572     // called more than once. jhrg 8/30/06
00573     if (d_bounding_box_set)
00574         throw
00575         InternalErr
00576         ("It is not possible to register more than one geographical constraint on a variable.");
00577 
00578     // Record the 'sense' of the latitude for use here and later on.
00579     d_latitude_sense = categorize_latitude();
00580 
00581     // Categorize the notation used by the bounding box (0/359 or -180/179).
00582     d_longitude_notation = categorize_notation(left, right);
00583 
00584     // If the notation uses -180/179, transform the request to 0/359 notation.
00585     if (d_longitude_notation == neg_pos)
00586         transform_constraint_to_pos_notation(left, right);
00587 
00588     // If the grid uses -180/179, transform it to 0/359 as well. This will make
00589     // subsequent logic easier and adds only a few extra operations, even with
00590     // large maps.
00591     Notation longitude_notation =
00592         categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
00593 
00594     if (longitude_notation == neg_pos)
00595         transform_longitude_to_pos_notation();
00596 
00597     if (!is_bounding_box_valid(left, top, right, bottom))
00598         throw Error("The bounding box does not intersect any data within this Grid or Array. The\ngeographical extent of these data are from latitude "
00599                     + double_to_string(d_lat[0]) + " to " 
00600                     + double_to_string(d_lat[d_lat_length-1]) 
00601                     + "\nand longitude " + double_to_string(d_lon[0]) 
00602                     + " to " + double_to_string(d_lon[d_lon_length-1])
00603                     + " while the bounding box provided was latitude "
00604                     + double_to_string(top) + " to " 
00605                     + double_to_string(bottom) + "\nand longitude " 
00606                     + double_to_string(left) + " to " 
00607                     + double_to_string(right));
00608 
00609     // This is simpler than the longitude case because there's no need to
00610     // test for several notations, no need to accommodate them in the return,
00611     // no modulo arithmetic for the axis and no need to account for a
00612     // constraint with two disconnected parts to be joined.
00613     find_latitude_indeces(top, bottom, d_latitude_sense,
00614                           d_latitude_index_top, d_latitude_index_bottom);
00615 
00616 
00617     // Find the longitude map indeces that correspond to the bounding box.
00618     find_longitude_indeces(left, right, d_longitude_index_left,
00619                            d_longitude_index_right);
00620 
00621     DBG(cerr << "Bounding box (tlbr): " << d_latitude_index_top << ", "
00622         << d_longitude_index_left << ", "
00623         << d_latitude_index_bottom << ", "
00624         << d_longitude_index_right << endl);
00625 
00626     d_bounding_box_set = true;
00627 }

Generated on Wed Nov 14 03:15:44 2007 for libdap++ by  doxygen 1.5.1