WSL/SLF GitLab Repository

FilterAlgorithms.cc 39.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
/***********************************************************************************/
/*  Copyright 2009 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/>.
*/
#include <meteoio/FilterAlgorithms.h>

using namespace std;

namespace mio {
 /**
 * @page filters Filters overview
 * The filtering infrastructure is described in FilterAlgorithms (for its API). The goal of this page is to give an overview of the available filters and their usage.
 *
 * @section filters_modes Modes of operation
 * It should be noted that filters often have two modes of operations: soft or hard. In soft mode, all value that is rejected is replaced by the filter parameter's value. This means that for a soft min filter set at 0.0, all values less than 0.0 will be replaced by 0.0. In hard mode, all rejected values are replaced by nodata.
 *
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
 * @section filtering_section Filtering section
 * The filters are specified for each parameter in the [Filters] section. This section contains
 * a list of the various meteo parameters with their associated choice of filtering algorithms and
 * optional parameters.The filters are applied serialy, in the order they are given in. An example of such section is given below:
 * @code
 * [Filters]
 * TA::filter1	= min_max
 * TA::arg1	= 230 330
 * 
 * RH::filter1	= min_max
 * RH::arg1	= -0.2 1.2
 * RH::filter2	= min_max
 * RH::arg2	= soft 0.0 1.0
 * 
 * HNW::filter1	= min
 * HNW::arg1	= -0.1
 * HNW::filter2	= min
 * HNW::arg2	= soft 0.
 * @endcode
 *
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
 * @section filters_available Available filters
 * The filters that are currently available are the following:
 * - rate: rate of change filter, see FilterAlgorithms::RateFilter
 * - min_max: range check filter, see FilterAlgorithms::MinMaxFilter
 * - min: minimum check filter, see FilterAlgorithms::MinValueFilter
 * - max: maximum check filter, see FilterAlgorithms::MaxValueFilter
 * - mad: median absolute deviation, see FilterAlgorithms::MedianAbsoluteDeviationFilter
 *
 * A few data transformations are also supported besides filtering:
 * - accumulate: data accumulates over a given period, see FilterAlgorithms::AccumulateProcess
 * - median_avg: running median average over a given window, see FilterAlgorithms::MedianAvgProcess
 * - mean_avg: running mean average over a given window, see FilterAlgorithms::MeanAvgProcess
 * - wind_avg: vector average over a given window, see FilterAlgorithms::WindAvgProcess
 */

std::map<std::string, FilterProperties> FilterAlgorithms::filterMap;
const bool FilterAlgorithms::__init = FilterAlgorithms::initStaticData();

bool FilterAlgorithms::initStaticData()
{
70
71
72
73
74
75
76
77
78
79
	filterMap["rate"]          = FilterProperties(true,  &FilterAlgorithms::RateFilter);
	filterMap["min_max"]       = FilterProperties(true,  &FilterAlgorithms::MinMaxFilter);
	filterMap["min"]           = FilterProperties(true,  &FilterAlgorithms::MinValueFilter);
	filterMap["max"]           = FilterProperties(true,  &FilterAlgorithms::MaxValueFilter);
	filterMap["mad"]           = FilterProperties(true,  &FilterAlgorithms::MedianAbsoluteDeviationFilter);
	filterMap["accumulate"]    = FilterProperties(false, &FilterAlgorithms::AccumulateProcess);
	filterMap["exp_smoothing"] = FilterProperties(false, &FilterAlgorithms::ExpSmoothingFilter);
	filterMap["median_avg"]    = FilterProperties(false, &FilterAlgorithms::MedianAvgProcess);
	filterMap["mean_avg"]      = FilterProperties(false, &FilterAlgorithms::MeanAvgProcess);
	filterMap["wind_avg"]      = FilterProperties(false, &FilterAlgorithms::WindAvgProcess);
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270

	return true;
}

const FilterProperties& FilterAlgorithms::filterProperties(const std::string& filtername)
{
	std::map<std::string, FilterProperties>::const_iterator it;
	it = filterMap.find(filtername);

	if (it==filterMap.end())
		throw UnknownValueException("Unknown filter called: " + filtername, AT);

	return it->second;
}

void FilterAlgorithms::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)
{
	isSoft = false;
	unsigned int argindex = 0;
	vecArgs_out.clear();
	double tmp;

	//Test for softness
	if (vecArgs_in.size() > 0)
		if (vecArgs_in[0] == "soft"){ isSoft = true; argindex=1;}

	try {
		for (unsigned int ii=argindex; ii<vecArgs_in.size(); ii++){
			IOUtils::convertString(tmp, vecArgs_in[ii]);
			vecArgs_out.push_back(tmp);
		}

		if ((vecArgs_out.size() < minArgs) || (vecArgs_out.size() > maxArgs))
			throw InvalidArgumentException("Wrong number of arguments for filter " + filtername, AT);
	} catch(std::exception& e){
		std::cerr << "[E] While processing arguments for filter " << filtername << std::endl;
		throw;
	}
}

