WSL/SLF GitLab Repository

Meteo2DInterpolator.cc 15.2 KB
Newer Older
1
/***********************************************************************************/
2
/*  Copyright 2014 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
/***********************************************************************************/
/* 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/Meteo2DInterpolator.h>

using namespace std;

namespace mio {
24

25
26
Meteo2DInterpolator::Meteo2DInterpolator(const Config& i_cfg, TimeSeriesManager& i_tsmanager, GridsManager& i_gridsmanager)
                    : cfg(i_cfg), tsmanager(i_tsmanager), gridsmanager(i_gridsmanager),
27
                      grid_buffer(0), mapAlgorithms(),
28
                      v_params(), v_coords(), v_stations(), virtual_point_cache(),
29
                      algorithms_ready(false), use_full_dem(false), downscaling(false), virtual_stations(false)
30
{
31
32
33
34
	size_t max_grids = 10; //default number of grids to keep in buffer
	cfg.getValue("BUFF_GRIDS", "Interpolations2D", max_grids, IOUtils::nothrow);
	grid_buffer.setMaxGrids(max_grids);
	
35
	setAlgorithms();
36
37
38
39
	cfg.getValue("Virtual_stations", "Input", virtual_stations, IOUtils::nothrow);
	if (virtual_stations) {
		initVirtualStations();
	}
40
41
}

42
Meteo2DInterpolator::~Meteo2DInterpolator()
43
{
44
45
46
	std::map<std::string, std::vector<InterpolationAlgorithm*> >::iterator iter;
	for (iter = mapAlgorithms.begin(); iter != mapAlgorithms.end(); ++iter) {
		const vector<InterpolationAlgorithm*>& vecAlgs = iter->second;
47
		for (size_t ii=0; ii<vecAlgs.size(); ++ii)
48
49
			delete vecAlgs[ii];
	}
50
51
}

52
53
54
55
56
/* By reading the Config object build up a list of user configured algorithms
* for each MeteoData::Parameters parameter (i.e. each member variable of MeteoData like ta, p, hnw, ...)
* Concept of this constructor: loop over all MeteoData::Parameters and then look
* for configuration of interpolation algorithms within the Config object.
*/
57
58
void Meteo2DInterpolator::setAlgorithms()
{
59
//HACK set callback to internal iomanager for virtual stations and downsampling!
60
	set<string> set_of_used_parameters;
61
62
63
64
65
	get_parameters(cfg, set_of_used_parameters);

	set<string>::const_iterator it;
	for (it = set_of_used_parameters.begin(); it != set_of_used_parameters.end(); ++it) {
		const std::string parname = *it;
66
		std::vector<std::string> tmpAlgorithms;
67
68
69
		const size_t nrOfAlgorithms = getAlgorithmsForParameter(cfg, parname, tmpAlgorithms);

		std::vector<InterpolationAlgorithm*> vecAlgorithms(nrOfAlgorithms);
70
		for (size_t jj=0; jj<nrOfAlgorithms; jj++) {
71
72
			std::vector<std::string> vecArgs;
			getArgumentsForAlgorithm(parname, tmpAlgorithms[jj], vecArgs);
73
			vecAlgorithms[jj] = AlgorithmFactory::getAlgorithm( tmpAlgorithms[jj], *this, vecArgs, tsmanager, gridsmanager);
74
75
		}

76
		if (nrOfAlgorithms>0) {
77
78
			mapAlgorithms[parname] = vecAlgorithms;
		}
79
80
	}
	algorithms_ready = true;
81
}
82

83
//get a list of all meteoparameters referenced in the Interpolations2D section
84
size_t Meteo2DInterpolator::get_parameters(const Config& cfg, std::set<std::string>& set_parameters)
85
{
86
87
88
89
90
91
92
93
	std::vector<std::string> vec_keys;
	cfg.findKeys(vec_keys, std::string(), "Interpolations2D");

	for (size_t ii=0; ii<vec_keys.size(); ii++) {
		const size_t found = vec_keys[ii].find_first_of(":");
		if (found != std::string::npos){
			const string tmp = vec_keys[ii].substr(0,found);
			set_parameters.insert( IOUtils::strToUpper(tmp) );
94
95
		}
	}
96
97

	return set_parameters.size();
98
99
}

100
101
void Meteo2DInterpolator::interpolate(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
                                      Grid2DObject& result)
102
103
{
	std::string InfoString;
104
	interpolate(date, dem, meteoparam, result, InfoString);
105
106
}

