libdap++  Updated for version 3.8.2
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 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 static char id[] not_used =
32  { "$Id: GridGeoConstraint.cc 24281 2011-03-09 00:22:31Z jimg $"
33  };
34 
35 #include <cmath>
36 
37 #include <iostream>
38 #include <sstream>
39 
40 //#define DODS_DEBUG
41 
42 #include "debug.h"
43 #include "dods-datatypes.h"
44 #include "GridGeoConstraint.h"
45 #include "Float64.h"
46 
47 #include "Error.h"
48 #include "InternalErr.h"
49 #include "ce_functions.h"
50 #include "util.h"
51 
52 using namespace std;
53 
54 namespace libdap {
55 
62 GridGeoConstraint::GridGeoConstraint(Grid *grid)
63  : GeoConstraint(), d_grid(grid), d_latitude(0), d_longitude(0)
64 {
65  if (d_grid->get_array()->dimensions() < 2
66  || d_grid->get_array()->dimensions() > 3)
67  throw Error("The geogrid() function works only with Grids of two or three dimensions.");
68 
69  // Is this Grid a geo-referenced grid? Throw Error if not.
70  if (!build_lat_lon_maps())
71  throw Error(string("The grid '") + d_grid->name()
72  + "' does not have identifiable latitude/longitude map vectors.");
73 
74  if (!lat_lon_dimensions_ok())
75  throw Error("The geogrid() function will only work when the Grid's Longitude and Latitude\nmaps are the rightmost dimensions.");
76 }
77 
79  : GeoConstraint(), d_grid(grid), d_latitude(0), d_longitude(0)
80 {
81  if (d_grid->get_array()->dimensions() < 2
82  || d_grid->get_array()->dimensions() > 3)
83  throw Error("The geogrid() function works only with Grids of two or three dimensions.");
84 
85  // Is this Grid a geo-referenced grid? Throw Error if not.
86  if (!build_lat_lon_maps(lat, lon))
87  throw Error(string("The grid '") + d_grid->name()
88  + "' does not have valid latitude/longitude map vectors.");
89 
90 
91  if (!lat_lon_dimensions_ok())
92  throw Error("The geogrid() function will only work when the Grid's Longitude and Latitude\nmaps are the rightmost dimensions.");
93 }
94 
110 bool GridGeoConstraint::build_lat_lon_maps()
111 {
112  Grid::Map_iter m = d_grid->map_begin();
113  // Assume that a Grid is correct and thus has exactly as many maps as its
114  // array part has dimensions. Thus don't bother to test the Grid's array
115  // dimension iterator for '!= dim_end()'.
116  Array::Dim_iter d = d_grid->get_array()->dim_begin();
117  // The fields d_latitude and d_longitude may be initialized to null or they
118  // may already contain pointers to the maps to use. In the latter case,
119  // skip the heuristics used in this code. However, given that all this
120  // method does is find the lat/lon maps, if they are given in the ctor,
121  // This method will likely not be called at all.
122  while (m != d_grid->map_end() && (!d_latitude || !d_longitude)) {
123  string units_value = (*m)->get_attr_table().get_attr("units");
124  units_value = remove_quotes(units_value);
125  string map_name = (*m)->name();
126 
127  // The 'units' attribute must match exactly; the name only needs to
128  // match a prefix.
129  if (!d_latitude
131  units_value, map_name)) {
132 
133  // Set both d_latitude (a pointer to the real map vector) and
134  // d_lat, a vector of the values represented as doubles. It's easier
135  // to work with d_lat, but it's d_latitude that needs to be set
136  // when constraining the grid. Also, record the grid variable's
137  // dimension iterator so that it's easier to set the Grid's Array
138  // (which also has to be constrained).
139  d_latitude = dynamic_cast < Array * >(*m);
140  if (!d_latitude)
141  throw InternalErr(__FILE__, __LINE__, "Expected an array.");
142  if (!d_latitude->read_p())
143  d_latitude->read();
144 
145  set_lat(extract_double_array(d_latitude)); // throws Error
146  set_lat_length(d_latitude->length());
147 
148  set_lat_dim(d);
149  }
150 
151  if (!d_longitude // && !units_value.empty()
153  units_value, map_name)) {
154 
155  d_longitude = dynamic_cast < Array * >(*m);
156  if (!d_longitude)
157  throw InternalErr(__FILE__, __LINE__, "Expected an array.");
158  if (!d_longitude->read_p())
159  d_longitude->read();
160 
161  set_lon(extract_double_array(d_longitude));
162  set_lon_length(d_longitude->length());
163 
164  set_lon_dim(d);
165 
166  if (m + 1 == d_grid->map_end())
168  }
169 
170  ++m;
171  ++d;
172  }
173 
174  return get_lat() && get_lon();
175 }
176 
184 bool GridGeoConstraint::build_lat_lon_maps(Array *lat, Array *lon)
185 {
186  Grid::Map_iter m = d_grid->map_begin();
187 
188  Array::Dim_iter d = d_grid->get_array()->dim_begin();
189 
190  while (m != d_grid->map_end() && (!d_latitude || !d_longitude)) {
191  // Look for the Grid map that matches the variable passed as 'lat'
192  if (!d_latitude && *m == lat) {
193 
194  d_latitude = lat;
195 
196  if (!d_latitude->read_p())
197  d_latitude->read();
198 
199  set_lat(extract_double_array(d_latitude)); // throws Error
200  set_lat_length(d_latitude->length());
201 
202  set_lat_dim(d);
203  }
204 
205  if (!d_longitude && *m == lon) {
206 
207  d_longitude = lon;
208 
209  if (!d_longitude->read_p())
210  d_longitude->read();
211 
212  set_lon(extract_double_array(d_longitude));
213  set_lon_length(d_longitude->length());
214 
215  set_lon_dim(d);
216 
217  if (m + 1 == d_grid->map_end())
219  }
220 
221  ++m;
222  ++d;
223  }
224 
225  return get_lat() && get_lon();
226 }
227 
238 bool
239 GridGeoConstraint::lat_lon_dimensions_ok()
240 {
241  // get the last two map iterators
242  Grid::Map_riter rightmost = d_grid->map_rbegin();
243  Grid::Map_riter next_rightmost = rightmost + 1;
244 
245  if (*rightmost == d_longitude && *next_rightmost == d_latitude)
247  else if (*rightmost == d_latitude && *next_rightmost == d_longitude)
249  else
250  return false;
251 
252  return true;
253 }
254 
277 {
278  if (!is_bounding_box_set())
279  throw InternalErr("The Latitude and Longitude constraints must be set before calling apply_constraint_to_data().");
280 
281  Array::Dim_iter fd = d_latitude->dim_begin();
282 
283  if (get_latitude_sense() == inverted) {
284  int tmp = get_latitude_index_top();
287  }
288 
289  // It's easy to flip the Latitude values; if the bottom index value
290  // is before/above the top index, return an error explaining that.
292  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.");
293 
294  // Constrain the lat vector and lat dim of the array
295  d_latitude->add_constraint(fd, get_latitude_index_top(), 1,
297  d_grid->get_array()->add_constraint(get_lat_dim(),
300 
301  // Does the longitude constraint cross the edge of the longitude vector?
302  // If so, reorder the grid's data (array), longitude map vector and the
303  // local vector of longitude data used for computation.
306 
307  // If the longitude constraint is 'split', join the two parts, reload
308  // the data into the Grid's Array and make sure the Array is marked as
309  // already read. This should be true for the whole Grid, but if some
310  // future modification changes that, the array will be covered here.
311  // Note that the following method only reads the data out and stores
312  // it in this object after joining the two parts. The method
313  // apply_constraint_to_data() transfers the data back from the this
314  // object to the DAP Grid variable.
316 
317  // Now that the data are all in local storage alter the indices; the
318  // left index has now been moved to 0, and the right index is now
319  // at lon_vector_length-left+right.
323  }
324 
325  // If the constraint used the -180/179 (neg_pos) notation, transform
326  // the longitude map so it uses the -180/179 notation. Note that at this
327  // point, d_longitude always uses the pos notation because of the earlier
328  // conditional transformation.
329 
330  // Do this _before_ applying the constraint since set_array_using_double()
331  // tests the array length using Vector::length() and that method returns
332  // the length _as constrained_. We want to move all of the longitude
333  // values from d_lon back into the map, not just the number that will be
334  // sent (although an optimization might do this, it's hard to imagine
335  // it would gain much).
336  if (get_longitude_notation() == neg_pos) {
338  }
339 
340  // Apply constraint; stride is always one and maps only have one dimension
341  fd = d_longitude->dim_begin();
342  d_longitude->add_constraint(fd, get_longitude_index_left(), 1,
344 
345  d_grid->get_array()->add_constraint(get_lon_dim(),
348 
349  // Transfer values from the local lat vector to the Grid's
350  // Here test the sense of the latitude vector and invert the vector if the
351  // sense is 'inverted' so that the top is always the northern-most value
352  if (get_latitude_sense() == inverted) {
353  DBG(cerr << "Inverted latitude sense" << endl);
356  // Now read the Array data and flip the latitudes.
360  }
361 
364 
367 
368  // Look for any non-lat/lon maps and make sure they are read correctly
369  Grid::Map_iter i = d_grid->map_begin();
370  Grid::Map_iter end = d_grid->map_end();
371  while (i != end) {
372  if (*i != d_latitude && *i != d_longitude) {
373  if ((*i)->send_p()) {
374  DBG(cerr << "reading grid map: " << (*i)->name() << endl);
375  //(*i)->set_read_p(false);
376  (*i)->read();
377  }
378  }
379  ++i;
380  }
381 
382  // ... and then the Grid's array if it has been read.
383  if (get_array_data()) {
384  int size = d_grid->get_array()->val2buf(get_array_data());
385 
386  if (size != get_array_data_size())
387  throw InternalErr(__FILE__, __LINE__, "Expected data size not copied to the Grid's buffer.");
388 
389  d_grid->set_read_p(true);
390  }
391  else {
392  d_grid->get_array()->read();
393  }
394 }
395 
396 } // namespace libdap