40 #include <dods-datatypes.h>
42 #include <InternalErr.h>
46 #include "GridGeoConstraint.h"
59 GridGeoConstraint::GridGeoConstraint(Grid *grid)
60 :
GeoConstraint(), d_grid(grid), d_latitude(0), d_longitude(0)
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.");
67 if (!build_lat_lon_maps())
68 throw Error(
string(
"The grid '") + d_grid->name()
69 +
"' does not have identifiable latitude/longitude map vectors.");
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).");
76 :
GeoConstraint(), d_grid(grid), d_latitude(0), d_longitude(0)
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.");
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.");
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).");
107 bool GridGeoConstraint::build_lat_lon_maps()
109 Grid::Map_iter m = d_grid->map_begin();
114 Array::Dim_iter d = d_grid->get_array()->dim_begin();
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();
129 && unit_or_name_match(get_coards_lat_units(), get_lat_names(),
130 units_value, map_name)) {
138 d_latitude =
dynamic_cast < Array *
>(*m);
140 throw InternalErr(__FILE__, __LINE__,
"Expected an array.");
141 if (!d_latitude->read_p())
144 set_lat(extract_double_array(d_latitude));
145 set_lat_length(d_latitude->length());
151 && unit_or_name_match(get_coards_lon_units(), get_lon_names(),
152 units_value, map_name)) {
154 d_longitude =
dynamic_cast < Array *
>(*m);
156 throw InternalErr(__FILE__, __LINE__,
"Expected an array.");
157 if (!d_longitude->read_p())
160 set_lon(extract_double_array(d_longitude));
161 set_lon_length(d_longitude->length());
165 if (m + 1 == d_grid->map_end())
166 set_longitude_rightmost(
true);
173 return get_lat() && get_lon();
183 bool GridGeoConstraint::build_lat_lon_maps(Array *lat, Array *lon)
185 Grid::Map_iter m = d_grid->map_begin();
187 Array::Dim_iter d = d_grid->get_array()->dim_begin();
189 while (m != d_grid->map_end() && (!d_latitude || !d_longitude)) {
191 if (!d_latitude && *m == lat) {
195 if (!d_latitude->read_p())
198 set_lat(extract_double_array(d_latitude));
199 set_lat_length(d_latitude->length());
204 if (!d_longitude && *m == lon) {
208 if (!d_longitude->read_p())
211 set_lon(extract_double_array(d_longitude));
212 set_lon_length(d_longitude->length());
216 if (m + 1 == d_grid->map_end())
217 set_longitude_rightmost(
true);
224 return get_lat() && get_lon();
238 GridGeoConstraint::lat_lon_dimensions_ok()
241 Grid::Map_riter rightmost = d_grid->map_rbegin();
242 Grid::Map_riter next_rightmost = rightmost + 1;
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);
277 if (!is_bounding_box_set())
278 throw InternalErr(
"The Latitude and Longitude constraints must be set before calling apply_constraint_to_data().");
280 Array::Dim_iter fd = d_latitude->dim_begin();
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);
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.");
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());
303 if (get_longitude_index_left() > get_longitude_index_right()) {
319 set_longitude_index_right(get_lon_length() - get_longitude_index_left()
320 + get_longitude_index_right());
321 set_longitude_index_left(0);
335 if (get_longitude_notation() == neg_pos) {
340 fd = d_longitude->dim_begin();
341 d_longitude->add_constraint(fd, get_longitude_index_left(), 1,
342 get_longitude_index_right());
344 d_grid->get_array()->add_constraint(get_lon_dim(),
345 get_longitude_index_left(),
346 1, get_longitude_index_right());
351 if (get_latitude_sense() == inverted) {
352 DBG(cerr <<
"Inverted latitude sense" << endl);
354 get_latitude_index_bottom() - get_latitude_index_top() + 1);
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);
361 set_array_using_double(d_latitude, get_lat() + get_latitude_index_top(),
362 get_latitude_index_bottom() - get_latitude_index_top() + 1);
364 set_array_using_double(d_longitude, get_lon() + get_longitude_index_left(),
365 get_longitude_index_right() - get_longitude_index_left() + 1);
368 Grid::Map_iter i = d_grid->map_begin();
369 Grid::Map_iter end = d_grid->map_end();
371 if (*i != d_latitude && *i != d_longitude) {
372 if ((*i)->send_p()) {
373 DBG(cerr <<
"reading grid map: " << (*i)->name() << endl);
382 if (get_array_data()) {
383 int size = d_grid->get_array()->val2buf(get_array_data());
385 if (size != get_array_data_size())
386 throw InternalErr(__FILE__, __LINE__,
"Expected data size not copied to the Grid's buffer.");
388 d_grid->set_read_p(
true);
391 d_grid->get_array()->read();