WSL/SLF GitLab Repository

ResamplingAlgorithms.cc 38.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
/***********************************************************************************/
/*  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/ResamplingAlgorithms.h>
19
#include <meteoio/MathOptim.h>
20
21
22
#include <meteoio/meteoLaws/Atmosphere.h>
#include <meteoio/meteoLaws/Sun.h>
#include <meteoio/meteoStats/libinterpol1D.h>
23
#include <meteoio/meteoFilters/ProcPSUMDistribute.h> //for the precipitation distribution
24

25
#include <cmath>
26
#include <algorithm>
27
28
29
30
31

using namespace std;

namespace mio {

32
ResamplingAlgorithms* ResamplingAlgorithmsFactory::getAlgorithm(const std::string& i_algoname, const std::string& parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs)
33
34
35
{
	const std::string algoname(IOUtils::strToUpper(i_algoname));

36
	if (algoname == "NONE" || algoname == "NO") {
37
		return new NoResampling(algoname, parname, dflt_window_size, vecArgs);
38
	} else if (algoname == "LINEAR") {
39
		return new LinearResampling(algoname, parname, dflt_window_size, vecArgs);
40
	} else if (algoname == "NEAREST") {
41
		return new NearestNeighbour(algoname, parname, dflt_window_size, vecArgs);
42
	} else if (algoname == "ACCUMULATE") {
43
		return new Accumulate(algoname, parname, dflt_window_size, vecArgs);
44
45
	} else if (algoname == "SOLAR") {
		return new Solar(algoname, parname, dflt_window_size, vecArgs);
46
	} else if (algoname == "DAILY_SOLAR") {
47
		return new Daily_solar(algoname, parname, dflt_window_size, vecArgs);
48
49
	} else if (algoname == "DAILY_AVG") {
		return new DailyAverage(algoname, parname, dflt_window_size, vecArgs);
50
51
52
53
	} else {
		throw IOException("The resampling algorithm '"+algoname+"' is not implemented" , AT);
	}
}
54

55
56
57
//compute the partial accumulation at the left of curr_date within a sampling interval
double ResamplingAlgorithms::partialAccumulateAtLeft(const std::vector<MeteoData>& vecM, const size_t& paramindex,
                                                     const size_t& pos, const Date& curr_date)
58
{
59
	const size_t end = pos+1;
60
	if (end>=vecM.size()) return IOUtils::nodata; //reaching the end of the input vector
61
62
63
64
65
66

	const double valend = vecM[end](paramindex);
	if (valend == IOUtils::nodata) return IOUtils::nodata;

	const double jul1 = vecM[pos].date.getJulian(true);
	const double jul2 = vecM[end].date.getJulian(true);
67

68
	const double left_accumulation = linearInterpolation(jul1, 0., jul2, valend, curr_date.getJulian(true));
69

70
71
72
73
74
75
76
	return left_accumulation;
}

//compute the partial accumulation at the right of curr_date within a sampling interval
double ResamplingAlgorithms::partialAccumulateAtRight(const std::vector<MeteoData>& vecM, const size_t& paramindex,
                                                      const size_t& pos, const Date& curr_date)
{
77
	const size_t end = pos+1;
78
	if (end>=vecM.size()) return IOUtils::nodata; //reaching the end of the input vector
79
80
81
82
83
84

	const double valend = vecM[end](paramindex);
	if (valend == IOUtils::nodata) return IOUtils::nodata;

	const double jul1 = vecM[pos].date.getJulian(true);
	const double jul2 = vecM[end].date.getJulian(true);
85

86
	const double left_accumulation = linearInterpolation(jul1, 0., jul2, valend, curr_date.getJulian(true));
87

88
	return valend - left_accumulation;
89
90
}

91
92
93
94
95
/**
 * @brief This function returns the last and next valid points around a given position
 * @param pos current position (index)
 * @param paramindex meteo parameter to use
 * @param vecM vector of MeteoData
96
 * @param resampling_date date to resample
97
98
99
100
101
 * @param window_size size of the search window
 * @param indexP1 index of point before the current position (IOUtils::npos if none could be found)
 * @param indexP2 index of point after the current position (IOUtils::npos if none could be found)
 */
void ResamplingAlgorithms::getNearestValidPts(const size_t& pos, const size_t& paramindex, const std::vector<MeteoData>& vecM, const Date& resampling_date,
102
                                              const double& window_size, size_t& indexP1, size_t& indexP2)
103
{
104
105
	indexP1=IOUtils::npos;
	indexP2=IOUtils::npos;
106

107
	const Date dateStart = resampling_date - window_size;
108
	for (size_t ii=pos; ii-- >0; ) {
109
		if (vecM[ii].date < dateStart) break;
110
		if (vecM[ii](paramindex) != IOUtils::nodata) {
111
112
113
114
115
116
117
			indexP1 = ii;
			break;
		}
	}

	//make sure the search window remains window_size
	const Date dateEnd = (indexP1 != IOUtils::npos)? vecM[indexP1].date+window_size : resampling_date+window_size;
118
	for (size_t ii=pos; ii<vecM.size(); ++ii) {
119
120
121
122
123
124
125
126
127
128
129
130
131
132
		if (vecM[ii].date > dateEnd) break;
		if (vecM[ii](paramindex) != IOUtils::nodata) {
			indexP2 = ii;
			break;
		}
	}
}

