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