void FilterAlgorithms::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)
{
	std::vector<std::string> vecArgs_new;
	bool alreadyDefined = false;

	for (unsigned int ii=0; ii<vecArgs_in.size(); ii++){
		//Check for "left", "right" or "center" keywords
		if ((vecArgs_in[ii] == "center") || (vecArgs_in[ii] == "left") || (vecArgs_in[ii] == "right")){
			if (alreadyDefined)
				throw InvalidArgumentException("Invalid or redefined arguments for filter " + filtername, AT);
			windowposition = vecArgs_in[ii];
			alreadyDefined = true;
		} else {
			vecArgs_new.push_back(vecArgs_in[ii]);
		}
	}

	parseFilterArguments(filtername, vecArgs_new, minArgs, maxArgs, isSoft, vecArgs_out);
}

bool FilterAlgorithms::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)
{
	vecWindow.clear();
	bool isSoft = false;
	std::string windowposition = "center"; //the default is a centered window
	std::vector<double> vecArgs;
	parseWindowFilterArguments(filtername, _vecArgs, 2, 2, isSoft, windowposition, vecArgs);

	if (vecArgs[0] < 1) //the window size has to be at least 1
		throw InvalidArgumentException("Number of data points in window of " +filtername+ " filter cannot be < 1", AT);
	if (vecArgs[1] < 0) //the time window has to be at least 0 seconds
		throw InvalidArgumentException("Time span of window for filter " +filtername+ " cannot be < 0 seconds", AT);

	//Deal with first parameter: minimal number of data points
	unsigned int windowSize = (unsigned int)vecArgs[0];
	unsigned int increment = 0;
	if (windowposition=="right"){
		increment = (unsigned int)windowSize - 1; //window reaches to the right
		if ((date != vecM[pos].date) && (increment>0)) increment--;
	} else if (windowposition=="left") {
		increment = 0;                        //window reaches to the left
	} else {
		increment = (unsigned int)(windowSize/2);                   //window is centered
		if (increment>0) increment--; //reason: the element at pos might be of greater date
	}

	unsigned int startposition = pos + increment;

	if (startposition > (vecM.size()-1)){
		if (!isSoft) return false; //if "soft" is not defined then there have to be enough data elements to the right
		startposition = (vecM.size()-1);
	}

	unsigned int endposition = 0;
	if (startposition < (windowSize-1)){
		return false; //not enough data points available
	} else {
		endposition = startposition - (windowSize - 1);
	}

	//Now deal with the second argument: the time window
	Date deltatime(vecArgs[1]/(24.*3600.)); //making a julian date out of the argument given in seconds
	while((deltatime > (vecM[startposition].date - vecM[endposition].date)) && (deltatime != 0.0)){
		//The time window is too small, so let's try enlargening it
		if ((windowposition == "right") && (!isSoft)){ //only expand to the right
			if (startposition<(vecM.size()-1)) startposition++;
			else return false; //time window not big enough
		} else if ((windowposition == "right") && (isSoft)){
			if (startposition<(vecM.size()-1)) startposition++;
			else if (endposition > 0) endposition--;
			else return false; //time window not big enough
		} else if ((windowposition == "left") && (!isSoft)){ //only expand to the left
			if (endposition > 0) endposition--;
			else return false; //time window not big enough
		} else if ((windowposition == "left") && (isSoft)){
			if (endposition > 0) endposition--;
			else if (startposition<(vecM.size()-1)) startposition++;
			else return false; //time window not big enough
		} else if ((windowposition == "center") && (!isSoft)){ //expand to left and right (equally)
			if ((endposition+startposition)%2 == 0) { //try to alternate when broadening window
				if (endposition > 0) endposition--;
				else return false; //time window not big enough
			} else {
				if (startposition<(vecM.size()-1)) startposition++;
				else return false; //time window not big enough
			}
		} else if ((windowposition == "center") && (isSoft)){ //expand to left and right (whereever possible)
			if ((endposition+startposition)%2 == 0) { //try to alternate when broadening window
				if (endposition > 0) endposition--;
				else if (startposition<(vecM.size()-1)) startposition++;
				else return false; //time window not big enough
			} else {
				if (startposition<(vecM.size()-1)) startposition++;
				else if (endposition > 0) endposition--;
				else return false; //time window not big enough
			}
		}
	}

	//Date gap(vecM[startposition].date - vecM[endposition].date);
	//cout << "The final gap is " << gap << "  windowposition: " << windowposition << endl;

	//Push all relevant data elements into the vector vecWindow
	//cout << "Start: " << startposition << "  End: " << endposition << endl;
	for (unsigned int ii=startposition+1; (ii--) > endposition; ){
		const double& tmp = vecM[ii].param(paramindex);
		if (tmp != IOUtils::nodata) {
			vecWindow.push_back(tmp);
			if (vecDate != NULL) (*vecDate).push_back(vecM[ii].date);
		}
		//cout << ii << ": pushed at vecM[" <<  ii << "] " << " : "  << endl;
	}

	return true;
}



/******************************************************************************
 * The following functions are implementations of different filter algorithms *
 ******************************************************************************/