/**
 * @brief This function solves the equation y = ax + b for two given points and returns y for a given x
 * @param x1 x-coordinate of first point
 * @param y1 y-coordinate of first point
 * @param x2 x-coordinate of second point
 * @param y2 y-coordinate of second point
133
 * @param x x-coordinate of desired point
134
135
136
 * @return y-coordinate of desired point
 */
double ResamplingAlgorithms::linearInterpolation(const double& x1, const double& y1,
137
                                       const double& x2, const double& y2, const double& x)
138
139
140
141
142
143
144
145
{
	if (x1 == x2)
		throw IOException("Attempted division by zero", AT);

	//Solving y = ax + b
	const double a = (y2 - y1) / (x2 - x1);
	const double b = y2 - a*x2;

146
	return (a*x + b);
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
/**
 * @brief For a given date, find the start of the day, considering that for midnight we return the day before!
 * (as is necessary for daily averages, sums, etc that can be provided at midnight for the day before)
 * @param resampling_date current date
 * @return start of the day or start of the day before in case of midnight
 */
Date ResamplingAlgorithms::getDailyStart(const Date& resampling_date)
{
	Date Start( resampling_date );
	Start.rnd(24*3600, Date::DOWN);
	if (Start==resampling_date) //if resampling_date=midnight GMT, the rounding lands on the exact same date
		Start -= 1.;
	
	return Start;
}

/**
 * @brief Find a unique value in a given time interval. 
 * This is useful for retrieving a unique daily average, daily sum, etc
 * @param vecM meteo data time series for the station
 * @param paramindex index of the meteo parameter to process
 * @param pos index of the current time step
 * @param intervalStart start of the interval within which we will look for a valid value (often, start of the day)
 * @param intervalEnd end of the interval within which we will look for a valid value (often, end of the day)
 * @return index of the parameter's valid value
 */
175
size_t ResamplingAlgorithms::getDailyValue(const std::vector<MeteoData>& vecM, const size_t& paramindex, size_t pos, const Date& intervalStart, const Date& intervalEnd)
176
177
178
{
	size_t indexP1=IOUtils::npos;
	size_t indexP2=IOUtils::npos;
179
180
181
	
	//if pos was not properly pre-positioned, do it
	if (vecM[pos].date<intervalStart) {
182
		for (; pos<vecM.size(); pos++)
183
184
185
			if (vecM[pos].date>=intervalStart) break;
	}
	if (vecM[pos].date>intervalEnd) {
186
187
		for (size_t ii = pos; ii-- >0; ) {
			pos=ii; //because ii gets corrupted at the final iteration if going all the way down
188
			if (vecM[pos].date<=intervalEnd) break;
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

	//look for daily sum before the current point
	for (size_t ii=pos; ii-- >0; ) {
		if (vecM[ii].date < intervalStart) break;
		if (vecM[ii](paramindex) != IOUtils::nodata) {
			indexP1 = ii;
			break;
		}
	}

	//look for daily sum after the current point
	for (size_t ii=pos; ii<vecM.size(); ++ii) {
		if (vecM[ii].date > intervalEnd) break;
		if (vecM[ii](paramindex) != IOUtils::nodata) {
			indexP2 = ii;
			break;
		}
	}
	
	if (indexP1!=IOUtils::npos && indexP2!=IOUtils::npos) {
		//if the data was provided at 00:00, a sum for each day has been found
		//we only keep the later one as being the sum of the past day
		int hour1, min1;
		vecM[indexP1].date.getTime(hour1, min1);
		if (hour1==0 && min1==0)
			indexP1=IOUtils::npos;
		else { //otherwise, this means multiple daily sums have been found for the same day
			const string msg = "More than one daily sum of solar radiation found between "+intervalStart.toString(Date::ISO)+" and "+intervalEnd.toString(Date::ISO);
			throw IOException(msg, AT);
		}
	}

	if (indexP1!=IOUtils::npos) return indexP1;
	if (indexP2!=IOUtils::npos) return indexP2;
	return IOUtils::npos;
}

228
/**********************************************************************************
229
 * The following functions are implementations of different resampling algorithms *
230
 **********************************************************************************/
231

232
233
234
NoResampling::NoResampling(const std::string& i_algoname, const std::string& i_parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs)
             : ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs)
{
235
	if (!vecArgs.empty()) //incorrect arguments, throw an exception
236
237
238
239
240
		throw InvalidArgumentException("Wrong number of arguments for \""+i_parname+"::"+i_algoname+"\"", AT);
}

std::string NoResampling::toString() const
{
241
	ostringstream ss;
242
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo << "[ ]";
243
244
245
	return ss.str();
}

246
void NoResampling::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
247
                            const std::vector<MeteoData>& vecM, MeteoData& md)
248
{
249
250
251
252
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);

	if (position == ResamplingAlgorithms::exact_match) {
253
		const double value = vecM[index](paramindex);
254
255
256
257
258
		if (value != IOUtils::nodata) {
			md(paramindex) = value; //propagate value
		}
	}

259
260
261
	return;
}

262
263

