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 21922 2010-01-07 18:23:57Z jimg $"
00033 };
00034
00035 #include <cstring>
00036 #include <cmath>
00037 #include <iostream>
00038 #include <sstream>
00039 #include <algorithm>
00040
00041
00042
00043
00044 #include "debug.h"
00045 #include "dods-datatypes.h"
00046 #include "GeoConstraint.h"
00047 #include "Float64.h"
00048
00049 #include "Error.h"
00050 #include "InternalErr.h"
00051 #include "ce_functions.h"
00052 #include "util.h"
00053
00054 using namespace std;
00055
00056 namespace libdap {
00057
00062 class is_prefix
00063 {
00064 private:
00065 string s;
00066 public:
00067 is_prefix(const string & in): s(in)
00068 {}
00069
00070 bool operator()(const string & prefix)
00071 {
00072 return s.find(prefix) == 0;
00073 }
00074 };
00075
00086 bool
00087 unit_or_name_match(set < string > units, set < string > names,
00088 const string & var_units, const string & var_name)
00089 {
00090 return (units.find(var_units) != units.end()
00091 || find_if(names.begin(), names.end(),
00092 is_prefix(var_name)) != names.end());
00093 }
00094
00109 GeoConstraint::Notation
00110 GeoConstraint::categorize_notation(const double left,
00111 const double right) const
00112 {
00113 return (left < 0 || right < 0) ? neg_pos : pos;
00114 }
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130 void
00131 GeoConstraint::transform_constraint_to_pos_notation(double &left,
00132 double &right) const
00133 {
00134 if (left < 0)
00135 left += 360 ;
00136
00137 if (right < 0)
00138 right += 360;
00139 }
00140
00149 void GeoConstraint::transform_longitude_to_pos_notation()
00150 {
00151
00152
00153
00154 for (int i = 0; i < d_lon_length; ++i)
00155 if (d_lon[i] < 0)
00156 d_lon[i] += 360;
00157 }
00158
00168 void GeoConstraint::transform_longitude_to_neg_pos_notation()
00169 {
00170 for (int i = 0; i < d_lon_length; ++i)
00171 if (d_lon[i] > 180)
00172 d_lon[i] -= 360;
00173 }
00174
00175 bool GeoConstraint::is_bounding_box_valid(const double left, const double top,
00176 const double right, const double bottom) const
00177 {
00178 if ((left < d_lon[0] && right < d_lon[0])
00179 || (left > d_lon[d_lon_length-1] && right > d_lon[d_lon_length-1]))
00180 return false;
00181
00182 if (d_latitude_sense == normal) {
00183
00184 if ((top > d_lat[0] && bottom > d_lat[0])
00185 || (top < d_lat[d_lat_length-1] && bottom < d_lat[d_lat_length-1]))
00186 return false;
00187 }
00188 else {
00189 if ((top < d_lat[0] && bottom < d_lat[0])
00190 || (top > d_lat[d_lat_length-1] && bottom > d_lat[d_lat_length-1]))
00191 return false;
00192 }
00193
00194 return true;
00195 }
00196
00207 void GeoConstraint::find_longitude_indeces(double left, double right,
00208 int &longitude_index_left,
00209 int &longitude_index_right) const
00210 {
00211
00212
00213
00214 double t_left = fmod(left, 360.0);
00215 double t_right = fmod(right, 360.0);
00216
00217
00218
00219
00220
00221 int i = 0;
00222 int lon_origin_index = 0;
00223 double smallest_lon = fmod(d_lon[0], 360.0);
00224 while (i < d_lon_length) {
00225 double curent_lon_value = fmod(d_lon[i], 360.0);
00226 if (smallest_lon > curent_lon_value) {
00227 smallest_lon = curent_lon_value;
00228 lon_origin_index = i;
00229 }
00230 ++i;
00231 }
00232
00233 DBG2(cerr << "lon_origin_index: " << lon_origin_index << endl);
00234
00235
00236
00237
00238 i = lon_origin_index;
00239 while (fmod(d_lon[i], 360.0) < t_left) {
00240 ++i;
00241 i = i % d_lon_length;
00242
00243
00244 if (i == lon_origin_index)
00245 throw Error("geogrid: Could not find an index for the longitude value '" + long_to_string(left) + "'");
00246 }
00247
00248 if (fmod(d_lon[i], 360.0) == t_left)
00249 longitude_index_left = i;
00250 else
00251 longitude_index_left = (i - 1) > 0 ? i - 1 : 0;
00252
00253 DBG2(cerr << "longitude_index_left: " << longitude_index_left << endl);
00254
00255
00256
00257 int largest_lon_index = (lon_origin_index - 1 + d_lon_length) % d_lon_length;
00258 i = largest_lon_index;
00259 while (fmod(d_lon[i], 360.0) > t_right) {
00260
00261 i = (i == 0) ? d_lon_length - 1 : i - 1;
00262 if (i == largest_lon_index)
00263 throw Error("geogrid: Could not find an index for the longitude value '" + long_to_string(right) + "'");
00264 }
00265
00266 if (fmod(d_lon[i], 360.0) == t_right)
00267 longitude_index_right = i;
00268 else
00269 longitude_index_right = (i + 1) < d_lon_length - 1 ? i + 1 : d_lon_length - 1;
00270
00271 DBG2(cerr << "longitude_index_right: " << longitude_index_right << endl);
00272 }
00273
00286 void GeoConstraint::find_latitude_indeces(double top, double bottom,
00287 LatitudeSense sense,
00288 int &latitude_index_top,
00289 int &latitude_index_bottom) const
00290 {
00291 int i, j;
00292
00293 if (sense == normal) {
00294 i = 0;
00295 while (i < d_lat_length - 1 && top < d_lat[i])
00296 ++i;
00297
00298 j = d_lat_length - 1;
00299 while (j > 0 && bottom > d_lat[j])
00300 --j;
00301
00302 if (d_lat[i] == top)
00303 latitude_index_top = i;
00304 else
00305 latitude_index_top = (i - 1) > 0 ? i - 1 : 0;
00306
00307 if (d_lat[j] == bottom)
00308 latitude_index_bottom = j;
00309 else
00310 latitude_index_bottom =
00311 (j + 1) < d_lat_length - 1 ? j + 1 : d_lat_length - 1;
00312 }
00313 else {
00314 i = d_lat_length - 1;
00315 while (i > 0 && d_lat[i] > top)
00316 --i;
00317
00318 j = 0;
00319 while (j < d_lat_length - 1 && d_lat[j] < bottom)
00320 ++j;
00321
00322 if (d_lat[i] == top)
00323 latitude_index_top = i;
00324 else
00325 latitude_index_top = (i + 1) < d_lat_length - 1 ? i + 1 : d_lat_length - 1;
00326
00327 if (d_lat[j] == bottom)
00328 latitude_index_bottom = j;
00329 else
00330 latitude_index_bottom = (j - 1) > 0 ? j - 1 : 0;
00331 }
00332 }
00333
00337 GeoConstraint::LatitudeSense GeoConstraint::categorize_latitude() const
00338 {
00339 return d_lat[0] >= d_lat[d_lat_length - 1] ? normal : inverted;
00340 }
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
00352 template<class T>
00353 static void transpose(std::vector<std::vector<T> > a,
00354 std::vector<std::vector<T> > b, int width, int height)
00355 {
00356 for (int i = 0; i < width; i++) {
00357 for (int j = 0; j < height; j++) {
00358 b[j][i] = a[i][j];
00359 }
00360 }
00361 }
00362
00370 void GeoConstraint::transpose_vector(double *src, const int length)
00371 {
00372 double *tmp = new double[length];
00373
00374 int i = 0, j = length-1;
00375 while (i < length)
00376 tmp[j--] = src[i++];
00377
00378 memcpy(src, tmp,length * sizeof(double));
00379
00380 delete[] tmp;
00381 }
00382
00383 static int
00384 count_size_except_latitude_and_longitude(Array &a)
00385 {
00386 if (a.dim_end() - a.dim_begin() <= 2)
00387 return 1;
00388
00389 int size = 1;
00390 for (Array::Dim_iter i = a.dim_begin(); (i + 2) != a.dim_end(); ++i)
00391 size *= a.dimension_size(i, true);
00392
00393 return size;
00394 }
00395
00396 void GeoConstraint::flip_latitude_within_array(Array &a, int lat_length,
00397 int lon_length)
00398 {
00399 if (!d_array_data) {
00400 a.read();
00401 d_array_data = static_cast<char*>(a.value());
00402 d_array_data_size = a.width();
00403 }
00404
00405 int size = count_size_except_latitude_and_longitude(a);
00406 char *tmp_data = new char[d_array_data_size];
00407 int array_elem_size = a.var()->width();
00408 int lat_lon_size = (d_array_data_size / size);
00409
00410 DBG(cerr << "lat, lon_length: " << lat_length << ", " << lon_length << endl);
00411 DBG(cerr << "size: " << size << endl);
00412 DBG(cerr << "d_array_data_size: " << d_array_data_size << endl);
00413 DBG(cerr << "array_elem_size: " << array_elem_size<< endl);
00414 DBG(cerr << "lat_lon_size: " << lat_lon_size<< endl);
00415
00416 for (int i = 0; i < size; ++i) {
00417 int lat = 0;
00418 int s_lat = lat_length - 1;
00419
00420 int lon_size = array_elem_size * lon_length;
00421 int offset = i * lat_lon_size;
00422 while (s_lat > -1)
00423 memcpy(tmp_data + offset + (lat++ * lon_size),
00424 d_array_data + offset + (s_lat-- * lon_size),
00425 lon_size);
00426 }
00427
00428 memcpy(d_array_data, tmp_data, d_array_data_size);
00429 delete [] tmp_data;
00430 }
00431
00440 void GeoConstraint::reorder_longitude_map(int longitude_index_left)
00441 {
00442 double *tmp_lon = new double[d_lon_length];
00443
00444 swap_vector_ends((char *) tmp_lon, (char *) d_lon, d_lon_length,
00445 longitude_index_left, sizeof(double));
00446
00447 memcpy(d_lon, tmp_lon, d_lon_length * sizeof(double));
00448
00449 delete[]tmp_lon;
00450 }
00451
00452 static int
00453 count_dimensions_except_longitude(Array &a)
00454 {
00455 int size = 1;
00456 for (Array::Dim_iter i = a.dim_begin(); (i + 1) != a.dim_end(); ++i)
00457 size *= a.dimension_size(i, true);
00458
00459 return size;
00460 }
00461
00479 void GeoConstraint::reorder_data_longitude_axis(Array &a, Array::Dim_iter lon_dim)
00480 {
00481
00482 if (!is_longitude_rightmost())
00483 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.");
00484
00485 DBG(cerr << "Constraint for the left half: " << get_longitude_index_left()
00486 << ", " << get_lon_length() - 1 << endl);
00487
00488
00489 a.add_constraint(lon_dim, get_longitude_index_left(), 1,
00490 get_lon_length() - 1);
00491 a.set_read_p(false);
00492 a.read();
00493 DBG2(a.print_val(stderr));
00494
00495
00496 int left_size = a.width();
00497 char *left_data = (char*)a.value();
00498
00499
00500
00501
00502
00503
00504
00505 DBG(cerr << "Constraint for the right half: " << 0
00506 << ", " << get_longitude_index_right() << endl);
00507
00508 a.add_constraint(lon_dim, 0, 1, get_longitude_index_right());
00509 a.set_read_p(false);
00510 a.read();
00511 DBG2(a.print_val(stderr));
00512
00513 char *right_data = (char*)a.value();
00514 int right_size = a.width();
00515
00516
00517 d_array_data_size = left_size + right_size;
00518 d_array_data = new char[d_array_data_size];
00519
00520
00521
00522
00523 int elem_size = a.var()->width();
00524 int left_row_size = (get_lon_length() - get_longitude_index_left()) * elem_size;
00525 int right_row_size = (get_longitude_index_right() + 1) * elem_size;
00526 int total_bytes_per_row = left_row_size + right_row_size;
00527
00528 DBG2(cerr << "elem_size: " << elem_size << "; left & right size: "
00529 << left_row_size << ", " << right_row_size << endl);
00530
00531
00532
00533 int rows_to_copy = count_dimensions_except_longitude(a);
00534 for (int i = 0; i < rows_to_copy; ++i) {
00535 DBG(cerr << "Copying " << i << "th row" << endl);
00536 DBG(cerr << "left memcpy: " << *(float *)(left_data + (left_row_size * i)) << endl);
00537
00538 memcpy(d_array_data + (total_bytes_per_row * i),
00539 left_data + (left_row_size * i),
00540 left_row_size);
00541
00542 DBG(cerr << "right memcpy: " << *(float *)(right_data + (right_row_size * i)) << endl);
00543
00544 memcpy(d_array_data + (total_bytes_per_row * i) + left_row_size,
00545 right_data + (right_row_size * i),
00546 right_row_size);
00547 }
00548
00549 delete[]left_data;
00550 delete[]right_data;
00551 }
00552
00555 GeoConstraint::GeoConstraint()
00556 : d_array_data(0), d_array_data_size(0),
00557 d_lat(0), d_lon(0),
00558 d_bounding_box_set(false),
00559 d_longitude_rightmost(false),
00560 d_longitude_notation(unknown_notation),
00561 d_latitude_sense(unknown_sense)
00562 {
00563
00564 d_coards_lat_units.insert("degrees_north");
00565 d_coards_lat_units.insert("degree_north");
00566 d_coards_lat_units.insert("degree_N");
00567 d_coards_lat_units.insert("degrees_N");
00568
00569 d_coards_lon_units.insert("degrees_east");
00570 d_coards_lon_units.insert("degree_east");
00571 d_coards_lon_units.insert("degrees_E");
00572 d_coards_lon_units.insert("degree_E");
00573
00574 d_lat_names.insert("COADSY");
00575 d_lat_names.insert("lat");
00576 d_lat_names.insert("Lat");
00577 d_lat_names.insert("LAT");
00578
00579 d_lon_names.insert("COADSX");
00580 d_lon_names.insert("lon");
00581 d_lon_names.insert("Lon");
00582 d_lon_names.insert("LON");
00583 }
00584
00595 void GeoConstraint::set_bounding_box(double top, double left,
00596 double bottom, double right)
00597 {
00598
00599
00600
00601 if (d_bounding_box_set)
00602 throw Error("It is not possible to register more than one geographical constraint on a variable.");
00603
00604
00605 d_latitude_sense = categorize_latitude();
00606 #if 0
00607 if (d_latitude_sense == inverted)
00608 throw Error("geogrid() does not currently work with inverted data (data where the north pole is at a negative latitude value).");
00609 #endif
00610
00611
00612 d_longitude_notation = categorize_notation(left, right);
00613
00614
00615 if (d_longitude_notation == neg_pos)
00616 transform_constraint_to_pos_notation(left, right);
00617
00618
00619
00620
00621 Notation longitude_notation =
00622 categorize_notation(d_lon[0], d_lon[d_lon_length - 1]);
00623
00624 if (longitude_notation == neg_pos)
00625 transform_longitude_to_pos_notation();
00626
00627 if (!is_bounding_box_valid(left, top, right, bottom))
00628 throw Error("The bounding box does not intersect any data within this Grid or Array. The\ngeographical extent of these data are from latitude "
00629 + double_to_string(d_lat[0]) + " to "
00630 + double_to_string(d_lat[d_lat_length-1])
00631 + "\nand longitude " + double_to_string(d_lon[0])
00632 + " to " + double_to_string(d_lon[d_lon_length-1])
00633 + " while the bounding box provided was latitude "
00634 + double_to_string(top) + " to "
00635 + double_to_string(bottom) + "\nand longitude "
00636 + double_to_string(left) + " to "
00637 + double_to_string(right));
00638
00639
00640
00641
00642
00643 find_latitude_indeces(top, bottom, d_latitude_sense,
00644 d_latitude_index_top, d_latitude_index_bottom);
00645
00646
00647
00648 find_longitude_indeces(left, right, d_longitude_index_left,
00649 d_longitude_index_right);
00650
00651 DBG(cerr << "Bounding box (tlbr): " << d_latitude_index_top << ", "
00652 << d_longitude_index_left << ", "
00653 << d_latitude_index_bottom << ", "
00654 << d_longitude_index_right << endl);
00655
00656 d_bounding_box_set = true;
00657 }
00658
00659 }