WSL/SLF GitLab Repository

ResamplingAlgorithms.cc 33 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
23
#include <meteoio/meteoLaws/Atmosphere.h>
#include <meteoio/meteoLaws/Sun.h>
#include <meteoio/meteoStats/libinterpol1D.h>
#include <meteoio/meteoFilters/ProcHNWDistribute.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
	} else if (algoname == "DAILY_SOLAR") {
45
		return new Daily_solar(algoname, parname, dflt_window_size, vecArgs);
46
47
	} else if (algoname == "DAILY_AVG") {
		return new DailyAverage(algoname, parname, dflt_window_size, vecArgs);
48
49
50
51
	} else {
		throw IOException("The resampling algorithm '"+algoname+"' is not implemented" , AT);
	}
}
52

53
54
55
//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)
56
{
57
	const size_t end = pos+1;
58
	if (end>=vecM.size()) return IOUtils::nodata; //reaching the end of the input vector
59
60
61
62
63
64

	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);
65

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

68
69
70
71
72
73
74
	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)
{
75
76
77
78
79
80
81
82
	const size_t end = pos+1;
	if(end>=vecM.size()) return IOUtils::nodata; //reaching the end of the input vector

	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);
83

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

86
	return valend - left_accumulation;
87
88
}

89
90
91
92
93
/**
 * @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
94
 * @param resampling_date date to resample
95
96
97
98
99
 * @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,
100
                                              const double& window_size, size_t& indexP1, size_t& indexP2) //HACK
101
{
102
103
	indexP1=IOUtils::npos;
	indexP2=IOUtils::npos;
104

105
	const Date dateStart = resampling_date - window_size;
106
	for (size_t ii=pos; ii-- >0; ) {
107
		if (vecM[ii].date < dateStart) break;
108
		if (vecM[ii](paramindex) != IOUtils::nodata) {
109
110
111
112
113
114
115
			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;
116
	for (size_t ii=pos; ii<vecM.size(); ++ii) {
117
118
119
120
121
122
123
124
125
126
127
128
129
130
		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
131
 * @param x x-coordinate of desired point
132
133
134
 * @return y-coordinate of desired point
 */
double ResamplingAlgorithms::linearInterpolation(const double& x1, const double& y1,
135
                                       const double& x2, const double& y2, const double& x)
