libdap++ Updated for version 3.8.2
|
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 21922 2010-01-07 18:23:57Z jimg $" 00033 }; 00034 00035 #include <cmath> 00036 #include <iostream> 00037 #include <sstream> 00038 00039 //#define DODS_DEBUG 00040 00041 #include "debug.h" 00042 #include "dods-datatypes.h" 00043 #include "ArrayGeoConstraint.h" 00044 #include "Float64.h" 00045 00046 #include "Error.h" 00047 #include "InternalErr.h" 00048 #include "ce_functions.h" 00049 00050 using namespace std; 00051 00052 namespace libdap { 00053 00055 void ArrayGeoConstraint::m_init() 00056 { 00057 if (d_array->dimensions() < 2 || d_array->dimensions() > 3) 00058 throw Error("The geoarray() function works only with Arrays of two or three dimensions."); 00059 00060 build_lat_lon_maps(); 00061 } 00062 00063 // In the other methods the ds_name parameter defaults to "" but that's not 00064 // possible here. Remove ds_name 00065 ArrayGeoConstraint::ArrayGeoConstraint(Array *array, double top, double left, 00066 double bottom, double right) 00067 : GeoConstraint(), d_array(array), 00068 d_extent(top, left, bottom, right), d_projection("plat-carre", "wgs84") 00069 00070 { 00071 m_init(); 00072 } 00073 00074 ArrayGeoConstraint::ArrayGeoConstraint(Array *array, 00075 double top, double left, double bottom, double right, 00076 const string &projection, const string &datum) 00077 : GeoConstraint(), d_array(array), 00078 d_extent(top, left, bottom, right), 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 (!is_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, get_lon_dim()); 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 00199 int size = d_array->val2buf(get_array_data()); 00200 00201 if (size != get_array_data_size()) 00202 throw InternalErr 00203 ("Expected data size not copied to the Grid's buffer."); 00204 d_array->set_read_p(true); 00205 } 00206 else { 00207 d_array->read(); 00208 } 00209 } 00210 00211 } // namespace libdap 00212