WSL/SLF GitLab Repository

InterpolationAlgorithms.cc 24.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
/***********************************************************************************/
/*  Copyright 2010 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/InterpolationAlgorithms.h>
19
#include <meteoio/meteolaws/Atmosphere.h>
20
21
#include <meteoio/Meteo2DInterpolator.h>
#include <meteoio/IOManager.h>
22
#include <meteoio/MathOptim.h>
23
24
25
26
27

using namespace std;

namespace mio {

28
InterpolationAlgorithm* AlgorithmFactory::getAlgorithm(const std::string& i_algoname,
29
                                                       Meteo2DInterpolator& i_mi,
30
                                                       const std::vector<std::string>& i_vecArgs, IOManager& iom)
31
{
32
	const std::string algoname( IOUtils::strToUpper(i_algoname) );
33

34
	if (algoname == "CST"){// constant fill
35
		return new ConstAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
36
	} else if (algoname == "STD_PRESS"){// standard air pressure interpolation
37
		return new StandardPressureAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
38
	} else if (algoname == "CST_LAPSE"){// constant fill with an elevation lapse rate
39
		return new ConstLapseRateAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
40
	} else if (algoname == "IDW"){// Inverse Distance Weighting fill
41
		return new IDWAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
42
	} else if (algoname == "IDW_LAPSE"){// Inverse Distance Weighting with an elevation lapse rate fill
43
		return new IDWLapseAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
44
	} /*else if (algoname == "LIDW_LAPSE"){// Inverse Distance Weighting with an elevation lapse rate fill, restricted to a local scale
45
		return new LocalIDWLapseAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
46
	}*/ else if (algoname == "RH"){// relative humidity interpolation
47
		return new RHAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
48
	} else if (algoname == "ILWR"){// long wave radiation interpolation
49
		return new ILWRAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
50
	} else if (algoname == "WIND_CURV"){// wind velocity interpolation (using a heuristic terrain effect)
51
		return new SimpleWindInterpolationAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
52
	} else if (algoname == "ODKRIG"){// ordinary kriging
53
		return new OrdinaryKrigingAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
54
	} else if (algoname == "ODKRIG_LAPSE"){// ordinary kriging with lapse rate
55
		return new LapseOrdinaryKrigingAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
56
	} else if (algoname == "USER"){// read user provided grid
57
		return new USERInterpolation(i_mi, i_vecArgs, i_algoname, iom);
58
	} else if (algoname == "HNW_SNOW"){// precipitation interpolation according to (Magnusson, 2010)
59
  		return new SnowHNWInterpolation(i_mi, i_vecArgs, i_algoname, iom);
60
61
62
63
64
	} else {
		throw IOException("The interpolation algorithm '"+algoname+"' is not implemented" , AT);
	}
}

65
66
size_t InterpolationAlgorithm::getData(const Date& i_date, const MeteoData::Parameters& i_param,
                                             std::vector<double>& o_vecData)
67
{
68
	iomanager.getMeteoData(i_date, vecMeteo);
69
70
71
72
73
74
75
76
77
	o_vecData.clear();
	for (size_t ii=0; ii<vecMeteo.size(); ii++){
		const double& val = vecMeteo[ii](i_param);
		if (val != IOUtils::nodata) {
			o_vecData.push_back(val);
		}
	}

	return o_vecData.size();
78
79
}

80
81
size_t InterpolationAlgorithm::getData(const Date& i_date, const MeteoData::Parameters& i_param,
                                             std::vector<double>& o_vecData, std::vector<StationData>& o_vecMeta)
82
{
83
	iomanager.getMeteoData(i_date, vecMeteo);
84
85
	o_vecData.clear();
	o_vecMeta.clear();
86
	for (size_t ii=0; ii<vecMeteo.size(); ii++){
87
		const double& val = vecMeteo[ii](i_param);
88
		if (val != IOUtils::nodata){
89
90
			o_vecData.push_back(val);
			o_vecMeta.push_back(vecMeteo[ii].meta);
91
92
93
		}
	}

94
	return o_vecData.size();
95
96
}

97
98
size_t InterpolationAlgorithm::getStationAltitudes(const std::vector<StationData>& i_vecMeta,
                                                         std::vector<double>& o_vecData) const