/**********************************************************************************/
264
265
266
267
NearestNeighbour::NearestNeighbour(const std::string& i_algoname, const std::string& i_parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs)
                 : ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs), extrapolate(false)
{
	const size_t nr_args = vecArgs.size();
268
269
270
	if (nr_args==0) return;
	if (nr_args==1) {
		if (vecArgs[0]=="extrapolate")
271
272
273
274
275
			extrapolate=true;
		else {
			IOUtils::convertString(window_size, vecArgs[0]);
			window_size /= 86400.; //user uses seconds, internally julian day is used
		}
276
	} else if (nr_args==2) {
277
278
		IOUtils::convertString(window_size, vecArgs[0]);
		window_size /= 86400.; //user uses seconds, internally julian day is used
279
		if (vecArgs[1]=="extrapolate")
280
281
282
283
284
285
286
287
288
289
			extrapolate=true;
		else
			throw InvalidArgumentException("Invalid argument \""+vecArgs[1]+"\" for \""+i_parname+"::"+i_algoname+"\"", AT);
	} else {
		throw InvalidArgumentException("Wrong number of arguments for \""+i_parname+"::"+i_algoname+"\"", AT);
	}
}

std::string NearestNeighbour::toString() const
{
290
	ostringstream ss;
291
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
292
	ss << "[ window_size=" << window_size << " extrapolate=" << boolalpha << extrapolate << noboolalpha << " ]";
293
294
295
	return ss.str();
}

296
void NearestNeighbour::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
297
                                const std::vector<MeteoData>& vecM, MeteoData& md)
298
{
299
300
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);
301

302
	if (position == ResamplingAlgorithms::exact_match) {
303
		const double value = vecM[index](paramindex);
304
		if (value != IOUtils::nodata) {
305
306
307
308
			md(paramindex) = value; //propagate value
			return;
		}
	}
309

310
311
	//if we are at the very beginning or end of vecM and !extrapolate, then there's nothing to do
	if (((!extrapolate) && (position == ResamplingAlgorithms::end))
312
	    || ((!extrapolate) && (position == ResamplingAlgorithms::begin)))
313
		return;
314

315
	const Date resampling_date = md.date;
316
	size_t indexP1=IOUtils::npos, indexP2=IOUtils::npos;
317
	getNearestValidPts(index, paramindex, vecM, resampling_date, window_size, indexP1, indexP2);
318
	const bool foundP1=(indexP1!=IOUtils::npos), foundP2=(indexP2!=IOUtils::npos);
319

320
	//Try to find the nearest neighbour, if there are two equally distant, then return the arithmetic mean
321
	if (foundP1 && foundP2) { //standard behavior
322
323
		const Duration diff1 = resampling_date - vecM[indexP1].date; //calculate time interval to element at index
		const Duration diff2 = vecM[indexP2].date - resampling_date; //calculate time interval to element at index
324
325
		const double val1 = vecM[indexP1](paramindex);
		const double val2 = vecM[indexP2](paramindex);
326

327
		if (IOUtils::checkEpsilonEquality(diff1.getJulian(true), diff2.getJulian(true), 0.1/1440.)) { //within 6 seconds
328
			md(paramindex) = Interpol1D::weightedMean(val1, val2, 0.5);
329
		} else if (diff1 < diff2) {
330
			md(paramindex) = val1;
331
		} else if (diff1 > diff2) {
332
			md(paramindex) = val2;
333
		}
334
	} else if (extrapolate) {
335
		if (foundP1 && !foundP2) { //nearest neighbour on found after index 'index'
336
			md(paramindex) = vecM[indexP1](paramindex);
337
		} else if (!foundP1 && foundP2) { //nearest neighbour on found before index 'index'
338
			md(paramindex) = vecM[indexP2](paramindex);
339
340
		} else { // no nearest neighbour with a value different from IOUtils::nodata
			return;
341
342
343
344
		}
	}
}

345
346

/**********************************************************************************/
347
348
349
350
LinearResampling::LinearResampling(const std::string& i_algoname, const std::string& i_parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs)
                 : ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs), extrapolate(false)
{
	const size_t nr_args = vecArgs.size();
351
352
353
	if (nr_args==0) return;
	if (nr_args==1) {
		if (vecArgs[0]=="extrapolate")
354
355
356
357
358
			extrapolate=true;
		else {
			IOUtils::convertString(window_size, vecArgs[0]);
			window_size /= 86400.; //user uses seconds, internally julian day is used
		}
359
	} else if (nr_args==2) {
360
361
		IOUtils::convertString(window_size, vecArgs[0]);
		window_size /= 86400.; //user uses seconds, internally julian day is used
362
		if (vecArgs[1]=="extrapolate")
363
364
365
366
367
368
369
370
371
372
			extrapolate=true;
		else
			throw InvalidArgumentException("Invalid argument \""+vecArgs[1]+"\" for \""+i_parname+"::"+i_algoname+"\"", AT);
	} else {
		throw InvalidArgumentException("Wrong number of arguments for \""+i_parname+"::"+i_algoname+"\"", AT);
	}
}

std::string LinearResampling::toString() const
{
373
	ostringstream ss;
374
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
375
	ss << "[ window_size=" << window_size << " extrapolate=" << boolalpha << extrapolate << noboolalpha << " ]";
376
377
378
	return ss.str();
}

379
void LinearResampling::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
380
                                const std::vector<MeteoData>& vecM, MeteoData& md)
381
{
382
383
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);
384

385
	if (position == ResamplingAlgorithms::exact_match) {
386
		const double value = vecM[index](paramindex);
387
		if (value != IOUtils::nodata) {
388
389
390
391
			md(paramindex) = value; //propagate value
			return;
		}
	}
