WSL/SLF GitLab Repository

FilterTukey.cc 4.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
38
/***********************************************************************************/
/*  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/FilterTukey.h>
#include <meteoio/IOUtils.h>
#include <cmath>

using namespace std;

namespace mio {

const double FilterTukey::k = 1.5; ///<How many times the stddev allowed as deviation to the smooth signal for valid points

FilterTukey::FilterTukey(const std::vector<std::string>& vec_args) : WindowedFilter("TUKEY")
{
	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;
}

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

46
		size_t start, end;
47
48
49
		if( get_window_specs(ii, ivec, start, end) ) {
			//Calculate std deviation
			const double std_dev  = getStdDev(ivec, param, start, end);
50

51
52
53
54
55
			const double u3 = getU3(ivec, ii, param);
			if(std_dev!=IOUtils::nodata && u3!=IOUtils::nodata) {
				if( abs(value-u3) > k*std_dev ) {
					value = IOUtils::nodata;
				}
56
57
			} else if(!is_soft) value = IOUtils::nodata;
		} else if(!is_soft) value = IOUtils::nodata;
58
59
60
	}
}

61
62
double FilterTukey::getStdDev(const std::vector<MeteoData>& ivec, const unsigned int& param, const size_t& start, const size_t& end)
{
63
	size_t count=0;
64
65
	double sum=0.;

66
67
	for(size_t ii=start; ii<=end; ii++) {
		const double& value = ivec[ii](param);
68
69
70
71
72
73
74
75
76
77
78
79
80
		if(value!=IOUtils::nodata) {
			sum += value;
			count++;
		}
	}

	if(count<=1) {
		return IOUtils::nodata;
	}

	//compensated variance algorithm, see https://secure.wikimedia.org/wikipedia/en/wiki/Algorithms_for_calculating_variance
	const double mean = sum/(double)count;
	double sum2=0., sum3=0.;
81
82
	for(size_t ii=start; ii<=end; ii++) {
		const double& value = ivec[ii](param);
83
		if(value!=IOUtils::nodata) {
84
85
86
			const double delta = value - mean;
			sum2 += delta*delta;
			sum3 += delta;
87
88
		}
	}
89
	const double variance = (sum2 - sum3*sum3/static_cast<double>(count)) / static_cast<double>(count - 1);
90
91
92
93

	return sqrt(variance);
}

94
double FilterTukey::getU3(const std::vector<MeteoData>& ivec, const size_t& i, const unsigned int& param)
95
{
96
	//exit if we don't have the required data points
97
	if( i<4 || i>=(ivec.size()-4) ) {
98
99
100
101
102
		return IOUtils::nodata;
	}

	//prepare intermediate variances
	std::vector<double> u2;
103
	for(char ii=-1; ii<=1; ii++) {
104
		std::vector<double> u1;
105
		for(char jj=-1; jj<=1; jj++) {
106
			std::vector<double> u;
107
108
			for(char kk=-2; kk<=2; kk++) {
				const size_t index = (i + (kk + jj + ii));
109
				const double value = ivec[index](param);
110
111
112
				if(value!=IOUtils::nodata)
					u.push_back( value );
			}
113
			if(!u.empty())
114
115
				u1.push_back( Interpol1D::getMedian(u) );
		}
116
		if(!u1.empty())
117
118
119
120
121
122
123
124
			u2.push_back( Interpol1D::getMedian(u1) );
		else
			u2.push_back( IOUtils::nodata );
	}

	//compute the variance u3
	//u3 = 1/4*( u2[0] + 2.*u2[1] + u2[2] )
	double u3=0.;
125
	size_t count=0;
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
	if(u2[0]!=IOUtils::nodata) {
		u3 += u2[0];
		count++;
	}
	if(u2[1]!=IOUtils::nodata) { //current timestep
		u3 += u2[1]*2.;
		count += 2;
	}
	if(u2[2]!=IOUtils::nodata) {
		u3 += u2[2];
		count++;
	}

	if(count>0)
		return u3/((double)count);
	else
		return IOUtils::nodata;
}

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

	if (vec_args.size() > 2){
150
		is_soft = ProcessingBlock::is_soft(vec_args);
151
152
153
154
	}

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

156
	convert_args(2, 2, vec_args, filter_args);
157
158

	if ((filter_args[0] < 1) || (filter_args[1] < 0)){
159
		throw InvalidArgumentException("Invalid window size configuration for filter " + getName(), AT);
160
161
	}

162
	min_data_points = (size_t)floor(filter_args[0]);
163
	min_time_span = Duration(filter_args[1] / 86400.0, 0.);
164
}
165
166

}