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) 2006 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: GeoConstraint.cc 19914 2008-11-25 17:26:22Z jimg $" 00033 }; 00034 00035 #include <cstring> 00036 #include <cmath> 00037 #include <iostream> 00038 #include <sstream> 00039 #include <algorithm> // for find_if 00040 00041 //#define DODS_DEBUG2 00042 00043 #include "debug.h" 00044 #include "dods-datatypes.h" 00045 #include "GeoConstraint.h" 00046 #include "Float64.h" 00047 00048 #include "Error.h" 00049 #include "InternalErr.h" 00050 #include "ce_functions.h" 00051 #include "util.h" 00052 00053 using namespace std; 00054 00055 namespace libdap { 00056 00061 class is_prefix 00062 { 00063 private: 00064 string s; 00065 public: 00066 is_prefix(const string & in): s(in) 00067 {} 00068 00069 bool operator()(const string & prefix) 00070 { 00071 return s.find(prefix) == 0; 00072 } 00073 }; 00074 00085 bool 00086 unit_or_name_match(set < string > units, set < string > names, 00087 const string & var_units, const string & var_name) 00088 { 00089 return (units.find(var_units) != units.end() 00090 || find_if(names.begin(), names.end(), 00091 is_prefix(var_name)) != names.end()); 00092 } 00093 00108 GeoConstraint::Notation 00109 GeoConstraint::categorize_notation(double left, 00110 double right) const 00111 { 00112 return (left < 0 || right < 0) ? neg_pos : pos; 00113 } 00114 00115 /* A private method to translate the longitude constraint from -180/179 00116 notation to 0/359 notation. 00117 @param left Value-result parameter; the left side of the bounding box 00118 @parm right Value-result parameter; the right side of the bounding box */ 00119 void 00120 GeoConstraint::transform_constraint_to_pos_notation(double &left, 00121 double &right) const 00122 { 00123 if (left < 0 || right < 0) { // double check 00124 left += 180; 00125 right += 180; 00126 } 00127 } 00128 00132 void GeoConstraint::transform_longitude_to_pos_notation() 00133 { 00134 // Assume earlier logic is correct (since the test is expensive) 00135 // for each value, add 180 00136 // Longitude could be represented using any of the numeric types 00137 for (int i = 0; i < d_lon_length; ++i) 00138 d_lon[i] += 180; 00139 } 00140 00144 void GeoConstraint::transform_longitude_to_neg_pos_notation() 00145 { 00146 for (int i = 0; i < d_lon_length; ++i) 00147 d_lon[i] -= 180; 00148 } 00149 00150 bool GeoConstraint::is_bounding_box_valid(double left, double top, 00151 double right, double bottom) const 00152 { 00153 if ((left < d_lon[0] && right < d_lon[0]) 00154 || (left > d_lon[d_lon_length-1] && right > d_lon[d_lon_length-1])) 00155 return false; 00156 00157 if (d_latitude_sense == normal) { 00158 // When sense is normal, the largest values are at the origin. 00159 if ((top > d_lat[0] && bottom > d_lat[0]) 00160 || (top < d_lat[d_lat_length-1] && bottom < d_lat[d_lat_length-1])) 00161 return false; 00162 } 00163 else { 00164 if ((top < d_lat[0] && bottom < d_lat[0]) 00165 || (top > d_lat[d_lat_length-1] && bottom > d_lat[d_lat_length-1])) 00166 return false; 00167 } 00168 00169 return true; 00170 } 00171 00182 void GeoConstraint::find_longitude_indeces(double left, double right, 00183 int &longitude_index_left, 00184 int &longitude_index_right) const 00185 { 00186 // Use all values 'modulo 360' to take into account the cases where the 00187 // constraint and/or the longitude vector are given using values greater 00188 // than 360 (i.e., when 380 is used to mean 20). 00189 double t_left = fmod(left, 360.0); 00190 double t_right = fmod(right, 360.0); 00191 00192 // Find the place where 'longitude starts.' That is, what value of the 00193 // index 'i' corresponds to the smallest value of d_lon. Why we do this: 00194 // Some data sources use offset longitude axes so that the 'seam' is 00195 // shifted to a place other than the date line. 00196 int i = 0; 00197 int smallest_lon_index = 0; 00198 double smallest_lon = fmod(d_lon[smallest_lon_index], 360.0); 00199 while (i < d_lon_length - 1) { 00200 if (smallest_lon > fmod(d_lon[i], 360.0)) { 00201 smallest_lon_index = i; 00202 smallest_lon = fmod(d_lon[smallest_lon_index], 360.0); 00203 } 00204 ++i; 00205 } 00206 DBG2(cerr << "smallest_lon_index: " << smallest_lon_index << endl); 00207 00208 // Scan from the index of the smallest value looking for the place where 00209 // the value is greater than or equal to the left most point of the bounding 00210 // box. 00211 i = smallest_lon_index; 00212 bool done = false; 00213 00214 DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) < t_left: " 00215 << fmod(d_lon[i], 360.0) << " < " << t_left << endl); 00216 00217 while (!done && fmod(d_lon[i], 360.0) < t_left) { 00218 00219 DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) < t_left: " 00220 << fmod(d_lon[i], 360.0) << " < " << t_left << endl); 00221 00222 ++i; 00223 i = i % d_lon_length; 00224 if (i == smallest_lon_index) 00225 done = true; 00226 } 00227 if (fmod(d_lon[i], 360.0) == t_left) 00228 longitude_index_left = i; 00229 else 00230 longitude_index_left = (i - 1) > 0 ? i - 1 : 0; 00231 00232 // Assume the vector is cirular --> the largest value is next to the 00233 // smallest. 00234 int largest_lon_index = 00235 (smallest_lon_index - 1 + d_lon_length) % d_lon_length; 00236 DBG2(cerr << "largest_lon_index: " << largest_lon_index << endl); 00237 i = largest_lon_index; 00238 done = false; 00239 00240 DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) > t_right: " 00241 << fmod(d_lon[i], 360.0) << " > " << t_right << endl); 00242 00243 while (!done && fmod(d_lon[i], 360.0) > t_right) { 00244 00245 DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) > t_right: " 00246 << fmod(d_lon[i], 360.0) << " > " << t_right << endl); 00247 00248 // This is like modulus but for 'counting down' 00249 i = (i == 0) ? d_lon_length - 1 : i - 1; 00250 if (i == largest_lon_index) 00251 done = true; 00252 } 00253 if (fmod(d_lon[i], 360.0) == t_right) 00254 longitude_index_right = i; 00255 else 00256 longitude_index_right = 00257 (i + 1) < d_lon_length - 1 ? i + 1 : d_lon_length - 1; 00258 } 00259 00272 void GeoConstraint::find_latitude_indeces(double top, double bottom, 00273 LatitudeSense sense, 00274 int &latitude_index_top, 00275 int &latitude_index_bottom) const 00276 { 00277 // Scan from the top down 00278 int i = 0; 00279 if (sense == normal) { 00280 while (i < d_lat_length - 1 && d_lat[i] > top) 00281 ++i; 00282 if (d_lat[i] == top) 00283 latitude_index_top = i; 00284 else 00285 latitude_index_top = (i - 1) > 0 ? i - 1 : 0; 00286 } 00287 else { 00288 while (i < d_lat_length - 1 && d_lat[i] < top) 00289 ++i; 00290 latitude_index_top = i; 00291 } 00292 00293 00294 // and from the bottom up 00295 i = d_lat_length - 1; 00296 if (sense == normal) { 00297 while (i > 0 && d_lat[i] < bottom) 00298 --i; 00299 if (d_lat[i] == bottom) 00300 latitude_index_bottom = i; 00301 else 00302 latitude_index_bottom = 00303 (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1; 00304 } 00305 else { 00306 while (i > 0 && d_lat[i] > bottom) 00307 --i; 00308 latitude_index_bottom = i; 00309 } 00310 } 00311 00315 GeoConstraint::LatitudeSense GeoConstraint::categorize_latitude() const 00316 { 00317 return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted; 00318 } 00319 00320 #if 0 00321 00328 void GeoConstraint::set_bounding_box_latitude(double top, double bottom) 00329 { 00330 // Record the 'sense' of the latitude for use here and later on. 00331 d_latitude_sense = categorize_latitude(); 00332 00333 // This is simpler than the longitude case because there's no need to test 00334 // for several notations, no need to accommodate them in the return, no 00335 // modulo arithmetic for the axis and no need to account for a constraint with 00336 // two disconnected parts to be joined. 00337 find_latitude_indeces(top, bottom, d_latitude_sense, 00338 d_latitude_index_top, d_latitude_index_bottom); 00339 } 00340 #endif 00341 00342 // Use 'index' as the pivot point. Move the points behind index to the front of 00343 // the vector and those points in front of and at index to the rear. 00344 static void 00345 swap_vector_ends(char *dest, char *src, int len, int index, int elem_sz) 00346 { 00347 memcpy(dest, src + index * elem_sz, (len - index) * elem_sz); 00348 00349 memcpy(dest + (len - index) * elem_sz, src, index * elem_sz); 00350 } 00351 00360 void GeoConstraint::reorder_longitude_map(int longitude_index_left) 00361 { 00362 double *tmp_lon = new double[d_lon_length]; 00363 00364 swap_vector_ends((char *) tmp_lon, (char *) d_lon, d_lon_length, 00365 longitude_index_left, sizeof(double)); 00366 00367 memcpy(d_lon, tmp_lon, d_lon_length * sizeof(double)); 00368 00369 delete[]tmp_lon; 00370 } 00371 00372 static int 00373 count_dimensions_except_longitude(Array &a) 00374 { 00375 int size = 1; 00376 for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i) 00377 size *= a.dimension_size(i, true); 00378 00379 return size; 00380 } 00381 00400 void GeoConstraint::reorder_data_longitude_axis(Array &a) 00401 { 00402 00403 if (!get_longitude_rightmost()) 00404 throw Error("This grid does not have Longitude as its rightmost dimension, the geogrid()\ndoes not support constraints that wrap around the edges of this type of grid."); 00405 00406 DBG(cerr << "Constraint for the left half: " << get_longitude_index_left() 00407 << ", " << get_lon_length() - 1 << endl); 00408 00409 a.add_constraint(d_lon_dim, get_longitude_index_left(), 1, 00410 get_lon_length() - 1); 00411 a.set_read_p(false); 00412 a.read(); 00413 DBG2(a.print_val(stderr)); 00414 #if 0 00415 char *left_data = 0; 00416 int left_size = a.buf2val((void **) & left_data); 00417 #endif 00418 char *left_data = (char*)a.value(); 00419 int left_size = a.length(); 00420 00421 // Build a constraint for the 'right' part, which 00422 // goes from the left edge of the array to the right index and read those 00423 // data. 00424 //a.clear_constraint(); 00425 00426 DBG(cerr << "Constraint for the right half: " << 0 00427 << ", " << get_longitude_index_right() << endl); 00428 00429 a.add_constraint(d_lon_dim, 0, 1, get_longitude_index_right()); 00430 a.set_read_p(false); 00431 a.read(); 00432 DBG2(a.print_val(stderr)); 00433 #if 0 00434 char *right_data = 0; 00435 int right_size = a.buf2val((void **) & right_data); 00436 #endif 00437 char *right_data = (char*)a.value(); 00438 int right_size = a.length(); 00439 00440 // Make one big lump O'data 00441 d_array_data_size = left_size + right_size; 00442 d_array_data = new char[d_array_data_size]; 00443 00444 // Assume COARDS conventions are being followed: lon varies fastest. 00445 // These *_elements variables are actually elements * bytes/element since 00446 // memcpy() uses bytes. 00447 int elem_width = a.var()->width(); 00448 int left_elements = (get_lon_length() - get_longitude_index_left()) * elem_width; 00449 int right_elements = (get_longitude_index_right() + 1) * elem_width; 00450 int total_elements_per_row = left_elements + right_elements; 00451 00452 // Interleave the left and right_data vectors. jhrg 8/31/06 00453 int rows_to_copy = count_dimensions_except_longitude(a); 00454 for (int i = 0; i < rows_to_copy; ++i) { 00455 memcpy(d_array_data + (total_elements_per_row * i), 00456 left_data + (left_elements * i), 00457 left_elements); 00458 memcpy(d_array_data + left_elements + (total_elements_per_row * i), 00459 right_data + (right_elements * i), 00460 right_elements); 00461 } 00462 00463 delete[]left_data; 00464 delete[]right_data; 00465 } 00466 00467 #if 0 00468 00481 void GeoConstraint::set_bounding_box_longitude(double left, double right) 00482 { 00483 // Categorize the notation used by the bounding box (0/359 or -180/179). 00484 d_longitude_notation = categorize_notation(left, right); 00485 00486 // If the notation uses -180/179, transform the request to 0/359 notation. 00487 if (d_longitude_notation == neg_pos) 00488 transform_constraint_to_pos_notation(left, right); 00489 00490 // If the grid uses -180/179, transform it to 0/359 as well. This will make 00491 // subsequent logic easier and adds only a few extra operations, even with 00492 // large maps. 00493 Notation longitude_notation = 00494 categorize_notation(d_lon[0], d_lon[d_lon_length - 1]); 00495 00496 if (longitude_notation == neg_pos) 00497 transform_longitude_to_pos_notation(); 00498 00499 // Find the longitude map indeces that correspond to the bounding box. 00500 find_longitude_indeces(left, right, d_longitude_index_left, 00501 d_longitude_index_right); 00502 } 00503 #endif 00504 00510 GeoConstraint::GeoConstraint() 00511 : d_array_data(0), d_array_data_size(0), 00512 d_lat(0), d_lon(0), 00513 d_bounding_box_set(false), 00514 d_longitude_notation(unknown_notation), 00515 d_latitude_sense(unknown_sense) 00516 { 00517 // Build sets of attribute values for easy searching. Maybe overkill??? 00518 d_coards_lat_units.insert("degrees_north"); 00519 d_coards_lat_units.insert("degree_north"); 00520 d_coards_lat_units.insert("degree_N"); 00521 d_coards_lat_units.insert("degrees_N"); 00522 00523 d_coards_lon_units.insert("degrees_east"); 00524 d_coards_lon_units.insert("degree_east"); 00525 d_coards_lon_units.insert("degrees_E"); 00526 d_coards_lon_units.insert("degree_E"); 00527 00528 d_lat_names.insert("COADSY"); 00529 d_lat_names.insert("lat"); 00530 d_lat_names.insert("Lat"); 00531 d_lat_names.insert("LAT"); 00532 00533 d_lon_names.insert("COADSX"); 00534 d_lon_names.insert("lon"); 00535 d_lon_names.insert("Lon"); 00536 d_lon_names.insert("LON"); 00537 } 00538 00549 void GeoConstraint::set_bounding_box(double left, double top, 00550 double right, double bottom) 00551 { 00552 // Ensure this method is called only once. What about pthreads? 00553 // The method Array::reset_constraint() might make this so it could be 00554 // called more than once. jhrg 8/30/06 00555 if (d_bounding_box_set) 00556 throw 00557 InternalErr 00558 ("It is not possible to register more than one geographical constraint on a variable."); 00559 00560 // Record the 'sense' of the latitude for use here and later on. 00561 d_latitude_sense = categorize_latitude(); 00562 00563 // Categorize the notation used by the bounding box (0/359 or -180/179). 00564 d_longitude_notation = categorize_notation(left, right); 00565 00566 // If the notation uses -180/179, transform the request to 0/359 notation. 00567 if (d_longitude_notation == neg_pos) 00568 transform_constraint_to_pos_notation(left, right); 00569 00570 // If the grid uses -180/179, transform it to 0/359 as well. This will make 00571 // subsequent logic easier and adds only a few extra operations, even with 00572 // large maps. 00573 Notation longitude_notation = 00574 categorize_notation(d_lon[0], d_lon[d_lon_length - 1]); 00575 00576 if (longitude_notation == neg_pos) 00577 transform_longitude_to_pos_notation(); 00578 00579 if (!is_bounding_box_valid(left, top, right, bottom)) 00580 throw Error("The bounding box does not intersect any data within this Grid or Array. The\ngeographical extent of these data are from latitude " 00581 + double_to_string(d_lat[0]) + " to " 00582 + double_to_string(d_lat[d_lat_length-1]) 00583 + "\nand longitude " + double_to_string(d_lon[0]) 00584 + " to " + double_to_string(d_lon[d_lon_length-1]) 00585 + " while the bounding box provided was latitude " 00586 + double_to_string(top) + " to " 00587 + double_to_string(bottom) + "\nand longitude " 00588 + double_to_string(left) + " to " 00589 + double_to_string(right)); 00590 00591 // This is simpler than the longitude case because there's no need to 00592 // test for several notations, no need to accommodate them in the return, 00593 // no modulo arithmetic for the axis and no need to account for a 00594 // constraint with two disconnected parts to be joined. 00595 find_latitude_indeces(top, bottom, d_latitude_sense, 00596 d_latitude_index_top, d_latitude_index_bottom); 00597 00598 00599 // Find the longitude map indeces that correspond to the bounding box. 00600 find_longitude_indeces(left, right, d_longitude_index_left, 00601 d_longitude_index_right); 00602 00603 DBG(cerr << "Bounding box (tlbr): " << d_latitude_index_top << ", " 00604 << d_longitude_index_left << ", " 00605 << d_latitude_index_bottom << ", " 00606 << d_longitude_index_right << endl); 00607 00608 d_bounding_box_set = true; 00609 } 00610 00611 } // namespace libdap