/**
 * @brief Exponential smooting filter, exponential moving average
 * s_0 = x_0
 * s_n = alpha*x_(t-1) + (1-alpha)*s_t-1
 * - http://en.wikipedia.org/wiki/Exponential_smoothing
 * - alpha needs to be provided as argument, as well as the window size and centering 
 * - nodata values are excluded from the moving average calculation
 * - Three arguments expected (all have to be present and valid for the filter to start operating):
 *   - minimal number of points in window
 *   - minimal time interval spanning the window (in seconds)
 *   - alpha value 0 < alpha < 1
 * - the arguments may be preceded by the keywords "left", "center" or "right", indicating the window position
 * - the keyword "soft" maybe added at the beginning, if the window position is allowed to be adjusted to the data present
 * @code
 *          TA::filter1 = exp_smoothing
 *          TA::arg1    = right 1 1800 0.8 (1800 seconds time span for the strictly right leaning window), alpha=0.8
 *          RH::filter1 = mean_avg
 *          RH::arg1    = 10 600 0.6 (strictly left window spanning 600 seconds and at least 10 points), alpha=0.6
 * @endcode
 */
271
272
273
void FilterAlgorithms::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)
274
{
275
276
	(void)vecS; (void)vecWindowS;
	if (vecArgs.size() < 3)
277
278
		throw InvalidArgumentException("Wrong number of arguments for ExpSmoothingFilter", AT);

279
	vector<string> myArgs = vecArgs;
280
281
282
283
284
285
286

	//Take away the last argument (alpha value) and check validity
	double alpha = 0.5;
	IOUtils::convertString(alpha, myArgs.at(myArgs.size()-1)); 
	if ((alpha <= 0) || (alpha >= 1))
		throw InvalidArgumentException("ExpSmoothingFilter: alpha must be in [0;1]", AT);

287
288
	myArgs.pop_back(); //delete alpha from myArgs

289
290
291

	bool isSoft = false;
	std::string windowposition = "center"; //the default is a centered window
292
293
	std::vector<double> doubleArgs;
	parseWindowFilterArguments("exp_smoothing", myArgs, 2, 2, isSoft, windowposition, doubleArgs);
294

295
296
297
298
	//for every element in vecWindowM, get the Window and perform exponential smoothing
	for (unsigned int ii=0; ii<vecWindowM.size(); ii++){
		vector<MeteoData> vecTmpWindow;
		unsigned int position = IOUtils::seek(vecWindowM[ii].date, vecM);
299
300

		if (position != IOUtils::npos){
301
302
303
			unsigned int posfind = getWindowData("exp_smoothing", vecM, position, myArgs, vecTmpWindow);
			if (posfind == IOUtils::npos)
				continue;
304
305

			if (windowposition == "left"){
306
307
308
				vecTmpWindow.erase(vecTmpWindow.begin()+posfind+1, vecTmpWindow.end()); //delete all after posfind
				vecWindowM[ii].param(paramindex) = ExpSmoothingAlgorithm(vecTmpWindow, paramindex, alpha);
				//cout << "ExpSmoothing: " << vecWindowM[ii].param(paramindex) << endl;
309
			} else if (windowposition == "right"){
310
311
312
313
				vecTmpWindow.erase(vecTmpWindow.begin(), vecTmpWindow.begin()+posfind); //delete all before posfind
				std::reverse(vecTmpWindow.begin(), vecTmpWindow.end()); //reverse the vector, posfind most significant
				vecWindowM[ii].param(paramindex) = ExpSmoothingAlgorithm(vecTmpWindow, paramindex, alpha);
				//cout << "ExpSmoothing: " << vecWindowM[ii].param(paramindex) << endl;
314
			} else { //centered window - regroup according to time difference with posfind
315
316
317
318
319
				for (unsigned int jj=0; jj<vecTmpWindow.size(); jj++)
					vecTmpWindow[jj].date=Date(abs(vecWindowM[ii].date.getJulianDate() - vecTmpWindow[jj].date.getJulianDate()));
				std::sort(vecTmpWindow.begin(), vecTmpWindow.end(), compareMeteoData);
				vecWindowM[ii].param(paramindex) = ExpSmoothingAlgorithm(vecTmpWindow, paramindex, alpha);
				//cout << "ExpSmoothing: " << vecWindowM[ii].param(paramindex) << endl;
320
			}
321
		}	
322
323
324
325
326
327
328
329
	}
}

bool FilterAlgorithms::compareMeteoData (const MeteoData& m1, const MeteoData& m2)
{ 
	return (m1.date>m2.date);
}

