WSL/SLF GitLab Repository

InterpolationAlgorithms.h 24 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
/***********************************************************************************/
/*  Copyright 2010 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    MeteoIO is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef __INTERPOLATIONALGORITHMS_H__
#define __INTERPOLATIONALGORITHMS_H__

21
22
#include <meteoio/DEMObject.h>
#include <meteoio/MeteoData.h>
23
#include <meteoio/meteostats/libinterpol1D.h>
24
#include <meteoio/meteostats/libinterpol2D.h>
25
#include <meteoio/meteostats/libfit1D.h>
26
27
28
29
30
31
32

#include <vector>
#include <string>
#include <set>

namespace mio {

33
class IOManager;
34
35
36
37
38
39
class Meteo2DInterpolator; // forward declaration, cyclic header include

/**
 * @page interpol2d Spatial interpolations
 * Using the vectors of MeteoData and StationData as filled by the IOInterface::readMeteoData call
 * as well as a grid of elevations (DEM, stored as a DEMObject), it is possible to get spatially
40
 * interpolated parameters.
41
42
43
44
45
46
47
48
 *
 * First, an interpolation method has to be selected for each variable which needs interpolation. Then the class computes
 * the interpolation for each 2D grid point, combining the inputs provided by the available data sources.
 * Any parameter of MeteoData can be interpolated, using the names given by \ref meteoparam. One has to keep
 * in mind that the interpolations are time-independent: each interpolation is done at a given time step and no
 * memory of (eventual) previous time steps is kept. This means that all parameters and variables that are
 * automatically calculated get recalculated anew for each time step.
 *
49
 * @section interpol2D_section Spatial interpolations section
50
51
 * Practically, the user
 * has to specify in his configuration file (typically io.ini), for each parameter to be interpolated, which
52
 * spatial interpolations algorithms should be considered, in the [Interpolations2D] section. This is provided as a space separated list of keywords
53
54
55
 * (one per interpolation algorithm). Please notice that some algorithms may require extra arguments.
 * Then, each algorithm will be evaluated (through the use of its rating method) and receive a grade (that might
 * depend on the number of available data, the quality of the data, etc). The algorithm that receives the higher
56
 * score within the user list, will be used for interpolating the selected variable. An example of such section is given below:
57
58
59
60
 * @code
 * [Interpolations2D]
 * TA::algorithms = IDW_LAPSE CST_LAPSE
 * TA::cst_lapse = -0.008
61
 *
62
 * RH::algorithms = RH IDW_LAPSE CST_LAPSE CST
63
 *
64
65
66
 * HNW::algorithms = HNW_SNOW IDW_LAPSE CST_LAPSE CST
 * HNW::hnw_snow = cst_lapse
 * HNW::cst_lapse = 0.0005 frac
67
 *
68
 * VW::algorithms = IDW_LAPSE CST_LAPSE
69
 *
70
71
72
 * P::algorithms = STD_PRESS
 * @endcode
 *
73
 * @section interpol2D_keywords Available algorithms
74
75
76
77
78
79
80
 * The keywords defining the algorithms are the following:
 * - STD_PRESS: standard atmospheric pressure as a function of the elevation of each cell (see StandardPressureAlgorithm)
 * - CST: constant value in each cell (see ConstAlgorithm)
 * - CST_LAPSE: constant value reprojected to the elevation of the cell (see ConstLapseRateAlgorithm)
 * - IDW: Inverse Distance Weighting averaging (see IDWAlgorithm)
 * - IDW_LAPSE: Inverse Distance Weighting averaging with reprojection to the elevation of the cell (see IDWLapseAlgorithm)
 * - RH: the dew point temperatures are interpolated using IDW_LAPSE, then reconverted locally to relative humidity (see RHAlgorithm)
81
 * - ILWR: the incoming long wave radiation is converted to emissivity and then interpolated (see ILWRAlgorithm)
82
 * - 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)
83
 * - HNW_SNOW: precipitation interpolation according to (Magnusson, 2011) (see SnowHNWInterpolation)
84
 * - ODKRIG: ordinary kriging (see OrdinaryKrigingAlgorithm)
85
 * - ODKRIG_LAPSE: ordinary kriging with lapse rate (see LapseOrdinaryKrigingAlgorithm)
86
 * - USER: user provided grids to be read from disk (if available, see USERInterpolation)
87
 * <!-- - LIDW_LAPSE: IDW_LAPSE restricted to a local scale (n neighbor stations, see LocalIDWLapseAlgorithm) -->
88
 *
89
 * @section interpol2D_lapse Lapse rates
90
91
92
93
94
95
96
 * Several algorithms use elevation trends, currently modelled as a linear relation. The slope of this linear relation can
 * sometimes be provided by the end user (through his io.ini configuration file), otherwise it is computed from the data.
 * In order to bring slightly more robustness, if the correlation between the input data and the computed linear regression
 * is not good enought (below 0.7, as defined in Interpol2D::LinRegression), the same regression will get re-calculated
 * with one point less (cycling throught all the points). The best result (ie: highest correlation coefficient) will be
 * kept. If the final correlation coefficient is less than 0.7, a warning is displayed.
 *
97
 * @section interpol2D_dev_use Developer usage
98
 * From the developer's point of view, all that has to be done is instantiate an IOManager object and call its
99
 * IOManager::interpolate method.
100
 * @code
101
 * 	Config cfg("io.ini");
102
 * 	IOManager io(cfg);
103
 *
104
 * 	//reading the dem (necessary for several spatial interpolations algoritms)
105
 * 	DEMObject dem;
106
 * 	io.readDEM(dem);
107
 *
108
 *	//performing spatial interpolations
109
 * 	Grid2DObject param;
110
 *	io.interpolate(date, dem, MeteoData::TA, param);
111
 *
112
113
 * @endcode
 *
114
 * @section interpol2D_biblio Bibliography
115
 * The interpolation algorithms have been inspired by the following papers:
116
117
118
119
120
 * - <i>"A Meteorological Distribution System for High-Resolution Terrestrial Modeling (MicroMet)"</i>, Liston and Elder, Journal of Hydrometeorology <b>7</b> (2006), 217-234.
 * - <i>"Simulating wind fields and snow redistribution using terrain-based parameters to model snow accumulation and melt over a semi-arid mountain catchment"</i>, Adam Winstral and Danny Marks, Hydrological Processes <b>16</b> (2002), 3585– 3603. DOI: 10.1002/hyp.1238 [NOT YET IMPLEMENTED]
 * - <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
121
 * - <i>"Statistics for spatial data"</i>, Noel A. C. Cressie, John Wiley & Sons, revised edition, 1993, 900 p.
122
 *
123
124
125
126
127
128
129
 * @author Mathias Bavay
 * @date   2010-04-12
 */