136
137
138
139
140
141
142
143
{
	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;

144
	return (a*x + b);
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
/**
 * @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
 */
173
size_t ResamplingAlgorithms::getDailyValue(const std::vector<MeteoData>& vecM, const size_t& paramindex, size_t pos, const Date& intervalStart, const Date& intervalEnd)
174
175
176
{
	size_t indexP1=IOUtils::npos;
	size_t indexP2=IOUtils::npos;
177
178
179
	
	//if pos was not properly pre-positioned, do it
	if (vecM[pos].date<intervalStart) {
180
		for (; pos<vecM.size(); pos++)
181
182
183
			if (vecM[pos].date>=intervalStart) break;
	}
	if (vecM[pos].date>intervalEnd) {
184
185
		for (size_t ii = pos; ii-- >0; ) {
			pos=ii; //because ii gets corrupted at the final iteration if going all the way down
186
			if (vecM[pos].date<=intervalEnd) break;
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

	//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;
}

226
/**********************************************************************************
227
 * The following functions are implementations of different resampling algorithms *
228
 **********************************************************************************/
229

230
231
232
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)
{
233
	if (!vecArgs.empty()) //incorrect arguments, throw an exception
234
235
236
237
238
		throw InvalidArgumentException("Wrong number of arguments for \""+i_parname+"::"+i_algoname+"\"", AT);
}

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

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

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

258
259
260
	return;
}

261
262

/**********************************************************************************/
263
264
265
266
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();
267
268
269
	if (nr_args==0) return;
	if (nr_args==1) {
		if (vecArgs[0]=="extrapolate")
270
271
272
273
274
			extrapolate=true;
		else {
			IOUtils::convertString(window_size, vecArgs[0]);
			window_size /= 86400.; //user uses seconds, internally julian day is used
		}
275
	} else if (nr_args==2) {
276
277
		IOUtils::convertString(window_size, vecArgs[0]);
		window_size /= 86400.; //user uses seconds, internally julian day is used
278
		if (vecArgs[1]=="extrapolate")
279
280
281
282
283
284
285
286
287
288
			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
{
289
	ostringstream ss;
290
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
291
	ss << "[ window_size=" << window_size << " extrapolate=" << boolalpha << extrapolate << noboolalpha << " ]";
292
293
294
	return ss.str();
}

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

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

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

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

319
	//Try to find the nearest neighbour, if there are two equally distant, then return the arithmetic mean
320
	if (foundP1 && foundP2) { //standard behavior
321
322
		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
323
324
		const double val1 = vecM[indexP1](paramindex);
		const double val2 = vecM[indexP2](paramindex);
325

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

344
345

/**********************************************************************************/
346
347
348
349
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();
350
351
352
	if (nr_args==0) return;
	if (nr_args==1) {
		if (vecArgs[0]=="extrapolate")
353
354
355
356
357
			extrapolate=true;
		else {
			IOUtils::convertString(window_size, vecArgs[0]);
			window_size /= 86400.; //user uses seconds, internally julian day is used
		}
358
	} else if (nr_args==2) {
359
360
		IOUtils::convertString(window_size, vecArgs[0]);
		window_size /= 86400.; //user uses seconds, internally julian day is used
361
		if (vecArgs[1]=="extrapolate")
362
363
364
365
366
367
368
369
370
371
			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
{
372
	ostringstream ss;
373
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
374
	ss << "[ window_size=" << window_size << " extrapolate=" << boolalpha << extrapolate << noboolalpha << " ]";
375
376
377
	return ss.str();
}

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

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

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

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

	//do nothing if we can't interpolate, and extrapolation is not explicitly activated
403
	if ((!extrapolate) && ((!foundP1) || (!foundP2)))
404
405
406
407
408
409
		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
410
411
412
	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) {
413
414
415
416
				indexP1 = ii;
				foundP1 = true;
				break;
			}
417
		}
418
419
	} else if (foundP1 && !foundP2) { //only nodata found after index, try looking before indexP1
		for (size_t ii=indexP1; (ii--) > 0; ){ 
420
			if (vecM[ii](paramindex) != IOUtils::nodata){
421
422
423
424
				indexP2=ii;
				foundP2 = true;
				break;
			}
425
		}
426
427
428
429
430
431
	}

	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
432
	const double val1 = vecM[indexP1](paramindex);
433
	const double jul1 = vecM[indexP1].date.getJulian(true);
434
	const double val2 = vecM[indexP2](paramindex);
435
	const double jul2 = vecM[indexP2].date.getJulian(true);
436

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

440

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

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

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

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

481
482
//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)
483
{
484
	size_t start_idx = IOUtils::npos;
485
	for (size_t idx=index; idx--> 0; ) {
486
		const Date curr_date = vecM[idx].date;
487
		if (curr_date <= dateStart) {
488
			start_idx = idx;
489
490
491
			break;
		}
	}
492

493
494
	return start_idx;
}
495

496
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
497
{//to keep in mind: start_idx is last index <= dateStart and index is first index >= resampling_date
498
499
500
	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);
501

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

505
506
507
	return sum;
}

508
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
509
{//to keep in mind: start_idx is last index <= dateStart and index is first index >= resampling_date
510
	double sum = IOUtils::nodata;
511
	//resample begining point, in the [start_idx ; start_idx+1] interval
512
	const double start_value = partialAccumulateAtRight(vecM, paramindex, start_idx, dateStart);
513
	if (start_value!=IOUtils::nodata)
514
		sum=start_value;
515
	else if (strict) return IOUtils::nodata;
516

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

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

533
534
535
	return sum;
}

536
537
538
//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)
539
{
540
541
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);
542
	if (position==ResamplingAlgorithms::begin || position==ResamplingAlgorithms::end)
543
		return;
544

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

548
	const Date dateStart(resampling_date.getJulian() - accumulate_period, resampling_date.getTimeZone());
549
	const size_t start_idx = findStartOfPeriod(vecM, index, dateStart);
550
551
552
553
554
555
	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;
	}

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

568

569
/**********************************************************************************/
570
571
572
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
573
574
const size_t Daily_solar::samples_per_day = 24*3; //every 20 minutes

575
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)
576
            : ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs), radiation(), station_index(), dateStart(), dateEnd(), loss_factor()