330
331
332
unsigned int FilterAlgorithms::getWindowData(const std::string& filtername, const std::vector<MeteoData>& vecM, 
						    const unsigned int& pos, 
						    const std::vector<std::string>& _vecArgs, std::vector<MeteoData>& vecResult)
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
{
	Date date(vecM.at(pos).date);
	vecResult.clear();
	bool isSoft = false;
	std::string windowposition = "center"; //the default is a centered window
	std::vector<double> vecArgs;
	parseWindowFilterArguments(filtername, _vecArgs, 2, 2, isSoft, windowposition, vecArgs);

	if (vecArgs[0] < 1) //the window size has to be at least 1
		throw InvalidArgumentException("Number of data points in window of " +filtername+ " filter cannot be < 1", AT);
	if (vecArgs[1] < 0) //the time window has to be at least 0 seconds
		throw InvalidArgumentException("Time span of window for filter " +filtername+ " cannot be < 0 seconds", AT);

	//Deal with first parameter: minimal number of data points
	unsigned int windowSize = (unsigned int)vecArgs[0];
	unsigned int increment = 0;
	if (windowposition=="right"){
		increment = (unsigned int)windowSize - 1; //window reaches to the right
		if ((date != vecM[pos].date) && (increment>0)) increment--;
	} else if (windowposition=="left") {
		increment = 0;                        //window reaches to the left
	} else {
		increment = (unsigned int)(windowSize/2);                   //window is centered
		if (increment>0) increment--; //reason: the element at pos might be of greater date
	}

	unsigned int startposition = pos + increment;

	if (startposition > (vecM.size()-1)){
362
		if (!isSoft) return IOUtils::npos; //if "!isSoft" then there have to be enough data elements to the right
363
364
365
366
367
		startposition = (vecM.size()-1);
	}

	unsigned int endposition = 0;
	if (startposition < (windowSize-1)){
368
		return IOUtils::npos; //not enough data points available
369
370
371
372
373
374
375
376
377
378
	} else {
		endposition = startposition - (windowSize - 1);
	}

	//Now deal with the second argument: the time window
	Date deltatime(vecArgs[1]/(24.*3600.)); //making a julian date out of the argument given in seconds
	while((deltatime > (vecM[startposition].date - vecM[endposition].date)) && (deltatime != 0.0)){
		//The time window is too small, so let's try enlargening it
		if ((windowposition == "right") && (!isSoft)){ //only expand to the right
			if (startposition<(vecM.size()-1)) startposition++;
379
			else return IOUtils::npos; //time window not big enough
380
381
382
		} else if ((windowposition == "right") && (isSoft)){
			if (startposition<(vecM.size()-1)) startposition++;
			else if (endposition > 0) endposition--;
383
			else return IOUtils::npos; //time window not big enough
384
385
		} else if ((windowposition == "left") && (!isSoft)){ //only expand to the left
			if (endposition > 0) endposition--;
386
			else return IOUtils::npos; //time window not big enough
387
388
389
		} else if ((windowposition == "left") && (isSoft)){
			if (endposition > 0) endposition--;
			else if (startposition<(vecM.size()-1)) startposition++;
390
			else return IOUtils::npos; //time window not big enough
391
392
393
		} else if ((windowposition == "center") && (!isSoft)){ //expand to left and right (equally)
			if ((endposition+startposition)%2 == 0) { //try to alternate when broadening window
				if (endposition > 0) endposition--;
394
				else return IOUtils::npos; //time window not big enough
395
396
			} else {
				if (startposition<(vecM.size()-1)) startposition++;
397
				else return IOUtils::npos; //time window not big enough
398
399
400
401
402
			}
		} else if ((windowposition == "center") && (isSoft)){ //expand to left and right (whereever possible)
			if ((endposition+startposition)%2 == 0) { //try to alternate when broadening window
				if (endposition > 0) endposition--;
				else if (startposition<(vecM.size()-1)) startposition++;
403
				else return IOUtils::npos; //time window not big enough
404
405
406
			} else {
				if (startposition<(vecM.size()-1)) startposition++;
				else if (endposition > 0) endposition--;
407
				else return IOUtils::npos; //time window not big enough
408
409
410
411
412
413
414
415
416
			}
		}
	}

	//Date gap(vecM[startposition].date - vecM[endposition].date);
	//cout << "The final gap is " << gap << "  windowposition: " << windowposition << endl;

	//Push all relevant data elements into the vector vecResult
	//cout << "POS" << pos << "  start: " << startposition << "  end: " << endposition << endl;
417
418
	unsigned int posofposition = IOUtils::npos;
	unsigned int counter = 0;
419
420
	for (unsigned int ii=endposition; ii<startposition+1; ii++){
		vecResult.push_back(vecM[ii]);
421
422
		if (date == vecM[ii].date) posofposition=counter;
		counter++;
423
424
425
		//cout << ii << ": pushed at vecM[" <<  ii << "] " << " : "  << endl;
	}

426
	return posofposition;
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
}

/**
 * @brief Actual implementation of the smoothing algorithm
 * @param vecMeteo A vector of MeteoData
 * @param paramindex The meteo data parameter to be smoothed
 * @param alpha The alpha for the exponential smoothing, the smoothing factor
 * @return The smoothing value for the specified meteo parameter of the last element in vecMeteo
 */