/**
 * @class InterpolationAlgorithm
 * @brief A class to perform 2D spatial interpolations. For more, see \ref interpol2d
130
131
 *
 * @ingroup stats
132
133
134
135
136
137
 * @author Thomas Egger
 * @date   2010-04-01
*/
class InterpolationAlgorithm {

	public:
138
139
		InterpolationAlgorithm(Meteo2DInterpolator& i_mi,
		                       const std::vector<std::string>& i_vecArgs,
140
		                       const std::string& i_algo, IOManager& iom) :
141
142
		                      algo(i_algo), mi(i_mi), date(0.), vecArgs(i_vecArgs), vecMeteo(), vecData(),
		                      vecMeta(), info(), param(MeteoData::firstparam), nrOfMeasurments(0), iomanager(iom) {};
143
144
		virtual ~InterpolationAlgorithm() {};
		//if anything is not ok (wrong parameter for this algo, insufficient data, etc) -> return zero
145
146
		virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param) = 0;
		virtual void calculate(const DEMObject& dem, Grid2DObject& grid) = 0;
147
		std::string getInfo() const;
148
149
		const std::string algo;

150
 	protected:
151
152
153
		size_t getData(const Date& i_date, const MeteoData::Parameters& i_param, std::vector<double>& o_vecData);
		size_t getData(const Date& i_date, const MeteoData::Parameters& i_param,
		               std::vector<double>& o_vecData, std::vector<StationData>& o_vecMeta);
