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 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 #ifdef WIN32
00041 #include <algorithm>
00042 #endif
00043
00044
00045
00046 #include "debug.h"
00047 #include "dods-datatypes.h"
00048 #include "GeoConstraint.h"
00049 #include "Float64.h"
00050
00051 #include "Error.h"
00052 #include "InternalErr.h"
00053 #include "ce_functions.h"
00054 #include "util.h"
00055
00056 using namespace std;
00057 using namespace libdap;
00058
00065 void remove_quotes(string & value)
00066 {
00067 if (!value.empty() && value[0] == '"'
00068 && value[value.length() - 1] == '"') {
00069 value.erase(0, 1);
00070 value.erase(value.length() - 1);
00071 }
00072
00073 return;
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
00289 void GeoConstraint::find_latitude_indeces(double top, double bottom,
00290 LatitudeSense sense,
00291 int &latitude_index_top,
00292 int &latitude_index_bottom) const
00293 {
00294
00295 int i = 0;
00296 if (sense == normal) {
00297 while (i < d_lat_length - 1 && d_lat[i] > top)
00298 ++i;
00299 if (d_lat[i] == top)
00300 latitude_index_top = i;
00301 else
00302 latitude_index_top = (i - 1) > 0 ? i - 1 : 0;
00303 }
00304 else {
00305 while (i < d_lat_length - 1 && d_lat[i] < top)
00306 ++i;
00307 latitude_index_top = i;
00308 }
00309
00310
00311
00312 i = d_lat_length - 1;
00313 if (sense == normal) {
00314 while (i > 0 && d_lat[i] < bottom)
00315 --i;
00316 if (d_lat[i] == bottom)
00317 latitude_index_bottom = i;
00318 else
00319 latitude_index_bottom =
00320 (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1;
00321 }
00322 else {
00323 while (i > 0 && d_lat[i] > bottom)
00324 --i;
00325 latitude_index_bottom = i;
00326 }
00327 }
00328
00332 GeoConstraint::LatitudeSense GeoConstraint::categorize_latitude() const
00333 {
00334 return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted;
00335 }
00336
00337 #if 0
00338
00345 void GeoConstraint::set_bounding_box_latitude(double top, double bottom)
00346 {
00347
00348 d_latitude_sense = categorize_latitude();
00349
00350
00351
00352
00353
00354 find_latitude_indeces(top, bottom, d_latitude_sense,
00355 d_latitude_index_top, d_latitude_index_bottom);
00356 }
00357 #endif
00358
00359
00360
00361 static void
00362 swap_vector_ends(char *dest, char *src, int len, int index, int elem_sz)
00363 {
00364 memcpy(dest, src + index * elem_sz, (len - index) * elem_sz);
00365
00366 memcpy(dest + (len - index) * elem_sz, src, index * elem_sz);
00367 }
00368
00377 void GeoConstraint::reorder_longitude_map(int longitude_index_left)
00378 {
00379 double *tmp_lon = new double[d_lon_length];
00380
00381 swap_vector_ends((char *) tmp_lon, (char *) d_lon, d_lon_length,
00382 longitude_index_left, sizeof(double));
00383
00384 memcpy(d_lon, tmp_lon, d_lon_length * sizeof(double));
00385
00386 delete[]tmp_lon;
00387 }
00388
00389 static int
00390 count_dimensions_except_longitude(Array &a)
00391 {
00392 int size = 1;
00393 for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i)
00394 size *= a.dimension_size(i, true);
00395
00396 return size;
00397 }
00398
00417 void GeoConstraint::reorder_data_longitude_axis(Array &a)
00418 {
00419
00420 if (!get_longitude_rightmost())
00421 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.");
00422
00423 DBG(cerr << "Constraint for the left half: " << get_longitude_index_left()
00424 << ", " << get_lon_length() - 1 << endl);
00425
00426 a.add_constraint(d_lon_dim, get_longitude_index_left(), 1,
00427 get_lon_length() - 1);
00428 a.set_read_p(false);
00429 a.read(get_dataset());
00430 DBG2(a.print_val(stderr));
00431 char *left_data = 0;
00432 int left_size = a.buf2val((void **) & left_data);
00433
00434
00435
00436
00437
00438
00439 DBG(cerr << "Constraint for the right half: " << 0
00440 << ", " << get_longitude_index_right() << endl);
00441
00442 a.add_constraint(d_lon_dim, 0, 1, get_longitude_index_right());
00443 a.set_read_p(false);
00444 a.read(get_dataset());
00445 DBG2(a.print_val(stderr));
00446 char *right_data = 0;
00447 int right_size = a.buf2val((void **) & right_data);
00448
00449
00450 d_array_data_size = left_size + right_size;
00451 d_array_data = new char[d_array_data_size];
00452
00453
00454
00455
00456 int elem_width = a.var()->width();
00457 int left_elements = (get_lon_length() - get_longitude_index_left()) * elem_width;
00458 int right_elements = (get_longitude_index_right() + 1) * elem_width;
00459 int total_elements_per_row = left_elements + right_elements;
00460
00461
00462 int rows_to_copy = count_dimensions_except_longitude(a);
00463 for (int i = 0; i < rows_to_copy; ++i) {
00464 memcpy(d_array_data + (total_elements_per_row * i),
00465 left_data + (left_elements * i),
00466 left_elements);
00467 memcpy(d_array_data + left_elements + (total_elements_per_row * i),
00468 right_data + (right_elements * i),
00469 right_elements);
00470 }
00471
00472 delete[]left_data;
00473 delete[]right_data;
00474 }
00475
00476 #if 0
00477
00490 void GeoConstraint::set_bounding_box_longitude(double left, double right)
00491 {
00492
00493 d_longitude_notation = categorize_notation(left, right);
00494
00495
00496 if (d_longitude_notation == neg_pos)
00497 transform_constraint_to_pos_notation(left, right);
00498
00499
00500
00501
00502 Notation longitude_notation =
00503 categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
00504
00505 if (longitude_notation == neg_pos)
00506 transform_longitude_to_pos_notation();
00507
00508
00509 find_longitude_indeces(left, right, d_longitude_index_left,
00510 d_longitude_index_right);
00511 }
00512 #endif
00513
00521 GeoConstraint::GeoConstraint(const string & ds_name)
00522 : d_dataset(ds_name), d_array_data(0), d_array_data_size(0),
00523 d_lat(0), d_lon(0),
00524 d_bounding_box_set(false),
00525 d_longitude_notation(unknown_notation),
00526 d_latitude_sense(unknown_sense)
00527 {
00528
00529 d_coards_lat_units.insert("degrees_north");
00530 d_coards_lat_units.insert("degree_north");
00531 d_coards_lat_units.insert("degree_N");
00532 d_coards_lat_units.insert("degrees_N");
00533
00534 d_coards_lon_units.insert("degrees_east");
00535 d_coards_lon_units.insert("degree_east");
00536 d_coards_lon_units.insert("degrees_E");
00537 d_coards_lon_units.insert("degree_E");
00538
00539 d_lat_names.insert("COADSY");
00540 d_lat_names.insert("lat");
00541 d_lat_names.insert("Lat");
00542 d_lat_names.insert("LAT");
00543
00544 d_lon_names.insert("COADSX");
00545 d_lon_names.insert("lon");
00546 d_lon_names.insert("Lon");
00547 d_lon_names.insert("LON");
00548 }
00549
00560 void GeoConstraint::set_bounding_box(double left, double top,
00561 double right, double bottom)
00562 {
00563
00564
00565
00566 if (d_bounding_box_set)
00567 throw
00568 InternalErr
00569 ("It is not possible to register more than one geographical constraint on a variable.");
00570
00571
00572 d_latitude_sense = categorize_latitude();
00573
00574
00575 d_longitude_notation = categorize_notation(left, right);
00576
00577
00578 if (d_longitude_notation == neg_pos)
00579 transform_constraint_to_pos_notation(left, right);
00580
00581
00582
00583
00584 Notation longitude_notation =
00585 categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
00586
00587 if (longitude_notation == neg_pos)
00588 transform_longitude_to_pos_notation();
00589
00590 if (!is_bounding_box_valid(left, top, right, bottom))
00591 throw Error("The bounding box does not intersect any data within this Grid or Array. The\ngeographical extent of these data are from latitude "
00592 + double_to_string(d_lat[0]) + " to "
00593 + double_to_string(d_lat[d_lat_length-1])
00594 + "\nand longitude " + double_to_string(d_lon[0])
00595 + " to " + double_to_string(d_lon[d_lon_length-1])
00596 + " while the bounding box provided was latitude "
00597 + double_to_string(top) + " to "
00598 + double_to_string(bottom) + "\nand longitude "
00599 + double_to_string(left) + " to "
00600 + double_to_string(right));
00601
00602
00603
00604
00605
00606 find_latitude_indeces(top, bottom, d_latitude_sense,
00607 d_latitude_index_top, d_latitude_index_bottom);
00608
00609
00610
00611 find_longitude_indeces(left, right, d_longitude_index_left,
00612 d_longitude_index_right);
00613
00614 DBG(cerr << "Bounding box (tlbr): " << d_latitude_index_top << ", "
00615 << d_longitude_index_left << ", "
00616 << d_latitude_index_bottom << ", "
00617 << d_longitude_index_right << endl);
00618
00619 d_bounding_box_set = true;
00620 }