99
{
100
	o_vecData.clear();
101
	for (size_t ii=0; ii<i_vecMeta.size(); ii++){
102
103
104
105
		const double& alt = i_vecMeta[ii].position.getAltitude();
		if (alt != IOUtils::nodata) {
			o_vecData.push_back(alt);
		}
106
107
	}

108
	return o_vecData.size();
109
110
}

111
112
113
114
115
/**
 * @brief Return an information string about the interpolation process
 * @return string containing some information (algorithm used, number of stations)
*/
std::string InterpolationAlgorithm::getInfo() const
116
{
117
	std::ostringstream os;
118
119
120
121
122
	os << algo << ", " << nrOfMeasurments;
	if(nrOfMeasurments==1)
		os << " station";
	else
		os << " stations";
123
	const std::string tmp( info.str() );
124
125
126
127
	if(!tmp.empty()) {
		os << ", " << tmp;
	}
	return os.str();
128
129
}

130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
/**
 * @brief Read the interpolation arguments and compute the trend accordingly
 *
 * @param vecAltitudes altitudes sorted similarly as the data in vecDat
 * @param vecDat data for the interpolated parameter
 * @param trend object containing the fitted trend to be used for detrending/retrending
*/
void InterpolationAlgorithm::getTrend(const std::vector<double>& vecAltitudes, const std::vector<double>& vecDat, Fit1D &trend) const
{
	bool status;
	if (vecArgs.empty()) {
		trend.setModel(Fit1D::NOISY_LINEAR, vecAltitudes, vecDat, false);
		status = trend.fit();
	} else if (vecArgs.size() == 1) {
		double lapse_rate;
		IOUtils::convertString(lapse_rate, vecArgs.front());
		trend.setModel(Fit1D::NOISY_LINEAR, vecAltitudes, vecDat, false);
		trend.setLapseRate(lapse_rate);
		status = trend.fit();
	} else if (vecArgs.size() == 2) {
		std::string extraArg;
		IOUtils::convertString(extraArg, vecArgs[1]);
		if(extraArg=="soft") { //soft
			trend.setModel(Fit1D::NOISY_LINEAR, vecAltitudes, vecDat, false);
			status = trend.fit();
			if(!status) {
				double lapse_rate;
				IOUtils::convertString(lapse_rate, vecArgs[0]);
				trend.setModel(Fit1D::NOISY_LINEAR, vecAltitudes, vecDat, false);
				trend.setLapseRate(lapse_rate);
				status = trend.fit();
			}
		} else if(extraArg=="frac") {
			double lapse_rate;
			IOUtils::convertString(lapse_rate, vecArgs[0]);
			trend.setModel(Fit1D::NOISY_LINEAR, vecAltitudes, vecDat, false);
			const double avgData = Interpol1D::arithmeticMean(vecDat);
			trend.setLapseRate(lapse_rate*avgData);
			status = trend.fit();
		} else {
			throw InvalidArgumentException("Unknown argument \""+extraArg+"\" supplied for the "+algo+" algorithm", AT);
		}
	} else { //incorrect arguments, throw an exception
		throw InvalidArgumentException("Wrong number of arguments supplied for the "+algo+" algorithm", AT);
	}

	if(!status)
		throw IOException("Interpolation FAILED for parameter " + MeteoData::getParameterName(param) + ": " + trend.getInfo(), AT);
}


void InterpolationAlgorithm::detrend(const Fit1D& trend, const std::vector<double>& vecAltitudes, std::vector<double> &vecDat) const
{
	if(vecDat.size() != vecAltitudes.size()) {
		std::ostringstream ss;
		ss << "Number of station data (" << vecDat.size() << ") and number of elevations (" << vecAltitudes.size() << ") don't match!";
		throw InvalidArgumentException(ss.str(), AT);
	}

	for(size_t ii=0; ii<vecAltitudes.size(); ii++) {
		double &val = vecDat[ii];
		if(val!=IOUtils::nodata)
			val -= trend( vecAltitudes[ii] );
	}
}

196
void InterpolationAlgorithm::retrend(const DEMObject& dem, const Fit1D& trend, Grid2DObject &grid) const
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
{
	const size_t nxy = grid.getNx()*grid.getNy();
	const size_t dem_nxy = dem.grid2D.getNx()*dem.grid2D.getNy();
	if(nxy != dem_nxy) {
		std::ostringstream ss;
		ss << "Dem size (" << dem.grid2D.getNx() << "," << dem.grid2D.getNy() << ") and";
		ss << "grid size (" << grid.getNx() << "," << grid.getNy() << ") don't match!";
		throw InvalidArgumentException(ss.str(), AT);
	}

	for(size_t ii=0; ii<nxy; ii++) {
		double &val = grid(ii);
		if(val!=IOUtils::nodata)
			val += trend.f( dem.grid2D(ii) );
	}
}

