00001
00005 #pragma once
00006
00007 #include <netcdfcpp.h>
00008 #include <usml/ublas/ublas.h>
00009 #include <usml/types/types.h>
00010
00011 namespace usml {
00012 namespace netcdf {
00013
00014 using namespace usml::ublas ;
00015 using namespace usml::types ;
00016
00019
00035 template< class DATA_TYPE, int NUM_DIMS > class netcdf_coards :
00036 public data_grid<DATA_TYPE,NUM_DIMS>
00037 {
00038 private:
00039
00050 seq_vector* make_axis( NcFile& file, NcDim* dimension ) {
00051
00052
00053
00054 NcVar* axis = file.get_var( dimension->name() ) ;
00055 if ( axis == 0 ) {
00056 throw std::invalid_argument("NetCDF variable not found");
00057 }
00058
00059
00060
00061 const int N = (int) axis->num_vals() ;
00062 bool isLinear = true ;
00063 bool isLog = true ;
00064 boost::numeric::ublas::vector<double> data(N) ;
00065
00066 double p1 = axis->as_double(0) ; data(0) = p1 ;
00067 double minValue = p1 ;
00068 double maxValue = p1 ;
00069
00070 if ( N > 1 ) {
00071 double p2 = axis->as_double(1) ; data(1) = p2 ;
00072 for ( int n=2 ; n < N ; ++n ) {
00073 double p3 = axis->as_double(n) ; data(n) = p3 ;
00074 maxValue = p3 ;
00075
00076 if ( abs( ((p3-p2)-(p2-p1)) / p2 ) > 1E-4 ) {
00077 isLinear = false ;
00078 }
00079 if ( p1 == 0.0 || p2 == 0.0 ||
00080 abs( (p3/p2)-(p2/p1) ) > 1E-5 )
00081 {
00082 isLog = false ;
00083 }
00084 if ( ! ( isLinear || isLog ) ) break ;
00085 p1 = p2 ;
00086 p2 = p3 ;
00087 }
00088 }
00089
00090
00091
00092 cout << axis->name() << " N=" << N << " minValue=" << minValue << " maxValue=" << maxValue << endl ;
00093 if ( isLinear ) {
00094 return new seq_linear( minValue, (maxValue-minValue)/(N-1), N ) ;
00095 } else if ( isLog ) {
00096 return new seq_log( minValue, (maxValue/minValue)/(N-1), N ) ;
00097 }
00098 return new seq_data( data ) ;
00099 }
00100
00101 public:
00102
00113 netcdf_coards( NcFile& file, NcToken name, bool read_fill=false ) {
00114
00115
00116
00117 NcVar* variable = file.get_var( name ) ;
00118 if ( variable == 0 ) {
00119 throw std::invalid_argument("NetCDF variable not found");
00120 }
00121
00122
00123
00124 size_t N = 1 ;
00125 for ( int n=0 ; n < NUM_DIMS ; ++n ) {
00126 seq_vector* ax = make_axis( file, variable->get_dim(n) ) ;
00127 this->_axis[n] = ax ;
00128 N *= this->_axis[n]->size() ;
00129 }
00130
00131
00132
00133 NcError nc_error( NcError::silent_nonfatal ) ;
00134
00135 DATA_TYPE missing = NAN ;
00136 NcAtt* att = variable->get_att("missing_value") ;
00137 if ( att ) {
00138 NcValues* values = att->values() ;
00139 missing = (DATA_TYPE) values->as_double(0) ;
00140 delete att ; delete values ;
00141 }
00142
00143 DATA_TYPE filling = NAN ;
00144 if ( read_fill ) {
00145 att = variable->get_att("_FillValue") ;
00146 if ( att ) {
00147 NcValues* values = att->values() ;
00148 filling = values->as_double(0) ;
00149 delete att ; delete values ;
00150 }
00151 }
00152
00153
00154
00155
00156 NcValues* values = variable->values() ;
00157 this->_data = new DATA_TYPE[N] ;
00158 for ( size_t n=0 ; n < N ; ++n ) {
00159 this->_data[n] = (DATA_TYPE) values->as_double((long) n) ;
00160 if ( ! isnan(missing) ) {
00161 if ( this->_data[n] == missing ) {
00162 this->_data[n] = filling ;
00163 }
00164 }
00165 }
00166 delete values ;
00167 }
00168 } ;
00169
00171 }
00172 }