OPeNDAP Hyrax Back End Server (BES)  Updated for version 3.8.3
GridGeoConstraint.cc
Go to the documentation of this file.
1 
2 // -*- mode: c++; c-basic-offset:4 -*-
3 
4 // This file is part of libdap, A C++ implementation of the OPeNDAP Data
5 // Access Protocol.
6 
7 // Copyright (c) 2002,2003 OPeNDAP, Inc.
8 // Author: James Gallagher <jgallagher@opendap.org>
9 //
10 // This library is free software; you can redistribute it and/or
11 // modify it under the terms of the GNU Lesser General Public
12 // License as published by the Free Software Foundation; either
13 // version 2.1 of the License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 //
24 // You can contact OPeNDAP, Inc. at PO Box 112, Saunderstown, RI. 02874-0112.
25 
26 // The Grid Selection Expression Clause class.
27 
28 
29 #include "config.h"
30 
31 #include <cmath>
32 
33 #include <iostream>
34 #include <sstream>
35 
36 //#define DODS_DEBUG
37 
38 #include <Float64.h>
39 #include <dods-datatypes.h>
40 #include <Error.h>
41 #include <InternalErr.h>
42 #include <util.h>
43 #include <debug.h>
44 
45 #include "GridGeoConstraint.h"
46 
47 using namespace std;
48 
49 namespace libdap {
50 
57 GridGeoConstraint::GridGeoConstraint(Grid *grid)
58  : GeoConstraint(), d_grid(grid), d_latitude(0), d_longitude(0)
59 {
60  if (d_grid->get_array()->dimensions() < 2
61  || d_grid->get_array()->dimensions() > 3)
62  throw Error("The geogrid() function works only with Grids of two or three dimensions.");
63 
64  // Is this Grid a geo-referenced grid? Throw Error if not.
65  if (!build_lat_lon_maps())
66  throw Error(string("The grid '") + d_grid->name()
67  + "' does not have identifiable latitude/longitude map vectors.");
68 
69  if (!lat_lon_dimensions_ok())
70  throw Error("The geogrid() function will only work when the Grid's Longitude and Latitude maps are the rightmost dimensions (grid: " + grid->name() + ", 1).");
71 }
72 
73 GridGeoConstraint::GridGeoConstraint(Grid *grid, Array *lat, Array *lon)
74  : GeoConstraint(), d_grid(grid), d_latitude(0), d_longitude(0)
75 {
76  if (d_grid->get_array()->dimensions() < 2
77  || d_grid->get_array()->dimensions() > 3)
78  throw Error("The geogrid() function works only with Grids of two or three dimensions.");
79 
80  // Is this Grid a geo-referenced grid? Throw Error if not.
81  if (!build_lat_lon_maps(lat, lon))
82  throw Error(string("The grid '") + d_grid->name()
83  + "' does not have valid latitude/longitude map vectors.");
84 
85 
86  if (!lat_lon_dimensions_ok())
87  throw Error("The geogrid() function will only work when the Grid's Longitude and Latitude maps are the rightmost dimensions (grid: " + grid->name() + ", 2).");
88 }
89 
105 bool GridGeoConstraint::build_lat_lon_maps()
106 {
107  Grid::Map_iter m = d_grid->map_begin();
108 
109  // Assume that a Grid is correct and thus has exactly as many maps as its
110  // array part has dimensions. Thus don't bother to test the Grid's array
111  // dimension iterator for '!= dim_end()'.
112  Array::Dim_iter d = d_grid->get_array()->dim_begin();
113 
114  // The fields d_latitude and d_longitude may be initialized to null or they
115  // may already contain pointers to the maps to use. In the latter case,
116  // skip the heuristics used in this code. However, given that all this
117  // method does is find the lat/lon maps, if they are given in the ctor,
118  // This method will likely not be called at all.
119  while (m != d_grid->map_end() && (!d_latitude || !d_longitude)) {
120  string units_value = (*m)->get_attr_table().get_attr("units");
121  units_value = remove_quotes(units_value);
122  string map_name = (*m)->name();
123 
124  // The 'units' attribute must match exactly; the name only needs to
125  // match a prefix.
126  if (!d_latitude
128  units_value, map_name)) {
129 
130  // Set both d_latitude (a pointer to the real map vector) and
131  // d_lat, a vector of the values represented as doubles. It's easier
132  // to work with d_lat, but it's d_latitude that needs to be set
133  // when constraining the grid. Also, record the grid variable's
134  // dimension iterator so that it's easier to set the Grid's Array
135  // (which also has to be constrained).
136  d_latitude = dynamic_cast < Array * >(*m);
137  if (!d_latitude)
138  throw InternalErr(__FILE__, __LINE__, "Expected an array.");
139  if (!d_latitude->read_p())
140  d_latitude->read();
141 
142  set_lat(extract_double_array(d_latitude)); // throws Error
143  set_lat_length(d_latitude->length());
144 
145  set_lat_dim(d);
146  }
147 
148  if (!d_longitude // && !units_value.empty()
150  units_value, map_name)) {
151 
152  d_longitude = dynamic_cast < Array * >(*m);
153  if (!d_longitude)
154  throw InternalErr(__FILE__, __LINE__, "Expected an array.");
155  if (!d_longitude->read_p())
156  d_longitude->read();
157 
158  set_lon(extract_double_array(d_longitude));
159  set_lon_length(d_longitude->length());
160 
161  set_lon_dim(d);
162 
163  if (m + 1 == d_grid->map_end())
165  }
166 
167  ++m;
168  ++d;
169  }
170 
171  return get_lat() && get_lon();
172 }
173 
181 bool GridGeoConstraint::build_lat_lon_maps(Array *lat, Array *lon)
182 {
183  Grid::Map_iter m = d_grid->map_begin();
184 
185  Array::Dim_iter d = d_grid->get_array()->dim_begin();
186 
187  while (m != d_grid->map_end() && (!d_latitude || !d_longitude)) {
188  // Look for the Grid map that matches the variable passed as 'lat'
189  if (!d_latitude && *m == lat) {
190 
191  d_latitude = lat;
192 
193  if (!d_latitude->read_p())
194  d_latitude->read();
195 
196  set_lat(extract_double_array(d_latitude)); // throws Error
197  set_lat_length(d_latitude->length());
198 
199  set_lat_dim(d);
200  }
201 
202  if (!d_longitude && *m == lon) {
203 
204  d_longitude = lon;
205 
206  if (!d_longitude->read_p())
207  d_longitude->read();
208 
209  set_lon(extract_double_array(d_longitude));
210  set_lon_length(d_longitude->length());
211 
212  set_lon_dim(d);
213 
214  if (m + 1 == d_grid->map_end())
216  }
217 
218  ++m;
219  ++d;
220  }
221 
222  return get_lat() && get_lon();
223 }
224 
235 bool
236 GridGeoConstraint::lat_lon_dimensions_ok()
237 {
238  // get the last two map iterators
239  Grid::Map_riter rightmost = d_grid->map_rbegin();
240  Grid::Map_riter next_rightmost = rightmost + 1;
241 
242  if (*rightmost == d_longitude && *next_rightmost == d_latitude)
244  else if (*rightmost == d_latitude && *next_rightmost == d_longitude)
246  else
247  return false;
248 
249  return true;
250 }
251 
274 {
275  if (!is_bounding_box_set())
276  throw InternalErr("The Latitude and Longitude constraints must be set before calling apply_constraint_to_data().");
277 
278  Array::Dim_iter fd = d_latitude->dim_begin();
279 
280  if (get_latitude_sense() == inverted) {
281  int tmp = get_latitude_index_top();
284  }
285 
286  // It's easy to flip the Latitude values; if the bottom index value
287  // is before/above the top index, return an error explaining that.
289  throw Error("The upper and lower latitude indices appear to be reversed. Please provide the latitude bounding box numbers giving the northern-most latitude first.");
290 
291  // Constrain the lat vector and lat dim of the array
292  d_latitude->add_constraint(fd, get_latitude_index_top(), 1,
294  d_grid->get_array()->add_constraint(get_lat_dim(),
297 
298  // Does the longitude constraint cross the edge of the longitude vector?
299  // If so, reorder the grid's data (array), longitude map vector and the
300  // local vector of longitude data used for computation.
303 
304  // If the longitude constraint is 'split', join the two parts, reload
305  // the data into the Grid's Array and make sure the Array is marked as
306  // already read. This should be true for the whole Grid, but if some
307  // future modification changes that, the array will be covered here.
308  // Note that the following method only reads the data out and stores
309  // it in this object after joining the two parts. The method
310  // apply_constraint_to_data() transfers the data back from the this
311  // object to the DAP Grid variable.
312  reorder_data_longitude_axis(*d_grid->get_array(), get_lon_dim());
313 
314  // Now that the data are all in local storage alter the indices; the
315  // left index has now been moved to 0, and the right index is now
316  // at lon_vector_length-left+right.
320  }
321 
322  // If the constraint used the -180/179 (neg_pos) notation, transform
323  // the longitude map so it uses the -180/179 notation. Note that at this
324  // point, d_longitude always uses the pos notation because of the earlier
325  // conditional transformation.
326 
327  // Do this _before_ applying the constraint since set_array_using_double()
328  // tests the array length using Vector::length() and that method returns
329  // the length _as constrained_. We want to move all of the longitude
330  // values from d_lon back into the map, not just the number that will be
331  // sent (although an optimization might do this, it's hard to imagine
332  // it would gain much).
333  if (get_longitude_notation() == neg_pos) {
335  }
336 
337  // Apply constraint; stride is always one and maps only have one dimension
338  fd = d_longitude->dim_begin();
339  d_longitude->add_constraint(fd, get_longitude_index_left(), 1,
341 
342  d_grid->get_array()->add_constraint(get_lon_dim(),
345 
346  // Transfer values from the local lat vector to the Grid's
347  // Here test the sense of the latitude vector and invert the vector if the
348  // sense is 'inverted' so that the top is always the northern-most value
349  if (get_latitude_sense() == inverted) {
350  DBG(cerr << "Inverted latitude sense" << endl);
353  // Now read the Array data and flip the latitudes.
354  flip_latitude_within_array(*d_grid->get_array(),
357  }
358 
359  set_array_using_double(d_latitude, get_lat() + get_latitude_index_top(),
361 
362  set_array_using_double(d_longitude, get_lon() + get_longitude_index_left(),
364 
365  // Look for any non-lat/lon maps and make sure they are read correctly
366  Grid::Map_iter i = d_grid->map_begin();
367  Grid::Map_iter end = d_grid->map_end();
368  while (i != end) {
369  if (*i != d_latitude && *i != d_longitude) {
370  if ((*i)->send_p()) {
371  DBG(cerr << "reading grid map: " << (*i)->name() << endl);
372  //(*i)->set_read_p(false);
373  (*i)->read();
374  }
375  }
376  ++i;
377  }
378 
379  // ... and then the Grid's array if it has been read.
380  if (get_array_data()) {
381  int size = d_grid->get_array()->val2buf(get_array_data());
382 
383  if (size != get_array_data_size())
384  throw InternalErr(__FILE__, __LINE__, "Expected data size not copied to the Grid's buffer.");
385 
386  d_grid->set_read_p(true);
387  }
388  else {
389  d_grid->get_array()->read();
390  }
391 }
392 
393 } // namespace libdap
void set_lat_dim(Array::Dim_iter lat)
char * get_array_data() const
void set_latitude_index_bottom(int bottom)
void set_lon(double *lon)
virtual void reorder_data_longitude_axis(Array &a, Array::Dim_iter lon_dim)
Reorder the data values relative to the longitude axis so that the reordered longitude map (see GeoCo...
void set_latitude_index_top(int top)
set< string > get_coards_lat_units() const
set< string > get_lon_names() const
void set_longitude_rightmost(bool state)
int get_longitude_index_left() const
Array::Dim_iter get_lon_dim() const
void set_lat(double *lat)
bool unit_or_name_match(set< string > units, set< string > names, const string &var_units, const string &var_name)
Look in the containers which hold the units attributes and variable name prefixes which are considere...
void set_lat_length(int len)
virtual void transpose_vector(double *src, const int length)
Given a vector of doubles, transpose the elements.
GridGeoConstraint(Grid *grid)
Initialize GeoConstraint with a Grid.
int get_longitude_index_right() const
int get_array_data_size() const
int get_latitude_index_bottom() const
set< string > get_coards_lon_units() const
virtual void flip_latitude_within_array(Array &a, int lat_length, int lon_length)
virtual void transform_longitude_to_neg_pos_notation()
Given that the Grid has a longitude map that uses the 'pos' notation, transform it to the 'neg_pos' n...
LatitudeSense get_latitude_sense() const
set< string > get_lat_names() const
void set_longitude_index_right(int right)
Notation get_longitude_notation() const
virtual void reorder_longitude_map(int longitude_index_left)
Reorder the elements in the longitude map so that the longitude constraint no longer crosses the edge...
int get_lon_length() const
void set_longitude_index_left(int left)
Encapsulate the logic needed to handle geographical constraints when they are applied to DAP Grid (an...
Definition: GeoConstraint.h:95
double * get_lon() const
Array::Dim_iter get_lat_dim() const
bool is_bounding_box_set() const
void set_lon_dim(Array::Dim_iter lon)
int get_latitude_index_top() const
double * get_lat() const
virtual void apply_constraint_to_data()
Once the bounding box is set use this method to apply the constraint.
void set_lon_length(int len)