WSL/SLF GitLab Repository

ResamplingAlgorithms.cc 30.8 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/GeneratorAlgorithms.h>
20
21
22
#include <meteoio/MathOptim.h>
#include <meteoio/meteolaws/Atmosphere.h>
#include <meteoio/meteolaws/Sun.h>
23
#include <meteoio/meteostats/libinterpol1D.h>
24
#include <cmath>
25
#include <algorithm>
26
27
28
29
30

using namespace std;

namespace mio {

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

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

50
51
52
//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)
53
{
54
55
56
57
58
59
60
61
	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);
62

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

65
66
67
68
69
70
71
	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)
{
72
73
74
75
76
77
78
79
	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);
80

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

83
	return valend - left_accumulation;
84
85
}

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

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

141
	return (a*x + b);
142
143
}

144
/**********************************************************************************
145
 * The following functions are implementations of different resampling algorithms *
146
 **********************************************************************************/
147

148
149
150
151
152
153
154
155
156
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)
{
	if(!vecArgs.empty()) //incorrect arguments, throw an exception
		throw InvalidArgumentException("Wrong number of arguments for \""+i_parname+"::"+i_algoname+"\"", AT);
}

std::string NoResampling::toString() const
{
157
	ostringstream ss;
158
159
160
161
162
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	ss << "[ ]";
	return ss.str();
}

163
void NoResampling::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
164
                            const std::vector<MeteoData>& vecM, MeteoData& md)
165
{
166
167
168
169
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);

	if (position == ResamplingAlgorithms::exact_match) {
170
		const double value = vecM[index](paramindex);
171
172
173
174
175
		if (value != IOUtils::nodata) {
			md(paramindex) = value; //propagate value
		}
	}

176
177
178
	return;
}

179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
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();
	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 NearestNeighbour::toString() const
{
205
	ostringstream ss;
206
207
208
209
210
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	ss << "[ window_size=" << window_size << " extrapolate=" << extrapolate << " ]";
	return ss.str();
}

211
void NearestNeighbour::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
212
                                const std::vector<MeteoData>& vecM, MeteoData& md)
213
{
214
215
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);
216

217
	if (position == ResamplingAlgorithms::exact_match) {
218
		const double value = vecM[index](paramindex);
219
		if (value != IOUtils::nodata) {
220
221
222
223
			md(paramindex) = value; //propagate value
			return;
		}
	}
224

225
226
	//if we are at the very beginning or end of vecM and !extrapolate, then there's nothing to do
	if (((!extrapolate) && (position == ResamplingAlgorithms::end))
227
	    || ((!extrapolate) && (position == ResamplingAlgorithms::begin)))
228
		return;
229

230
	const Date resampling_date = md.date;
231
	size_t indexP1=IOUtils::npos, indexP2=IOUtils::npos;
232
	getNearestValidPts(index, paramindex, vecM, resampling_date, window_size, indexP1, indexP2);
233
	const bool foundP1=(indexP1!=IOUtils::npos), foundP2=(indexP2!=IOUtils::npos);
234

235
	//Try to find the nearest neighbour, if there are two equally distant, then return the arithmetic mean
236
	if (foundP1 && foundP2) { //standard behavior
237
238
		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
239
240
		const double val1 = vecM[indexP1](paramindex);
		const double val2 = vecM[indexP2](paramindex);
241

242
		if (IOUtils::checkEpsilonEquality(diff1.getJulian(true), diff2.getJulian(true), 0.1/1440.)){ //within 6 seconds
243
			md(paramindex) = Interpol1D::weightedMean(val1, val2, 0.5);
244
		} else if (diff1 < diff2){
245
			md(paramindex) = val1;
246
		} else if (diff1 > diff2){
247
			md(paramindex) = val2;
248
		}
249
	} else if (extrapolate) {
250
251
252
253
		if(foundP1 && !foundP2){ //nearest neighbour on found after index 'index'
			md(paramindex) = vecM[indexP1](paramindex);
		} else if (!foundP1 && foundP2){ //nearest neighbour on found before index 'index'
			md(paramindex) = vecM[indexP2](paramindex);
254
255
		} else { // no nearest neighbour with a value different from IOUtils::nodata
			return;
256
257
258
259
		}
	}
}

