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 }