/**********************************************************************************/
/*                    Implementation of the various algorithms                    */
216
/**********************************************************************************/
217
double StandardPressureAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
218
219
{
	date = i_date;
220
	param = in_param;
221
	nrOfMeasurments = getData(date, param, vecData, vecMeta);
222

223
224
225
226
227
228
229
230
231
	if (param != MeteoData::P)
		return 0.0;

	if (nrOfMeasurments == 0)
		return 1.0;

	return 0.1;
}

232
void StandardPressureAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid) {
233
	Interpol2D::stdPressure(dem, grid);
234
235
}

236
double ConstAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
237
238
{
	date = i_date;
239
	param = in_param;
240
	nrOfMeasurments = getData(date, param, vecData, vecMeta);
241
242
243
244
245
246
247
248
249
250
251
252

	if (nrOfMeasurments == 0){
		return 0.0;
	} else if (nrOfMeasurments == 1){
		return 0.8;
	} else if (nrOfMeasurments > 1){
		return 0.2;
	}

	return 0.0;
}

253
void ConstAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid) {
254
	Interpol2D::constant(Interpol1D::arithmeticMean(vecData), dem, grid);
255
256
}

257
double ConstLapseRateAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
258
259
{
	date = i_date;
260
	param = in_param;
261
	nrOfMeasurments = getData(date, param, vecData, vecMeta);
262
263
264
265

	if (nrOfMeasurments == 0){
		return 0.0;
	} else if (nrOfMeasurments == 1){
266
		if (!vecArgs.empty())
267
268
269
270
271
272
273
			return 0.9; //the laspe rate is provided
		else
			return 0.0; //no lapse rate is provided and it can not be computed
	} else if (nrOfMeasurments == 2){
		//in any case, we can do at least as good as IDW_LAPSE
		return 0.71;
	} else if (nrOfMeasurments>2){
274
275
276
277
278
279
		return 0.2;
	}

	return 0.2;
}

280
void ConstLapseRateAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
281
{
282
	info.clear(); info.str("");
283
	vector<double> vecAltitudes;
284
	getStationAltitudes(vecMeta, vecAltitudes);
285
	if (vecAltitudes.empty())
286
		throw IOException("Not enough data for spatially interpolating parameter " + MeteoData::getParameterName(param), AT);
287

288
289
290
291
	Fit1D trend;
	getTrend(vecAltitudes, vecData, trend);
	info << trend.getInfo();
	detrend(trend, vecAltitudes, vecData);
292
	Interpol2D::constant(Interpol1D::arithmeticMean(vecData), dem, grid);
293
	retrend(dem, trend, grid);
294
295
}

296
double IDWAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
297
298
{
	date = i_date;
299
	param = in_param;
300
	nrOfMeasurments = getData(date, param, vecData, vecMeta);
301
302
303
304
305
306
307
308
309
310
311
312

	if (nrOfMeasurments == 0){
		return 0.0;
	} else if (nrOfMeasurments == 1){
		return 0.3;
	} else if (nrOfMeasurments > 1){
		return 0.5;
	}

	return 0.2;
}

313
void IDWAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
314
{
315
316
317
	Interpol2D::IDW(vecData, vecMeta, dem, grid);
}

318
double IDWLapseAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
319
320
{
	date = i_date;
321
	param = in_param;
322
	nrOfMeasurments = getData(date, param, vecData, vecMeta);
323

324
325
326
327
328
329
	if (nrOfMeasurments == 0)
		return 0.0;

	return 0.7;
}

330
void IDWLapseAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
331
{
332
	info.clear(); info.str("");
333
	vector<double> vecAltitudes;
334
	getStationAltitudes(vecMeta, vecAltitudes);
335
	if (vecAltitudes.empty())
336
		throw IOException("Not enough data for spatially interpolating parameter " + MeteoData::getParameterName(param), AT);
337

338
339
340
341
342
	Fit1D trend;
	getTrend(vecAltitudes, vecData, trend);
	info << trend.getInfo();
	detrend(trend, vecAltitudes, vecData);
	Interpol2D::IDW(vecData, vecMeta, dem, grid); //the meta should NOT be used for elevations!
343
	retrend(dem, trend, grid);
344
345
}

