First (non working) implementation of the ordinary kriging. The variogram fit...

First (non working) implementation of the ordinary kriging. The variogram fit still has to be done... But a few const issues have been fixed in Matrix and a new method added (to read a 1x1 matrix as a scalar).

Fixed some doxygen warnings, added a howto for writing spatial interpolations.

Better Oracle OCCI detection, renammed into FindOCCI.cmake and used by the ImisIO plugin compilation.

meteoio/FilterAlgorithms.cc

@@ -516,7 +516,7 @@ double FilterAlgorithms::ExpSmoothingAlgorithm(const std::vector<MeteoData>& vec
 
 /**
  * @brief Actual implementation of the smoothing algorithm
- * @param vecMeteo A vector of MeteoData
+ * @param vecMyMeteo A vector of MeteoData
  * @param paramindex The meteo data parameter to be smoothed
  * @return The wma smoothing value for the specified meteo parameter of the last element in vecMeteo
  */

meteoio/FilterAlgorithms.h

@@ -76,30 +76,31 @@ class FilterAlgorithms {
 		static void ExpSmoothingFilter(const std::vector<MeteoData>& vecM, const std::vector<StationData>& vecS, 
                          const std::vector<std::string>& vecArgs, const MeteoData::Parameters& paramindex,
                          std::vector<MeteoData>& vecWindowM, std::vector<StationData>& vecWindowS);
-		static void WMASmoothingFilter(const std::vector<MeteoData>& vecM, const std::vector<StationData>& vecS, 
-					const std::vector<std::string>& vecArgs, const MeteoData::Parameters& paramindex,
+		static void WMASmoothingFilter(const std::vector<MeteoData>& vecM, const std::vector<StationData>& vecS,
+		         const std::vector<std::string>& vecArgs, const MeteoData::Parameters& paramindex,
                          std::vector<MeteoData>& vecWindowM, std::vector<StationData>& vecWindowS);
 
-		static double ExpSmoothingAlgorithm(const std::vector<MeteoData>& vecMeteo, 
-									 const unsigned int& paramindex, const double& alpha);
-		static double WMASmoothingAlgorithm(const std::vector<MeteoData>& vecMeteo, const unsigned int& paramindex);
-
  	private:
 		static bool compareMeteoData (const MeteoData& m1, const MeteoData& m2);
 		static void parseFilterArguments(const std::string& filtername, const std::vector<std::string>& vecArgs_in,
-							const unsigned int& minArgs, const unsigned int& maxArgs,
-							bool& isSoft, std::vector<double>& vecArgs_out);
+		                    const unsigned int& minArgs, const unsigned int& maxArgs,
+		                    bool& isSoft, std::vector<double>& vecArgs_out);
 		static void parseWindowFilterArguments(const std::string& filtername, const std::vector<std::string>& vecArgs_in,
-								const unsigned int& minArgs, const unsigned int& maxArgs,
-								bool& isSoft, std::string& windowposition, std::vector<double>& vecArgs_out);
+		                    const unsigned int& minArgs, const unsigned int& maxArgs,
+		                    bool& isSoft, std::string& windowposition, std::vector<double>& vecArgs_out);
 		static unsigned int getWindowData(const std::string& filtername, const std::vector<MeteoData>& vecM,
                                     const unsigned int& pos, 
                                     const std::vector<std::string>& _vecArgs, std::vector<MeteoData>& vecResult);
 		static bool getWindowData(const std::string& filtername, const std::vector<MeteoData>& vecM,
-						const unsigned int& pos,
-						const Date& date, const std::vector<std::string>& _vecArgs,
-						const unsigned int& paramindex, std::vector<double>& vecWindow,
-						std::vector<Date> *vecDate = NULL);
+		                          const unsigned int& pos,
+		                          const Date& date, const std::vector<std::string>& _vecArgs,
+		                          const unsigned int& paramindex, std::vector<double>& vecWindow,
+		                          std::vector<Date> *vecDate = NULL);
+
+		static double ExpSmoothingAlgorithm(const std::vector<MeteoData>& vecMeteo, 
+		         const unsigned int& paramindex, const double& alpha);
+		static double WMASmoothingAlgorithm(const std::vector<MeteoData>& vecMyMeteo, const unsigned int& paramindex);
+
 		static std::map<std::string, FilterProperties> filterMap;
 		static const bool __init;    ///<helper variable to enable the init of static collection data
 		static bool initStaticData();///<initialize the static map filterMap

