WSL/SLF GitLab Repository

Meteo2DInterpolator.cc 22.1 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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
/**
 * @page virtual_stations Virtual stations handling
 * It is possible to use spatially interpolated meteorological fields or time series of 2D grids to extract meteorological time series for a set of points.
 * This is handled as "virtual stations" since the data **will seem to originate from these virtual stations points** where no station is present. 
 * 
 * *Virtual stations are currently in a "beta" stage and their implementation might significantly change in the future.*
 *
 * @section virtual_stations_from_interpolation From spatial interpolations
 * The data from real input stations (as read by the plugin defined with the METEO key in the [input] section) is filtered/processed, temporally interpolated and 
 * spatially interpolated as defined in the configuration file. Then time series are reconstructed from these grids at a set of defined points. This behavior is configured
 * by the following keys (in the [Input] section):
 *    + VIRTUAL_STATIONS set to *true*;
 *    + VSTATION# : provide the lat, lon and (optionally) the epsg code for a virtual station;
 *    + VIRTUAL_PARAMETERS: list of MeteoData::Parameters that have to be interpolated to populate the virtual stations;
 *    + VSTATIONS_REFRESH_RATE: how often to rebuild the spatial interpolations, in seconds (otherwise, only IOUtils::nodata will be returned);
 *    + INTERPOL_USE_FULL_DEM: should the spatial interpolations be performed on the whole DEM? (this is necessary for some algorithms, for example WINSTAL);
 * 
 * Currently, a DEM also has to be provided since this will be used to retrieve the elevation, slope and azimuth of the virtual stations.
 * 
 * @code
 * DEM = ARC
 * DEMFILE = ./input/surface-grids/davos.asc
 *
 * #here, the real data as measured by some stations 
 * METEO		= IMIS
 * DBNAME		= sdbo
 * DBUSER		= xxx
 * DBPASS		= xxx
 * STATION1	= *WFJ
 * STATION2	= STB2
 * STATION3	= WFJ2
 * STATION4	= *DAV
 * 
 * #here the locations where the data will be generated. The caller will only see these stations!
 * Virtual_stations = true
 * VSTATION1 = 46.793029 9.821343
 * VSTATION2 = 46.793031 9.831572
 * Virtual_parameters = TA RH PSUM ILWR P VW RSWR
 * VSTATIONS_REFRESH_RATE = 21600
 * @endcode
 * 
* @section virtual_stations_from_grids From gridded data
* The meteorological time series are extracted from time series of user-provided grids. therefore a plugin for 2D grids must have been defined (with the GRID2D key in
* the [Input] section). The following keys control this downscaling process:
*    + DOWNSCALING set to *true*;
*    + VSTATION# : provide the lat, lon and (optionally) the epsg code for a virtual station;
*    + VIRTUAL_PARAMETERS: list of MeteoData::Parameters that have to be interpolated to populate the virtual stations; 
* 
* Currently, a DEM has to be provided in order to check the position of the stations and the consistency of the grids.
 * 
 */

77
78
Meteo2DInterpolator::Meteo2DInterpolator(const Config& i_cfg, TimeSeriesManager& i_tsmanager, GridsManager& i_gridsmanager)
                    : cfg(i_cfg), tsmanager(i_tsmanager), gridsmanager(i_gridsmanager),
79
                      grid_buffer(0), mapAlgorithms(),
80
                      v_params(), v_coords(), v_stations(), virtual_point_cache(),
81
                      vstations_refresh_rate(IOUtils::unodata), algorithms_ready(false), use_full_dem(false)
82
{
83
84
85
86
	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);
	
87
	setAlgorithms();
88
	bool virtual_stations = false; ///< compute the meteo values at virtual stations
89
	cfg.getValue("Virtual_stations", "Input", virtual_stations, IOUtils::nothrow);