260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
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();
	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 LinearResampling::toString() const
{
286
	ostringstream ss;
287
288
289
290
291
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	ss << "[ window_size=" << window_size << " extrapolate=" << extrapolate << " ]";
	return ss.str();
}

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

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

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

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

	//do nothing if we can't interpolate, and extrapolation is not explicitly activated
317
	if ((!extrapolate) && ((!foundP1) || (!foundP2)))
318
319
320
321
322
323
		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
324
	if (!foundP1 && foundP2){ //only nodata values found before index, try looking after indexP2
325
		for (size_t ii=indexP2+1; ii<vecM.size(); ii++){
326
			if (vecM[ii](paramindex) != IOUtils::nodata){
327
328
329
330
				indexP1 = ii;
				foundP1 = true;
				break;
			}
331
		}
332
	} else if (foundP1 && !foundP2){ //only nodata found after index, try looking before indexP1
333
		for (size_t ii=indexP1; (ii--) > 0; ){
334
			if (vecM[ii](paramindex) != IOUtils::nodata){
335
336
337
338
				indexP2=ii;
				foundP2 = true;
				break;
			}
339
		}
340
341
342
343
344
345
	}

	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
346
	const double val1 = vecM[indexP1](paramindex);
347
	const double jul1 = vecM[indexP1].date.getJulian(true);
348
	const double val2 = vecM[indexP2](paramindex);
349
	const double jul2 = vecM[indexP2].date.getJulian(true);
350

351
	md(paramindex) = linearInterpolation(jul1, val1, jul2, val2, resampling_date.getJulian(true));
352
353
}

354

355
Accumulate::Accumulate(const std::string& i_algoname, const std::string& i_parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs)
356
357
           : ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs),
           vecCache(), station_index(),
358
           measured_period(IOUtils::nodata), accumulate_period(IOUtils::nodata), strict(false)
359
{
360
	const size_t nr_args = vecArgs.size();
361
362
	if(nr_args<1 || nr_args>3)
		throw InvalidArgumentException("Please at least provide accumulation period (in seconds) for \""+i_parname+"::"+i_algoname+"\"", AT);
363

364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
	vector<double> numArgs;
	for(size_t ii=0; ii<nr_args; ii++) {
		if(IOUtils::isNumeric(vecArgs[ii])) {
			double tmp;
			IOUtils::convertString(tmp, vecArgs[ii]);
			tmp /= 86400.; //user uses seconds, internally julian day is used
			if(tmp<=0.) {
				std::ostringstream ss;
				ss << "Invalid accumulation period (" << accumulate_period << ") for \"" << i_parname << "::" << i_algoname << "\"";
				throw InvalidArgumentException(ss.str(), AT);
			}
			numArgs.push_back( tmp );
		} else if (vecArgs[ii]=="strict" && !strict) {
			if(strict) //do not set strict more than once!
				throw InvalidArgumentException("Do not provide \"strict\" more than once for \""+i_parname+"::"+i_algoname+"\"", AT);
			strict = true;
		} else throw InvalidArgumentException("Invalid argument \""+vecArgs[ii]+"\" for \""+i_parname+"::"+i_algoname+"\"", AT);
381
	}
382

383
384
385
386
387
388
389
	if(numArgs.size()==1) {
		accumulate_period = numArgs[0];
	} else if(numArgs.size()==2) {
		measured_period = numArgs[0];
		accumulate_period = numArgs[1];
	} else //this should never happen
		throw InvalidArgumentException("Too many arguments provided for \""+i_parname+"::"+i_algoname+"\"", AT);
390
391
}