154
		size_t getStationAltitudes(const std::vector<StationData>& i_vecMeta, std::vector<double>& o_vecData) const;
155
156
		void getTrend(const std::vector<double>& vecAltitudes, const std::vector<double>& vecDat, Fit1D &trend) const;
		void detrend(const Fit1D& trend, const std::vector<double>& vecAltitudes, std::vector<double> &vecDat) const;
157
		void retrend(const DEMObject& dem, const Fit1D& trend, Grid2DObject &grid) const;
158

159
		Meteo2DInterpolator& mi;
160
		Date date;
161
		const std::vector<std::string> vecArgs; //we must keep our own copy, it is different for each algorithm!
162
163

		std::vector<MeteoData> vecMeteo;
164
165
		std::vector<double> vecData; ///<store the measurement for the given parameter
		std::vector<StationData> vecMeta; ///<store the station data for the given parameter
166
		std::ostringstream info; ///<to store some extra information about the interplation process
167
		MeteoData::Parameters param; ///<the parameter that we will interpolate
168
		size_t nrOfMeasurments; ///<the available number of measurements
169
		IOManager& iomanager;
170
171
172
173
};

class AlgorithmFactory {
	public:
174
		static InterpolationAlgorithm* getAlgorithm(const std::string& i_algoname,
175
		                                            Meteo2DInterpolator& i_mi,
176
		                                            const std::vector<std::string>& i_vecArgs, IOManager& iom);
177
178
179
180
};

/**
 * @class ConstAlgorithm
181
 * @brief Constant filling interpolation algorithm.
182
183
184
185
 * Fill the grid with the average of the inputs for this parameter.
 */
class ConstAlgorithm : public InterpolationAlgorithm {
	public:
186
187
188
		ConstAlgorithm(Meteo2DInterpolator& i_mi,
					const std::vector<std::string>& i_vecArgs,
					const std::string& i_algo, IOManager& iom)
189
			: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom) {}
190
191
		virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
		virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
192
193
194
195
196
197
198
199
200
};

/**
 * @class StandardPressureAlgorithm
 * @brief Standard atmospheric pressure interpolation algorithm.
 * Fill the grid with the standard atmosphere's pressure, depending on the local elevation.
 */
class StandardPressureAlgorithm : public InterpolationAlgorithm {
	public:
201
202
203
		StandardPressureAlgorithm(Meteo2DInterpolator& i_mi,
					const std::vector<std::string>& i_vecArgs,
					const std::string& i_algo, IOManager& iom)
204
			: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom) {}
205
206
		virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
		virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
207
208
209
210
211
212
};

/**
 * @class ConstLapseRateAlgorithm
 * @brief Constant filling with elevation lapse rate interpolation algorithm.
 * Assuming that average values occured at the average of the elevations, the grid is filled with average values
213
214
215
 * reprojected to real grid elevation according to a lapse rate. The lapse rate is either calculated from the data
 * (if no extra argument is provided), or given by the user-provided the optional argument <i>"cst_lapse"</i>.
 * If followed by <i>"soft"</i>, then an attempt to calculate the lapse rate from the data is made, any only if
216
217
 * unsuccessful, then user provided lapse rate is used as a fallback. If the optional user given lapse rate is
 * followed by <i>"frac"</i>, then the lapse rate is understood as a fractional lapse rate, that is a relative change
218
 * of the value as a function of the elevation (for example, +0.05% per meters given as 0.0005). In this case, no attempt to calculate
219
 * the fractional lapse from the data is made.
220
221
222
 */
class ConstLapseRateAlgorithm : public InterpolationAlgorithm {
	public:
223
224
225
		ConstLapseRateAlgorithm(Meteo2DInterpolator& i_mi,
					const std::vector<std::string>& i_vecArgs,
					const std::string& i_algo, IOManager& iom)
226
			: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom) {}
