WSL/SLF GitLab Repository

ProcessingBlock.cc 12.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/***********************************************************************************/
/*  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/>.
*/
18
19
20
21
22
23
#include <fstream>
#include <sstream>
#include <errno.h>
#include <cstring>

#include <meteoio/FileUtils.h>
24
25
26
27
28
29
30
#include <meteoio/meteoFilters/ProcessingBlock.h>
#include <meteoio/meteoFilters/FilterSuppr.h>
#include <meteoio/meteoFilters/FilterMin.h>
#include <meteoio/meteoFilters/FilterMax.h>
#include <meteoio/meteoFilters/FilterMinMax.h>
#include <meteoio/meteoFilters/FilterStdDev.h>
#include <meteoio/meteoFilters/FilterRate.h>
31
#include <meteoio/meteoFilters/FilterUnheatedPSUM.h>
32
33
#include <meteoio/meteoFilters/FilterTukey.h>
#include <meteoio/meteoFilters/FilterMAD.h>
34
#include <meteoio/meteoFilters/ProcAggregate.h>
35
#include <meteoio/meteoFilters/ProcIIR.h>
36
37
38
#include <meteoio/meteoFilters/ProcUndercatch_WMO.h>
#include <meteoio/meteoFilters/ProcUndercatch_Forland.h>
#include <meteoio/meteoFilters/ProcUndercatch_Hamon.h>
39
#include <meteoio/meteoFilters/ProcPSUMDistribute.h>
40
#include <meteoio/meteoFilters/ProcUnventilatedT.h>
41
#include <meteoio/meteoFilters/ProcShade.h>
42
43
#include <meteoio/meteoFilters/ProcAdd.h>
#include <meteoio/meteoFilters/ProcMult.h>
44
#include <meteoio/meteoFilters/ProcNoise.h>
45
46
#include <meteoio/meteoFilters/ProcExpSmoothing.h>
#include <meteoio/meteoFilters/ProcWMASmoothing.h>
47
#include <meteoio/meteoFilters/FilterNoChange.h>
48
49
#include <meteoio/meteoFilters/FilterTimeconsistency.h>
#include <meteoio/meteoFilters/FilterOffsetsnowdepth.h>
50
#include <meteoio/meteoFilters/FilterDeGrass.h>
51
52

