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