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 19914 2008-11-25 17:26:22Z 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
00061 class is_prefix
00062 {
00063 private:
00064 string s;
00065 public:
00066 is_prefix(const string & in): s(in)
00067 {}
00068
00069 bool operator()(const string & prefix)
00070 {
00071 return s.find(prefix) == 0;
00072 }
00073 };
00074
00085 bool
00086 unit_or_name_match(set < string > units, set < string > names,
00087 const string & var_units, const string & var_name)
00088 {
00089 return (units.find(var_units) != units.end()
00090 || find_if(names.begin(), names.end(),
00091 is_prefix(var_name)) != names.end());
00092 }
00093
00108 GeoConstraint::Notation
00109 GeoConstraint::categorize_notation(double left,
00110 double right) const
00111 {
00112 return (left < 0 || right < 0) ? neg_pos : pos;
00113 }
00114
00115
00116
00117
00118
00119 void
00120 GeoConstraint::transform_constraint_to_pos_notation(double &left,
00121 double &right) const
00122 {
00123 if (left < 0 || right < 0) {
00124 left += 180;
00125 right += 180;
00126 }
00127 }
00128
00132 void GeoConstraint::transform_longitude_to_pos_notation()
00133 {
00134
00135
00136
00137 for (int i = 0; i < d_lon_length; ++i)
00138 d_lon[i] += 180;
00139 }
00140
00144 void GeoConstraint::transform_longitude_to_neg_pos_notation()
00145 {
00146 for (int i = 0; i < d_lon_length; ++i)
00147 d_lon[i] -= 180;
00148 }
00149
00150 bool GeoConstraint::is_bounding_box_valid(double left, double top,
00151 double right, double bottom) const
00152 {
00153 if ((left < d_lon[0] && right < d_lon[0])
00154 || (left > d_lon[d_lon_length-1] && right > d_lon[d_lon_length-1]))
00155 return false;
00156
00157 if (d_latitude_sense == normal) {
00158
00159 if ((top > d_lat[0] && bottom > d_lat[0])
00160 || (top < d_lat[d_lat_length-1] && bottom < d_lat[d_lat_length-1]))
00161 return false;
00162 }
00163 else {
00164 if ((top < d_lat[0] && bottom < d_lat[0])
00165 || (top > d_lat[d_lat_length-1] && bottom > d_lat[d_lat_length-1]))
00166 return false;
00167 }
00168
00169 return true;
00170 }
00171
00182 void GeoConstraint::find_longitude_indeces(double left, double right,
00183 int &longitude_index_left,
00184 int &longitude_index_right) const
00185 {
00186
00187
00188
00189 double t_left = fmod(left, 360.0);
00190 double t_right = fmod(right, 360.0);
00191
00192
00193
00194
00195
00196 int i = 0;
00197 int smallest_lon_index = 0;
00198 double smallest_lon = fmod(d_lon[smallest_lon_index], 360.0);
00199 while (i < d_lon_length - 1) {
00200 if (smallest_lon > fmod(d_lon[i], 360.0)) {
00201 smallest_lon_index = i;
00202 smallest_lon = fmod(d_lon[smallest_lon_index], 360.0);
00203 }
00204 ++i;
00205 }
00206 DBG2(cerr << "smallest_lon_index: " << smallest_lon_index << endl);
00207
00208
00209
00210
00211 i = smallest_lon_index;
00212 bool done = false;
00213
00214 DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) < t_left: "
00215 << fmod(d_lon[i], 360.0) << " < " << t_left << endl);
00216
00217 while (!done && fmod(d_lon[i], 360.0) < t_left) {
00218
00219 DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) < t_left: "
00220 << fmod(d_lon[i], 360.0) << " < " << t_left << endl);
00221
00222 ++i;
00223 i = i % d_lon_length;
00224 if (i == smallest_lon_index)
00225 done = true;
00226 }
00227 if (fmod(d_lon[i], 360.0) == t_left)
00228 longitude_index_left = i;
00229 else
00230 longitude_index_left = (i - 1) > 0 ? i - 1 : 0;
00231
00232
00233
00234 int largest_lon_index =
00235 (smallest_lon_index - 1 + d_lon_length) % d_lon_length;
00236 DBG2(cerr << "largest_lon_index: " << largest_lon_index << endl);
00237 i = largest_lon_index;
00238 done = false;
00239
00240 DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) > t_right: "
00241 << fmod(d_lon[i], 360.0) << " > " << t_right << endl);
00242
00243 while (!done && fmod(d_lon[i], 360.0) > t_right) {
00244
00245 DBG2(cerr << "fmod(d_lon[" << i << "], 360.0) > t_right: "
00246 << fmod(d_lon[i], 360.0) << " > " << t_right << endl);
00247
00248
00249 i = (i == 0) ? d_lon_length - 1 : i - 1;
00250 if (i == largest_lon_index)
00251 done = true;
00252 }
00253 if (fmod(d_lon[i], 360.0) == t_right)
00254 longitude_index_right = i;
00255 else
00256 longitude_index_right =
00257 (i + 1) < d_lon_length - 1 ? i + 1 : d_lon_length - 1;
00258 }
00259
00272 void GeoConstraint::find_latitude_indeces(double top, double bottom,
00273 LatitudeSense sense,
00274 int &latitude_index_top,
00275 int &latitude_index_bottom) const
00276 {
00277
00278 int i = 0;
00279 if (sense == normal) {
00280 while (i < d_lat_length - 1 && d_lat[i] > top)
00281 ++i;
00282 if (d_lat[i] == top)
00283 latitude_index_top = i;
00284 else
00285 latitude_index_top = (i - 1) > 0 ? i - 1 : 0;
00286 }
00287 else {
00288 while (i < d_lat_length - 1 && d_lat[i] < top)
00289 ++i;
00290 latitude_index_top = i;
00291 }
00292
00293
00294
00295 i = d_lat_length - 1;
00296 if (sense == normal) {
00297 while (i > 0 && d_lat[i] < bottom)
00298 --i;
00299 if (d_lat[i] == bottom)
00300 latitude_index_bottom = i;
00301 else
00302 latitude_index_bottom =
00303 (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1;
00304 }
00305 else {
00306 while (i > 0 && d_lat[i] > bottom)
00307 --i;
00308 latitude_index_bottom = i;
00309 }
00310 }
00311
00315 GeoConstraint::LatitudeSense GeoConstraint::categorize_latitude() const
00316 {
00317 return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted;
00318 }
00319
00320 #if 0
00321
00328 void GeoConstraint::set_bounding_box_latitude(double top, double bottom)
00329 {
00330
00331 d_latitude_sense = categorize_latitude();
00332
00333
00334
00335
00336
00337 find_latitude_indeces(top, bottom, d_latitude_sense,
00338 d_latitude_index_top, d_latitude_index_bottom);
00339 }
00340 #endif
00341
00342
00343
00344 static void
00345 swap_vector_ends(char *dest, char *src, int len, int index, int elem_sz)
00346 {
00347 memcpy(dest, src + index * elem_sz, (len - index) * elem_sz);
00348
00349 memcpy(dest + (len - index) * elem_sz, src, index * elem_sz);
00350 }
00351
00360 void GeoConstraint::reorder_longitude_map(int longitude_index_left)
00361 {
00362 double *tmp_lon = new double[d_lon_length];
00363
00364 swap_vector_ends((char *) tmp_lon, (char *) d_lon, d_lon_length,
00365 longitude_index_left, sizeof(double));
00366
00367 memcpy(d_lon, tmp_lon, d_lon_length * sizeof(double));
00368
00369 delete[]tmp_lon;
00370 }
00371
00372 static int
00373 count_dimensions_except_longitude(Array &a)
00374 {
00375 int size = 1;
00376 for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i)
00377 size *= a.dimension_size(i, true);
00378
00379 return size;
00380 }
00381
00400 void GeoConstraint::reorder_data_longitude_axis(Array &a)
00401 {
00402
00403 if (!get_longitude_rightmost())
00404 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.");
00405
00406 DBG(cerr << "Constraint for the left half: " << get_longitude_index_left()
00407 << ", " << get_lon_length() - 1 << endl);
00408
00409 a.add_constraint(d_lon_dim, get_longitude_index_left(), 1,
00410 get_lon_length() - 1);
00411 a.set_read_p(false);
00412 a.read();
00413 DBG2(a.print_val(stderr));
00414 #if 0
00415 char *left_data = 0;
00416 int left_size = a.buf2val((void **) & left_data);
00417 #endif
00418 char *left_data = (char*)a.value();
00419 int left_size = a.length();
00420
00421
00422
00423
00424
00425
00426 DBG(cerr << "Constraint for the right half: " << 0
00427 << ", " << get_longitude_index_right() << endl);
00428
00429 a.add_constraint(d_lon_dim, 0, 1, get_longitude_index_right());
00430 a.set_read_p(false);
00431 a.read();
00432 DBG2(a.print_val(stderr));
00433 #if 0
00434 char *right_data = 0;
00435 int right_size = a.buf2val((void **) & right_data);
00436 #endif
00437 char *right_data = (char*)a.value();
00438 int right_size = a.length();
00439
00440
00441 d_array_data_size = left_size + right_size;
00442 d_array_data = new char[d_array_data_size];
00443
00444
00445
00446
00447 int elem_width = a.var()->width();
00448 int left_elements = (get_lon_length() - get_longitude_index_left()) * elem_width;
00449 int right_elements = (get_longitude_index_right() + 1) * elem_width;
00450 int total_elements_per_row = left_elements + right_elements;
00451
00452
00453 int rows_to_copy = count_dimensions_except_longitude(a);
00454 for (int i = 0; i < rows_to_copy; ++i) {
00455 memcpy(d_array_data + (total_elements_per_row * i),
00456 left_data + (left_elements * i),
00457 left_elements);
00458 memcpy(d_array_data + left_elements + (total_elements_per_row * i),
00459 right_data + (right_elements * i),
00460 right_elements);
00461 }
00462
00463 delete[]left_data;
00464 delete[]right_data;
00465 }
00466
00467 #if 0
00468
00481 void GeoConstraint::set_bounding_box_longitude(double left, double right)
00482 {
00483
00484 d_longitude_notation = categorize_notation(left, right);
00485
00486
00487 if (d_longitude_notation == neg_pos)
00488 transform_constraint_to_pos_notation(left, right);
00489
00490
00491
00492
00493 Notation longitude_notation =
00494 categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
00495
00496 if (longitude_notation == neg_pos)
00497 transform_longitude_to_pos_notation();
00498
00499
00500 find_longitude_indeces(left, right, d_longitude_index_left,
00501 d_longitude_index_right);
00502 }
00503 #endif
00504
00510 GeoConstraint::GeoConstraint()
00511 : d_array_data(0), d_array_data_size(0),
00512 d_lat(0), d_lon(0),
00513 d_bounding_box_set(false),
00514 d_longitude_notation(unknown_notation),
00515 d_latitude_sense(unknown_sense)
00516 {
00517
00518 d_coards_lat_units.insert("degrees_north");
00519 d_coards_lat_units.insert("degree_north");
00520 d_coards_lat_units.insert("degree_N");
00521 d_coards_lat_units.insert("degrees_N");
00522
00523 d_coards_lon_units.insert("degrees_east");
00524 d_coards_lon_units.insert("degree_east");
00525 d_coards_lon_units.insert("degrees_E");
00526 d_coards_lon_units.insert("degree_E");
00527
00528 d_lat_names.insert("COADSY");
00529 d_lat_names.insert("lat");
00530 d_lat_names.insert("Lat");
00531 d_lat_names.insert("LAT");
00532
00533 d_lon_names.insert("COADSX");
00534 d_lon_names.insert("lon");
00535 d_lon_names.insert("Lon");
00536 d_lon_names.insert("LON");
00537 }
00538
00549 void GeoConstraint::set_bounding_box(double left, double top,
00550 double right, double bottom)
00551 {
00552
00553
00554
00555 if (d_bounding_box_set)
00556 throw
00557 InternalErr
00558 ("It is not possible to register more than one geographical constraint on a variable.");
00559
00560
00561 d_latitude_sense = categorize_latitude();
00562
00563
00564 d_longitude_notation = categorize_notation(left, right);
00565
00566
00567 if (d_longitude_notation == neg_pos)
00568 transform_constraint_to_pos_notation(left, right);
00569
00570
00571
00572
00573 Notation longitude_notation =
00574 categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
00575
00576 if (longitude_notation == neg_pos)
00577 transform_longitude_to_pos_notation();
00578
00579 if (!is_bounding_box_valid(left, top, right, bottom))
00580 throw Error("The bounding box does not intersect any data within this Grid or Array. The\ngeographical extent of these data are from latitude "
00581 + double_to_string(d_lat[0]) + " to "
00582 + double_to_string(d_lat[d_lat_length-1])
00583 + "\nand longitude " + double_to_string(d_lon[0])
00584 + " to " + double_to_string(d_lon[d_lon_length-1])
00585 + " while the bounding box provided was latitude "
00586 + double_to_string(top) + " to "
00587 + double_to_string(bottom) + "\nand longitude "
00588 + double_to_string(left) + " to "
00589 + double_to_string(right));
00590
00591
00592
00593
00594
00595 find_latitude_indeces(top, bottom, d_latitude_sense,
00596 d_latitude_index_top, d_latitude_index_bottom);
00597
00598
00599
00600 find_longitude_indeces(left, right, d_longitude_index_left,
00601 d_longitude_index_right);
00602
00603 DBG(cerr << "Bounding box (tlbr): " << d_latitude_index_top << ", "
00604 << d_longitude_index_left << ", "
00605 << d_latitude_index_bottom << ", "
00606 << d_longitude_index_right << endl);
00607
00608 d_bounding_box_set = true;
00609 }
00610
00611 }