libdap++ Updated for version 3.8.2
|
00001 // -*- mode: c++; c-basic-offset:4 -*- 00002 00003 // This file is part of libdap, A C++ implementation of the OPeNDAP Data 00004 // Access Protocol. 00005 00006 // Copyright (c) 2002,2003 OPeNDAP, Inc. 00007 // Author: James Gallagher <jgallagher@opendap.org> 00008 // 00009 // This library is free software; you can redistribute it and/or 00010 // modify it under the terms of the GNU Lesser General Public 00011 // License as published by the Free Software Foundation; either 00012 // version 2.1 of the License, or (at your option) any later version. 00013 // 00014 // This library is distributed in the hope that it will be useful, 00015 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 // Lesser General Public License for more details. 00018 // 00019 // You should have received a copy of the GNU Lesser General Public 00020 // License along with this library; if not, write to the Free Software 00021 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 // 00023 // You can contact OPeNDAP, Inc. at PO Box 112, Saunderstown, RI. 02874-0112. 00024 00025 // (c) COPYRIGHT URI/MIT 1999 00026 // Please read the full copyright statement in the file COPYRIGHT_URI. 00027 // 00028 // Authors: 00029 // jhrg,jimg James Gallagher <jgallagher@gso.uri.edu> 00030 00031 00032 // These functions are used by the CE evaluator 00033 // 00034 // 1/15/99 jhrg 00035 00036 #include "config.h" 00037 00038 static char rcsid[]not_used = 00039 { "$Id: ce_functions.cc 20518 2009-03-05 23:39:46Z jimg $" 00040 }; 00041 00042 #include <limits.h> 00043 00044 #include <cstdlib> // used by strtod() 00045 #include <cerrno> 00046 #include <cmath> 00047 #include <iostream> 00048 #include <vector> 00049 #include <algorithm> 00050 00051 // #define DODS_DEBUG 00052 00053 #include "BaseType.h" 00054 #include "Array.h" 00055 #include "Sequence.h" 00056 #include "Grid.h" 00057 #include "Error.h" 00058 #include "RValue.h" 00059 00060 #include "GSEClause.h" 00061 #include "GridGeoConstraint.h" 00062 #include "ArrayGeoConstraint.h" 00063 00064 #include "ce_functions.h" 00065 #include "gse_parser.h" 00066 #include "gse.tab.hh" 00067 #include "debug.h" 00068 #include "util.h" 00069 00070 // We wrapped VC++ 6.x strtod() to account for a short coming 00071 // in that function in regards to "NaN". I don't know if this 00072 // still applies in more recent versions of that product. 00073 // ROM - 12/2007 00074 #ifdef WIN32 00075 #include <limits> 00076 double w32strtod(const char *, char **); 00077 #endif 00078 00079 using namespace std; 00080 00081 int gse_parse(void *arg); 00082 void gse_restart(FILE * in); 00083 00084 // Glue routines declared in gse.lex 00085 void gse_switch_to_buffer(void *new_buffer); 00086 void gse_delete_buffer(void *buffer); 00087 void *gse_string(const char *yy_str); 00088 00089 namespace libdap { 00090 00092 inline bool double_eq(double lhs, double rhs, double epsilon = 1.0e-5) 00093 { 00094 if (lhs > rhs) 00095 return (lhs - rhs) < ((lhs + rhs) / epsilon); 00096 else 00097 return (rhs - lhs) < ((lhs + rhs) / epsilon); 00098 } 00099 00107 string extract_string_argument(BaseType * arg) 00108 { 00109 if (arg->type() != dods_str_c) 00110 throw Error(malformed_expr, 00111 "The function requires a DAP string argument."); 00112 00113 if (!arg->read_p()) 00114 throw InternalErr(__FILE__, __LINE__, 00115 "The CE Evaluator built an argument list where some constants held no values."); 00116 00117 string s = dynamic_cast<Str&>(*arg).value(); 00118 00119 DBG(cerr << "s: " << s << endl); 00120 00121 return s; 00122 } 00123 template<class T> static void set_array_using_double_helper(Array * a, 00124 double *src, int src_len) 00125 { 00126 T *values = new T[src_len]; 00127 for (int i = 0; i < src_len; ++i) 00128 values[i] = (T) src[i]; 00129 00130 #ifdef VAL2BUF 00131 a->val2buf(values, true); 00132 #else 00133 a->set_value(values, src_len); 00134 #endif 00135 00136 delete[]values; 00137 } 00138 00156 void set_array_using_double(Array * dest, double *src, int src_len) 00157 { 00158 // Simple types are Byte, ..., Float64, String and Url. 00159 if (dest->type() == dods_array_c && !dest->var()->is_simple_type() || dest->var()->type() == dods_str_c || dest->var()->type() == dods_url_c) 00160 throw InternalErr(__FILE__, __LINE__, 00161 "The function requires a DAP numeric-type array argument."); 00162 00163 // Test sizes. Note that Array::length() takes any constraint into account 00164 // when it returns the length. Even if this was removed, the 'helper' 00165 // function this uses calls Vector::val2buf() which uses Vector::width() 00166 // which in turn uses length(). 00167 if (dest->length() != src_len) 00168 throw InternalErr(__FILE__, __LINE__, 00169 "The source and destination array sizes don't match (" 00170 + long_to_string(src_len) + " versus " 00171 + long_to_string(dest->length()) + ")."); 00172 00173 // The types of arguments that the CE Parser will build for numeric 00174 // constants are limited to Uint32, Int32 and Float64. See ce_expr.y. 00175 // Expanded to work for any numeric type so it can be used for more than 00176 // just arguments. 00177 switch (dest->var()->type()) { 00178 case dods_byte_c: 00179 set_array_using_double_helper<dods_byte>(dest, src, src_len); 00180 break; 00181 case dods_uint16_c: 00182 set_array_using_double_helper<dods_uint16>(dest, src, src_len); 00183 break; 00184 case dods_int16_c: 00185 set_array_using_double_helper<dods_int16>(dest, src, src_len); 00186 break; 00187 case dods_uint32_c: 00188 set_array_using_double_helper<dods_uint32>(dest, src, src_len); 00189 break; 00190 case dods_int32_c: 00191 set_array_using_double_helper<dods_int32>(dest, src, src_len); 00192 break; 00193 case dods_float32_c: 00194 set_array_using_double_helper<dods_float32>(dest, src, src_len); 00195 break; 00196 case dods_float64_c: 00197 set_array_using_double_helper<dods_float64>(dest, src, src_len); 00198 break; 00199 default: 00200 throw InternalErr(__FILE__, __LINE__, 00201 "The argument list built by the CE parser contained an unsupported numeric type."); 00202 } 00203 00204 // Set the read_p property. 00205 dest->set_read_p(true); 00206 } 00207 00208 template<class T> static double *extract_double_array_helper(Array * a) 00209 { 00210 int length = a->length(); 00211 00212 T *b = new T[length]; 00213 a->value(b); 00214 00215 double *dest = new double[length]; 00216 for (int i = 0; i < length; ++i) 00217 dest[i] = (double) b[i]; 00218 delete[]b; 00219 00220 return dest; 00221 } 00222 00227 double *extract_double_array(Array * a) 00228 { 00229 // Simple types are Byte, ..., Float64, String and Url. 00230 if (a->type() == dods_array_c && !a->var()->is_simple_type() || a->var()->type() == dods_str_c || a->var()->type() == dods_url_c) 00231 throw Error(malformed_expr, 00232 "The function requires a DAP numeric-type array argument."); 00233 00234 if (!a->read_p()) 00235 throw InternalErr(__FILE__, __LINE__, 00236 string("The Array '") + a->name() + 00237 "'does not contain values."); 00238 00239 // The types of arguments that the CE Parser will build for numeric 00240 // constants are limited to Uint32, Int32 and Float64. See ce_expr.y. 00241 // Expanded to work for any numeric type so it can be used for more than 00242 // just arguments. 00243 switch (a->var()->type()) { 00244 case dods_byte_c: 00245 return extract_double_array_helper<dods_byte>(a); 00246 case dods_uint16_c: 00247 return extract_double_array_helper<dods_uint16>(a); 00248 case dods_int16_c: 00249 return extract_double_array_helper<dods_int16>(a); 00250 case dods_uint32_c: 00251 return extract_double_array_helper<dods_uint32>(a); 00252 case dods_int32_c: 00253 return extract_double_array_helper<dods_int32>(a); 00254 case dods_float32_c: 00255 return extract_double_array_helper<dods_float32>(a); 00256 case dods_float64_c: 00257 return extract_double_array_helper<dods_float64>(a); 00258 default: 00259 throw InternalErr(__FILE__, __LINE__, 00260 "The argument list built by the CE parser contained an unsupported numeric type."); 00261 } 00262 } 00263 00271 double extract_double_value(BaseType * arg) 00272 { 00273 // Simple types are Byte, ..., Float64, String and Url. 00274 if (!arg->is_simple_type() || arg->type() == dods_str_c || arg->type() 00275 == dods_url_c) 00276 throw Error(malformed_expr, 00277 "The function requires a DAP numeric-type argument."); 00278 00279 if (!arg->read_p()) 00280 throw InternalErr(__FILE__, __LINE__, 00281 "The CE Evaluator built an argument list where some constants held no values."); 00282 00283 // The types of arguments that the CE Parser will build for numeric 00284 // constants are limited to Uint32, Int32 and Float64. See ce_expr.y. 00285 // Expanded to work for any numeric type so it can be used for more than 00286 // just arguments. 00287 switch (arg->type()) { 00288 case dods_byte_c: 00289 return (double)(dynamic_cast<Byte&>(*arg).value()); 00290 case dods_uint16_c: 00291 return (double)(dynamic_cast<UInt16&>(*arg).value()); 00292 case dods_int16_c: 00293 return (double)(dynamic_cast<Int16&>(*arg).value()); 00294 case dods_uint32_c: 00295 return (double)(dynamic_cast<UInt32&>(*arg).value()); 00296 case dods_int32_c: 00297 return (double)(dynamic_cast<Int32&>(*arg).value()); 00298 case dods_float32_c: 00299 return (double)(dynamic_cast<Float32&>(*arg).value()); 00300 case dods_float64_c: 00301 return dynamic_cast<Float64&>(*arg).value(); 00302 default: 00303 throw InternalErr(__FILE__, __LINE__, 00304 "The argument list built by the CE parser contained an unsupported numeric type."); 00305 } 00306 } 00307 00310 void 00311 function_version(int, BaseType *[], DDS &, BaseType **btpp) 00312 { 00313 string 00314 xml_value = 00315 "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\ 00316 <functions>\ 00317 <function name=\"version\" version=\"1.0\"/>\ 00318 <function name=\"grid\" version=\"1.0\"/>\ 00319 <function name=\"geogrid\" version=\"1.0b2\"/>\ 00320 <function name=\"geoarray\" version=\"0.9b1\"/>\ 00321 <function name=\"linear_scale\" version=\"1.0b1\"/>\ 00322 </functions>"; 00323 00324 Str *response = new Str("version"); 00325 00326 response->set_value(xml_value); 00327 *btpp = response; 00328 return; 00329 } 00330 00331 static void parse_gse_expression(gse_arg * arg, BaseType * expr) 00332 { 00333 gse_restart(0); // Restart the scanner. 00334 void *cls = gse_string(extract_string_argument(expr).c_str()); 00335 // gse_switch_to_buffer(cls); // Get set to scan the string. 00336 bool status = gse_parse((void *) arg) == 0; 00337 gse_delete_buffer(cls); 00338 if (!status) 00339 throw Error(malformed_expr, "Error parsing grid selection."); 00340 } 00341 00342 static void apply_grid_selection_expr(Grid * grid, GSEClause * clause) 00343 { 00344 // Basic plan: For each map, look at each clause and set start and stop 00345 // to be the intersection of the ranges in those clauses. 00346 Grid::Map_iter map_i = grid->map_begin(); 00347 while (map_i != grid->map_end() && (*map_i)->name() != clause->get_map_name()) 00348 ++map_i; 00349 00350 if (map_i == grid->map_end()) 00351 throw Error(malformed_expr,"The map vector '" + clause->get_map_name() 00352 + "' is not in the grid '" + grid->name() + "'."); 00353 00354 // Use pointer arith & the rule that map order must match array dim order 00355 Array::Dim_iter grid_dim = (grid->get_array()->dim_begin() + (map_i - grid->map_begin())); 00356 00357 Array *map = dynamic_cast < Array * >((*map_i)); 00358 if (!map) 00359 throw InternalErr(__FILE__, __LINE__, "Expected an Array"); 00360 int start = max(map->dimension_start(map->dim_begin()), clause->get_start()); 00361 int stop = min(map->dimension_stop(map->dim_begin()), clause->get_stop()); 00362 00363 if (start > stop) { 00364 ostringstream msg; 00365 msg 00366 << "The expressions passed to grid() do not result in an inclusive \n" 00367 << "subset of '" << clause->get_map_name() 00368 << "'. The map's values range " << "from " 00369 << clause->get_map_min_value() << " to " 00370 << clause->get_map_max_value() << "."; 00371 throw Error(malformed_expr,msg.str()); 00372 } 00373 00374 DBG(cerr << "Setting constraint on " << map->name() 00375 << "[" << start << ":" << stop << "]" << endl); 00376 00377 // Stride is always one. 00378 map->add_constraint(map->dim_begin(), start, 1, stop); 00379 grid->get_array()->add_constraint(grid_dim, start, 1, stop); 00380 } 00381 00382 static void apply_grid_selection_expressions(Grid * grid, 00383 vector < GSEClause * >clauses) 00384 { 00385 vector < GSEClause * >::iterator clause_i = clauses.begin(); 00386 while (clause_i != clauses.end()) 00387 apply_grid_selection_expr(grid, *clause_i++); 00388 00389 grid->set_read_p(false); 00390 } 00391 00428 void 00429 function_grid(int argc, BaseType * argv[], DDS &, BaseType **btpp) 00430 { 00431 DBG(cerr << "Entering function_grid..." << endl); 00432 00433 string 00434 info = 00435 "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\ 00436 <function name=\"grid\" version=\"1.0\"/>\ 00437 The grid() function takes a grid variable and zero or more relational\ 00438 expressions. Each relational expression is applied to the grid using the\ 00439 server's constraint evaluator and the resulting grid is returned. The\ 00440 expressions may use constants and the grid's map vectors but may not use\ 00441 any other variables. Two forms of expression are provide: \"var relop const\"\ 00442 and \"const relop var relop const\". For example: grid(sst, \"10<=TIME<20\")\ 00443 and grid(sst, \"10<=TIME\", \"TIME<20\") are both legal and, in this case,\ 00444 also equivalent.\ 00445 </function>"; 00446 00447 if (argc == 0) { 00448 Str *response = new Str("info"); 00449 response->set_value(info); 00450 *btpp = response; 00451 return; 00452 } 00453 00454 Grid *original_grid = dynamic_cast < Grid * >(argv[0]); 00455 if (!original_grid) 00456 throw Error(malformed_expr,"The first argument to grid() must be a Grid variable!"); 00457 00458 // Duplicate the grid; DODSFilter::send_data() will delete the variable 00459 // after serializing it. 00460 Grid *l_grid = dynamic_cast < Grid * >(original_grid->ptr_duplicate()); 00461 if (!l_grid) 00462 throw InternalErr(__FILE__, __LINE__, "Expected a Grid."); 00463 00464 DBG(cerr << "grid: past initialization code" << endl); 00465 00466 // Read the maps. Do this before calling parse_gse_expression(). Avoid 00467 // reading the array until the constraints have been applied because it 00468 // might be really large. 00469 00470 // This version makes sure to set the send_p flags which is needed for 00471 // the hdf4 handler (and is what should be done in general). 00472 Grid::Map_iter i = l_grid->map_begin(); 00473 while (i != l_grid->map_end()) 00474 (*i++)->set_send_p(true); 00475 l_grid->read(); 00476 00477 DBG(cerr << "grid: past map read" << endl); 00478 00479 // argv[1..n] holds strings; each are little expressions to be parsed. 00480 // When each expression is parsed, the parser makes a new instance of 00481 // GSEClause. GSEClause checks to make sure the named map really exists 00482 // in the Grid and that the range of values given makes sense. 00483 vector < GSEClause * > clauses; 00484 gse_arg *arg = new gse_arg(l_grid); 00485 for (int i = 1; i < argc; ++i) { 00486 parse_gse_expression(arg, argv[i]); 00487 clauses.push_back(arg->get_gsec()); 00488 } 00489 delete arg; 00490 arg = 0; 00491 00492 apply_grid_selection_expressions(l_grid, clauses); 00493 00494 DBG(cerr << "grid: past gse application" << endl); 00495 00496 l_grid->get_array()->set_send_p(true); 00497 00498 l_grid->read(); 00499 00500 *btpp = l_grid; 00501 return; 00502 } 00503 00539 void 00540 function_geogrid(int argc, BaseType * argv[], DDS &, BaseType **btpp) 00541 { 00542 string 00543 info = 00544 "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\ 00545 <function name=\"geogrid\" version=\"1.0b2\"/>\ 00546 geogrid() applies a constraint given in latitude and longitude to a\ 00547 DAP Grid variable. The arguments to the function are:\ 00548 geogrid(<grid variable>, <upper latitude>, <left longitude>,\ 00549 <lower latitude>, <right longitude> [selection expressions - see grid()])\ 00550 geogrid() returns the version of the function.\ 00551 The function will always return a single Grid variable whose values\ 00552 completely cover the given region, although there may be cases when\ 00553 some additional data is also returned. If the longitude values 'wrap\ 00554 around' the right edge of the data, then the function will make two\ 00555 requests and return those joined together as a single Grid.\ 00556 </function>"; 00557 00558 if (argc == 0) { 00559 Str *response = new Str("version"); 00560 response->set_value(info); 00561 *btpp = response; 00562 return ; 00563 } 00564 00565 if (argc < 5) 00566 throw Error(malformed_expr,"Wrong number of arguments to geogrid(). See geogrid() for more information."); 00567 00568 Grid *l_grid = dynamic_cast < Grid * >(argv[0]->ptr_duplicate()); 00569 if (!l_grid) 00570 throw Error(malformed_expr,"The first argument to geogrid() must be a Grid variable!"); 00571 00572 // Read the maps. Do this before calling parse_gse_expression(). Avoid 00573 // reading the array until the constraints have been applied because it 00574 // might be really large. 00575 // 00576 // Trick: Some handlers build Grids from a combination of Array 00577 // variables and attributes. Those handlers (e.g., hdf4) use the send_p 00578 // property to determine which parts of the Grid to read *but they can 00579 // only read the maps from within Grid::read(), not the map's read()*. 00580 // Since the Grid's array does not have send_p set, it will not be read 00581 // by the call below to Grid::read(). 00582 Grid::Map_iter i = l_grid->map_begin(); 00583 while (i != l_grid->map_end()) 00584 (*i++)->set_send_p(true); 00585 l_grid->read(); 00586 // Calling read() above sets the read_p flag for the entire grid; clear it 00587 // for the grid's array so that later on the code will be sure to read it 00588 // under all circumstances. 00589 l_grid->get_array()->set_read_p(false);DBG(cerr << "geogrid: past map read" << endl); 00590 00591 // Look for Grid Selection Expressions tacked onto the end of the BB 00592 // specification. If there are any, evaluate them before evaluating the BB. 00593 if (argc > 5) { 00594 // argv[5..n] holds strings; each are little Grid Selection Expressions 00595 // to be parsed and evaluated. 00596 vector < GSEClause * > clauses; 00597 gse_arg *arg = new gse_arg(l_grid); 00598 for (int i = 5; i < argc; ++i) { 00599 parse_gse_expression(arg, argv[i]); 00600 clauses.push_back(arg->get_gsec()); 00601 } 00602 delete arg; 00603 arg = 0; 00604 00605 apply_grid_selection_expressions(l_grid, clauses); 00606 } 00607 00608 try { 00609 // Build a GeoConstraint object. If there are no longitude/latitude 00610 // maps then this constructor throws Error. 00611 GridGeoConstraint gc(l_grid); 00612 00613 // This sets the bounding box and modifies the maps to match the 00614 // notation of the box (0/359 or -180/179) 00615 double top = extract_double_value(argv[1]); 00616 double left = extract_double_value(argv[2]); 00617 double bottom = extract_double_value(argv[3]); 00618 double right = extract_double_value(argv[4]); 00619 gc.set_bounding_box(left, top, right, bottom); 00620 DBG(cerr << "geogrid: past bounding box set" << endl); 00621 00622 // This also reads all of the data into the grid variable 00623 gc.apply_constraint_to_data(); 00624 DBG(cerr << "geogrid: past apply constraint" << endl); 00625 00626 // In this function the l_grid pointer is the same as the pointer returned 00627 // by this call. The caller of the function must free the pointer. 00628 *btpp = gc.get_constrained_grid(); 00629 return; 00630 } 00631 catch (Error &e) { 00632 throw e; 00633 } 00634 catch (exception & e) { 00635 throw 00636 InternalErr(string 00637 ("A C++ exception was thrown from inside geogrid(): ") 00638 + e.what()); 00639 } 00640 } 00641 00642 // These static functions could be moved to a class that provides a more 00643 // general interface for COARDS/CF someday. Assume each BaseType comes bundled 00644 // with an attribute table. 00645 00646 // This was ripped from parser-util.cc 00647 static double string_to_double(const char *val) 00648 { 00649 char *ptr; 00650 errno = 0; 00651 // Clear previous value. 5/21/2001 jhrg 00652 00653 #ifdef WIN32 00654 double v = w32strtod(val, &ptr); 00655 #else 00656 double v = strtod(val, &ptr); 00657 #endif 00658 00659 if ((v == 0.0 && (val == ptr || errno == HUGE_VAL || errno == ERANGE)) 00660 || *ptr != '\0') { 00661 throw Error(malformed_expr,string("Could not convert the string '") + val + "' to a double."); 00662 } 00663 00664 double abs_val = fabs(v); 00665 if (abs_val > DODS_DBL_MAX || (abs_val != 0.0 && abs_val < DODS_DBL_MIN)) 00666 throw Error(malformed_expr,string("Could not convert the string '") + val + "' to a double."); 00667 00668 return v; 00669 } 00670 00674 static double get_attribute_double_value(BaseType *var, 00675 vector<string> &attributes) 00676 { 00677 AttrTable &attr = var->get_attr_table(); 00678 string attribute_value = ""; 00679 string values = ""; 00680 vector<string>::iterator i = attributes.begin(); 00681 while (attribute_value == "" && i != attributes.end()) { 00682 values += *i; 00683 if (!values.empty()) 00684 values += ", "; 00685 attribute_value = attr.get_attr(*i++); 00686 } 00687 00688 // If the value string is empty, then look at the grid's array (if it's a 00689 // grid or throw an Error. 00690 if (attribute_value.empty()) { 00691 if (var->type() == dods_grid_c) 00692 return get_attribute_double_value(dynamic_cast<Grid&>(*var).get_array(), attributes); 00693 else 00694 throw Error(malformed_expr,string("No COARDS '") + values.substr(0, values.length() - 2) 00695 + "' attribute was found for the variable '" 00696 + var->name() + "'."); 00697 } 00698 00699 return string_to_double(remove_quotes(attribute_value).c_str()); 00700 } 00701 00702 static double get_attribute_double_value(BaseType *var, const string &attribute) 00703 { 00704 AttrTable &attr = var->get_attr_table(); 00705 string attribute_value = attr.get_attr(attribute); 00706 00707 // If the value string is empty, then look at the grid's array (if it's a 00708 // grid or throw an Error. 00709 if (attribute_value.empty()) { 00710 if (var->type() == dods_grid_c) 00711 return get_attribute_double_value(dynamic_cast<Grid&>(*var).get_array(), attribute); 00712 else 00713 throw Error(malformed_expr,string("No COARDS '") + attribute 00714 + "' attribute was found for the variable '" 00715 + var->name() + "'."); 00716 } 00717 00718 return string_to_double(remove_quotes(attribute_value).c_str()); 00719 } 00720 00721 static double get_y_intercept(BaseType *var) 00722 { 00723 vector<string> attributes; 00724 attributes.push_back("add_offset"); 00725 attributes.push_back("add_off"); 00726 return get_attribute_double_value(var, attributes); 00727 } 00728 00729 static double get_slope(BaseType *var) 00730 { 00731 return get_attribute_double_value(var, "scale_factor"); 00732 } 00733 00734 static double get_missing_value(BaseType *var) 00735 { 00736 return get_attribute_double_value(var, "missing_value"); 00737 } 00738 00751 void 00752 function_linear_scale(int argc, BaseType * argv[], DDS &, BaseType **btpp) 00753 { 00754 string 00755 info = 00756 "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\ 00757 <function name=\"linear_scale\" version=\"1.0b1\">\ 00758 The linear_scale() function applies the familiar y=mx+b equation to data.\ 00759 If only the name of a variable is given, the function looks for the COARDS\ 00760 scale_factor, add_offset and missing_value attributes. In the equation,\ 00761 'm' is scale_factor, 'b' is add_offset and data values which match\ 00762 missing_value are not scaled. If add_offset cannot be found, it defaults to\ 00763 zero; if missing_value cannot be found, the test for it is not performed.\ 00764 The caller can also provide values for these and, if provided, those values\ 00765 are used and the dataset's attributes are ignored. If only scale_factor is\ 00766 provided, the defaults for 'b' and the missing value flag are used.\ 00767 Similarly, if only 'm' and 'b' are provided the missing value test is not\ 00768 performed. The values are always returned in a Float64 array.\ 00769 </function>"; 00770 00771 if (argc == 0) { 00772 Str *response = new Str("info"); 00773 response->set_value(info); 00774 *btpp = response; 00775 return; 00776 } 00777 00778 // Check for 1 or 3 arguments: 1 --> use attributes; 3 --> m & b supplied 00779 DBG(cerr << "argc = " << argc << endl); 00780 if (!(argc == 1 || argc == 3 || argc == 4)) 00781 throw Error(malformed_expr,"Wrong number of arguments to linear_scale(). See linear_scale() for more information"); 00782 00783 // Get m & b 00784 bool use_missing = false; 00785 double m, b, missing = 0.0; 00786 if (argc == 4) { 00787 m = extract_double_value(argv[1]); 00788 b = extract_double_value(argv[2]); 00789 missing = extract_double_value(argv[3]); 00790 use_missing = true; 00791 } else if (argc == 3) { 00792 m = extract_double_value(argv[1]); 00793 b = extract_double_value(argv[2]); 00794 use_missing = false; 00795 } else { 00796 m = get_slope(argv[0]); 00797 00798 // This is really a hack; on a fair number of datasets, the y intercept 00799 // is not given and is assumed to be 0. Here the function looks and 00800 // catches the error if a y intercept is not found. 00801 try { 00802 b = get_y_intercept(argv[0]); 00803 } 00804 catch (Error &e) { 00805 b = 0.0; 00806 } 00807 00808 // This is not the best plan; the get_missing_value() function should 00809 // do something other than throw, but to do that would require mayor 00810 // surgery on get_attribute_double_value(). 00811 try { 00812 missing = get_missing_value(argv[0]); 00813 use_missing = true; 00814 } 00815 catch (Error &e) { 00816 use_missing = false; 00817 } 00818 } 00819 00820 DBG(cerr << "m: " << m << ", b: " << b << endl);DBG(cerr << "use_missing: " << use_missing << ", missing: " << missing << endl); 00821 00822 // Read the data, scale and return the result. Must replace the new data 00823 // in a constructor (i.e., Array part of a Grid). 00824 BaseType *dest = 0; 00825 double *data; 00826 if (argv[0]->type() == dods_grid_c) { 00827 Array &source = *dynamic_cast<Grid&>(*argv[0]).get_array(); 00828 source.set_send_p(true); 00829 source.read(); 00830 data = extract_double_array(&source); 00831 int length = source.length(); 00832 int i = 0; 00833 while (i < length) { 00834 DBG2(cerr << "data[" << i << "]: " << data[i] << endl); 00835 if (!use_missing || !double_eq(data[i], missing)) 00836 data[i] = data[i] * m + b; 00837 DBG2(cerr << " >> data[" << i << "]: " << data[i] << endl); 00838 ++i; 00839 } 00840 00841 // Vector::add_var will delete the existing 'template' variable 00842 Float64 *temp_f = new Float64(source.name()); 00843 source.add_var(temp_f); 00844 #ifdef VAL2BUF 00845 source.val2buf(static_cast<void*>(data), false); 00846 #else 00847 source.set_value(data, i); 00848 #endif 00849 delete [] data; // val2buf copies. 00850 delete temp_f; // add_var copies and then adds. 00851 dest = argv[0]; 00852 } else if (argv[0]->is_vector_type()) { 00853 Array &source = dynamic_cast<Array&>(*argv[0]); 00854 source.set_send_p(true); 00855 // If the array is really a map, make sure to read using the Grid 00856 // because of the HDF4 handler's odd behavior WRT dimensions. 00857 if (source.get_parent() && source.get_parent()->type() == dods_grid_c) 00858 source.get_parent()->read(); 00859 else 00860 source.read(); 00861 00862 data = extract_double_array(&source); 00863 int length = source.length(); 00864 int i = 0; 00865 while (i < length) { 00866 if (!use_missing || !double_eq(data[i], missing)) 00867 data[i] = data[i] * m + b; 00868 ++i; 00869 } 00870 00871 Float64 *temp_f = new Float64(source.name()); 00872 source.add_var(temp_f); 00873 00874 source.val2buf(static_cast<void*>(data), false); 00875 00876 delete [] data; // val2buf copies. 00877 delete temp_f; // add_var copies and then adds. 00878 00879 dest = argv[0]; 00880 } else if (argv[0]->is_simple_type() && !(argv[0]->type() == dods_str_c 00881 || argv[0]->type() == dods_url_c)) { 00882 double data = extract_double_value(argv[0]); 00883 if (!use_missing || !double_eq(data, missing)) 00884 data = data * m + b; 00885 00886 dest = new Float64(argv[0]->name()); 00887 00888 dest->val2buf(static_cast<void*>(&data)); 00889 00890 } else { 00891 throw Error(malformed_expr,"The linear_scale() function works only for numeric Grids, Arrays and scalars."); 00892 } 00893 00894 *btpp = dest; 00895 return; 00896 } 00897 00914 void 00915 function_geoarray(int argc, BaseType * argv[], DDS &, BaseType **btpp) 00916 { 00917 string 00918 info = 00919 "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\ 00920 <function name=\"geoarray\" version=\"0.9b1\"/>\ 00921 The geoarray() function supports two different sets of arguments:\ 00922 geoarray(var,left,top,right,bottom)\ 00923 geoarray(var,left,top,right,bottom,var_left,var_top,var_right,var_bottom)\ 00924 In the first version 'var' is the target of the selection and 'left', 'top',\ 00925 'right' and 'bottom' are the corners of a longitude-latitude box that defines\ 00926 the selection. In the second version the 'var_left', ..., parameters give the\ 00927 longitude and latitude extent of the entire array. The projection and datum are\ 00928 assumed to be Plat-Carre and WGS84.\ 00929 </function>"; 00930 00931 if (argc == 0) { 00932 Str *response = new Str("version"); 00933 response->set_value(info); 00934 *btpp = response; 00935 return; 00936 } 00937 00938 DBG(cerr << "argc = " << argc << endl); 00939 if (!(argc == 5 || argc == 9 || argc == 11)) 00940 throw Error(malformed_expr,"Wrong number of arguments to geoarray(). See geoarray() for more information."); 00941 00942 // Check the Array (and dup because the caller will free the variable). 00943 Array *l_array = dynamic_cast < Array * >(argv[0]->ptr_duplicate()); 00944 if (!l_array) 00945 throw Error(malformed_expr,"The first argument to geoarray() must be an Array variable!"); 00946 00947 try { 00948 00949 // Read the bounding box and variable extents from the params 00950 double bb_top = extract_double_value(argv[1]); 00951 double bb_left = extract_double_value(argv[2]); 00952 double bb_bottom = extract_double_value(argv[3]); 00953 double bb_right = extract_double_value(argv[4]); 00954 00955 switch (argc) { 00956 case 5: { 00957 ArrayGeoConstraint agc(l_array); 00958 00959 agc.set_bounding_box(bb_left, bb_top, bb_right, bb_bottom); 00960 // This also reads all of the data into the grid variable 00961 agc.apply_constraint_to_data(); 00962 *btpp = agc.get_constrained_array(); 00963 return; 00964 break; 00965 } 00966 case 9: { 00967 double var_top = extract_double_value(argv[5]); 00968 double var_left = extract_double_value(argv[6]); 00969 double var_bottom = extract_double_value(argv[7]); 00970 double var_right = extract_double_value(argv[8]); 00971 ArrayGeoConstraint agc (l_array, var_left, var_top, var_right, var_bottom); 00972 00973 agc.set_bounding_box(bb_left, bb_top, bb_right, bb_bottom); 00974 // This also reads all of the data into the grid variable 00975 agc.apply_constraint_to_data(); 00976 *btpp = agc.get_constrained_array(); 00977 return; 00978 break; 00979 } 00980 case 11: { 00981 double var_top = extract_double_value(argv[5]); 00982 double var_left = extract_double_value(argv[6]); 00983 double var_bottom = extract_double_value(argv[7]); 00984 double var_right = extract_double_value(argv[8]); 00985 string projection = extract_string_argument(argv[9]); 00986 string datum = extract_string_argument(argv[10]); 00987 ArrayGeoConstraint agc(l_array, 00988 var_left, var_top, var_right, var_bottom, 00989 projection, datum); 00990 00991 agc.set_bounding_box(bb_left, bb_top, bb_right, bb_bottom); 00992 // This also reads all of the data into the grid variable 00993 agc.apply_constraint_to_data(); 00994 *btpp = agc.get_constrained_array(); 00995 return; 00996 break; 00997 } 00998 default: 00999 throw InternalErr(__FILE__, __LINE__, "Wrong number of args to geoarray."); 01000 } 01001 } 01002 catch (Error & e) { 01003 throw e; 01004 } 01005 catch (exception & e) { 01006 throw 01007 InternalErr(string 01008 ("A C++ exception was thrown from inside geoarray(): ") 01009 + e.what()); 01010 01011 } 01012 01013 throw InternalErr(__FILE__, __LINE__, "Impossible condition in geoarray."); 01014 } 01015 01016 void register_functions(ConstraintEvaluator & ce) 01017 { 01018 ce.add_function("grid", function_grid); 01019 ce.add_function("geogrid", function_geogrid); 01020 ce.add_function("linear_scale", function_linear_scale); 01021 ce.add_function("geoarray", function_geoarray); 01022 ce.add_function("version", function_version); 01023 } 01024 01025 } // namespace libdap