227
228
		virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
		virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
229
230
231
232
233
234
235
236
237
238
239
};

/**
 * @class IDWAlgorithm
 * @brief Inverse Distance Weighting interpolation algorithm.
 * Each cell receives the weighted average of the whole data set with weights being 1/r²
 * (r being the distance of the current cell to the contributing station) and renormalized
 * (so that the sum of the weights is equal to 1.0).
 */
class IDWAlgorithm : public InterpolationAlgorithm {
	public:
240
241
242
		IDWAlgorithm(Meteo2DInterpolator& i_mi,
					const std::vector<std::string>& i_vecArgs,
					const std::string& i_algo, IOManager& iom)
243
			: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom) {}
244
245
		virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
		virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
246
247
248
249
250
251
252
};

/**
 * @class IDWLapseAlgorithm
 * @brief Inverse Distance Weighting interpolation algorithm with elevation detrending/reprojection.
 * The input data is projected to a reference elevation and spatially interpolated using an Inverse Distance
 * Weighting interpolation algorithm (see IDWAlgorithm). Then, each value is reprojected to the real
253
254
255
 * elevation of the relative cell. The lapse rate is either calculated from the data
 * (if no extra argument is provided), or given by the user-provided the optional argument <i>"idw_lapse"</i>.
 * If followed by <i>"soft"</i>, then an attempt to calculate the lapse rate from the data is made, any only if
256
257
258
 * unsuccessful or too bad (r^2<0.6), then the user provided lapse rate is used as a fallback.
 * If the optional user given lapse rate is
 * followed by <i>"frac"</i>, then the lapse rate is understood as a fractional lapse rate, that is a relative change
259
 * of the value as a function of the elevation (for example, +0.05% per meters given as 0.0005). In this case, no attempt to calculate
260
 * the fractional lapse from the data is made.
261
262
263
 */
class IDWLapseAlgorithm : public InterpolationAlgorithm {
	public:
264
265
266
		IDWLapseAlgorithm(Meteo2DInterpolator& i_mi,
					const std::vector<std::string>& i_vecArgs,
					const std::string& i_algo, IOManager& iom)
267
			: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom) {}
268
269
		virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
		virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
270
271
272
273
274
275
276
277
};


/**
 * @class LocalIDWLapseAlgorithm
 * @brief Inverse Distance Weighting interpolation algorithm with elevation detrending/reprojection.
 * The closest n stations (n being given as an extra argument of <i>"idw_lapse"</i>) to each pixel are
 * used to compute the local lapse rate, allowing to project the contributions of these n stations to the
278
 * local pixel with an inverse distance weight. We also assume a two segments regression for altitude detrending with
279
 * a fixed 1200m above sea level inflection point.
280
281
282
 */
class LocalIDWLapseAlgorithm : public InterpolationAlgorithm {
	public:
283
284
285
		LocalIDWLapseAlgorithm(Meteo2DInterpolator& i_mi,
					const std::vector<std::string>& i_vecArgs,
					const std::string& i_algo, IOManager& iom)
286
			: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom) {}
287
288
		virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
		virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
289
290
291
292
293
294
295
296
297
298
299
300
301
302
};

/**
 * @class RHAlgorithm
 * @brief Relative humidity interpolation algorithm.
 * This is an implementation of the method described in (Liston & Elder, 2006): for each input point, the dew
 * point temperature is calculated. Then, the dew point temperatures are spatially interpolated using IDWLapseAlgorithm.
 * Finally, each local dew point temperature is converted back to a local relative humidity.
 *
 * As a side effect, the user must have defined algorithms to be used for air temperature (since this is needed for dew
 * point to RH conversion)
 */
class RHAlgorithm : public InterpolationAlgorithm {
	public:
303
304
305
		RHAlgorithm(Meteo2DInterpolator& i_mi,
					const std::vector<std::string>& i_vecArgs,
					const std::string& i_algo, IOManager& iom)
306
			: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom), vecDataTA(), vecDataRH() {}