392

393
394
	//if we are at the very beginning or end of vecM and !extrapolate, then there's nothing to do
	if (((!extrapolate) && (position == ResamplingAlgorithms::end))
395
	    || ((!extrapolate) && (position == ResamplingAlgorithms::begin)))
396
		return;
397

398
	const Date resampling_date = md.date;
399
	size_t indexP1=IOUtils::npos, indexP2=IOUtils::npos;
400
	getNearestValidPts(index, paramindex, vecM, resampling_date, window_size, indexP1, indexP2);
401
	bool foundP1=(indexP1!=IOUtils::npos), foundP2=(indexP2!=IOUtils::npos);
402
403

	//do nothing if we can't interpolate, and extrapolation is not explicitly activated
404
	if ((!extrapolate) && ((!foundP1) || (!foundP2)))
405
406
407
408
409
410
		return;
	//do nothing if not at least one value different from IOUtils::nodata has been found
	if (!foundP1 && !foundP2)
		return;

	//At this point we either have a valid indexP1 or indexP2 and we can at least try to extrapolate
411
412
413
	if (!foundP1 && foundP2) { //only nodata values found before index, try looking after indexP2
		for (size_t ii=indexP2+1; ii<vecM.size(); ii++) {
			if (vecM[ii](paramindex) != IOUtils::nodata) {
414
415
416
417
				indexP1 = ii;
				foundP1 = true;
				break;
			}
418
		}
419
420
	} else if (foundP1 && !foundP2) { //only nodata found after index, try looking before indexP1
		for (size_t ii=indexP1; (ii--) > 0; ){ 
421
			if (vecM[ii](paramindex) != IOUtils::nodata){
422
423
424
425
				indexP2=ii;
				foundP2 = true;
				break;
			}
426
		}
427
428
429
430
431
432
	}

	if (!foundP1 || !foundP2) //now at least two points need to be present
		return;

	//At this point indexP1 and indexP2 point to values that are different from IOUtils::nodata
433
	const double val1 = vecM[indexP1](paramindex);
434
	const double jul1 = vecM[indexP1].date.getJulian(true);
435
	const double val2 = vecM[indexP2](paramindex);
436
	const double jul2 = vecM[indexP2].date.getJulian(true);
437

438
	md(paramindex) = linearInterpolation(jul1, val1, jul2, val2, resampling_date.getJulian(true));
439
440
}

441

442
/**********************************************************************************/
443
Accumulate::Accumulate(const std::string& i_algoname, const std::string& i_parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs)
444
           : ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs),
445
             accumulate_period(IOUtils::nodata), strict(false)
446
{
447
	const size_t nr_args = vecArgs.size();
448
	if (nr_args<1 || nr_args>2)
449
		throw InvalidArgumentException("Please at least provide accumulation period (in seconds) for \""+i_parname+"::"+i_algoname+"\"", AT);
450

451
	bool period_read = false;
452
453
454
	for (size_t ii=0; ii<nr_args; ii++) {
		if (IOUtils::isNumeric(vecArgs[ii])) {
			if (period_read==true)
455
				throw InvalidArgumentException("Two arguments "+i_algoname+" resampling has been deprecated! Please use the \"PSUM_Distribute\" Processing Element instead!", AT);
456

457
458
			IOUtils::convertString(accumulate_period, vecArgs[ii]);
			accumulate_period /= 86400.; //user uses seconds, internally julian day is used
459
			if (accumulate_period<=0.) {
460
461
462
463
				std::ostringstream ss;
				ss << "Invalid accumulation period (" << accumulate_period << ") for \"" << i_parname << "::" << i_algoname << "\"";
				throw InvalidArgumentException(ss.str(), AT);
			}
464
			period_read = true;
465
		} else if (vecArgs[ii]=="strict" && !strict) {
466
			if (strict) //do not set strict more than once!
467
468
				throw InvalidArgumentException("Do not provide \"strict\" more than once for \""+i_parname+"::"+i_algoname+"\"", AT);
			strict = true;
469
470
		} else 
			throw InvalidArgumentException("Invalid argument \""+vecArgs[ii]+"\" for \""+i_parname+"::"+i_algoname+"\"", AT);
471
	}
472
473
}

474
475
std::string Accumulate::toString() const
{
476
	ostringstream ss;
477
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
478
	ss << "[ period=" << accumulate_period << " strict=" << boolalpha << strict << noboolalpha << " ]";
479
480
481
	return ss.str();
}

482
483
//find the index just before the start of accumulation period
size_t Accumulate::findStartOfPeriod(const std::vector<MeteoData>& vecM, const size_t& index, const Date& dateStart)
484
{
485
	size_t start_idx = IOUtils::npos;
486
	for (size_t idx=index; idx--> 0; ) {
487
		const Date curr_date = vecM[idx].date;
488
		if (curr_date <= dateStart) {
489
			start_idx = idx;
490
491
492
			break;
		}
	}
493

494
495
	return start_idx;
}
496

