ce_functions.cc

Go to the documentation of this file.
00001 
00002 // -*- mode: c++; c-basic-offset:4 -*-
00003 
00004 // This file is part of libdap, A C++ implementation of the OPeNDAP Data
00005 // Access Protocol.
00006 
00007 // Copyright (c) 2002,2003 OPeNDAP, Inc.
00008 // Author: James Gallagher <jgallagher@opendap.org>
00009 //
00010 // This library is free software; you can redistribute it and/or
00011 // modify it under the terms of the GNU Lesser General Public
00012 // License as published by the Free Software Foundation; either
00013 // version 2.1 of the License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful,
00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00023 //
00024 // You can contact OPeNDAP, Inc. at PO Box 112, Saunderstown, RI. 02874-0112.
00025 
00026 // (c) COPYRIGHT URI/MIT 1999
00027 // Please read the full copyright statement in the file COPYRIGHT_URI.
00028 //
00029 // Authors:
00030 //      jhrg,jimg       James Gallagher <jgallagher@gso.uri.edu>
00031 
00032 
00033 // These functions are used by the CE evaluator
00034 //
00035 // 1/15/99 jhrg
00036 
00037 #include "config.h"
00038 
00039 static char rcsid[] not_used =
00040     { "$Id: ce_functions.cc 16612 2007-06-04 22:08:57Z jimg $"
00041     };
00042 
00043 #include <errno.h>      // used by strtod()
00044 #include <limits.h>
00045 #include <math.h>
00046 
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.h"
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 using std::vector;
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 {
00090 
00092 inline bool
00093 double_eq(double lhs, double rhs, double epsilon = 1.0e-5)
00094 {
00095     if (lhs > rhs)
00096         return (lhs - rhs) < ((lhs + rhs) / epsilon);
00097     else
00098         return (rhs - lhs) < ((lhs + rhs) / epsilon);
00099 }
00100 
00108 string extract_string_argument(BaseType * arg)
00109 {
00110     if (arg->type() != dods_str_c)
00111         throw Error(malformed_expr,
00112                     "The function requires a DAP string-type argument.");
00113 
00114     if (!arg->read_p())
00115         throw InternalErr(__FILE__, __LINE__,
00116                           "The CE Evaluator built an argument list where some constants held no values.");
00117 
00118     string *sp = 0;
00119     arg->buf2val((void **) &sp);
00120     string s = *sp;
00121     delete sp;
00122     sp = 0;
00123 
00124     DBG(cerr << "s: " << s << endl);
00125 
00126     return s;
00127 }
00128 template < class T >
00129 static void
00130 set_array_using_double_helper(Array * a, double *src,
00131                               int src_len)
00132 {
00133     T *values = new T[src_len];
00134     for (int i = 0; i < src_len; ++i)
00135         values[i] = (T) src[i];
00136     a->val2buf(values, true);
00137     delete[]values;
00138 }
00139 
00157 void set_array_using_double(Array * dest, double *src, int src_len)
00158 {
00159     // Simple types are Byte, ..., Float64, String and Url.
00160     if (dest->type() == dods_array_c && !dest->var()->is_simple_type()
00161         || dest->var()->type() == dods_str_c
00162         || dest->var()->type() == dods_url_c)
00163         throw InternalErr(__FILE__, __LINE__,
00164                           "The function requires a DAP numeric-type array argument.");
00165 
00166     // Test sizes. Note that Array::length() takes any constraint into account
00167     // when it returns the length. Even if this was removed, the 'helper'
00168     // function this uses calls Vector::val2buf() which uses Vector::width()
00169     // which in turn uses length().
00170     if (dest->length() != src_len)
00171         throw InternalErr(__FILE__, __LINE__,
00172                           "The source and destination array sizes don't match ("
00173                           + long_to_string(src_len) + " versus "
00174                           + long_to_string(dest->length()) + ").");
00175 
00176     // The types of arguments that the CE Parser will build for numeric
00177     // constants are limited to Uint32, Int32 and Float64. See ce_expr.y.
00178     // Expanded to work for any numeric type so it can be used for more than
00179     // just arguments.
00180     switch (dest->var()->type()) {
00181     case dods_byte_c:
00182         set_array_using_double_helper < dods_byte > (dest, src,
00183                 src_len);
00184         break;
00185     case dods_uint16_c:
00186         set_array_using_double_helper < dods_uint16 > (dest, src,
00187                 src_len);
00188         break;
00189     case dods_int16_c:
00190         set_array_using_double_helper < dods_int16 > (dest, src,
00191                 src_len);
00192         break;
00193     case dods_uint32_c:
00194         set_array_using_double_helper < dods_uint32 > (dest, src,
00195                 src_len);
00196         break;
00197     case dods_int32_c:
00198         set_array_using_double_helper < dods_int32 > (dest, src,
00199                 src_len);
00200         break;
00201     case dods_float32_c:
00202         set_array_using_double_helper < dods_float32 > (dest, src,
00203                 src_len);
00204         break;
00205     case dods_float64_c:
00206         set_array_using_double_helper < dods_float64 > (dest, src,
00207                 src_len);
00208         break;
00209     default:
00210         throw InternalErr(__FILE__, __LINE__,
00211                           "The argument list built by the CE parser contained an unsupported numeric type.");
00212     }
00213 
00214     // Set the read_p property.
00215     dest->set_read_p(true);
00216 }
00217 
00218 template < class T >
00219 static double *extract_double_array_helper(Array * a)
00220 {
00221     int length = a->length();
00222     double *dest = new double[length];
00223     T *b = new T[length];
00224     a->buf2val((void **) &b);
00225     for (int i = 0; i < length; ++i)
00226         dest[i] = (double) b[i];
00227     delete[]b;
00228     return dest;
00229 }
00230 
00235 double *extract_double_array(Array * a)
00236 {
00237     // Simple types are Byte, ..., Float64, String and Url.
00238     if (a->type() == dods_array_c && !a->var()->is_simple_type()
00239         || a->var()->type() == dods_str_c
00240         || a->var()->type() == dods_url_c)
00241         throw Error(malformed_expr,
00242                     "The function requires a DAP numeric-type array argument.");
00243 
00244     if (!a->read_p())
00245         throw InternalErr(__FILE__, __LINE__,
00246                           string("The Array '") + a->name() +
00247                           "'does not contain values.");
00248 
00249     // The types of arguments that the CE Parser will build for numeric
00250     // constants are limited to Uint32, Int32 and Float64. See ce_expr.y.
00251     // Expanded to work for any numeric type so it can be used for more than
00252     // just arguments.
00253     switch (a->var()->type()) {
00254     case dods_byte_c:
00255         return extract_double_array_helper < dods_byte > (a);
00256     case dods_uint16_c:
00257         return extract_double_array_helper < dods_uint16 > (a);
00258     case dods_int16_c:
00259         return extract_double_array_helper < dods_int16 > (a);
00260     case dods_uint32_c:
00261         return extract_double_array_helper < dods_uint32 > (a);
00262     case dods_int32_c:
00263         return extract_double_array_helper < dods_int32 > (a);
00264     case dods_float32_c:
00265         return extract_double_array_helper < dods_float32 > (a);
00266     case dods_float64_c:
00267         return extract_double_array_helper < dods_float64 > (a);
00268     default:
00269         throw InternalErr(__FILE__, __LINE__,
00270                           "The argument list built by the CE parser contained an unsupported numeric type.");
00271     }
00272 }
00273 
00281 double extract_double_value(BaseType * arg)
00282 {
00283     // Simple types are Byte, ..., Float64, String and Url.
00284     if (!arg->is_simple_type()
00285         || arg->type() == dods_str_c || arg->type() == dods_url_c)
00286         throw Error(malformed_expr,
00287                     "The function requires a DAP numeric-type argument.");
00288 
00289     if (!arg->read_p())
00290         throw InternalErr(__FILE__, __LINE__,
00291                           "The CE Evaluator built an argument list where some constants held no values.");
00292 
00293     // The types of arguments that the CE Parser will build for numeric
00294     // constants are limited to Uint32, Int32 and Float64. See ce_expr.y.
00295     // Expanded to work for any numeric type so it can be used for more than
00296     // just arguments.
00297     switch (arg->type()) {
00298     case dods_byte_c: {
00299             dods_byte i;
00300             dods_byte *pi = &i;
00301             arg->buf2val((void **) &pi);
00302             return (double)(i);
00303         }
00304     case dods_uint16_c: {
00305             dods_uint16 i;
00306             dods_uint16 *pi = &i;
00307             arg->buf2val((void **) &pi);
00308             return (double)(i);
00309         }
00310     case dods_int16_c: {
00311             dods_int16 i;
00312             dods_int16 *pi = &i;
00313             arg->buf2val((void **) &pi);
00314             return (double)(i);
00315         }
00316     case dods_uint32_c: {
00317             dods_uint32 i;
00318             dods_uint32 *pi = &i;
00319             arg->buf2val((void **) &pi);
00320             return (double)(i);
00321         }
00322     case dods_int32_c: {
00323             dods_int32 i;
00324             dods_int32 *pi = &i;
00325             arg->buf2val((void **) &pi);
00326             return (double)(i);
00327         }
00328     case dods_float32_c: {
00329             dods_float32 i;
00330             dods_float32 *pi = &i;
00331             arg->buf2val((void **) &pi);
00332             return (double) i;
00333         }
00334     case dods_float64_c: {
00335             DBG(cerr << "arg->value(): " << dynamic_cast <
00336                 Float64 * >(arg)->value() << endl);
00337             dods_float64 i;
00338             dods_float64 *pi = &i;
00339             arg->buf2val((void **) &pi);
00340             DBG(cerr << "i: " << i << endl);
00341             return i;
00342         }
00343     default:
00344         throw InternalErr(__FILE__, __LINE__,
00345                           "The argument list built by the CE parser contained an unsupported numeric type.");
00346     }
00347 }
00348 
00349 #if 0
00350 // In reality no server implements this; it _should_ be removed. 03/28/05 jhrg
00351 BaseType *func_length(int argc, BaseType * argv[], DDS & dds)
00352 {
00353     if (argc != 1) {
00354         throw Error("Wrong number of arguments to length().");
00355     }
00356 
00357     switch (argv[0]->type()) {
00358 
00359     case dods_sequence_c: {
00360             Sequence *var = dynamic_cast < Sequence * >(argv[0]);
00361             if (!var)
00362                 throw
00363                 Error("Expected a Sequence variable in length()");
00364             dods_int32 result = var->length();
00365 
00366             BaseType *ret = dds.get_factory()->NewInt32("constant");
00367             ret->val2buf(&result);
00368             ret->set_read_p(true);
00369             ret->set_send_p(true);
00370 
00371             return ret;
00372         }
00373 
00374     default:
00375         throw Error("Wrong type argument to length()");
00376     }
00377 
00378     return 0;
00379 }
00380 #endif
00381 
00384 BaseType *function_version(int, BaseType *[], DDS &, const string &)
00385 {
00386     string xml_value = "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\
00387                        <functions>\
00388                        <function name=\"version\" version=\"1.0\"/>\
00389                        <function name=\"grid\" version=\"1.0\"/>\
00390                        <function name=\"geogrid\" version=\"1.0b2\"/>\
00391                        <function name=\"geoarray\" version=\"0.9b1\"/>\
00392                        <function name=\"linear_scale\" version=\"1.0b1\"/>\
00393                        </functions>";
00394 
00395     Str *response = new Str("version");
00396 
00397     response->set_value(xml_value);
00398 
00399     return response;
00400 }
00401 
00402 static void parse_gse_expression(gse_arg * arg, BaseType * expr)
00403 {
00404     gse_restart(0);         // Restart the scanner.
00405     void *cls = gse_string(extract_string_argument(expr).c_str());
00406     // gse_switch_to_buffer(cls); // Get set to scan the string.
00407     bool status = gse_parse((void *) arg) == 0;
00408     gse_delete_buffer(cls);
00409     if (!status)
00410         throw Error(malformed_expr, "Error parsing grid selection.");
00411 }
00412 
00413 static void apply_grid_selection_expr(Grid * grid, GSEClause * clause)
00414 {
00415     // Basic plan: For each map, look at each clause and set start and stop
00416     // to be the intersection of the ranges in those clauses.
00417     Grid::Map_iter map_i = grid->map_begin();
00418     while (map_i != grid->map_end()
00419            && (*map_i)->name() != clause->get_map_name())
00420         ++map_i;
00421 
00422     if (map_i == grid->map_end())
00423         throw Error("The map vector '" + clause->get_map_name()
00424                     + "' is not in the grid '" + grid->name() + "'.");
00425 
00426     // Use pointer arith & the rule that map order must match array dim order
00427     Array::Dim_iter grid_dim = (grid->get_array()->dim_begin()
00428                                 + (map_i - grid->map_begin()));
00429 
00430     Array *map = dynamic_cast < Array * >((*map_i));
00431     int start =
00432         max(map->dimension_start(map->dim_begin()),
00433             clause->get_start());
00434     int stop =
00435         min(map->dimension_stop(map->dim_begin()), clause->get_stop());
00436 
00437     if (start > stop) {
00438         ostringstream msg;
00439         msg <<
00440         "The expressions passed to grid() do not result in an inclusive \n"
00441         << "subset of '" << clause->get_map_name()
00442         << "'. The map's values range " << "from " <<
00443         clause->get_map_min_value() << " to "
00444         << clause->get_map_max_value() << ".";
00445         throw Error(msg.str());
00446     }
00447 
00448     DBG(cerr << "Setting constraint on " << map->name()
00449         << "[" << start << ":" << stop << "]" << endl);
00450 
00451     // Stride is always one.
00452     map->add_constraint(map->dim_begin(), start, 1, stop);
00453     grid->get_array()->add_constraint(grid_dim, start, 1, stop);
00454 }
00455 
00456 static void
00457 apply_grid_selection_expressions(Grid * grid,
00458                                  vector < GSEClause * >clauses)
00459 {
00460     vector < GSEClause * >::iterator clause_i = clauses.begin();
00461     while (clause_i != clauses.end())
00462         apply_grid_selection_expr(grid, *clause_i++);
00463 
00464     grid->set_read_p(false);
00465 }
00466 
00503 BaseType *function_grid(int argc, BaseType * argv[], DDS &,
00504                         const string & dataset)
00505 {
00506     DBG(cerr << "Entering function_grid..." << endl);
00507 
00508     string info =
00509         "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\
00510         <function name=\"grid\" version=\"1.0\"/>\
00511         The grid() function takes a grid variable and zero or more relational\
00512         expressions. Each relational expression is applied to the grid using the\
00513         server's constraint evaluator and the resulting grid is returned. The\
00514         expressions may use constants and the grid's map vectors but may not use\
00515         any other variables. Two forms of expression are provide: \"var relop const\"\
00516         and \"const relop var relop const\". For example: grid(sst, \"10<=TIME<20\")\
00517         and grid(sst, \"10<=TIME\", \"TIME<20\") are both legal and, in this case,\
00518         also equivalent.\
00519         </function>";
00520 
00521     if (argc == 0) {
00522         Str *response = new Str("info");
00523         response->set_value(info);
00524         return response;
00525     }
00526 
00527     Grid *original_grid = dynamic_cast < Grid * >(argv[0]);
00528     if (!original_grid)
00529         throw Error("The first argument to grid() must be a Grid variable!");
00530 
00531     // Duplicate the grid; DODSFilter::send_data() will delete the variable
00532     // after serializing it.
00533     Grid *l_grid =
00534         dynamic_cast < Grid * >(original_grid->ptr_duplicate());
00535 
00536     DBG(cerr << "grid: past initialization code" << endl);
00537 
00538     // Read the maps. Do this before calling parse_gse_expression(). Avoid
00539     // reading the array until the constraints have been applied because it
00540     // might be really large.
00541 
00542     // This version makes sure to set the send_p flags which is needed for
00543     // the hdf4 handler (and is what should be done in general).
00544     Grid::Map_iter i = l_grid->map_begin();
00545     while (i != l_grid->map_end())
00546         (*i++)->set_send_p(true);
00547     l_grid->read(dataset);
00548 
00549 #if 0
00550     // Old, pre-hdf_handler compatible version.
00551     Grid::Map_iter i = l_grid->map_begin();
00552     while (i != l_grid->map_end())
00553         (*i++)->read(dataset);
00554 #endif
00555     DBG(cerr << "grid: past map read" << endl);
00556 
00557     // argv[1..n] holds strings; each are little expressions to be parsed.
00558     // When each expression is parsed, the parser makes a new instance of
00559     // GSEClause. GSEClause checks to make sure the named map really exists
00560     // in the Grid and that the range of values given makes sense.
00561     vector < GSEClause * >clauses;
00562     gse_arg *arg = new gse_arg(l_grid);
00563     for (int i = 1; i < argc; ++i) {
00564         parse_gse_expression(arg, argv[i]);
00565         clauses.push_back(arg->get_gsec());
00566     }
00567     delete arg;
00568     arg = 0;
00569 
00570     apply_grid_selection_expressions(l_grid, clauses);
00571 
00572     DBG(cerr << "grid: past gse application" << endl);
00573 
00574     l_grid->get_array()->set_send_p(true);
00575 
00576     l_grid->read(dataset);
00577 
00578     return l_grid;
00579 }
00580 
00616 BaseType *function_geogrid(int argc, BaseType * argv[], DDS &,
00617                            const string & dataset)
00618 {
00619     string info =
00620         "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\
00621         <function name=\"geogrid\" version=\"1.0b2\"/>\
00622         geogrid() applies a constraint given in latitude and longitude to a\
00623         DAP Grid variable. The arguments to the function are:\
00624         geogrid(<grid variable>, <upper latitude>, <left longitude>,\
00625         <lower latitude>, <right longitude> [selection expressions - see grid()])\
00626         geogrid(\"version\") returns the version of the function.\
00627         The function will always return a single Grid variable whose values\
00628         completely cover the given region, although there may be cases when\
00629         some additional data is also returned. If the longitude values 'wrap\
00630         around' the right edge of the data, then the function will make two\
00631         requests and return those joined together as a single Grid.\
00632         </function>";
00633 
00634     if (argc == 0) {
00635         Str *response = new Str("version");
00636         response->set_value(info);
00637         return response;
00638     }
00639 
00640     if (argc < 5)
00641         throw Error("Wrong number of arguments to geogrid(). See geogrid() for more information.");
00642 
00643     Grid *l_grid = dynamic_cast < Grid * >(argv[0]->ptr_duplicate());
00644     if (!l_grid)
00645         throw Error("The first argument to geogrid() must be a Grid variable!");
00646 
00647     // Read the maps. Do this before calling parse_gse_expression(). Avoid
00648     // reading the array until the constraints have been applied because it
00649     // might be really large.
00650     //
00651     // Trick: Some handlers build Grids from a combination of Array
00652     // variables and attributes. Those handlers (e.g., hdf4) use the send_p
00653     // property to determine which parts of the Grid to read *but they can
00654     // only read the maps from within Grid::read(), not the map's read()*.
00655     // Since the Grid's array does not have send_p set, it will not be read
00656     // by the call below to Grid::read().
00657     Grid::Map_iter i = l_grid->map_begin();
00658     while (i != l_grid->map_end())
00659         (*i++)->set_send_p(true);
00660     l_grid->read(dataset);
00661     // Calling read() above sets the read_p flag for the entire grid; clear it
00662     // for the grid's array so that later on the code will be sure to read it
00663     // under all circumstances.
00664     l_grid->get_array()->set_read_p(false);
00665     DBG(cerr << "geogrid: past map read" << endl);
00666 
00667     // Look for Grid Selection Expressions tacked onto the end of the BB
00668     // specification. If there are any, evaluate them before evaluating the BB.
00669     if (argc > 5) {
00670         // argv[5..n] holds strings; each are little Grid Selection Expressions
00671         // to be parsed and evaluated.
00672         vector < GSEClause * >clauses;
00673         gse_arg *arg = new gse_arg(l_grid);
00674         for (int i = 5; i < argc; ++i) {
00675             parse_gse_expression(arg, argv[i]);
00676             clauses.push_back(arg->get_gsec());
00677         }
00678         delete arg;
00679         arg = 0;
00680 
00681         apply_grid_selection_expressions(l_grid, clauses);
00682     }
00683 
00684     try {
00685         // Build a GeoConstraint object. If there are no longitude/latitude
00686         // maps then this constructor throws Error.
00687         GridGeoConstraint gc(l_grid, dataset);
00688 
00689         // This sets the bounding box and modifies the maps to match the
00690         // notation of the box (0/359 or -180/179)
00691         double top = extract_double_value(argv[1]);
00692         double left = extract_double_value(argv[2]);
00693         double bottom = extract_double_value(argv[3]);
00694         double right = extract_double_value(argv[4]);
00695         gc.set_bounding_box(left, top, right, bottom);
00696         DBG(cerr << "geogrid: past bounding box set" << endl);
00697 
00698         // This also reads all of the data into the grid variable
00699         gc.apply_constraint_to_data();
00700         DBG(cerr << "geogrid: past apply constraint" << endl);
00701 
00702         // In this function the l_grid pointer is the same as the pointer returned
00703         // by this call. The caller of the function must free the pointer.
00704         return gc.get_constrained_grid();
00705     }
00706     catch (Error & e) {
00707         throw;
00708     }
00709     catch (exception & e) {
00710         throw
00711         InternalErr(string
00712                     ("A C++ exception was thrown from inside geogrid(): ")
00713                     + e.what());
00714     }
00715 }
00716 
00717 // These static functions could be moved to a class that provides a more
00718 // general interface for COARDS/CF someday. Assume each BaseType comes bundled
00719 // with an attribute table.
00720 
00721 // This was ripped from parser-util.cc
00722 static double
00723 string_to_double(const char *val)
00724 {
00725     char *ptr;
00726     errno = 0;                  // Clear previous value. 5/21/2001 jhrg
00727 
00728 #ifdef WIN32
00729     double v = w32strtod(val, &ptr);
00730 #else
00731     double v = strtod(val, &ptr);
00732 #endif
00733 
00734     if ((v == 0.0 && (val == ptr || errno == HUGE_VAL || errno == ERANGE))
00735         || *ptr != '\0') {
00736         throw Error(string("Could not convert the string '") + val + "' to a double.");
00737     }
00738 
00739     double abs_val = fabs(v);
00740     if (abs_val > DODS_DBL_MAX || (abs_val != 0.0 && abs_val < DODS_DBL_MIN))
00741         throw Error(string("Could not convert the string '") + val + "' to a double.");
00742 
00743     return v;
00744 }
00745 
00746 static string
00747 remove_quotes(const string &s)
00748 {
00749     if (s[0] == '\"' && s[s.length()-1] == '\"')
00750         return s.substr(1, s.length() - 2);
00751     else
00752         return s;
00753 }
00754 
00758 static double
00759 get_attribute_double_value(BaseType *var, vector<string> &attributes)
00760 {
00761     AttrTable &attr = var->get_attr_table();
00762     string attribute_value = "";
00763     string values = "";
00764     vector<string>::iterator i = attributes.begin();
00765     while (attribute_value == "" && i != attributes.end()) {
00766         values += *i; if (!values.empty()) values += ", ";
00767         attribute_value = attr.get_attr(*i++);
00768     }
00769 
00770     // If the value string is empty, then look at the grid's array (if it's a
00771     // grid or throw an Error.
00772     if (attribute_value.empty()) {
00773         if (var->type() == dods_grid_c)
00774             return get_attribute_double_value(dynamic_cast<Grid&>(*var).get_array(), attributes);
00775         else
00776             throw Error(string("No COARDS '") + values.substr(0, values.length() - 2)
00777                         + "' attribute was found for the variable '"
00778                         + var->name() + "'.");
00779     }
00780 
00781     return string_to_double(remove_quotes(attribute_value).c_str());
00782 }
00783 
00784 static double
00785 get_attribute_double_value(BaseType *var, const string &attribute)
00786 {
00787     AttrTable &attr = var->get_attr_table();
00788     string attribute_value = attr.get_attr(attribute);
00789 
00790     // If the value string is empty, then look at the grid's array (if it's a
00791     // grid or throw an Error.
00792     if (attribute_value.empty()) {
00793         if (var->type() == dods_grid_c)
00794             return get_attribute_double_value(dynamic_cast<Grid&>(*var).get_array(), attribute);
00795         else
00796             throw Error(string("No COARDS '") + attribute
00797                         + "' attribute was found for the variable '"
00798                         + var->name() + "'.");
00799     }
00800 
00801     return string_to_double(remove_quotes(attribute_value).c_str());
00802 }
00803 
00804 static double
00805 get_y_intercept(BaseType *var)
00806 {
00807     vector<string> attributes;
00808     attributes.push_back("add_offset");
00809     attributes.push_back("add_off");
00810     return get_attribute_double_value(var, attributes);
00811 }
00812 
00813 static double
00814 get_slope(BaseType *var)
00815 {
00816     return get_attribute_double_value(var, "scale_factor");
00817 }
00818 
00819 static double
00820 get_missing_value(BaseType *var)
00821 {
00822     return get_attribute_double_value(var, "missing_value");
00823 }
00824 
00837 BaseType *
00838 function_linear_scale(int argc, BaseType * argv[], DDS &, const string & dataset)
00839 {
00840     string info =
00841         "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\
00842         <function name=\"linear_scale\" version=\"1.0b1\">\
00843         The linear_scale() function applies the familiar y=mx+b equation to data.\
00844         If only the name of a variable is given, the function looks for the COARDS\
00845         scale_factor, add_offset and missing_value attributes. In the equation,\
00846         'm' is scale_factor, 'b' is add_offset and data values which match\
00847         missing_value are not scaled. If add_offset cannot be found, it defaults to\
00848         zero; if missing_value cannot be found, the test for it is not performed.\
00849         The caller can also provide values for these and, if provided, those values\
00850         are used and the dataset's attributes are ignored. If only scale_factor is\
00851         provided, the defaults for 'b' and the missing value flag are used.\
00852         Similarly, if only 'm' and 'b' are provided the missing value test is not\
00853         performed. The values are always returned in a Float64 array.\
00854         </function>";
00855 
00856     if (argc == 0) {
00857         Str *response = new Str("info");
00858         response->set_value(info);
00859         return response;
00860     }
00861 
00862     // Check for 1 or 3 arguments: 1 --> use attributes; 3 --> m & b supplied
00863     DBG(cerr << "argc = " << argc << endl);
00864     if (!(argc == 1 || argc == 3 || argc == 4))
00865         throw Error("Wrong number of arguments to linear_scale(). See linear_scale() for more information");
00866 
00867     // Get m & b
00868     bool use_missing;
00869     double m, b, missing;
00870     if (argc == 4) {
00871         m = extract_double_value(argv[1]);
00872         b = extract_double_value(argv[2]);
00873         missing = extract_double_value(argv[3]);
00874         use_missing = true;
00875     }
00876     else if (argc == 3) {
00877         m = extract_double_value(argv[1]);
00878         b = extract_double_value(argv[2]);
00879         use_missing = false;
00880     }
00881     else {
00882         m = get_slope(argv[0]);
00883 
00884         // This is really a hack; on a fair number of datasets, the y intercept
00885         // is not given and is assumed to be 0. Here the function looks and
00886         // catches the error if a y intercept is not found.
00887         try {
00888             b = get_y_intercept(argv[0]);
00889         }
00890         catch (Error &e) {
00891             b = 0.0;
00892         }
00893 
00894         // This is not the best plan; the get_missing_value() function should
00895         // do something other than throw, but to do that would require mayor
00896         // surgery on get_attribute_double_value().
00897         try {
00898             missing = get_missing_value(argv[0]);
00899             use_missing = true;
00900         }
00901         catch (Error &e) {
00902             use_missing = false;
00903         }
00904     }
00905 
00906     DBG(cerr << "m: " << m << ", b: " << b << endl);
00907     DBG(cerr << "use_missing: " << use_missing << ", missing: " << missing << endl);
00908 
00909     // Read the data, scale and return the result. Must replace the new data
00910     // in a constructor (i.e., Array part of a Grid).
00911     BaseType *dest = 0;
00912     double *data;
00913     if (argv[0]->type() == dods_grid_c) {
00914         Array &source = *dynamic_cast<Grid&>(*argv[0]).get_array();
00915         source.set_send_p(true);
00916         source.read(dataset);
00917         data = extract_double_array(&source);
00918         int length = source.length();
00919         int i = 0;
00920         while (i < length) {
00921             DBG2(cerr << "data[" << i << "]: " << data[i] << endl);
00922             if (!use_missing || !double_eq(data[i], missing))
00923                 data[i] = data[i] * m + b;
00924             DBG2(cerr << " >> data[" << i << "]: " << data[i] << endl);
00925             ++i;
00926         }
00927 
00928         // Vector::add_var will delete the existing 'template' variable
00929         Float64 *temp_f = new Float64(source.name());
00930         source.add_var(temp_f);
00931         source.val2buf(static_cast<void*>(data), false);
00932         delete temp_f;              // add_var copies and then adds.
00933         dest = argv[0];
00934     }
00935     else if (argv[0]->is_vector_type()) {
00936         Array &source = dynamic_cast<Array&>(*argv[0]);
00937         source.set_send_p(true);
00938         // If the array is really a map, make sure to read using the Grid
00939         // because of the HDF4 handler's odd behavior WRT dimensions.
00940         if (source.get_parent()
00941             && source.get_parent()->type() == dods_grid_c)
00942             source.get_parent()->read(dataset);
00943         else
00944             source.read(dataset);
00945 
00946         data = extract_double_array(&source);
00947         int length = source.length();
00948         int i = 0;
00949         while (i < length) {
00950             if (!use_missing || !double_eq(data[i], missing))
00951                 data[i] = data[i] * m + b;
00952             ++i;
00953         }
00954 
00955         Float64 *temp_f = new Float64(source.name());
00956         source.add_var(temp_f);
00957         source.val2buf(static_cast<void*>(data), false);
00958         delete temp_f;              // add_var copies and then adds.
00959 
00960         dest = argv[0];
00961     }
00962     else if (argv[0]->is_simple_type()
00963              && !(argv[0]->type() == dods_str_c
00964                   || argv[0]->type() == dods_url_c)) {
00965         double data = extract_double_value(argv[0]);
00966         if (!use_missing || !double_eq(data, missing))
00967             data = data * m + b;
00968 
00969         dest = new Float64(argv[0]->name());
00970         dest->val2buf(static_cast<void*>(&data));
00971     }
00972     else {
00973         throw Error("The linear_scale() function works only for numeric Grids, Arrays and scalars.");
00974     }
00975 
00976     return dest;
00977 }
00978 
00995 BaseType *
00996 function_geoarray(int argc, BaseType * argv[], DDS &, const string & dataset)
00997 {
00998     string info =
00999         "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\
01000         <function name=\"geoarray\" version=\"0.9b1\"/>\
01001         The geoarray() function supports two different sets of arguments:\
01002         geoarray(var,left,top,right,bottom)\
01003         geoarray(var,left,top,right,bottom,var_left,var_top,var_right,var_bottom)\
01004         In the first version 'var' is the target of the selection and 'left', 'top',\
01005         'right' and 'bottom' are the corners of a longitude-latitude box that defines\
01006         the selection. In the second version the 'var_left', ..., parameters give the\
01007         longitude and latitude extent of the entire array. The projection and datum are\
01008         assumed to be Plat-Carre and WGS84.\
01009         </function>";
01010 
01011     if (argc == 0) {
01012         Str *response = new Str("version");
01013         response->set_value(info);
01014         return response;
01015     }
01016 
01017     DBG(cerr << "argc = " << argc << endl);
01018     if (!(argc == 5 || argc == 9 || argc == 11))
01019         throw Error("Wrong number of arguments to geoarray(). See geoarray() for more information.");
01020 
01021     // Check the Array (and dup because the caller will free the variable).
01022     Array *l_array = dynamic_cast < Array * >(argv[0]->ptr_duplicate());
01023     if (!l_array)
01024         throw Error("The first argument to geoarray() must be an Array variable!");
01025 
01026     try {
01027 
01028         // Read the bounding box and variable extents from the params
01029         double bb_top = extract_double_value(argv[1]);
01030         double bb_left = extract_double_value(argv[2]);
01031         double bb_bottom = extract_double_value(argv[3]);
01032         double bb_right = extract_double_value(argv[4]);
01033 
01034         ArrayGeoConstraint *agc = 0;
01035         switch (argc) {
01036         case 5:
01037             agc = new ArrayGeoConstraint(l_array, dataset);
01038             break;
01039         case 9: {
01040                 double var_top = extract_double_value(argv[5]);
01041                 double var_left = extract_double_value(argv[6]);
01042                 double var_bottom = extract_double_value(argv[7]);
01043                 double var_right = extract_double_value(argv[8]);
01044                 agc = new ArrayGeoConstraint(l_array, dataset,
01045                                              var_left, var_top, var_right, var_bottom);
01046                 break;
01047             }
01048         case 11: {
01049                 double var_top = extract_double_value(argv[5]);
01050                 double var_left = extract_double_value(argv[6]);
01051                 double var_bottom = extract_double_value(argv[7]);
01052                 double var_right = extract_double_value(argv[8]);
01053                 string projection = extract_string_argument(argv[9]);
01054                 string datum = extract_string_argument(argv[10]);
01055                 agc = new ArrayGeoConstraint(l_array, dataset,
01056                                              var_left, var_top, var_right, var_bottom,
01057                                              projection, datum);
01058                 break;
01059             }
01060         default:
01061             throw InternalErr(__FILE__, __LINE__, "Wrong number of args to geoarray.");
01062         }
01063 
01064         agc->set_bounding_box(bb_left, bb_top, bb_right, bb_bottom);
01065 
01066         // This also reads all of the data into the grid variable
01067         agc->apply_constraint_to_data();
01068 
01069         // In this function the l_grid pointer is the same as the pointer returned
01070         // by this call. The caller of the function must free the pointer.
01071         return agc->get_constrained_array();
01072     }
01073     catch (Error & e) {
01074         throw;
01075     }
01076     catch (exception & e) {
01077         throw
01078         InternalErr(string
01079                     ("A C++ exception was thrown from inside geoarray(): ")
01080                     + e.what());
01081 
01082     }
01083 }
01084 
01085 void register_functions(ConstraintEvaluator & ce)
01086 {
01087     ce.add_function("grid", function_grid);
01088     ce.add_function("geogrid", function_geogrid);
01089     ce.add_function("linear_scale", function_linear_scale);
01090     ce.add_function("geoarray", function_geoarray);
01091     ce.add_function("version", function_version);
01092 }
01093 
01094 }                               // namespace libdap

Generated on Wed Jun 27 12:56:38 2007 for libdap++ by  doxygen 1.4.7