double FilterAlgorithms::ExpSmoothingAlgorithm(const std::vector<MeteoData>& vecMeteo, const unsigned int& paramindex,
                                               const double& alpha)
{
	if (vecMeteo.size() == 0)
		return IOUtils::nodata;

	bool initCompleted = false;
	double expavg = IOUtils::nodata;
	for (unsigned int ii=0; ii<vecMeteo.size(); ii++){
		const double& currentval = vecMeteo[ii].param(paramindex);
		if (currentval != IOUtils::nodata){
			if (!initCompleted){
				expavg = currentval;
				initCompleted = true;
			} else {
				expavg = alpha*currentval + (1-alpha)*expavg;
			}
		}
	}

	return expavg;
}

/**
 * @brief Rate of change filter.
 * Calculate the change rate (ie: slope) between two points, if it is above a user given value, reject the point. Remarks:
 * - the maximum permissible rate of change (per seconds) has to be provided as an argument
 * @code
 * TA::filter1	= rate
 * TA::arg1	= 0.01
 * @endcode
 */
468
469
470
void FilterAlgorithms::RateFilter(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)
471
{
472
473
	(void)vecS; (void)vecWindowS;

474
475
	//parse arguments and check whether they are valid
	bool isSoft = false;
476
477
	std::vector<double> doubleArgs;
	parseFilterArguments("rate", vecArgs, 1, 1, isSoft, doubleArgs);
478
479
480

	//disable the warnings triggered by the commented out section that follows
	(void)vecM; (void)paramindex; (void)vecWindowM;
481
482
	/*
	const double& maxRateOfChange = doubleArgs[0];
483
484

	//Run actual Rate filter over all relevant meteo data
485
486
487
488
489
490
491
492
493
	for (unsigned int ii=2; ii<vecWindowM.size(); ii++){
		MeteoData& current     = vecWindowM[ii];
		MeteoData& previous    = vecWindowM[ii-1];
		MeteoData& preprevious = vecWindowM[ii-2];
	}


	
	if(vecWindowM.size()>1) {
494
		//the request time step is NOT part of the data, it will have to be resampled
495
496
		const double prev_value = vecM[pos-vecWindowM.size()].param(paramindex);
		const double prev_time = vecM[pos-vecWindowM.size()].date.getJulianDate()*24.*3600.; //in seconds
497

498
499
500
501
		for(unsigned int ii=0; ii<vecWindowM.size(); ii++){
			double& value = vecWindowM[ii].param(paramindex);
			const double curr_value = vecWindowM[ii].param(paramindex);
			const double curr_time = vecWindowM[ii].date.getJulianDate()*24.*3600.; //in seconds
502
503
			const double local_rate = abs((curr_value-prev_value)/(curr_time-prev_time));

504
			if( local_rate > doubleArgs[0] ) {
505
506
507
508
				value = IOUtils::nodata;
			}
		}
	} else {
509
		//const double tmp = vecWindowM.size();
510
511
512
		//std::cout << tmp << std::endl;
		//the request time step is part of the data
		if(pos>1) { //if we are at the start of the data set, we can not apply the filter...
513
			double& value = vecWindowM[0].param(paramindex);
514
515
516
517
518
519
			const double curr_value = vecM[pos].param(paramindex);
			const double curr_time = vecM[pos].date.getJulianDate()*24.*3600.; //in seconds
			const double prev_value = vecM[pos-1].param(paramindex);
			const double prev_time = vecM[pos-1].date.getJulianDate()*24.*3600.; //in seconds
			const double local_rate = abs((curr_value-prev_value)/(curr_time-prev_time));

520
			if( local_rate > doubleArgs[0] ) {
521
522
523
524
525
526
				value = IOUtils::nodata;
			}
		} else {
			//the filter can not be applied
		}
	}
527
	*/
528
529
530
531
532
533
534
535
536
537
538
539
540
}

/**
 * @brief Min/Max range filter.
 * Reject all values greater than the max or smaller than the min. Remarks:
 * - two arguments have to be provided, min and max (in SI)
 * - the keyword "soft" maybe added, in such a case all data greater than the max would be assigned
 * the maximum permissible value and all data smaller than the min would be assigned the minimum permissible value
 * @code
 * TA::filter1	= min_max
 * TA::arg1	= 230 330
 * @endcode
 */
541
542
543
void FilterAlgorithms::MinMaxFilter(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)
544
{
545
	(void)vecM; (void)vecS; (void)vecWindowS;
546
547
	//parse arguments and check whether they are valid
	bool isSoft = false;
548
549
	std::vector<double> doubleArgs;
	parseFilterArguments("min_max", vecArgs, 2, 2, isSoft, doubleArgs);
550

551
	sort(doubleArgs.begin(), doubleArgs.end()); //the two parameters are sorted ascending
552
553

	//Run actual MinMax filter over all relevant meteo data
554
555
	for(unsigned int ii=0; ii<vecWindowM.size(); ii++){
		double& value = vecWindowM[ii].param(paramindex);
556
557
558

		if (value == IOUtils::nodata) continue;

559
560
		if (value<doubleArgs[0]){
			if (isSoft) value=doubleArgs[0];
561
562
563
			else value=IOUtils::nodata;
			//cout << "Changed: " << value << endl;
		}
564
565
		if (value>doubleArgs[1]){
			if (isSoft) value=doubleArgs[1];
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
			else value=IOUtils::nodata;
			//cout << "Changed: " << value << endl;
		}
	}
}