392
393
std::string Accumulate::toString() const
{
394
	ostringstream ss;
395
396
397
398
399
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	ss << "[ accumulate_period=" << accumulate_period << " strict=" << strict << " ]";
	return ss.str();
}

400
401
//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)
402
{
403
	size_t start_idx = IOUtils::npos;
404
	for (size_t idx=index; idx--> 0; ) {
405
406
407
		const Date curr_date = vecM[idx].date;
		if(curr_date <= dateStart) {
			start_idx = idx;
408
409
410
			break;
		}
	}
411

412
413
	return start_idx;
}
414

415
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
416
{//to keep in mind: start_idx is last index <= dateStart and index is first index >= resampling_date
417
418
419
	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);
420
421

	if(start_val!=IOUtils::nodata && end_val!=IOUtils::nodata)
422
		sum = end_val - start_val;
423

424
425
426
	return sum;
}

427
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
428
{//to keep in mind: start_idx is last index <= dateStart and index is first index >= resampling_date
429
	double sum = IOUtils::nodata;
430
	//resample begining point, in the [start_idx ; start_idx+1] interval
431
432
433
434
	const double start_value = partialAccumulateAtRight(vecM, paramindex, start_idx, dateStart);
	if(start_value!=IOUtils::nodata)
		sum=start_value;
	else if(strict) return IOUtils::nodata;
435

436
437
	//sum all whole periods AFTER the begining point, in the [start_idx+2 ; index-1] interval
	for(size_t idx=(start_idx+2); idx<index; idx++) {
438
439
440
441
		const double curr_value = vecM[idx](paramindex);
		if(curr_value!=IOUtils::nodata) {
			if(sum!=IOUtils::nodata) sum += curr_value;
			else sum = curr_value;
442
		} else if(strict) return IOUtils::nodata;
443
	}
444

445
	//resample end point, in the [index-1 ; index] interval
446
	const double end_val = partialAccumulateAtLeft(vecM, paramindex, index-1, resampling_date);
447
448
449
	if(end_val!=IOUtils::nodata) {
		if(sum!=IOUtils::nodata) sum += end_val;
		else sum = end_val;
450
	} else if(strict) return IOUtils::nodata;
451

452
453
454
	return sum;
}

455
456
//do we already have this station?
size_t Accumulate::getStationIndex(const std::string& key)
457
{
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
	const size_t nr_stations = station_index.size();
	for(size_t ii=0; ii<nr_stations; ++ii) {
		if(station_index[ii]==key)
			return ii;
	}

	return IOUtils::npos;
}

