Go to the documentation of this file.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