346
double LocalIDWLapseAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
347
348
{
	date = i_date;
349
	param = in_param;
350
	nrOfMeasurments = getData(date, param, vecData, vecMeta);
351
352

	if (nrOfMeasurments == 0)
353
354
		return 0.0;

355
356
	return 0.7;
}
357

358
void LocalIDWLapseAlgorithm::calculate(const DEMObject& /*dem*/, Grid2DObject& /*grid*/)
359
{
360
361
	/*info.clear(); info.str("");
	if (nrOfMeasurments == 0)
362
363
364
		throw IOException("Interpolation FAILED for parameter " + MeteoData::getParameterName(param), AT);

	//Set regression coefficients
365
	std::vector<double> vecCoefficients(4, 0.0);
366

367
	size_t nrOfNeighbors=0;
368
369
370
371
372
373
374
375
	//Get the optional arguments for the algorithm: lapse rate, lapse rate usage
	if (vecArgs.size() == 1) { //compute lapse rate on a reduced data set
		IOUtils::convertString(nrOfNeighbors, vecArgs[0]);
	} else { //incorrect arguments, throw an exception
		throw InvalidArgumentException("Wrong number of arguments supplied for the IDW_LAPSE algorithm", AT);
	}

	double r2=0.;
Mathias Bavay's avatar
Mathias Bavay committed
376
	Interpol2D::LocalLapseIDW(vecData, vecMeta, dem, nrOfNeighbors, grid, r2);
377
	info << "r^2=" << Optim::pow2( r2 );*/
378
379
}

380
double RHAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
381
382
383
384
{
	//This algorithm is only valid for RH
	if (in_param != MeteoData::RH)
		return 0.0;
385

386
	date = i_date;
387
	param = in_param;
388
389
	vecData.clear(); vecMeta.clear();
	vecDataTA.clear(); vecDataRH.clear();
390

391
	nrOfMeasurments = 0;
392
	iomanager.getMeteoData(date, vecMeteo);
393
	for (size_t ii=0; ii<vecMeteo.size(); ii++){
394
395
396
		if ((vecMeteo[ii](MeteoData::RH) != IOUtils::nodata) && (vecMeteo[ii](MeteoData::TA) != IOUtils::nodata)){
			vecDataTA.push_back(vecMeteo[ii](MeteoData::TA));
			vecDataRH.push_back(vecMeteo[ii](MeteoData::RH));
397
			vecMeta.push_back(vecMeteo[ii].meta);
398
399
400
401
			nrOfMeasurments++;
		}
	}

402
	if (nrOfMeasurments==0)
403
		return 0.0;
404
	if( (nrOfMeasurments<vecDataRH.size()/2) || ( nrOfMeasurments<2 ) )
405
406
407
408
409
		return 0.6;

	return 0.9;
}

410
void RHAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
411
{
412
	info.clear(); info.str("");
413
	vector<double> vecAltitudes;
414
	getStationAltitudes(vecMeta, vecAltitudes);
415
	if (vecAltitudes.empty())
416
		throw IOException("Not enough data for spatially interpolating parameter " + MeteoData::getParameterName(param), AT);
417
418

	Grid2DObject ta;
419
	mi.interpolate(date, dem, MeteoData::TA, ta); //get TA interpolation from call back to Meteo2DInterpolator
420

421
	//RH->Td, interpolations, Td->RH
422
	//Compute dew point temperatures at stations
423
	std::vector<double> vecTd(vecDataRH.size());
424
	for (size_t ii=0; ii<vecDataRH.size(); ii++){
425
		vecTd[ii] = Atmosphere::RhtoDewPoint(vecDataRH[ii], vecDataTA[ii], 1);
426
	}
427

428
429
430
431
432
	Fit1D trend;
	getTrend(vecAltitudes, vecTd, trend);
	info << trend.getInfo();
	detrend(trend, vecAltitudes, vecTd);
	Interpol2D::IDW(vecTd, vecMeta, dem, grid); //the meta should NOT be used for elevations!
433
	retrend(dem, trend, grid);
434
435

	//Recompute Rh from the interpolated td
436
437
	for (size_t jj=0; jj<grid.nrows; jj++) {
		for (size_t ii=0; ii<grid.ncols; ii++) {
438
439
440
			double &value = grid.grid2D(ii,jj);
			if(value!=IOUtils::nodata)
				value = Atmosphere::DewPointtoRh(value, ta.grid2D(ii,jj), 1);
441
442
443
444
		}
	}
}

