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 18315 2008-03-03 20:14:44Z 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 ArrayGeoConstraint::ArrayGeoConstraint(Array *array, const string &ds_name,
00064 double left, double top, double right, double bottom)
00065 : GeoConstraint(ds_name), d_array(array),
00066 d_extent(left, top, right, bottom), d_projection("plat-carre", "wgs84")
00067
00068 {
00069 m_init();
00070 }
00071
00072 ArrayGeoConstraint::ArrayGeoConstraint(Array *array, const string &ds_name,
00073 double left, double top, double right, double bottom,
00074 const string &projection, const string &datum)
00075 : GeoConstraint(ds_name), d_array(array),
00076 d_extent(left, top, right, bottom), d_projection(projection, datum)
00077
00078 {
00079 m_init();
00080 }
00081
00093 bool ArrayGeoConstraint::build_lat_lon_maps()
00094 {
00095
00096 set_longitude_rightmost(true);
00097 set_lon_dim(d_array->dim_begin() + (d_array->dimensions() - 1));
00098
00099 int number_elements_longitude = d_array->dimension_size(get_lon_dim());
00100 double *lon_map = new double[number_elements_longitude];
00101 for (int i = 0; i < number_elements_longitude; ++i) {
00102 lon_map[i] = ((d_extent.d_right - d_extent.d_left) / (number_elements_longitude - 1)) * i + d_extent.d_left;
00103 }
00104 set_lon(lon_map);
00105 set_lon_length(number_elements_longitude);
00106
00107
00108 set_lat_dim(d_array->dim_begin() + (d_array->dimensions() - 2));
00109
00110 int number_elements_latitude = d_array->dimension_size(get_lat_dim());
00111 double *lat_map = new double[number_elements_latitude];
00112 for (int i = 0; i < number_elements_latitude; ++i) {
00113 lat_map[i] = ((d_extent.d_bottom - d_extent.d_top) / (number_elements_latitude - 1)) * i + d_extent.d_top;
00114 }
00115 set_lat(lat_map);
00116 set_lat_length(number_elements_latitude);
00117
00118 return get_lat() && get_lon();
00119 }
00120
00128 bool
00129 ArrayGeoConstraint::lat_lon_dimensions_ok()
00130 {
00131 return true;
00132 }
00133
00148 void ArrayGeoConstraint::apply_constraint_to_data()
00149 {
00150 if (!get_bounding_box_set())
00151 throw InternalErr(
00152 "The Latitude and Longitude constraints must be set before calling\n\
00153 apply_constraint_to_data().");
00154
00155 if (get_latitude_sense() == inverted) {
00156 int tmp = get_latitude_index_top();
00157 set_latitude_index_top(get_latitude_index_bottom());
00158 set_latitude_index_bottom(tmp);
00159 }
00160
00161
00162
00163 if (get_latitude_index_top() > get_latitude_index_bottom())
00164 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.");
00165
00166 d_array->add_constraint(get_lat_dim(),
00167 get_latitude_index_top(), 1,
00168 get_latitude_index_bottom());
00169
00170
00171
00172 if (get_longitude_index_left() > get_longitude_index_right()) {
00173 reorder_data_longitude_axis(*d_array);
00174
00175
00176
00177
00178
00179 set_longitude_index_right(get_lon_length() - get_longitude_index_left()
00180 + get_longitude_index_right());
00181 set_longitude_index_left(0);
00182 }
00183
00184
00185
00186
00187
00188
00189 d_array->add_constraint(get_lon_dim(),
00190 get_longitude_index_left(), 1,
00191 get_longitude_index_right());
00192
00193
00194
00195 if (get_array_data()) {
00196 int size = d_array->val2buf(get_array_data());
00197 if (size != get_array_data_size())
00198 throw InternalErr
00199 ("Expected data size not copied to the Grid's buffer.");
00200 d_array->set_read_p(true);
00201 }
00202 else {
00203 d_array->read(get_dataset());
00204 }
00205 }
00206
00207 }
00208