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