/**
 * @brief Min range filter.
 * Reject all values smaller than the min. Remarks:
 * - the minimum permissible value has to be provided has an argument (in SI)
 * - the keyword "soft" maybe added, in such a case all data smaller than the min would be assigned
 * the minimum permissible value
 * @code
 * TA::filter1	= min
 * TA::arg1	= 230
 * @endcode
 */
583
584
585
void FilterAlgorithms::MinValueFilter(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)
586
{
587
	(void)vecM; (void)vecS; (void)vecWindowS;
588
589
	//parse arguments and check whether they are valid
	bool isSoft = false;
590
591
	std::vector<double> doubleArgs;
	parseFilterArguments("min", vecArgs, 1, 1, isSoft, doubleArgs);
592
593

	//Run actual MinValue filter over all relevant meteo data
594
595
596
	for(unsigned int ii=0; ii<vecWindowM.size(); ii++){
		double& value = vecWindowM[ii].param(paramindex);

597
		if (value == IOUtils::nodata) continue;
598
599
600

		if (value<doubleArgs[0]){
			if (isSoft) value=doubleArgs[0];
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
			else value=IOUtils::nodata;
		}
	}
}

/**
 * @brief Max range filter.
 * Reject all values greater than the max. Remarks:
 * - the maximum permissible value has to be provided has an argument (in SI)
 * - the keyword "soft" maybe added, in such a case all data greater than the max would be assigned
 * the maximum permissible value
 * @code
 * TA::filter1	= max
 * TA::arg1	= 330
 * @endcode
 */
617
618
619
void FilterAlgorithms::MaxValueFilter(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)
620
{
621
	(void)vecM; (void)vecS; (void)vecWindowS;
622
623
	//parse arguments and check whether they are valid
	bool isSoft = false;
624
625
	std::vector<double> doubleArgs;
	parseFilterArguments("max", vecArgs, 1, 1, isSoft, doubleArgs);
626
627

	//Run actual MaxValue filter over all relevant meteo data
628
629
630
	for(unsigned int ii=0; ii<vecWindowM.size(); ii++){
		double& value = vecWindowM[ii].param(paramindex);

631
		if (value == IOUtils::nodata) continue;
632
633
634

		if (value>doubleArgs[0]){
			if (isSoft) value=doubleArgs[0];
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
			else value=IOUtils::nodata;
		}
	}
}

/**
 * @brief Median Absolute Deviation.
 * Values outside of median ± 3 σ_MAD are rejected. The σ_MAD is calculated as follow:\n
 * <center>\f$ \sigma_{MAD} = K \cdot \mathop{median_i} \left( \left| X_i - \mathop{median_j} ( X_j ) \right| \right) \f$ with \f$ K = \Phi^{-1}; \Phi = 0.6745 \f$ </center>\n
 * See http://en.wikipedia.org/wiki/Median_absolute_deviation
 * for more information about the Mean Absolute Deviation.
 * @code
 * Valid examples for the io.ini file:
 *          TA::filter1 = mad
 *          TA::arg1    = soft left 1 1800  (1800 seconds time span for the left leaning window)
 *          RH::filter1 = mad
 *          RH::arg1    = 10 600            (strictly centered window spanning 600 seconds and at least 10 points)
 * @endcode
 */
654
655
656
657
void FilterAlgorithms::MedianAbsoluteDeviationFilter(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)
658
{
659
660
661
	(void)vecS; (void)vecWindowS;

	for (unsigned int ii=0; ii<vecWindowM.size(); ii++){
662
663
664
665
		double& value = vecWindowM[ii].param(paramindex);
		if (value == IOUtils::nodata) //No need to run filter for nodata points
			continue;

666
		unsigned int pos = IOUtils::seek(vecWindowM[ii].date, vecM, false);
667

668
669
670
		std::vector<double> vecWindow;
		if (!getWindowData("mad", vecM, pos, vecWindowM[ii].date, vecArgs, paramindex, vecWindow))
			return; //Not enough data to meet user configuration
671
	
672
673
674
675
		//Calculate MAD
		const double K = 1. / 0.6745;
		double mad     = IOUtils::nodata;
		double median  = IOUtils::nodata;
676

677
678
679
680
681
682
		try {
			median = Interpol1D::getMedian(vecWindow);
			mad    = Interpol1D::getMedianAverageDeviation(vecWindow);
		} catch(exception& e){
			return;
		}
683
	
684
		double sigma = mad * K;
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699

		if( (value>(median + 3.*sigma)) || (value<(median - 3.*sigma)) ) {
			value = IOUtils::nodata;
		}
	}
}


/**
 * @brief Accumulation over a user given period.
 * The input data is accumulated over a given time interval (given as filter argument, in minutes).
 * This is for example needed for converting rain gauges measurements read every 10 minutes to
 * hourly precipitation measurements. Remarks:
 * - the accumulation period has to be provided as an argument (in seconds)
 * @code
700
701
 * HNW::filter1 = accumulate
 * HNW::arg1	 = 3600
702
703
 * @endcode
 */