meteoio/InterpolationAlgorithms.cc

@@ -42,6 +42,7 @@ bool AlgorithmFactory::initStaticData()
 	setAlgorithms.insert("RH");        // relative humidity interpolation
 	setAlgorithms.insert("WIND_CURV"); // wind velocity interpolation (using a heuristic terrain effect)
 	setAlgorithms.insert("HNW_SNOW"); // precipitation interpolation according to (Magnusson, 2010)
+	setAlgorithms.insert("ODKRIG"); // ordinary kriging
 	setAlgorithms.insert("USER"); // read user provided grid
 
 	return true;
@@ -76,9 +77,9 @@ InterpolationAlgorithm* AlgorithmFactory::getAlgorithm(const std::string& _algon
 		return new RHAlgorithm(_mi, _dem, _vecMeteo, _vecStation, _vecArgs, _algoname);
 	} else if (algoname == "WIND_CURV"){
 		return new SimpleWindInterpolationAlgorithm(_mi, _dem, _vecMeteo, _vecStation, _vecArgs, _algoname);
-	/*} else if (algoname == "ODKRIG"){
+	} else if (algoname == "ODKRIG"){
 		return new OrdinaryKrigingAlgorithm(_mi, _dem, _vecMeteo, _vecStation, _vecArgs, _algoname);
-	*/} else if (algoname == "USER"){
+	} else if (algoname == "USER"){
 		return new USERInterpolation(_mi, _dem, _vecMeteo, _vecStation, _vecArgs, _algoname);
 	} else if (algoname == "HNW_SNOW"){
 		return new SnowHNWInterpolation(_mi, _dem, _vecMeteo, _vecStation, _vecArgs, _algoname);
@@ -565,5 +566,37 @@ void SnowHNWInterpolation::calculate(const MeteoData::Parameters& param, Grid2DO
 	grid.grid2D *= orig_mean/new_mean;
 }
 
+double OrdinaryKrigingAlgorithm::getQualityRating(const MeteoData::Parameters& param)
+{
+	throw IOException("ODKRIG interpolation algorithm not yet implemented...", AT);
+	vector<double> vecData;
+	unsigned int nrOfMeasurments = getData(param, vecData);
+
+	if(nrOfMeasurments>=20) return 0.9;
+	return 0.;
+}
+
+
+void OrdinaryKrigingAlgorithm::calculate(const MeteoData::Parameters& param, Grid2DObject& grid)
+{
+//optimization: getrange (from variogram fit -> exclude stations that are at distances > range (-> smaller matrix)
+	throw IOException("ODKRIG interpolation algorithm not yet implemented...", AT);
+
+	vector<double> vecData;
+	vector<StationData> vecMeta;
+	getData(param, vecData, vecMeta);
+
+	Interpol2D::ODKriging(vecData, vecMeta, dem, grid);
+	
+}
+
+double OrdinaryKrigingAlgorithm::computeVariogram(const std::vector<StationData>& /*vecStations*/) const
+{
+	//return variogramm fit of covariance between stations i and j
+	//HACK: todo!
+	
+	return 1.;
+}
+
 } //namespace
 

meteoio/InterpolationAlgorithms.h