307
308
		virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
		virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
309
310
	private:
		std::vector<double> vecDataTA, vecDataRH; ///<vectors of extracted TA and RH
311
312
};

313
314
315
316
317
318
/**
 * @class ILWRAlgorithm
 * @brief Incoming Long Wave Radiation interpolation algorithm.
 * Each ILWR is converted to an emissivity (using the local air temperature), interpolated using CST_LAPSE or IDW_LAPSE with
 * a fixed lapse rate and reconverted to ILWR.
 *
319
 * As a side effect, the user must have defined algorithms to be used for air temperature (since this is needed for
320
321
322
323
324
325
326
 * emissivity to ILWR conversion)
 */
class ILWRAlgorithm : public InterpolationAlgorithm {
	public:
		ILWRAlgorithm(Meteo2DInterpolator& i_mi,
					const std::vector<std::string>& i_vecArgs,
					const std::string& i_algo, IOManager& iom)
327
			: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom), vecDataEA() {}
328
329
		virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
		virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
330
331
332
333
	private:
		std::vector<double> vecDataEA; ///<vectors of extracted emissivities
};

334
335
336
337
338
339
340
341
342
343
344
/**
 * @class SimpleWindInterpolationAlgorithm
 * @brief Curvature/slope influenced  wind interpolation algorithm.
 * This is an implementation of the method described in (Liston & Elder, 2006): the wind speed and direction are
 * spatially interpolated using IDWLapseAlgorithm for the wind speed and using the user defined wind direction
 * interpolation algorithm. Then, the wind speed and direction fields are altered by wind weighting factors
 * and wind diverting factors (respectively) calculated from the local curvature and slope
 * (as taken from the DEM, see DEMObject).
 */
class SimpleWindInterpolationAlgorithm : public InterpolationAlgorithm {
	public:
345
346
347
		SimpleWindInterpolationAlgorithm(Meteo2DInterpolator& i_mi,
					const std::vector<std::string>& i_vecArgs,
					const std::string& i_algo, IOManager& iom)
348
			: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom), vecDataVW(), vecDataDW() {}
349
350
		virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
		virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
351
352
	private:
		std::vector<double> vecDataVW, vecDataDW; ///<vectors of extracted VW and DW
353
354
355
};

/**
356
 * @class USERInterpolation
357
 * @brief Reads user provided gridded data on the disk.
358
359
360
361
362
363
364
365
366
367
 * The grids are all in a directory that is given as the algorithm's argument. The files must be named
 * according to the following schema:
 * - {numeric date}_{capitalized meteo parameter}.asc, for example 200812011500_TA.asc
 * - Default_{capitalized meteo parameter}.asc for the grid to use when no measurements exist (which prevents
 * retrieving the date for the interpolation)
 * The meteo parameters can be found in \ref meteoparam "MeteoData". Example of use:
 * @code
 * TA::algorithms = USER
 * TA::user = ./meteo_grids
 * @endcode
368
 *
369
 */
370
class USERInterpolation : public InterpolationAlgorithm {
371
	public:
372
373
374
		USERInterpolation(Meteo2DInterpolator& i_mi,
					const std::vector<std::string>& i_vecArgs,
					const std::string& i_algo, IOManager& iom)
375
			: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom), filename() {nrOfMeasurments=0;}
376
377
		virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
		virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
378
	private:
379
		std::string getGridFileName() const;
380
		std::string filename;
381
382
};