497
double Accumulate::easySampling(const std::vector<MeteoData>& vecM, const size_t& paramindex, const size_t& /*index*/, const size_t& start_idx, const Date& dateStart, const Date& resampling_date) const
498
{//to keep in mind: start_idx is last index <= dateStart and index is first index >= resampling_date
499
500
501
	double sum = IOUtils::nodata;
	const double start_val = partialAccumulateAtLeft(vecM, paramindex, start_idx, dateStart);
	const double end_val = partialAccumulateAtLeft(vecM, paramindex, start_idx, resampling_date);
502

503
	if (start_val!=IOUtils::nodata && end_val!=IOUtils::nodata)
504
		sum = end_val - start_val;
505

506
507
508
	return sum;
}

509
double Accumulate::complexSampling(const std::vector<MeteoData>& vecM, const size_t& paramindex, const size_t& index, const size_t& start_idx, const Date& dateStart, const Date& resampling_date) const
510
{//to keep in mind: start_idx is last index <= dateStart and index is first index >= resampling_date
511
	double sum = IOUtils::nodata;
512
	//resample begining point, in the [start_idx ; start_idx+1] interval
513
	const double start_value = partialAccumulateAtRight(vecM, paramindex, start_idx, dateStart);
514
	if (start_value!=IOUtils::nodata)
515
		sum=start_value;
516
	else if (strict) return IOUtils::nodata;
517

518
	//sum all whole periods AFTER the begining point, in the [start_idx+2 ; index-1] interval
519
	for (size_t idx=(start_idx+2); idx<index; idx++) {
520
		const double curr_value = vecM[idx](paramindex);
521
522
		if (curr_value!=IOUtils::nodata) {
			if (sum!=IOUtils::nodata) sum += curr_value;
523
			else sum = curr_value;
524
		} else if (strict) return IOUtils::nodata;
525
	}
526

527
	//resample end point, in the [index-1 ; index] interval
528
	const double end_val = partialAccumulateAtLeft(vecM, paramindex, index-1, resampling_date);
529
530
	if (end_val!=IOUtils::nodata) {
		if (sum!=IOUtils::nodata) sum += end_val;
531
		else sum = end_val;
532
	} else if (strict) return IOUtils::nodata;
533

534
535
536
	return sum;
}

537
538
539
//index is the first element AFTER the resampling_date
void Accumulate::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
                          const std::vector<MeteoData>& vecM, MeteoData& md)
540
{
541
542
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);
543
	if (position==ResamplingAlgorithms::begin || position==ResamplingAlgorithms::end)
544
		return;
545

546
547
	md(paramindex) = IOUtils::nodata;
	const Date resampling_date = md.date;
548

549
	const Date dateStart(resampling_date.getJulian() - accumulate_period, resampling_date.getTimeZone());
550
	const size_t start_idx = findStartOfPeriod(vecM, index, dateStart);
551
552
553
554
555
556
	if (start_idx==IOUtils::npos) {//No acceptable starting point found
		cerr << "[W] Could not accumulate " << vecM.at(0).getNameForParameter(paramindex) << ": ";
		cerr << "not enough data for accumulation period at date " << resampling_date.toString(Date::ISO) << "\n";
		return;
	}

557
	if ((index - start_idx) <= 1) {//easy upsampling when start & stop are in the same input time step
558
		//upsampling (for example, generate 15min values from hourly data)
559
		const double sum = easySampling(vecM, paramindex, index, start_idx, dateStart, resampling_date);
560
		md(paramindex) = sum; //if resampling was unsuccesful, sum==IOUtils::nodata
561
	} else {
562
		//downsampling (for example, generate daily values from hourly data)
563
564
		//and upsampling when resampled period falls accross a measurement timestamp
		const double sum = complexSampling(vecM, paramindex, index, start_idx, dateStart, resampling_date);
565
566
		md(paramindex) = sum; //if resampling was unsuccesful, sum==IOUtils::nodata
	}
567
568
}

569

570
571
572
573
574
575
/**********************************************************************************/
const double Solar::soil_albedo = .23; //grass
const double Solar::snow_albedo = .85; //snow
const double Solar::snow_thresh = .1; //if snow height greater than this threshold -> snow albedo

Solar::Solar(const std::string& i_algoname, const std::string& i_parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs)
576
            : ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs), cache_losses(), extrapolate(false)
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
{
	const size_t nr_args = vecArgs.size();
	if (nr_args==0) return;
	if (nr_args==1) {
		if (vecArgs[0]=="extrapolate")
			extrapolate=true;
		else {
			IOUtils::convertString(window_size, vecArgs[0]);
			window_size /= 86400.; //user uses seconds, internally julian day is used
		}
	} else if (nr_args==2) {
		IOUtils::convertString(window_size, vecArgs[0]);
		window_size /= 86400.; //user uses seconds, internally julian day is used
		if (vecArgs[1]=="extrapolate")
			extrapolate=true;
		else
			throw InvalidArgumentException("Invalid argument \""+vecArgs[1]+"\" for \""+i_parname+"::"+i_algoname+"\"", AT);
	} else {
		throw InvalidArgumentException("Wrong number of arguments for \""+i_parname+"::"+i_algoname+"\"", AT);
	}
}

std::string Solar::toString() const
{
	ostringstream ss;
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo << "[ ]";
	return ss.str();
}

