libdap++  Updated for version 3.8.2
ce_functions.cc
Go to the documentation of this file.
1 // -*- mode: c++; c-basic-offset:4 -*-
2 
3 // This file is part of libdap, A C++ implementation of the OPeNDAP Data
4 // Access Protocol.
5 
6 // Copyright (c) 2002,2003 OPeNDAP, Inc.
7 // Author: James Gallagher <jgallagher@opendap.org>
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the GNU Lesser General Public
11 // License as published by the Free Software Foundation; either
12 // version 2.1 of the License, or (at your option) any later version.
13 //
14 // This library is distributed in the hope that it will be useful,
15 // but WITHOUT ANY WARRANTY; without even the implied warranty of
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 // Lesser General Public License for more details.
18 //
19 // You should have received a copy of the GNU Lesser General Public
20 // License along with this library; if not, write to the Free Software
21 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 //
23 // You can contact OPeNDAP, Inc. at PO Box 112, Saunderstown, RI. 02874-0112.
24 
25 // (c) COPYRIGHT URI/MIT 1999
26 // Please read the full copyright statement in the file COPYRIGHT_URI.
27 //
28 // Authors:
29 // jhrg,jimg James Gallagher <jgallagher@gso.uri.edu>
30 
31 
32 // These functions are used by the CE evaluator
33 //
34 // 1/15/99 jhrg
35 
36 #include "config.h"
37 
38 static char rcsid[]not_used = { "$Id: ce_functions.cc 25066 2011-11-30 18:39:37Z jimg $" };
39 
40 #include <limits.h>
41 
42 #include <cstdlib> // used by strtod()
43 #include <cerrno>
44 #include <cmath>
45 #include <iostream>
46 #include <vector>
47 #include <algorithm>
48 
49 //#define DODS_DEBUG
50 #undef FUNCTION_DAP // undef so the dap() function always returns an error;
51 // use keywords instead.
52 
53 #include "BaseType.h"
54 #include "Byte.h"
55 #include "Int16.h"
56 #include "UInt16.h"
57 #include "Int32.h"
58 #include "UInt32.h"
59 #include "Float32.h"
60 #include "Float64.h"
61 #include "Str.h"
62 #include "Url.h"
63 #include "Array.h"
64 #include "Structure.h"
65 #include "Sequence.h"
66 #include "Grid.h"
67 #include "Error.h"
68 #include "RValue.h"
69 
70 #include "GSEClause.h"
71 #include "GridGeoConstraint.h"
72 #include "ArrayGeoConstraint.h"
73 
74 #include "ce_functions.h"
75 #include "gse_parser.h"
76 #include "gse.tab.hh"
77 #include "debug.h"
78 #include "util.h"
79 
80 // We wrapped VC++ 6.x strtod() to account for a short coming
81 // in that function in regards to "NaN". I don't know if this
82 // still applies in more recent versions of that product.
83 // ROM - 12/2007
84 #ifdef WIN32
85 #include <limits>
86 double w32strtod(const char *, char **);
87 #endif
88 
89 using namespace std;
90 
91 int gse_parse(void *arg);
92 void gse_restart(FILE * in);
93 
94 // Glue routines declared in gse.lex
95 void gse_switch_to_buffer(void *new_buffer);
96 void gse_delete_buffer(void *buffer);
97 void *gse_string(const char *yy_str);
98 
99 namespace libdap {
100 
102 inline bool double_eq(double lhs, double rhs, double epsilon = 1.0e-5)
103 {
104  if (lhs > rhs)
105  return (lhs - rhs) < ((lhs + rhs) / epsilon);
106  else
107  return (rhs - lhs) < ((lhs + rhs) / epsilon);
108 }
109 
118 {
119  if (arg->type() != dods_str_c)
120  throw Error(malformed_expr, "The function requires a DAP string argument.");
121 
122  if (!arg->read_p())
123  throw InternalErr(__FILE__, __LINE__,
124  "The CE Evaluator built an argument list where some constants held no values.");
125 
126  string s = dynamic_cast<Str&> (*arg).value();
127 
128  DBG(cerr << "s: " << s << endl);
129 
130  return s;
131 }
132 
133 template<class T> static void set_array_using_double_helper(Array * a, double *src, int src_len)
134 {
135  T *values = new T[src_len];
136  for (int i = 0; i < src_len; ++i)
137  values[i] = (T) src[i];
138 
139 #ifdef VAL2BUF
140  a->val2buf(values, true);
141 #else
142  a->set_value(values, src_len);
143 #endif
144 
145  delete[] values;
146 }
147 
165 void set_array_using_double(Array * dest, double *src, int src_len)
166 {
167  // Simple types are Byte, ..., Float64, String and Url.
168  if ((dest->type() == dods_array_c && !dest->var()->is_simple_type()) || dest->var()->type() == dods_str_c
169  || dest->var()->type() == dods_url_c)
170  throw InternalErr(__FILE__, __LINE__, "The function requires a DAP numeric-type array argument.");
171 
172  // Test sizes. Note that Array::length() takes any constraint into account
173  // when it returns the length. Even if this was removed, the 'helper'
174  // function this uses calls Vector::val2buf() which uses Vector::width()
175  // which in turn uses length().
176  if (dest->length() != src_len)
177  throw InternalErr(
178  __FILE__,
179  __LINE__,
180  "The source and destination array sizes don't match (" + long_to_string(src_len) + " versus "
181  + long_to_string(dest->length()) + ").");
182 
183  // The types of arguments that the CE Parser will build for numeric
184  // constants are limited to Uint32, Int32 and Float64. See ce_expr.y.
185  // Expanded to work for any numeric type so it can be used for more than
186  // just arguments.
187  switch (dest->var()->type()) {
188  case dods_byte_c:
189  set_array_using_double_helper<dods_byte> (dest, src, src_len);
190  break;
191  case dods_uint16_c:
192  set_array_using_double_helper<dods_uint16> (dest, src, src_len);
193  break;
194  case dods_int16_c:
195  set_array_using_double_helper<dods_int16> (dest, src, src_len);
196  break;
197  case dods_uint32_c:
198  set_array_using_double_helper<dods_uint32> (dest, src, src_len);
199  break;
200  case dods_int32_c:
201  set_array_using_double_helper<dods_int32> (dest, src, src_len);
202  break;
203  case dods_float32_c:
204  set_array_using_double_helper<dods_float32> (dest, src, src_len);
205  break;
206  case dods_float64_c:
207  set_array_using_double_helper<dods_float64> (dest, src, src_len);
208  break;
209  default:
210  throw InternalErr(__FILE__, __LINE__,
211  "The argument list built by the CE parser contained an unsupported numeric type.");
212  }
213 
214  // Set the read_p property.
215  dest->set_read_p(true);
216 }
217 
218 template<class T> static double *extract_double_array_helper(Array * a)
219 {
220  int length = a->length();
221 
222  T *b = new T[length];
223  a->value(b);
224 
225  double *dest = new double[length];
226  for (int i = 0; i < length; ++i)
227  dest[i] = (double) b[i];
228  delete[] b;
229 
230  return dest;
231 }
232 
238 {
239  // Simple types are Byte, ..., Float64, String and Url.
240  if ((a->type() == dods_array_c && !a->var()->is_simple_type()) || a->var()->type() == dods_str_c
241  || a->var()->type() == dods_url_c)
242  throw Error(malformed_expr, "The function requires a DAP numeric-type array argument.");
243 
244  if (!a->read_p())
245  throw InternalErr(__FILE__, __LINE__, string("The Array '") + a->name() + "'does not contain values.");
246 
247  // The types of arguments that the CE Parser will build for numeric
248  // constants are limited to Uint32, Int32 and Float64. See ce_expr.y.
249  // Expanded to work for any numeric type so it can be used for more than
250  // just arguments.
251  switch (a->var()->type()) {
252  case dods_byte_c:
253  return extract_double_array_helper<dods_byte> (a);
254  case dods_uint16_c:
255  return extract_double_array_helper<dods_uint16> (a);
256  case dods_int16_c:
257  return extract_double_array_helper<dods_int16> (a);
258  case dods_uint32_c:
259  return extract_double_array_helper<dods_uint32> (a);
260  case dods_int32_c:
261  return extract_double_array_helper<dods_int32> (a);
262  case dods_float32_c:
263  return extract_double_array_helper<dods_float32> (a);
264  case dods_float64_c:
265  return extract_double_array_helper<dods_float64> (a);
266  default:
267  throw InternalErr(__FILE__, __LINE__,
268  "The argument list built by the CE parser contained an unsupported numeric type.");
269  }
270 }
271 
280 {
281  // Simple types are Byte, ..., Float64, String and Url.
282  if (!arg->is_simple_type() || arg->type() == dods_str_c || arg->type() == dods_url_c)
283  throw Error(malformed_expr, "The function requires a DAP numeric-type argument.");
284 
285  if (!arg->read_p())
286  throw InternalErr(__FILE__, __LINE__,
287  "The CE Evaluator built an argument list where some constants held no values.");
288 
289  // The types of arguments that the CE Parser will build for numeric
290  // constants are limited to Uint32, Int32 and Float64. See ce_expr.y.
291  // Expanded to work for any numeric type so it can be used for more than
292  // just arguments.
293  switch (arg->type()) {
294  case dods_byte_c:
295  return (double) (dynamic_cast<Byte&> (*arg).value());
296  case dods_uint16_c:
297  return (double) (dynamic_cast<UInt16&> (*arg).value());
298  case dods_int16_c:
299  return (double) (dynamic_cast<Int16&> (*arg).value());
300  case dods_uint32_c:
301  return (double) (dynamic_cast<UInt32&> (*arg).value());
302  case dods_int32_c:
303  return (double) (dynamic_cast<Int32&> (*arg).value());
304  case dods_float32_c:
305  return (double) (dynamic_cast<Float32&> (*arg).value());
306  case dods_float64_c:
307  return dynamic_cast<Float64&> (*arg).value();
308  default:
309  throw InternalErr(__FILE__, __LINE__,
310  "The argument list built by the CE parser contained an unsupported numeric type.");
311  }
312 }
313 
320 void function_version(int, BaseType *[], DDS &, BaseType **btpp)
321 {
322  string
323  xml_value =
324  "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\
325  <functions>\
326  <function name=\"geogrid\" version=\"1.2\"/>\
327  <function name=\"grid\" version=\"1.0\"/>\
328  <function name=\"linear_scale\" version=\"1.0b1\"/>\
329  <function name=\"version\" version=\"1.0\"/>\
330  <function name=\"dap\" version=\"1.0\"/>\
331  </functions>";
332 
333  // <function name=\"geoarray\" version=\"0.9b1\"/>
334 
335  Str *response = new Str("version");
336 
337  response->set_value(xml_value);
338  *btpp = response;
339  return;
340 }
341 
343 {
344 #ifdef FUNCTION_DAP
345  if (argc != 1) {
346  throw Error("The 'dap' function must be called with a version number.\n\
347  see http://docs.opendap.org/index.php/Server_Side_Processing_Functions#dap");
348  }
349 
350  double pv = extract_double_value(argv[0]);
351  dds.set_dap_version(pv);
352 #else
353  throw Error(
354  "The 'dap' function is not supported in lieu of Constraint expression 'keywords.'\n\
355 see http://docs.opendap.org/index.php/Server_Side_Processing_Functions#keywords");
356 #endif
357 }
358 
359 static void parse_gse_expression(gse_arg * arg, BaseType * expr)
360 {
361  gse_restart(0); // Restart the scanner.
362  void *cls = gse_string(extract_string_argument(expr).c_str());
363  // gse_switch_to_buffer(cls); // Get set to scan the string.
364  bool status = gse_parse((void *) arg) == 0;
365  gse_delete_buffer(cls);
366  if (!status)
367  throw Error(malformed_expr, "Error parsing grid selection.");
368 }
369 
370 static void apply_grid_selection_expr(Grid * grid, GSEClause * clause)
371 {
372  // Basic plan: For each map, look at each clause and set start and stop
373  // to be the intersection of the ranges in those clauses.
374  Grid::Map_iter map_i = grid->map_begin();
375  while (map_i != grid->map_end() && (*map_i)->name() != clause->get_map_name())
376  ++map_i;
377 
378  if (map_i == grid->map_end())
379  throw Error(malformed_expr,
380  "The map vector '" + clause->get_map_name() + "' is not in the grid '" + grid->name() + "'.");
381 
382  // Use pointer arith & the rule that map order must match array dim order
383  Array::Dim_iter grid_dim = (grid->get_array()->dim_begin() + (map_i - grid->map_begin()));
384 
385  Array *map = dynamic_cast<Array *> ((*map_i));
386  if (!map)
387  throw InternalErr(__FILE__, __LINE__, "Expected an Array");
388  int start = max(map->dimension_start(map->dim_begin()), clause->get_start());
389  int stop = min(map->dimension_stop(map->dim_begin()), clause->get_stop());
390 
391  if (start > stop) {
392  ostringstream msg;
393  msg << "The expressions passed to grid() do not result in an inclusive \n" << "subset of '"
394  << clause->get_map_name() << "'. The map's values range " << "from " << clause->get_map_min_value()
395  << " to " << clause->get_map_max_value() << ".";
396  throw Error(malformed_expr, msg.str());
397  }
398 
399  DBG(cerr << "Setting constraint on " << map->name()
400  << "[" << start << ":" << stop << "]" << endl);
401 
402  // Stride is always one.
403  map->add_constraint(map->dim_begin(), start, 1, stop);
404  grid->get_array()->add_constraint(grid_dim, start, 1, stop);
405 }
406 
407 static void apply_grid_selection_expressions(Grid * grid, vector<GSEClause *> clauses)
408 {
409  vector<GSEClause *>::iterator clause_i = clauses.begin();
410  while (clause_i != clauses.end())
411  apply_grid_selection_expr(grid, *clause_i++);
412 
413  grid->set_read_p(false);
414 }
415 
452 void function_grid(int argc, BaseType * argv[], DDS &, BaseType **btpp)
453 {
454  DBG(cerr << "Entering function_grid..." << endl);
455 
456  string
457  info =
458  string("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n")
459  + "<function name=\"grid\" version=\"1.0\" href=\"http://docs.opendap.org/index.php/Server_Side_Processing_Functions#grid\">\n"
460  + "</function>\n";
461 
462  if (argc == 0) {
463  Str *response = new Str("info");
464  response->set_value(info);
465  *btpp = response;
466  return;
467  }
468 
469  Grid *original_grid = dynamic_cast<Grid *> (argv[0]);
470  if (!original_grid)
471  throw Error(malformed_expr, "The first argument to grid() must be a Grid variable!");
472 
473  // Duplicate the grid; ResponseBuilder::send_data() will delete the variable
474  // after serializing it.
475  BaseType *btp = original_grid->ptr_duplicate();
476  Grid *l_grid = dynamic_cast<Grid *> (btp);
477  if (!l_grid) {
478  delete btp;
479  throw InternalErr(__FILE__, __LINE__, "Expected a Grid.");
480  }
481 
482  DBG(cerr << "grid: past initialization code" << endl);
483 
484  // Read the maps. Do this before calling parse_gse_expression(). Avoid
485  // reading the array until the constraints have been applied because it
486  // might be really large.
487 
488  // This version makes sure to set the send_p flags which is needed for
489  // the hdf4 handler (and is what should be done in general).
490  Grid::Map_iter i = l_grid->map_begin();
491  while (i != l_grid->map_end())
492  (*i++)->set_send_p(true);
493  l_grid->read();
494 
495  DBG(cerr << "grid: past map read" << endl);
496 
497  // argv[1..n] holds strings; each are little expressions to be parsed.
498  // When each expression is parsed, the parser makes a new instance of
499  // GSEClause. GSEClause checks to make sure the named map really exists
500  // in the Grid and that the range of values given makes sense.
501  vector<GSEClause *> clauses;
502  gse_arg *arg = new gse_arg(l_grid);
503  for (int i = 1; i < argc; ++i) {
504  parse_gse_expression(arg, argv[i]);
505  clauses.push_back(arg->get_gsec());
506  }
507  delete arg;
508  arg = 0;
509 
510  apply_grid_selection_expressions(l_grid, clauses);
511 
512  DBG(cerr << "grid: past gse application" << endl);
513 
514  l_grid->get_array()->set_send_p(true);
515 
516  l_grid->read();
517 
518  *btpp = l_grid;
519  return;
520 }
521 
556 void function_geogrid(int argc, BaseType * argv[], DDS &, BaseType **btpp)
557 {
558  string
559  info =
560  string("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n")
561  + "<function name=\"geogrid\" version=\"1.2\" href=\"http://docs.opendap.org/index.php/Server_Side_Processing_Functions#geogrid\">\n"
562  + "</function>";
563 
564  if (argc == 0) {
565  Str *response = new Str("version");
566  response->set_value(info);
567  *btpp = response;
568  return;
569  }
570 
571  // There are two main forms of this function, one that takes a Grid and one
572  // that takes a Grid and two Arrays. The latter provides a way to explicitly
573  // tell the function which maps contain lat and lon data. The remaining
574  // arguments are the same for both versions, although that includes a
575  // varying argument list.
576 
577  // Look at the types of the first three arguments to determine which of the
578  // two forms were used to call this function.
579  Grid *l_grid = 0;
580  if (argc < 1 || !(l_grid = dynamic_cast<Grid *> (argv[0]->ptr_duplicate())))
581  throw Error(malformed_expr, "The first argument to geogrid() must be a Grid variable!");
582 
583  // Both forms require at least this many args
584  if (argc < 5)
585  throw Error(malformed_expr,
586  "Wrong number of arguments to geogrid() (expected at least 5 args). See geogrid() for more information.");
587 
588  bool grid_lat_lon_form;
589  Array *l_lat = 0;
590  Array *l_lon = 0;
591  if (!(l_lat = dynamic_cast<Array *> (argv[1]))) //->ptr_duplicate())))
592  grid_lat_lon_form = false;
593  else if (!(l_lon = dynamic_cast<Array *> (argv[2]))) //->ptr_duplicate())))
594  throw Error(malformed_expr,
595  "When using the Grid, Lat, Lon form of geogrid() both the lat and lon maps must be given (lon map missing)!");
596  else
597  grid_lat_lon_form = true;
598 
599  if (grid_lat_lon_form && argc < 7)
600  throw Error(malformed_expr,
601  "Wrong number of arguments to geogrid() (expected at least 7 args). See geogrid() for more information.");
602 
603 #if 0
604  Grid *l_grid = dynamic_cast < Grid * >(argv[0]->ptr_duplicate());
605  if (!l_grid)
606  throw Error(malformed_expr,"The first argument to geogrid() must be a Grid variable!");
607 #endif
608  // Read the maps. Do this before calling parse_gse_expression(). Avoid
609  // reading the array until the constraints have been applied because it
610  // might be really large.
611  //
612  // Trick: Some handlers build Grids from a combination of Array
613  // variables and attributes. Those handlers (e.g., hdf4) use the send_p
614  // property to determine which parts of the Grid to read *but they can
615  // only read the maps from within Grid::read(), not the map's read()*.
616  // Since the Grid's array does not have send_p set, it will not be read
617  // by the call below to Grid::read().
618  Grid::Map_iter i = l_grid->map_begin();
619  while (i != l_grid->map_end())
620  (*i++)->set_send_p(true);
621 
622  l_grid->read();
623  // Calling read() above sets the read_p flag for the entire grid; clear it
624  // for the grid's array so that later on the code will be sure to read it
625  // under all circumstances.
626  l_grid->get_array()->set_read_p(false);
627  DBG(cerr << "geogrid: past map read" << endl);
628 
629  // Look for Grid Selection Expressions tacked onto the end of the BB
630  // specification. If there are any, evaluate them before evaluating the BB.
631  int min_arg_count = (grid_lat_lon_form) ? 7 : 5;
632  if (argc > min_arg_count) {
633  // argv[5..n] holds strings; each are little Grid Selection Expressions
634  // to be parsed and evaluated.
635  vector<GSEClause *> clauses;
636  gse_arg *arg = new gse_arg(l_grid);
637  for (int i = min_arg_count; i < argc; ++i) {
638  parse_gse_expression(arg, argv[i]);
639  clauses.push_back(arg->get_gsec());
640  }
641  delete arg;
642  arg = 0;
643 
644  apply_grid_selection_expressions(l_grid, clauses);
645  }
646 
647  try {
648  // Build a GeoConstraint object. If there are no longitude/latitude
649  // maps then this constructor throws Error.
650  GridGeoConstraint gc(l_grid);
651 
652  // This sets the bounding box and modifies the maps to match the
653  // notation of the box (0/359 or -180/179)
654  int box_index_offset = (grid_lat_lon_form) ? 3 : 1;
655  double top = extract_double_value(argv[box_index_offset]);
656  double left = extract_double_value(argv[box_index_offset + 1]);
657  double bottom = extract_double_value(argv[box_index_offset + 2]);
658  double right = extract_double_value(argv[box_index_offset + 3]);
659  gc.set_bounding_box(top, left, bottom, right);
660  DBG(cerr << "geogrid: past bounding box set" << endl);
661 
662  // This also reads all of the data into the grid variable
664  DBG(cerr << "geogrid: past apply constraint" << endl);
665 
666  // In this function the l_grid pointer is the same as the pointer returned
667  // by this call. The caller of the function must free the pointer.
668  *btpp = gc.get_constrained_grid();
669  return;
670  } catch (Error &e) {
671  throw e;
672  } catch (exception & e) {
673  throw InternalErr(string("A C++ exception was thrown from inside geogrid(): ") + e.what());
674  }
675 }
676 
677 // These static functions could be moved to a class that provides a more
678 // general interface for COARDS/CF someday. Assume each BaseType comes bundled
679 // with an attribute table.
680 
681 // This was ripped from parser-util.cc
682 static double string_to_double(const char *val)
683 {
684  char *ptr;
685  errno = 0;
686  // Clear previous value. 5/21/2001 jhrg
687 
688 #ifdef WIN32
689  double v = w32strtod(val, &ptr);
690 #else
691  double v = strtod(val, &ptr);
692 #endif
693 
694  if ((v == 0.0 && (val == ptr || errno == HUGE_VAL || errno == ERANGE)) || *ptr != '\0') {
695  throw Error(malformed_expr, string("Could not convert the string '") + val + "' to a double.");
696  }
697 
698  double abs_val = fabs(v);
699  if (abs_val > DODS_DBL_MAX || (abs_val != 0.0 && abs_val < DODS_DBL_MIN))
700  throw Error(malformed_expr, string("Could not convert the string '") + val + "' to a double.");
701 
702  return v;
703 }
704 
714 static double get_attribute_double_value(BaseType *var, vector<string> &attributes)
715 {
716  // This code also builds a list of the attribute values that have been
717  // passed in but not found so that an informative message can be returned.
718  AttrTable &attr = var->get_attr_table();
719  string attribute_value = "";
720  string values = "";
721  vector<string>::iterator i = attributes.begin();
722  while (attribute_value == "" && i != attributes.end()) {
723  values += *i;
724  if (!values.empty())
725  values += ", ";
726  attribute_value = attr.get_attr(*i++);
727  }
728 
729  // If the value string is empty, then look at the grid's array (if it's a
730  // grid) or throw an Error.
731  if (attribute_value.empty()) {
732  if (var->type() == dods_grid_c)
733  return get_attribute_double_value(dynamic_cast<Grid&> (*var).get_array(), attributes);
734  else
735  throw Error(
737  string("No COARDS/CF '") + values.substr(0, values.length() - 2)
738  + "' attribute was found for the variable '" + var->name() + "'.");
739  }
740 
741  return string_to_double(remove_quotes(attribute_value).c_str());
742 }
743 
744 static double get_attribute_double_value(BaseType *var, const string &attribute)
745 {
746  AttrTable &attr = var->get_attr_table();
747  string attribute_value = attr.get_attr(attribute);
748 
749  // If the value string is empty, then look at the grid's array (if it's a
750  // grid or throw an Error.
751  if (attribute_value.empty()) {
752  if (var->type() == dods_grid_c)
753  return get_attribute_double_value(dynamic_cast<Grid&> (*var).get_array(), attribute);
754  else
755  throw Error(malformed_expr,
756  string("No COARDS '") + attribute + "' attribute was found for the variable '" + var->name() + "'.");
757  }
758 
759  return string_to_double(remove_quotes(attribute_value).c_str());
760 }
761 
762 static double get_y_intercept(BaseType *var)
763 {
764  vector<string> attributes;
765  attributes.push_back("add_offset");
766  attributes.push_back("add_off");
767  return get_attribute_double_value(var, attributes);
768 }
769 
770 static double get_slope(BaseType *var)
771 {
772  return get_attribute_double_value(var, "scale_factor");
773 }
774 
775 static double get_missing_value(BaseType *var)
776 {
777  return get_attribute_double_value(var, "missing_value");
778 }
779 
792 void function_linear_scale(int argc, BaseType * argv[], DDS &, BaseType **btpp)
793 {
794  string
795  info =
796  string("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n")
797  + "<function name=\"linear_scale\" version=\"1.0b1\" href=\"http://docs.opendap.org/index.php/Server_Side_Processing_Functions#linear_scale\">\n"
798  + "</function>";
799 
800  if (argc == 0) {
801  Str *response = new Str("info");
802  response->set_value(info);
803  *btpp = response;
804  return;
805  }
806 
807  // Check for 1 or 3 arguments: 1 --> use attributes; 3 --> m & b supplied
808  DBG(cerr << "argc = " << argc << endl);
809  if (!(argc == 1 || argc == 3 || argc == 4))
810  throw Error(malformed_expr,
811  "Wrong number of arguments to linear_scale(). See linear_scale() for more information");
812 
813  // Get m & b
814  bool use_missing = false;
815  double m, b, missing = 0.0;
816  if (argc == 4) {
817  m = extract_double_value(argv[1]);
818  b = extract_double_value(argv[2]);
819  missing = extract_double_value(argv[3]);
820  use_missing = true;
821  }
822  else if (argc == 3) {
823  m = extract_double_value(argv[1]);
824  b = extract_double_value(argv[2]);
825  use_missing = false;
826  }
827  else {
828  m = get_slope(argv[0]);
829 
830  // This is really a hack; on a fair number of datasets, the y intercept
831  // is not given and is assumed to be 0. Here the function looks and
832  // catches the error if a y intercept is not found.
833  try {
834  b = get_y_intercept(argv[0]);
835  } catch (Error &e) {
836  b = 0.0;
837  }
838 
839  // This is not the best plan; the get_missing_value() function should
840  // do something other than throw, but to do that would require mayor
841  // surgery on get_attribute_double_value().
842  try {
843  missing = get_missing_value(argv[0]);
844  use_missing = true;
845  } catch (Error &e) {
846  use_missing = false;
847  }
848  }
849 
850  DBG(cerr << "m: " << m << ", b: " << b << endl);DBG(cerr << "use_missing: " << use_missing << ", missing: " << missing << endl);
851 
852  // Read the data, scale and return the result. Must replace the new data
853  // in a constructor (i.e., Array part of a Grid).
854  BaseType *dest = 0;
855  double *data;
856  if (argv[0]->type() == dods_grid_c) {
857  Array &source = *dynamic_cast<Grid&> (*argv[0]).get_array();
858  source.set_send_p(true);
859  source.read();
860  data = extract_double_array(&source);
861  int length = source.length();
862  int i = 0;
863  while (i < length) {
864  DBG2(cerr << "data[" << i << "]: " << data[i] << endl);
865  if (!use_missing || !double_eq(data[i], missing))
866  data[i] = data[i] * m + b;
867  DBG2(cerr << " >> data[" << i << "]: " << data[i] << endl);
868  ++i;
869  }
870 
871  // Vector::add_var will delete the existing 'template' variable
872  Float64 *temp_f = new Float64(source.name());
873  source.add_var(temp_f);
874 #ifdef VAL2BUF
875  source.val2buf(static_cast<void*>(data), false);
876 #else
877  source.set_value(data, i);
878 #endif
879  delete[] data; // val2buf copies.
880  delete temp_f; // add_var copies and then adds.
881  dest = argv[0];
882  }
883  else if (argv[0]->is_vector_type()) {
884  Array &source = dynamic_cast<Array&> (*argv[0]);
885  source.set_send_p(true);
886  // If the array is really a map, make sure to read using the Grid
887  // because of the HDF4 handler's odd behavior WRT dimensions.
888  if (source.get_parent() && source.get_parent()->type() == dods_grid_c)
889  source.get_parent()->read();
890  else
891  source.read();
892 
893  data = extract_double_array(&source);
894  int length = source.length();
895  int i = 0;
896  while (i < length) {
897  if (!use_missing || !double_eq(data[i], missing))
898  data[i] = data[i] * m + b;
899  ++i;
900  }
901 
902  Float64 *temp_f = new Float64(source.name());
903  source.add_var(temp_f);
904 
905  source.val2buf(static_cast<void*> (data), false);
906 
907  delete[] data; // val2buf copies.
908  delete temp_f; // add_var copies and then adds.
909 
910  dest = argv[0];
911  }
912  else if (argv[0]->is_simple_type() && !(argv[0]->type() == dods_str_c || argv[0]->type() == dods_url_c)) {
913  double data = extract_double_value(argv[0]);
914  if (!use_missing || !double_eq(data, missing))
915  data = data * m + b;
916 
917  dest = new Float64(argv[0]->name());
918 
919  dest->val2buf(static_cast<void*> (&data));
920 
921  }
922  else {
923  throw Error(malformed_expr, "The linear_scale() function works only for numeric Grids, Arrays and scalars.");
924  }
925 
926  *btpp = dest;
927  return;
928 }
929 
930 #if 0
931 
948 void
949 function_geoarray(int argc, BaseType * argv[], DDS &, BaseType **btpp)
950 {
951  string info =
952  string("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n") +
953  "<function name=\"geoarray\" version=\"0.9b1\" href=\"http://docs.opendap.org/index.php/Server_Side_Processing_Functions#geoarray\">\n" +
954  "</function>";
955 
956  if (argc == 0) {
957  Str *response = new Str("version");
958  response->set_value(info);
959  *btpp = response;
960  return;
961  }
962 
963  DBG(cerr << "argc = " << argc << endl);
964  if (!(argc == 5 || argc == 9 || argc == 11))
965  throw Error(malformed_expr,"Wrong number of arguments to geoarray(). See geoarray() for more information.");
966 
967  // Check the Array (and dup because the caller will free the variable).
968  Array *l_array = dynamic_cast < Array * >(argv[0]->ptr_duplicate());
969  if (!l_array)
970  throw Error(malformed_expr,"The first argument to geoarray() must be an Array variable!");
971 
972  try {
973 
974  // Read the bounding box and variable extents from the params
975  double bb_top = extract_double_value(argv[1]);
976  double bb_left = extract_double_value(argv[2]);
977  double bb_bottom = extract_double_value(argv[3]);
978  double bb_right = extract_double_value(argv[4]);
979 
980  switch (argc) {
981  case 5: {
982  ArrayGeoConstraint agc(l_array);
983 
984  agc.set_bounding_box(bb_left, bb_top, bb_right, bb_bottom);
985  // This also reads all of the data into the grid variable
986  agc.apply_constraint_to_data();
987  *btpp = agc.get_constrained_array();
988  return;
989  break;
990  }
991  case 9: {
992  double var_top = extract_double_value(argv[5]);
993  double var_left = extract_double_value(argv[6]);
994  double var_bottom = extract_double_value(argv[7]);
995  double var_right = extract_double_value(argv[8]);
996  ArrayGeoConstraint agc (l_array, var_left, var_top, var_right, var_bottom);
997 
998  agc.set_bounding_box(bb_left, bb_top, bb_right, bb_bottom);
999  // This also reads all of the data into the grid variable
1000  agc.apply_constraint_to_data();
1001  *btpp = agc.get_constrained_array();
1002  return;
1003  break;
1004  }
1005  case 11: {
1006  double var_top = extract_double_value(argv[5]);
1007  double var_left = extract_double_value(argv[6]);
1008  double var_bottom = extract_double_value(argv[7]);
1009  double var_right = extract_double_value(argv[8]);
1010  string projection = extract_string_argument(argv[9]);
1011  string datum = extract_string_argument(argv[10]);
1012  ArrayGeoConstraint agc(l_array,
1013  var_left, var_top, var_right, var_bottom,
1014  projection, datum);
1015 
1016  agc.set_bounding_box(bb_left, bb_top, bb_right, bb_bottom);
1017  // This also reads all of the data into the grid variable
1018  agc.apply_constraint_to_data();
1019  *btpp = agc.get_constrained_array();
1020  return;
1021  break;
1022  }
1023  default:
1024  throw InternalErr(__FILE__, __LINE__, "Wrong number of args to geoarray.");
1025  }
1026  }
1027  catch (Error & e) {
1028  throw e;
1029  }
1030  catch (exception & e) {
1031  throw
1032  InternalErr(string
1033  ("A C++ exception was thrown from inside geoarray(): ")
1034  + e.what());
1035 
1036  }
1037 
1038  throw InternalErr(__FILE__, __LINE__, "Impossible condition in geoarray.");
1039 }
1040 #endif
1041 
1043 {
1044  ce.add_function("grid", function_grid);
1045  ce.add_function("geogrid", function_geogrid);
1046  ce.add_function("linear_scale", function_linear_scale);
1047 #if 0
1048  ce.add_function("geoarray", function_geoarray);
1049 #endif
1050  ce.add_function("version", function_version);
1051 
1052  ce.add_function("dap", function_dap);
1053 }
1054 
1055 } // namespace libdap