ArrayGeoConstraint.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: ArrayGeoConstraint.cc 16612 2007-06-04 22:08:57Z 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 "ArrayGeoConstraint.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 
00057 void ArrayGeoConstraint::m_init()
00058 {
00059     if (d_array->dimensions() < 2 || d_array->dimensions() > 3)
00060         throw Error("The geoarray() function works only with Arrays of two or three dimensions.");
00061 
00062     build_lat_lon_maps();
00063 }
00064 
00065 ArrayGeoConstraint::ArrayGeoConstraint(Array *array, const string &ds_name,
00066                                        double left, double top, double right, double bottom)
00067         : GeoConstraint(ds_name), d_array(array),
00068         d_extent(left, top, right, bottom), d_projection("plat-carre", "wgs84")
00069 
00070 {
00071     m_init();
00072 }
00073 
00074 ArrayGeoConstraint::ArrayGeoConstraint(Array *array, const string &ds_name,
00075                                        double left, double top, double right, double bottom,
00076                                        const string &projection, const string &datum)
00077         : GeoConstraint(ds_name), d_array(array),
00078         d_extent(left, top, right, bottom), d_projection(projection, datum)
00079 
00080 {
00081     m_init();
00082 }
00083 
00095 bool ArrayGeoConstraint::build_lat_lon_maps()
00096 {
00097     // Find the longitude dimension: Assume it is the rightmost
00098     set_longitude_rightmost(true);
00099     set_lon_dim(d_array->dim_begin() + (d_array->dimensions() - 1));
00100 
00101     int number_elements_longitude = d_array->dimension_size(get_lon_dim());
00102     double *lon_map = new double[number_elements_longitude];
00103     for (int i = 0; i < number_elements_longitude; ++i) {
00104         lon_map[i] = ((d_extent.d_right - d_extent.d_left) / (number_elements_longitude - 1)) * i + d_extent.d_left;
00105     }
00106     set_lon(lon_map);
00107     set_lon_length(number_elements_longitude);
00108 
00109     // Find the latitude dimension: Assume it is the next-rightmost
00110     set_lat_dim(d_array->dim_begin() + (d_array->dimensions() - 2));
00111 
00112     int number_elements_latitude = d_array->dimension_size(get_lat_dim());
00113     double *lat_map = new double[number_elements_latitude];
00114     for (int i = 0; i < number_elements_latitude; ++i) {
00115         lat_map[i] = ((d_extent.d_bottom - d_extent.d_top) / (number_elements_latitude - 1)) * i + d_extent.d_top;
00116     }
00117     set_lat(lat_map);
00118     set_lat_length(number_elements_latitude);
00119 
00120     return get_lat() && get_lon();
00121 }
00122 
00130 bool
00131 ArrayGeoConstraint::lat_lon_dimensions_ok()
00132 {
00133     return true;
00134 }
00135 
00150 void ArrayGeoConstraint::apply_constraint_to_data()
00151 {
00152     if (!get_bounding_box_set())
00153         throw InternalErr(
00154             "The Latitude and Longitude constraints must be set before calling\n\
00155             apply_constraint_to_data().");
00156 
00157     if (get_latitude_sense() == inverted) {
00158         int tmp = get_latitude_index_top();
00159         set_latitude_index_top(get_latitude_index_bottom());
00160         set_latitude_index_bottom(tmp);
00161     }
00162 
00163     // It's easy to flip the Latitude values; if the bottom index value
00164     // is before/above the top index, return an error explaining that.
00165     if (get_latitude_index_top() > get_latitude_index_bottom())
00166         throw Error("The upper and lower latitude indexes appear to be reversed. Please provide\nthe latitude bounding box numbers giving the northern-most latitude first.");
00167 
00168     d_array->add_constraint(get_lat_dim(),
00169                             get_latitude_index_top(), 1,
00170                             get_latitude_index_bottom());
00171 
00172     // Does the longitude constraint cross the edge of the longitude vector?
00173     // If so, reorder the data (array).
00174     if (get_longitude_index_left() > get_longitude_index_right()) {
00175         reorder_data_longitude_axis(*d_array);
00176 
00177         // Now the data are all in local storage
00178 
00179         // alter the indexes; the left index has now been moved to 0, and the right
00180         // index is now at lon_vector_length-left+right.
00181         set_longitude_index_right(get_lon_length() - get_longitude_index_left()
00182                                   + get_longitude_index_right());
00183         set_longitude_index_left(0);
00184     }
00185     // If the constraint used the -180/179 (neg_pos) notation, transform
00186     // the longitude map s it uses the -180/179 notation. Note that at this
00187     // point, d_longitude always uses the pos notation because of the earlier
00188     // conditional transformation.
00189 
00190     // Apply constraint; stride is always one
00191     d_array->add_constraint(get_lon_dim(),
00192                             get_longitude_index_left(), 1,
00193                             get_longitude_index_right());
00194 
00195     // Load the array if it has been read, which will be the case if
00196     // reorder_data_longitude_axis() has been called.
00197     if (get_array_data()) {
00198         int size = d_array->val2buf(get_array_data());
00199         if (size != get_array_data_size())
00200             throw InternalErr
00201             ("Expected data size not copied to the Grid's buffer.");
00202         d_array->set_read_p(true);
00203     }
00204     else {
00205         d_array->read(get_dataset());
00206     }
00207 }
00208 

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