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

Generated on Tue Jun 10 18:00:29 2008 for libdap++ by  doxygen 1.5.4