90
91
	if (virtual_stations) {
		cfg.getValue("VSTATIONS_REFRESH_RATE", "Input", vstations_refresh_rate, IOUtils::nothrow);
92
		cfg.getValue("Interpol_Use_Full_DEM", "Input", use_full_dem, IOUtils::nothrow);
93
94
	}
	
95
	bool downscaling = false; ///< Are we downscaling meteo grids instead of interpolating stations' data?
96
	cfg.getValue("Downscaling", "Input", downscaling, IOUtils::nothrow);
97
98
99
	if (virtual_stations && downscaling)
		throw InvalidArgumentException("It is not possible to use both Virtual_stations and Downscaling!", AT);
		
100
	if (virtual_stations || downscaling) {
101
		initVirtualStations(downscaling); //adjust the coordinates if downscaling
102
	}
103
104
}

105
Meteo2DInterpolator::~Meteo2DInterpolator()
106
{
107
108
109
	std::map<std::string, std::vector<InterpolationAlgorithm*> >::iterator iter;
	for (iter = mapAlgorithms.begin(); iter != mapAlgorithms.end(); ++iter) {
		const vector<InterpolationAlgorithm*>& vecAlgs = iter->second;
110
		for (size_t ii=0; ii<vecAlgs.size(); ++ii)
111
112
			delete vecAlgs[ii];
	}
113
114
}

115
/* By reading the Config object build up a list of user configured algorithms
116
* for each MeteoData::Parameters parameter (i.e. each member variable of MeteoData like ta, p, psum, ...)
117
118
119
* Concept of this constructor: loop over all MeteoData::Parameters and then look
* for configuration of interpolation algorithms within the Config object.
*/
120
121
void Meteo2DInterpolator::setAlgorithms()
{
122
	std::set<std::string> set_of_used_parameters;
123
124
	get_parameters(cfg, set_of_used_parameters);

125
	std::set<std::string>::const_iterator it;
126
	for (it = set_of_used_parameters.begin(); it != set_of_used_parameters.end(); ++it) {
127
		const std::string parname( *it );
128
		std::vector<std::string> tmpAlgorithms;
129
130
131
		const size_t nrOfAlgorithms = getAlgorithmsForParameter(cfg, parname, tmpAlgorithms);

		std::vector<InterpolationAlgorithm*> vecAlgorithms(nrOfAlgorithms);
132
		for (size_t jj=0; jj<nrOfAlgorithms; jj++) {
133
134
			std::vector<std::string> vecArgs;
			getArgumentsForAlgorithm(parname, tmpAlgorithms[jj], vecArgs);
135
			vecAlgorithms[jj] = AlgorithmFactory::getAlgorithm( tmpAlgorithms[jj], *this, vecArgs, tsmanager, gridsmanager);
136
137
		}

138
		if (nrOfAlgorithms>0) {
139
140
			mapAlgorithms[parname] = vecAlgorithms;
		}
141
142
	}
	algorithms_ready = true;
143
}
144

145
//get a list of all meteoparameters referenced in the Interpolations2D section
146
size_t Meteo2DInterpolator::get_parameters(const Config& cfg, std::set<std::string>& set_parameters)
147
{
148
	std::vector<std::string> vec_keys;
149
	cfg.findKeys(vec_keys, "::algorithms", "Interpolations2D", true);
150
151
152
153

	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){
154
			const std::string tmp( vec_keys[ii].substr(0,found) );
155
			set_parameters.insert( IOUtils::strToUpper(tmp) );
156
157
		}
	}
158
159

	return set_parameters.size();
160
161
}

162
163
void Meteo2DInterpolator::interpolate(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
                                      Grid2DObject& result)
164
165
{
	std::string InfoString;
166
	interpolate(date, dem, meteoparam, result, InfoString);
167
168
}

169
170
void Meteo2DInterpolator::interpolate(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
                                      Grid2DObject& result, std::string& InfoString)
