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 16088 2007-03-28 21:42:19Z jimg $"
00033 };
00034
00035 #include <math.h>
00036 #include <string.h>
00037
00038 #include <iostream>
00039 #include <sstream>
00040
00041
00042
00043 #include "debug.h"
00044 #include "dods-datatypes.h"
00045 #include "GridGeoConstraint.h"
00046 #include "Float64.h"
00047
00048 #include "Error.h"
00049 #include "InternalErr.h"
00050 #include "ce_functions.h"
00051
00052 using namespace std;
00053 using namespace libdap;
00054
00055
00063 GridGeoConstraint::GridGeoConstraint(Grid *grid, const string &ds_name)
00064 : GeoConstraint(ds_name), d_grid(grid), d_latitude(0), d_longitude(0)
00065 {
00066 if (d_grid->get_array()->dimensions() < 2
00067 || d_grid->get_array()->dimensions() > 3)
00068 throw Error("The geogrid() function works only with Grids of two or three dimensions.");
00069
00070
00071 if (!build_lat_lon_maps())
00072 throw Error(string("The grid '") + d_grid->name()
00073 +
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
00095 bool GridGeoConstraint::build_lat_lon_maps()
00096 {
00097 Grid::Map_iter m = d_grid->map_begin();
00098
00099
00100
00101 Array::Dim_iter d = d_grid->get_array()->dim_begin();
00102
00103 while (m != d_grid->map_end() && (!d_latitude || !d_longitude)) {
00104 string units_value = (*m)->get_attr_table().get_attr("units");
00105 remove_quotes(units_value);
00106 string map_name = (*m)->name();
00107
00108
00109
00110 if (!d_latitude
00111 && unit_or_name_match(get_coards_lat_units(), get_lat_names(),
00112 units_value, map_name)) {
00113
00114
00115
00116
00117
00118
00119
00120 d_latitude = dynamic_cast < Array * >(*m);
00121 if (!d_latitude->read_p())
00122 d_latitude->read(get_dataset());
00123
00124 set_lat(extract_double_array(d_latitude));
00125 set_lat_length(d_latitude->length());
00126
00127 set_lat_dim(d);
00128 }
00129
00130 if (!d_longitude
00131 && unit_or_name_match(get_coards_lon_units(), get_lon_names(),
00132 units_value, map_name)) {
00133
00134 d_longitude = dynamic_cast < Array * >(*m);
00135 if (!d_longitude->read_p())
00136 d_longitude->read(get_dataset());
00137
00138 set_lon(extract_double_array(d_longitude));
00139 set_lon_length(d_longitude->length());
00140
00141 set_lon_dim(d);
00142 }
00143
00144 ++m;
00145 ++d;
00146 }
00147
00148 return get_lat() && get_lon();
00149 }
00150
00161 bool
00162 GridGeoConstraint::lat_lon_dimensions_ok()
00163 {
00164
00165 Grid::Map_riter rightmost = d_grid->map_rbegin();
00166 Grid::Map_riter next_rightmost = rightmost + 1;
00167
00168 if (*rightmost == d_longitude && *next_rightmost == d_latitude)
00169 set_longitude_rightmost(true);
00170 else if (*rightmost == d_latitude && *next_rightmost == d_longitude)
00171 set_longitude_rightmost(false);
00172 else
00173 return false;
00174
00175 return true;
00176 }
00177
00199 void GridGeoConstraint::apply_constraint_to_data()
00200 {
00201 if (!get_bounding_box_set())
00202 throw
00203 InternalErr
00204 ("The Latitude and Longitude constraints must be set before calling apply_constraint_to_data().");
00205
00206 Array::Dim_iter fd = d_latitude->dim_begin();
00207 if (get_latitude_sense() == inverted) {
00208 int tmp = get_latitude_index_top();
00209 set_latitude_index_top(get_latitude_index_bottom());
00210 set_latitude_index_bottom(tmp);
00211 }
00212
00213
00214
00215 if (get_latitude_index_top() > get_latitude_index_bottom())
00216 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.");
00217
00218
00219 d_latitude->add_constraint(fd, get_latitude_index_top(), 1,
00220 get_latitude_index_bottom());
00221 d_grid->get_array()->add_constraint(get_lat_dim(),
00222 get_latitude_index_top(), 1,
00223 get_latitude_index_bottom());
00224
00225
00226
00227
00228 if (get_longitude_index_left() > get_longitude_index_right()) {
00229 reorder_longitude_map(get_longitude_index_left());
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239 reorder_data_longitude_axis(*d_grid->get_array());
00240
00241
00242
00243
00244
00245 set_longitude_index_right(get_lon_length() - get_longitude_index_left()
00246 + get_longitude_index_right());
00247 set_longitude_index_left(0);
00248 }
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260 if (get_longitude_notation() == neg_pos) {
00261 transform_longitude_to_neg_pos_notation();
00262 }
00263
00264 fd = d_longitude->dim_begin();
00265 d_longitude->add_constraint(fd, get_longitude_index_left(), 1,
00266 get_longitude_index_right());
00267
00268 d_grid->get_array()->add_constraint(get_lon_dim(),
00269 get_longitude_index_left(),
00270 1, get_longitude_index_right());
00271
00272
00273 set_array_using_double(d_latitude, get_lat() + get_latitude_index_top(),
00274 get_latitude_index_bottom() - get_latitude_index_top() + 1);
00275
00276 set_array_using_double(d_longitude, get_lon() + get_longitude_index_left(),
00277 get_longitude_index_right() - get_longitude_index_left() + 1);
00278
00279
00280 if (get_array_data()) {
00281 int size = d_grid->get_array()->val2buf(get_array_data());
00282 if (size != get_array_data_size())
00283 throw
00284 InternalErr
00285 ("Expected data size not copied to the Grid's buffer.");
00286 d_grid->set_read_p(true);
00287 }
00288 else {
00289 d_grid->get_array()->read(get_dataset());
00290 }
00291 }