WSL/SLF GitLab Repository

MeteoData.cc 15.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/***********************************************************************************/
/*  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/>.
*/
18
19
#include <meteoio/dataClasses/MeteoData.h>
#include <meteoio/dataClasses/StationData.h>
20

21
22
23
#include <cmath>
#include <limits>

24
25
26
using namespace std;
namespace mio {

27
28
29
30
31
32
33
34
35
36
37
38
/************************************************************
 * static section                                           *
 ************************************************************/
const size_t MeteoGrids::nrOfParameters =  MeteoGrids::lastparam - MeteoGrids::firstparam + 1;
std::vector<std::string> MeteoGrids::paramname;
const bool MeteoGrids::__init = MeteoGrids::initStaticData();

bool MeteoGrids::initStaticData()
{
	//the order must be the same as in the enum
	paramname.push_back("TA");
	paramname.push_back("RH");
39
40
	paramname.push_back("QI");
	paramname.push_back("TD");
41
42
43
44
45
	paramname.push_back("VW");
	paramname.push_back("DW");
	paramname.push_back("VW_MAX");
	paramname.push_back("ISWR");
	paramname.push_back("RSWR");
46
47
	paramname.push_back("ISWR_DIFF");
	paramname.push_back("ISWR_DIR");
48
	paramname.push_back("ILWR");
49
	paramname.push_back("TAU_CLD");
50
	paramname.push_back("HS");
51
52
53
54
	paramname.push_back("PSUM");
	paramname.push_back("PSUM_PH");
	paramname.push_back("PSUM_L");
	paramname.push_back("PSUM_S");
55
56
57
	paramname.push_back("TSG");
	paramname.push_back("TSS");
	paramname.push_back("P");
58
	paramname.push_back("P_SEA");
59
60
61
	paramname.push_back("U");
	paramname.push_back("V");
	paramname.push_back("W");
62
	paramname.push_back("SWE");
63
	paramname.push_back("ROT");
64
	paramname.push_back("ALB");
65
	paramname.push_back("DEM");
66
	paramname.push_back("SHADE");
67
68
	paramname.push_back("SLOPE");
	paramname.push_back("AZI");
69
70
71
72
73
74
75
76
77
78
79
80

	return true;
}

const std::string& MeteoGrids::getParameterName(const size_t& parindex)
{
	if (parindex >= MeteoGrids::nrOfParameters)
		throw IndexOutOfBoundsException("Trying to get name for parameter that does not exist", AT);

	return paramname[parindex];
}

81
82
83
84
85
86
87
88
89
90
size_t MeteoGrids::getParameterIndex(const std::string& parname)
{
	for (size_t ii=0; ii<MeteoGrids::nrOfParameters; ii++) {
		if (paramname[ii] == parname)
			return ii;
	}

	return IOUtils::npos; //parameter not a part of MeteoGrids
}

91
92
93
/************************************************************
 * static section                                           *
 ************************************************************/
94
const double MeteoData::epsilon = 1e-5;
95
const size_t MeteoData::nrOfParameters =  MeteoData::lastparam - MeteoData::firstparam + 1;
96
std::vector<std::string> MeteoData::s_default_paramname;
97
98
99
100
const bool MeteoData::__init = MeteoData::initStaticData();

bool MeteoData::initStaticData()
{
101
	//Since the parameters enum starts at 0, this is enough to associate an index with its name
102
	s_default_paramname.push_back("P");
103
104
	s_default_paramname.push_back("TA");
	s_default_paramname.push_back("RH");
105
106
107
	s_default_paramname.push_back("TSG");
	s_default_paramname.push_back("TSS");
	s_default_paramname.push_back("HS");
108
109
110
111
	s_default_paramname.push_back("VW");
	s_default_paramname.push_back("DW");
	s_default_paramname.push_back("VW_MAX");
	s_default_paramname.push_back("RSWR");
112
	s_default_paramname.push_back("ISWR");
113
	s_default_paramname.push_back("ILWR");
114
	s_default_paramname.push_back("TAU_CLD");
115
	s_default_paramname.push_back("PSUM");
116
	s_default_paramname.push_back("PSUM_PH");
117

118
119
120
	return true;
}

121
const std::string& MeteoData::getParameterName(const size_t& parindex)
122
123
124
125
{
	if (parindex >= MeteoData::nrOfParameters)
		throw IndexOutOfBoundsException("Trying to access meteo parameter that does not exist", AT);

126
	return MeteoData::s_default_paramname[parindex];
127
128
129
130
131
132
}

/************************************************************
 * non-static section                                       *
 ************************************************************/

133
const std::string& MeteoData::getNameForParameter(const size_t& parindex) const
134
{
135
	if (parindex >= MeteoData::nrOfAllParameters)
136
		throw IndexOutOfBoundsException("Trying to get name for parameter that does not exist", AT);
137

138
139
	if (parindex<MeteoData::nrOfParameters) return MeteoData::s_default_paramname[parindex];
	return extra_param_name[parindex-MeteoData::nrOfParameters];
140
141
}

142
bool MeteoData::param_exists(const std::string& i_paramname) const
143
{
144
145
146
147
148
149
	for (size_t ii = 0; ii<MeteoData::nrOfParameters; ii++) {
		if (s_default_paramname[ii] == i_paramname)
			return true;
	}
	
	const size_t current_size = extra_param_name.size();
150
	for (size_t ii = 0; ii<current_size; ii++) {
151
		if (extra_param_name[ii] == i_paramname)
152
153
			return true;
	}
154

155
	return false;
156
157
}

158
size_t MeteoData::addParameter(const std::string& i_paramname)
159
{
160
	//check if name is already taken
161
	const size_t current_index = getParameterIndex(i_paramname);
162
163
	if (current_index != IOUtils::npos)
		return current_index; //do nothing, because parameter is already present
164
165

	//add parameter
166
	extra_param_name.push_back(i_paramname);
167
	data.push_back(IOUtils::nodata);
168
169

	//Increase nrOfAllParameters
170
	nrOfAllParameters++;
171
172

	return (nrOfAllParameters - 1);
173
174
}

175
size_t MeteoData::getNrOfParameters() const
176
177
{
	return nrOfAllParameters;
178
179
}

180
MeteoData::MeteoData()
181
         : date(0.0, 0.), meta(), extra_param_name(), data(MeteoData::nrOfParameters, IOUtils::nodata), nrOfAllParameters(MeteoData::nrOfParameters), resampled(false)
182
{ }
183

184
MeteoData::MeteoData(const Date& date_in)
185
         : date(date_in), meta(), extra_param_name(), data(MeteoData::nrOfParameters, IOUtils::nodata), nrOfAllParameters(MeteoData::nrOfParameters), resampled(false)
186
{ }
187

188
MeteoData::MeteoData(const Date& date_in, const StationData& meta_in)
189
         : date(date_in), meta(meta_in), extra_param_name(), data(MeteoData::nrOfParameters, IOUtils::nodata), nrOfAllParameters(MeteoData::nrOfParameters), resampled(false)
190
191
{ }

192
void MeteoData::setDate(const Date& in_date)
193
{
194
	date = in_date;
195
196
}

197
198
void MeteoData::reset()
{
199
	std::fill(data.begin(), data.end(), IOUtils::nodata);
200
201
}

202
203
204
205
206
/**
* @brief Standardize nodata values
* The plugin-specific nodata values are replaced by MeteoIO's internal nodata value
*/
void MeteoData::standardizeNodata(const double& plugin_nodata) {
207
	for (size_t ii=0; ii<nrOfAllParameters; ii++) {
208
		//loop through all meteo params and check whether they're nodata values
209
210
		if (data[ii] <= plugin_nodata)
			data[ii] = IOUtils::nodata;
211
212
213
	}
}

214
bool MeteoData::isResampled() const
215
216
217
218
{
	return resampled;
}

219
void MeteoData::setResampled(const bool& in_resampled)
220
{
221
	resampled = in_resampled;
222
223
224
225
226
}

bool MeteoData::operator==(const MeteoData& in) const
{
	//An object is equal if the date is equal and all meteo parameters are equal
227
	if (date != in.date) {
228
		return false;
229
	}
230

231
	if (nrOfAllParameters != in.nrOfAllParameters) { //the number of meteo parameters has to be consistent
232
		return false;
233
	}
234

235
	for (size_t ii=0; ii<nrOfAllParameters; ii++) {
236
237
		//const double epsilon_rel = (fabs(data[ii]) < fabs(in.data[ii]) ? fabs(in.data[ii]) : fabs(data[ii])) * std::numeric_limits<double>::epsilon(); // Hack not working...
		//const double epsilon_rel = (fabs(data[ii]) < fabs(in.data[ii]) ? fabs(in.data[ii]) : fabs(data[ii])) * 0.0000001; // Hack not working with 0 == 0 ....
238
		if( !IOUtils::checkEpsilonEquality(data[ii], in.data[ii], epsilon) ){
239
			return false;
240
		}
241
242
	}

243
	return true;
244
245
246
247
248
249
250
}

bool MeteoData::operator!=(const MeteoData& in) const
{
	return !(*this==in);
}

251
double& MeteoData::operator()(const size_t& parindex)
252
253
{
#ifndef NOSAFECHECKS
254
	if (parindex >= nrOfAllParameters)//getNrOfParameters())
255
		throw IndexOutOfBoundsException("Trying to access meteo parameter that does not exist", AT);
256
#endif
257
	return data[parindex];
258
259
}

260
const double& MeteoData::operator()(const size_t& parindex) const
261
262
{
#ifndef NOSAFECHECKS
263
	if (parindex >= nrOfAllParameters)//getNrOfParameters())
264
265
		throw IndexOutOfBoundsException("Trying to access meteo parameter that does not exist", AT);
#endif
266
	return data[parindex];
267
268
}

269
double& MeteoData::operator()(const std::string& parname)
270
{
271
	const size_t index = getParameterIndex(parname);
272
273
	if (index == IOUtils::npos)
		throw IndexOutOfBoundsException("Trying to access meteo parameter that does not exist: " + parname, AT);
274

275
	return operator()(index);
276
277
}

278
const double& MeteoData::operator()(const std::string& parname) const
279
{
280
	const size_t index = getParameterIndex(parname);
281
282
	if (index == IOUtils::npos)
		throw IndexOutOfBoundsException("Trying to access meteo parameter that does not exist: " + parname, AT);
283

284
	return operator()(index);
285
286
}

287
size_t MeteoData::getParameterIndex(const std::string& parname) const
288
{
289
290
291
292
293
294
295
296
	for (size_t ii = 0; ii<MeteoData::nrOfParameters; ii++) {
		if (s_default_paramname[ii] == parname)
			return ii;
	}
	
	const size_t current_size = extra_param_name.size();
	for (size_t ii=0; ii<current_size; ii++) {
		if (extra_param_name[ii] == parname)
297
			return ii;
298
	}
299
300
301
302

	return IOUtils::npos; //parameter not a part of MeteoData
}

303
const std::string MeteoData::toString() const {
304
	std::ostringstream os;
305
	os << "<meteo>\n";
306
307
	os << meta.toString();
	os << date.toString(Date::FULL) << "\n";
308

309
310
	for (size_t ii=0; ii<getNrOfParameters(); ii++) {
		const double& value = operator()(ii);
311
		if (value != IOUtils::nodata)
312
			os << setw(8) << getNameForParameter(ii) << ":" << setw(15) << value << endl;
313
	}
314

315
	os << "</meteo>\n";
316
317
318
319
320
321
	return os.str();
}

std::iostream& operator<<(std::iostream& os, const MeteoData& data) {
	os << data.date;
	os << data.meta;
322
	const size_t s_vector = data.extra_param_name.size();
323
	os.write(reinterpret_cast<const char*>(&s_vector), sizeof(size_t));
324
	for (size_t ii=0; ii<s_vector; ii++) {
325
		const size_t s_string = data.extra_param_name[ii].size();
326
		os.write(reinterpret_cast<const char*>(&s_string), sizeof(size_t));
327
		os.write(reinterpret_cast<const char*>(&data.extra_param_name[ii][0]), s_string*sizeof(data.extra_param_name[ii][0]));
328
	}
329

330
331
332
333
334
335
	const size_t s_data = data.data.size();
	os.write(reinterpret_cast<const char*>(&s_data), sizeof(s_data));
	os.write(reinterpret_cast<const char*>(&data.data[0]), s_data*sizeof(data.data[0]));

	os.write(reinterpret_cast<const char*>(&data.nrOfAllParameters), sizeof(data.nrOfAllParameters));
	os.write(reinterpret_cast<const char*>(&data.resampled), sizeof(data.resampled));
336
337
338
	return os;
}

339
340
341
342
343
std::iostream& operator>>(std::iostream& is, MeteoData& data) {
	is >> data.date;
	is >> data.meta;
	size_t s_vector;
	is.read(reinterpret_cast<char*>(&s_vector), sizeof(size_t));
344
	data.extra_param_name.resize(s_vector);
345
	for (size_t ii=0; ii<s_vector; ii++) {
346
347
		size_t s_string;
		is.read(reinterpret_cast<char*>(&s_string), sizeof(size_t));
348
349
		data.extra_param_name[ii].resize(s_string);
		is.read(reinterpret_cast<char*>(&data.extra_param_name[ii][0]), s_string*sizeof(data.extra_param_name[ii][0]));
350
351
352
353
354
355
356
357
358
	}

	size_t s_data;
	is.read(reinterpret_cast<char*>(&s_data), sizeof(size_t));
	data.data.resize(s_data);
	is.read(reinterpret_cast<char*>(&data.data[0]), s_data*sizeof(data.data[0]));

	is.read(reinterpret_cast<char*>(&data.nrOfAllParameters), sizeof(data.nrOfAllParameters));
	is.read(reinterpret_cast<char*>(&data.resampled), sizeof(data.resampled));
359
	return is;
360
361
}

362
363
void MeteoData::merge(std::vector<MeteoData>& vec1, const std::vector<MeteoData>& vec2, const bool& simple_merge)
{
364
	if (vec2.empty()) return;
365

366
	if (simple_merge || vec1.empty()) {
367
		vec1.reserve( vec1.size()+vec2.size() );
368
369
		for(size_t ii=0; ii<vec2.size(); ii++) vec1.push_back( vec2[ii] );
	} else {
370
		for (size_t ii=0; ii<vec2.size(); ii++) merge(vec1, vec2[ii]);
371
372
373
374
375
	}
}

void MeteoData::merge(std::vector<MeteoData>& vec, const MeteoData& meteo2, const bool& simple_merge)
{
376
377
	if (!simple_merge) {
		for (size_t ii=0; ii<vec.size(); ii++) {
378
			//two stations are considered the same if they point to the same 3D position
379
			if (vec[ii].meta.position==meteo2.meta.position) {
380
381
382
383
384
385
386
387
388
389
				vec[ii].merge(meteo2);
				return;
			}
		}
	}

	//the station was not found in the vector -> adding it
	vec.push_back( meteo2 );
}

390
391
392
393
394
395
396
397
398
399
400
void MeteoData::merge(std::vector<MeteoData>& vec)
{
	const size_t nElems = vec.size();
	if (nElems<2) return;
	
	std::vector<MeteoData> vecResult;
	std::vector<size_t> mergeIdx(nElems);
	for(size_t ii=0; ii<nElems; ii++) mergeIdx[ii] = ii;
	
	for (size_t ii=0; ii<nElems; ii++) {
		if (mergeIdx[ii]==IOUtils::npos) continue; //this element has already been merged, skip
401
		for (size_t jj=ii+1; jj<nElems; jj++) {
402
403
404
405
406
407
408
409
410
411
412
			if (vec[ii].meta.position==vec[jj].meta.position) {
				vec[ii].merge( vec[jj] );
				mergeIdx[jj]=IOUtils::npos; //this element will be skipped in the next loops
			}
		}
		vecResult.push_back( vec[ii] );
	}
	
	vec.swap( vecResult );
}

413
MeteoData MeteoData::merge(MeteoData meteo1, const MeteoData& meteo2)
414
{
415
416
	meteo1.merge(meteo2);
	return meteo1;
417
418
419
420
}

void MeteoData::merge(const MeteoData& meteo2)
{
421
	if (!date.isUndef() && !meteo2.date.isUndef() && date!=meteo2.date) {
422
		//the data must be time synchronized!
423
		std::ostringstream ss;
424
425
426
427
428
		ss << "Trying to merge MeteoData at " << date.toString(Date::ISO);
		ss << " with MeteoData at " << meteo2.date.toString(Date::ISO);
		throw InvalidArgumentException(ss.str(), AT);
	}

429
	if (date.isUndef()) date=meteo2.date;
430
431
	meta.merge(meteo2.meta);

432
	if (meteo2.resampled==true ) resampled=true;
433
434

	//merge standard parameters
435
436
	for (size_t ii=0; ii<nrOfParameters; ii++) {
		if (data[ii]==IOUtils::nodata) data[ii]=meteo2.data[ii];
437
438
439
440
441
442
443
	}

	//merge extra parameters
	const size_t nrExtra1 = nrOfAllParameters - nrOfParameters;
	const size_t nrExtra2 = meteo2.nrOfAllParameters - nrOfParameters;

	//no extra parameters to add -> return
444
	if (nrExtra2==0) return;
445
446

	//extra parameters only in meteo2 -> add them
447
448
	if (nrExtra1==0) {
		for (size_t ii=0; ii<nrExtra2; ii++) {
449
			const size_t new_idx = addParameter( meteo2.extra_param_name[nrOfParameters+ii] );
450
451
452
453
454
455
456
457
			data[new_idx] = meteo2.data[nrOfParameters+ii];
		}
		return;
	}

	//extra parameters in both the current meteodata and meteo2 -> tedious merge...
	std::vector<bool> meteo2_flags(nrExtra2, false); //to keep track of which elements have been copied
	//merge the extra field of the current meteodata
458
	for (size_t ii=0; ii<nrExtra1; ii++) {
459
		const std::string curr_name = extra_param_name[nrOfParameters+ii];
460
		//look for this parameter in meteo2
461
		for (size_t jj=0; jj<nrExtra2; jj++) {
462
			if (meteo2.extra_param_name[nrOfParameters+jj]==curr_name ) {
463
				meteo2_flags[jj] = true;
464
				if (data[nrOfParameters+ii]==IOUtils::nodata) data[nrOfParameters+ii]=meteo2.data[nrOfParameters+jj];
465
466
467
468
469
				break;
			}
		}
	}
	//merge the extra fields of meteo2 that were NOT in the current meteodata
470
471
	for (size_t ii=0; ii<nrExtra2; ii++) {
		if (meteo2_flags[ii]==false) {
472
			const size_t new_idx = addParameter( meteo2.extra_param_name[nrOfParameters+ii] );
473
474
475
476
477
			data[new_idx] = meteo2.data[nrOfParameters+ii];
		}
	}
}

478
} //namespace