GridGeoConstraint.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) 2002,2003 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: GridGeoConstraint.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 
00041 //#define DODS_DEBUG
00042 
00043 #include "debug.h"
00044 #include "dods-datatypes.h"
00045 #include "GridGeoConstraint.h"
00046 #include "Float64.h"
00047 
00048 #include "Error.h"
00049 #include "InternalErr.h"
00050 #include "ce_functions.h"
00051 
00052 using namespace std;
00053 using namespace libdap;
00054 
00055 
00063 GridGeoConstraint::GridGeoConstraint(Grid *grid, const string &ds_name)
00064         : GeoConstraint(ds_name), d_grid(grid), d_latitude(0), d_longitude(0)
00065 {
00066     if (d_grid->get_array()->dimensions() < 2
00067         || d_grid->get_array()->dimensions() > 3)
00068         throw Error("The geogrid() function works only with Grids of two or three dimensions.");
00069 
00070     // Is this Grid a geo-referenced grid? Throw Error if not.
00071     if (!build_lat_lon_maps())
00072         throw Error(string("The grid '") + d_grid->name()
00073                     +
00074                     "' does not have identifiable latitude/longitude map vectors.");
00075 
00076     if (!lat_lon_dimensions_ok())
00077         throw Error("The geogrid() function will only work when the Grid's Longitude and Latitude\nmaps are the rightmost dimensions.");
00078 }
00079 
00095 bool GridGeoConstraint::build_lat_lon_maps()
00096 {
00097     Grid::Map_iter m = d_grid->map_begin();
00098     // Assume that a Grid is correct and thus has exactly as many maps at its
00099     // array part has dimensions. Thus don't bother to test the Grid's array
00100     // dimension iterator for '!= dim_end()'.
00101     Array::Dim_iter d = d_grid->get_array()->dim_begin();
00102     // The fields d_latitude and d_longitude are initialized to null
00103     while (m != d_grid->map_end() && (!d_latitude || !d_longitude)) {
00104         string units_value = (*m)->get_attr_table().get_attr("units");
00105         remove_quotes(units_value);
00106         string map_name = (*m)->name();
00107 
00108         // The 'units' attribute must match exactly; the name only needs to
00109         // match a prefix.
00110         if (!d_latitude
00111             && unit_or_name_match(get_coards_lat_units(), get_lat_names(),
00112                                   units_value, map_name)) {
00113 
00114             // Set both d_latitude (a pointer to the real map vector) and
00115             // d_lat, a vector of the values represented as doubles. It's easier
00116             // to work with d_lat, but it's d_latitude that needs to be set
00117             // when constraining the grid. Also, record the grid variable's
00118             // dimension iterator so that it's easier to set the Grid's Array
00119             // (which also has to be constrained).
00120             d_latitude = dynamic_cast < Array * >(*m);
00121             if (!d_latitude->read_p())
00122                 d_latitude->read(get_dataset());
00123 
00124             set_lat(extract_double_array(d_latitude));   // throws Error
00125             set_lat_length(d_latitude->length());
00126 
00127             set_lat_dim(d);
00128         }
00129 
00130         if (!d_longitude        // && !units_value.empty()
00131             && unit_or_name_match(get_coards_lon_units(), get_lon_names(),
00132                                   units_value, map_name)) {
00133 
00134             d_longitude = dynamic_cast < Array * >(*m);
00135             if (!d_longitude->read_p())
00136                 d_longitude->read(get_dataset());
00137 
00138             set_lon(extract_double_array(d_longitude));
00139             set_lon_length(d_longitude->length());
00140 
00141             set_lon_dim(d);
00142         }
00143 
00144         ++m;
00145         ++d;
00146     }
00147 
00148     return get_lat() && get_lon();
00149 }
00150 
00161 bool
00162 GridGeoConstraint::lat_lon_dimensions_ok()
00163 {
00164     // get the last two map iterators
00165     Grid::Map_riter rightmost = d_grid->map_rbegin();
00166     Grid::Map_riter next_rightmost = rightmost + 1;
00167 
00168     if (*rightmost == d_longitude && *next_rightmost == d_latitude)
00169         set_longitude_rightmost(true);
00170     else if (*rightmost == d_latitude && *next_rightmost == d_longitude)
00171         set_longitude_rightmost(false);
00172     else
00173         return false;
00174 
00175     return true;
00176 }
00177 
00199 void GridGeoConstraint::apply_constraint_to_data()
00200 {
00201     if (!get_bounding_box_set())
00202         throw
00203         InternalErr
00204         ("The Latitude and Longitude constraints must be set before calling apply_constraint_to_data().");
00205 
00206     Array::Dim_iter fd = d_latitude->dim_begin();
00207     if (get_latitude_sense() == inverted) {
00208         int tmp = get_latitude_index_top();
00209         set_latitude_index_top(get_latitude_index_bottom());
00210         set_latitude_index_bottom(tmp);
00211     }
00212 
00213     // It's esy to flip the Latitude values; if the bottom index value
00214     // is before/above the top index, return an error explaining that.
00215     if (get_latitude_index_top() > get_latitude_index_bottom())
00216         throw Error("The upper and lower latitude indices appear to be reversed. Please provide the latitude bounding box numbers giving the northern-most latitude first.");
00217 
00218     // Constrain the lat vector and lat dim of the array
00219     d_latitude->add_constraint(fd, get_latitude_index_top(), 1,
00220                                get_latitude_index_bottom());
00221     d_grid->get_array()->add_constraint(get_lat_dim(),
00222                                         get_latitude_index_top(), 1,
00223                                         get_latitude_index_bottom());
00224 
00225     // Does the longitude constraint cross the edge of the longitude vector?
00226     // If so, reorder the grid's data (array), longitude map vector and the
00227     // local vector of longitude data used for computation.
00228     if (get_longitude_index_left() > get_longitude_index_right()) {
00229         reorder_longitude_map(get_longitude_index_left());
00230 
00231         // If the longitude constraint is 'split', join the two parts, reload
00232         // the data into the Grid's Array and make sure the Array is marked as
00233         // already read. This should be true for the whole Grid, but if some
00234         // future modification changes that, the array will be covered here.
00235         // Note that the following method only reads the data out and stores
00236         // it in this object after joining the two parts. The method
00237         // apply_constraint_to_data() transfers the data back from the this
00238         // objet to the DAP Grid variable.
00239         reorder_data_longitude_axis(*d_grid->get_array());
00240 
00241         // Now the data are all in local storage
00242 
00243         // alter the indices; the left index has now been moved to 0, and the right
00244         // index is now at lon_vector_length-left+right.
00245         set_longitude_index_right(get_lon_length() - get_longitude_index_left()
00246                                   + get_longitude_index_right());
00247         set_longitude_index_left(0);
00248     }
00249     // If the constraint used the -180/179 (neg_pos) notation, transform
00250     // the longitude map s it uses the -180/179 notation. Note that at this
00251     // point, d_longitude always uses the pos notation because of the earlier
00252     // conditional transforamtion.
00253 
00254     // Do this _before_ applying the constraint since set_array_using_double()
00255     // tests the array length using Vector::length() and that method returns
00256     // the length _as constrained_. We want to move all of the longitude
00257     // values from d_lon back into the map, not just the number that will be
00258     // sent (although an optimization might do this, it's hard to imagine
00259     // it would gain much).
00260     if (get_longitude_notation() == neg_pos) {
00261         transform_longitude_to_neg_pos_notation();
00262     }
00263     // Apply constraint; stride is always one and maps only have one dimension
00264     fd = d_longitude->dim_begin();
00265     d_longitude->add_constraint(fd, get_longitude_index_left(), 1,
00266                                 get_longitude_index_right());
00267 
00268     d_grid->get_array()->add_constraint(get_lon_dim(),
00269                                         get_longitude_index_left(),
00270                                         1, get_longitude_index_right());
00271 
00272     // Transfer values from the local lat vector to the Grid's
00273     set_array_using_double(d_latitude, get_lat() + get_latitude_index_top(),
00274                            get_latitude_index_bottom() - get_latitude_index_top() + 1);
00275 
00276     set_array_using_double(d_longitude, get_lon() + get_longitude_index_left(),
00277                            get_longitude_index_right() - get_longitude_index_left() + 1);
00278 
00279     // ... and then the Grid's array if it has been read.
00280     if (get_array_data()) {
00281         int size = d_grid->get_array()->val2buf(get_array_data());
00282         if (size != get_array_data_size())
00283             throw
00284             InternalErr
00285             ("Expected data size not copied to the Grid's buffer.");
00286         d_grid->set_read_p(true);
00287     }
00288     else {
00289         d_grid->get_array()->read(get_dataset());
00290     }
00291 }

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