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

Generated on Wed Nov 14 03:15:43 2007 for libdap++ by  doxygen 1.5.1