606
double Solar::getPotentialH(const MeteoData& md)
607
608
609
610
611
612
613
614
{
	const double lat = md.meta.position.getLat();
	const double lon = md.meta.position.getLon();
	const double alt = md.meta.position.getAltitude();
	if (lat==IOUtils::nodata || lon==IOUtils::nodata || alt==IOUtils::nodata) return IOUtils::nodata;
	SunObject sun(lat, lon, alt);
	
	const double P = (md(MeteoData::P)!=IOUtils::nodata)? md(MeteoData::P) : Atmosphere::stdAirPressure(alt);
615
	const double HS = md(MeteoData::HS);
616
617
618
619
620
621
622
623
624
625
	double albedo = 0.5;
	if (HS!=IOUtils::nodata) //no big deal if we can not adapt the albedo
		albedo = (HS>=snow_thresh)? snow_albedo : soil_albedo;
	//if we don't have TA and RH, set them so the reduced precipitable water will get an average value
	double TA = md(MeteoData::TA);
	double RH = md(MeteoData::RH);
	if (TA==IOUtils::nodata || RH==IOUtils::nodata) {
		TA = 274.98;
		RH = 0.666;
	}
626
	
627
628
629
630
631
632
633
634
635
	sun.setDate(md.date.getJulian(), md.date.getTimeZone());
	sun.calculateRadiation(TA, RH, P, albedo);
	double toa, direct, diffuse;
	sun.getHorizontalRadiation(toa, direct, diffuse);
	const double global_h = direct+diffuse;
	
	return global_h;
}

636
bool Solar::computeLossFactor(const size_t& index, const size_t& paramindex,
637
                           const std::vector<MeteoData>& vecM, const Date& resampling_date, Points &pts)
638
639
640
641
642
643
644
{
	size_t indexP1=IOUtils::npos, indexP2=IOUtils::npos;
	getNearestValidPts(index, paramindex, vecM, resampling_date, window_size, indexP1, indexP2);
	bool foundP1=(indexP1!=IOUtils::npos), foundP2=(indexP2!=IOUtils::npos);
	
	if (!extrapolate && (!foundP1 || !foundP2)) return false;
	
645
646
647
648
	const double jul1 = vecM[indexP1].date.getJulian();
	if (jul1==pts.jul2) { //if the new start point is the old end point
		pts.jul1 = pts.jul2;
		pts.loss1 = pts.loss2;
649
650
	} else {
		const double pot1 = (foundP1)? getPotentialH( vecM[indexP1] ) : IOUtils::nodata;
651
652
		pts.loss1 = (pot1!=IOUtils::nodata && pot1>0.)? vecM[indexP1](paramindex) / pot1 : IOUtils::nodata;
		pts.jul1 = jul1;
653
654
655
	}
	
	const double pot2 = (foundP2)? getPotentialH( vecM[indexP2] ) : IOUtils::nodata;
656
657
	pts.loss2 = (pot2!=IOUtils::nodata && pot2>0.)? vecM[indexP2](paramindex) / pot2 : IOUtils::nodata;
	pts.jul2 = vecM[indexP2].date.getJulian();
658
659
660
	return true;
}

661
double Solar::interpolateLossFactor(const double& resampling_jul, const Points &pts)
662
{
663
664
665
666
667
668
669
	if (pts.loss1!=IOUtils::nodata && pts.loss2!=IOUtils::nodata) {
		const double weight = (resampling_jul - pts.jul1) / (pts.jul2 - pts.jul1);
		return Interpol1D::weightedMean(pts.loss1, pts.loss2, weight);
	} else if (pts.loss1==IOUtils::nodata && pts.loss2!=IOUtils::nodata){
		return pts.loss2;
	} else if (pts.loss2==IOUtils::nodata && pts.loss1!=IOUtils::nodata) {
		return pts.loss1;
670
671
672
673
674
	}
	
	return 1.;
}

675
676
677
678
679
680
681
682
683
684
685
686
void Solar::resample(const size_t& index, const ResamplingPosition& /*position*/, const size_t& paramindex,
                           const std::vector<MeteoData>& vecM, MeteoData& md)
{
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);

	if (paramindex!=MeteoData::ISWR && paramindex!=MeteoData::RSWR)
		throw IOException("This method only applies to short wave radiation! (either ISWR or RSWR)", AT);

	const double pot_pt = getPotentialH( md );
	if (pot_pt==IOUtils::nodata) return;
	
687
688
689
690
	const double resampling_jul = md.date.getJulian();
	Points pts( cache_losses[ vecM[0].meta.getHash() ] );
	if (pts.jul1==0. || (resampling_jul<pts.jul1 || resampling_jul>pts.jul2)) {
		const bool status = computeLossFactor(index, paramindex, vecM, md.date, pts);
691
		if (!status) return;
692
		cache_losses[ vecM[0].meta.getHash() ] = pts;
693
694
	}
	
695
	const double loss = interpolateLossFactor(resampling_jul, pts);
696
697
698
699
700
	md(paramindex) = loss * pot_pt;
	md.setResampled(true);
}


701
/**********************************************************************************/
702
703
704
const double Daily_solar::soil_albedo = .23; //grass
const double Daily_solar::snow_albedo = .85; //snow
const double Daily_solar::snow_thresh = .1; //if snow height greater than this threshold -> snow albedo
705
706
const size_t Daily_solar::samples_per_day = 24*3; //every 20 minutes

707
Daily_solar::Daily_solar(const std::string& i_algoname, const std::string& i_parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs)
708
            : ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs), radiation(), station_index(), dateStart(), dateEnd(), loss_factor()