//when an accumulation was measured on a period greater than the sampling period of the input data,
//for example when accumulating over 1 day but reporting values every hours, this has to be called
//to distribute the data over the extra timesteps (marked as "nodata" but this does not really
//means "nodata")
size_t Accumulate::distributeMeasurements(const std::vector<MeteoData>& vecM, const size_t& paramindex, const size_t& index, const Date& resampling_date, std::vector<MeteoData> &vecResult)
{
	vecResult.resize(0);
474

475
476
477
478
479
480
	//find the next accumulated value
	size_t end_idx = IOUtils::npos; //first index that has data
	const Date dateMaxEnd = resampling_date + measured_period;
	for(size_t idx=index; idx<vecM.size(); idx++) {
		const Date curr_date = vecM[idx].date;
		if(curr_date > dateMaxEnd) //current point would fall outside the accumulation period anyway
481
			break;
482
483
484
485
486
487
488
489
490
491
492
		if(vecM[idx](paramindex)!=IOUtils::nodata) {
			end_idx = idx;
			break;
		}
	}
	if(end_idx==IOUtils::npos) {
		std::ostringstream ss;
		ss << "Redistribution of precipitation before reaccumulation failed: precipitation value required ";
		ss << "in the " << resampling_date.toString(Date::ISO) << " - " << dateMaxEnd.toString(Date::ISO) << " interval!\n";
		throw NoAvailableDataException(ss.str(), AT);
	}
493

494
495
	const Date dateEnd = vecM[end_idx].date;
	const double precip = vecM[end_idx](paramindex);
496

497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
	//find the start date of the measured accumulation period
	const Date dateStart = vecM[end_idx].date - measured_period;
	size_t start_idx = IOUtils::npos; // last index <= (dateEnd - measured_period)
	for (size_t idx=index; idx--> 0; ) {
		const Date curr_date = vecM[idx].date;
		if(curr_date <= dateStart) {
			start_idx = idx;
			break;
		}
		if(vecM[idx](paramindex)!=IOUtils::nodata) { //this means the measured period is wrong
			const double interval = dateEnd.getJulian(true) - curr_date.getJulian(true);
			std::ostringstream ss;
			ss << "Accumulation period for \"" << parname;
			ss << "\" found to be " << interval*86400. << " for " << dateEnd.toString(Date::ISO);
			ss << " when " << measured_period*86400. << " was given as arguments of \"" << algo << "\"";
			throw IOException(ss.str(), AT);
		}
	}
	if(precip==0. && start_idx==IOUtils::npos) {
		//even if we don't know the start of the period, we can redistribute 0!
		start_idx = 0;
	}

	if(start_idx==IOUtils::npos) {
		std::ostringstream ss;
		ss << "Redistribution of precipitation before reaccumulation failed: precipitation value required ";
		ss << "at " << dateStart.toString(Date::ISO) << "!\n";
		throw NoAvailableDataException(ss.str(), AT);
	}
	if(precip!=0 && dateStart!=vecM[start_idx].date) {
		throw IOException("Currently, only multiples of the data sampling rate are allowed for accumulation periods!", AT);
	}
529

530
531
	vecResult.assign(vecM.begin()+start_idx, vecM.begin()+end_idx+1);
	//HSSweGenerator::CstDistributeHNW(precip, 1, vecResult.size()-1, paramindex, vecResult);
532
	HSSweGenerator::SmartDistributeHNW(precip, 1, vecResult.size()-1, paramindex, vecResult);
533

534
535
536
537
538
539
	return index-(start_idx);
}

void Accumulate::coreResample(const size_t& index, const size_t& paramindex,
                              const std::vector<MeteoData>& vecM, const Date& resampling_date, MeteoData& md) const
{
540
	const Date dateStart(resampling_date.getJulian() - accumulate_period, resampling_date.getTimeZone());
541
	const size_t start_idx = findStartOfPeriod(vecM, index, dateStart);
542
543
544
545
546
547
	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;
	}

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

560
561
562
563
564
565
566
567
//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)
{
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);
	if(position==ResamplingAlgorithms::begin || position==ResamplingAlgorithms::end)
		return;
568

569
570
	md(paramindex) = IOUtils::nodata;
	const Date resampling_date = md.date;
571

572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
	if(measured_period==IOUtils::nodata) {
		coreResample(index, paramindex, vecM, resampling_date, md);
	} else {
		size_t new_index = IOUtils::npos;
		const string key = md.meta.stationID;
		size_t stat_idx = getStationIndex(key);
		if(stat_idx!=IOUtils::npos) { //station is known
			const Date cache_start = vecCache[stat_idx].front().date;
			const Date cache_end = vecCache[stat_idx].back().date;
			if(resampling_date<cache_start || resampling_date>cache_end) {
				new_index = distributeMeasurements(vecM, paramindex, index, resampling_date, vecCache[stat_idx]);
			} else
				new_index = IOUtils::seek(resampling_date, vecCache[stat_idx], false);
		} else { //a new station is created
			station_index.push_back( key );
			vector<MeteoData> tmp;
			new_index = distributeMeasurements(vecM, paramindex, index, resampling_date, tmp);
			vecCache.push_back( tmp );
			stat_idx = vecCache.size()-1;
		}
592

593
594
595
596
		coreResample(new_index, paramindex, vecCache[stat_idx], resampling_date, md);
	}
}