704
705
706
void FilterAlgorithms::AccumulateProcess(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)
707
{
708
709
	(void)vecS; (void)vecWindowS;
	
710
711
	//parse arguments and check whether they are valid
	bool isSoft = false;
712
713
	std::vector<double> doubleArgs;
	parseFilterArguments("accumulate", vecArgs, 1, 1, isSoft, doubleArgs);
714

715
	Date deltatime(doubleArgs[0]/(24.*3600.)); //making a julian date out of the argument given in seconds
716

717
718
	for (unsigned int ii=0; ii<vecWindowM.size(); ii++){
		unsigned int pos = IOUtils::seek(vecWindowM[ii].date, vecM);
719

720
721
722
723
724
		if (pos != IOUtils::npos){
			unsigned int startpos = pos;
			double sum = 0.0;
			while((vecM[pos].date + deltatime) > vecM[startpos].date){
				const double& val = vecM[pos].param(paramindex);
725

726
727
				if (val != IOUtils::nodata)
					sum += val;
728

729
730
				if (pos > 0) pos--;
				else break;
731
			}
732
733
734
735
			vecWindowM[ii].param(paramindex) = sum;
			//cout << "sum: " << vecWindowM[ii].param(paramindex) << endl;
		} else {
			throw IOException("Could not find an element to start accumulation", AT);
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
		}
	}
}

/**
 * @brief Median averaging.
 * The median average filter returns the median value of all values within a user given time window. Remarks:
 * - nodata values are excluded from the median
 * - if there is an even number of window elements the arithmetic mean of the two central elements is used to calculate the median
 * - Two arguments expected (both have to be fullfilled for the filter to start operating):
 *   - minimal number of points in window
 *   - minimal time interval spanning the window (in seconds)
 * - the two arguments may be preceded by the keywords "left", "center" or "right", indicating the window position
 * - the keyword "soft" maybe added, if the window position is allowed to be adjusted to the data present
 *
 * @code
 * Valid examples for the io.ini file:
 *          TA::filter1 = median_avg
 *          TA::arg1    = soft left 1 1800  (1800 seconds time span for the left leaning window)
 *          RH::filter1 = median_avg
 *          RH::arg1    = 10 600            (strictly centered window spanning 600 seconds and at least 10 points)
 * @endcode
 */
759
760
761
void FilterAlgorithms::MedianAvgProcess(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)
762
{
763
	(void)vecS; (void)vecWindowS;
764
765
766
767
768
769
770
771
772
773
774
	//Remarks:
	//1) nodata values are not considered when calculating the median
	//2) if there is an even number of window elements the arithmetic mean is used to calculate the median
	//3) Two arguments expected (both have to be fullfilled for the filter to start operating):
	//   1. minimal number of points in window
	//   2. minimal time interval spanning the window
	//4) the two arguments may be preceded by the keywords "left", "center" or "right", indicating the window
	//   position
	//5) the keyword "soft" maybe added, if the window position is allowed to be adjusted to the data present

	//for every element in vecFilteredM, get the Window and calculate median average
775
	for (unsigned int ii=0; ii<vecWindowM.size(); ii++){
776
		std::vector<double> vecWindow;
777
		unsigned int position = IOUtils::seek(vecWindowM[ii].date, vecM);
778

779
		if (!getWindowData("median_avg", vecM, position, vecWindowM[ii].date, vecArgs, paramindex, vecWindow))
780
781
782
783
784
			continue; //Not enough data to meet user configuration

		double median = IOUtils::nodata;
		try {
			median = Interpol1D::getMedian(vecWindow);
785
			vecWindowM[ii].param(paramindex) = median;
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
		} catch(exception& e){
			continue; //the median calculation did not work out, filter is not applied, value unchanged
		}
	}
}

/**
 * @brief Mean averaging.
 * The mean average filter returns the mean value of all values within a user given time window. Remarks:
 * - nodata values are excluded from the mean
 * - Two arguments expected (both have to be fullfilled for the filter to start operating):
 *   - minimal number of points in window
 *   - minimal time interval spanning the window (in seconds)
 * - the two arguments may be preceded by the keywords "left", "center" or "right", indicating the window position
 * - the keyword "soft" maybe added, if the window position is allowed to be adjusted to the data present
 *
 * @code
 * Valid examples for the io.ini file:
 *          TA::filter1 = mean_avg
 *          TA::arg1    = soft left 1 1800 (1800 seconds time span for the left leaning window)
 *          RH::filter1 = mean_avg
 *          RH::arg1    = 10 600          (strictly centered window spanning 600 seconds and at least 10 points)
 * @endcode
 */
