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: ArrayGeoConstraint.cc 21922 2010-01-07 18:23:57Z jimg $"
00033 };
00034
00035 #include <cmath>
00036 #include <iostream>
00037 #include <sstream>
00038
00039
00040
00041 #include "debug.h"
00042 #include "dods-datatypes.h"
00043 #include "ArrayGeoConstraint.h"
00044 #include "Float64.h"
00045
00046 #include "Error.h"
00047 #include "InternalErr.h"
00048 #include "ce_functions.h"
00049
00050 using namespace std;
00051
00052 namespace libdap {
00053
00055 void ArrayGeoConstraint::m_init()
00056 {
00057 if (d_array->dimensions() < 2 || d_array->dimensions() > 3)
00058 throw Error("The geoarray() function works only with Arrays of two or three dimensions.");
00059
00060 build_lat_lon_maps();
00061 }
00062
00063
00064
00065 ArrayGeoConstraint::ArrayGeoConstraint(Array *array, double top, double left,
00066 double bottom, double right)
00067 : GeoConstraint(), d_array(array),
00068 d_extent(top, left, bottom, right), d_projection("plat-carre", "wgs84")
00069
00070 {
00071 m_init();
00072 }
00073
00074 ArrayGeoConstraint::ArrayGeoConstraint(Array *array,
00075 double top, double left, double bottom, double right,
00076 const string &projection, const string &datum)
00077 : GeoConstraint(), d_array(array),
00078 d_extent(top, left, bottom, right), d_projection(projection, datum)
00079
00080 {
00081 m_init();
00082 }
00083
00095 bool ArrayGeoConstraint::build_lat_lon_maps()
00096 {
00097
00098 set_longitude_rightmost(true);
00099 set_lon_dim(d_array->dim_begin() + (d_array->dimensions() - 1));
00100
00101 int number_elements_longitude = d_array->dimension_size(get_lon_dim());
00102 double *lon_map = new double[number_elements_longitude];
00103 for (int i = 0; i < number_elements_longitude; ++i) {
00104 lon_map[i] = ((d_extent.d_right - d_extent.d_left) / (number_elements_longitude - 1)) * i + d_extent.d_left;
00105 }
00106 set_lon(lon_map);
00107 set_lon_length(number_elements_longitude);
00108
00109
00110 set_lat_dim(d_array->dim_begin() + (d_array->dimensions() - 2));
00111
00112 int number_elements_latitude = d_array->dimension_size(get_lat_dim());
00113 double *lat_map = new double[number_elements_latitude];
00114 for (int i = 0; i < number_elements_latitude; ++i) {
00115 lat_map[i] = ((d_extent.d_bottom - d_extent.d_top) / (number_elements_latitude - 1)) * i + d_extent.d_top;
00116 }
00117 set_lat(lat_map);
00118 set_lat_length(number_elements_latitude);
00119
00120 return get_lat() && get_lon();
00121 }
00122
00130 bool
00131 ArrayGeoConstraint::lat_lon_dimensions_ok()
00132 {
00133 return true;
00134 }
00135
00150 void ArrayGeoConstraint::apply_constraint_to_data()
00151 {
00152 if (!is_bounding_box_set())
00153 throw InternalErr(
00154 "The Latitude and Longitude constraints must be set before calling\n\
00155 apply_constraint_to_data().");
00156
00157 if (get_latitude_sense() == inverted) {
00158 int tmp = get_latitude_index_top();
00159 set_latitude_index_top(get_latitude_index_bottom());
00160 set_latitude_index_bottom(tmp);
00161 }
00162
00163
00164
00165 if (get_latitude_index_top() > get_latitude_index_bottom())
00166 throw Error("The upper and lower latitude indexes appear to be reversed. Please provide\nthe latitude bounding box numbers giving the northern-most latitude first.");
00167
00168 d_array->add_constraint(get_lat_dim(),
00169 get_latitude_index_top(), 1,
00170 get_latitude_index_bottom());
00171
00172
00173
00174 if (get_longitude_index_left() > get_longitude_index_right()) {
00175 reorder_data_longitude_axis(*d_array, get_lon_dim());
00176
00177
00178
00179
00180
00181 set_longitude_index_right(get_lon_length() - get_longitude_index_left()
00182 + get_longitude_index_right());
00183 set_longitude_index_left(0);
00184 }
00185
00186
00187
00188
00189
00190
00191 d_array->add_constraint(get_lon_dim(),
00192 get_longitude_index_left(), 1,
00193 get_longitude_index_right());
00194
00195
00196
00197 if (get_array_data()) {
00198
00199 int size = d_array->val2buf(get_array_data());
00200
00201 if (size != get_array_data_size())
00202 throw InternalErr
00203 ("Expected data size not copied to the Grid's buffer.");
00204 d_array->set_read_p(true);
00205 }
00206 else {
00207 d_array->read();
00208 }
00209 }
00210
00211 }
00212