WSL/SLF GitLab Repository

Meteo2DInterpolator.cc 6.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
/***********************************************************************************/
/*  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/Meteo2DInterpolator.h>

using namespace std;

namespace mio {

25
Meteo2DInterpolator::Meteo2DInterpolator(const Config& _cfg, const DEMObject& _dem, 
26
27
28
29
30
								 const std::vector<MeteoData>& _vecMeteo, 
								 const std::vector<StationData>& _vecStation) 
	: cfg(_cfg), dem(_dem), vecMeteo(_vecMeteo), vecStation(_vecStation)
{
	/*
31
	 * By reading the Config object build up a list of user configured algorithms
32
33
	 * 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
34
	 * for configuration of interpolation algorithms within the Config object.
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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
	 */
	for (unsigned int ii=0; ii < MeteoData::nrOfParameters; ii++){ //loop over all MeteoData member variables
		std::vector<std::string> tmpAlgorithms;
		const std::string& parname = MeteoData::getParameterName(ii); //Current parameter name
		unsigned int nrOfAlgorithms = getAlgorithmsForParameter(parname, tmpAlgorithms);

		if (nrOfAlgorithms > 0)
			mapAlgorithms[parname] = tmpAlgorithms;
	}

	//check whether the size of the two vectors is equal
	if (vecMeteo.size() != vecStation.size())
		throw IOException("Size of vector<MeteoData> and vector<StationData> are not equal", AT);

	//check that the stations are using the same projection as the dem
	for (unsigned int i=0; i<(unsigned int)vecStation.size(); i++) {
		if(!vecStation[i].position.isSameProj(dem.llcorner)) {
			throw IOException("Some stations are not using the same geographic projection as the DEM", AT);
		}
	}
}

void Meteo2DInterpolator::interpolate(const MeteoData::Parameters& meteoparam, Grid2DObject& result) const
{
	//Show algorithms to be used for this parameter
	map<string, vector<string> >::const_iterator it = mapAlgorithms.find(MeteoData::getParameterName(meteoparam));
	
	if (it != mapAlgorithms.end()){
		double maxQualityRating = 0.0;
		auto_ptr<InterpolationAlgorithm> bestalgorithm(NULL);
		vector<string> vecArgs;

		for (unsigned int ii=0; ii < it->second.size(); ii++){
			const string& algoname = it->second.at(ii);
			getArgumentsForAlgorithm(meteoparam, algoname, vecArgs);
			
			//Get the configured algorithm
			auto_ptr<InterpolationAlgorithm> algorithm(AlgorithmFactory::getAlgorithm(algoname, *this, dem, 
			                                            vecMeteo, vecStation, vecArgs));
			//Get the quality rating and compare to previously computed quality ratings
			const double rating = algorithm->getQualityRating(meteoparam);
			if ((rating != 0.0) && (rating > maxQualityRating)){
				//we use ">" so that in case of equality, the first choice will be kept
				bestalgorithm = algorithm; //remember this algorithm: ownership belongs to bestalgorithm
				maxQualityRating = rating;
			}
		}

		//finally execute the algorithm with the best quality rating or throw an exception
		if (bestalgorithm.get() == NULL) {
			throw IOException("No interpolation algorithm with quality rating >0 found", AT);
		}
		bestalgorithm->calculate(meteoparam, result);
	} else {
		//Some default message, that interpolation for this parameter needs configuration
		throw IOException("You need to configure the interpolation algorithms for parameter " + 
					   MeteoData::getParameterName(meteoparam), AT);
		
	}

	//check that the output grid is using the same projection as the dem
	if(!result.llcorner.isSameProj(dem.llcorner)) {
		throw IOException("The output grid is not using the same geographic projection as the DEM", AT);
	}
99
100
101
102
103
104
105

	//Run soft min/max filter for RH and HNW
	if (meteoparam == MeteoData::RH){
		Meteo2DInterpolator::checkMinMax(0.0, 1.0, result);
	} else if (meteoparam == MeteoData::HNW){
		Meteo2DInterpolator::checkMinMax(0.0, 10000.0, result);
	}
106
107
108
109
110
111
}

unsigned int Meteo2DInterpolator::getAlgorithmsForParameter(const std::string& parname, std::vector<std::string>& vecAlgorithms)
{
	/* 
	 * This function retrieves the user defined interpolation algorithms for 
112
	 * parameter 'parname' by querying the Config object
113
114
115
116
117
118
119
120
121
122
123
124
125
	 */
	std::vector<std::string> vecKeys;

	vecAlgorithms.clear();
	cfg.findKeys(vecKeys, parname+"::algorithms", "Interpolations2D");

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

	if (vecKeys.size() == 0)
		return 0;


126
	cfg.getValue(vecKeys.at(0), "Interpolations2D", vecAlgorithms, Config::nothrow);
127
128
129
130
131
132
133
134
135
136

	return vecAlgorithms.size();
}

unsigned int Meteo2DInterpolator::getArgumentsForAlgorithm(const MeteoData::Parameters& param, 
                                                           const std::string& algorithm,
                                                           std::vector<std::string>& vecArgs) const
{
	vecArgs.clear();
	const string keyname = MeteoData::getParameterName(param) +"::"+ algorithm;
137
	cfg.getValue(keyname, "Interpolations2D", vecArgs, Config::nothrow);
138
139
140
141

	return vecArgs.size();
}

142
143
144
145
146
147
void Meteo2DInterpolator::checkMinMax(const double& minval, const double& maxval, Grid2DObject& gridobj)
{
	unsigned int nx=0, ny=0;

	gridobj.grid2D.size(nx, ny);

148
149
150
151
	for (unsigned int jj=0; jj<ny; jj++){
		for (unsigned int ii=0; ii<nx; ii++){
			double& value = gridobj.grid2D(ii,jj);
			if (value == IOUtils::nodata){
152
				//Do nothing
153
154
155
156
			} else if (value < minval){
				value = minval;
			} else if (value > maxval){
				value = maxval;
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
#ifdef _POPC_
#include "marshal_meteoio.h"
void Meteo2DInterpolator::Serialize(POPBuffer &buf, bool pack)
{
//TODO: check this serialization!! It seems dubious that it would work at all...
	if (pack)
	{
		//buf.Pack(&cfg,1);
		/*buf.Pack(&dem,1);
		marshal_METEO_DATASET(buf, vecMeteo, 0, FLAG_MARSHAL, NULL);
		marshal_STATION_DATASET(buf, vecStation, 0, FLAG_MARSHAL, NULL);*/
		marshal_map_str_vecstr(buf, mapAlgorithms, 0, FLAG_MARSHAL, NULL);
	}
	else
	{
		//buf.UnPack(&cfg,1);
		/*buf.UnPack(&dem,1);
		marshal_METEO_DATASET(buf, vecMeteo, 0, !FLAG_MARSHAL, NULL);
		marshal_STATION_DATASET(buf, vecStation, 0, !FLAG_MARSHAL, NULL);*/
		marshal_map_str_vecstr(buf, mapAlgorithms, 0, !FLAG_MARSHAL, NULL);
	}
}
#endif

} //namespace