107
108
void Meteo2DInterpolator::interpolate(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
                                      Grid2DObject& result, std::string& InfoString)
109
{
110
	if (!algorithms_ready)
111
		setAlgorithms();
112

113
114
	const string param_name = MeteoData::getParameterName(meteoparam);
	
115
	//Get grid from buffer if it exists
116
117
118
	std::ostringstream grid_hash;
	grid_hash << dem.llcorner.printLatLon() << " " << dem.getNx() << "x" << dem.getNy() << " @" << dem.cellsize << " " << date.toString(Date::ISO) << " " << param_name;
	if (grid_buffer.get(result, grid_hash.str(), InfoString))
119
120
		return;

121
	//Show algorithms to be used for this parameter
122
	const map<string, vector<InterpolationAlgorithm*> >::iterator it = mapAlgorithms.find(param_name);
123
	if (it==mapAlgorithms.end()) {
124
125
		throw IOException("No interpolation algorithms configured for parameter "+param_name, AT);
	}
126

127
128
129
130
131
132
133
134
135
136
	//look for algorithm with the highest quality rating
	const vector<InterpolationAlgorithm*>& vecAlgs = it->second;
	double maxQualityRating = -1.;
	size_t bestalgorithm = 0;
	for (size_t ii=0; ii < vecAlgs.size(); ++ii){
		const double rating = vecAlgs[ii]->getQualityRating(date, meteoparam);
		if ((rating != 0.0) && (rating > maxQualityRating)) {
			//we use ">" so that in case of equality, the first choice will be kept
			bestalgorithm = ii;
			maxQualityRating = rating;
137
138
139
		}
	}

140
	//finally execute the algorithm with the best quality rating or throw an exception
141
	if (maxQualityRating<=0.0)
142
		throw IOException("No interpolation algorithm with quality rating >0 found for parameter "+param_name+" on "+date.toString(Date::ISO_TZ), AT);
143
144
145
	vecAlgs[bestalgorithm]->calculate(dem, result);
	InfoString = vecAlgs[bestalgorithm]->getInfo();

146
	//Run soft min/max filter for RH, HNW and HS
147
148
149
150
	if (meteoparam == MeteoData::RH){
		Meteo2DInterpolator::checkMinMax(0.0, 1.0, result);
	} else if (meteoparam == MeteoData::HNW){
		Meteo2DInterpolator::checkMinMax(0.0, 10000.0, result);
151
152
	} else if (meteoparam == MeteoData::HS){
		Meteo2DInterpolator::checkMinMax(0.0, 10000.0, result);
153
154
	} else if (meteoparam == MeteoData::VW){
		Meteo2DInterpolator::checkMinMax(0.0, 10000.0, result);
155
	}
156

157
158
	//save grid in buffer
	grid_buffer.push(result, grid_hash.str(), InfoString);
159
160
161
162
}

//HACK make sure that skip_virtual_stations = true before calling this method when using virtual stations!
void Meteo2DInterpolator::interpolate(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
163
                            const std::vector<Coords>& in_coords, std::vector<double>& result, std::string& info_string)
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
196
197
198
199
200
201
202
203
204
205
206
{
	result.clear();
	vector<Coords> vec_coords(in_coords);

	if (use_full_dem) {
		Grid2DObject result_grid;
		interpolate(date, dem, meteoparam, result_grid, info_string);
		const bool gridify_success = dem.gridify(vec_coords);
		if (!gridify_success)
			throw InvalidArgumentException("Coordinate given to interpolate is outside of dem", AT);

		for (size_t ii=0; ii<vec_coords.size(); ii++) {
			//we know the i,j are positive because of gridify_success
			const size_t pt_i = static_cast<size_t>( vec_coords[ii].getGridI() );
			const size_t pt_j = static_cast<size_t>( vec_coords[ii].getGridJ() );
			result.push_back( result_grid(pt_i,pt_j) );
		}
	} else {
		for (size_t ii=0; ii<vec_coords.size(); ii++) {
			const bool gridify_success = dem.gridify(vec_coords[ii]);
			if (!gridify_success)
				throw InvalidArgumentException("Coordinate given to interpolate is outside of dem", AT);

			//we know the i,j are positive because of gridify_success
			const size_t pt_i = static_cast<size_t>( vec_coords[ii].getGridI() );
			const size_t pt_j = static_cast<size_t>( vec_coords[ii].getGridJ() );

			//Make new DEM with just one point, namely the one specified by vec_coord[ii]
			//Copy all other properties of the big DEM into the new one
			DEMObject one_point_dem(dem, pt_i, pt_j, 1, 1, false);

			one_point_dem.min_altitude = dem.min_altitude;
			one_point_dem.max_altitude = dem.max_altitude;
			one_point_dem.min_slope = dem.min_slope;
			one_point_dem.max_slope = dem.max_slope;
			one_point_dem.min_curvature = dem.min_curvature;
			one_point_dem.max_curvature = dem.max_curvature;

			Grid2DObject result_grid;
			interpolate(date, one_point_dem, meteoparam, result_grid, info_string);
			result.push_back(result_grid(0,0));
		}
	}
207
208
}

