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