USML
netcdf_coards.h
1 
5 #ifndef USML_NETCDF_COARDS_H
6 #define USML_NETCDF_COARDS_H
7 
8 #include <netcdfcpp.h>
9 #include <usml/ublas/ublas.h>
10 #include <usml/types/types.h>
11 
12 namespace usml {
13 namespace netcdf {
14 
15 using namespace usml::ublas ;
16 using namespace usml::types ;
17 
20 
36 template< class DATA_TYPE, int NUM_DIMS > class netcdf_coards :
37  public data_grid<DATA_TYPE,NUM_DIMS>
38 {
39  private:
40 
51  seq_vector* make_axis( NcFile& file, NcDim* dimension ) {
52 
53  // search for this axis in the NetCDF file
54 
55  NcVar* axis = file.get_var( dimension->name() ) ;
56  if ( axis == 0 ) {
57  throw std::invalid_argument("NetCDF variable not found");
58  }
59 
60  // search for linear or logrithmic pattern in this data
61 
62  const int N = (int) axis->num_vals() ;
63  bool isLinear = true ;
64  bool isLog = true ;
65  boost::numeric::ublas::vector<double> data(N) ;
66 
67  double p1 = axis->as_double(0) ; data(0) = p1 ;
68  double minValue = p1 ;
69  double maxValue = p1 ;
70 
71  if ( N > 1 ) {
72  double p2 = axis->as_double(1) ; data(1) = p2 ;
73  for ( int n=2 ; n < N ; ++n ) {
74  double p3 = axis->as_double(n) ; data(n) = p3 ;
75  maxValue = p3 ;
76  // printf("diff[%d] = %f\n",n,p3-p2);
77  if ( abs( ((p3-p2)-(p2-p1)) / p2 ) > 1E-4 ) {
78  isLinear = false ;
79  }
80  if ( p1 == 0.0 || p2 == 0.0 ||
81  abs( (p3/p2)-(p2/p1) ) > 1E-5 )
82  {
83  isLog = false ;
84  }
85  if ( ! ( isLinear || isLog ) ) break ;
86  p1 = p2 ;
87  p2 = p3 ;
88  }
89  }
90 
91  // build a new sequence for this axis
92 
93  cout << axis->name() << " N=" << N << " minValue=" << minValue << " maxValue=" << maxValue << endl ;
94  if ( isLinear ) {
95  return new seq_linear( minValue, (maxValue-minValue)/(N-1), N ) ;
96  } else if ( isLog ) {
97  return new seq_log( minValue, (maxValue/minValue)/(N-1), N ) ;
98  }
99  return new seq_data( data ) ;
100  }
101 
102  public:
103 
114  netcdf_coards( NcFile& file, NcToken name, bool read_fill=false ) {
115 
116  // search for this grid in the NetCDF file
117 
118  NcVar* variable = file.get_var( name ) ;
119  if ( variable == 0 ) {
120  throw std::invalid_argument("NetCDF variable not found");
121  }
122 
123  // read axis data from NetCDF file.
124 
125  size_t N = 1 ;
126  for ( int n=0 ; n < NUM_DIMS ; ++n ) {
127  seq_vector* ax = make_axis( file, variable->get_dim(n) ) ;
128  this->_axis[n] = ax ;
129  N *= this->_axis[n]->size() ;
130  }
131 
132  // extract missing attribute information
133 
134  NcError nc_error( NcError::silent_nonfatal ) ;
135 
136  DATA_TYPE missing = NAN ; // default value for missing data
137  NcAtt* att = variable->get_att("missing_value") ;
138  if ( att ) {
139  NcValues* values = att->values() ;
140  missing = (DATA_TYPE) values->as_double(0) ;
141  delete att ; delete values ;
142  }
143 
144  DATA_TYPE filling = NAN ; // default for fill value
145  if ( read_fill ) {
146  att = variable->get_att("_FillValue") ;
147  if ( att ) {
148  NcValues* values = att->values() ;
149  filling = values->as_double(0) ;
150  delete att ; delete values ;
151  }
152  }
153 
154  // copy interpolant data from the NetCDF file into local memory.
155  // replace missing data with fill value
156 
157  NcValues* values = variable->values() ;
158  this->_data = new DATA_TYPE[N] ;
159  for ( unsigned n=0 ; n < N ; ++n ) {
160  this->_data[n] = (DATA_TYPE) values->as_double(n) ;
161  if ( ! isnan(missing) ) {
162  if ( this->_data[n] == missing ) {
163  this->_data[n] = filling ;
164  }
165  }
166  }
167  delete values ;
168  }
169 } ;
170 
172 } // end of namespace netcdf
173 } // end of namespace usml
174 
175 #endif
Reads a single COARDS data grid from a netCDF file.
Definition: netcdf_coards.h:36
seq_vector * make_axis(NcFile &file, NcDim *dimension)
Construct a seq_vector from NetCDF dimension object.
Definition: netcdf_coards.h:51
BOOST_UBLAS_INLINE matrix_unary1_traits< E, scalar_abs< typename E::value_type > >::result_type abs(const matrix_expression< E > &e)
Magnitude of a complex matrix.
Definition: matrix_math.h:257
Sequence defined by an evenly spaced grid of points.
Definition: seq_linear.h:19
Sequence defined by a logarithmically spaced grid of points.
Definition: seq_log.h:20
netcdf_coards(NcFile &file, NcToken name, bool read_fill=false)
Extract a named data grid from an open NetCDF file.
Definition: netcdf_coards.h:114
Sequence defined by an unevenly spaced vector of points.
Definition: seq_data.h:28
size_type size() const
Returns the number of elements in this sequence.
Definition: seq_vector.h:138
N-dimensional data set and its associated axes.
Definition: data_grid.h:64
A read-only, monotonic sequence of values.
Definition: seq_vector.h:36