209
size_t Meteo2DInterpolator::getAlgorithmsForParameter(const Config& cfg, const std::string& parname, std::vector<std::string>& vecAlgorithms)
210
{
211
212
	// This function retrieves the user defined interpolation algorithms for
	// parameter 'parname' by querying the Config object
213
	vecAlgorithms.clear();
214
	std::vector<std::string> vecKeys;
215
216
217
218
219
	cfg.findKeys(vecKeys, parname+"::algorithms", "Interpolations2D");

	if (vecKeys.size() > 1)
		throw IOException("Multiple definitions of " + parname + "::algorithms in config file", AT);;

220
	if (vecKeys.empty())
221
222
		return 0;

223
	cfg.getValue(vecKeys[0], "Interpolations2D", vecAlgorithms, IOUtils::nothrow);
224
225
226
	return vecAlgorithms.size();
}

227
size_t Meteo2DInterpolator::getArgumentsForAlgorithm(const std::string& param,
228
229
                                                     const std::string& algorithm,
                                                     std::vector<std::string>& vecArgs) const
230
231
{
	vecArgs.clear();
232
	const string keyname = param +"::"+ algorithm;
233
	cfg.getValue(keyname, "Interpolations2D", vecArgs, IOUtils::nothrow);
234
235
236
237

	return vecArgs.size();
}

238
239
void Meteo2DInterpolator::checkMinMax(const double& minval, const double& maxval, Grid2DObject& gridobj)
{
Mathias Bavay's avatar
Mathias Bavay committed
240
	const size_t nxy = gridobj.getNx() * gridobj.getNy();
241

Mathias Bavay's avatar
Mathias Bavay committed
242
	for (size_t ii=0; ii<nxy; ii++){
243
		double& value = gridobj(ii);
244
245
246
247
248
249
250
		if (value == IOUtils::nodata){
			continue;
		}
		if (value < minval) {
			value = minval;
		} else if (value > maxval) {
			value = maxval;
251
252
253
254
		}
	}
}

255
256
257
258
259
void Meteo2DInterpolator::check_projections(const DEMObject& dem, const std::vector<MeteoData>& vec_meteo)
{
	//check that the stations are using the same projection as the dem
	for (size_t ii=0; ii<vec_meteo.size(); ii++) {
		const StationData& meta = vec_meteo[ii].meta;
260
		if (!meta.position.isSameProj(dem.llcorner)) {
261
			std::ostringstream os;
262
263
264
265
266
267
268
269
270
271
			std::string type, args;
			meta.position.getProj(type, args);
			os << "Station " << meta.stationID << " is using projection (" << type << " " << args << ") ";
			dem.llcorner.getProj(type, args);
			os << "while DEM is using projection ("<< type << " " << args << ") ";
			throw IOException(os.str(), AT);
		}
	}
}

272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
//get the stations' data to use for downscaling (=true measurements)
size_t Meteo2DInterpolator::getVirtualMeteoData(const vstations_policy& strategy, const Date& i_date, METEO_SET& vecMeteo)
{
	if (strategy==VSTATIONS) {
		return getVirtualStationsData(i_date, vecMeteo);
	} else if (strategy==DOWNSCALING) {
		//extract all grid points
		return 0; //hack
	} else if (strategy==SMART_DOWNSCALING) {
		//for each parameter:
		//call iomanager.read2DGrid()
		//extract relevant points and fill vecMeteo
		//and loop over all grids!
		//
		//Questions/tricks:
		//   throw exception if grids don't match between parameters
		//   should we recompute which points to take between each time steps? (ie more flexibility)
		//   should the generated virtual stations data be filtered? (useful for applying corrections)
		//          if so: have an own iomanager and use iomanager.push_meteo_data()
		return 0; //hack
	}

	throw UnknownValueException("Unknown virtual station strategy", AT);
}