709
710
{
	const size_t nr_args = vecArgs.size();
711
	if (nr_args>0) {
712
713
714
715
716
717
		throw InvalidArgumentException("Too many arguments for \""+i_parname+"::"+i_algoname+"\"", AT);
	}
}

std::string Daily_solar::toString() const
{
718
	ostringstream ss;
719
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo << "[ ]";
720
721
722
	return ss.str();
}

723
724
725
726
727
728
729
730
731
732
733
734
//compute the daily sum of solar radiation as well as fill a vector containing the solar radiation at a regular sampling rate
double Daily_solar::compRadiation(const double& lat, const double& lon, const double& alt, const double& HS, const size_t& stat_idx)
{
	const double P = Atmosphere::stdAirPressure(alt);
	double albedo = 0.5;
	if (HS!=IOUtils::nodata) //no big deal if we can not adapt the albedo
		albedo = (HS>=snow_thresh)? snow_albedo : soil_albedo;
	const double TA=274.98, RH=0.666; //the reduced precipitable water will get an average value

	SunObject sun(lat, lon, alt);
	double sum = 0.;
	size_t index=0;
735
	for (Date date(dateStart[stat_idx]); date<dateEnd[stat_idx]; date += 1./double(samples_per_day)) {
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
		//compute potential solar radiation at this time step
		sun.setDate(date.getJulian(), date.getTimeZone());
		sun.calculateRadiation(TA, RH, P, albedo);
		double toa, direct, diffuse;
		sun.getHorizontalRadiation(toa, direct, diffuse);
		const double global_h = direct+diffuse;

		//compute the integral by a simple triangle method
		sum += global_h * static_cast<double>(24*3600/samples_per_day);

		//store the radiation for later reuse
		radiation[stat_idx][index++] = global_h;
	}

	return sum;
}

//interpolate the solar radiation from the vector containing regularly sampled solar radiation
double Daily_solar::getSolarInterpol(const Date& resampling_date, const size_t& stat_idx) const
{
	const double in_day = (resampling_date.getJulian() - dateStart[stat_idx].getJulian()) / (dateEnd[stat_idx].getJulian() - dateStart[stat_idx].getJulian());
	const size_t vec_index = static_cast<size_t>( Optim::floor(in_day*samples_per_day) ); //so the sample will be between vec_index and vec_index+1
	const double weight = in_day*static_cast<double>(samples_per_day) - static_cast<double>(vec_index);

	return Interpol1D::weightedMean(radiation[stat_idx][vec_index], radiation[stat_idx][vec_index+1], weight);
}

//a new, previously unknown station has been found, allocate the memory
size_t Daily_solar::getStationIndex(const std::string& key)
{
	const size_t nr_stations = station_index.size();
767
768
	for (size_t ii=0; ii<nr_stations; ++ii) {
		if (station_index[ii]==key)
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
			return ii;
	}

	radiation.push_back( vector<double>(samples_per_day, 0.) );
	loss_factor.push_back( 0. );

	Date null_date(0., 0.);
	dateStart.push_back( null_date );
	dateEnd.push_back( null_date );

	station_index.push_back( key );
	return nr_stations; //the new index is the old nr_stations
}

void Daily_solar::resample(const size_t& index, const ResamplingPosition& /*position*/, const size_t& paramindex,
                           const std::vector<MeteoData>& vecM, MeteoData& md)
785
786
787
788
{
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);

789
	if (paramindex!=MeteoData::ISWR && paramindex!=MeteoData::RSWR)
790
791
792
793
794
		throw IOException("This method only applies to short wave radiation! (either ISWR or RSWR)", AT);

	const double lat = md.meta.position.getLat();
	const double lon = md.meta.position.getLon();
	const double alt = md.meta.position.getAltitude();
795
	if (lat==IOUtils::nodata || lon==IOUtils::nodata || alt==IOUtils::nodata) return;
796
	const double HS = md(MeteoData::HS);
797

798
799
	//get station index
	const size_t stat_idx = getStationIndex(md.meta.stationID);
800

801
	//has the radiation already been calculated for this day and station?
802
	if (md.date<dateStart[stat_idx] || md.date>=dateEnd[stat_idx]) {
803
804
805
806
		dateStart[stat_idx] = getDailyStart(md.date);
		dateEnd[stat_idx] = dateStart[stat_idx]+1.;
		//const size_t indexP = getNearestValidPt(vecM, paramindex, stat_idx, index);
		const size_t indexP = getDailyValue(vecM, paramindex, index, dateStart[stat_idx], dateEnd[stat_idx]);
807
		if (indexP==IOUtils::npos) { //no daily sum found for the current day
808
809
810
811
812
813
			loss_factor[stat_idx] = IOUtils::nodata;
			return;
		}

		const double daily_sum = compRadiation(lat, lon, alt, HS, stat_idx);
		loss_factor[stat_idx] = (daily_sum>0)? vecM[indexP](paramindex) / daily_sum : 0.; //in case of polar night...
814
815
	}

816
	if (loss_factor[stat_idx]==IOUtils::nodata) //the station could not be calculated for this day
817
818
819
820
		return;

	//interpolate radiation for this timestep and write it out
	const double rad = getSolarInterpol(md.date, stat_idx);
