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 24370 2011-03-28 16:21: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 vector<char> tmp_data(d_array_data_size); 00407 int array_elem_size = a.var()->width(); 00408 int lat_lon_size = (d_array_data_size / size); 00409 00410 DBG(cerr << "lat, lon_length: " << lat_length << ", " << lon_length << endl); 00411 DBG(cerr << "size: " << size << endl); 00412 DBG(cerr << "d_array_data_size: " << d_array_data_size << endl); 00413 DBG(cerr << "array_elem_size: " << array_elem_size<< endl); 00414 DBG(cerr << "lat_lon_size: " << lat_lon_size<< endl); 00415 00416 for (int i = 0; i < size; ++i) { 00417 int lat = 0; 00418 int s_lat = lat_length - 1; 00419 // lon_length is the element size; memcpy() needs the number of bytes 00420 int lon_size = array_elem_size * lon_length; 00421 int offset = i * lat_lon_size; 00422 while (s_lat > -1) 00423 memcpy(&tmp_data[0] + offset + (lat++ * lon_size), 00424 d_array_data + offset + (s_lat-- * lon_size), 00425 lon_size); 00426 } 00427 00428 memcpy(d_array_data, &tmp_data[0], d_array_data_size); 00429 // delete [] tmp_data; 00430 } 00431 00440 void GeoConstraint::reorder_longitude_map(int longitude_index_left) 00441 { 00442 double *tmp_lon = new double[d_lon_length]; 00443 00444 swap_vector_ends((char *) tmp_lon, (char *) d_lon, d_lon_length, 00445 longitude_index_left, sizeof(double)); 00446 00447 memcpy(d_lon, tmp_lon, d_lon_length * sizeof(double)); 00448 00449 delete[]tmp_lon; 00450 } 00451 00452 static int 00453 count_dimensions_except_longitude(Array &a) 00454 { 00455 int size = 1; 00456 for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i) 00457 size *= a.dimension_size(i, true); 00458 00459 return size; 00460 } 00461 00479 void GeoConstraint::reorder_data_longitude_axis(Array &a, Array::Dim_iter lon_dim) 00480 { 00481 00482 if (!is_longitude_rightmost()) 00483 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."); 00484 00485 DBG(cerr << "Constraint for the left half: " << get_longitude_index_left() 00486 << ", " << get_lon_length() - 1 << endl); 00487 00488 // Build a constraint for the left part and get those values 00489 a.add_constraint(lon_dim, get_longitude_index_left(), 1, 00490 get_lon_length() - 1); 00491 a.set_read_p(false); 00492 a.read(); 00493 DBG2(a.print_val(stderr)); 00494 00495 // Save the left-hand data to local storage 00496 int left_size = a.width(); // width() == length() * element size 00497 char *left_data = (char*)a.value(); // value() allocates and copies 00498 00499 // Build a constraint for the 'right' part, which goes from the left edge 00500 // of the array to the right index and read those data. 00501 // (Don't call a.clear_constraint() since that will clear the constraint on 00502 // all the dimensions while add_constraint() will replace a constraint on 00503 // the given dimension). 00504 00505 DBG(cerr << "Constraint for the right half: " << 0 00506 << ", " << get_longitude_index_right() << endl); 00507 00508 a.add_constraint(lon_dim, 0, 1, get_longitude_index_right()); 00509 a.set_read_p(false); 00510 a.read(); 00511 DBG2(a.print_val(stderr)); 00512 00513 char *right_data = (char*)a.value(); 00514 int right_size = a.width(); 00515 00516 // Make one big lump O'data 00517 d_array_data_size = left_size + right_size; 00518 d_array_data = new char[d_array_data_size]; 00519 00520 // Assume COARDS conventions are being followed: lon varies fastest. 00521 // These *_elements variables are actually elements * bytes/element since 00522 // memcpy() uses bytes. 00523 int elem_size = a.var()->width(); 00524 int left_row_size = (get_lon_length() - get_longitude_index_left()) * elem_size; 00525 int right_row_size = (get_longitude_index_right() + 1) * elem_size; 00526 int total_bytes_per_row = left_row_size + right_row_size; 00527 00528 DBG2(cerr << "elem_size: " << elem_size << "; left & right size: " 00529 << left_row_size << ", " << right_row_size << endl); 00530 00531 // This will work for any number of dimension so long as longitude is the 00532 // right-most array dimension. 00533 int rows_to_copy = count_dimensions_except_longitude(a); 00534 for (int i = 0; i < rows_to_copy; ++i) { 00535 DBG(cerr << "Copying " << i << "th row" << endl); 00536 DBG(cerr << "left memcpy: " << *(float *)(left_data + (left_row_size * i)) << endl); 00537 00538 memcpy(d_array_data + (total_bytes_per_row * i), 00539 left_data + (left_row_size * i), 00540 left_row_size); 00541 00542 DBG(cerr << "right memcpy: " << *(float *)(right_data + (right_row_size * i)) << endl); 00543 00544 memcpy(d_array_data + (total_bytes_per_row * i) + left_row_size, 00545 right_data + (right_row_size * i), 00546 right_row_size); 00547 } 00548 00549 delete[]left_data; 00550 delete[]right_data; 00551 } 00552 00555 GeoConstraint::GeoConstraint() 00556 : d_array_data(0), d_array_data_size(0), 00557 d_lat(0), d_lon(0), 00558 d_bounding_box_set(false), 00559 d_longitude_rightmost(false), 00560 d_longitude_notation(unknown_notation), 00561 d_latitude_sense(unknown_sense) 00562 { 00563 // Build sets of attribute values for easy searching. Maybe overkill??? 00564 d_coards_lat_units.insert("degrees_north"); 00565 d_coards_lat_units.insert("degree_north"); 00566 d_coards_lat_units.insert("degree_N"); 00567 d_coards_lat_units.insert("degrees_N"); 00568 00569 d_coards_lon_units.insert("degrees_east"); 00570 d_coards_lon_units.insert("degree_east"); 00571 d_coards_lon_units.insert("degrees_E"); 00572 d_coards_lon_units.insert("degree_E"); 00573 00574 d_lat_names.insert("COADSY"); 00575 d_lat_names.insert("lat"); 00576 d_lat_names.insert("Lat"); 00577 d_lat_names.insert("LAT"); 00578 00579 d_lon_names.insert("COADSX"); 00580 d_lon_names.insert("lon"); 00581 d_lon_names.insert("Lon"); 00582 d_lon_names.insert("LON"); 00583 } 00584 00595 void GeoConstraint::set_bounding_box(double top, double left, 00596 double bottom, double right) 00597 { 00598 // Ensure this method is called only once. What about pthreads? 00599 // The method Array::reset_constraint() might make this so it could be 00600 // called more than once. jhrg 8/30/06 00601 if (d_bounding_box_set) 00602 throw Error("It is not possible to register more than one geographical constraint on a variable."); 00603 00604 // Record the 'sense' of the latitude for use here and later on. 00605 d_latitude_sense = categorize_latitude(); 00606 #if 0 00607 if (d_latitude_sense == inverted) 00608 throw Error("geogrid() does not currently work with inverted data (data where the north pole is at a negative latitude value)."); 00609 #endif 00610 00611 // Categorize the notation used by the bounding box (0/359 or -180/179). 00612 d_longitude_notation = categorize_notation(left, right); 00613 00614 // If the notation uses -180/179, transform the request to 0/359 notation. 00615 if (d_longitude_notation == neg_pos) 00616 transform_constraint_to_pos_notation(left, right); 00617 00618 // If the grid uses -180/179, transform it to 0/359 as well. This will make 00619 // subsequent logic easier and adds only a few extra operations, even with 00620 // large maps. 00621 Notation longitude_notation = 00622 categorize_notation(d_lon[0], d_lon[d_lon_length - 1]); 00623 00624 if (longitude_notation == neg_pos) 00625 transform_longitude_to_pos_notation(); 00626 00627 if (!is_bounding_box_valid(left, top, right, bottom)) 00628 throw Error("The bounding box does not intersect any data within this Grid or Array. The\ngeographical extent of these data are from latitude " 00629 + double_to_string(d_lat[0]) + " to " 00630 + double_to_string(d_lat[d_lat_length-1]) 00631 + "\nand longitude " + double_to_string(d_lon[0]) 00632 + " to " + double_to_string(d_lon[d_lon_length-1]) 00633 + " while the bounding box provided was latitude " 00634 + double_to_string(top) + " to " 00635 + double_to_string(bottom) + "\nand longitude " 00636 + double_to_string(left) + " to " 00637 + double_to_string(right)); 00638 00639 // This is simpler than the longitude case because there's no need to 00640 // test for several notations, no need to accommodate them in the return, 00641 // no modulo arithmetic for the axis and no need to account for a 00642 // constraint with two disconnected parts to be joined. 00643 find_latitude_indeces(top, bottom, d_latitude_sense, 00644 d_latitude_index_top, d_latitude_index_bottom); 00645 00646 00647 // Find the longitude map indexes that correspond to the bounding box. 00648 find_longitude_indeces(left, right, d_longitude_index_left, 00649 d_longitude_index_right); 00650 00651 DBG(cerr << "Bounding box (tlbr): " << d_latitude_index_top << ", " 00652 << d_longitude_index_left << ", " 00653 << d_latitude_index_bottom << ", " 00654 << d_longitude_index_right << endl); 00655 00656 d_bounding_box_set = true; 00657 } 00658 00659 } // namespace libdap