171
{
172
	if (!algorithms_ready)
173
		setAlgorithms();
174

175
	const std::string param_name( MeteoData::getParameterName(meteoparam) );
176
	
177
	//Get grid from buffer if it exists
178
	std::ostringstream grid_hash;
179
	grid_hash << dem.llcorner.toString(Coords::LATLON) << " " << dem.getNx() << "x" << dem.getNy() << " @" << dem.cellsize << " " << date.toString(Date::ISO) << " " << param_name;
180
	if (grid_buffer.get(result, grid_hash.str(), InfoString))
181
182
		return;

183
	//Show algorithms to be used for this parameter
184
	const std::map<string, vector<InterpolationAlgorithm*> >::iterator it = mapAlgorithms.find(param_name);
185
	if (it==mapAlgorithms.end())
186
		throw IOException("No interpolation algorithms configured for parameter "+param_name, AT);
187

188
189
190
191
192
193
194
195
196
197
	//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;
198
199
200
		}
	}

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

207
	//Run soft min/max filter for RH, PSUM and HS
208
209
	if (meteoparam == MeteoData::RH){
		Meteo2DInterpolator::checkMinMax(0.0, 1.0, result);
210
	} else if (meteoparam == MeteoData::PSUM){
211
		Meteo2DInterpolator::checkMinMax(0.0, 10000.0, result);
212
213
	} else if (meteoparam == MeteoData::HS){
		Meteo2DInterpolator::checkMinMax(0.0, 10000.0, result);
214
215
	} else if (meteoparam == MeteoData::VW){
		Meteo2DInterpolator::checkMinMax(0.0, 10000.0, result);
216
	}
217

218
219
	//save grid in buffer
	grid_buffer.push(result, grid_hash.str(), InfoString);
220
221
222
223
}

//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,
224
                            const std::vector<Coords>& in_coords, std::vector<double>& result, std::string& info_string)
225
226
{
	result.clear();
227
	std::vector<Coords> vec_coords(in_coords);
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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

	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));
		}
	}
268
269
}

270
size_t Meteo2DInterpolator::getAlgorithmsForParameter(const Config& cfg, const std::string& parname, std::vector<std::string>& vecAlgorithms)
271
{
272
273
	// This function retrieves the user defined interpolation algorithms for
	// parameter 'parname' by querying the Config object
274
	vecAlgorithms.clear();
275
	std::vector<std::string> vecKeys;
276
277
278
279
280
	cfg.findKeys(vecKeys, parname+"::algorithms", "Interpolations2D");

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

281
	if (vecKeys.empty()) return 0;
282

283
	cfg.getValue(vecKeys[0], "Interpolations2D", vecAlgorithms, IOUtils::nothrow);
284
285
286
	return vecAlgorithms.size();
}

287
size_t Meteo2DInterpolator::getArgumentsForAlgorithm(const std::string& param,
288
289
                                                     const std::string& algorithm,
                                                     std::vector<std::string>& vecArgs) const
290
291
{
	vecArgs.clear();
292
	const std::string keyname = param +"::"+ algorithm;
293
	cfg.getValue(keyname, "Interpolations2D", vecArgs, IOUtils::nothrow);
294
295
296
297

	return vecArgs.size();
}

298
299
void Meteo2DInterpolator::checkMinMax(const double& minval, const double& maxval, Grid2DObject& gridobj)
{
300
	for (size_t ii=0; ii<gridobj.size(); ii++){
301
		double& value = gridobj(ii);
302
		if (value == IOUtils::nodata)
303
			continue;
304

305
306
307
308
		if (value < minval) {
			value = minval;
		} else if (value > maxval) {
			value = maxval;
309
310
311
312
		}
	}
}

313
314
315
316
317
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;
318
		if (!meta.position.isSameProj(dem.llcorner)) {
319
			std::ostringstream os;
320
321
322
323
324
325
326
327
328
329
			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);
		}
	}
}

330
331
332
333
//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) {
334
		//this reads station data, interpolates the stations and extract points from the interpolated grids