445
double ILWRAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
446
447
448
449
450
451
{
	//This algorithm is only valid for ILWR
	if (in_param != MeteoData::ILWR)
		return 0.0;

	date = i_date;
452
	param = in_param;
453
454
	vecData.clear(); vecMeta.clear();
	vecDataEA.clear();
455
456

	nrOfMeasurments = 0;
457
	for (size_t ii=0; ii<vecMeteo.size(); ii++){
458
459
		if ((vecMeteo[ii](MeteoData::ILWR) != IOUtils::nodata) && (vecMeteo[ii](MeteoData::TA) != IOUtils::nodata)){
			vecDataEA.push_back( Atmosphere::blkBody_Emissivity( vecMeteo[ii](MeteoData::ILWR), vecMeteo[ii](MeteoData::TA)) );
460
461
462
463
464
			vecMeta.push_back(vecMeteo[ii].meta);
			nrOfMeasurments++;
		}
	}

465
	if (nrOfMeasurments==0)
466
467
468
469
470
		return 0.0;

	return 0.9;
}

471
void ILWRAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
472
{
473
	info.clear(); info.str("");
474
475
	vector<double> vecAltitudes;
	getStationAltitudes(vecMeta, vecAltitudes);
476
	if (vecAltitudes.empty())
477
		throw IOException("Not enough data for spatially interpolating parameter " + MeteoData::getParameterName(param), AT);
478
479
480
481

	Grid2DObject ta;
	mi.interpolate(date, dem, MeteoData::TA, ta); //get TA interpolation from call back to Meteo2DInterpolator

482
483
484
485
486
	Fit1D trend;
	getTrend(vecAltitudes, vecDataEA, trend);
	info << trend.getInfo();
	detrend(trend, vecAltitudes, vecDataEA);
	Interpol2D::IDW(vecDataEA, vecMeta, dem, grid); //the meta should NOT be used for elevations!
487
	retrend(dem, trend, grid);
488
489

	//Recompute Rh from the interpolated td
490
491
	for (size_t jj=0; jj<grid.nrows; jj++) {
		for (size_t ii=0; ii<grid.ncols; ii++) {
492
493
			double &value = grid.grid2D(ii,jj);
			value = Atmosphere::blkBody_Radiation(value, ta.grid2D(ii,jj));
494
495
496
497
		}
	}
}

498
double SimpleWindInterpolationAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
499
500
501
502
503
504
{
	//This algorithm is only valid for VW
	if (in_param != MeteoData::VW)
		return 0.0;

	date = i_date;
505
	param = in_param;
506
507
	vecData.clear(); vecMeta.clear();
	vecDataVW.clear(); vecDataDW.clear();
508

509
	nrOfMeasurments = 0;
510
	iomanager.getMeteoData(date, vecMeteo);
511
	for (size_t ii=0; ii<vecMeteo.size(); ii++){
512
513
514
		if ((vecMeteo[ii](MeteoData::VW) != IOUtils::nodata) && (vecMeteo[ii](MeteoData::DW) != IOUtils::nodata)){
			vecDataVW.push_back(vecMeteo[ii](MeteoData::VW));
			vecDataDW.push_back(vecMeteo[ii](MeteoData::DW));
515
			vecMeta.push_back(vecMeteo[ii].meta);
516
517
518
			nrOfMeasurments++;
		}
	}
519

520
	if (nrOfMeasurments==0)
521
		return 0.0;
522
	if( (nrOfMeasurments<vecDataVW.size()/2) || ( nrOfMeasurments<2 ) )
523
524
525
526
527
		return 0.6;

	return 0.9;
}

528
void SimpleWindInterpolationAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
529
{
530
	info.clear(); info.str("");
531
	vector<double> vecAltitudes;
532
	getStationAltitudes(vecMeta, vecAltitudes);
533
	if (vecAltitudes.empty())
534
		throw IOException("Not enough data for spatially interpolating parameter " + MeteoData::getParameterName(param), AT);
535
536

	Grid2DObject dw;
537
	mi.interpolate(date, dem, MeteoData::DW, dw); //get DW interpolation from call back to Meteo2DInterpolator
538

539
540
541
542
543
	Fit1D trend;
	getTrend(vecAltitudes, vecDataVW, trend);
	info << trend.getInfo();
	detrend(trend, vecAltitudes, vecDataVW);
	Interpol2D::IDW(vecDataVW, vecMeta, dem, grid); //the meta should NOT be used for elevations!
544
	retrend(dem, trend, grid);
545
546
547
	Interpol2D::SimpleDEMWindInterpolate(dem, grid, dw);
}

