libdap++ Updated for version 3.8.2

GeoConstraint.cc

Go to the documentation of this file.
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