WSL/SLF GitLab Repository

ResamplingAlgorithms.cc 23.6 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
#include <cmath>
24
#include <algorithm>
25
26
27
28
29

using namespace std;

namespace mio {

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

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

49
50
double ResamplingAlgorithms::funcval(size_t pos, const size_t& paramindex, const std::vector<MeteoData>& vecM,
                                     const Date& date, const bool& start_pt)
51
{
52
53
54
55
56
57
58
59
60
61
62
63
	if(!start_pt) pos--;
	const double valstart = vecM[pos](paramindex);
	if (vecM[pos].date == date) return valstart;

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

65
66
67
68
	if(start_pt)
		return valend - linearInterpolation(jul1, 0., jul2, valend, date.getJulian(true));
	else
		return linearInterpolation(jul1, 0., jul2, valend, date.getJulian(true));
69
70
}

71
72
73
74
75
/**
 * @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
76
 * @param resampling_date date to resample
77
78
79
80
81
 * @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,
82
                                              const double& window_size, size_t& indexP1, size_t& indexP2) //HACK
83
{
84
85
	indexP1=IOUtils::npos;
	indexP2=IOUtils::npos;
86

87
88
	const Date dateStart = resampling_date - window_size;

89
	for (size_t ii=pos; ii-- >0; ) {
90
91
92
93
94
95
96
97
98
		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;
99

100
	for (size_t ii=pos; ii<vecM.size(); ++ii) {
101
102
103
104
105
106
107
108
109
110
111
112
113
114
		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
115
 * @param x x-coordinate of desired point
116
117
118
 * @return y-coordinate of desired point
 */
double ResamplingAlgorithms::linearInterpolation(const double& x1, const double& y1,
119
                                       const double& x2, const double& y2, const double& x)
120
121
122
123
124
125
126
127
{
	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;

128
	return (a*x + b);
129
130
}

131
/**********************************************************************************
132
 * The following functions are implementations of different resampling algorithms *
133
 **********************************************************************************/
134

135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
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
{
	stringstream ss;
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	ss << "[ ]";
	return ss.str();
}

150
void NoResampling::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
151
                            const std::vector<MeteoData>& vecM, MeteoData& md)
152
{
153
154
155
156
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);

	if (position == ResamplingAlgorithms::exact_match) {
157
		const double value = vecM[index](paramindex);
158
159
160
161
162
		if (value != IOUtils::nodata) {
			md(paramindex) = value; //propagate value
		}
	}

163
164
165
	return;
}

166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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
{
	stringstream ss;
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	ss << "[ window_size=" << window_size << " extrapolate=" << extrapolate << " ]";
	return ss.str();
}

198
void NearestNeighbour::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
199
                                const std::vector<MeteoData>& vecM, MeteoData& md)
200
{
201
202
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);
203

204
	if (position == ResamplingAlgorithms::exact_match) {
205
		const double value = vecM[index](paramindex);
206
		if (value != IOUtils::nodata) {
207
208
209
210
			md(paramindex) = value; //propagate value
			return;
		}
	}
211

212
213
	//if we are at the very beginning or end of vecM and !extrapolate, then there's nothing to do
	if (((!extrapolate) && (position == ResamplingAlgorithms::end))
214
	    || ((!extrapolate) && (position == ResamplingAlgorithms::begin)))
215
		return;
216

217
	const Date resampling_date = md.date;
218
	size_t indexP1=IOUtils::npos, indexP2=IOUtils::npos;
219
	getNearestValidPts(index, paramindex, vecM, resampling_date, window_size, indexP1, indexP2);
220
	const bool foundP1=(indexP1!=IOUtils::npos), foundP2=(indexP2!=IOUtils::npos);
221

222
	//Try to find the nearest neighbour, if there are two equally distant, then return the arithmetic mean
223
	if (foundP1 && foundP2) { //standard behavior
224
225
		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
226
227
		const double val1 = vecM[indexP1](paramindex);
		const double val2 = vecM[indexP2](paramindex);
228

229
		if (IOUtils::checkEpsilonEquality(diff1.getJulian(true), diff2.getJulian(true), 0.1/1440.)){ //within 6 seconds
230
			md(paramindex) = Interpol1D::weightedMean(val1, val2, 0.5);
231
		} else if (diff1 < diff2){
232
			md(paramindex) = val1;
233
		} else if (diff1 > diff2){
234
			md(paramindex) = val2;
235
		}
236
	} else if (extrapolate) {
237
238
239
240
		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);