383
/**
384
 * @class SnowHNWInterpolation
385
 * @brief Precipitation distribution according to the local slope and curvature.
386
387
388
 * The precipitation distribution is initialized using a specified algorithm (IDW_LAPSE by default, see IDWLapseAlgorithm).
 * An optional parameter can be given to specify which algorithm has to be used for initializing the grid.
 * Please do not forget to provide the arguments of the chosen algorithm itself if necessary!
389
 *
390
391
 * After this initialization, the pixels whose air temperatures are below or at freezing are modified according
 * to the method described in <i>"Quantitative evaluation of different hydrological modelling approaches
392
 * in a partly glacierized Swiss watershed"</i>, Magnusson et Al., Hydrological Processes, <b>25</b>, 2071-2084, 2011 and
393
 * <i>"Modelling runoff from highly glacierized alpine catchments in a changing climate"</i>, Huss et All., Hydrological Processes, <b>22</b>, 3888-3902, 2008.
394
395
396
397
398
399
400
401
 *
 * An example using this algorithm, initializing the grid with a constant lapse rate fill using +0.05% precipitation increase per meter of elevation, is given below:
 * @code
 * HNW::algorithms = HNW_SNOW
 * HNW::hnw_snow = cst_lapse
 * HNW::cst_lapse = 0.0005 frac
 * @endcode
 *
402
 * @author Florian Kobierska, Jan Magnusson and Mathias Bavay
403
 */
404
class SnowHNWInterpolation : public InterpolationAlgorithm {
405
	public:
406
407
408
		SnowHNWInterpolation(Meteo2DInterpolator& i_mi,
					const std::vector<std::string>& i_vecArgs,
					const std::string& i_algo, IOManager& iom)
409
  			: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom), internal_dem() {}
410
		virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
411
412
413
		virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
    private:
        DEMObject internal_dem;
414
};
415

416
417
/**
 * @class OrdinaryKrigingAlgorithm
418
 * @brief Ordinary kriging.
419
 * This implements ordinary krigging (see https://secure.wikimedia.org/wikipedia/en/wiki/Kriging)
420
421
 * with user-selectable variogram model (see https://secure.wikimedia.org/wikipedia/en/wiki/Variogram).
 * The variogram and krigging
422
 * coefficients are re-computed fresh for each new grid (or time step).
423
424
425
426
427
 *
 * The variogram is currently computed with the current data (as 1/2*(X1-X2)^2), which makes it quite
 * uninteresting... The next improvement will consist in calculating the covariances (used to build the
 * variogram) from time series (thus reflecting the time-correlation between stations).
 *
428
429
430
 * The available variogram models are found in Fit1D::regression and given as optional arguments
 * (by default, LINVARIO is used). Several models can be given, the first that can fit the data will be used
 * for the current timestep:
431
432
 * @code
 * TA::algorithms = ODKRIG
433
 * TA::odkrig = SPHERICVARIO linvario
434
435
 * @endcode
 *
436
437
438
439
 * @author Mathias Bavay
 */
class OrdinaryKrigingAlgorithm : public InterpolationAlgorithm {
	public:
440
441
442
		OrdinaryKrigingAlgorithm(Meteo2DInterpolator& i_mi,
					const std::vector<std::string>& i_vecArgs,
					const std::string& i_algo, IOManager& iom)
443
			: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom), variogram() {}
444
445
		virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
		virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
446
447
	protected:
		void getDataForVariogram(std::vector<double> &distData, std::vector<double> &variData);
448
		bool computeVariogram();
449
		Fit1D variogram;
450
451
};

452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468

/**
 * @class LapseOrdinaryKrigingAlgorithm
 * @brief Ordinary kriging with detrending.
 * This is very similar to OrdinaryKrigingAlgorithm but performs detrending on the data.
 * @code
 * TA::algorithms = ODKRIG_LAPSE
 * TA::odkrig = SPHERICVARIO
 * @endcode
 *
 * @author Mathias Bavay
 */
class LapseOrdinaryKrigingAlgorithm : public OrdinaryKrigingAlgorithm {
	public:
		LapseOrdinaryKrigingAlgorithm(Meteo2DInterpolator& i_mi,
					const std::vector<std::string>& i_vecArgs,
					const std::string& i_algo, IOManager& iom)
469
			: OrdinaryKrigingAlgorithm(i_mi, i_vecArgs, i_algo, iom) {}
470
		virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
471
472
};

473
474
475
} //end namespace mio

#endif