335
336
		return getVirtualStationsData(i_date, vecMeteo);
	} else if (strategy==SMART_DOWNSCALING) {
337
338
		//This reads already gridded data and extract points from the grids
		
339
340
		//Questions/tricks:
		//   should we recompute which points to take between each time steps? (ie more flexibility)
341
		//   move this to the IOHandler so the data could be filtered, corrected, interpolated
342
343
344
		return getVirtualStationsFromGrid(i_date, vecMeteo);
	} else if (strategy==DOWNSCALING) {
		//extract all grid points
345
346
347
348
349
350
		return 0; //hack
	}

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

351
352
353
354
/** @brief read the list of virtual stations
 * @param adjust_coordinates should the coordinates be recomputed to match DEM cells?
 */
void Meteo2DInterpolator::initVirtualStations(const bool& adjust_coordinates)
355
{
356
	if (!cfg.keyExists("DEM", "Input"))
357
		throw NoDataException("In order to use virtual stations, please provide a DEM!", AT);
358
	DEMObject dem;
359
	dem.setUpdatePpt( DEMObject::SLOPE ); //we only need the elevation
360
361
362
363
364
	gridsmanager.readDEM(dem);

	//get virtual stations coordinates
	std::string coordin, coordinparam, coordout, coordoutparam;
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
365
	dem.llcorner.setProj(coordin, coordinparam); //make sure the DEM and the VStations are in the same projection
366
367
	const double dem_easting = dem.llcorner.getEasting();
	const double dem_northing = dem.llcorner.getNorthing();
368

369
370
	std::vector<std::string> vecStation, vecKeys;
	cfg.getValues("Vstation", "INPUT", vecStation, vecKeys);
371
	for (size_t ii=0; ii<vecStation.size(); ii++) {
372
		//The coordinate specification is given as either: "easting northing epsg" or "lat lon"
373
		Coords curr_point(coordin, coordinparam, vecStation[ii]);
374

375
376
377
378
379
380
381
382
		if (!curr_point.isNodata()) {
			v_coords.push_back( curr_point ); //so we can check for duplicates

			if (!dem.gridify(curr_point)) {
				ostringstream ss;
				ss << "Virtual station \"" << vecStation[ii] << "\" is not contained is provided DEM";
				throw NoDataException(ss.str(), AT);
			}
383

384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
			const size_t i = curr_point.getGridI(), j = curr_point.getGridJ();
			if (adjust_coordinates) { //adjust coordinates to match the chosen cell
				const double easting = dem_easting + dem.cellsize*static_cast<double>(i);
				const double northing = dem_northing + dem.cellsize*static_cast<double>(j);
				curr_point.setXY(easting, northing, dem(i,j));
				curr_point.setGridIndex(static_cast<int>(i), static_cast<int>(j), IOUtils::inodata, true);
			} else {
				curr_point.setAltitude(dem(i,j), false);
			}

			//remove duplicate stations, ie stations that have same easting,northing,altitude
			bool is_duplicate = false;
			for (size_t jj=0; jj<ii; jj++) {
				if (curr_point==v_coords[jj]) {
					std::cout << "[W] removing VSTATION" << ii+1 << " as a duplicate of VSTATION" << jj+1 << "\n";
					is_duplicate = true;
					break;
				}
			}
			if (!is_duplicate) {
				//extract vstation number, build the station name and station ID
				const std::string id_num = vecKeys[ii].substr(string("Vstation").length());
				StationData sd(curr_point, "VIR"+id_num, "Virtual_Station_"+id_num);
				sd.setSlope(dem.slope(i,j), dem.azi(i,j));
				v_stations.push_back( sd );
409
410
			}
		}
411
412
413
414
	}
}

size_t Meteo2DInterpolator::getVirtualStationsData(const Date& i_date, METEO_SET& vecMeteo)
415
{
416
417
418
	vecMeteo.clear();

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

	//get data from real input stations
	METEO_SET vecTrueMeteo;
427
	tsmanager.getMeteoData(i_date, vecTrueMeteo);
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
	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 );
	}