@@ -74,6 +74,7 @@ class Meteo2DInterpolator; // forward declaration, cyclic header include
  * - RH: the dew point temperatures are interpolated using IDW_LAPSE, then reconverted locally to relative humidity (see RHAlgorithm)
  * - WIND_CURV: the wind field (VW and DW) is interpolated using IDW_LAPSE and then altered depending on the local curvature and slope (taken from the DEM, see SimpleWindInterpolationAlgorithm)
  * - HNW_SNOW: precipitation interpolation according to (Magnusson, 2010) (see SnowHNWInterpolation)
+ * - ODKRIG: ordinary kriging THIS IS NOT YET USABLE
  * - USER: user provided grids to be read from disk (if available, see USERinterpolation) THIS IS NOT YET USABLE
  *
  * @section lapse Lapse rates
@@ -108,6 +109,7 @@ class Meteo2DInterpolator; // forward declaration, cyclic header include
  * - <i>"Quantitative evaluation of different hydrological modelling approaches in a partly glacierized Swiss watershed"</i>, Jan Magnusson, Daniel Farinotti, Tobias Jonas and Mathias Bavay, Hydrological Processes, 2010, under review.
  * - <i>"Modelling runoff from highly glacierized alpine catchments in a changing climate"</i>, Matthias Huss, Daniel Farinotti, Andreas Bauder and Martin Funk, Hydrological Processes, <b>22</b>, 3888-3902, 2008.
  * - <i>"Geostatistics for Natural Resources Evaluation"</i>, Pierre Goovaerts, Oxford University Press, Applied Geostatistics Series, 1997, 483 p., ISBN 0-19-511538-4
+ * - <i>"Statistics for spatial data"</i>, Noel A. C. Cressie, John Wiley & Sons, revised edition, 1993, 900 p.
  * 
  * @author Mathias Bavay
  * @date   2010-04-12
@@ -354,6 +356,27 @@ class SnowHNWInterpolation : public InterpolationAlgorithm {
 		virtual void calculate(const MeteoData::Parameters& param, Grid2DObject& grid);
 };
 
+/**
+ * @class OrdinaryKrigingAlgorithm
+ * @brief Ordinary kriging THIS IS NOT USABLE YET
+ * @author Mathias Bavay
+ */
+class OrdinaryKrigingAlgorithm : public InterpolationAlgorithm {
+	public:
+		OrdinaryKrigingAlgorithm(const Meteo2DInterpolator& _mi,
+		               const DEMObject& _dem,
+		               const std::vector<MeteoData>& _vecMeteo,
+		               const std::vector<StationData>& _vecStation,
+		               const std::vector<std::string>& _vecArgs,
+		               const std::string _algo)
+			: InterpolationAlgorithm(_mi, _dem, _vecMeteo, _vecStation, _vecArgs, _algo) {}
+		virtual double getQualityRating(const MeteoData::Parameters& param);
+		virtual void calculate(const MeteoData::Parameters& param, Grid2DObject& grid);
+
+	private:
+		double computeVariogram(const std::vector<StationData>& vecStations) const;
+};
+
 } //end namespace mio
 
 #endif

meteoio/MainPage.h