241
242
		} else { // no nearest neighbour with a value different from IOUtils::nodata
			return;
243
244
245
246
		}
	}
}

247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
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
{
	stringstream ss;
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	ss << "[ window_size=" << window_size << " extrapolate=" << extrapolate << " ]";
	return ss.str();
}

279
void LinearResampling::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
280
                                const std::vector<MeteoData>& vecM, MeteoData& md)
281
{
282
283
	if (index >= vecM.size())
		throw IOException("The index of the element to be resampled is out of bounds", AT);
284

285
	if (position == ResamplingAlgorithms::exact_match) {
286
		const double value = vecM[index](paramindex);
287
		if (value != IOUtils::nodata) {
288
289
290
291
			md(paramindex) = value; //propagate value
			return;
		}
	}
292

293
294
	//if we are at the very beginning or end of vecM and !extrapolate, then there's nothing to do
	if (((!extrapolate) && (position == ResamplingAlgorithms::end))
295
	    || ((!extrapolate) && (position == ResamplingAlgorithms::begin)))
296
		return;
297

298
	const Date resampling_date = md.date;
299
	size_t indexP1=IOUtils::npos, indexP2=IOUtils::npos;
300
	getNearestValidPts(index, paramindex, vecM, resampling_date, window_size, indexP1, indexP2);
301
	bool foundP1=(indexP1!=IOUtils::npos), foundP2=(indexP2!=IOUtils::npos);
302
303

	//do nothing if we can't interpolate, and extrapolation is not explicitly activated
304
	if ((!extrapolate) && ((!foundP1) || (!foundP2)))
305
306
307
308
309
310
		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
311
	if (!foundP1 && foundP2){ //only nodata values found before index, try looking after indexP2
312
		for (size_t ii=indexP2+1; ii<vecM.size(); ii++){
313
			if (vecM[ii](paramindex) != IOUtils::nodata){
314
315
316
317
				indexP1 = ii;
				foundP1 = true;
				break;
			}
318
		}
319
	} else if (foundP1 && !foundP2){ //only nodata found after index, try looking before indexP1
320
		for (size_t ii=indexP1; (ii--) > 0; ){
321
			if (vecM[ii](paramindex) != IOUtils::nodata){
322
323
324
325
				indexP2=ii;
				foundP2 = true;
				break;
			}
326
		}
327
328
329
330
331
332
	}

	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
333
	const double val1 = vecM[indexP1](paramindex);
334
	const double jul1 = vecM[indexP1].date.getJulian(true);
335
	const double val2 = vecM[indexP2](paramindex);
336
	const double jul2 = vecM[indexP2].date.getJulian(true);
337

338
	md(paramindex) = linearInterpolation(jul1, val1, jul2, val2, resampling_date.getJulian(true));
339
340
}

341
342
Accumulate::Accumulate(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), accumulate_period(0.), strict(false)
343
{
344
345
	const size_t nr_args = vecArgs.size();
	if (nr_args==1 || nr_args==2) {
346
		IOUtils::convertString(accumulate_period, vecArgs[0]);
347
		accumulate_period /= 86400.; //user uses seconds, internally julian day is used
348
349
		if(accumulate_period<=0.) {
			std::stringstream tmp;
350
			tmp << "Invalid accumulation period (" << accumulate_period << ") for \"" << i_parname << "::" << i_algoname << "\"";
351
352
			throw InvalidArgumentException(tmp.str(), AT);
		}
353
354
355
356
357
		if(nr_args==2) {
			if(vecArgs[1]=="strict")
				strict=true;
			else
				throw InvalidArgumentException("Invalid argument \""+vecArgs[1]+"\" for \""+i_parname+"::"+i_algoname+"\"", AT);
358
		}
359
	} else {
360
		throw InvalidArgumentException("Please provide accumulation period (in seconds) for \""+i_parname+"::"+i_algoname+"\"", AT);
361
	}
362
363
}

364
365
366
367
368
369
370
371
std::string Accumulate::toString() const
{
	stringstream ss;
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	ss << "[ accumulate_period=" << accumulate_period << " strict=" << strict << " ]";
	return ss.str();
}

