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