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: GeoConstraint.cc 18315 2008-03-03 20:14:44Z jimg $"
00033 };
00034
00035 #include <cstring>
00036 #include <cmath>
00037 #include <iostream>
00038 #include <sstream>
00039 #include <algorithm>
00040
00041
00042
00043 #include "debug.h"
00044 #include "dods-datatypes.h"
00045 #include "GeoConstraint.h"
00046 #include "Float64.h"
00047
00048 #include "Error.h"
00049 #include "InternalErr.h"
00050 #include "ce_functions.h"
00051 #include "util.h"
00052
00053 using namespace std;
00054
00055 namespace libdap {
00056
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00080 class is_prefix
00081 {
00082 private:
00083 string s;
00084 public:
00085 is_prefix(const string & in): s(in)
00086 {}
00087
00088 bool operator()(const string & prefix)
00089 {
00090 return s.find(prefix) == 0;
00091 }
00092 };
00093
00104 bool
00105 unit_or_name_match(set < string > units, set < string > names,
00106 const string & var_units, const string & var_name)
00107 {
00108 return (units.find(var_units) != units.end()
00109 || find_if(names.begin(), names.end(),
00110 is_prefix(var_name)) != names.end());
00111 }
00112
00127 GeoConstraint::Notation
00128 GeoConstraint::categorize_notation(double left,
00129 double right) const
00130 {
00131 return (left < 0 || right < 0) ? neg_pos : pos;
00132 }
00133
00134
00135
00136
00137
00138 void
00139 GeoConstraint::transform_constraint_to_pos_notation(double &left,
00140 double &right) const
00141 {
00142 if (left < 0 || right < 0) {
00143 left += 180;
00144 right += 180;
00145 }
00146 }
00147
00151 void GeoConstraint::transform_longitude_to_pos_notation()
00152 {
00153
00154
00155
00156 for (int i = 0; i < d_lon_length; ++i)
00157 d_lon[i] += 180;
00158 }
00159
00163 void GeoConstraint::transform_longitude_to_neg_pos_notation()
00164 {
00165 for (int i = 0; i < d_lon_length; ++i)
00166 d_lon[i] -= 180;
00167 }
00168
00169 bool GeoConstraint::is_bounding_box_valid(double left, double top,
00170 double right, double bottom) const
00171 {
00172 if ((left < d_lon[0] && right < d_lon[0])
00173 || (left > d_lon[d_lon_length-1] && right > d_lon[d_lon_length-1]))
00174 return false;
00175
00176 if (d_latitude_sense == normal) {
00177
00178 if ((top > d_lat[0] && bottom > d_lat[0])
00179 || (top < d_lat[d_lat_length-1] && bottom < d_lat[d_lat_length-1]))
00180 return false;
00181 }
00182 else {
00183 if ((top < d_lat[0] && bottom < d_lat[0])
00184 || (top > d_lat[d_lat_length-1] && bottom > d_lat[d_lat_length-1]))
00185 return false;
00186 }
00187
00188 return true;
00189 }
00190
00201 void GeoConstraint::find_longitude_indeces(double left, double right,
00202 int &longitude_index_left,
00203 int &longitude_index_right) const
00204 {
00205
00206
00207
00208 double t_left = fmod(left, 360.0);
00209 double t_right = fmod(right, 360.0);
00210
00211
00212
00213
00214
00215 int i = 0;
00216 int smallest_lon_index = 0;
00217 double smallest_lon = fmod(d_lon[smallest_lon_index], 360.0);
00218 while (i < d_lon_length - 1) {
00219 if (smallest_lon > fmod(d_lon[i], 360.0)) {
00220 smallest_lon_index = i;
00221 smallest_lon = fmod(d_lon[smallest_lon_index], 360.0);
00222 }
00223 ++i;
00224 }
00225 DBG2(cerr << "smallest_lon_index: " << smallest_lon_index << endl);
00226
00227
00228
00229
00230 i = smallest_lon_index;
00231 bool done = false;
00232
00233 DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) < t_left: "
00234 << fmod(d_lon[i], 360.0) << " < " << t_left << endl);
00235
00236 while (!done && fmod(d_lon[i], 360.0) < t_left) {
00237
00238 DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) < t_left: "
00239 << fmod(d_lon[i], 360.0) << " < " << t_left << endl);
00240
00241 ++i;
00242 i = i % d_lon_length;
00243 if (i == smallest_lon_index)
00244 done = true;
00245 }
00246 if (fmod(d_lon[i], 360.0) == t_left)
00247 longitude_index_left = i;
00248 else
00249 longitude_index_left = (i - 1) > 0 ? i - 1 : 0;
00250
00251
00252
00253 int largest_lon_index =
00254 (smallest_lon_index - 1 + d_lon_length) % d_lon_length;
00255 DBG2(cerr << "largest_lon_index: " << largest_lon_index << endl);
00256 i = largest_lon_index;
00257 done = false;
00258
00259 DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) > t_right: "
00260 << fmod(d_lon[i], 360.0) << " > " << t_right << endl);
00261
00262 while (!done && fmod(d_lon[i], 360.0) > t_right) {
00263
00264 DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) > t_right: "
00265 << fmod(d_lon[i], 360.0) << " > " << t_right << endl);
00266
00267
00268 i = (i == 0) ? d_lon_length - 1 : i - 1;
00269 if (i == largest_lon_index)
00270 done = true;
00271 }
00272 if (fmod(d_lon[i], 360.0) == t_right)
00273 longitude_index_right = i;
00274 else
00275 longitude_index_right =
00276 (i + 1) < d_lon_length - 1 ? i + 1 : d_lon_length - 1;
00277 }
00278
00291 void GeoConstraint::find_latitude_indeces(double top, double bottom,
00292 LatitudeSense sense,
00293 int &latitude_index_top,
00294 int &latitude_index_bottom) const
00295 {
00296
00297 int i = 0;
00298 if (sense == normal) {
00299 while (i < d_lat_length - 1 && d_lat[i] > top)
00300 ++i;
00301 if (d_lat[i] == top)
00302 latitude_index_top = i;
00303 else
00304 latitude_index_top = (i - 1) > 0 ? i - 1 : 0;
00305 }
00306 else {
00307 while (i < d_lat_length - 1 && d_lat[i] < top)
00308 ++i;
00309 latitude_index_top = i;
00310 }
00311
00312
00313
00314 i = d_lat_length - 1;
00315 if (sense == normal) {
00316 while (i > 0 && d_lat[i] < bottom)
00317 --i;
00318 if (d_lat[i] == bottom)
00319 latitude_index_bottom = i;
00320 else
00321 latitude_index_bottom =
00322 (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1;
00323 }
00324 else {
00325 while (i > 0 && d_lat[i] > bottom)
00326 --i;
00327 latitude_index_bottom = i;
00328 }
00329 }
00330
00334 GeoConstraint::LatitudeSense GeoConstraint::categorize_latitude() const
00335 {
00336 return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted;
00337 }
00338
00339 #if 0
00340
00347 void GeoConstraint::set_bounding_box_latitude(double top, double bottom)
00348 {
00349
00350 d_latitude_sense = categorize_latitude();
00351
00352
00353
00354
00355
00356 find_latitude_indeces(top, bottom, d_latitude_sense,
00357 d_latitude_index_top, d_latitude_index_bottom);
00358 }
00359 #endif
00360
00361
00362
00363 static void
00364 swap_vector_ends(char *dest, char *src, int len, int index, int elem_sz)
00365 {
00366 memcpy(dest, src + index * elem_sz, (len - index) * elem_sz);
00367
00368 memcpy(dest + (len - index) * elem_sz, src, index * elem_sz);
00369 }
00370
00379 void GeoConstraint::reorder_longitude_map(int longitude_index_left)
00380 {
00381 double *tmp_lon = new double[d_lon_length];
00382
00383 swap_vector_ends((char *) tmp_lon, (char *) d_lon, d_lon_length,
00384 longitude_index_left, sizeof(double));
00385
00386 memcpy(d_lon, tmp_lon, d_lon_length * sizeof(double));
00387
00388 delete[]tmp_lon;
00389 }
00390
00391 static int
00392 count_dimensions_except_longitude(Array &a)
00393 {
00394 int size = 1;
00395 for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i)
00396 size *= a.dimension_size(i, true);
00397
00398 return size;
00399 }
00400
00419 void GeoConstraint::reorder_data_longitude_axis(Array &a)
00420 {
00421
00422 if (!get_longitude_rightmost())
00423 throw Error("This grid does not have Longitude as its rightmost dimension, the geogrid()\ndoes not support constraints that wrap around the edges of this type of grid.");
00424
00425 DBG(cerr << "Constraint for the left half: " << get_longitude_index_left()
00426 << ", " << get_lon_length() - 1 << endl);
00427
00428 a.add_constraint(d_lon_dim, get_longitude_index_left(), 1,
00429 get_lon_length() - 1);
00430 a.set_read_p(false);
00431 a.read(get_dataset());
00432 DBG2(a.print_val(stderr));
00433 #if 0
00434 char *left_data = 0;
00435 int left_size = a.buf2val((void **) & left_data);
00436 #endif
00437 char *left_data = (char*)a.value();
00438 int left_size = a.length();
00439
00440
00441
00442
00443
00444
00445 DBG(cerr << "Constraint for the right half: " << 0
00446 << ", " << get_longitude_index_right() << endl);
00447
00448 a.add_constraint(d_lon_dim, 0, 1, get_longitude_index_right());
00449 a.set_read_p(false);
00450 a.read(get_dataset());
00451 DBG2(a.print_val(stderr));
00452 #if 0
00453 char *right_data = 0;
00454 int right_size = a.buf2val((void **) & right_data);
00455 #endif
00456 char *right_data = (char*)a.value();
00457 int right_size = a.length();
00458
00459
00460 d_array_data_size = left_size + right_size;
00461 d_array_data = new char[d_array_data_size];
00462
00463
00464
00465
00466 int elem_width = a.var()->width();
00467 int left_elements = (get_lon_length() - get_longitude_index_left()) * elem_width;
00468 int right_elements = (get_longitude_index_right() + 1) * elem_width;
00469 int total_elements_per_row = left_elements + right_elements;
00470
00471
00472 int rows_to_copy = count_dimensions_except_longitude(a);
00473 for (int i = 0; i < rows_to_copy; ++i) {
00474 memcpy(d_array_data + (total_elements_per_row * i),
00475 left_data + (left_elements * i),
00476 left_elements);
00477 memcpy(d_array_data + left_elements + (total_elements_per_row * i),
00478 right_data + (right_elements * i),
00479 right_elements);
00480 }
00481
00482 delete[]left_data;
00483 delete[]right_data;
00484 }
00485
00486 #if 0
00487
00500 void GeoConstraint::set_bounding_box_longitude(double left, double right)
00501 {
00502
00503 d_longitude_notation = categorize_notation(left, right);
00504
00505
00506 if (d_longitude_notation == neg_pos)
00507 transform_constraint_to_pos_notation(left, right);
00508
00509
00510
00511
00512 Notation longitude_notation =
00513 categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
00514
00515 if (longitude_notation == neg_pos)
00516 transform_longitude_to_pos_notation();
00517
00518
00519 find_longitude_indeces(left, right, d_longitude_index_left,
00520 d_longitude_index_right);
00521 }
00522 #endif
00523
00528 GeoConstraint::GeoConstraint(const string & ds_name)
00529 : d_dataset(ds_name), d_array_data(0), d_array_data_size(0),
00530 d_lat(0), d_lon(0),
00531 d_bounding_box_set(false),
00532 d_longitude_notation(unknown_notation),
00533 d_latitude_sense(unknown_sense)
00534 {
00535
00536 d_coards_lat_units.insert("degrees_north");
00537 d_coards_lat_units.insert("degree_north");
00538 d_coards_lat_units.insert("degree_N");
00539 d_coards_lat_units.insert("degrees_N");
00540
00541 d_coards_lon_units.insert("degrees_east");
00542 d_coards_lon_units.insert("degree_east");
00543 d_coards_lon_units.insert("degrees_E");
00544 d_coards_lon_units.insert("degree_E");
00545
00546 d_lat_names.insert("COADSY");
00547 d_lat_names.insert("lat");
00548 d_lat_names.insert("Lat");
00549 d_lat_names.insert("LAT");
00550
00551 d_lon_names.insert("COADSX");
00552 d_lon_names.insert("lon");
00553 d_lon_names.insert("Lon");
00554 d_lon_names.insert("LON");
00555 }
00556
00567 void GeoConstraint::set_bounding_box(double left, double top,
00568 double right, double bottom)
00569 {
00570
00571
00572
00573 if (d_bounding_box_set)
00574 throw
00575 InternalErr
00576 ("It is not possible to register more than one geographical constraint on a variable.");
00577
00578
00579 d_latitude_sense = categorize_latitude();
00580
00581
00582 d_longitude_notation = categorize_notation(left, right);
00583
00584
00585 if (d_longitude_notation == neg_pos)
00586 transform_constraint_to_pos_notation(left, right);
00587
00588
00589
00590
00591 Notation longitude_notation =
00592 categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
00593
00594 if (longitude_notation == neg_pos)
00595 transform_longitude_to_pos_notation();
00596
00597 if (!is_bounding_box_valid(left, top, right, bottom))
00598 throw Error("The bounding box does not intersect any data within this Grid or Array. The\ngeographical extent of these data are from latitude "
00599 + double_to_string(d_lat[0]) + " to "
00600 + double_to_string(d_lat[d_lat_length-1])
00601 + "\nand longitude " + double_to_string(d_lon[0])
00602 + " to " + double_to_string(d_lon[d_lon_length-1])
00603 + " while the bounding box provided was latitude "
00604 + double_to_string(top) + " to "
00605 + double_to_string(bottom) + "\nand longitude "
00606 + double_to_string(left) + " to "
00607 + double_to_string(right));
00608
00609
00610
00611
00612
00613 find_latitude_indeces(top, bottom, d_latitude_sense,
00614 d_latitude_index_top, d_latitude_index_bottom);
00615
00616
00617
00618 find_longitude_indeces(left, right, d_longitude_index_left,
00619 d_longitude_index_right);
00620
00621 DBG(cerr << "Bounding box (tlbr): " << d_latitude_index_top << ", "
00622 << d_longitude_index_left << ", "
00623 << d_latitude_index_bottom << ", "
00624 << d_longitude_index_right << endl);
00625
00626 d_bounding_box_set = true;
00627 }
00628
00629 }