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 23563 2010-09-14 17:43:32Z 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_DEBUG 00042 //#define DODS_DEBUG2 00043 00044 #include "debug.h" 00045 #include "dods-datatypes.h" 00046 #include "GeoConstraint.h" 00047 #include "Float64.h" 00048 00049 #include "Error.h" 00050 #include "InternalErr.h" 00051 #include "ce_functions.h" 00052 #include "util.h" 00053 00054 using namespace std; 00055 00056 namespace libdap { 00057 00062 class is_prefix 00063 { 00064 private: 00065 string s; 00066 public: 00067 is_prefix(const string & in): s(in) 00068 {} 00069 00070 bool operator()(const string & prefix) 00071 { 00072 return s.find(prefix) == 0; 00073 } 00074 }; 00075 00086 bool 00087 unit_or_name_match(set < string > units, set < string > names, 00088 const string & var_units, const string & var_name) 00089 { 00090 return (units.find(var_units) != units.end() 00091 || find_if(names.begin(), names.end(), 00092 is_prefix(var_name)) != names.end()); 00093 } 00094 00109 GeoConstraint::Notation 00110 GeoConstraint::categorize_notation(const double left, 00111 const double right) const 00112 { 00113 return (left < 0 || right < 0) ? neg_pos : pos; 00114 } 00115 00116 /* A private method to translate the longitude constraint from -180/179 00117 notation to 0/359 notation. 00118 00119 About the two notations: 00120 0 180 360 00121 |---------------------------|-------------------------| 00122 0 180,-180 0 00123 00124 so in the neg-pos notation (using the name I give it in this class) both 00125 180 and -180 are the same longitude. And in the pos notation 0 and 360 are 00126 the same. 00127 00128 @param left Value-result parameter; the left side of the bounding box 00129 @parm right Value-result parameter; the right side of the bounding box */ 00130 void 00131 GeoConstraint::transform_constraint_to_pos_notation(double &left, 00132 double &right) const 00133 { 00134 if (left < 0) 00135 left += 360 ; 00136 00137 if (right < 0) 00138 right += 360; 00139 } 00140 00149 void GeoConstraint::transform_longitude_to_pos_notation() 00150 { 00151 // Assume earlier logic is correct (since the test is expensive) 00152 // for each value, add 180 00153 // Longitude could be represented using any of the numeric types 00154 for (int i = 0; i < d_lon_length; ++i) 00155 if (d_lon[i] < 0) 00156 d_lon[i] += 360; 00157 } 00158 00168 void GeoConstraint::transform_longitude_to_neg_pos_notation() 00169 { 00170 for (int i = 0; i < d_lon_length; ++i) 00171 if (d_lon[i] > 180) 00172 d_lon[i] -= 360; 00173 } 00174 00175 bool GeoConstraint::is_bounding_box_valid(const double left, const double top, 00176 const double right, const double bottom) const 00177 { 00178 if ((left < d_lon[0] && right < d_lon[0]) 00179 || (left > d_lon[d_lon_length-1] && right > d_lon[d_lon_length-1])) 00180 return false; 00181 00182 if (d_latitude_sense == normal) { 00183 // When sense is normal, the largest values are at the origin. 00184 if ((top > d_lat[0] && bottom > d_lat[0]) 00185 || (top < d_lat[d_lat_length-1] && bottom < d_lat[d_lat_length-1])) 00186 return false; 00187 } 00188 else { 00189 if ((top < d_lat[0] && bottom < d_lat[0]) 00190 || (top > d_lat[d_lat_length-1] && bottom > d_lat[d_lat_length-1])) 00191 return false; 00192 } 00193 00194 return true; 00195 } 00196 00207 void GeoConstraint::find_longitude_indeces(double left, double right, 00208 int &longitude_index_left, int &longitude_index_right) const 00209 { 00210 // Use all values 'modulo 360' to take into account the cases where the 00211 // constraint and/or the longitude vector are given using values greater 00212 // than 360 (i.e., when 380 is used to mean 20). 00213 double t_left = fmod(left, 360.0); 00214 double t_right = fmod(right, 360.0); 00215 00216 // Find the place where 'longitude starts.' That is, what value of the 00217 // index 'i' corresponds to the smallest value of d_lon. Why we do this: 00218 // Some data sources use offset longitude axes so that the 'seam' is 00219 // shifted to a place other than the date line. 00220 int i = 0; 00221 int lon_origin_index = 0; 00222 double smallest_lon = fmod(d_lon[0], 360.0); 00223 while (i < d_lon_length) { 00224 double curent_lon_value = fmod(d_lon[i], 360.0); 00225 if (smallest_lon > curent_lon_value) { 00226 smallest_lon = curent_lon_value; 00227 lon_origin_index = i; 00228 } 00229 ++i; 00230 } 00231 00232 DBG2(cerr << "lon_origin_index: " << lon_origin_index << endl); 00233 00234 // Scan from the index of the smallest value looking for the place where 00235 // the value is greater than or equal to the left most point of the bounding 00236 // box. 00237 i = lon_origin_index; 00238 while (fmod(d_lon[i], 360.0) < t_left) { 00239 ++i; 00240 i = i % d_lon_length; 00241 00242 // If we cycle completely through all the values/indices, throw 00243 if (i == lon_origin_index) 00244 throw Error("geogrid: Could not find an index for the longitude value '" + double_to_string(left) + "'"); 00245 } 00246 00247 if (fmod(d_lon[i], 360.0) == t_left) 00248 longitude_index_left = i; 00249 else 00250 longitude_index_left = (i - 1) > 0 ? i - 1 : 0; 00251 00252 DBG2(cerr << "longitude_index_left: " << longitude_index_left << endl); 00253 00254 // Assume the vector is circular --> the largest value is next to the 00255 // smallest. 00256 int largest_lon_index = (lon_origin_index - 1 + d_lon_length) % d_lon_length; 00257 i = largest_lon_index; 00258 while (fmod(d_lon[i], 360.0) > t_right) { 00259 // This is like modulus but for 'counting down' 00260 i = (i == 0) ? d_lon_length - 1 : i - 1; 00261 if (i == largest_lon_index) 00262 throw Error("geogrid: Could not find an index for the longitude value '" + double_to_string(right) + "'"); 00263 } 00264 00265 if (fmod(d_lon[i], 360.0) == t_right) 00266 longitude_index_right = i; 00267 else 00268 longitude_index_right = (i + 1) < d_lon_length - 1 ? i + 1 : d_lon_length - 1; 00269 00270 DBG2(cerr << "longitude_index_right: " << longitude_index_right << endl); 00271 } 00272 00285 void GeoConstraint::find_latitude_indeces(double top, double bottom, 00286 LatitudeSense sense, 00287 int &latitude_index_top, 00288 int &latitude_index_bottom) const 00289 { 00290 int i, j; 00291 00292 if (sense == normal) { 00293 i = 0; 00294 while (i < d_lat_length - 1 && top < d_lat[i]) 00295 ++i; 00296 00297 j = d_lat_length - 1; 00298 while (j > 0 && bottom > d_lat[j]) 00299 --j; 00300 00301 if (d_lat[i] == top) 00302 latitude_index_top = i; 00303 else 00304 latitude_index_top = (i - 1) > 0 ? i - 1 : 0; 00305 00306 if (d_lat[j] == bottom) 00307 latitude_index_bottom = j; 00308 else 00309 latitude_index_bottom = 00310 (j + 1) < d_lat_length - 1 ? j + 1 : d_lat_length - 1; 00311 } 00312 else { 00313 i = d_lat_length - 1; 00314 while (i > 0 && d_lat[i] > top) 00315 --i; 00316 00317 j = 0; 00318 while (j < d_lat_length - 1 && d_lat[j] < bottom) 00319 ++j; 00320 00321 if (d_lat[i] == top) 00322 latitude_index_top = i; 00323 else 00324 latitude_index_top = (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1; 00325 00326 if (d_lat[j] == bottom) 00327 latitude_index_bottom = j; 00328 else 00329 latitude_index_bottom = (j - 1) > 0 ? j - 1 : 0; 00330 } 00331 } 00332 00336 GeoConstraint::LatitudeSense GeoConstraint::categorize_latitude() const 00337 { 00338 return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted; 00339 } 00340 00341 // Use 'index' as the pivot point. Move the points behind index to the front of 00342 // the vector and those points in front of and at index to the rear. 00343 static void 00344 swap_vector_ends(char *dest, char *src, int len, int index, int elem_sz) 00345 { 00346 memcpy(dest, src + index * elem_sz, (len - index) * elem_sz); 00347 00348 memcpy(dest + (len - index) * elem_sz, src, index * elem_sz); 00349 } 00350 00351 template<class T> 00352 static void transpose(std::vector<std::vector<T> > a, 00353 std::vector<std::vector<T> > b, int width, int height) 00354 { 00355 for (int i = 0; i < width; i++) { 00356 for (int j = 0; j < height; j++) { 00357 b[j][i] = a[i][j]; 00358 } 00359 } 00360 } 00361 00369 void GeoConstraint::transpose_vector(double *src, const int length) 00370 { 00371 double *tmp = new double[length]; 00372 00373 int i = 0, j = length-1; 00374 while (i < length) 00375 tmp[j--] = src[i++]; 00376 00377 memcpy(src, tmp,length * sizeof(double)); 00378 00379 delete[] tmp; 00380 } 00381 00382 static int 00383 count_size_except_latitude_and_longitude(Array &a) 00384 { 00385 if (a.dim_end() - a.dim_begin() <= 2) // < 2 is really an error... 00386 return 1; 00387 00388 int size = 1; 00389 for (Array::Dim_iter i = a.dim_begin(); (i + 2) != a.dim_end(); ++i) 00390 size *= a.dimension_size(i, true); 00391 00392 return size; 00393 } 00394 00395 void GeoConstraint::flip_latitude_within_array(Array &a, int lat_length, 00396 int lon_length) 00397 { 00398 if (!d_array_data) { 00399 a.read(); 00400 d_array_data = static_cast<char*>(a.value()); 00401 d_array_data_size = a.width(); // Bytes not elements 00402 } 00403 00404 int size = count_size_except_latitude_and_longitude(a); 00405 char *tmp_data = new char[d_array_data_size]; 00406 int array_elem_size = a.var()->width(); 00407 int lat_lon_size = (d_array_data_size / size); 00408 00409 DBG(cerr << "lat, lon_length: " << lat_length << ", " << lon_length << endl); 00410 DBG(cerr << "size: " << size << endl); 00411 DBG(cerr << "d_array_data_size: " << d_array_data_size << endl); 00412 DBG(cerr << "array_elem_size: " << array_elem_size<< endl); 00413 DBG(cerr << "lat_lon_size: " << lat_lon_size<< endl); 00414 00415 for (int i = 0; i < size; ++i) { 00416 int lat = 0; 00417 int s_lat = lat_length - 1; 00418 // lon_length is the element size; memcpy() needs the number of bytes 00419 int lon_size = array_elem_size * lon_length; 00420 int offset = i * lat_lon_size; 00421 while (s_lat > -1) 00422 memcpy(tmp_data + offset + (lat++ * lon_size), 00423 d_array_data + offset + (s_lat-- * lon_size), 00424 lon_size); 00425 } 00426 00427 memcpy(d_array_data, tmp_data, d_array_data_size); 00428 delete [] tmp_data; 00429 } 00430 00439 void GeoConstraint::reorder_longitude_map(int longitude_index_left) 00440 { 00441 double *tmp_lon = new double[d_lon_length]; 00442 00443 swap_vector_ends((char *) tmp_lon, (char *) d_lon, d_lon_length, 00444 longitude_index_left, sizeof(double)); 00445 00446 memcpy(d_lon, tmp_lon, d_lon_length * sizeof(double)); 00447 00448 delete[]tmp_lon; 00449 } 00450 00451 static int 00452 count_dimensions_except_longitude(Array &a) 00453 { 00454 int size = 1; 00455 for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i) 00456 size *= a.dimension_size(i, true); 00457 00458 return size; 00459 } 00460 00478 void GeoConstraint::reorder_data_longitude_axis(Array &a, Array::Dim_iter lon_dim) 00479 { 00480 00481 if (!is_longitude_rightmost()) 00482 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."); 00483 00484 DBG(cerr << "Constraint for the left half: " << get_longitude_index_left() 00485 << ", " << get_lon_length() - 1 << endl); 00486 00487 // Build a constraint for the left part and get those values 00488 a.add_constraint(lon_dim, get_longitude_index_left(), 1, 00489 get_lon_length() - 1); 00490 a.set_read_p(false); 00491 a.read(); 00492 DBG2(a.print_val(stderr)); 00493 00494 // Save the left-hand data to local storage 00495 int left_size = a.width(); // width() == length() * element size 00496 char *left_data = (char*)a.value(); // value() allocates and copies 00497 00498 // Build a constraint for the 'right' part, which goes from the left edge 00499 // of the array to the right index and read those data. 00500 // (Don't call a.clear_constraint() since that will clear the constraint on 00501 // all the dimensions while add_constraint() will replace a constraint on 00502 // the given dimension). 00503 00504 DBG(cerr << "Constraint for the right half: " << 0 00505 << ", " << get_longitude_index_right() << endl); 00506 00507 a.add_constraint(lon_dim, 0, 1, get_longitude_index_right()); 00508 a.set_read_p(false); 00509 a.read(); 00510 DBG2(a.print_val(stderr)); 00511 00512 char *right_data = (char*)a.value(); 00513 int right_size = a.width(); 00514 00515 // Make one big lump O'data 00516 d_array_data_size = left_size + right_size; 00517 d_array_data = new char[d_array_data_size]; 00518 00519 // Assume COARDS conventions are being followed: lon varies fastest. 00520 // These *_elements variables are actually elements * bytes/element since 00521 // memcpy() uses bytes. 00522 int elem_size = a.var()->width(); 00523 int left_row_size = (get_lon_length() - get_longitude_index_left()) * elem_size; 00524 int right_row_size = (get_longitude_index_right() + 1) * elem_size; 00525 int total_bytes_per_row = left_row_size + right_row_size; 00526 00527 DBG2(cerr << "elem_size: " << elem_size << "; left & right size: " 00528 << left_row_size << ", " << right_row_size << endl); 00529 00530 // This will work for any number of dimension so long as longitude is the 00531 // right-most array dimension. 00532 int rows_to_copy = count_dimensions_except_longitude(a); 00533 for (int i = 0; i < rows_to_copy; ++i) { 00534 DBG(cerr << "Copying " << i << "th row" << endl); 00535 DBG(cerr << "left memcpy: " << *(float *)(left_data + (left_row_size * i)) << endl); 00536 00537 memcpy(d_array_data + (total_bytes_per_row * i), 00538 left_data + (left_row_size * i), 00539 left_row_size); 00540 00541 DBG(cerr << "right memcpy: " << *(float *)(right_data + (right_row_size * i)) << endl); 00542 00543 memcpy(d_array_data + (total_bytes_per_row * i) + left_row_size, 00544 right_data + (right_row_size * i), 00545 right_row_size); 00546 } 00547 00548 delete[]left_data; 00549 delete[]right_data; 00550 } 00551 00554 GeoConstraint::GeoConstraint() 00555 : d_array_data(0), d_array_data_size(0), 00556 d_lat(0), d_lon(0), 00557 d_bounding_box_set(false), 00558 d_longitude_rightmost(false), 00559 d_longitude_notation(unknown_notation), 00560 d_latitude_sense(unknown_sense) 00561 { 00562 // Build sets of attribute values for easy searching. Maybe overkill??? 00563 d_coards_lat_units.insert("degrees_north"); 00564 d_coards_lat_units.insert("degree_north"); 00565 d_coards_lat_units.insert("degree_N"); 00566 d_coards_lat_units.insert("degrees_N"); 00567 00568 d_coards_lon_units.insert("degrees_east"); 00569 d_coards_lon_units.insert("degree_east"); 00570 d_coards_lon_units.insert("degrees_E"); 00571 d_coards_lon_units.insert("degree_E"); 00572 00573 d_lat_names.insert("COADSY"); 00574 d_lat_names.insert("lat"); 00575 d_lat_names.insert("Lat"); 00576 d_lat_names.insert("LAT"); 00577 00578 d_lon_names.insert("COADSX"); 00579 d_lon_names.insert("lon"); 00580 d_lon_names.insert("Lon"); 00581 d_lon_names.insert("LON"); 00582 } 00583 00594 void GeoConstraint::set_bounding_box(double top, double left, 00595 double bottom, double right) 00596 { 00597 // Ensure this method is called only once. What about pthreads? 00598 // The method Array::reset_constraint() might make this so it could be 00599 // called more than once. jhrg 8/30/06 00600 if (d_bounding_box_set) 00601 throw Error("It is not possible to register more than one geographical constraint on a variable."); 00602 00603 // Record the 'sense' of the latitude for use here and later on. 00604 d_latitude_sense = categorize_latitude(); 00605 #if 0 00606 if (d_latitude_sense == inverted) 00607 throw Error("geogrid() does not currently work with inverted data (data where the north pole is at a negative latitude value)."); 00608 #endif 00609 00610 // Categorize the notation used by the bounding box (0/359 or -180/179). 00611 d_longitude_notation = categorize_notation(left, right); 00612 00613 // If the notation uses -180/179, transform the request to 0/359 notation. 00614 if (d_longitude_notation == neg_pos) 00615 transform_constraint_to_pos_notation(left, right); 00616 00617 // If the grid uses -180/179, transform it to 0/359 as well. This will make 00618 // subsequent logic easier and adds only a few extra operations, even with 00619 // large maps. 00620 Notation longitude_notation = 00621 categorize_notation(d_lon[0], d_lon[d_lon_length - 1]); 00622 00623 if (longitude_notation == neg_pos) 00624 transform_longitude_to_pos_notation(); 00625 00626 if (!is_bounding_box_valid(left, top, right, bottom)) 00627 throw Error("The bounding box does not intersect any data within this Grid or Array. The\ngeographical extent of these data are from latitude " 00628 + double_to_string(d_lat[0]) + " to " 00629 + double_to_string(d_lat[d_lat_length-1]) 00630 + "\nand longitude " + double_to_string(d_lon[0]) 00631 + " to " + double_to_string(d_lon[d_lon_length-1]) 00632 + " while the bounding box provided was latitude " 00633 + double_to_string(top) + " to " 00634 + double_to_string(bottom) + "\nand longitude " 00635 + double_to_string(left) + " to " 00636 + double_to_string(right)); 00637 00638 // This is simpler than the longitude case because there's no need to 00639 // test for several notations, no need to accommodate them in the return, 00640 // no modulo arithmetic for the axis and no need to account for a 00641 // constraint with two disconnected parts to be joined. 00642 find_latitude_indeces(top, bottom, d_latitude_sense, 00643 d_latitude_index_top, d_latitude_index_bottom); 00644 00645 00646 // Find the longitude map indexes that correspond to the bounding box. 00647 find_longitude_indeces(left, right, d_longitude_index_left, 00648 d_longitude_index_right); 00649 00650 DBG(cerr << "Bounding box (tlbr): " << d_latitude_index_top << ", " 00651 << d_longitude_index_left << ", " 00652 << d_latitude_index_bottom << ", " 00653 << d_longitude_index_right << endl); 00654 00655 d_bounding_box_set = true; 00656 } 00657 00658 } // namespace libdap