libdap++  Updated for version 3.11.7
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 25915 2012-10-24 00:14:58Z 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 
51 #undef FUNCTION_DAP // undef so the dap() function always returns an error;
52 // use keywords instead.
53 
54 #include "BaseType.h"
55 #include "Byte.h"
56 #include "Int16.h"
57 #include "UInt16.h"
58 #include "Int32.h"
59 #include "UInt32.h"
60 #include "Float32.h"
61 #include "Float64.h"
62 #include "Str.h"
63 #include "Url.h"
64 #include "Array.h"
65 #include "Structure.h"
66 #include "Sequence.h"
67 #include "Grid.h"
68 #include "Error.h"
69 #include "RValue.h"
70 
71 #include "GSEClause.h"
72 #include "GridGeoConstraint.h"
73 #include "ArrayGeoConstraint.h"
74 
75 #include "ce_functions.h"
76 #include "gse_parser.h"
77 #include "gse.tab.hh"
78 #include "debug.h"
79 #include "util.h"
80 
81 // We wrapped VC++ 6.x strtod() to account for a short coming
82 // in that function in regards to "NaN". I don't know if this
83 // still applies in more recent versions of that product.
84 // ROM - 12/2007
85 #ifdef WIN32
86 #include <limits>
87 double w32strtod(const char *, char **);
88 #endif
89 
90 using namespace std;
91 
92 int gse_parse(void *arg);
93 void gse_restart(FILE * in);
94 
95 // Glue routines declared in gse.lex
96 void gse_switch_to_buffer(void *new_buffer);
97 void gse_delete_buffer(void *buffer);
98 void *gse_string(const char *yy_str);
99 
100 namespace libdap {
101 
103 inline bool double_eq(double lhs, double rhs, double epsilon = 1.0e-5)
104 {
105  if (lhs > rhs)
106  return (lhs - rhs) < ((lhs + rhs) / epsilon);
107  else
108  return (rhs - lhs) < ((lhs + rhs) / epsilon);
109 }
110 
119 {
120  if (arg->type() != dods_str_c)
121  throw Error(malformed_expr, "The function requires a DAP string argument.");
122 
123  if (!arg->read_p())
124  throw InternalErr(__FILE__, __LINE__,
125  "The CE Evaluator built an argument list where some constants held no values.");
126 
127  string s = dynamic_cast<Str&> (*arg).value();
128 
129  DBG(cerr << "s: " << s << endl);
130 
131  return s;
132 }
133 
134 // @todo Replace new with vector<T> (vector<T> values(src_len);)
135 template<class T> static void set_array_using_double_helper(Array * a, double *src, int src_len)
136 {
137  T *values = new T[src_len];
138  for (int i = 0; i < src_len; ++i)
139  values[i] = (T) src[i];
140 
141 #ifdef VAL2BUF
142  a->val2buf(values, true);
143 #else
144  a->set_value(values, src_len);
145 #endif
146 
147  delete[] values;
148 }
149 
167 void set_array_using_double(Array * dest, double *src, int src_len)
168 {
169  // Simple types are Byte, ..., Float64, String and Url.
170  if ((dest->type() == dods_array_c && !dest->var()->is_simple_type()) || dest->var()->type() == dods_str_c
171  || dest->var()->type() == dods_url_c)
172  throw InternalErr(__FILE__, __LINE__, "The function requires a DAP numeric-type array argument.");
173 
174  // Test sizes. Note that Array::length() takes any constraint into account
175  // when it returns the length. Even if this was removed, the 'helper'
176  // function this uses calls Vector::val2buf() which uses Vector::width()
177  // which in turn uses length().
178  if (dest->length() != src_len)
179  throw InternalErr(
180  __FILE__,
181  __LINE__,
182  "The source and destination array sizes don't match (" + long_to_string(src_len) + " versus "
183  + long_to_string(dest->length()) + ").");
184 
185  // The types of arguments that the CE Parser will build for numeric
186  // constants are limited to Uint32, Int32 and Float64. See ce_expr.y.
187  // Expanded to work for any numeric type so it can be used for more than
188  // just arguments.
189  switch (dest->var()->type()) {
190  case dods_byte_c:
191  set_array_using_double_helper<dods_byte> (dest, src, src_len);
192  break;
193  case dods_uint16_c:
194  set_array_using_double_helper<dods_uint16> (dest, src, src_len);
195  break;
196  case dods_int16_c:
197  set_array_using_double_helper<dods_int16> (dest, src, src_len);
198  break;
199  case dods_uint32_c:
200  set_array_using_double_helper<dods_uint32> (dest, src, src_len);
201  break;
202  case dods_int32_c:
203  set_array_using_double_helper<dods_int32> (dest, src, src_len);
204  break;
205  case dods_float32_c:
206  set_array_using_double_helper<dods_float32> (dest, src, src_len);
207  break;
208  case dods_float64_c:
209  set_array_using_double_helper<dods_float64> (dest, src, src_len);
210  break;
211  default:
212  throw InternalErr(__FILE__, __LINE__,
213  "The argument list built by the CE parser contained an unsupported numeric type.");
214  }
215 
216  // Set the read_p property.
217  dest->set_read_p(true);
218 }
219 
220 template<class T> static double *extract_double_array_helper(Array * a)
221 {
222  int length = a->length();
223  // Could improve this using vector<T>. jhrg
224  T *b = new T[length];
225  a->value(b);
226 
227  double *dest = new double[length];
228  for (int i = 0; i < length; ++i)
229  dest[i] = (double) b[i];
230  delete[] b;
231 
232  return dest;
233 }
234 
240 {
241  // Simple types are Byte, ..., Float64, String and Url.
242  if ((a->type() == dods_array_c && !a->var()->is_simple_type()) || a->var()->type() == dods_str_c
243  || a->var()->type() == dods_url_c)
244  throw Error(malformed_expr, "The function requires a DAP numeric-type array argument.");
245 
246  if (!a->read_p())
247  throw InternalErr(__FILE__, __LINE__, string("The Array '") + a->name() + "'does not contain values.");
248 
249  // The types of arguments that the CE Parser will build for numeric
250  // constants are limited to Uint32, Int32 and Float64. See ce_expr.y.
251  // Expanded to work for any numeric type so it can be used for more than
252  // just arguments.
253  switch (a->var()->type()) {
254  case dods_byte_c:
255  return extract_double_array_helper<dods_byte> (a);
256  case dods_uint16_c:
257  return extract_double_array_helper<dods_uint16> (a);
258  case dods_int16_c:
259  return extract_double_array_helper<dods_int16> (a);
260  case dods_uint32_c:
261  return extract_double_array_helper<dods_uint32> (a);
262  case dods_int32_c:
263  return extract_double_array_helper<dods_int32> (a);
264  case dods_float32_c:
265  return extract_double_array_helper<dods_float32> (a);
266  case dods_float64_c:
267  return extract_double_array_helper<dods_float64> (a);
268  default:
269  throw InternalErr(__FILE__, __LINE__,
270  "The argument list built by the CE parser contained an unsupported numeric type.");
271  }
272 }
273 
282 {
283  // Simple types are Byte, ..., Float64, String and Url.
284  if (!arg->is_simple_type() || arg->type() == dods_str_c || arg->type() == dods_url_c)
285  throw Error(malformed_expr, "The function requires a DAP numeric-type argument.");
286 
287  if (!arg->read_p())
288  throw InternalErr(__FILE__, __LINE__,
289  "The CE Evaluator built an argument list where some constants held no values.");
290 
291  // The types of arguments that the CE Parser will build for numeric
292  // constants are limited to Uint32, Int32 and Float64. See ce_expr.y.
293  // Expanded to work for any numeric type so it can be used for more than
294  // just arguments.
295  switch (arg->type()) {
296  case dods_byte_c:
297  return (double) (dynamic_cast<Byte&> (*arg).value());
298  case dods_uint16_c:
299  return (double) (dynamic_cast<UInt16&> (*arg).value());
300  case dods_int16_c:
301  return (double) (dynamic_cast<Int16&> (*arg).value());
302  case dods_uint32_c:
303  return (double) (dynamic_cast<UInt32&> (*arg).value());
304  case dods_int32_c:
305  return (double) (dynamic_cast<Int32&> (*arg).value());
306  case dods_float32_c:
307  return (double) (dynamic_cast<Float32&> (*arg).value());
308  case dods_float64_c:
309  return dynamic_cast<Float64&> (*arg).value();
310  default:
311  throw InternalErr(__FILE__, __LINE__,
312  "The argument list built by the CE parser contained an unsupported numeric type.");
313  }
314 }
315 
322 void function_version(int, BaseType *[], DDS &, BaseType **btpp)
323 {
324  string
325  xml_value =
326  "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\
327  <functions>\
328  <function name=\"geogrid\" version=\"1.2\"/>\
329  <function name=\"grid\" version=\"1.0\"/>\
330  <function name=\"linear_scale\" version=\"1.0b1\"/>\
331  <function name=\"version\" version=\"1.0\"/>\
332  <function name=\"dap\" version=\"1.0\"/>\
333  </functions>";
334 
335  // <function name=\"geoarray\" version=\"0.9b1\"/>
336 
337  Str *response = new Str("version");
338 
339  response->set_value(xml_value);
340  *btpp = response;
341  return;
342 }
343 
345 {
346 #ifdef FUNCTION_DAP
347  if (argc != 1) {
348  throw Error("The 'dap' function must be called with a version number.\n\
349  see http://docs.opendap.org/index.php/Server_Side_Processing_Functions#dap");
350  }
351 
352  double pv = extract_double_value(argv[0]);
353  dds.set_dap_version(pv);
354 #else
355  throw Error(
356  "The 'dap' function is not supported in lieu of Constraint expression 'keywords.'\n\
357 see http://docs.opendap.org/index.php/Server_Side_Processing_Functions#keywords");
358 #endif
359 }
360 
361 static void parse_gse_expression(gse_arg * arg, BaseType * expr)
362 {
363  gse_restart(0); // Restart the scanner.
364  void *cls = gse_string(extract_string_argument(expr).c_str());
365  // gse_switch_to_buffer(cls); // Get set to scan the string.
366  bool status = gse_parse((void *) arg) == 0;
367  gse_delete_buffer(cls);
368  if (!status)
369  throw Error(malformed_expr, "Error parsing grid selection.");
370 }
371 
372 static void apply_grid_selection_expr(Grid * grid, GSEClause * clause)
373 {
374  // Basic plan: For each map, look at each clause and set start and stop
375  // to be the intersection of the ranges in those clauses.
376  Grid::Map_iter map_i = grid->map_begin();
377  while (map_i != grid->map_end() && (*map_i)->name() != clause->get_map_name())
378  ++map_i;
379 
380  if (map_i == grid->map_end())
381  throw Error(malformed_expr,
382  "The map vector '" + clause->get_map_name() + "' is not in the grid '" + grid->name() + "'.");
383 
384  // Use pointer arith & the rule that map order must match array dim order
385  Array::Dim_iter grid_dim = (grid->get_array()->dim_begin() + (map_i - grid->map_begin()));
386 
387  Array *map = dynamic_cast<Array *> ((*map_i));
388  if (!map)
389  throw InternalErr(__FILE__, __LINE__, "Expected an Array");
390  int start = max(map->dimension_start(map->dim_begin()), clause->get_start());
391  int stop = min(map->dimension_stop(map->dim_begin()), clause->get_stop());
392 
393  if (start > stop) {
394  ostringstream msg;
395  msg << "The expressions passed to grid() do not result in an inclusive \n" << "subset of '"
396  << clause->get_map_name() << "'. The map's values range " << "from " << clause->get_map_min_value()
397  << " to " << clause->get_map_max_value() << ".";
398  throw Error(malformed_expr, msg.str());
399  }
400 
401  DBG(cerr << "Setting constraint on " << map->name()
402  << "[" << start << ":" << stop << "]" << endl);
403 
404  // Stride is always one.
405  map->add_constraint(map->dim_begin(), start, 1, stop);
406  grid->get_array()->add_constraint(grid_dim, start, 1, stop);
407 }
408 
409 static void apply_grid_selection_expressions(Grid * grid, vector<GSEClause *> clauses)
410 {
411  vector<GSEClause *>::iterator clause_i = clauses.begin();
412  while (clause_i != clauses.end())
413  apply_grid_selection_expr(grid, *clause_i++);
414 
415  grid->set_read_p(false);
416 }
417 
427 void function_miic_ex2(int argc, BaseType * argv[], DDS &dds, BaseType **btpp)
428 {
429  Array *l_lat = 0;
430  Array *l_lon = 0;
431  switch (argc) {
432  case 0: {
433  // First find the latitude and longitude variables. This assumes the file is
434  // CF. Note that the names of these are not passed into the function so it
435  // looks up the variables in the DDS. If the names were passed in, this step
436  // would be skipped.
437  if (!(l_lat = dynamic_cast<Array *>(dds.var("Latitude")))
438  || !(l_lon = dynamic_cast<Array *>(dds.var("Latitude"))))
439  throw Error(malformed_expr, "Could not find the Latitude or Longitude data!");
440  break;
441  }
442  case 2: {
443  l_lat = dynamic_cast<Array*>(argv[0]);
444  l_lon = dynamic_cast<Array*>(argv[1]);
445  if (!l_lat || !l_lon)
446  throw Error(malformed_expr, "Expected two Array variables as arguments.");
447  break;
448  }
449  default:
450  throw Error(malformed_expr, "Expected either zero or two arguments.");
451  }
452 
453  // Now read the data values into C arrays the function can use. The length of the
454  // data is l_lat->length() and l_lon->length() resp. Use delete[] to release the
455  // storage. Also note that the Array* must be used to determine the number of
456  // dimensions of the arrays - extract_double_array() returns a simple vector.
457  l_lat->read();
458  double *lat = extract_double_array(l_lat);
459  l_lon->read();
460  double *lon = extract_double_array(l_lon);
461 
462  // For this example, make a Structure, add two new variables to it and stuff
463  // these values in them. Make the new variable one-dimensional arrays (vectors)
464  // just to keep the code simple.
465  Structure *dest = new Structure("MODIS_Geo_information");
466 
467  Array *new_lat = new Array("MODIS_Latitude", new Float64("MODIS_Latitude"));
468  new_lat->append_dim(l_lat->length());
469  new_lat->set_value(lat, l_lat->length());
470  dest->add_var(new_lat);
471 
472  Array *new_lon = new Array("MODIS_Longtude", new Float64("MODIS_Longtude"));
473  new_lon->append_dim(l_lon->length());
474  new_lon->set_value(lon, l_lon->length());
475  dest->add_var(new_lon);
476 
477  *btpp = dest;
478 }
479 
516 void function_grid(int argc, BaseType * argv[], DDS &, BaseType **btpp)
517 {
518  DBG(cerr << "Entering function_grid..." << endl);
519 
520  string
521  info =
522  string("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n")
523  + "<function name=\"grid\" version=\"1.0\" href=\"http://docs.opendap.org/index.php/Server_Side_Processing_Functions#grid\">\n"
524  + "</function>\n";
525 
526  if (argc == 0) {
527  Str *response = new Str("info");
528  response->set_value(info);
529  *btpp = response;
530  return;
531  }
532 
533  Grid *original_grid = dynamic_cast<Grid *> (argv[0]);
534  if (!original_grid)
535  throw Error(malformed_expr, "The first argument to grid() must be a Grid variable!");
536 
537  // Duplicate the grid; ResponseBuilder::send_data() will delete the variable
538  // after serializing it.
539  BaseType *btp = original_grid->ptr_duplicate();
540  Grid *l_grid = dynamic_cast<Grid *> (btp);
541  if (!l_grid) {
542  delete btp;
543  throw InternalErr(__FILE__, __LINE__, "Expected a Grid.");
544  }
545 
546  DBG(cerr << "grid: past initialization code" << endl);
547 
548  // Read the maps. Do this before calling parse_gse_expression(). Avoid
549  // reading the array until the constraints have been applied because it
550  // might be really large.
551 
552  // This version makes sure to set the send_p flags which is needed for
553  // the hdf4 handler (and is what should be done in general).
554  Grid::Map_iter i = l_grid->map_begin();
555  while (i != l_grid->map_end())
556  (*i++)->set_send_p(true);
557 
558  l_grid->read();
559 
560  DBG(cerr << "grid: past map read" << endl);
561 
562  // argv[1..n] holds strings; each are little expressions to be parsed.
563  // When each expression is parsed, the parser makes a new instance of
564  // GSEClause. GSEClause checks to make sure the named map really exists
565  // in the Grid and that the range of values given makes sense.
566  vector<GSEClause *> clauses;
567  gse_arg *arg = new gse_arg(l_grid);
568  for (int i = 1; i < argc; ++i) {
569  parse_gse_expression(arg, argv[i]);
570  clauses.push_back(arg->get_gsec());
571  }
572  delete arg;
573  arg = 0;
574 
575  apply_grid_selection_expressions(l_grid, clauses);
576 
577  DBG(cerr << "grid: past gse application" << endl);
578 
579  l_grid->get_array()->set_send_p(true);
580 
581  l_grid->read();
582 
583  // Make a new grid here and copy just the parts of the Grid
584  // that are in the current projection - this means reading
585  // the array slicing information, extracting the correct
586  // values and building destination arrays with just those
587  // values.
588 
589  *btpp = l_grid;
590  return;
591 }
592 
627 void function_geogrid(int argc, BaseType * argv[], DDS &, BaseType **btpp)
628 {
629  string
630  info =
631  string("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n")
632  + "<function name=\"geogrid\" version=\"1.2\" href=\"http://docs.opendap.org/index.php/Server_Side_Processing_Functions#geogrid\">\n"
633  + "</function>";
634 
635  if (argc == 0) {
636  Str *response = new Str("version");
637  response->set_value(info);
638  *btpp = response;
639  return;
640  }
641 
642  // There are two main forms of this function, one that takes a Grid and one
643  // that takes a Grid and two Arrays. The latter provides a way to explicitly
644  // tell the function which maps contain lat and lon data. The remaining
645  // arguments are the same for both versions, although that includes a
646  // varying argument list.
647 
648  // Look at the types of the first three arguments to determine which of the
649  // two forms were used to call this function.
650  Grid *l_grid = 0;
651  if (argc < 1 || !(l_grid = dynamic_cast<Grid *> (argv[0]->ptr_duplicate())))
652  throw Error(malformed_expr, "The first argument to geogrid() must be a Grid variable!");
653 
654  // Both forms require at least this many args
655  if (argc < 5)
656  throw Error(malformed_expr,
657  "Wrong number of arguments to geogrid() (expected at least 5 args). See geogrid() for more information.");
658 
659  bool grid_lat_lon_form;
660  Array *l_lat = 0;
661  Array *l_lon = 0;
662  if (!(l_lat = dynamic_cast<Array *> (argv[1]))) //->ptr_duplicate())))
663  grid_lat_lon_form = false;
664  else if (!(l_lon = dynamic_cast<Array *> (argv[2]))) //->ptr_duplicate())))
665  throw Error(malformed_expr,
666  "When using the Grid, Lat, Lon form of geogrid() both the lat and lon maps must be given (lon map missing)!");
667  else
668  grid_lat_lon_form = true;
669 
670  if (grid_lat_lon_form && argc < 7)
671  throw Error(malformed_expr,
672  "Wrong number of arguments to geogrid() (expected at least 7 args). See geogrid() for more information.");
673 
674 #if 0
675  Grid *l_grid = dynamic_cast < Grid * >(argv[0]->ptr_duplicate());
676  if (!l_grid)
677  throw Error(malformed_expr,"The first argument to geogrid() must be a Grid variable!");
678 #endif
679  // Read the maps. Do this before calling parse_gse_expression(). Avoid
680  // reading the array until the constraints have been applied because it
681  // might be really large.
682  //
683  // Trick: Some handlers build Grids from a combination of Array
684  // variables and attributes. Those handlers (e.g., hdf4) use the send_p
685  // property to determine which parts of the Grid to read *but they can
686  // only read the maps from within Grid::read(), not the map's read()*.
687  // Since the Grid's array does not have send_p set, it will not be read
688  // by the call below to Grid::read().
689  Grid::Map_iter i = l_grid->map_begin();
690  while (i != l_grid->map_end())
691  (*i++)->set_send_p(true);
692 
693  l_grid->read();
694  // Calling read() above sets the read_p flag for the entire grid; clear it
695  // for the grid's array so that later on the code will be sure to read it
696  // under all circumstances.
697  l_grid->get_array()->set_read_p(false);
698  DBG(cerr << "geogrid: past map read" << endl);
699 
700  // Look for Grid Selection Expressions tacked onto the end of the BB
701  // specification. If there are any, evaluate them before evaluating the BB.
702  int min_arg_count = (grid_lat_lon_form) ? 7 : 5;
703  if (argc > min_arg_count) {
704  // argv[5..n] holds strings; each are little Grid Selection Expressions
705  // to be parsed and evaluated.
706  vector<GSEClause *> clauses;
707  gse_arg *arg = new gse_arg(l_grid);
708  for (int i = min_arg_count; i < argc; ++i) {
709  parse_gse_expression(arg, argv[i]);
710  clauses.push_back(arg->get_gsec());
711  }
712  delete arg;
713  arg = 0;
714 
715  apply_grid_selection_expressions(l_grid, clauses);
716  }
717 
718  try {
719  // Build a GeoConstraint object. If there are no longitude/latitude
720  // maps then this constructor throws Error.
721  GridGeoConstraint gc(l_grid);
722 
723  // This sets the bounding box and modifies the maps to match the
724  // notation of the box (0/359 or -180/179)
725  int box_index_offset = (grid_lat_lon_form) ? 3 : 1;
726  double top = extract_double_value(argv[box_index_offset]);
727  double left = extract_double_value(argv[box_index_offset + 1]);
728  double bottom = extract_double_value(argv[box_index_offset + 2]);
729  double right = extract_double_value(argv[box_index_offset + 3]);
730  gc.set_bounding_box(top, left, bottom, right);
731  DBG(cerr << "geogrid: past bounding box set" << endl);
732 
733  // This also reads all of the data into the grid variable
735  DBG(cerr << "geogrid: past apply constraint" << endl);
736 
737  // In this function the l_grid pointer is the same as the pointer returned
738  // by this call. The caller of the function must free the pointer.
739  *btpp = gc.get_constrained_grid();
740  return;
741  } catch (Error &e) {
742  throw e;
743  } catch (exception & e) {
744  throw InternalErr(string("A C++ exception was thrown from inside geogrid(): ") + e.what());
745  }
746 }
747 
748 // These static functions could be moved to a class that provides a more
749 // general interface for COARDS/CF someday. Assume each BaseType comes bundled
750 // with an attribute table.
751 
752 // This was ripped from parser-util.cc
753 static double string_to_double(const char *val)
754 {
755  char *ptr;
756  errno = 0;
757  // Clear previous value. 5/21/2001 jhrg
758 
759 #ifdef WIN32
760  double v = w32strtod(val, &ptr);
761 #else
762  double v = strtod(val, &ptr);
763 #endif
764 
765  if ((v == 0.0 && (val == ptr || errno == HUGE_VAL || errno == ERANGE)) || *ptr != '\0') {
766  throw Error(malformed_expr, string("Could not convert the string '") + val + "' to a double.");
767  }
768 
769  double abs_val = fabs(v);
770  if (abs_val > DODS_DBL_MAX || (abs_val != 0.0 && abs_val < DODS_DBL_MIN))
771  throw Error(malformed_expr, string("Could not convert the string '") + val + "' to a double.");
772 
773  return v;
774 }
775 
785 static double get_attribute_double_value(BaseType *var, vector<string> &attributes)
786 {
787  // This code also builds a list of the attribute values that have been
788  // passed in but not found so that an informative message can be returned.
789  AttrTable &attr = var->get_attr_table();
790  string attribute_value = "";
791  string values = "";
792  vector<string>::iterator i = attributes.begin();
793  while (attribute_value == "" && i != attributes.end()) {
794  values += *i;
795  if (!values.empty())
796  values += ", ";
797  attribute_value = attr.get_attr(*i++);
798  }
799 
800  // If the value string is empty, then look at the grid's array (if it's a
801  // grid) or throw an Error.
802  if (attribute_value.empty()) {
803  if (var->type() == dods_grid_c)
804  return get_attribute_double_value(dynamic_cast<Grid&> (*var).get_array(), attributes);
805  else
806  throw Error(
808  string("No COARDS/CF '") + values.substr(0, values.length() - 2)
809  + "' attribute was found for the variable '" + var->name() + "'.");
810  }
811 
812  return string_to_double(remove_quotes(attribute_value).c_str());
813 }
814 
815 static double get_attribute_double_value(BaseType *var, const string &attribute)
816 {
817  AttrTable &attr = var->get_attr_table();
818  string attribute_value = attr.get_attr(attribute);
819 
820  // If the value string is empty, then look at the grid's array (if it's a
821  // grid or throw an Error.
822  if (attribute_value.empty()) {
823  if (var->type() == dods_grid_c)
824  return get_attribute_double_value(dynamic_cast<Grid&> (*var).get_array(), attribute);
825  else
826  throw Error(malformed_expr,
827  string("No COARDS '") + attribute + "' attribute was found for the variable '" + var->name() + "'.");
828  }
829 
830  return string_to_double(remove_quotes(attribute_value).c_str());
831 }
832 
833 static double get_y_intercept(BaseType *var)
834 {
835  vector<string> attributes;
836  attributes.push_back("add_offset");
837  attributes.push_back("add_off");
838  return get_attribute_double_value(var, attributes);
839 }
840 
841 static double get_slope(BaseType *var)
842 {
843  return get_attribute_double_value(var, "scale_factor");
844 }
845 
846 static double get_missing_value(BaseType *var)
847 {
848  return get_attribute_double_value(var, "missing_value");
849 }
850 
863 void function_linear_scale(int argc, BaseType * argv[], DDS &, BaseType **btpp)
864 {
865  string
866  info =
867  string("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n")
868  + "<function name=\"linear_scale\" version=\"1.0b1\" href=\"http://docs.opendap.org/index.php/Server_Side_Processing_Functions#linear_scale\">\n"
869  + "</function>";
870 
871  if (argc == 0) {
872  Str *response = new Str("info");
873  response->set_value(info);
874  *btpp = response;
875  return;
876  }
877 
878  // Check for 1 or 3 arguments: 1 --> use attributes; 3 --> m & b supplied
879  DBG(cerr << "argc = " << argc << endl);
880  if (!(argc == 1 || argc == 3 || argc == 4))
881  throw Error(malformed_expr,
882  "Wrong number of arguments to linear_scale(). See linear_scale() for more information");
883 
884  // Get m & b
885  bool use_missing = false;
886  double m, b, missing = 0.0;
887  if (argc == 4) {
888  m = extract_double_value(argv[1]);
889  b = extract_double_value(argv[2]);
890  missing = extract_double_value(argv[3]);
891  use_missing = true;
892  }
893  else if (argc == 3) {
894  m = extract_double_value(argv[1]);
895  b = extract_double_value(argv[2]);
896  use_missing = false;
897  }
898  else {
899  m = get_slope(argv[0]);
900 
901  // This is really a hack; on a fair number of datasets, the y intercept
902  // is not given and is assumed to be 0. Here the function looks and
903  // catches the error if a y intercept is not found.
904  try {
905  b = get_y_intercept(argv[0]);
906  } catch (Error &e) {
907  b = 0.0;
908  }
909 
910  // This is not the best plan; the get_missing_value() function should
911  // do something other than throw, but to do that would require mayor
912  // surgery on get_attribute_double_value().
913  try {
914  missing = get_missing_value(argv[0]);
915  use_missing = true;
916  } catch (Error &e) {
917  use_missing = false;
918  }
919  }
920 
921  DBG(cerr << "m: " << m << ", b: " << b << endl);DBG(cerr << "use_missing: " << use_missing << ", missing: " << missing << endl);
922 
923  // Read the data, scale and return the result. Must replace the new data
924  // in a constructor (i.e., Array part of a Grid).
925  BaseType *dest = 0;
926  double *data;
927  if (argv[0]->type() == dods_grid_c) {
928 #if 0
929  // For a Grid, the function scales only the Array part.
930  Array *source = dynamic_cast<Grid*>(argv[0])->get_array();
931  //argv[0]->set_send_p(true);
932  //source->set_send_p(true);
933  source->read();
934  data = extract_double_array(source);
935  int length = source->length();
936  for (int i = 0; i < length; ++i)
937  data[i] = data[i] * m + b;
938 #if 0
939  int i = 0;
940  while (i < length) {
941  DBG2(cerr << "data[" << i << "]: " << data[i] << endl);
942  if (!use_missing || !double_eq(data[i], missing))
943  data[i] = data[i] * m + b;
944  DBG2(cerr << " >> data[" << i << "]: " << data[i] << endl);
945  ++i;
946  }
947 #endif
948  // Vector::add_var will delete the existing 'template' variable
949  Float64 *temp_f = new Float64(source->name());
950  source->add_var(temp_f);
951 
952 #ifdef VAL2BUF
953  source.val2buf(static_cast<void*>(data), false);
954 #else
955  source->set_value(data, length);
956 #endif
957  delete[] data; // val2buf copies.
958  delete temp_f; // add_var copies and then adds.
959  dest = argv[0];
960  dest->set_send_p(true);
961 #endif
962  // Grab the whole Grid; note that the scaling is done only on the array part
963  Grid &source = dynamic_cast<Grid&>(*argv[0]);
964 
965  DBG(cerr << "Grid send_p: " << source.send_p() << endl);
966  DBG(cerr << "Grid Array send_p: " << source.get_array()->send_p() << endl);
967 
968  // Read the grid; set send_p since Grid is a kind of constructor and
969  // read will only be called on it's fields if their send_p flag is set
970  source.set_send_p(true);
971  source.read();
972 
973  // Get the Array part and read the values
974  Array *a = source.get_array();
975  //a->read();
976  data = extract_double_array(a);
977 
978  // Now scale the data.
979  int length = a->length();
980  for (int i = 0; i < length; ++i)
981  data[i] = data[i] * m + b;
982 #if 0
983  // read the maps so that those values will be copied when the source Grid
984  // is copied to the dest Grid
985  Grid::Map_iter s = source.map_begin();
986  while (s != source.map_end()) {
987  static_cast<Array*>(*s)->read();
988  ++s;
989  }
990 #endif
991  // Copy source Grid to result Grid. Could improve on this by not using this
992  // trick since it copies all of 'source' to 'dest', including the main Array.
993  // The next bit of code will replace those values with the newly scaled ones.
994  Grid *result = new Grid(source);
995 
996  // Now load the transferred values; use Float64 as the new type of the result
997  // Grid Array.
998  result->get_array()->add_var_nocopy(new Float64(source.name()));
999  result->get_array()->set_value(data, length);
1000  delete[] data;
1001 
1002 #if 0
1003  // Now set the maps (NB: the copy constructor does not copy data)
1004  Grid::Map_iter s = source.map_begin();
1005  Grid::Map_iter d = result->map_begin();
1006  while (s != source.map_end()) {
1007  Array *a = static_cast<Array*>(*s);
1008  a->read();
1009  switch(a->var()->type()) {
1010  case dods_byte_c: {
1011  vector<dods_byte> v(a->length());
1012  a->value(&v[0]);
1013  static_cast<Array*>(*d)->set_value(v, v.size());
1014  break;
1015  }
1016  case dods_float32_c: {
1017  vector<dods_float32> v(a->length());
1018  a->value(&v[0]);
1019  static_cast<Array*>(*d)->set_value(v, a->length());
1020  break;
1021  }
1022  default:
1023  throw Error("Non-numeric Grid Map not supported by linear_scale().");
1024  }
1025  ++s; ++d;
1026  }
1027 #endif
1028 
1029  // FIXME result->set_send_p(true);
1030  DBG(cerr << "Grid send_p: " << result->send_p() << endl);
1031  DBG(cerr << "Grid Array send_p: " << result->get_array()->send_p() << endl);
1032 
1033  dest = result;
1034  }
1035  else if (argv[0]->is_vector_type()) {
1036 #if 0
1037  Array &source = dynamic_cast<Array&> (*argv[0]);
1038  source.set_send_p(true);
1039  // If the array is really a map, make sure to read using the Grid
1040  // because of the HDF4 handler's odd behavior WRT dimensions.
1041  if (source.get_parent() && source.get_parent()->type() == dods_grid_c)
1042  source.get_parent()->read();
1043  else
1044  source.read();
1045 
1046  data = extract_double_array(&source);
1047  int length = source.length();
1048  int i = 0;
1049  while (i < length) {
1050  if (!use_missing || !double_eq(data[i], missing))
1051  data[i] = data[i] * m + b;
1052  ++i;
1053  }
1054 
1055  Float64 *temp_f = new Float64(source.name());
1056  source.add_var(temp_f);
1057 
1058  source.val2buf(static_cast<void*> (data), false);
1059 
1060  delete[] data; // val2buf copies.
1061  delete temp_f; // add_var copies and then adds.
1062 
1063  dest = argv[0];
1064 #endif
1065  Array &source = dynamic_cast<Array&>(*argv[0]);
1066  // If the array is really a map, make sure to read using the Grid
1067  // because of the HDF4 handler's odd behavior WRT dimensions.
1068  if (source.get_parent() && source.get_parent()->type() == dods_grid_c) {
1069  source.get_parent()->set_send_p(true);
1070  source.get_parent()->read();
1071  }
1072  else
1073  source.read();
1074 
1075  data = extract_double_array(&source);
1076  int length = source.length();
1077  for (int i = 0; i < length; ++i)
1078  data[i] = data[i] * m + b;
1079 
1080  Array *result = new Array(source);
1081 
1082  result->add_var_nocopy(new Float64(source.name()));
1083  result->set_value(data, length);
1084 
1085  delete[] data; // val2buf copies.
1086 
1087  dest = result;
1088  }
1089  else if (argv[0]->is_simple_type() && !(argv[0]->type() == dods_str_c || argv[0]->type() == dods_url_c)) {
1090  double data = extract_double_value(argv[0]);
1091  if (!use_missing || !double_eq(data, missing))
1092  data = data * m + b;
1093 
1094  Float64 *fdest = new Float64(argv[0]->name());
1095 
1096  fdest->set_value(data);
1097  // dest->val2buf(static_cast<void*> (&data));
1098  dest = fdest;
1099  }
1100  else {
1101  throw Error(malformed_expr, "The linear_scale() function works only for numeric Grids, Arrays and scalars.");
1102  }
1103 
1104  *btpp = dest;
1105  return;
1106 }
1107 
1108 #if 0
1109 
1126 void
1127 function_geoarray(int argc, BaseType * argv[], DDS &, BaseType **btpp)
1128 {
1129  string info =
1130  string("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n") +
1131  "<function name=\"geoarray\" version=\"0.9b1\" href=\"http://docs.opendap.org/index.php/Server_Side_Processing_Functions#geoarray\">\n" +
1132  "</function>";
1133 
1134  if (argc == 0) {
1135  Str *response = new Str("version");
1136  response->set_value(info);
1137  *btpp = response;
1138  return;
1139  }
1140 
1141  DBG(cerr << "argc = " << argc << endl);
1142  if (!(argc == 5 || argc == 9 || argc == 11))
1143  throw Error(malformed_expr,"Wrong number of arguments to geoarray(). See geoarray() for more information.");
1144 
1145  // Check the Array (and dup because the caller will free the variable).
1146  Array *l_array = dynamic_cast < Array * >(argv[0]->ptr_duplicate());
1147  if (!l_array)
1148  throw Error(malformed_expr,"The first argument to geoarray() must be an Array variable!");
1149 
1150  try {
1151 
1152  // Read the bounding box and variable extents from the params
1153  double bb_top = extract_double_value(argv[1]);
1154  double bb_left = extract_double_value(argv[2]);
1155  double bb_bottom = extract_double_value(argv[3]);
1156  double bb_right = extract_double_value(argv[4]);
1157 
1158  switch (argc) {
1159  case 5: {
1160  ArrayGeoConstraint agc(l_array);
1161 
1162  agc.set_bounding_box(bb_left, bb_top, bb_right, bb_bottom);
1163  // This also reads all of the data into the grid variable
1164  agc.apply_constraint_to_data();
1165  *btpp = agc.get_constrained_array();
1166  return;
1167  break;
1168  }
1169  case 9: {
1170  double var_top = extract_double_value(argv[5]);
1171  double var_left = extract_double_value(argv[6]);
1172  double var_bottom = extract_double_value(argv[7]);
1173  double var_right = extract_double_value(argv[8]);
1174  ArrayGeoConstraint agc (l_array, var_left, var_top, var_right, var_bottom);
1175 
1176  agc.set_bounding_box(bb_left, bb_top, bb_right, bb_bottom);
1177  // This also reads all of the data into the grid variable
1178  agc.apply_constraint_to_data();
1179  *btpp = agc.get_constrained_array();
1180  return;
1181  break;
1182  }
1183  case 11: {
1184  double var_top = extract_double_value(argv[5]);
1185  double var_left = extract_double_value(argv[6]);
1186  double var_bottom = extract_double_value(argv[7]);
1187  double var_right = extract_double_value(argv[8]);
1188  string projection = extract_string_argument(argv[9]);
1189  string datum = extract_string_argument(argv[10]);
1190  ArrayGeoConstraint agc(l_array,
1191  var_left, var_top, var_right, var_bottom,
1192  projection, datum);
1193 
1194  agc.set_bounding_box(bb_left, bb_top, bb_right, bb_bottom);
1195  // This also reads all of the data into the grid variable
1196  agc.apply_constraint_to_data();
1197  *btpp = agc.get_constrained_array();
1198  return;
1199  break;
1200  }
1201  default:
1202  throw InternalErr(__FILE__, __LINE__, "Wrong number of args to geoarray.");
1203  }
1204  }
1205  catch (Error & e) {
1206  throw e;
1207  }
1208  catch (exception & e) {
1209  throw
1210  InternalErr(string
1211  ("A C++ exception was thrown from inside geoarray(): ")
1212  + e.what());
1213 
1214  }
1215 
1216  throw InternalErr(__FILE__, __LINE__, "Impossible condition in geoarray.");
1217 }
1218 #endif
1219 
1221 {
1222  ce.add_function("grid", function_grid);
1223  ce.add_function("geogrid", function_geogrid);
1224  ce.add_function("linear_scale", function_linear_scale);
1225  ce.add_function("version", function_version);
1226 
1227  ce.add_function("miic_ex2", function_miic_ex2);
1228 
1229  ce.add_function("dap", function_dap);
1230 }
1231 
1232 } // namespace libdap