577
578
{
	const size_t nr_args = vecArgs.size();
579
	if (nr_args>0) {
580
581
582
583
584
585
		throw InvalidArgumentException("Too many arguments for \""+i_parname+"::"+i_algoname+"\"", AT);
	}
}

std::string Daily_solar::toString() const
{
586
	ostringstream ss;
587
588
589
590
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	return ss.str();
}

591
592
593
594
595
596
597
598
599
600
601
602
//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;
603
	for (Date date(dateStart[stat_idx]); date<dateEnd[stat_idx]; date += 1./double(samples_per_day)) {
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
		//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();
635
636
	for (size_t ii=0; ii<nr_stations; ++ii) {
		if (station_index[ii]==key)
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
			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)
653
654
655
656
{
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);

657
	if (paramindex!=MeteoData::ISWR && paramindex!=MeteoData::RSWR)
658
659
660
661
662
		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();
663
	if (lat==IOUtils::nodata || lon==IOUtils::nodata || alt==IOUtils::nodata) return;
664
	const double HS = md(MeteoData::HS);
665

666
667
	//get station index
	const size_t stat_idx = getStationIndex(md.meta.stationID);
668

669
	//has the radiation already been calculated for this day and station?
670
	if (md.date<dateStart[stat_idx] || md.date>=dateEnd[stat_idx]) {
671
672
673
674
		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]);
675
		if (indexP==IOUtils::npos) { //no daily sum found for the current day
676
677
678
679
680
681
			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...
682
683
	}

684
	if (loss_factor[stat_idx]==IOUtils::nodata) //the station could not be calculated for this day
685
686
687
688
		return;

	//interpolate radiation for this timestep and write it out
	const double rad = getSolarInterpol(md.date, stat_idx);
689
	if (paramindex==MeteoData::ISWR) {
690
691
692
693
694
695
696
		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;
	}
697

698
	md.setResampled(true);
699
700
}

701
702
703
704
705
706
707
708
709
710
711
712

/**********************************************************************************/
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]);
713
		phase *= -1; //shift the minimum *later* in the day
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
	} else {
		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();
}

732
733
734
735
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
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
	
	double A=range;
	if (A==IOUtils::nodata) { //ie no "range" was provided by the user
		//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;
		
		if (!(min!=IOUtils::nodata) != !(max!=IOUtils::nodata)) { //cheap implementation of XOR
			if  (min!=IOUtils::nodata)
				A = (avg -  min) * 2.;
			else //max!=IOUtils::nodata
				A =(max - avg) * 2.;
		} else {
			if (min==IOUtils::nodata && max==IOUtils::nodata)
				return IOUtils::nodata;
			else
				throw InvalidArgumentException("Providing both AVG, MIN and MAX for \'"+algo+"\' is not supported!", AT);
		}
	}

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

766
767
768
769
770
771
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);

772
	//define various time parameters
773
	const double julian = md.date.getJulian();
774
775
776
	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
	
777
778
779
780
	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 
781
	
782
783
784
	//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);
785
	
786
787
788
789
790
	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;
791
		
792
793
794
795
		if (val0!=IOUtils::nodata) {
			const double w0 = 1./Optim::pow2( frac_day+0.5 );
			val += w0*val0;
			norm += w0;
796
		}
797
798
799
800
		if (val2!=IOUtils::nodata) {
			const double w2 = 1./Optim::pow2( 1.5-frac_day );
			val += w2*val2;
			norm += w2;
801
		}
802
		val /= norm;
803
804
	}
	
805
	md(paramindex) = val;
806
807
	//very basic range checking
	//HACK this means that this should be implemented as a data creator so it could be filtered afterward
808
	//HACK or we could one more pass of filtering *after* the resampling
809
	if (paramindex==MeteoData::RH) {
810
811
		if (md(paramindex)<0.01) md(paramindex)=0.01;
		if (md(paramindex)>1.) md(paramindex)=1.;
812
813
814
815
816
817
818
819
	} else if (paramindex==MeteoData::TA ||  paramindex==MeteoData::VW ||  paramindex==MeteoData::VW_MAX || paramindex==MeteoData::ISWR || paramindex==MeteoData::RSWR || paramindex==MeteoData::ILWR || paramindex==MeteoData::TSG || paramindex==MeteoData::TSS || paramindex==MeteoData::TAU_CLD) {
		if (md(paramindex)<0.)
			md(paramindex)=0.;
	}
		
	md.setResampled(true);
}

820
821
} //namespace