810
811
812
void FilterAlgorithms::MeanAvgProcess(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)
813
{
814
	(void)vecS; (void)vecWindowS;
815
816
817
818
819
820
821
822
823
	//Remarks:
	//1) nodata values are not considered when calculating the mean
	//2) Two arguments expected (both have to be fullfilled for the filter to start operating):
	//   1. minimal number of points in window
	//   2. minimal time interval spanning the window
	//3) the two arguments may be preceded by the keywords "left", "center" or "right", indicating the window
	//   position
	//4) the keyword "soft" maybe added, if the window position is allowed to be adjusted to the data present

824
825
	//for every element in vecWindowM, get the Window and calculate mean average
	for (unsigned int ii=0; ii<vecWindowM.size(); ii++){
826
		std::vector<double> vecWindow;
827
		unsigned int position = IOUtils::seek(vecWindowM[ii].date, vecM);
828

829
		if (!getWindowData("mean_avg", vecM, position, vecWindowM[ii].date, vecArgs, paramindex, vecWindow))
830
831
832
833
834
			continue; //Not enough data to meet user configuration

		double mean = IOUtils::nodata;
		try {
			mean = Interpol1D::arithmeticMean(vecWindow);
835
			vecWindowM[ii].param(paramindex) = mean;
836
837
838
839
840
841
842
843
844
		} catch(exception& e){
			continue; //the mean calculation did not work out, filter is not applied, value unchanged
		}
	}
}

/**
 * @brief Wind vector averaging.
 * This calculates the vector average over a user given time period. Each wind vector within this period
845
 * is added and the final sum is normalized by the number of vectors that have been added. Important information:
846
847
848
849
850
851
852
853
854
855
856
857
858
859
 * - nodata values are excluded from the mean
 * - Two arguments expected (both have to be fullfilled for the filter to start operating):
 *   - minimal number of points in window
 *   - minimal time interval spanning the window (in seconds)
 * - the two arguments may be preceded by the keywords "left", "center" or "right", indicating the window position
 * - the keyword "soft" maybe added, if the window position is allowed to be adjusted to the data present
 * @code
 * Valid examples for the io.ini file:
 *          VW::filter1 = wind_avg
 *          VW::arg1    = soft left 1 1800 (1800 seconds time span for the left leaning window)
 *          VW::filter1 = wind_avg
 *          VW::arg1    = 10 600          (strictly centered window spanning 600 seconds and at least 10 points)
 * @endcode
 */
860
861
862
void FilterAlgorithms::WindAvgProcess(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)
863
{
864
	(void)vecS; (void)vecWindowS; (void)paramindex;
865
866
867
868
869
870
871
872
873
	//Remarks:
	//1) nodata values are not considered when calculating the mean
	//2) Two arguments expected (both have to be fullfilled for the filter to start operating):
	//   1. minimal number of points in window
	//   2. minimal time interval spanning the window
	//3) the two arguments may be preceded by the keywords "left", "center" or "right", indicating the window
	//   position
	//4) the keyword "soft" maybe added, if the window position is allowed to be adjusted to the data present
	//
874
875
876
877
878
	for (unsigned int ii=0; ii<vecWindowM.size(); ii++){
		unsigned int pos = IOUtils::seek(vecWindowM[ii].date, vecM);
		
		if (pos == IOUtils::npos)
			continue;
879

880
881
882
883
		std::vector<double> vecWindowVW, vecWindowDW;
		std::vector<Date> vecDateVW, vecDateDW;
		if (!getWindowData("wind_avg", vecM, pos, vecWindowM[ii].date, vecArgs, MeteoData::VW, vecWindowVW, &vecDateVW))
			continue; //Not enough data to meet user configuration
884

885
886
		if (!getWindowData("wind_avg", vecM, pos, vecWindowM[ii].date, vecArgs, MeteoData::DW, vecWindowDW, &vecDateDW))
			continue; //Not enough data to meet user configuration
887

888
889
		if (vecWindowVW.size() != vecWindowDW.size()) //same amount of data points necessary
			continue;
890

891
892
		for (unsigned int jj=0; jj<vecDateVW.size(); jj++){ //The VW and DW data points have to correlate
			if (vecDateVW[jj] != vecDateDW[jj]) continue;
893
894
		}

895
896
897
898
		//Calculate mean
		double meanspeed     = IOUtils::nodata;
		double meandirection = IOUtils::nodata;
		unsigned int vecSize = vecWindowVW.size();
899

900
901
902
903
904
905
906
907
908
909
910
		if (vecSize == 0){
			continue; //only nodata values detected or other problem
		} else {
			//calculate ve and vn
			double ve=0.0, vn=0.0;
			for (unsigned int jj=0; jj<vecSize; jj++){
				ve += vecWindowVW[jj] * sin(vecWindowDW[jj] * M_PI / 180.); //turn into radians
				vn += vecWindowVW[jj] * cos(vecWindowDW[jj] * M_PI / 180.); //turn into radians
			}
			ve /= vecSize;
			vn /= vecSize;
911

912
913
914
915
916
917
918
			meanspeed = sqrt(ve*ve + vn*vn);
			meandirection = fmod( atan2(ve,vn) * 180. / M_PI + 360. , 360.); // turn into degrees [0;360)
		}

		vecWindowM[ii].param(MeteoData::VW) = meanspeed;
		vecWindowM[ii].param(MeteoData::DW) = meandirection;
	}
919
920
921
922
}

} //namespace