@@ -48,6 +48,7 @@ namespace mio {
  * -# Expanding MeteoIO
  *    -# How to \subpage dev_plugins "write a Plugin"
  *    -# How to \subpage dev_filters "write a Filter"
+ *    -# How to \subpage dev_2Dinterpol "write a spatial interpolation algorithm"
  */
 
  /**
@@ -196,6 +197,7 @@ namespace mio {
 //Plugins overview given in IOHandler.cc
 //Filters overview given in FilterAlgorithms.cc
 //Plugin development given in IOInterface.h
+//2D interpolation development given in Meteo2DInterpolator.h
 
  /**
  * @page dev_filters Filter developer's guide

meteoio/Matrix.cc

@@ -165,7 +165,7 @@ Matrix& Matrix::operator+=(const Matrix& rhs) {
 	return *this;
 }
 
-const Matrix Matrix::operator+(const Matrix& rhs) {
+const Matrix Matrix::operator+(const Matrix& rhs) const {
 	Matrix result = *this;
 	result += rhs; //already implemented
 
@@ -183,7 +183,7 @@ Matrix& Matrix::operator+=(const double& rhs) {
 	return *this;
 }
 
-const Matrix Matrix::operator+(const double& rhs) {
+const Matrix Matrix::operator+(const double& rhs) const {
 	Matrix result = *this;
 	result += rhs; //already implemented
 
@@ -210,7 +210,7 @@ Matrix& Matrix::operator-=(const Matrix& rhs) {
 	return *this;
 }
 
-const Matrix Matrix::operator-(const Matrix& rhs) {
+const Matrix Matrix::operator-(const Matrix& rhs) const {
 	Matrix result = *this;
 	result -= rhs; //already implemented
 
@@ -223,7 +223,7 @@ Matrix& Matrix::operator-=(const double& rhs) {
 	return *this;
 }
 
-const Matrix Matrix::operator-(const double& rhs) {
+const Matrix Matrix::operator-(const double& rhs) const {
 	Matrix result = *this;
 	result += -rhs; //already implemented
 
@@ -258,7 +258,7 @@ Matrix& Matrix::operator*=(const Matrix& rhs) {
 	return *this;
 }
 
-const Matrix Matrix::operator*(const Matrix& rhs) {
+const Matrix Matrix::operator*(const Matrix& rhs) const {
 	Matrix result = *this;
 	result *= rhs; //already implemented
 
@@ -275,7 +275,7 @@ Matrix& Matrix::operator*=(const double& rhs) {
 	return *this;
 }
 
-const Matrix Matrix::operator*(const double& rhs) {
+const Matrix Matrix::operator*(const double& rhs) const {
 	Matrix result = *this;
 	result *= rhs; //already implemented
 
@@ -287,18 +287,36 @@ Matrix& Matrix::operator/=(const double& rhs) {
 	return *this;
 }
 
-const Matrix Matrix::operator/(const double& rhs) {
+const Matrix Matrix::operator/(const double& rhs) const {
 	Matrix result = *this;
 	result *= 1./rhs; //already implemented
 
 	return result;
 }
 
+double Matrix::scalar(const Matrix& m) {
+	return m.scalar();
+}
+
+const double Matrix::scalar() const {
+	if(ncols!=1 || nrows!=1) {
+		std::stringstream tmp;
+		tmp << "Trying to get scalar value of a non (1x1) matrix ";
+		tmp << "(" << nrows << "," << ncols << ") !";
+		throw IOException(tmp.str(), AT);
+	}
+	return operator()(1,1);
+}
+
+Matrix Matrix::T(const Matrix& m) {
+	return m.T();
+}
+
 /*void Matrix::T() {
 
 }*/
 
-const Matrix Matrix::T() {
+const Matrix Matrix::T() const {
 //other possibility: create a "transpose" flag that simply swaps the data reading...
 	Matrix result(ncols, nrows);
 	for(unsigned int i=1; i<=result.nrows; i++) {
@@ -366,7 +384,7 @@ void Matrix::LU(Matrix& L, Matrix& U) const {
 	
 }*/
 
-const Matrix Matrix::inv() {
+const Matrix Matrix::inv() const {
 //This uses an LU decomposition followed by backward and forward solving for the inverse
 	if(nrows!=ncols) {
 		std::stringstream tmp;

meteoio/Matrix.h

@@ -92,11 +92,19 @@ class Matrix {
 		double& operator ()(const unsigned int& x, const unsigned int& y);
 		double operator ()(const unsigned int& x, const unsigned int& y) const;
 
+		/**
+		* @brief Return a 1x1 matrix as a scalar
+		* @return scalar value
+		*/
+		const double scalar() const;
+		static double scalar(const Matrix& m);
+
 		/**
 		* @brief matrix transpose
 		* @return transposed matrix
 		*/
-		const Matrix T();
+		const Matrix T() const;
+		static Matrix T(const Matrix& m);
 		//void T();
 
 		/**
@@ -105,7 +113,7 @@ class Matrix {
 		* backward and forward solving of LU * A-1 = I
 		* @return inversed matrix
 		*/
-		const Matrix inv();
+		const Matrix inv() const;
 		//void inv();
 
 		/**
@@ -127,22 +135,22 @@ class Matrix {
 		friend std::ostream& operator<<(std::ostream& os, const Matrix& data);
 
 		Matrix& operator+=(const Matrix& rhs);
-		const Matrix operator+(const Matrix& rhs);
+		const Matrix operator+(const Matrix& rhs) const;
 		Matrix& operator+=(const double& rhs);
-		const Matrix operator+(const double& rhs);
+		const Matrix operator+(const double& rhs) const;
 
 		Matrix& operator-=(const Matrix& rhs);
-		const Matrix operator-(const Matrix& rhs);
+		const Matrix operator-(const Matrix& rhs) const;
 		Matrix& operator-=(const double& rhs);
-		const Matrix operator-(const double& rhs);
+		const Matrix operator-(const double& rhs) const;
 
 		Matrix& operator*=(const Matrix& rhs);
-		const Matrix operator*(const Matrix& rhs);
+		const Matrix operator*(const Matrix& rhs) const;
 		Matrix& operator*=(const double& rhs);
-		const Matrix operator*(const double& rhs);
+		const Matrix operator*(const double& rhs) const;
 
 		Matrix& operator/=(const double& rhs);
-		const Matrix operator/(const double& rhs);
+		const Matrix operator/(const double& rhs) const;
 
 		bool operator==(const Matrix&) const; ///<Operator that tests for equality
 		bool operator!=(const Matrix&) const; ///<Operator that tests for inequality

meteoio/Meteo2DInterpolator.h

@@ -27,6 +27,51 @@
 
 namespace mio {
 
+/**
+ * @page dev_2Dinterpol How to write a spatial interpolation algorithm
+ * Point measurements can be spatially interpolated by MeteoIO, through the use of interpolation
+ * algorithms. The user will then be able to choose for each meteorological parameter which
+ * interpolations could be applicable and the system will choose (through a heuristic) which one
+ * to apply at each time step (depending on the conditions of the moment, like the number of measurements).
+ *
+ * @section structure_2Dinterpol Structure
+ * The selection of which interpolation algorithm to use at any given time step, for a given parameter is
+ * performed by the Meteo2DInterpolator class. This class provides the interface for the spatial
+ * interpolations. The interpolation algorithms themselves derive from the
+ * InterpolationAlgorithm class that standardizes their public API (which would not be used by anything
+ * but a Meteo2DInterpolator object). It contains a getQualityRating() method that must return a quality
+ * index between 0 (algorithm not applicable) and 1 (perfect result if using this algorithm). This is currently
+ * only based on extremely simple heuristics, using general knowledge about the applicability of the various
+ * spatial interpolation methods depending on some obvious factors (number of measurement points, etc). The
+ * Meteo2DInterpolator object will call this method from all the algorithms listed by the user (in his io.ini
+ * configuration file) and keep the one that gets the highest score for interpolating the current parameter
+ * at the current time step. The interpolation is then done calling the algorithm's calculate method.
+ *
+ * @section implementation_2Dinterpol Implementation
+ * It is therefore necessary to create in InterpolationAlgorithms.cc (and declared in the .h) a new class,
+ * nammed after the algorithm that will be implemented and inheriting InterpolationAlgorithm. Two methods need
+ * to be implemented (the constructor being inherited from InterpolationAlgorithm and automatically called
+ * by an object factory):
+ * - double getQualityRating(const MeteoData::Parameters& param)
+ * - void calculate(const MeteoData::Parameters& param, Grid2DObject& grid)
+ * 
+ * The calculate method must properly erase and reste the grid that it receives before filling it. If necessary,
+ * (as is the case for precipitation and relative humidity, for example) the grid can be checked for min/max by
+ * calling checkMinMax() at the end of Meteo2DInterpolator::interpolate.
+ *
+ * The new class and its associated end user key must be declared in AlgorithmFactory::initStaticData() and
+ * the constructor called in AlgorithmFactory::getAlgorithm. It is recommended that any generic statistical
+ * spatial processing be implemented as a static class in libinterpol2D.cc so that it could be reused by other
+ * algorithms (see for example Interpol2D::IDW and IDWCore). In any case, proper doxygen documentation
+ * must be written alongside the implementation.
+ *
+ * @section doc_2Dinterpol Documentation
+ * The newly added interpolation algorithm must be added to the list of available algorithms in
+ * InterpolationAlgorithms.h with a proper description. An example can also be given in the example section
+ * of the same file. Please feel free to add necessary bibliographic references to the bibliographic section!
+ * 
+*/
+
 /**
  * @class Meteo2DInterpolator
  * @brief A class to spatially interpolate meteo parameters. For more, see \ref interpol2d

meteoio/ResamplingAlgorithms.cc

@@ -299,7 +299,7 @@ void ResamplingAlgorithms::Accumulate(const unsigned int& pos, const MeteoData::
           vecM[pos].param(paramindex) = sum;
           return;
 	}
-          
+	
 	if ((interval_end+1) == pos){
 		valend = funcval(vecM, interval_end, vecM[pos].date, paramindex);
           if (valend != IOUtils::nodata)
@@ -324,8 +324,8 @@ void ResamplingAlgorithms::Accumulate(const unsigned int& pos, const MeteoData::
 	//TODO:check if at least one point has been summed. If not -> nodata
 }
 
-double ResamplingAlgorithms::funcval(const vector<MeteoData>& vecM, const unsigned int& index,
-							  const Date& date, const unsigned int& paramindex)
+double ResamplingAlgorithms::funcval(const std::vector<MeteoData>& vecM, const unsigned int& index,
+                                     const Date& date, const unsigned int& paramindex)
 {
 	unsigned int start = index;
 	if (vecM[start].isResampled()){
@@ -353,8 +353,8 @@ double ResamplingAlgorithms::funcval(const vector<MeteoData>& vecM, const unsign
 			return IOUtils::nodata;
 
 		return Interpol1D::linearInterpolation(vecM[start].date.getJulianDate(), 0.0,
-									    vecM[end].date.getJulianDate(), valend,
-									    date.getJulianDate());
+		                                       vecM[end].date.getJulianDate(), valend,
+		                                       date.getJulianDate());
 	}
 
 	return IOUtils::nodata;

meteoio/libinterpol2D.cc

@@ -653,4 +653,92 @@ void Interpol2D::PrecipSnow(const DEMObject& dem, const Grid2DObject& ta, Grid2D
 	}
 }
 
+/**
+* @brief Ordinary Kriging matrix formulation
+* This implements the matrix formulation of Ordinary Kriging, as shown (for example) in 
+* <i>"Statistics for spatial data"</i>, Noel A. C. Cressie, John Wiley & Sons, revised edition, 1993, pp122.
+* @param vecData vector containing the values as measured at the stations
+* @param vecStations vector of stations
+* @param dem digital elevation model
+* @param grid 2D array of precipitation to fill
+* @author Mathias Bavay
+*/
+void Interpol2D::ODKriging(const std::vector<double>& vecData, const std::vector<StationData>& vecStations, const DEMObject& dem, Grid2DObject& grid)
+{
+	grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);
+	unsigned int nrOfMeasurments = vecStations.size();
+	Matrix G(nrOfMeasurments, nrOfMeasurments);
+	Matrix gamma((unsigned int)1, nrOfMeasurments);
+	const Matrix One((unsigned int)1, nrOfMeasurments, 1.);
+	const Matrix One_T = One.T();
+
+	//fill the G matrix
+	for(unsigned int i=1; i<=nrOfMeasurments; i++) {
+		const Coords& st1 = vecStations[i].position;
+		const double x1 = st1.getEasting();
+		const double y1 = st1.getNorthing();
+
+		for(unsigned int j=1; j<=nrOfMeasurments; j++) {
+			//compute distance between stations
+			const Coords& st2 = vecStations[j].position;
+			const double DX = x1-st2.getEasting();
+			const double DY = y1-st2.getNorthing();
+			const double distance = fastSqrt_Q3(DX*DX + DY*DY);
+
+			G(i,j) = varioFit(distance);
+		}
+	}
+
+	//G inverse matrix
+	const Matrix Ginv = G.inv();
+
+	//calculate constant denominator
+	const Matrix OneT_Ginv = One_T * Ginv;
+	const double denom = Matrix::scalar( OneT_Ginv * One );
+
+	//now, calculate each point
+	for(unsigned int i=1; i<=grid.ncols; i++) {
+		for(unsigned int j=1; j<=grid.nrows; j++) {
+			const double x = grid.llcorner.getEasting()+i*grid.cellsize;
+			const double y = grid.llcorner.getNorthing()+j*grid.cellsize;
+
+			//fill gamma
+			for(unsigned int st=1; st<=nrOfMeasurments; st++) {
+				//compute distance between cell and each station
+				const Coords& position = vecStations[st].position;
+				const double DX = x-position.getEasting();
+				const double DY = y-position.getNorthing();
+				const double distance = fastSqrt_Q3(DX*DX + DY*DY);
+
+				gamma(1,st) = varioFit(distance);
+			}
+
+			const Matrix lambda = Matrix::T(gamma + One * ((1. - Matrix::scalar(OneT_Ginv*gamma)) / denom) ) * Ginv;
+
+			//calculate local parameter interpolation
+			double p = 0.;
+			for(unsigned int st=1; st<=nrOfMeasurments; st++) {
+				p += lambda(1,st) * vecData[st-1]; //because the vector starts at 0
+			}
+			grid.grid2D(i,j) = p;
+		}
+	}
+}
+
+/**
+* @brief Return the fitted variogram value for a given distance
+* The various variogram models can be found in
+* <i>"Statistics for spatial data"</i>, Noel A. C. Cressie, John Wiley & Sons, revised edition, 1993, pp63.
+* @param distance distance to the measured value
+* @return variogram fit
+* @author Mathias Bavay
+*/
+double Interpol2D::varioFit(const double& /*distance*/)
+{
+	//return variogramm fit of covariance between stations i and j
+	//HACK: todo!
+	
+	return 1.;
+}
+
 } //namespace

meteoio/libinterpol2D.h

@@ -22,14 +22,7 @@
 #ifndef INTERPOL2D_H
 #define INTERPOL2D_H
 
-#include <meteoio/StationData.h>
-#include <meteoio/MeteoData.h>
-#include <meteoio/Grid2DObject.h>
-#include <meteoio/DEMObject.h>
-#include <meteoio/Date.h>
-#include <meteoio/IOExceptions.h>
-#include <meteoio/IOUtils.h>
-#include <meteoio/libinterpol1D.h>
+#include <meteoio/MeteoIO.h>
 
 #include <cmath>
 #include <vector>
@@ -37,8 +30,8 @@
 #include <assert.h>
 
 #define MAX_INPUT_STATIONS 255
-#define GRAVITY	9.80665		     // (m s-2)
-#define GAS_CONSTANT_AIR 287.	     // (J kg-1 K-1)
+#define GRAVITY	9.80665		     // (m s-2) //HACK: move to libmeteolaws
+#define GAS_CONSTANT_AIR 287.	     // (J kg-1 K-1) //HACK: move to libmeteolaws
 
 namespace mio {
 
@@ -76,6 +69,9 @@ class Interpol2D {
                                 const DEMObject& dem, Grid2DObject& grid);
 		static void SimpleDEMWindInterpolate(const DEMObject& dem, Grid2DObject& VW, Grid2DObject& DW);
 		static void PrecipSnow(const DEMObject& dem, const Grid2DObject& ta, Grid2DObject& grid);
+		static void ODKriging(const std::vector<double>& vecData,
+		                      const std::vector<StationData>& vecStations,
+		                      const DEMObject& dem, Grid2DObject& grid);
 
 		//projections functions
 		static double ConstProject(const double& value, const double& altitude, const double& new_altitude,
@@ -113,6 +109,9 @@ class Interpol2D {
 		double weightInvDistN(const double& d2);
 		double dist_pow; //power for the weighting method weightInvDistN
 
+		//variogram methods
+		static double varioFit(const double& distance);
+
 	private:
 		//static members
 		const static double wind_ys; ///coefficient for wind dependency on slope

meteoio/plugins/CMakeLists.txt

@@ -124,13 +124,10 @@ IF(PLUGIN_GSNIO)
 ENDIF(PLUGIN_GSNIO)
 
 IF(PLUGIN_IMISIO)
-  set(ORACLE $ENV{ORACLE_HOME})
-  find_path(ORACLE_INCLUDE_DIR occi.h "${ORACLE}/client_1/include")
-  find_library(ORACLE_OCCI NAMES occi PATHS ${ORACLE}/client_1/lib)
-  find_library(ORACLE_CLNTSH NAMES clntsh PATHS ${ORACLE}/client_1/lib)
-  INCLUDE_DIRECTORIES(${INCLUDE_DIRECTORIES} ${ORACLE_INCLUDE_DIR})
+  FIND_PACKAGE(OCCI)
+  INCLUDE_DIRECTORIES(${INCLUDE_DIRECTORIES} ${OCCI_INCLUDE_DIR})
   ADD_LIBRARY(${IMISIONAME} ImisIO.cc)
-  TARGET_LINK_LIBRARIES(${IMISIONAME} ${ORACLE_OCCI} ${ORACLE_CLNTSH} ${SHAREDNAME})
+  TARGET_LINK_LIBRARIES(${IMISIONAME} ${OCCI_SHARED} ${SHAREDNAME})
   SET_TARGET_PROPERTIES(${IMISIONAME} PROPERTIES 
     LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/lib/plugins
     VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH} 

meteoio/plugins/ImisIO.cc

@@ -436,8 +436,8 @@ void ImisIO::assimilateAnetzData(const unsigned int& indexStart, const unsigned
 	}
 }
 
-double ImisIO::getHNW(const vector<MeteoData>& vecAnetz, const AnetzData& ad, const unsigned int& index, 
-				  const std::map<std::string, unsigned int>& mapAnetzNames)
+double ImisIO::getHNW(const std::vector<MeteoData>& vecAnetz, const AnetzData& ad, const unsigned int& index, 
+                      const std::map<std::string, unsigned int>& mapAnetzNames)
 {
 	double hnw = 0.0;
 	map<string, unsigned int>::const_iterator it;

tools/cmake/FindOCCI.cmake

+# - Find Oracle OCCI library
+# Find the Oracle OCCI library and includes
+#
+# It defines the following variables:
+# OCCI_LIBRARY - the library
+# OCCI_INCLUDE_DIR - the header
+# OCCI_FOUND - Set to false, or undefined, if we haven't found
+
+include(LibFindMacros)
+
+# Finally the library itself
+find_library(OCCI_LIBRARY
+  NAMES occi clntsh
+  HINTS $ENV{ORACLE_HOME}/client_1/lib $ENV{ORACLE_BASE}/client_1/lib
+  PATHS "/usr/lib" "/usr/local/lib" "~/usr/lib" "/opt/lib" ENV LD_LIBRARY_PATH
+  PATH_SUFFIXES "oracle/*/client_1/lib"