libdap++  Updated for version 3.8.2
ce_functions.cc
Go to the documentation of this file.
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 24370 2011-03-28 16:21:32Z 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 #undef FUNCTION_DAP     // undef so the dap() function always returns an error;
00053                         // use keywords instead.
00054 
00055 #include "BaseType.h"
00056 #include "Byte.h"
00057 #include "Int16.h"
00058 #include "UInt16.h"
00059 #include "Int32.h"
00060 #include "UInt32.h"
00061 #include "Float32.h"
00062 #include "Float64.h"
00063 #include "Str.h"
00064 #include "Url.h"
00065 #include "Array.h"
00066 #include "Structure.h"
00067 #include "Sequence.h"
00068 #include "Grid.h"
00069 #include "Error.h"
00070 #include "RValue.h"
00071 
00072 #include "GSEClause.h"
00073 #include "GridGeoConstraint.h"
00074 #include "ArrayGeoConstraint.h"
00075 
00076 #include "ce_functions.h"
00077 #include "gse_parser.h"
00078 #include "gse.tab.hh"
00079 #include "debug.h"
00080 #include "util.h"
00081 
00082 //  We wrapped VC++ 6.x strtod() to account for a short coming
00083 //  in that function in regards to "NaN".  I don't know if this
00084 //  still applies in more recent versions of that product.
00085 //  ROM - 12/2007
00086 #ifdef WIN32
00087 #include <limits>
00088 double w32strtod(const char *, char **);
00089 #endif
00090 
00091 using namespace std;
00092 
00093 int gse_parse(void *arg);
00094 void gse_restart(FILE * in);
00095 
00096 // Glue routines declared in gse.lex
00097 void gse_switch_to_buffer(void *new_buffer);
00098 void gse_delete_buffer(void *buffer);
00099 void *gse_string(const char *yy_str);
00100 
00101 namespace libdap {
00102 
00104 inline bool double_eq(double lhs, double rhs, double epsilon = 1.0e-5)
00105 {
00106     if (lhs > rhs)
00107         return (lhs - rhs) < ((lhs + rhs) / epsilon);
00108     else
00109         return (rhs - lhs) < ((lhs + rhs) / epsilon);
00110 }
00111 
00119 string extract_string_argument(BaseType * arg)
00120 {
00121     if (arg->type() != dods_str_c)
00122         throw Error(malformed_expr,
00123                 "The function requires a DAP string argument.");
00124 
00125     if (!arg->read_p())
00126         throw InternalErr(__FILE__, __LINE__,
00127                 "The CE Evaluator built an argument list where some constants held no values.");
00128 
00129     string s = dynamic_cast<Str&>(*arg).value();
00130 
00131     DBG(cerr << "s: " << s << endl);
00132 
00133     return s;
00134 }
00135 
00136 template<class T> static void set_array_using_double_helper(Array * a,
00137         double *src, int src_len)
00138 {
00139     T *values = new T[src_len];
00140     for (int i = 0; i < src_len; ++i)
00141         values[i] = (T) src[i];
00142 
00143 #ifdef VAL2BUF
00144     a->val2buf(values, true);
00145 #else
00146     a->set_value(values, src_len);
00147 #endif
00148 
00149     delete[]values;
00150 }
00151 
00169 void set_array_using_double(Array * dest, double *src, int src_len)
00170 {
00171     // Simple types are Byte, ..., Float64, String and Url.
00172     if ((dest->type() == dods_array_c && !dest->var()->is_simple_type()) 
00173         || dest->var()->type() == dods_str_c 
00174         || dest->var()->type() == dods_url_c)
00175         throw InternalErr(__FILE__, __LINE__,
00176                 "The function requires a DAP numeric-type array argument.");
00177 
00178     // Test sizes. Note that Array::length() takes any constraint into account
00179     // when it returns the length. Even if this was removed, the 'helper'
00180     // function this uses calls Vector::val2buf() which uses Vector::width()
00181     // which in turn uses length().
00182     if (dest->length() != src_len)
00183         throw InternalErr(__FILE__, __LINE__,
00184                 "The source and destination array sizes don't match ("
00185                 + long_to_string(src_len) + " versus "
00186                 + long_to_string(dest->length()) + ").");
00187 
00188     // The types of arguments that the CE Parser will build for numeric
00189     // constants are limited to Uint32, Int32 and Float64. See ce_expr.y.
00190     // Expanded to work for any numeric type so it can be used for more than
00191     // just arguments.
00192     switch (dest->var()->type()) {
00193     case dods_byte_c:
00194         set_array_using_double_helper<dods_byte>(dest, src, src_len);
00195         break;
00196     case dods_uint16_c:
00197         set_array_using_double_helper<dods_uint16>(dest, src, src_len);
00198         break;
00199     case dods_int16_c:
00200         set_array_using_double_helper<dods_int16>(dest, src, src_len);
00201         break;
00202     case dods_uint32_c:
00203         set_array_using_double_helper<dods_uint32>(dest, src, src_len);
00204         break;
00205     case dods_int32_c:
00206         set_array_using_double_helper<dods_int32>(dest, src, src_len);
00207         break;
00208     case dods_float32_c:
00209         set_array_using_double_helper<dods_float32>(dest, src, src_len);
00210         break;
00211     case dods_float64_c:
00212         set_array_using_double_helper<dods_float64>(dest, src, src_len);
00213         break;
00214     default:
00215         throw InternalErr(__FILE__, __LINE__,
00216                 "The argument list built by the CE parser contained an unsupported numeric type.");
00217     }
00218 
00219     // Set the read_p property.
00220     dest->set_read_p(true);
00221 }
00222 
00223 template<class T> static double *extract_double_array_helper(Array * a)
00224 {
00225     int length = a->length();
00226 
00227     T *b = new T[length];
00228     a->value(b);
00229 
00230     double *dest = new double[length];
00231     for (int i = 0; i < length; ++i)
00232         dest[i] = (double) b[i];
00233     delete[]b;
00234 
00235     return dest;
00236 }
00237 
00242 double *extract_double_array(Array * a)
00243 {
00244     // Simple types are Byte, ..., Float64, String and Url.
00245     if ((a->type() == dods_array_c && !a->var()->is_simple_type())
00246         || a->var()->type() == dods_str_c || a->var()->type() == dods_url_c)
00247         throw Error(malformed_expr,
00248                 "The function requires a DAP numeric-type array argument.");
00249 
00250     if (!a->read_p())
00251         throw InternalErr(__FILE__, __LINE__,
00252                 string("The Array '") + a->name() +
00253                 "'does not contain values.");
00254 
00255     // The types of arguments that the CE Parser will build for numeric
00256     // constants are limited to Uint32, Int32 and Float64. See ce_expr.y.
00257     // Expanded to work for any numeric type so it can be used for more than
00258     // just arguments.
00259     switch (a->var()->type()) {
00260     case dods_byte_c:
00261         return extract_double_array_helper<dods_byte>(a);
00262     case dods_uint16_c:
00263         return extract_double_array_helper<dods_uint16>(a);
00264     case dods_int16_c:
00265         return extract_double_array_helper<dods_int16>(a);
00266     case dods_uint32_c:
00267         return extract_double_array_helper<dods_uint32>(a);
00268     case dods_int32_c:
00269         return extract_double_array_helper<dods_int32>(a);
00270     case dods_float32_c:
00271         return extract_double_array_helper<dods_float32>(a);
00272     case dods_float64_c:
00273         return extract_double_array_helper<dods_float64>(a);
00274     default:
00275         throw InternalErr(__FILE__, __LINE__,
00276                 "The argument list built by the CE parser contained an unsupported numeric type.");
00277     }
00278 }
00279 
00287 double extract_double_value(BaseType * arg)
00288 {
00289     // Simple types are Byte, ..., Float64, String and Url.
00290     if (!arg->is_simple_type() || arg->type() == dods_str_c || arg->type()
00291             == dods_url_c)
00292         throw Error(malformed_expr,
00293                 "The function requires a DAP numeric-type argument.");
00294 
00295     if (!arg->read_p())
00296         throw InternalErr(__FILE__, __LINE__,
00297                 "The CE Evaluator built an argument list where some constants held no values.");
00298 
00299     // The types of arguments that the CE Parser will build for numeric
00300     // constants are limited to Uint32, Int32 and Float64. See ce_expr.y.
00301     // Expanded to work for any numeric type so it can be used for more than
00302     // just arguments.
00303     switch (arg->type()) {
00304     case dods_byte_c:
00305         return (double)(dynamic_cast<Byte&>(*arg).value());
00306     case dods_uint16_c:
00307         return (double)(dynamic_cast<UInt16&>(*arg).value());
00308     case dods_int16_c:
00309         return (double)(dynamic_cast<Int16&>(*arg).value());
00310     case dods_uint32_c:
00311         return (double)(dynamic_cast<UInt32&>(*arg).value());
00312     case dods_int32_c:
00313         return (double)(dynamic_cast<Int32&>(*arg).value());
00314     case dods_float32_c:
00315         return (double)(dynamic_cast<Float32&>(*arg).value());
00316     case dods_float64_c:
00317         return dynamic_cast<Float64&>(*arg).value();
00318     default:
00319         throw InternalErr(__FILE__, __LINE__,
00320                 "The argument list built by the CE parser contained an unsupported numeric type.");
00321     }
00322 }
00323 
00330 void
00331 function_version(int, BaseType *[], DDS &, BaseType **btpp)
00332 {
00333     string
00334             xml_value =
00335                     "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\
00336                        <functions>\
00337                        <function name=\"geogrid\" version=\"1.2\"/>\
00338                        <function name=\"grid\" version=\"1.0\"/>\
00339                        <function name=\"linear_scale\" version=\"1.0b1\"/>\
00340                        <function name=\"version\" version=\"1.0\"/>\
00341                        <function name=\"dap\" version=\"1.0\"/>\
00342                      </functions>";
00343 
00344     //                        <function name=\"geoarray\" version=\"0.9b1\"/>
00345 
00346     Str *response = new Str("version");
00347 
00348     response->set_value(xml_value);
00349     *btpp = response;
00350     return;
00351 }
00352 
00353 void
00354 function_dap(int, BaseType *[], DDS &, ConstraintEvaluator &)
00355 {
00356 #ifdef FUNCTION_DAP
00357     if (argc != 1) {
00358         throw Error("The 'dap' function must be called with a version number.\n\
00359         see http://docs.opendap.org/index.php/Server_Side_Processing_Functions#dap");
00360     }
00361 
00362     double pv = extract_double_value(argv[0]);
00363     dds.set_dap_version(pv);
00364 #else
00365     throw Error("The 'dap' function is not supported in lieu of Constraint expression 'keywords.'\n\
00366 see http://docs.opendap.org/index.php/Server_Side_Processing_Functions#keywords");
00367 #endif
00368 }
00369 
00370 static void parse_gse_expression(gse_arg * arg, BaseType * expr)
00371 {
00372     gse_restart(0); // Restart the scanner.
00373     void *cls = gse_string(extract_string_argument(expr).c_str());
00374     // gse_switch_to_buffer(cls); // Get set to scan the string.
00375     bool status = gse_parse((void *) arg) == 0;
00376     gse_delete_buffer(cls);
00377     if (!status)
00378         throw Error(malformed_expr, "Error parsing grid selection.");
00379 }
00380 
00381 static void apply_grid_selection_expr(Grid * grid, GSEClause * clause)
00382 {
00383     // Basic plan: For each map, look at each clause and set start and stop
00384     // to be the intersection of the ranges in those clauses.
00385     Grid::Map_iter map_i = grid->map_begin();
00386     while (map_i != grid->map_end() && (*map_i)->name() != clause->get_map_name())
00387         ++map_i;
00388 
00389     if (map_i == grid->map_end())
00390         throw Error(malformed_expr,"The map vector '" + clause->get_map_name()
00391                 + "' is not in the grid '" + grid->name() + "'.");
00392 
00393     // Use pointer arith & the rule that map order must match array dim order
00394     Array::Dim_iter grid_dim = (grid->get_array()->dim_begin() + (map_i - grid->map_begin()));
00395 
00396     Array *map = dynamic_cast < Array * >((*map_i));
00397     if (!map)
00398         throw InternalErr(__FILE__, __LINE__, "Expected an Array");
00399     int start = max(map->dimension_start(map->dim_begin()), clause->get_start());
00400     int stop = min(map->dimension_stop(map->dim_begin()), clause->get_stop());
00401 
00402     if (start > stop) {
00403         ostringstream msg;
00404         msg
00405                 << "The expressions passed to grid() do not result in an inclusive \n"
00406                 << "subset of '" << clause->get_map_name()
00407                 << "'. The map's values range " << "from "
00408                 << clause->get_map_min_value() << " to "
00409                 << clause->get_map_max_value() << ".";
00410         throw Error(malformed_expr,msg.str());
00411     }
00412 
00413     DBG(cerr << "Setting constraint on " << map->name()
00414             << "[" << start << ":" << stop << "]" << endl);
00415 
00416     // Stride is always one.
00417     map->add_constraint(map->dim_begin(), start, 1, stop);
00418     grid->get_array()->add_constraint(grid_dim, start, 1, stop);
00419 }
00420 
00421 static void apply_grid_selection_expressions(Grid * grid,
00422         vector < GSEClause * >clauses)
00423 {
00424     vector < GSEClause * >::iterator clause_i = clauses.begin();
00425     while (clause_i != clauses.end())
00426         apply_grid_selection_expr(grid, *clause_i++);
00427 
00428     grid->set_read_p(false);
00429 }
00430 
00467 void
00468 function_grid(int argc, BaseType * argv[], DDS &, BaseType **btpp)
00469 {
00470     DBG(cerr << "Entering function_grid..." << endl);
00471 
00472     string info =
00473     string("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n") +
00474     "<function name=\"grid\" version=\"1.0\" href=\"http://docs.opendap.org/index.php/Server_Side_Processing_Functions#grid\">\n" +
00475     "</function>\n";
00476 
00477     if (argc == 0) {
00478         Str *response = new Str("info");
00479         response->set_value(info);
00480         *btpp = response;
00481         return;
00482     }
00483 
00484     Grid *original_grid = dynamic_cast < Grid * >(argv[0]);
00485     if (!original_grid)
00486         throw Error(malformed_expr,"The first argument to grid() must be a Grid variable!");
00487 
00488     // Duplicate the grid; ResponseBuilder::send_data() will delete the variable
00489     // after serializing it.
00490     BaseType *btp = original_grid->ptr_duplicate();
00491     Grid *l_grid = dynamic_cast < Grid * >(btp);
00492     if (!l_grid) {
00493         delete btp;
00494         throw InternalErr(__FILE__, __LINE__, "Expected a Grid.");
00495     }
00496     
00497     DBG(cerr << "grid: past initialization code" << endl);
00498 
00499     // Read the maps. Do this before calling parse_gse_expression(). Avoid
00500     // reading the array until the constraints have been applied because it
00501     // might be really large.
00502 
00503     // This version makes sure to set the send_p flags which is needed for
00504     // the hdf4 handler (and is what should be done in general).
00505     Grid::Map_iter i = l_grid->map_begin();
00506     while (i != l_grid->map_end())
00507         (*i++)->set_send_p(true);
00508     l_grid->read();
00509 
00510     DBG(cerr << "grid: past map read" << endl);
00511 
00512     // argv[1..n] holds strings; each are little expressions to be parsed.
00513     // When each expression is parsed, the parser makes a new instance of
00514     // GSEClause. GSEClause checks to make sure the named map really exists
00515     // in the Grid and that the range of values given makes sense.
00516     vector < GSEClause * > clauses;
00517     gse_arg *arg = new gse_arg(l_grid);
00518     for (int i = 1; i < argc; ++i) {
00519         parse_gse_expression(arg, argv[i]);
00520         clauses.push_back(arg->get_gsec());
00521     }
00522     delete arg;
00523     arg = 0;
00524 
00525     apply_grid_selection_expressions(l_grid, clauses);
00526 
00527     DBG(cerr << "grid: past gse application" << endl);
00528 
00529     l_grid->get_array()->set_send_p(true);
00530 
00531     l_grid->read();
00532 
00533     *btpp = l_grid;
00534     return;
00535 }
00536 
00571 void
00572 function_geogrid(int argc, BaseType * argv[], DDS &, BaseType **btpp)
00573 {
00574     string info =
00575     string("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n") +
00576     "<function name=\"geogrid\" version=\"1.2\" href=\"http://docs.opendap.org/index.php/Server_Side_Processing_Functions#geogrid\">\n"+
00577     "</function>";
00578 
00579     if (argc == 0) {
00580         Str *response = new Str("version");
00581         response->set_value(info);
00582         *btpp = response;
00583         return ;
00584     }
00585 
00586     // There are two main forms of this function, one that takes a Grid and one
00587     // that takes a Grid and two Arrays. The latter provides a way to explicitly
00588     // tell the function which maps contain lat and lon data. The remaining
00589     // arguments are the same for both versions, although that includes a
00590     // varying argument list.
00591 
00592     // Look at the types of the first three arguments to determine which of the
00593     // two forms were used to call this function.
00594     Grid *l_grid = 0;
00595     if (argc < 1 || !(l_grid = dynamic_cast < Grid * >(argv[0]->ptr_duplicate())))
00596         throw Error(malformed_expr,"The first argument to geogrid() must be a Grid variable!");
00597 
00598     // Both forms require at least this many args
00599     if (argc < 5)
00600         throw Error(malformed_expr,"Wrong number of arguments to geogrid() (expected at least 5 args). See geogrid() for more information.");
00601 
00602     bool grid_lat_lon_form;
00603     Array *l_lat = 0;
00604     Array *l_lon = 0;
00605     if (!(l_lat = dynamic_cast < Array * >(argv[1]))) //->ptr_duplicate())))
00606         grid_lat_lon_form = false;
00607     else if (!(l_lon = dynamic_cast < Array * >(argv[2]))) //->ptr_duplicate())))
00608         throw Error(malformed_expr,"When using the Grid, Lat, Lon form of geogrid() both the lat and lon maps must be given (lon map missing)!");
00609     else
00610         grid_lat_lon_form = true;
00611 
00612     if (grid_lat_lon_form && argc < 7)
00613         throw Error(malformed_expr,"Wrong number of arguments to geogrid() (expected at least 7 args). See geogrid() for more information.");
00614 
00615 #if 0
00616     Grid *l_grid = dynamic_cast < Grid * >(argv[0]->ptr_duplicate());
00617     if (!l_grid)
00618         throw Error(malformed_expr,"The first argument to geogrid() must be a Grid variable!");
00619 #endif
00620     // Read the maps. Do this before calling parse_gse_expression(). Avoid
00621     // reading the array until the constraints have been applied because it
00622     // might be really large.
00623     //
00624     // Trick: Some handlers build Grids from a combination of Array
00625     // variables and attributes. Those handlers (e.g., hdf4) use the send_p
00626     // property to determine which parts of the Grid to read *but they can
00627     // only read the maps from within Grid::read(), not the map's read()*.
00628     // Since the Grid's array does not have send_p set, it will not be read
00629     // by the call below to Grid::read().
00630     Grid::Map_iter i = l_grid->map_begin();
00631     while (i != l_grid->map_end())
00632         (*i++)->set_send_p(true);
00633 
00634     l_grid->read();
00635     // Calling read() above sets the read_p flag for the entire grid; clear it
00636     // for the grid's array so that later on the code will be sure to read it
00637     // under all circumstances.
00638     l_grid->get_array()->set_read_p(false);
00639     DBG(cerr << "geogrid: past map read" << endl);
00640 
00641     // Look for Grid Selection Expressions tacked onto the end of the BB
00642     // specification. If there are any, evaluate them before evaluating the BB.
00643     int min_arg_count = (grid_lat_lon_form) ? 7 : 5;
00644     if (argc > min_arg_count) {
00645         // argv[5..n] holds strings; each are little Grid Selection Expressions
00646         // to be parsed and evaluated.
00647         vector < GSEClause * > clauses;
00648         gse_arg *arg = new gse_arg(l_grid);
00649         for (int i = min_arg_count; i < argc; ++i) {
00650             parse_gse_expression(arg, argv[i]);
00651             clauses.push_back(arg->get_gsec());
00652         }
00653         delete arg;
00654         arg = 0;
00655 
00656         apply_grid_selection_expressions(l_grid, clauses);
00657     }
00658 
00659     try {
00660         // Build a GeoConstraint object. If there are no longitude/latitude
00661         // maps then this constructor throws Error.
00662         GridGeoConstraint gc(l_grid);
00663 
00664         // This sets the bounding box and modifies the maps to match the
00665         // notation of the box (0/359 or -180/179)
00666         int box_index_offset = (grid_lat_lon_form) ? 3 : 1;
00667         double top = extract_double_value(argv[box_index_offset]);
00668         double left = extract_double_value(argv[box_index_offset + 1]);
00669         double bottom = extract_double_value(argv[box_index_offset + 2]);
00670         double right = extract_double_value(argv[box_index_offset + 3]);
00671         gc.set_bounding_box(top, left, bottom, right);
00672         DBG(cerr << "geogrid: past bounding box set" << endl);
00673 
00674         // This also reads all of the data into the grid variable
00675         gc.apply_constraint_to_data();
00676         DBG(cerr << "geogrid: past apply constraint" << endl);
00677 
00678         // In this function the l_grid pointer is the same as the pointer returned
00679         // by this call. The caller of the function must free the pointer.
00680         *btpp = gc.get_constrained_grid();
00681         return;
00682     }
00683     catch (Error &e) {
00684         throw e;
00685     }
00686     catch (exception & e) {
00687         throw
00688         InternalErr(string
00689                 ("A C++ exception was thrown from inside geogrid(): ")
00690                 + e.what());
00691     }
00692 }
00693 
00694 // These static functions could be moved to a class that provides a more
00695 // general interface for COARDS/CF someday. Assume each BaseType comes bundled
00696 // with an attribute table.
00697 
00698 // This was ripped from parser-util.cc
00699 static double string_to_double(const char *val)
00700 {
00701     char *ptr;
00702     errno = 0;
00703     // Clear previous value. 5/21/2001 jhrg
00704 
00705 #ifdef WIN32
00706     double v = w32strtod(val, &ptr);
00707 #else
00708     double v = strtod(val, &ptr);
00709 #endif
00710 
00711     if ((v == 0.0 && (val == ptr || errno == HUGE_VAL || errno == ERANGE))
00712             || *ptr != '\0') {
00713         throw Error(malformed_expr,string("Could not convert the string '") + val + "' to a double.");
00714     }
00715 
00716     double abs_val = fabs(v);
00717     if (abs_val > DODS_DBL_MAX || (abs_val != 0.0 && abs_val < DODS_DBL_MIN))
00718         throw Error(malformed_expr,string("Could not convert the string '") + val + "' to a double.");
00719 
00720     return v;
00721 }
00722 
00732 static double get_attribute_double_value(BaseType *var,
00733         vector<string> &attributes)
00734 {
00735     // This code also builds a list of the attribute values that have been
00736     // passed in but not found so that an informative message can be returned.
00737     AttrTable &attr = var->get_attr_table();
00738     string attribute_value = "";
00739     string values = "";
00740     vector<string>::iterator i = attributes.begin();
00741     while (attribute_value == "" && i != attributes.end()) {
00742         values += *i;
00743         if (!values.empty())
00744             values += ", ";
00745         attribute_value = attr.get_attr(*i++);
00746     }
00747 
00748     // If the value string is empty, then look at the grid's array (if it's a
00749     // grid) or throw an Error.
00750     if (attribute_value.empty()) {
00751         if (var->type() == dods_grid_c)
00752             return get_attribute_double_value(dynamic_cast<Grid&>(*var).get_array(), attributes);
00753         else
00754             throw Error(malformed_expr,string("No COARDS/CF '") + values.substr(0, values.length() - 2)
00755                     + "' attribute was found for the variable '"
00756                     + var->name() + "'.");
00757     }
00758 
00759     return string_to_double(remove_quotes(attribute_value).c_str());
00760 }
00761 
00762 static double get_attribute_double_value(BaseType *var, const string &attribute)
00763 {
00764     AttrTable &attr = var->get_attr_table();
00765     string attribute_value = attr.get_attr(attribute);
00766 
00767     // If the value string is empty, then look at the grid's array (if it's a
00768     // grid or throw an Error.
00769     if (attribute_value.empty()) {
00770         if (var->type() == dods_grid_c)
00771             return get_attribute_double_value(dynamic_cast<Grid&>(*var).get_array(), attribute);
00772         else
00773             throw Error(malformed_expr,string("No COARDS '") + attribute
00774                     + "' attribute was found for the variable '"
00775                     + var->name() + "'.");
00776     }
00777 
00778     return string_to_double(remove_quotes(attribute_value).c_str());
00779 }
00780 
00781 static double get_y_intercept(BaseType *var)
00782 {
00783     vector<string> attributes;
00784     attributes.push_back("add_offset");
00785     attributes.push_back("add_off");
00786     return get_attribute_double_value(var, attributes);
00787 }
00788 
00789 static double get_slope(BaseType *var)
00790 {
00791     return get_attribute_double_value(var, "scale_factor");
00792 }
00793 
00794 static double get_missing_value(BaseType *var)
00795 {
00796     return get_attribute_double_value(var, "missing_value");
00797 }
00798 
00811 void
00812 function_linear_scale(int argc, BaseType * argv[], DDS &, BaseType **btpp)
00813 {
00814     string info =
00815     string("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n") +
00816     "<function name=\"linear_scale\" version=\"1.0b1\" href=\"http://docs.opendap.org/index.php/Server_Side_Processing_Functions#linear_scale\">\n" +
00817     "</function>";
00818 
00819     if (argc == 0) {
00820         Str *response = new Str("info");
00821         response->set_value(info);
00822         *btpp = response;
00823         return;
00824     }
00825 
00826     // Check for 1 or 3 arguments: 1 --> use attributes; 3 --> m & b supplied
00827     DBG(cerr << "argc = " << argc << endl);
00828     if (!(argc == 1 || argc == 3 || argc == 4))
00829         throw Error(malformed_expr,"Wrong number of arguments to linear_scale(). See linear_scale() for more information");
00830 
00831     // Get m & b
00832     bool use_missing = false;
00833     double m, b, missing = 0.0;
00834     if (argc == 4) {
00835         m = extract_double_value(argv[1]);
00836         b = extract_double_value(argv[2]);
00837         missing = extract_double_value(argv[3]);
00838         use_missing = true;
00839     } else if (argc == 3) {
00840         m = extract_double_value(argv[1]);
00841         b = extract_double_value(argv[2]);
00842         use_missing = false;
00843     } else {
00844         m = get_slope(argv[0]);
00845 
00846         // This is really a hack; on a fair number of datasets, the y intercept
00847         // is not given and is assumed to be 0. Here the function looks and
00848         // catches the error if a y intercept is not found.
00849         try {
00850             b = get_y_intercept(argv[0]);
00851         }
00852         catch (Error &e) {
00853             b = 0.0;
00854         }
00855 
00856         // This is not the best plan; the get_missing_value() function should
00857         // do something other than throw, but to do that would require mayor
00858         // surgery on get_attribute_double_value().
00859         try {
00860             missing = get_missing_value(argv[0]);
00861             use_missing = true;
00862         }
00863         catch (Error &e) {
00864             use_missing = false;
00865         }
00866     }
00867 
00868     DBG(cerr << "m: " << m << ", b: " << b << endl);DBG(cerr << "use_missing: " << use_missing << ", missing: " << missing << endl);
00869 
00870     // Read the data, scale and return the result. Must replace the new data
00871     // in a constructor (i.e., Array part of a Grid).
00872     BaseType *dest = 0;
00873     double *data;
00874     if (argv[0]->type() == dods_grid_c) {
00875         Array &source = *dynamic_cast<Grid&>(*argv[0]).get_array();
00876         source.set_send_p(true);
00877         source.read();
00878         data = extract_double_array(&source);
00879         int length = source.length();
00880         int i = 0;
00881         while (i < length) {
00882             DBG2(cerr << "data[" << i << "]: " << data[i] << endl);
00883             if (!use_missing || !double_eq(data[i], missing))
00884                 data[i] = data[i] * m + b;
00885             DBG2(cerr << " >> data[" << i << "]: " << data[i] << endl);
00886             ++i;
00887         }
00888 
00889         // Vector::add_var will delete the existing 'template' variable
00890         Float64 *temp_f = new Float64(source.name());
00891         source.add_var(temp_f);
00892 #ifdef VAL2BUF
00893         source.val2buf(static_cast<void*>(data), false);
00894 #else
00895         source.set_value(data, i);
00896 #endif
00897         delete [] data; // val2buf copies.
00898         delete temp_f; // add_var copies and then adds.
00899         dest = argv[0];
00900     } else if (argv[0]->is_vector_type()) {
00901         Array &source = dynamic_cast<Array&>(*argv[0]);
00902         source.set_send_p(true);
00903         // If the array is really a map, make sure to read using the Grid
00904         // because of the HDF4 handler's odd behavior WRT dimensions.
00905         if (source.get_parent() && source.get_parent()->type() == dods_grid_c)
00906             source.get_parent()->read();
00907         else
00908             source.read();
00909 
00910         data = extract_double_array(&source);
00911         int length = source.length();
00912         int i = 0;
00913         while (i < length) {
00914             if (!use_missing || !double_eq(data[i], missing))
00915                 data[i] = data[i] * m + b;
00916             ++i;
00917         }
00918 
00919         Float64 *temp_f = new Float64(source.name());
00920         source.add_var(temp_f);
00921 
00922         source.val2buf(static_cast<void*>(data), false);
00923 
00924         delete [] data; // val2buf copies.
00925         delete temp_f; // add_var copies and then adds.
00926 
00927         dest = argv[0];
00928     } else if (argv[0]->is_simple_type() && !(argv[0]->type() == dods_str_c
00929             || argv[0]->type() == dods_url_c)) {
00930         double data = extract_double_value(argv[0]);
00931         if (!use_missing || !double_eq(data, missing))
00932             data = data * m + b;
00933 
00934         dest = new Float64(argv[0]->name());
00935 
00936         dest->val2buf(static_cast<void*>(&data));
00937 
00938     } else {
00939         throw Error(malformed_expr,"The linear_scale() function works only for numeric Grids, Arrays and scalars.");
00940     }
00941 
00942     *btpp = dest;
00943     return;
00944 }
00945 
00946 #if 0
00947 
00964 void
00965 function_geoarray(int argc, BaseType * argv[], DDS &, BaseType **btpp)
00966 {
00967     string info =
00968     string("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n") +
00969     "<function name=\"geoarray\" version=\"0.9b1\" href=\"http://docs.opendap.org/index.php/Server_Side_Processing_Functions#geoarray\">\n" +
00970     "</function>";
00971 
00972     if (argc == 0) {
00973         Str *response = new Str("version");
00974         response->set_value(info);
00975         *btpp = response;
00976         return;
00977     }
00978 
00979     DBG(cerr << "argc = " << argc << endl);
00980     if (!(argc == 5 || argc == 9 || argc == 11))
00981         throw Error(malformed_expr,"Wrong number of arguments to geoarray(). See geoarray() for more information.");
00982 
00983     // Check the Array (and dup because the caller will free the variable).
00984     Array *l_array = dynamic_cast < Array * >(argv[0]->ptr_duplicate());
00985     if (!l_array)
00986         throw Error(malformed_expr,"The first argument to geoarray() must be an Array variable!");
00987 
00988     try {
00989 
00990         // Read the bounding box and variable extents from the params
00991         double bb_top = extract_double_value(argv[1]);
00992         double bb_left = extract_double_value(argv[2]);
00993         double bb_bottom = extract_double_value(argv[3]);
00994         double bb_right = extract_double_value(argv[4]);
00995 
00996         switch (argc) {
00997             case 5: {
00998                 ArrayGeoConstraint agc(l_array);
00999 
01000                         agc.set_bounding_box(bb_left, bb_top, bb_right, bb_bottom);
01001                                 // This also reads all of the data into the grid variable
01002                         agc.apply_constraint_to_data();
01003                         *btpp = agc.get_constrained_array();
01004                         return;
01005                 break;
01006             }
01007             case 9: {
01008                 double var_top = extract_double_value(argv[5]);
01009                 double var_left = extract_double_value(argv[6]);
01010                 double var_bottom = extract_double_value(argv[7]);
01011                 double var_right = extract_double_value(argv[8]);
01012                 ArrayGeoConstraint agc (l_array, var_left, var_top, var_right, var_bottom);
01013 
01014                         agc.set_bounding_box(bb_left, bb_top, bb_right, bb_bottom);
01015                                 // This also reads all of the data into the grid variable
01016                         agc.apply_constraint_to_data();
01017                         *btpp =  agc.get_constrained_array();
01018                         return;
01019                 break;
01020             }
01021             case 11: {
01022                 double var_top = extract_double_value(argv[5]);
01023                 double var_left = extract_double_value(argv[6]);
01024                 double var_bottom = extract_double_value(argv[7]);
01025                 double var_right = extract_double_value(argv[8]);
01026                 string projection = extract_string_argument(argv[9]);
01027                 string datum = extract_string_argument(argv[10]);
01028                 ArrayGeoConstraint agc(l_array,
01029                         var_left, var_top, var_right, var_bottom,
01030                         projection, datum);
01031 
01032                         agc.set_bounding_box(bb_left, bb_top, bb_right, bb_bottom);
01033                                 // This also reads all of the data into the grid variable
01034                         agc.apply_constraint_to_data();
01035                         *btpp = agc.get_constrained_array();
01036                         return;
01037                 break;
01038             }
01039             default:
01040                 throw InternalErr(__FILE__, __LINE__, "Wrong number of args to geoarray.");
01041         }
01042     }
01043     catch (Error & e) {
01044         throw e;
01045     }
01046     catch (exception & e) {
01047         throw
01048         InternalErr(string
01049                 ("A C++ exception was thrown from inside geoarray(): ")
01050                 + e.what());
01051 
01052     }
01053 
01054     throw InternalErr(__FILE__, __LINE__, "Impossible condition in geoarray.");
01055 }
01056 #endif
01057 
01058 void register_functions(ConstraintEvaluator & ce)
01059 {
01060     ce.add_function("grid", function_grid);
01061     ce.add_function("geogrid", function_geogrid);
01062     ce.add_function("linear_scale", function_linear_scale);
01063 #if 0
01064     ce.add_function("geoarray", function_geoarray);
01065 #endif
01066     ce.add_function("version", function_version);
01067 
01068     ce.add_function("dap", function_dap);
01069 }
01070 
01071 } // namespace libdap