548
std::string USERInterpolation::getGridFileName() const
549
{
550
//HACK: use read2DGrid(grid, MeteoGrid::Parameters, Date) instead?
551
552
553
	if (vecArgs.size() != 1){
		throw InvalidArgumentException("Please provide the path to the grids for the USER interpolation algorithm", AT);
	}
554
	const std::string ext(".asc");
555
	const std::string& grid_path = vecArgs[0];
556
	std::string gridname = grid_path + "/";
557

558
	if(!vecMeteo.empty()) {
559
		const Date& timestep = vecMeteo.at(0).date;
560
		gridname =  gridname + timestep.toString(Date::NUM) + "_" + MeteoData::getParameterName(param) + ext;
561
	} else {
562
		gridname = gridname + "Default" + "_" + MeteoData::getParameterName(param) + ext;
563
	}
564

565
566
567
	return gridname;
}

568
double USERInterpolation::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
569
{
570
571
	date = i_date;
	param = in_param;
572
	filename = getGridFileName();
573
574

	if (!IOUtils::validFileName(filename)) {
575
		cerr << "[E] Invalid grid filename for USER interpolation algorithm: " << filename << "\n";
576
577
578
579
580
581
582
583
584
		return 0.0;
	}
	if(IOUtils::fileExists(filename)) {
		return 1.0;
	} else {
		return 0.0;
	}
}

585
void USERInterpolation::calculate(const DEMObject& dem, Grid2DObject& grid)
586
{
587
	iomanager.read2DGrid(grid, filename);
588
	if(!grid.isSameGeolocalization(dem)) {
589
590
		throw InvalidArgumentException("[E] trying to load a grid(" + filename + ") that has not the same georeferencing as the DEM!", AT);
	}
591
592
}

593
double SnowHNWInterpolation::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
594
595
{
	date = i_date;
596
	param = in_param;
597
	nrOfMeasurments = getData(date, param, vecData, vecMeta);
598

599
600
601
	if (nrOfMeasurments == 0)
		return 0.0;

602
	return 0.9;
603
604
}

605
void SnowHNWInterpolation::calculate(const DEMObject& dem, Grid2DObject& grid)
606
607
608
609
{
    info.clear(); info.str("");
    if(internal_dem.isEmpty())
        internal_dem=dem;
610
	info.clear(); info.str("");
611
612
	//retrieve optional arguments
	std::string base_algo;
613
	if (vecArgs.empty()){
614
615
616
617
		base_algo=std::string("IDW_LAPSE");
	} else if (vecArgs.size() == 1){
		IOUtils::convertString(base_algo, vecArgs[0]);
	} else { //incorrect arguments, throw an exception
618
		throw InvalidArgumentException("Wrong number of arguments supplied for the HNW_SNOW algorithm", AT);
619
620
621
622
	}

	//initialize precipitation grid with user supplied algorithm (IDW_LAPSE by default)
	IOUtils::toUpper(base_algo);
623
	vector<string> vecArgs2;
Mathias Bavay's avatar
Mathias Bavay committed
624
	mi.getArgumentsForAlgorithm(MeteoData::getParameterName(param), base_algo, vecArgs2);
625
	auto_ptr<InterpolationAlgorithm> algorithm(AlgorithmFactory::getAlgorithm(base_algo, mi, vecArgs2, iomanager));
626
	algorithm->getQualityRating(date, param);
627
	algorithm->calculate(internal_dem, grid);
628
	info << algorithm->getInfo();
629
630
631
632
	const double orig_mean = grid.grid2D.getMean();

	 //get TA interpolation from call back to Meteo2DInterpolator
	Grid2DObject ta;
633
	mi.interpolate(date, internal_dem, MeteoData::TA, ta);
634
635

	//slope/curvature correction for solid precipitation
636
637
638
639
640
641
642
	Interpol2D::SteepSlopeRedistribution(internal_dem, ta, grid);
    //Interpol2D::CurvatureCorrection(internal_dem, ta, grid);
    //add "virtual snow height" to the internal dem
    for(size_t ii=0; ii<dem.ncols*dem.nrows; ++ii)
        internal_dem.grid2D(ii) += 3.*grid.grid2D(ii);

    iomanager.write2DGrid(internal_dem, MeteoGrids::DEM, date);
643
644
}

