libdap++ Updated for version 3.8.2

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 23447 2010-08-27 19:38:04Z jimg $"
00033     };
00034 
00035 #include <cmath>
00036 
00037 #include <iostream>
00038 #include <sstream>
00039 
00040 //#define DODS_DEBUG
00041 
00042 #include "debug.h"
00043 #include "dods-datatypes.h"
00044 #include "GridGeoConstraint.h"
00045 #include "Float64.h"
00046 
00047 #include "Error.h"
00048 #include "InternalErr.h"
00049 #include "ce_functions.h"
00050 #include "util.h"
00051 
00052 using namespace std;
00053 
00054 namespace libdap {
00055 
00064 GridGeoConstraint::GridGeoConstraint(Grid *grid)
00065         : GeoConstraint(), d_grid(grid), d_latitude(0), d_longitude(0)
00066 {
00067     if (d_grid->get_array()->dimensions() < 2
00068         || d_grid->get_array()->dimensions() > 3)
00069         throw Error("The geogrid() function works only with Grids of two or three dimensions.");
00070 
00071     // Is this Grid a geo-referenced grid? Throw Error if not.
00072     if (!build_lat_lon_maps())
00073         throw Error(string("The grid '") + d_grid->name()
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 
00080 GridGeoConstraint::GridGeoConstraint(Grid *grid, Array *lat, Array *lon)
00081         : GeoConstraint(), d_grid(grid), d_latitude(0), d_longitude(0)
00082 {
00083     if (d_grid->get_array()->dimensions() < 2
00084         || d_grid->get_array()->dimensions() > 3)
00085         throw Error("The geogrid() function works only with Grids of two or three dimensions.");
00086 
00087     // Is this Grid a geo-referenced grid? Throw Error if not.
00088     if (!build_lat_lon_maps(lat, lon))
00089         throw Error(string("The grid '") + d_grid->name()
00090                     + "' does not have valid latitude/longitude map vectors.");
00091 
00092 
00093     if (!lat_lon_dimensions_ok())
00094         throw Error("The geogrid() function will only work when the Grid's Longitude and Latitude\nmaps are the rightmost dimensions.");
00095 }
00096 
00112 bool GridGeoConstraint::build_lat_lon_maps()
00113 {
00114     Grid::Map_iter m = d_grid->map_begin();
00115     // Assume that a Grid is correct and thus has exactly as many maps as its
00116     // array part has dimensions. Thus don't bother to test the Grid's array
00117     // dimension iterator for '!= dim_end()'.
00118     Array::Dim_iter d = d_grid->get_array()->dim_begin();
00119     // The fields d_latitude and d_longitude may be initialized to null or they
00120     // may already contain pointers to the maps to use. In the latter case,
00121     // skip the heuristics used in this code. However, given that all this
00122     // method does is find the lat/lon maps, if they are given in the ctor,
00123     // This method will likely not be called at all.
00124     while (m != d_grid->map_end() && (!d_latitude || !d_longitude)) {
00125         string units_value = (*m)->get_attr_table().get_attr("units");
00126         units_value = remove_quotes(units_value);
00127         string map_name = (*m)->name();
00128 
00129         // The 'units' attribute must match exactly; the name only needs to
00130         // match a prefix.
00131         if (!d_latitude
00132             && unit_or_name_match(get_coards_lat_units(), get_lat_names(),
00133                                   units_value, map_name)) {
00134 
00135             // Set both d_latitude (a pointer to the real map vector) and
00136             // d_lat, a vector of the values represented as doubles. It's easier
00137             // to work with d_lat, but it's d_latitude that needs to be set
00138             // when constraining the grid. Also, record the grid variable's
00139             // dimension iterator so that it's easier to set the Grid's Array
00140             // (which also has to be constrained).
00141             d_latitude = dynamic_cast < Array * >(*m);
00142             if (!d_latitude)
00143                 throw InternalErr(__FILE__, __LINE__, "Expected an array.");
00144             if (!d_latitude->read_p())
00145                 d_latitude->read();
00146 
00147             set_lat(extract_double_array(d_latitude));   // throws Error
00148             set_lat_length(d_latitude->length());
00149 
00150             set_lat_dim(d);
00151         }
00152 
00153         if (!d_longitude        // && !units_value.empty()
00154             && unit_or_name_match(get_coards_lon_units(), get_lon_names(),
00155                                   units_value, map_name)) {
00156 
00157             d_longitude = dynamic_cast < Array * >(*m);
00158             if (!d_longitude)
00159                 throw InternalErr(__FILE__, __LINE__, "Expected an array.");
00160             if (!d_longitude->read_p())
00161                 d_longitude->read();
00162 
00163             set_lon(extract_double_array(d_longitude));
00164             set_lon_length(d_longitude->length());
00165 
00166             set_lon_dim(d);
00167 
00168             if (m + 1 == d_grid->map_end())
00169                 set_longitude_rightmost(true);
00170         }
00171 
00172         ++m;
00173         ++d;
00174     }
00175 
00176     return get_lat() && get_lon();
00177 }
00178 
00186 bool GridGeoConstraint::build_lat_lon_maps(Array *lat, Array *lon)
00187 {
00188     Grid::Map_iter m = d_grid->map_begin();
00189 
00190     Array::Dim_iter d = d_grid->get_array()->dim_begin();
00191 
00192     while (m != d_grid->map_end() && (!d_latitude || !d_longitude)) {
00193         // Look for the Grid map that matches the variable passed as 'lat'
00194         if (!d_latitude && *m == lat) {
00195 
00196             d_latitude = lat;
00197 
00198             if (!d_latitude->read_p())
00199                 d_latitude->read();
00200 
00201             set_lat(extract_double_array(d_latitude));   // throws Error
00202             set_lat_length(d_latitude->length());
00203 
00204             set_lat_dim(d);
00205         }
00206 
00207         if (!d_longitude && *m == lon) {
00208 
00209             d_longitude = lon;
00210 
00211             if (!d_longitude->read_p())
00212                 d_longitude->read();
00213 
00214             set_lon(extract_double_array(d_longitude));
00215             set_lon_length(d_longitude->length());
00216 
00217             set_lon_dim(d);
00218 
00219             if (m + 1 == d_grid->map_end())
00220                 set_longitude_rightmost(true);
00221         }
00222 
00223         ++m;
00224         ++d;
00225     }
00226 
00227     return get_lat() && get_lon();
00228 }
00229 
00240 bool
00241 GridGeoConstraint::lat_lon_dimensions_ok()
00242 {
00243     // get the last two map iterators
00244     Grid::Map_riter rightmost = d_grid->map_rbegin();
00245     Grid::Map_riter next_rightmost = rightmost + 1;
00246 
00247     if (*rightmost == d_longitude && *next_rightmost == d_latitude)
00248         set_longitude_rightmost(true);
00249     else if (*rightmost == d_latitude && *next_rightmost == d_longitude)
00250         set_longitude_rightmost(false);
00251     else
00252         return false;
00253 
00254     return true;
00255 }
00256 
00278 void GridGeoConstraint::apply_constraint_to_data()
00279 {
00280     if (!is_bounding_box_set())
00281         throw InternalErr("The Latitude and Longitude constraints must be set before calling apply_constraint_to_data().");
00282 
00283     Array::Dim_iter fd = d_latitude->dim_begin();
00284 
00285     if (get_latitude_sense() == inverted) {
00286         int tmp = get_latitude_index_top();
00287         set_latitude_index_top(get_latitude_index_bottom());
00288         set_latitude_index_bottom(tmp);
00289     }
00290 
00291     // It's easy to flip the Latitude values; if the bottom index value
00292     // is before/above the top index, return an error explaining that.
00293     if (get_latitude_index_top() > get_latitude_index_bottom())
00294         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.");
00295 
00296     // Constrain the lat vector and lat dim of the array
00297     d_latitude->add_constraint(fd, get_latitude_index_top(), 1,
00298                                get_latitude_index_bottom());
00299     d_grid->get_array()->add_constraint(get_lat_dim(),
00300                                         get_latitude_index_top(), 1,
00301                                         get_latitude_index_bottom());
00302 
00303     // Does the longitude constraint cross the edge of the longitude vector?
00304     // If so, reorder the grid's data (array), longitude map vector and the
00305     // local vector of longitude data used for computation.
00306     if (get_longitude_index_left() > get_longitude_index_right()) {
00307         reorder_longitude_map(get_longitude_index_left());
00308 
00309         // If the longitude constraint is 'split', join the two parts, reload
00310         // the data into the Grid's Array and make sure the Array is marked as
00311         // already read. This should be true for the whole Grid, but if some
00312         // future modification changes that, the array will be covered here.
00313         // Note that the following method only reads the data out and stores
00314         // it in this object after joining the two parts. The method
00315         // apply_constraint_to_data() transfers the data back from the this
00316         // object to the DAP Grid variable.
00317         reorder_data_longitude_axis(*d_grid->get_array(), get_lon_dim());
00318 
00319         // Now that the data are all in local storage alter the indices; the
00320         // left index has now been moved to 0, and the right index is now
00321         // at lon_vector_length-left+right.
00322         set_longitude_index_right(get_lon_length() - get_longitude_index_left()
00323                                   + get_longitude_index_right());
00324         set_longitude_index_left(0);
00325     }
00326 
00327     // If the constraint used the -180/179 (neg_pos) notation, transform
00328     // the longitude map so it uses the -180/179 notation. Note that at this
00329     // point, d_longitude always uses the pos notation because of the earlier
00330     // conditional transformation.
00331 
00332     // Do this _before_ applying the constraint since set_array_using_double()
00333     // tests the array length using Vector::length() and that method returns
00334     // the length _as constrained_. We want to move all of the longitude
00335     // values from d_lon back into the map, not just the number that will be
00336     // sent (although an optimization might do this, it's hard to imagine
00337     // it would gain much).
00338     if (get_longitude_notation() == neg_pos) {
00339         transform_longitude_to_neg_pos_notation();
00340     }
00341 
00342     // Apply constraint; stride is always one and maps only have one dimension
00343     fd = d_longitude->dim_begin();
00344     d_longitude->add_constraint(fd, get_longitude_index_left(), 1,
00345                                 get_longitude_index_right());
00346 
00347     d_grid->get_array()->add_constraint(get_lon_dim(),
00348                                         get_longitude_index_left(),
00349                                         1, get_longitude_index_right());
00350 
00351     // Transfer values from the local lat vector to the Grid's
00352     // Here test the sense of the latitude vector and invert the vector if the
00353     // sense is 'inverted' so that the top is always the northern-most value
00354     if (get_latitude_sense() == inverted) {
00355         DBG(cerr << "Inverted latitude sense" << endl);
00356         transpose_vector(get_lat() + get_latitude_index_top(),
00357                 get_latitude_index_bottom() - get_latitude_index_top() + 1);
00358         // Now read the Array data and flip the latitudes.
00359         flip_latitude_within_array(*d_grid->get_array(),
00360                 get_latitude_index_bottom() - get_latitude_index_top() + 1,
00361                 get_longitude_index_right() - get_longitude_index_left() + 1);
00362     }
00363 
00364     set_array_using_double(d_latitude, get_lat() + get_latitude_index_top(),
00365                            get_latitude_index_bottom() - get_latitude_index_top() + 1);
00366 
00367     set_array_using_double(d_longitude, get_lon() + get_longitude_index_left(),
00368                            get_longitude_index_right() - get_longitude_index_left() + 1);
00369 
00370     // Look for any non-lat/lon maps and make sure they are read correctly
00371     Grid::Map_iter i = d_grid->map_begin();
00372     Grid::Map_iter end = d_grid->map_end();
00373     while (i != end) {
00374         if (*i != d_latitude && *i != d_longitude) {
00375             if ((*i)->send_p()) {
00376                 DBG(cerr << "reading grid map: " << (*i)->name() << endl);
00377                 //(*i)->set_read_p(false);
00378                 (*i)->read();
00379             }
00380         }
00381         ++i;
00382     }
00383 
00384     // ... and then the Grid's array if it has been read.
00385     if (get_array_data()) {
00386         int size = d_grid->get_array()->val2buf(get_array_data());
00387 
00388         if (size != get_array_data_size())
00389             throw InternalErr(__FILE__, __LINE__, "Expected data size not copied to the Grid's buffer.");
00390 
00391         d_grid->set_read_p(true);
00392     }
00393     else {
00394         d_grid->get_array()->read();
00395     }
00396 }
00397 
00398 } // namespace libdap