597

598
599
600
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
601
602
const size_t Daily_solar::samples_per_day = 24*3; //every 20 minutes

603
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)
604
            : ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs), radiation(), station_index(), dateStart(), dateEnd(), loss_factor()
605
606
607
608
609
610
611
612
613
{
	const size_t nr_args = vecArgs.size();
	if(nr_args>0) {
		throw InvalidArgumentException("Too many arguments for \""+i_parname+"::"+i_algoname+"\"", AT);
	}
}

std::string Daily_solar::toString() const
{
614
	ostringstream ss;
615
616
617
618
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	return ss.str();
}

619
620
//look for the daily sum of solar radiation for the current day
size_t Daily_solar::getNearestValidPt(const std::vector<MeteoData>& vecM, const size_t& paramindex,  const size_t& stat_idx, const size_t& pos) const
621
622
623
624
625
626
{
	size_t indexP1=IOUtils::npos;
	size_t indexP2=IOUtils::npos;

	//look for daily sum before the current point
	for (size_t ii=pos; ii-- >0; ) {
627
		if (vecM[ii].date < dateStart[stat_idx]) break;
628
629
630
631
632
633
634
635
		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) {
636
		if (vecM[ii].date > dateEnd[stat_idx]) break;
637
638
639
640
641
642
		if (vecM[ii](paramindex) != IOUtils::nodata) {
			indexP2 = ii;
			break;
		}
	}

643
644
645
646
	if(indexP1!=IOUtils::npos && indexP2!=IOUtils::npos) {
		const string msg = "More than one daily sum of solar radiation found between "+dateStart[stat_idx].toString(Date::ISO)+" and "+dateEnd[stat_idx].toString(Date::ISO);
		throw IOException(msg, AT);
	}
647
648
649
650
651
652

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

653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
//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;
	for(Date date(dateStart[stat_idx]); date<dateEnd[stat_idx]; date += 1./double(samples_per_day)) {
		//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();
	for(size_t ii=0; ii<nr_stations; ++ii) {
		if(station_index[ii]==key)
			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::setDayStartAndEnd(const Date& resampling_date, const size_t& stat_idx)
{
	dateStart[stat_idx] = resampling_date;
	dateStart[stat_idx].rnd(24*3600, Date::DOWN);
	if(dateStart[stat_idx]==resampling_date) //if resampling_date=midnight GMT, the rounding lands on the exact same date
		dateStart[stat_idx] -= 1.;

	dateEnd[stat_idx] = resampling_date;
	dateEnd[stat_idx].rnd(24*3600, Date::UP);
}

void Daily_solar::resample(const size_t& index, const ResamplingPosition& /*position*/, const size_t& paramindex,
                           const std::vector<MeteoData>& vecM, MeteoData& md)
726
727
728
729
730
731
732
733
734
735
736
{
	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 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;
737
	const double HS = md(MeteoData::HS);
738

739
740
	//get station index
	const size_t stat_idx = getStationIndex(md.meta.stationID);
741

742
743
744
745
746
747
748
749
750
751
752
	//has the radiation already been calculated for this day and station?
	if(md.date<dateStart[stat_idx] || md.date>=dateEnd[stat_idx]) {
		setDayStartAndEnd(md.date, stat_idx);
		const size_t indexP = getNearestValidPt(vecM, paramindex, stat_idx, index);
		if(indexP==IOUtils::npos) { //no daily sum found for the current day
			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...
753
754
	}

755
756
757
758
759
760
761
762
763
764
765
766
767
	if(loss_factor[stat_idx]==IOUtils::nodata) //the station could not be calculated for this day
		return;

	//interpolate radiation for this timestep and write it out
	const double rad = getSolarInterpol(md.date, stat_idx);
	if(paramindex==MeteoData::ISWR) {
		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;
	}
768

769
	md.setResampled(true);
770
771
}

772
773
} //namespace