645
646
647
648
649
void OrdinaryKrigingAlgorithm::getDataForVariogram(std::vector<double> &distData, std::vector<double> &variData)
{
	distData.clear();
	variData.clear();

650
	for(size_t j=0; j<nrOfMeasurments; j++) {
651
652
653
654
655
		const Coords& st1 = vecMeta[j].position;
		const double x1 = st1.getEasting();
		const double y1 = st1.getNorthing();
		const double val1 = vecData[j];

656
		for(size_t i=0; i<j; i++) {
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
			//compute distance between stations
			const Coords& st2 = vecMeta[i].position;
			const double val2 = vecData[i];
			const double DX = x1-st2.getEasting();
			const double DY = y1-st2.getNorthing();
			const double distance = Optim::fastSqrt_Q3(Optim::pow2(DX) + Optim::pow2(DY));

			distData.push_back(distance);
			variData.push_back( 0.5*Optim::pow2(val1-val2) );
		}
	}
}

bool OrdinaryKrigingAlgorithm::computeVariogram()
{//return variogram fit of covariance between stations i and j
	std::vector<double> distData, variData;
	getDataForVariogram(distData, variData);

	std::vector<string> vario_types( vecArgs );
	if(vario_types.empty()) vario_types.push_back("LINVARIO");
	/*std::vector<string> vario_types;
	vario_types.push_back("SPHERICVARIO");
	vario_types.push_back("EXPVARIO");
	vario_types.push_back("LINVARIO");*/

	size_t args_index=0;
	do {
		const string vario_model = IOUtils::strToUpper( vario_types[args_index] );
		const bool status = variogram.setModel(vario_model, distData, variData);
		if(status) {
			info << " - " << vario_model;
			return true;
		}

		args_index++;
	} while (args_index<vario_types.size());

	return false;
}

697
double OrdinaryKrigingAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
698
699
{
	date = i_date;
700
	param = in_param;
701
	nrOfMeasurments = getData(date, param, vecData, vecMeta);
702

703
	if(nrOfMeasurments==0) return 0.;
704
	if(nrOfMeasurments>=7) return 0.9;
705
	return 0.1;
706
707
}

708
void OrdinaryKrigingAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
709
{
710
	info.clear(); info.str("");
711
712
713
714
715
716
717
	//optimization: getRange (from variogram fit -> exclude stations that are at distances > range (-> smaller matrix)
	//or, get max range from io.ini, build variogram from this user defined max range
	if(!computeVariogram()) //only refresh once a month, or once a week, etc
		throw IOException("The variogram for parameter " + MeteoData::getParameterName(param) + " could not be computed!", AT);
	Interpol2D::ODKriging(vecData, vecMeta, dem, variogram, grid);
}

718
void LapseOrdinaryKrigingAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
719
{
720
	info.clear(); info.str("");
721
722
	//optimization: getRange (from variogram fit -> exclude stations that are at distances > range (-> smaller matrix)
	//or, get max range from io.ini, build variogram from this user defined max range
723
724
	vector<double> vecAltitudes;
	getStationAltitudes(vecMeta, vecAltitudes);
725
	if (vecAltitudes.empty())
726
727
		throw IOException("Not enough data for spatially interpolating parameter " + MeteoData::getParameterName(param), AT);

728
729
730
	Fit1D trend(Fit1D::NOISY_LINEAR, vecAltitudes, vecData, false);
	if(!trend.fit())
		throw IOException("Interpolation FAILED for parameter " + MeteoData::getParameterName(param) + ": " + trend.getInfo(), AT);
731
732
733
	info << trend.getInfo();
	detrend(trend, vecAltitudes, vecData);

734
735
	if(!computeVariogram()) //only refresh once a month, or once a week, etc
		throw IOException("The variogram for parameter " + MeteoData::getParameterName(param) + " could not be computed!", AT);
736
	Interpol2D::ODKriging(vecData, vecMeta, dem, variogram, grid);
737

738
	retrend(dem, trend, grid);
739
}
740

741
} //namespace