372
void Accumulate::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
373
                          const std::vector<MeteoData>& vecM, MeteoData& md)
374
375
376
377
378
379
380
381
{
	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;

	const Date resampling_date = md.date;
	md(paramindex) = IOUtils::nodata;
382

383
	//find start of accumulation period and initialize the sum
384
	double sum = IOUtils::nodata;
385
	const Date dateStart(resampling_date.getJulian() - accumulate_period, resampling_date.getTimeZone());
386
	bool found_start=false;
387
	size_t start_idx; //this is the index of the first data point before the window
388
	for (start_idx=index+1; (start_idx--) > 0; ) {
389
		const Date date = vecM[start_idx].date;
390
391
392
		if(date <= dateStart) {
			if(date<dateStart) {
				const double start_value = funcval(start_idx, paramindex, vecM, dateStart, true); //resampling the starting point
393
394
395
				if(start_value!=IOUtils::nodata)
					sum=start_value;
				else if(strict) return;
396
			}
397
398
399
400
			found_start=true;
			break;
		}
	}
401
402
	if (!found_start) {
		cerr << "[W] Could not accumulate " << vecM.at(0).getNameForParameter(paramindex) << ": ";
403
404
		cerr << "not enough data for accumulation period at date " << resampling_date.toString(Date::ISO) << "\n";
		md(paramindex) = IOUtils::nodata;
405
		return;
406
	}
407
408
409
410
411
412
413
414
	if(vecM[start_idx].date != dateStart) start_idx++; //we need to skip the first point that was already used in the interpolation
	//if up-sampling, take a quicker path (for example, generate 15min values from hourly data)
	if(start_idx==index) {
		const double start_val = funcval(start_idx, paramindex, vecM, dateStart, false);
		const double end_val = funcval(index, paramindex, vecM, resampling_date, false);
		if(start_val!=IOUtils::nodata && end_val!=IOUtils::nodata) md(paramindex) = end_val - start_val;
		return;
	}
415

416
417
418
419
420
421
	 //sum all whole periods
	for(size_t idx=(start_idx+1); idx<index; idx++) {
		const double curr_value = vecM[idx](paramindex);
		if(curr_value!=IOUtils::nodata) {
			if(sum!=IOUtils::nodata) sum += curr_value;
			else sum = curr_value;
422
		} else if(strict) return;
423
	}
424

425
426
427
428
429
	//resample end point
	const double end_val = funcval(index, paramindex, vecM, resampling_date, false);
	if(end_val!=IOUtils::nodata) {
		if(sum!=IOUtils::nodata) sum += end_val;
		else sum = end_val;
430
	} else if(strict) return;
431

432
	//write out sum
433
	md(paramindex) = sum;
434
435
}

436
437
438
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
439
440
const size_t Daily_solar::samples_per_day = 24*3; //every 20 minutes

441
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)
442
            : ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs), radiation(), station_index(), dateStart(), dateEnd(), loss_factor()
443
444
445
446
447
448
449
450
451
452
453
454
455
456
{
	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
{
	stringstream ss;
	ss << right << setw(10) << parname << "::"  << left << setw(15) << algo;
	return ss.str();
}

457
458
//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
459
460
461
462
463
464
{
	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; ) {
465
		if (vecM[ii].date < dateStart[stat_idx]) break;
466
467
468
469
470
471
472
473
		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) {
474
		if (vecM[ii].date > dateEnd[stat_idx]) break;
475
476
477
478
479
480
		if (vecM[ii](paramindex) != IOUtils::nodata) {
			indexP2 = ii;
			break;
		}
	}

481
482
483
484
	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);
	}
485
486
487
488
489
490

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

491
492
493
494
495
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
529
530
531
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
//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)
564
565
566
567
568
569
570
571
572
573
574
{
	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;
575
	const double HS = md(MeteoData::HS);
576

577
578
	//get station index
	const size_t stat_idx = getStationIndex(md.meta.stationID);
579

580
581
582
583
584
585
586
587
588
589
590
	//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...
591
592
	}

593
594
595
596
597
598
599
600
601
602
603
604
605
	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;
	}
606

607
	md.setResampled(true);
608
609
}

610
611
} //namespace