void Meteo2DInterpolator::initVirtualStations()
{
	if(!cfg.keyExists("DEM", "Input"))
		throw NoAvailableDataException("In order to use virtual stations, please provide a DEM!", AT);
	DEMObject dem;
	gridsmanager.readDEM(dem);

	//get virtual stations coordinates
	std::string coordin, coordinparam, coordout, coordoutparam;
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);

	std::vector<std::string> vecStation;
	cfg.getValues("Vstation", "INPUT", vecStation);
	for(size_t ii=0; ii<vecStation.size(); ii++) {
311
		//The coordinate specification is given as either: "easting northing epsg" or "lat lon"
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
		Coords tmp(coordin, coordinparam, vecStation[ii]);
		if(!tmp.isNodata())
			v_coords.push_back( tmp );
	}

	//create stations' metadata
	for(size_t ii=0; ii<v_coords.size(); ii++) {
		if(!dem.gridify(v_coords[ii])) {
			ostringstream ss;
			ss << "Virtual station \"" << vecStation[ii] << "\" is not contained is provided DEM";
			throw NoAvailableDataException(ss.str(), AT);
		}

		const size_t i = v_coords[ii].getGridI(), j = v_coords[ii].getGridJ();
		v_coords[ii].setAltitude(dem(i,j), false);

		ostringstream name;
		name << "Virtual_Station_" << ii+1;
		ostringstream id;
		id << "VIR" << ii+1;
		StationData sd(v_coords[ii], id.str(), name.str());
		sd.setSlope(dem.slope(i,j), dem.azi(i,j));

		v_stations.push_back( sd );
	}

	cfg.getValue("Interpol_Use_Full_DEM", "Input", use_full_dem, IOUtils::nothrow);
}

size_t Meteo2DInterpolator::getVirtualStationsData(const Date& i_date, METEO_SET& vecMeteo)
{ //HACK use own private iomanager and do caching in its private one
	vecMeteo.clear();

	// Check if data is available in cache
	const map<Date, vector<MeteoData> >::const_iterator it = virtual_point_cache.find(i_date);
	if (it != virtual_point_cache.end()){
		vecMeteo = it->second;
		return vecMeteo.size();
	}

	//get data from real input stations
	METEO_SET vecTrueMeteo;
354
	tsmanager.getMeteoData(i_date, vecTrueMeteo);
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
	if (vecTrueMeteo.empty()) return 0;

	if (v_params.empty()) {
		//get parameters to interpolate if not already done
		//we need valid data in order to handle extra parameters
		std::vector<std::string> vecStr;
		cfg.getValue("Virtual_parameters", "Input", vecStr);
		for (size_t ii=0; ii<vecStr.size(); ii++) {
			v_params.push_back( vecTrueMeteo[0].getParameterIndex(vecStr[ii]) );
		}
	}

	//create stations without measurements
	for (size_t ii=0; ii<v_stations.size(); ii++) {
		MeteoData md(i_date, v_stations[ii]);
		vecMeteo.push_back( md );
	}

	//fill meteo parameters
	DEMObject dem;
	gridsmanager.readDEM(dem); //this is not a big deal since it will be in the buffer
	string info_string;
	for (size_t param=0; param<v_params.size(); param++) {
		std::vector<double> result;
		interpolate(i_date, dem, static_cast<MeteoData::Parameters>(v_params[param]), v_coords, result, info_string);
		for (size_t ii=0; ii<v_coords.size(); ii++)
			vecMeteo[ii](v_params[param]) = result[ii];
	}

	return vecMeteo.size();
}

387

388
const std::string Meteo2DInterpolator::toString() const {
389
	ostringstream os;
390
	os << "<Meteo2DInterpolator>\n";
391
	os << "Config& cfg = " << hex << &cfg << dec << "\n";
392
393
	os << "TimeSeriesManager& tsmanager = "  << hex << &tsmanager << dec << "\n";
	os << "GridsManager& gridsmanager = "  << hex << &gridsmanager << dec << "\n";
394

395
	os << "Spatial resampling algorithms:\n";
396
	std::map<std::string, std::vector<InterpolationAlgorithm*> >::const_iterator iter;
397
	for (iter = mapAlgorithms.begin(); iter != mapAlgorithms.end(); ++iter) {
398
		os << setw(10) << iter->first << "::";
399
		for (size_t jj=0; jj<iter->second.size(); jj++) {
400
			os << iter->second[jj]->algo << " ";
401
402
403
404
		}
		os << "\n";
	}

405
	//cache content
406
	os << grid_buffer.toString();
407
	os << "</Meteo2DInterpolator>\n";
408
	return os.str();
409
410
}

411
} //namespace