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 16612 2007-06-04 22:08:57Z jimg $"
00033 };
00034
00035 #include <math.h>
00036 #include <string.h>
00037
00038 #include <iostream>
00039 #include <sstream>
00040
00041
00042
00043 #include "debug.h"
00044 #include "dods-datatypes.h"
00045 #include "ArrayGeoConstraint.h"
00046 #include "Float64.h"
00047
00048 #include "Error.h"
00049 #include "InternalErr.h"
00050 #include "ce_functions.h"
00051
00052 using namespace std;
00053 using namespace libdap;
00054
00055
00057 void ArrayGeoConstraint::m_init()
00058 {
00059 if (d_array->dimensions() < 2 || d_array->dimensions() > 3)
00060 throw Error("The geoarray() function works only with Arrays of two or three dimensions.");
00061
00062 build_lat_lon_maps();
00063 }
00064
00065 ArrayGeoConstraint::ArrayGeoConstraint(Array *array, const string &ds_name,
00066 double left, double top, double right, double bottom)
00067 : GeoConstraint(ds_name), d_array(array),
00068 d_extent(left, top, right, bottom), d_projection("plat-carre", "wgs84")
00069
00070 {
00071 m_init();
00072 }
00073
00074 ArrayGeoConstraint::ArrayGeoConstraint(Array *array, const string &ds_name,
00075 double left, double top, double right, double bottom,
00076 const string &projection, const string &datum)
00077 : GeoConstraint(ds_name), d_array(array),
00078 d_extent(left, top, right, bottom), 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 (!get_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);
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 int size = d_array->val2buf(get_array_data());
00199 if (size != get_array_data_size())
00200 throw InternalErr
00201 ("Expected data size not copied to the Grid's buffer.");
00202 d_array->set_read_p(true);
00203 }
00204 else {
00205 d_array->read(get_dataset());
00206 }
00207 }
00208