445
	if ( vstations_refresh_rate!=IOUtils::unodata && Date::mod(i_date, vstations_refresh_rate) != 0) return vecMeteo.size();
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460

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

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
size_t Meteo2DInterpolator::getVirtualStationsFromGrid(const Date& i_date, METEO_SET& vecMeteo)
{
	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();
	}
	
	if (v_params.empty()) { //get parameters to interpolate if not already done
		std::vector<std::string> vecStr;
		cfg.getValue("Virtual_parameters", "Input", vecStr);		
		for (size_t ii=0; ii<vecStr.size(); ii++) {
			const size_t param_idx = MeteoGrids::getParameterIndex( vecStr[ii] );
			if (param_idx==IOUtils::npos)
				throw InvalidArgumentException("Invalid parameter '" + vecStr[ii] + "', only standard parameters can be extracted from grids for virtual stations! ", AT);
			v_params.push_back( param_idx );
		}
	}

	//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 );
	}
	
489
490
491
492
	DEMObject dem;
	dem.setUpdatePpt( DEMObject::NO_UPDATE ); //we only need the elevation
	gridsmanager.readDEM(dem); //this is not a big deal since it will be in the buffer
	
493
494
495
496
497
	for (size_t param=0; param<v_params.size(); param++) { //loop over required parameters
		const MeteoGrids::Parameters grid_param = static_cast<MeteoGrids::Parameters>(v_params[param]);
		Grid2DObject grid;
		gridsmanager.read2DGrid(grid, grid_param, i_date);
		
498
499
500
		if (!grid.isSameGeolocalization(dem))
			throw InvalidArgumentException("In SMART_DOWNSCALING, the DEM and the source grid don't match for '"+MeteoGrids::getParameterName(grid_param)+"' on "+i_date.toString(Date::ISO));
		
501
502
503
504
505
		for (size_t ii=0; ii<v_stations.size(); ii++) { //loop over all virtual stations
			const size_t grid_i = v_stations[ii].position.getGridI(); //this should work since invalid stations have been removed in init
			const size_t grid_j = v_stations[ii].position.getGridJ();
			
			//check if this is a standard MeteoData parameter
Mathias Bavay's avatar
Mathias Bavay committed
506
			const  size_t meteo_param = vecMeteo[ii].getParameterIndex( MeteoGrids::getParameterName(grid_param) ); //is this name also a meteoparameter?
507
			if (meteo_param!=IOUtils::npos)
Mathias Bavay's avatar
Mathias Bavay committed
508
				vecMeteo[ii]( static_cast<MeteoData::Parameters>(meteo_param) ) = grid(grid_i, grid_j);
509
510
511
512
513
		}
	}
	
	return vecMeteo.size();
}
514

515
516
const std::string Meteo2DInterpolator::toString() const
{
517
	ostringstream os;
518
	os << "<Meteo2DInterpolator>\n";
519
	os << "Config& cfg = " << hex << &cfg << dec << "\n";
520
521
	os << "TimeSeriesManager& tsmanager = "  << hex << &tsmanager << dec << "\n";
	os << "GridsManager& gridsmanager = "  << hex << &gridsmanager << dec << "\n";
522

523
	os << "Spatial resampling algorithms:\n";
524
	std::map<std::string, std::vector<InterpolationAlgorithm*> >::const_iterator iter;
525
	for (iter = mapAlgorithms.begin(); iter != mapAlgorithms.end(); ++iter) {
526
		os << setw(10) << iter->first << "::";
527
		for (size_t jj=0; jj<iter->second.size(); jj++) {
528
			os << iter->second[jj]->algo << " ";
529
530
531
532
		}
		os << "\n";
	}

533
	//cache content
534
	os << grid_buffer.toString();
535
	os << "</Meteo2DInterpolator>\n";
536
	return os.str();
537
538
}

539
} //namespace