• Main Page
  • Related Pages
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

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 21922 2010-01-07 18:23:57Z 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,
00209         int &longitude_index_right) const
00210 {
00211     // Use all values 'modulo 360' to take into account the cases where the
00212     // constraint and/or the longitude vector are given using values greater
00213     // than 360 (i.e., when 380 is used to mean 20).
00214     double t_left = fmod(left, 360.0);
00215     double t_right = fmod(right, 360.0);
00216 
00217     // Find the place where 'longitude starts.' That is, what value of the
00218     // index 'i' corresponds to the smallest value of d_lon. Why we do this:
00219     // Some data sources use offset longitude axes so that the 'seam' is
00220     // shifted to a place other than the date line.
00221     int i = 0;
00222     int lon_origin_index = 0;
00223     double smallest_lon = fmod(d_lon[0], 360.0);
00224     while (i < d_lon_length) {
00225         double curent_lon_value = fmod(d_lon[i], 360.0);
00226         if (smallest_lon > curent_lon_value) {
00227             smallest_lon = curent_lon_value;
00228             lon_origin_index = i;
00229         }
00230         ++i;
00231     }
00232 
00233     DBG2(cerr << "lon_origin_index: " << lon_origin_index << endl);
00234 
00235     // Scan from the index of the smallest value looking for the place where
00236     // the value is greater than or equal to the left most point of the bounding
00237     // box.
00238     i = lon_origin_index;
00239     while (fmod(d_lon[i], 360.0) < t_left) {
00240         ++i;
00241         i = i % d_lon_length;
00242 
00243         // If we cycle completely through all the values/indices, throw
00244         if (i == lon_origin_index)
00245             throw Error("geogrid: Could not find an index for the longitude value '" + long_to_string(left) + "'");
00246     }
00247 
00248     if (fmod(d_lon[i], 360.0) == t_left)
00249         longitude_index_left = i;
00250     else
00251         longitude_index_left = (i - 1) > 0 ? i - 1 : 0;
00252 
00253     DBG2(cerr << "longitude_index_left: " << longitude_index_left << endl);
00254 
00255     // Assume the vector is circular --> the largest value is next to the
00256     // smallest.
00257     int largest_lon_index = (lon_origin_index - 1 + d_lon_length) % d_lon_length;
00258     i = largest_lon_index;
00259     while (fmod(d_lon[i], 360.0) > t_right) {
00260         // This is like modulus but for 'counting down'
00261         i = (i == 0) ? d_lon_length - 1 : i - 1;
00262         if (i == largest_lon_index)
00263             throw Error("geogrid: Could not find an index for the longitude value '" + long_to_string(right) + "'");
00264     }
00265 
00266     if (fmod(d_lon[i], 360.0) == t_right)
00267         longitude_index_right = i;
00268     else
00269         longitude_index_right = (i + 1) < d_lon_length - 1 ? i + 1 : d_lon_length - 1;
00270 
00271     DBG2(cerr << "longitude_index_right: " << longitude_index_right << endl);
00272 }
00273 
00286 void GeoConstraint::find_latitude_indeces(double top, double bottom,
00287         LatitudeSense sense,
00288         int &latitude_index_top,
00289         int &latitude_index_bottom) const
00290 {
00291     int i, j;
00292 
00293     if (sense == normal) {
00294         i = 0;
00295         while (i < d_lat_length - 1 && top < d_lat[i])
00296             ++i;
00297 
00298         j = d_lat_length - 1;
00299         while (j > 0 && bottom > d_lat[j])
00300             --j;
00301 
00302         if (d_lat[i] == top)
00303             latitude_index_top = i;
00304         else
00305             latitude_index_top = (i - 1) > 0 ? i - 1 : 0;
00306 
00307         if (d_lat[j] == bottom)
00308             latitude_index_bottom = j;
00309         else
00310             latitude_index_bottom =
00311                 (j + 1) < d_lat_length - 1 ? j + 1 : d_lat_length - 1;
00312     }
00313     else {
00314         i = d_lat_length - 1;
00315         while (i > 0 && d_lat[i] > top)
00316             --i;
00317 
00318         j = 0;
00319         while (j < d_lat_length - 1 && d_lat[j] < bottom)
00320             ++j;
00321 
00322         if (d_lat[i] == top)
00323             latitude_index_top = i;
00324         else
00325             latitude_index_top = (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1;
00326 
00327         if (d_lat[j] == bottom)
00328             latitude_index_bottom = j;
00329         else
00330             latitude_index_bottom = (j - 1) > 0 ? j - 1 : 0;
00331     }
00332 }
00333 
00337 GeoConstraint::LatitudeSense GeoConstraint::categorize_latitude() const
00338 {
00339     return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted;
00340 }
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 
00352 template<class T>
00353 static void transpose(std::vector<std::vector<T> > a,
00354         std::vector<std::vector<T> > b, int width, int height)
00355 {
00356     for (int i = 0; i < width; i++) {
00357         for (int j = 0; j < height; j++) {
00358             b[j][i] = a[i][j];
00359         }
00360     }
00361 }
00362 
00370 void GeoConstraint::transpose_vector(double *src, const int length)
00371 {
00372     double *tmp = new double[length];
00373 
00374     int i = 0, j = length-1;
00375     while (i < length)
00376         tmp[j--] = src[i++];
00377 
00378     memcpy(src, tmp,length * sizeof(double));
00379 
00380     delete[] tmp;
00381 }
00382 
00383 static int
00384 count_size_except_latitude_and_longitude(Array &a)
00385 {
00386     if (a.dim_end() - a.dim_begin() <= 2)       // < 2 is really an error...
00387         return 1;
00388 
00389     int size = 1;
00390     for (Array::Dim_iter i = a.dim_begin(); (i + 2) != a.dim_end(); ++i)
00391         size *= a.dimension_size(i, true);
00392 
00393     return size;
00394 }
00395 
00396 void GeoConstraint::flip_latitude_within_array(Array &a, int lat_length,
00397         int lon_length)
00398 {
00399     if (!d_array_data) {
00400         a.read();
00401         d_array_data = static_cast<char*>(a.value());
00402         d_array_data_size = a.width();  // Bytes not elements
00403     }
00404 
00405     int size = count_size_except_latitude_and_longitude(a);
00406     char *tmp_data = new char[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 + 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, 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 indeces 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

Generated on Mon Sep 6 2010 18:58:17 for libdap++ by  doxygen 1.7.1