821
	if (paramindex==MeteoData::ISWR) {
822
823
824
825
826
827
828
		md(paramindex) = loss_factor[stat_idx] * rad;
	} else {
		double albedo = 0.5;
		if (HS!=IOUtils::nodata) //no big deal if we can not adapt the albedo
			albedo = (HS>=snow_thresh)? snow_albedo : soil_albedo;
		md(paramindex) = loss_factor[stat_idx] * rad * albedo;
	}
829

830
	md.setResampled(true);
831
832
}

833
834
835
836
837
838
839
840
841
842
843
844

/**********************************************************************************/
DailyAverage::DailyAverage(const std::string& i_algoname, const std::string& i_parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs)
                 : ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs), range(IOUtils::nodata), phase(0.)
{
	const size_t nr_args = vecArgs.size();
	
	if (nr_args==1) {
		IOUtils::convertString(range, vecArgs[0]);
	} else if (nr_args==2) {
		IOUtils::convertString(range, vecArgs[0]);
		IOUtils::convertString(phase, vecArgs[1]);
845
		phase *= -1; //shift the minimum *later* in the day
846
	} else if (nr_args>2) {
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
		throw InvalidArgumentException("Wrong number of arguments for \""+i_parname+"::"+i_algoname+"\"", AT);
	}
}

std::string DailyAverage::toString() const
{
	ostringstream ss;
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	ss << "[ window_size=" << window_size;
	if (range!=IOUtils::nodata) 
		ss << " range=" << range;
	if (phase!=0.)
		ss << " phase=" << phase;
	ss << " ]";
	return ss.str();
}

864
865
866
867
868
869
870
871
872
double DailyAverage::getValue(const std::vector<MeteoData>& vecM, const size_t& paramindex, const size_t& index, const Date& dayStart, const double& frac_day) const
{
	const Date dayEnd = dayStart+1.;
	const size_t indexP = getDailyValue(vecM, paramindex, index, dayStart, dayEnd);
	if (indexP==IOUtils::npos)
		return IOUtils::nodata;
	
	const double avg = vecM[indexP](paramindex); //from getDailyValue we are guaranteed to have a value
	
873
874
875
876
877
878
879
880
881
882
883
884
885
	//search for min or max or both
	const string param_name = vecM[indexP].getNameForParameter(paramindex);
	const string min_name = param_name+"_MIN";
	const string max_name = param_name+"_MAX";
	const double min = (vecM[indexP].param_exists( min_name ))?  vecM[indexP](min_name) : IOUtils::nodata;
	const double max = (vecM[indexP].param_exists( max_name ))?  vecM[indexP](max_name) : IOUtils::nodata;
	double A=range; //fallback: user provided range
	
	if (min==IOUtils::nodata && max==IOUtils::nodata) {
		if (A==IOUtils::nodata) {
			ostringstream ss;
			ss << "Missing parameters for " << algo << " on " << dayStart.toString(Date::ISO_DATE) << ". Please consider providing a fallback range argument!";
			throw InvalidArgumentException(ss.str(), AT);
886
		}
887
	} else if (min!=IOUtils::nodata && max==IOUtils::nodata) {
888
		A = (avg -  min);
889
	} else if (min==IOUtils::nodata && max!=IOUtils::nodata) {
890
		A =(max - avg);
891
892
	} else
		throw InvalidArgumentException("Providing both AVG, MIN and MAX for \'"+algo+"\' is not supported!", AT);
893
894
895
896

	return A * sin( 2.*Cst::PI * (frac_day-.25+phase) ) + avg;
}

897
898
899
900
901
902
void DailyAverage::resample(const size_t& index, const ResamplingPosition& /*position*/, const size_t& paramindex,
                                const std::vector<MeteoData>& vecM, MeteoData& md)
{
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);

903
	//define various time parameters
904
	const double julian = md.date.getJulian();
905
906
907
	double frac_day = (julian+.5 - Optim::intPart(julian+.5)); // julian day starts at noon and we want the fraction of the day, from 0 to 1!
	if (frac_day==0.) frac_day = 1.; //because we considere here that midnight belongs to the day before
	
908
909
910
911
	const Date dayStart = getDailyStart(md.date);
	const double val1 = getValue(vecM, paramindex, index, dayStart, frac_day);
	if (val1==IOUtils::nodata)
		return; //we could also decide to use val0 & val2 
912
	
913
914
915
	//get the other values (day before, day after)
	const double val0 = getValue(vecM, paramindex, index, dayStart-1., frac_day);
	const double val2 = getValue(vecM, paramindex, index, dayStart+1., frac_day);
916
	
917
918
919
920
921
	double val = val1;
	if (frac_day!=0.5) { //at noon, we purely take val1
		const double w1 = 1./Optim::pow2( abs(.5-frac_day) );
		double norm = w1;
		val *= w1;
922
		
923
924
925
926
		if (val0!=IOUtils::nodata) {
			const double w0 = 1./Optim::pow2( frac_day+0.5 );
			val += w0*val0;
			norm += w0;
927
		}
928
929
930
931
		if (val2!=IOUtils::nodata) {
			const double w2 = 1./Optim::pow2( 1.5-frac_day );
			val += w2*val2;
			norm += w2;
932
		}
933
		val /= norm;
934
935
	}
	
936
	md(paramindex) = val;
937
938
	//very basic range checking
	//HACK this means that this should be implemented as a data creator so it could be filtered afterward
939
	//HACK or we could one more pass of filtering *after* the resampling
940
	if (paramindex==MeteoData::RH) {
941
942