WSL/SLF GitLab Repository

FilterStdDev.cc 3.78 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
30
31
32
33
34
35
36
37
/***********************************************************************************/
/*  Copyright 2011 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/meteofilters/FilterStdDev.h>
#include <cmath>

using namespace std;

namespace mio {

const double FilterStdDev::sigma = 2.; ///<How many times the stddev allowed for valid points

FilterStdDev::FilterStdDev(const std::vector<std::string>& vec_args) : WindowedFilter("STD_DEV")
{
	parse_args(vec_args);

	//This is safe, but maybe too imprecise:
	properties.time_before = min_time_span;
	properties.time_after  = min_time_span;
	properties.points_before = min_data_points;
	properties.points_after = min_data_points;
}

38
void FilterStdDev::process(const unsigned int& param, const std::vector<MeteoData>& ivec,
39
40
                           std::vector<MeteoData>& ovec)
{
41
42
	ovec = ivec;
	for (size_t ii=0; ii<ovec.size(); ii++){ //for every element in ivec, get a window
43
44
		double& value = ovec[ii](param);
		if(value==IOUtils::nodata) continue;
45
46
47
48

		//Calculate deviation
		double mean     = IOUtils::nodata;
		double std_dev  = IOUtils::nodata;
49

50
		size_t start, end;
51
52
		if( get_window_specs(ii, ivec, start, end) ) {
			getStat(ivec, param, start, end, std_dev, mean);
53
54
55
		}  else if(!is_soft) {
			value = IOUtils::nodata;
			continue;
56
57
		}

58
59
60
61
		if(mean==IOUtils::nodata) {
			if(!is_soft) value = IOUtils::nodata;
			continue;
		}
62

63
64
		if( value!=IOUtils::nodata && abs(value-mean)>sigma*std_dev) {
			value = IOUtils::nodata;
65
66
67
68
		}
	}
}

69
70
71
void FilterStdDev::getStat(const std::vector<MeteoData>& ivec, const unsigned int& param,
                           const size_t& start, const size_t& end, double& stddev, double& mean)
{
72
	size_t count=0;
73
	double sum=0.;
74

75
76
	for(size_t ii=start; ii<=end; ii++) {
		const double& value = ivec[ii](param);
77
78
79
80
81
82
83
84
85
86
		if(value!=IOUtils::nodata) {
			sum += value;
			count++;
		}
	}

	if(count<=1) {
		mean = IOUtils::nodata;
		stddev = IOUtils::nodata;
	} else {
87
		//compensated variance algorithm, see https://secure.wikimedia.org/wikipedia/en/wiki/Algorithms_for_calculating_variance
88
		mean = sum/(double)count;
89
		double sum2=0., sum3=0.;
90
91
		for(size_t ii=start; ii<=end; ii++) {
			const double& value = ivec[ii](param);
92
			if(value!=IOUtils::nodata) {
93
94
95
				const double delta = value - mean;
				sum2 += delta*delta;
				sum3 += delta;
96
97
			}
		}
98
		const double variance = (sum2 - sum3*sum3/static_cast<double>(count)) / static_cast<double>(count - 1);
99
100
101
102
103
104
105
106
107
		stddev = sqrt(variance);
	}
}

void FilterStdDev::parse_args(std::vector<std::string> vec_args)
{
	vector<double> filter_args;

	if (vec_args.size() > 2){
108
		is_soft = ProcessingBlock::is_soft(vec_args);
109
110
111
112
	}

	if (vec_args.size() > 2)
		centering = (WindowedFilter::Centering)WindowedFilter::get_centering(vec_args);
113

114
	convert_args(2, 2, vec_args, filter_args);
115
116

	if ((filter_args[0] < 1) || (filter_args[1] < 0)){
117
		throw InvalidArgumentException("Invalid window size configuration for filter " + getName(), AT);
118
119
120
121
	}

	min_data_points = (unsigned int)floor(filter_args[0]);
	min_time_span = Duration(filter_args[1] / 86400.0, 0.);
122
}
123
124

}