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 25066 2011-11-30 18:39:37Z 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 
114  // Assume that a Grid is correct and thus has exactly as many maps as its
115  // array part has dimensions. Thus don't bother to test the Grid's array
116  // dimension iterator for '!= dim_end()'.
117  Array::Dim_iter d = d_grid->get_array()->dim_begin();
118 
119  // The fields d_latitude and d_longitude may be initialized to null or they
120  // may already contain pointers to the maps to use. In the latter case,
121  // skip the heuristics used in this code. However, given that all this
122  // method does is find the lat/lon maps, if they are given in the ctor,
123  // This method will likely not be called at all.
124  while (m != d_grid->map_end() && (!d_latitude || !d_longitude)) {
125  string units_value = (*m)->get_attr_table().get_attr("units");
126  units_value = remove_quotes(units_value);
127  string map_name = (*m)->name();
128 
129  // The 'units' attribute must match exactly; the name only needs to
130  // match a prefix.
131  if (!d_latitude
133  units_value, map_name)) {
134 
135  // Set both d_latitude (a pointer to the real map vector) and
136  // d_lat, a vector of the values represented as doubles. It's easier
137  // to work with d_lat, but it's d_latitude that needs to be set
138  // when constraining the grid. Also, record the grid variable's
139  // dimension iterator so that it's easier to set the Grid's Array
140  // (which also has to be constrained).
141  d_latitude = dynamic_cast < Array * >(*m);
142  if (!d_latitude)
143  throw InternalErr(__FILE__, __LINE__, "Expected an array.");
144  if (!d_latitude->read_p())
145  d_latitude->read();
146 
147  set_lat(extract_double_array(d_latitude)); // throws Error
148  set_lat_length(d_latitude->length());
149 
150  set_lat_dim(d);
151  }
152 
153  if (!d_longitude // && !units_value.empty()
155  units_value, map_name)) {
156 
157  d_longitude = dynamic_cast < Array * >(*m);
158  if (!d_longitude)
159  throw InternalErr(__FILE__, __LINE__, "Expected an array.");
160  if (!d_longitude->read_p())
161  d_longitude->read();
162 
163  set_lon(extract_double_array(d_longitude));
164  set_lon_length(d_longitude->length());
165 
166  set_lon_dim(d);
167 
168  if (m + 1 == d_grid->map_end())
170  }
171 
172  ++m;
173  ++d;
174  }
175 
176  return get_lat() && get_lon();
177 }
178 
186 bool GridGeoConstraint::build_lat_lon_maps(Array *lat, Array *lon)
187 {
188  Grid::Map_iter m = d_grid->map_begin();
189 
190  Array::Dim_iter d = d_grid->get_array()->dim_begin();
191 
192  while (m != d_grid->map_end() && (!d_latitude || !d_longitude)) {
193  // Look for the Grid map that matches the variable passed as 'lat'
194  if (!d_latitude && *m == lat) {
195 
196  d_latitude = lat;
197 
198  if (!d_latitude->read_p())
199  d_latitude->read();
200 
201  set_lat(extract_double_array(d_latitude)); // throws Error
202  set_lat_length(d_latitude->length());
203 
204  set_lat_dim(d);
205  }
206 
207  if (!d_longitude && *m == lon) {
208 
209  d_longitude = lon;
210 
211  if (!d_longitude->read_p())
212  d_longitude->read();
213 
214  set_lon(extract_double_array(d_longitude));
215  set_lon_length(d_longitude->length());
216 
217  set_lon_dim(d);
218 
219  if (m + 1 == d_grid->map_end())
221  }
222 
223  ++m;
224  ++d;
225  }
226 
227  return get_lat() && get_lon();
228 }
229 
240 bool
241 GridGeoConstraint::lat_lon_dimensions_ok()
242 {
243  // get the last two map iterators
244  Grid::Map_riter rightmost = d_grid->map_rbegin();
245  Grid::Map_riter next_rightmost = rightmost + 1;
246 
247  if (*rightmost == d_longitude && *next_rightmost == d_latitude)
249  else if (*rightmost == d_latitude && *next_rightmost == d_longitude)
251  else
252  return false;
253 
254  return true;
255 }
256 
279 {
280  if (!is_bounding_box_set())
281  throw InternalErr("The Latitude and Longitude constraints must be set before calling apply_constraint_to_data().");
282 
283  Array::Dim_iter fd = d_latitude->dim_begin();
284 
285  if (get_latitude_sense() == inverted) {
286  int tmp = get_latitude_index_top();
289  }
290 
291  // It's easy to flip the Latitude values; if the bottom index value
292  // is before/above the top index, return an error explaining that.
294  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.");
295 
296  // Constrain the lat vector and lat dim of the array
297  d_latitude->add_constraint(fd, get_latitude_index_top(), 1,
299  d_grid->get_array()->add_constraint(get_lat_dim(),
302 
303  // Does the longitude constraint cross the edge of the longitude vector?
304  // If so, reorder the grid's data (array), longitude map vector and the
305  // local vector of longitude data used for computation.
308 
309  // If the longitude constraint is 'split', join the two parts, reload
310  // the data into the Grid's Array and make sure the Array is marked as
311  // already read. This should be true for the whole Grid, but if some
312  // future modification changes that, the array will be covered here.
313  // Note that the following method only reads the data out and stores
314  // it in this object after joining the two parts. The method
315  // apply_constraint_to_data() transfers the data back from the this
316  // object to the DAP Grid variable.
318 
319  // Now that the data are all in local storage alter the indices; the
320  // left index has now been moved to 0, and the right index is now
321  // at lon_vector_length-left+right.
325  }
326 
327  // If the constraint used the -180/179 (neg_pos) notation, transform
328  // the longitude map so it uses the -180/179 notation. Note that at this
329  // point, d_longitude always uses the pos notation because of the earlier
330  // conditional transformation.
331 
332  // Do this _before_ applying the constraint since set_array_using_double()
333  // tests the array length using Vector::length() and that method returns
334  // the length _as constrained_. We want to move all of the longitude
335  // values from d_lon back into the map, not just the number that will be
336  // sent (although an optimization might do this, it's hard to imagine
337  // it would gain much).
338  if (get_longitude_notation() == neg_pos) {
340  }
341 
342  // Apply constraint; stride is always one and maps only have one dimension
343  fd = d_longitude->dim_begin();
344  d_longitude->add_constraint(fd, get_longitude_index_left(), 1,
346 
347  d_grid->get_array()->add_constraint(get_lon_dim(),
350 
351  // Transfer values from the local lat vector to the Grid's
352  // Here test the sense of the latitude vector and invert the vector if the
353  // sense is 'inverted' so that the top is always the northern-most value
354  if (get_latitude_sense() == inverted) {
355  DBG(cerr << "Inverted latitude sense" << endl);
358  // Now read the Array data and flip the latitudes.
362  }
363 
366 
369 
370  // Look for any non-lat/lon maps and make sure they are read correctly
371  Grid::Map_iter i = d_grid->map_begin();
372  Grid::Map_iter end = d_grid->map_end();
373  while (i != end) {
374  if (*i != d_latitude && *i != d_longitude) {
375  if ((*i)->send_p()) {
376  DBG(cerr << "reading grid map: " << (*i)->name() << endl);
377  //(*i)->set_read_p(false);
378  (*i)->read();
379  }
380  }
381  ++i;
382  }
383 
384  // ... and then the Grid's array if it has been read.
385  if (get_array_data()) {
386  int size = d_grid->get_array()->val2buf(get_array_data());
387 
388  if (size != get_array_data_size())
389  throw InternalErr(__FILE__, __LINE__, "Expected data size not copied to the Grid's buffer.");
390 
391  d_grid->set_read_p(true);
392  }
393  else {
394  d_grid->get_array()->read();
395  }
396 }
397 
398 } // namespace libdap