namespace mio {
53
54
55
56
57
58
59
/**
 * @page processing Processing overview
 * The pre-processing infrastructure is described in ProcessingBlock (for its API). The goal of this page is to give an overview of the available filters and processing elements and their usage.
 *
 * @section processing_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.
 *
60
 * Several filter take arguments describing a processing window (for example, FilterStdDev). In such a case, two values are given:
61
 * - the minimum time span of the window
62
63
 * - the minimum number of points in the window
 *
64
65
 * The ProcessingBlock will walk through the data, starting at the current point and adding points to the processing window.  When both of these criterias are met,
 * the window is accepted. This means that a window defined as "6 21600" is a window that contains 6 points minimum AND spans at least 21600 seconds.
66
 *
67
68
69
70
71
72
73
74
 * @section processing_section Filtering section
 * The filters are specified for each parameter in the [Filters] section. This section contains
 * a list of the various meteo parameters (see MeteoData) 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
Mathias Bavay's avatar
Mathias Bavay committed
75
 *
76
77
78
79
 * RH::filter1	= min_max
 * RH::arg1	= -0.2 1.2
 * RH::filter2	= min_max
 * RH::arg2	= soft 0.0 1.0
Mathias Bavay's avatar
Mathias Bavay committed
80
 *
81
82
83
84
 * PSUM::filter1	= min
 * PSUM::arg1	= -0.1
 * PSUM::filter2	= min
 * PSUM::arg2	= soft 0.
85
86
87
 * @endcode
 *
 * @section processing_available Available processing elements
88
 * New filters can easily be developed. The filters that are currently available are the following:
89
90
 * - MIN: minimum check filter, see FilterMin
 * - MAX: maximum check filter, see FilterMax
91
92
 * - MIN_MAX: range check filter, see FilterMinMax
 * - RATE: rate of change filter, see FilterRate
Mathias Bavay's avatar
Mathias Bavay committed
93
 * - STD_DEV: reject data outside mean +/- k*stddev, see FilterStdDev
94
 * - MAD: median absolute deviation, see FilterMAD
95
 * - TUKEY: Tukey53H spike detection, based on median, see FilterTukey
96
 * - UNHEATED_RAINGAUGE: detection of snow melting in a rain gauge, see FilterUnheatedPSUM
97
 * - NO_CHANGE: reject data that changes too little (low variance), see FilterNoChange
98
99
 * - TIME_CONSISTENCY: reject data that changes too much , see FilterTimeconsistency
 * - DETECT_GRASS: detection of grass growing under the snow height sensor, see FilterDeGrass
100
 *
Mathias Bavay's avatar
Mathias Bavay committed
101
 * Some data transformations are also supported besides filtering, both very basic and generic data transformations:
102
 * - SUPPR: delete all or some data, see FilterSuppr
Mathias Bavay's avatar
Mathias Bavay committed
103
104
 * - ADD: adds a given offset to the data, see ProcAdd
 * - MULT: multiply the data by a given factor, see ProcMult
105
 * - NOISE: add to (or multiply by) randomly distributed noise, see ProcNoise
Mathias Bavay's avatar
Mathias Bavay committed
106
107
 *
 * As well as more specific data transformations:
108
 * - EXP_SMOOTHING: exponential smoothing of data, see ProcExpSmoothing
109
 * - WMA_SMOOTHING: weighted moving average smoothing of data, see ProcWMASmoothing
110
 * - IIR: Low Pass or High Pass critically damped filter, see ProcIIR
111
 * - AGGREGATE: various data aggregation algorithms, see ProcAggregate
112
 * - UNDERCATCH_WMO: WMO rain gauge correction for undercatch, using various correction models, see ProcUndercatch_WMO
113
 * - UNDERCATCH_FORLAND: Forland1996 rain gauge correction for solid and liquid undercatch, using various correction models, see ProcUndercatch_Forland
114
 * - UNDERCATCH_HAMON: Hamon1973 rain gauge correction for undercatch, see ProcUndercatch_Hamon
115
 * - UNVENTILATED_T: unventilated temperature sensor correction, see ProcUnventilatedT
116
 * - PSUM_DISTRIBUTE: distribute accumulated precipitation over preceeding timesteps, see ProcPSUMDistribute
117
 * - SHADE: apply a shading mask to the Incoming or Reflected Short Wave Radiation, see ProcShade
118
 *
119
 */
120

121
ProcessingBlock* BlockFactory::getBlock(const std::string& blockname, const std::vector<std::string>& vec_args, const Config& cfg)
122
{
123
124
125
126
	//the indenting is a little weird, this is in order to show the same groups as in the documentation above
	
	//normal filters
	if (blockname == "MIN"){
127
		return new FilterMin(vec_args, blockname);
128
	} else if (blockname == "MAX"){
129
		return new FilterMax(vec_args, blockname);
130
	} else if (blockname == "MIN_MAX"){
131
		return new FilterMinMax(vec_args, blockname);
132
	} else if (blockname == "RATE"){
133
		return new FilterRate(vec_args, blockname);
134
135
	} else if (blockname == "STD_DEV"){
		return new FilterStdDev(vec_args, blockname);
136
	} else if (blockname == "MAD"){
137
		return new FilterMAD(vec_args, blockname);
138
139
	} else if (blockname == "TUKEY"){
		return new FilterTukey(vec_args, blockname);
140
	} else if (blockname == "UNHEATED_RAINGAUGE"){
141
		return new FilterUnheatedPSUM(vec_args, blockname);
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
	} else if (blockname == "NO_CHANGE"){
		return new FilterNoChange(vec_args, blockname);
	} else if (blockname == "TIME_CONSISTENCY"){
		return new FilterTimeconsistency(vec_args, blockname);
	} else if (blockname == "OFFSNOW"){
		return new FilterOffsetsnowdepth(vec_args, blockname);
	} else if (blockname == "DETECT_GRASS"){
		return new FilterDeGrass(vec_args, blockname);
	} 
	
	//general data transformations
	else if (blockname == "SUPPR"){
		return new FilterSuppr(vec_args, blockname, cfg.getConfigRootDir(), cfg.get("TIME_ZONE", "Input"));
	} else if (blockname == "ADD"){
		return new ProcAdd(vec_args, blockname, cfg.getConfigRootDir());
	} else if (blockname == "MULT"){
		return new ProcMult(vec_args, blockname, cfg.getConfigRootDir());
	} else if (blockname == "NOISE"){
		return new ProcNoise(vec_args, blockname);
	}
	
	//more specific data transformations
	else if (blockname == "EXP_SMOOTHING"){
		return new ProcExpSmoothing(vec_args, blockname);
	} else if (blockname == "WMA_SMOOTHING"){
		return new ProcWMASmoothing(vec_args, blockname);
168
169
	} else if (blockname == "IIR"){
		return new ProcIIR(vec_args, blockname);
170
171
	} else if (blockname == "AGGREGATE"){
		return new ProcAggregate(vec_args, blockname);
172
	} else if (blockname == "UNDERCATCH_WMO"){
173
		return new ProcUndercatch_WMO(vec_args, blockname);
174
175
	} else if (blockname == "UNDERCATCH_FORLAND"){
		return new ProcUndercatch_Forland(vec_args, blockname);
176
	} else if (blockname == "UNDERCATCH_HAMON"){
177
		return new ProcUndercatch_Hamon(vec_args, blockname);
178
	} else if (blockname == "UNVENTILATED_T"){
179
		return new ProcUnventilatedT(vec_args, blockname);
180
181
	} else if (blockname == "PSUM_DISTRIBUTE"){
		return new ProcPSUMDistribute(vec_args, blockname);
182
	} else if (blockname == "SHADE"){
183
		return new ProcShade(vec_args, blockname, cfg);
184
	} else {
185
		throw IOException("The processing block '"+blockname+"' does not exist! " , AT);
186
187
188
189
	}

}

190
ProcessingBlock::ProcessingBlock(const std::string& name) : properties(), block_name(name)
191
{}
192

Mathias Bavay's avatar
Mathias Bavay committed
193
void ProcessingBlock::convert_args(const size_t& min_nargs, const size_t& max_nargs,
194
                               const std::vector<std::string>& vec_args, std::vector<double>& dbl_args) const
195
{
Mathias Bavay's avatar
Mathias Bavay committed
196
197
	const size_t nr_of_args = vec_args.size();
	if ((nr_of_args < min_nargs) || (nr_of_args > max_nargs))
198
199
		throw InvalidArgumentException("Wrong number of arguments for filter/processing element \"" + getName() + "\"", AT);

Mathias Bavay's avatar
Mathias Bavay committed
200
	for (size_t ii=0; ii<nr_of_args; ii++){
201
202
203
204
205
206
		double tmp = IOUtils::nodata;
		IOUtils::convertString(tmp, vec_args[ii]);
		dbl_args.push_back(tmp);
	}
}

207
bool ProcessingBlock::is_soft(std::vector<std::string>& vec_args) {
208
209
210
	for (size_t ii=0; ii<vec_args.size(); ++ii) {
		if (IOUtils::strToUpper( vec_args[ii] ) == "SOFT") {
			vec_args.erase( vec_args.begin() + ii );
211
212
213
214
215
216
217
			return true;
		}
	}

	return false;
}

218
std::string ProcessingBlock::getName() const {
219
220
221
	return block_name;
}

222
223
void ProcessingBlock::readCorrections(const std::string& filter, const std::string& filename, const char& c_type, const double& init, std::vector<double> &corrections)
{
224
	std::ifstream fin(filename.c_str());
225
226
227
228
	if (fin.fail()) {
		std::ostringstream ss;
		ss << "Filter " << filter << ": ";
		ss << "error opening file \"" << filename << "\", possible reason: " << std::strerror(errno);
229
		throw AccessException(ss.str(), AT);
230
231
232
233
234
235
236
237
	}

	if (c_type=='m') corrections.resize(12, init);
	else if (c_type=='d') corrections.resize(366, init);
	else if (c_type=='h') corrections.resize(24, init);
	const size_t maxIndex = corrections.size();
	const size_t minIndex = (c_type=='h')? 0 : 1;

238
	const char eoln = FileUtils::getEoln(fin); //get the end of line character for the file
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
271
272
273
274
275
276
277

	try {
		size_t index, lcount=0;
		double value;
		do {
			lcount++;
			std::string line;
			getline(fin, line, eoln); //read complete line
			IOUtils::stripComments(line);
			IOUtils::trim(line);
			if (line.empty()) continue;

			std::istringstream iss(line);
			iss.setf(std::ios::fixed);
			iss.precision(std::numeric_limits<double>::digits10);
			iss >> std::skipws >> index;
			if ( !iss || index<minIndex || (index-minIndex)>=maxIndex) {
				std::ostringstream ss;
				ss << "Invalid index in file " << filename << " at line " << lcount;
				throw InvalidArgumentException(ss.str(), AT);
			}
			iss >> std::skipws >> value;
			if ( !iss ){
				std::ostringstream ss;
				ss << "Invalid value in file " << filename << " at line " << lcount;
				throw InvalidArgumentException(ss.str(), AT);
			}

			corrections.at( index-minIndex ) = value;
		} while (!fin.eof());
		fin.close();
	} catch (const std::exception&){
		if (fin.is_open()) {//close fin if open
			fin.close();
		}
		throw;
	}
}

278
279
ProcessingBlock::~ProcessingBlock() {}

280
const std::string ProcessingBlock::toString() const {
281
	std::ostringstream os;
282
283
	os << "[" << block_name << " ";
	os << properties.toString();
284
	os << "]";
285
	return os.str();
286
287
}

288
const ProcessingProperties& ProcessingBlock::getProperties() const {
289
290
291
	return properties;
}

292
293
const std::string ProcessingProperties::toString() const
{
294
	std::ostringstream os;
295
296
	const double h_before = time_before.getJulian()*24.;
	const double h_after = time_after.getJulian()*24.;
Mathias Bavay's avatar
Mathias Bavay committed
297
298
	const size_t p_before = points_before;
	const size_t p_after = points_after;
299
300

	os << "{";
301
302
303
	if(h_before>0. || h_after>0.) os << "-" << h_before << " +" << h_after << " h; ";
	if(p_before>0 || p_after>0) os << "-" << p_before << " +" << p_after << " pts; ";
	if(stage==ProcessingProperties::first)
304
		os << "p¹";
305
	if(stage==ProcessingProperties::second)
306
		os << "p²";
307
	if(stage==ProcessingProperties::both)
308
		os << "p½";
309
	os << "}";
310
	return os.str();
311
312
}

313
} //end namespace