WSL/SLF GitLab Repository

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

using namespace std;

namespace mio {

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

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

48
49
double ResamplingAlgorithms::funcval(size_t pos, const size_t& paramindex, const std::vector<MeteoData>& vecM,
                                     const Date& date, const bool& start_pt)
50
{
51
52
53
54
55
56
57
58
59
60
61
62
	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);
63

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

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

86
87
	const Date dateStart = resampling_date - window_size;

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

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

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

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

134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
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();
}

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

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

162
163
164
	return;
}

165
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
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();
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

363
364
365
366
367
368
369
370
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();
}

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

382
	//find start of accumulation period and initialize the sum
383
	double sum = IOUtils::nodata;
384
	const Date dateStart(resampling_date.getJulian() - accumulate_period, resampling_date.getTimeZone());
385
	bool found_start=false;
386
	size_t start_idx; //this is the index of the first data point before the window
387
	for (start_idx=index+1; (start_idx--) > 0; ) {
388
		const Date date = vecM[start_idx].date;
389
390
391
		if(date <= dateStart) {
			if(date<dateStart) {
				const double start_value = funcval(start_idx, paramindex, vecM, dateStart, true); //resampling the starting point
392
393
394
				if(start_value!=IOUtils::nodata)
					sum=start_value;
				else if(strict) return;
395
			}
396
397
398
399
			found_start=true;
			break;
		}
	}
400
401
	if (!found_start) {
		cerr << "[W] Could not accumulate " << vecM.at(0).getNameForParameter(paramindex) << ": ";
402
403
		cerr << "not enough data for accumulation period at date " << resampling_date.toString(Date::ISO) << "\n";
		md(paramindex) = IOUtils::nodata;
404
		return;
405
	}
406
407
408
409
410
411
412
413
	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;
	}
414

415
416
417
418
419
420
	 //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;
421
		} else if(strict) return;
422
	}
423

424
425
426
427
428
	//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;
429
	} else if(strict) return;
430

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

435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
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
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
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)
            : ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs)
{
	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();
}

size_t Daily_solar::getNearestValidPt(const size_t& pos, const size_t& paramindex, const std::vector<MeteoData>& vecM, const Date& resampling_date)
{
	size_t indexP1=IOUtils::npos;
	size_t indexP2=IOUtils::npos;

	//look for daily sum before the current point
	Date dateStart( resampling_date );
	dateStart.rnd(24*3600, Date::DOWN);
	if(dateStart==resampling_date) //if resampling_date=midnight GMT, the rounding lands on the exact same date
		dateStart -= 1.;
	for (size_t ii=pos; ii-- >0; ) {
		if (vecM[ii].date < dateStart) break;
		if (vecM[ii](paramindex) != IOUtils::nodata){
			indexP1 = ii;
			break;
		}
	}

	//look for daily sum after the current point
	Date dateEnd( resampling_date );
	dateEnd.rnd(24*3600, Date::UP);
	for (size_t ii=pos; ii<vecM.size(); ++ii) {
		if (vecM[ii].date > dateEnd) break;
		if (vecM[ii](paramindex) != IOUtils::nodata) {
			indexP2 = ii;
			break;
		}
	}

	if(indexP1!=IOUtils::npos && indexP2!=IOUtils::npos)
		throw IOException("More than one daily sum of solar radiation found for "+resampling_date.toString(Date::ISO)+" !", AT);

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

void Daily_solar::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
                           const std::vector<MeteoData>& vecM, MeteoData& md) const
{
	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);

	md.setResampled(true);

	const size_t indexP = getNearestValidPt(index, paramindex, vecM, md.date);
	if(indexP==IOUtils::npos) //no daily sum found for the current day
		return;

	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;

	double P = md(MeteoData::P);
	double TA = md(MeteoData::TA);
	double RH = md(MeteoData::RH);
	double HS = md(MeteoData::HS);
	double albedo = 0.5;

	if (P==IOUtils::nodata) P = Atmosphere::stdAirPressure(alt);
	if (HS!=IOUtils::nodata) //no big deal if we can not adapt the albedo
		albedo = (HS>=snow_thresh)? snow_albedo : soil_albedo;
	if(TA==IOUtils::nodata || RH==IOUtils::nodata) {
		//set TA & RH so the reduced precipitable water will get an average value
		TA=274.98;
		RH=0.666;
	}

	SunObject sun(lat, lon, alt, md.date.getJulian(), md.date.getTimeZone());
	sun.calculateRadiation(TA, RH, P, albedo);
	double toa, direct, diffuse;
	sun.getHorizontalRadiation(toa, direct, diffuse);
	const double sum_toa_h = sun.approxTOADailySum();
	//const double loss_factor = (sum_toa_h>0.)? vecM[indexP](paramindex) / sum_toa_h : 0.;
	const double loss_factor = 1.;

	if(paramindex==MeteoData::ISWR)
		md(paramindex) = loss_factor * (toa);
	else
		md(paramindex) = loss_factor * (toa) * albedo;
}

540
541
} //namespace