WSL/SLF GitLab Repository

MeteoData.cc 15.8 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
	paramname.push_back("TSG");
	paramname.push_back("TSS");
57
	paramname.push_back("TSOIL");
58
	paramname.push_back("P");
59
	paramname.push_back("P_SEA");
60
61
62
	paramname.push_back("U");
	paramname.push_back("V");
	paramname.push_back("W");
63
	paramname.push_back("SWE");
64
	paramname.push_back("ROT");
65
	paramname.push_back("ALB");
66
	paramname.push_back("DEM");
67
	paramname.push_back("SHADE");
68
69
	paramname.push_back("SLOPE");
	paramname.push_back("AZI");
70
71
72
73
74
75
76
77
78
79
80
81

	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];
}

82
83
84
85
86
87
88
89
90
91
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
}

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

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

119
120
121
	return true;
}

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

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

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

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

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

143
bool MeteoData::param_exists(const std::string& i_paramname) const
144
{
145
146
147
148
149
150
	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();
151
	for (size_t ii = 0; ii<current_size; ii++) {
152
		if (extra_param_name[ii] == i_paramname)
153
154
			return true;
	}
155

156
	return false;
157
158
}

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

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

	//Increase nrOfAllParameters
171
	nrOfAllParameters++;
172
173

	return (nrOfAllParameters - 1);
174
175
}

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

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

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

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

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

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

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

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

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

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

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

236
	for (size_t ii=0; ii<nrOfAllParameters; ii++) {
237
238
		//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 ....
239
		if ( !IOUtils::checkEpsilonEquality(data[ii], in.data[ii], epsilon) ){
240
			return false;
241
		}
242
243
	}

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

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

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

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

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

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

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

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

288
size_t MeteoData::getParameterIndex(const std::string& parname) const
289
{
290
291
292
293
294
295
296
297
	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)
298
			return ii+MeteoData::nrOfParameters;
299
	}
300
301
302
303

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

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

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

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

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

331
332
333
334
335
336
	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));
337
338
339
	return os;
}

340
341
342
343
344
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));
345
	data.extra_param_name.resize(s_vector);
346
	for (size_t ii=0; ii<s_vector; ii++) {
347
348
		size_t s_string;
		is.read(reinterpret_cast<char*>(&s_string), sizeof(size_t));
349
350
		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]));
351
352
353
354
355
356
357
358
359
	}

	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));
360
	return is;
361
362
}

363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
void MeteoData::mergeTimeSeries(std::vector<MeteoData>& vec1, const std::vector<MeteoData>& vec2)
{
	if (vec1.empty() || vec2.empty()) return;
	if (vec1.back().date<vec2.front().date) return; //vec1 is before vec2
	if (vec1.front().date>vec2.back().date) return; //vec1 is after vec2
	
	size_t idx2 = 0;
	for (size_t ii=0; ii<vec1.size(); ii++) {
		while ((vec1[ii].date>vec2[idx2].date) && (idx2<vec2.size()))
			idx2++;
		
		if (idx2==vec2.size()) return; //no more chances of common timestamps
		
		if (vec1[ii].date==vec2[idx2].date) {//we found a common timestamp
			vec1[ii].merge( vec2[idx2] );
		}
	}
}

382
383
void MeteoData::merge(std::vector<MeteoData>& vec1, const std::vector<MeteoData>& vec2, const bool& simple_merge)
{
384
	if (vec2.empty()) return;
385

386
	if (simple_merge || vec1.empty()) {
387
		vec1.reserve( vec1.size()+vec2.size() );
388
		for (size_t ii=0; ii<vec2.size(); ii++) vec1.push_back( vec2[ii] );
389
	} else {
390
		for (size_t ii=0; ii<vec2.size(); ii++) merge(vec1, vec2[ii]);
391
392
393
394
395
	}
}

void MeteoData::merge(std::vector<MeteoData>& vec, const MeteoData& meteo2, const bool& simple_merge)
{
396
397
	if (!simple_merge) {
		for (size_t ii=0; ii<vec.size(); ii++) {
398
			//two stations are considered the same if they point to the same 3D position
399
			if (vec[ii].meta.position==meteo2.meta.position) {
400
401
402
403
404
405
406
407
408
409
				vec[ii].merge(meteo2);
				return;
			}
		}
	}

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

410
411
412
413
414
415
416
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);
417
	for (size_t ii=0; ii<nElems; ii++) mergeIdx[ii] = ii;
418
419
420
	
	for (size_t ii=0; ii<nElems; ii++) {
		if (mergeIdx[ii]==IOUtils::npos) continue; //this element has already been merged, skip
421
		for (size_t jj=ii+1; jj<nElems; jj++) {
422
423
424
425
426
427
428
429
430
431
432
			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 );
}

433
MeteoData MeteoData::merge(MeteoData meteo1, const MeteoData& meteo2)
434
{
435
436
	meteo1.merge(meteo2);
	return meteo1;
437
438
439
440
}

void MeteoData::merge(const MeteoData& meteo2)
{
441
	if (!date.isUndef() && !meteo2.date.isUndef() && date!=meteo2.date) {
442
		//the data must be time synchronized!
443
		std::ostringstream ss;
444
445
446
447
448
		ss << "Trying to merge MeteoData at " << date.toString(Date::ISO);
		ss << " with MeteoData at " << meteo2.date.toString(Date::ISO);
		throw InvalidArgumentException(ss.str(), AT);
	}

449
	if (date.isUndef()) date=meteo2.date;
450
451
	meta.merge(meteo2.meta);

452
	if (meteo2.resampled==true ) resampled=true;
453
454

	//merge standard parameters
455
456
	for (size_t ii=0; ii<nrOfParameters; ii++) {
		if (data[ii]==IOUtils::nodata) data[ii]=meteo2.data[ii];
457
458
459
460
461
462
463
	}

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

	//no extra parameters to add -> return
464
	if (nrExtra2==0) return;
465
466

	//extra parameters only in meteo2 -> add them
467
468
	if (nrExtra1==0) {
		for (size_t ii=0; ii<nrExtra2; ii++) {
469
			const size_t new_idx = addParameter( meteo2.extra_param_name[ii] );
470
471
472
473
474
475
476
477
			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
478
	for (size_t ii=0; ii<nrExtra1; ii++) {
479
		const std::string curr_name = extra_param_name[ii];
480
		//look for this parameter in meteo2
481
		for (size_t jj=0; jj<nrExtra2; jj++) {
482
			if (meteo2.extra_param_name[jj]==curr_name ) {
483
				meteo2_flags[jj] = true;
484
				if (data[nrOfParameters+ii]==IOUtils::nodata) data[nrOfParameters+ii]=meteo2.data[nrOfParameters+jj];
485
486
487
488
489
				break;
			}
		}
	}
	//merge the extra fields of meteo2 that were NOT in the current meteodata
490
491
	for (size_t ii=0; ii<nrExtra2; ii++) {
		if (meteo2_flags[ii]==false) {
492
			const size_t new_idx = addParameter( meteo2.extra_param_name[ii] );
493
494
495
496
497
			data[new_idx] = meteo2.data[nrOfParameters+ii];
		}
	}
}

498
} //namespace