WSL/SLF GitLab Repository

ResamplingAlgorithms.cc 26 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
20
21
#include <meteoio/MathOptim.h>
#include <meteoio/meteolaws/Atmosphere.h>
#include <meteoio/meteolaws/Sun.h>
22
#include <meteoio/meteostats/libinterpol1D.h>
23
24
#include <meteoio/meteofilters/ProcHNWDistribute.h> //for the precipitation distribution

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
36
{
	const std::string algoname(IOUtils::strToUpper(i_algoname));

	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 == "DAILY_SOLAR"){
		return new Daily_solar(algoname, parname, dflt_window_size, vecArgs);
46
47
48
49
	} else {
		throw IOException("The resampling algorithm '"+algoname+"' is not implemented" , AT);
	}
}
50

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

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

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

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

84
	return valend - left_accumulation;
85
86
}

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

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

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

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

149
150
151
152
153
154
155
156
157
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
{
158
	ostringstream ss;
159
160
161
162
163
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	ss << "[ ]";
	return ss.str();
}

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

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

177
178
179
	return;
}

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
205
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
{
206
	ostringstream ss;
207
208
209
210
211
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	ss << "[ window_size=" << window_size << " extrapolate=" << extrapolate << " ]";
	return ss.str();
}

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

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

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

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

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

243
		if (IOUtils::checkEpsilonEquality(diff1.getJulian(true), diff2.getJulian(true), 0.1/1440.)){ //within 6 seconds
244
			md(paramindex) = Interpol1D::weightedMean(val1, val2, 0.5);
245
		} else if (diff1 < diff2){
246
			md(paramindex) = val1;
247
		} else if (diff1 > diff2){
248
			md(paramindex) = val2;
249
		}
250
	} else if (extrapolate) {
251
252
253
254
		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);
255
256
		} else { // no nearest neighbour with a value different from IOUtils::nodata
			return;
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
286
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
{
287
	ostringstream ss;
288
289
290
291
292
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	ss << "[ window_size=" << window_size << " extrapolate=" << extrapolate << " ]";
	return ss.str();
}

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

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

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

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

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

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

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

355

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

364
365
	for(size_t ii=0; ii<nr_args; ii++) {
		if(IOUtils::isNumeric(vecArgs[ii])) {
366
367
368
			IOUtils::convertString(accumulate_period, vecArgs[ii]);
			accumulate_period /= 86400.; //user uses seconds, internally julian day is used
			if(accumulate_period<=0.) {
369
370
371
372
373
374
375
376
377
				std::ostringstream ss;
				ss << "Invalid accumulation period (" << accumulate_period << ") for \"" << i_parname << "::" << i_algoname << "\"";
				throw InvalidArgumentException(ss.str(), AT);
			}
		} 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);
378
	}
379
380
}

381
382
std::string Accumulate::toString() const
{
383
	ostringstream ss;
384
385
386
387
388
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	ss << "[ accumulate_period=" << accumulate_period << " strict=" << strict << " ]";
	return ss.str();
}

389
390
//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)
391
{
392
	size_t start_idx = IOUtils::npos;
393
	for (size_t idx=index; idx--> 0; ) {
394
395
396
		const Date curr_date = vecM[idx].date;
		if(curr_date <= dateStart) {
			start_idx = idx;
397
398
399
			break;
		}
	}
400

401
402
	return start_idx;
}
403

404
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
405
{//to keep in mind: start_idx is last index <= dateStart and index is first index >= resampling_date
406
407
408
	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);
409
410

	if(start_val!=IOUtils::nodata && end_val!=IOUtils::nodata)
411
		sum = end_val - start_val;
412

413
414
415
	return sum;
}

416
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
417
{//to keep in mind: start_idx is last index <= dateStart and index is first index >= resampling_date
418
	double sum = IOUtils::nodata;
419
	//resample begining point, in the [start_idx ; start_idx+1] interval
420
421
422
423
	const double start_value = partialAccumulateAtRight(vecM, paramindex, start_idx, dateStart);
	if(start_value!=IOUtils::nodata)
		sum=start_value;
	else if(strict) return IOUtils::nodata;
424

425
426
	//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++) {
427
428
429
430
		const double curr_value = vecM[idx](paramindex);
		if(curr_value!=IOUtils::nodata) {
			if(sum!=IOUtils::nodata) sum += curr_value;
			else sum = curr_value;
431
		} else if(strict) return IOUtils::nodata;
432
	}
433

434
	//resample end point, in the [index-1 ; index] interval
435
	const double end_val = partialAccumulateAtLeft(vecM, paramindex, index-1, resampling_date);
436
437
438
	if(end_val!=IOUtils::nodata) {
		if(sum!=IOUtils::nodata) sum += end_val;
		else sum = end_val;
439
	} else if(strict) return IOUtils::nodata;
440

441
442
443
	return sum;
}

444
445
446
//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)
447
{
448
449
450
451
	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;
452

453
454
	md(paramindex) = IOUtils::nodata;
	const Date resampling_date = md.date;
455

456
	const Date dateStart(resampling_date.getJulian() - accumulate_period, resampling_date.getTimeZone());
457
	const size_t start_idx = findStartOfPeriod(vecM, index, dateStart);
458
459
460
461
462
463
	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;
	}

464
	if((index - start_idx) <= 1) {//easy upsampling when start & stop are in the same input time step
465
		//upsampling (for example, generate 15min values from hourly data)
466
		const double sum = easySampling(vecM, paramindex, index, start_idx, dateStart, resampling_date);
467
		md(paramindex) = sum; //if resampling was unsuccesful, sum==IOUtils::nodata
468
	} else {
469
		//downsampling (for example, generate daily values from hourly data)
470
471
		//and upsampling when resampled period falls accross a measurement timestamp
		const double sum = complexSampling(vecM, paramindex, index, start_idx, dateStart, resampling_date);
472
473
		md(paramindex) = sum; //if resampling was unsuccesful, sum==IOUtils::nodata
	}
474
475
}

476

477
478
479
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
480
481
const size_t Daily_solar::samples_per_day = 24*3; //every 20 minutes

482
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)
483
            : ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs), radiation(), station_index(), dateStart(), dateEnd(), loss_factor()
484
485
486
487
488
489
490
491
492
{
	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
{
493
	ostringstream ss;
494
495
496
497
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	return ss.str();
}

498
499
//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
500
501
502
503
504
505
{
	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; ) {
506
		if (vecM[ii].date < dateStart[stat_idx]) break;
507
508
509
510
511
512
513
514
		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) {
515
		if (vecM[ii].date > dateEnd[stat_idx]) break;
516
517
518
519
520
521
		if (vecM[ii](paramindex) != IOUtils::nodata) {
			indexP2 = ii;
			break;
		}
	}

522
523
524
525
	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);
	}
526
527
528
529
530
531

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

532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
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
//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)
605
606
607
608
609
610
611
612
613
614
615
{
	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;
616
	const double HS = md(MeteoData::HS);
617

618
619
	//get station index
	const size_t stat_idx = getStationIndex(md.meta.stationID);
620

621
622
623
624
625
626
627
628
629
630
631
	//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...
632
633
	}

634
635
636
637
638
639
640
641
642
643
644
645
646
	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;
	}
647

648
	md.setResampled(true);
649
650
}

651
652
} //namespace