Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "config.h"
00030
00031 static char id[] not_used =
00032 { "$Id: GridGeoConstraint.cc 21922 2010-01-07 18:23:57Z jimg $"
00033 };
00034
00035 #include <cmath>
00036
00037 #include <iostream>
00038 #include <sstream>
00039
00040
00041
00042 #include "debug.h"
00043 #include "dods-datatypes.h"
00044 #include "GridGeoConstraint.h"
00045 #include "Float64.h"
00046
00047 #include "Error.h"
00048 #include "InternalErr.h"
00049 #include "ce_functions.h"
00050 #include "util.h"
00051
00052 using namespace std;
00053
00054 namespace libdap {
00055
00064 GridGeoConstraint::GridGeoConstraint(Grid *grid)
00065 : GeoConstraint(), d_grid(grid), d_latitude(0), d_longitude(0)
00066 {
00067 if (d_grid->get_array()->dimensions() < 2
00068 || d_grid->get_array()->dimensions() > 3)
00069 throw Error("The geogrid() function works only with Grids of two or three dimensions.");
00070
00071
00072 if (!build_lat_lon_maps())
00073 throw Error(string("The grid '") + d_grid->name()
00074 +
00075 "' does not have identifiable latitude/longitude map vectors.");
00076
00077 if (!lat_lon_dimensions_ok())
00078 throw Error("The geogrid() function will only work when the Grid's Longitude and Latitude\nmaps are the rightmost dimensions.");
00079 }
00080
00096 bool GridGeoConstraint::build_lat_lon_maps()
00097 {
00098 Grid::Map_iter m = d_grid->map_begin();
00099
00100
00101
00102 Array::Dim_iter d = d_grid->get_array()->dim_begin();
00103
00104 while (m != d_grid->map_end() && (!d_latitude || !d_longitude)) {
00105 string units_value = (*m)->get_attr_table().get_attr("units");
00106 units_value = remove_quotes(units_value);
00107 string map_name = (*m)->name();
00108
00109
00110
00111 if (!d_latitude
00112 && unit_or_name_match(get_coards_lat_units(), get_lat_names(),
00113 units_value, map_name)) {
00114
00115
00116
00117
00118
00119
00120
00121 d_latitude = dynamic_cast < Array * >(*m);
00122 if (!d_latitude)
00123 throw InternalErr(__FILE__, __LINE__, "Expected an array.");
00124 if (!d_latitude->read_p())
00125 d_latitude->read();
00126
00127 set_lat(extract_double_array(d_latitude));
00128 set_lat_length(d_latitude->length());
00129
00130 set_lat_dim(d);
00131 }
00132
00133 if (!d_longitude
00134 && unit_or_name_match(get_coards_lon_units(), get_lon_names(),
00135 units_value, map_name)) {
00136
00137 d_longitude = dynamic_cast < Array * >(*m);
00138 if (!d_longitude)
00139 throw InternalErr(__FILE__, __LINE__, "Expected an array.");
00140 if (!d_longitude->read_p())
00141 d_longitude->read();
00142
00143 set_lon(extract_double_array(d_longitude));
00144 set_lon_length(d_longitude->length());
00145
00146 set_lon_dim(d);
00147
00148 if (m + 1 == d_grid->map_end())
00149 set_longitude_rightmost(true);
00150 }
00151
00152 ++m;
00153 ++d;
00154 }
00155
00156 return get_lat() && get_lon();
00157 }
00158
00169 bool
00170 GridGeoConstraint::lat_lon_dimensions_ok()
00171 {
00172
00173 Grid::Map_riter rightmost = d_grid->map_rbegin();
00174 Grid::Map_riter next_rightmost = rightmost + 1;
00175
00176 if (*rightmost == d_longitude && *next_rightmost == d_latitude)
00177 set_longitude_rightmost(true);
00178 else if (*rightmost == d_latitude && *next_rightmost == d_longitude)
00179 set_longitude_rightmost(false);
00180 else
00181 return false;
00182
00183 return true;
00184 }
00185
00207 void GridGeoConstraint::apply_constraint_to_data()
00208 {
00209 if (!is_bounding_box_set())
00210 throw InternalErr("The Latitude and Longitude constraints must be set before calling apply_constraint_to_data().");
00211
00212 Array::Dim_iter fd = d_latitude->dim_begin();
00213
00214 if (get_latitude_sense() == inverted) {
00215 int tmp = get_latitude_index_top();
00216 set_latitude_index_top(get_latitude_index_bottom());
00217 set_latitude_index_bottom(tmp);
00218 }
00219
00220
00221
00222 if (get_latitude_index_top() > get_latitude_index_bottom())
00223 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.");
00224
00225
00226 d_latitude->add_constraint(fd, get_latitude_index_top(), 1,
00227 get_latitude_index_bottom());
00228 d_grid->get_array()->add_constraint(get_lat_dim(),
00229 get_latitude_index_top(), 1,
00230 get_latitude_index_bottom());
00231
00232
00233
00234
00235 if (get_longitude_index_left() > get_longitude_index_right()) {
00236 reorder_longitude_map(get_longitude_index_left());
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246 reorder_data_longitude_axis(*d_grid->get_array(), get_lon_dim());
00247
00248
00249
00250
00251 set_longitude_index_right(get_lon_length() - get_longitude_index_left()
00252 + get_longitude_index_right());
00253 set_longitude_index_left(0);
00254 }
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267 if (get_longitude_notation() == neg_pos) {
00268 transform_longitude_to_neg_pos_notation();
00269 }
00270
00271
00272 fd = d_longitude->dim_begin();
00273 d_longitude->add_constraint(fd, get_longitude_index_left(), 1,
00274 get_longitude_index_right());
00275
00276 d_grid->get_array()->add_constraint(get_lon_dim(),
00277 get_longitude_index_left(),
00278 1, get_longitude_index_right());
00279
00280
00281
00282
00283 if (get_latitude_sense() == inverted) {
00284 DBG(cerr << "Inverted latitude sense" << endl);
00285 transpose_vector(get_lat() + get_latitude_index_top(),
00286 get_latitude_index_bottom() - get_latitude_index_top() + 1);
00287
00288 flip_latitude_within_array(*d_grid->get_array(),
00289 get_latitude_index_bottom() - get_latitude_index_top() + 1,
00290 get_longitude_index_right() - get_longitude_index_left() + 1);
00291 }
00292
00293 set_array_using_double(d_latitude, get_lat() + get_latitude_index_top(),
00294 get_latitude_index_bottom() - get_latitude_index_top() + 1);
00295
00296 set_array_using_double(d_longitude, get_lon() + get_longitude_index_left(),
00297 get_longitude_index_right() - get_longitude_index_left() + 1);
00298
00299
00300 Grid::Map_iter i = d_grid->map_begin();
00301 Grid::Map_iter end = d_grid->map_end();
00302 while (i != end) {
00303 if (*i != d_latitude && *i != d_longitude) {
00304 if ((*i)->send_p()) {
00305 DBG(cerr << "reading grid map: " << (*i)->name() << endl);
00306
00307 (*i)->read();
00308 }
00309 }
00310 ++i;
00311 }
00312
00313
00314 if (get_array_data()) {
00315 int size = d_grid->get_array()->val2buf(get_array_data());
00316
00317 if (size != get_array_data_size())
00318 throw InternalErr(__FILE__, __LINE__, "Expected data size not copied to the Grid's buffer.");
00319
00320 d_grid->set_read_p(true);
00321 }
00322 else {
00323 d_grid->get_array()->read();
00324 }
00325 }
00326
00327 }