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