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 23447 2010-08-27 19:38:04Z 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 + "' does not have identifiable latitude/longitude map vectors.");
00075
00076 if (!lat_lon_dimensions_ok())
00077 throw Error("The geogrid() function will only work when the Grid's Longitude and Latitude\nmaps are the rightmost dimensions.");
00078 }
00079
00080 GridGeoConstraint::GridGeoConstraint(Grid *grid, Array *lat, Array *lon)
00081 : GeoConstraint(), d_grid(grid), d_latitude(0), d_longitude(0)
00082 {
00083 if (d_grid->get_array()->dimensions() < 2
00084 || d_grid->get_array()->dimensions() > 3)
00085 throw Error("The geogrid() function works only with Grids of two or three dimensions.");
00086
00087
00088 if (!build_lat_lon_maps(lat, lon))
00089 throw Error(string("The grid '") + d_grid->name()
00090 + "' does not have valid latitude/longitude map vectors.");
00091
00092
00093 if (!lat_lon_dimensions_ok())
00094 throw Error("The geogrid() function will only work when the Grid's Longitude and Latitude\nmaps are the rightmost dimensions.");
00095 }
00096
00112 bool GridGeoConstraint::build_lat_lon_maps()
00113 {
00114 Grid::Map_iter m = d_grid->map_begin();
00115
00116
00117
00118 Array::Dim_iter d = d_grid->get_array()->dim_begin();
00119
00120
00121
00122
00123
00124 while (m != d_grid->map_end() && (!d_latitude || !d_longitude)) {
00125 string units_value = (*m)->get_attr_table().get_attr("units");
00126 units_value = remove_quotes(units_value);
00127 string map_name = (*m)->name();
00128
00129
00130
00131 if (!d_latitude
00132 && unit_or_name_match(get_coards_lat_units(), get_lat_names(),
00133 units_value, map_name)) {
00134
00135
00136
00137
00138
00139
00140
00141 d_latitude = dynamic_cast < Array * >(*m);
00142 if (!d_latitude)
00143 throw InternalErr(__FILE__, __LINE__, "Expected an array.");
00144 if (!d_latitude->read_p())
00145 d_latitude->read();
00146
00147 set_lat(extract_double_array(d_latitude));
00148 set_lat_length(d_latitude->length());
00149
00150 set_lat_dim(d);
00151 }
00152
00153 if (!d_longitude
00154 && unit_or_name_match(get_coards_lon_units(), get_lon_names(),
00155 units_value, map_name)) {
00156
00157 d_longitude = dynamic_cast < Array * >(*m);
00158 if (!d_longitude)
00159 throw InternalErr(__FILE__, __LINE__, "Expected an array.");
00160 if (!d_longitude->read_p())
00161 d_longitude->read();
00162
00163 set_lon(extract_double_array(d_longitude));
00164 set_lon_length(d_longitude->length());
00165
00166 set_lon_dim(d);
00167
00168 if (m + 1 == d_grid->map_end())
00169 set_longitude_rightmost(true);
00170 }
00171
00172 ++m;
00173 ++d;
00174 }
00175
00176 return get_lat() && get_lon();
00177 }
00178
00186 bool GridGeoConstraint::build_lat_lon_maps(Array *lat, Array *lon)
00187 {
00188 Grid::Map_iter m = d_grid->map_begin();
00189
00190 Array::Dim_iter d = d_grid->get_array()->dim_begin();
00191
00192 while (m != d_grid->map_end() && (!d_latitude || !d_longitude)) {
00193
00194 if (!d_latitude && *m == lat) {
00195
00196 d_latitude = lat;
00197
00198 if (!d_latitude->read_p())
00199 d_latitude->read();
00200
00201 set_lat(extract_double_array(d_latitude));
00202 set_lat_length(d_latitude->length());
00203
00204 set_lat_dim(d);
00205 }
00206
00207 if (!d_longitude && *m == lon) {
00208
00209 d_longitude = lon;
00210
00211 if (!d_longitude->read_p())
00212 d_longitude->read();
00213
00214 set_lon(extract_double_array(d_longitude));
00215 set_lon_length(d_longitude->length());
00216
00217 set_lon_dim(d);
00218
00219 if (m + 1 == d_grid->map_end())
00220 set_longitude_rightmost(true);
00221 }
00222
00223 ++m;
00224 ++d;
00225 }
00226
00227 return get_lat() && get_lon();
00228 }
00229
00240 bool
00241 GridGeoConstraint::lat_lon_dimensions_ok()
00242 {
00243
00244 Grid::Map_riter rightmost = d_grid->map_rbegin();
00245 Grid::Map_riter next_rightmost = rightmost + 1;
00246
00247 if (*rightmost == d_longitude && *next_rightmost == d_latitude)
00248 set_longitude_rightmost(true);
00249 else if (*rightmost == d_latitude && *next_rightmost == d_longitude)
00250 set_longitude_rightmost(false);
00251 else
00252 return false;
00253
00254 return true;
00255 }
00256
00278 void GridGeoConstraint::apply_constraint_to_data()
00279 {
00280 if (!is_bounding_box_set())
00281 throw InternalErr("The Latitude and Longitude constraints must be set before calling apply_constraint_to_data().");
00282
00283 Array::Dim_iter fd = d_latitude->dim_begin();
00284
00285 if (get_latitude_sense() == inverted) {
00286 int tmp = get_latitude_index_top();
00287 set_latitude_index_top(get_latitude_index_bottom());
00288 set_latitude_index_bottom(tmp);
00289 }
00290
00291
00292
00293 if (get_latitude_index_top() > get_latitude_index_bottom())
00294 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.");
00295
00296
00297 d_latitude->add_constraint(fd, get_latitude_index_top(), 1,
00298 get_latitude_index_bottom());
00299 d_grid->get_array()->add_constraint(get_lat_dim(),
00300 get_latitude_index_top(), 1,
00301 get_latitude_index_bottom());
00302
00303
00304
00305
00306 if (get_longitude_index_left() > get_longitude_index_right()) {
00307 reorder_longitude_map(get_longitude_index_left());
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317 reorder_data_longitude_axis(*d_grid->get_array(), get_lon_dim());
00318
00319
00320
00321
00322 set_longitude_index_right(get_lon_length() - get_longitude_index_left()
00323 + get_longitude_index_right());
00324 set_longitude_index_left(0);
00325 }
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338 if (get_longitude_notation() == neg_pos) {
00339 transform_longitude_to_neg_pos_notation();
00340 }
00341
00342
00343 fd = d_longitude->dim_begin();
00344 d_longitude->add_constraint(fd, get_longitude_index_left(), 1,
00345 get_longitude_index_right());
00346
00347 d_grid->get_array()->add_constraint(get_lon_dim(),
00348 get_longitude_index_left(),
00349 1, get_longitude_index_right());
00350
00351
00352
00353
00354 if (get_latitude_sense() == inverted) {
00355 DBG(cerr << "Inverted latitude sense" << endl);
00356 transpose_vector(get_lat() + get_latitude_index_top(),
00357 get_latitude_index_bottom() - get_latitude_index_top() + 1);
00358
00359 flip_latitude_within_array(*d_grid->get_array(),
00360 get_latitude_index_bottom() - get_latitude_index_top() + 1,
00361 get_longitude_index_right() - get_longitude_index_left() + 1);
00362 }
00363
00364 set_array_using_double(d_latitude, get_lat() + get_latitude_index_top(),
00365 get_latitude_index_bottom() - get_latitude_index_top() + 1);
00366
00367 set_array_using_double(d_longitude, get_lon() + get_longitude_index_left(),
00368 get_longitude_index_right() - get_longitude_index_left() + 1);
00369
00370
00371 Grid::Map_iter i = d_grid->map_begin();
00372 Grid::Map_iter end = d_grid->map_end();
00373 while (i != end) {
00374 if (*i != d_latitude && *i != d_longitude) {
00375 if ((*i)->send_p()) {
00376 DBG(cerr << "reading grid map: " << (*i)->name() << endl);
00377
00378 (*i)->read();
00379 }
00380 }
00381 ++i;
00382 }
00383
00384
00385 if (get_array_data()) {
00386 int size = d_grid->get_array()->val2buf(get_array_data());
00387
00388 if (size != get_array_data_size())
00389 throw InternalErr(__FILE__, __LINE__, "Expected data size not copied to the Grid's buffer.");
00390
00391 d_grid->set_read_p(true);
00392 }
00393 else {
00394 d_grid->get_array()->read();
00395 }
00396 }
00397
00398 }