WSL/SLF GitLab Repository

OshdIO.cc 24.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
#include <meteoio/plugins/OshdIO.h>
19
#include <meteoio/meteoLaws/Meteoconst.h>
20
21
#include <meteoio/FileUtils.h>

22
#include <matio.h>
23
#include <algorithm>
24
25
26
27
28
29
30
31
32
33
34
35
36
37

using namespace std;

namespace mio {
/**
 * @page oshd OshdIO
 * This plugin reads the meteorological forecast data downscaled for each of the Swiss meteorological networks IMIS/ANETZ stations
 * as preprocessed by the <A HREF="www.wsl.ch/fe/gebirgshydrologie/schnee_hydro/oshd/index_EN">Operational Snow-Hydrological Service</A> 
 * of the <A HREF="www.wsl.ch">WSL</A>. The data is written as Matlab 
 * <A HREF="http://www.mathworks.com/help/pdf_doc/matlab/matfile_format.pdf">binary files (.mat)</A>, one per meteorological parameter and per timestep, 
 * available on an access-controlled server after each new <A HREF="www.cosmo-model.org/">COSMO</A> run. It therefore requires a third party 
 * library to read this file format: the Open Source <A HREF="https://sourceforge.net/projects/matio/">MatIO</A> library. This can be installed directly from
 * the default repositories under Linux or installed by downloading the proper package for Windows or OsX.
 * 
38
39
40
41
42
 * \note If non-ascii characters have been used and the file has been created under Windows, some of the strings might end up using the UTF-16 encoding.
 * This requires a recent version of libmatio (see <A HREF="https://github.com/tbeu/matio/issues/34">this issue</A>). Another option would be to 
 * add at the begining of the Matlab routine a call to *feature('DefaultCharacterSet', 'UTF8')* in order to switch from the current default (which can be read by the same call, 
 * ommitting the 'UTF8' option) to the (<A HREF="http://blog.omega-prime.co.uk/?p=150">partial</A>) UTF-8  encoding of Matlab.
 * 
43
 * @section oshd_data_structure Data structure
44
 * The files are named with the following schema: <i>{parameter}_{timestep}_{cosmo model version}_F_{run time}.mat</i> with the following possible values:
45
 *     + *parameter* is one of idif, idir, albd, ilwr, pair, prec, rhum, tcor, wcor, wdir;
46
 *     + *timestep* is written as purely numeric ISO with minute resolution;
47
 *     + *cosmo model version* could be any of cosmo7, cosmo2, cosmo1, cosmoE;
48
49
50
51
52
53
54
55
56
57
58
59
60
 *     + *run time* is the purely numeric ISO date and time of when COSMO produced the dataset.
 * 
 * The files have the following internal data structure (represented as "name {data type}"):
 * @verbatim
      stat {1x1 struct}
        ├── time {1x1 array of doubles}
        ├── data {1x623 array of doubles}
        ├── acro {1x623 array of arrays of char}
        ├── dunit {array of char}
        ├── type {array of char}
        └── name {array of char}
  @endverbatim
 * 
61
62
63
64
65
66
67
68
69
70
71
 * The stations' acronyms follow a fixed order but their coordinates must be provided in a separate file, given as *METAFILE* key (see below). This file
 * must have the following structure (the *x* and *y* coordinates being the CH1903 easting and northing, respectively): 
 * @verbatim
      statlist {1x1 struct}
        ├── acro {1x623 array of arrays of char}
        ├── name {1x623 array of arrays of char}
        ├── x {1x623 array of doubles}
        ├── y {1x623 array of doubles}
        └── z {1x623 array of doubles}
  @endverbatim
 *
72
73
74
75
76
 *
 * @section oshd_keywords Keywords
 * This plugin uses the following keywords:
 * - COORDSYS: coordinate system (see Coords); [Input] and [Output] section
 * - COORDPARAM: extra coordinates parameters (see Coords); [Input] and [Output] section
77
 * - METEOPATH: directory containing all the data files with the proper file naming schema; [Input] section
78
 * - METEOPATH_RECURSIVE: should *meteopath* be searched recursively for files? (default: no); [Input] section
79
 * - STATION#: input stations' IDs (in METEOPATH). As many meteofiles as needed may be specified
80
81
 * - METAFILE: file containing the stations' IDs, names and location; [Input] section (either within METEOPATH if not path is 
 provided or within the provided path)
82
 * - OSHD_DEBUG: write out extra information to better show what is in the files
83
84
85
86
87
88
 *
 * @section oshd_example Example use
 * @code
 * [Input]
 * METEO = OSHD
 * METEOPATH = /local/LATEST_03h_RUN
89
 * METEOPATH_RECURSIVE = true
90
 * METAFILE  = STAT_LIST.mat ;another possibility could be /local/metadata/STAT_LIST.mat
91
92
93
94
 * STATION1  = ATT2
 * STATION2  = WFJ2
 * @endcode
 *
95
96
97
98
99
100
101
102
103
104
 */

/**********************************************************************************
 * Here we define some wrappers around libmatio. These are not declared as class methods in
 * order to avoid having to expose matio.h when including OshdIO.h
 **********************************************************************************/
void listFields(matvar_t *matvar)
{
	const unsigned int nrFields = Mat_VarGetNumberOfFields(matvar);
	char * const *fields = Mat_VarGetStructFieldnames(matvar);
105
106
	for (unsigned int ii=0; ii<nrFields; ii++) 
		printf("field[%d] = %s\n", ii, fields[ii]);
107
108
109
110
111
112
}

void printStructure(matvar_t *matvar)
{
	//Mat_VarPrint(field, 0);
	printf("name=%s class_type=%d data_type=%d rank=%d", matvar->name, matvar->class_type, matvar->data_type, matvar->rank);
113
114
	for (int ii=0; ii<matvar->rank; ii++) 
		printf("\tdims[%d]=%d", ii, (int)matvar->dims[ii]);
115
116
117
	printf("\n");
}

118

119
120
121
122
123
124
125
126
std::string readString(const std::string &filename, const std::string &fieldname, mat_t *matfp, matvar_t *matvar)
{
	matvar_t *field = Mat_VarGetStructFieldByName(matvar, fieldname.c_str(), 0);
	if (matvar==NULL) throw NotFoundException("could not read field '"+fieldname+"' in file '"+filename+"'", AT);
	if (Mat_VarReadDataAll(matfp, field)) 
		throw InvalidFormatException("could not read field '"+fieldname+"' in file '"+filename+"'", AT);
	if (field->class_type!=MAT_C_CHAR) throw InvalidFormatException("field '"+fieldname+"' in file '"+filename+"' is not a type string", AT);
	
127
	return std::string( static_cast<char*>(field->data) );
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
}

void readStringVector(const std::string &filename, const std::string &fieldname, mat_t *matfp, matvar_t *matvar, std::vector<std::string> &vecString)
{
	vecString.clear();
	
	matvar_t *field = Mat_VarGetStructFieldByName(matvar, fieldname.c_str(), 0);
	if (matvar==NULL) 	throw NotFoundException("could not read field '"+fieldname+"' in file '"+filename+"'", AT);
	if (Mat_VarReadDataAll(matfp, field)) 
		throw InvalidFormatException("could not read field '"+fieldname+"' in file '"+filename+"'", AT);
	
	if (field->class_type!=MAT_C_CELL) throw InvalidFormatException("field '"+fieldname+"' in file '"+filename+"' is not a cell type", AT);
	if (field->data_type!=MAT_T_CELL) throw InvalidFormatException("field '"+fieldname+"' in file '"+filename+"' is not a cell array data type", AT);
	
	matvar_t *cell = Mat_VarGetCell(matvar, 0);
	if (cell==NULL) throw InvalidFormatException("could not read data in field '"+fieldname+"' in file '"+filename+"'", AT);
	if (field->rank!=2) throw InvalidFormatException("invalid rank for field '"+fieldname+"' in file '"+filename+"'", AT);
	
	const size_t nrows = field->dims[0];
	const size_t ncols = field->dims[1];
	if (nrows!=1) throw InvalidFormatException("invalid nrows for field '"+fieldname+"' in file '"+filename+"'", AT);
	
	vecString.resize( ncols );
	for (size_t ii=0; ii<ncols; ii++) {
		cell = Mat_VarGetCell(field, static_cast<int>(ii));
		if (cell->rank!=2) throw InvalidFormatException("invalid cell rank in file '"+filename+"'", AT);
		if (cell->class_type!=MAT_C_CHAR) throw InvalidFormatException("field '"+fieldname+"' in file '"+filename+"' is not a type string", AT);
		vecString[ii] = static_cast<char*>(cell->data);
	}
}

void readDoubleVector(const std::string &filename, const std::string &fieldname, mat_t *matfp, matvar_t *matvar, std::vector<double> &vecData)
{
	vecData.clear();
	
	matvar_t *field = Mat_VarGetStructFieldByName(matvar, fieldname.c_str(), 0);
	if (matvar==NULL) 	throw NotFoundException("could not read field '"+fieldname+"' in file '"+filename+"'", AT);
	if (Mat_VarReadDataAll(matfp, field)) 
		throw InvalidFormatException("could not read field '"+fieldname+"' in file '"+filename+"'", AT);
	
	if (field->class_type!=MAT_C_DOUBLE) throw InvalidFormatException("field '"+fieldname+"' in file '"+filename+"' is not a double type", AT);
	if (field->rank!=2) throw InvalidFormatException("invalid rank for field '"+fieldname+"' in file '"+filename+"'", AT);
	
	const size_t nrows = field->dims[0];
	const size_t ncols = field->dims[1];
	if (nrows!=1) throw InvalidFormatException("invalid nrows for field '"+fieldname+"' in file '"+filename+"'", AT);
	
175
	const double* matData( static_cast<double*>( field->data ) );
176
177
178
179
180
181
	vecData.resize( ncols );
	for (size_t ii=0; ii<ncols; ii++) {
		vecData[ii] = matData[ii];
	}
}

182
183
184
185
186
void printFileStructure(const std::string& filename)
{
	mat_t *matfp = Mat_Open(filename.c_str(), MAT_ACC_RDONLY);
	if ( NULL == matfp ) throw AccessException(filename, AT);

Mathias Bavay's avatar
Mathias Bavay committed
187
	std::cout << "<" << FileUtils::getFilename( filename ) << ">\n";
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
	matvar_t *matvar;
	while ( (matvar = Mat_VarReadNextInfo(matfp)) != NULL ) {
		std::cout << "\t" << matvar->name << " [";
		for (int ii=0; ii<matvar->rank; ii++) {
			std::cout << (int)matvar->dims[ii];
			if (ii<(matvar->rank-1)) std::cout << "x";
		}
		std::cout << "]\n";

		const unsigned int nrFields = Mat_VarGetNumberOfFields(matvar);
		char * const *fields = Mat_VarGetStructFieldnames(matvar);
		for (unsigned int ii=0; ii<nrFields; ii++) {
			const std::string field_name( fields[ii] );
			matvar_t *field = Mat_VarGetStructFieldByName(matvar, field_name.c_str(), 0);
			const std::string prefix = (ii<(nrFields-1))? "├──" : "└──";
			std::cout << "\t" << prefix << field_name;
			if (field->class_type==MAT_C_CHAR)
Mathias Bavay's avatar
Mathias Bavay committed
205
				std::cout << " = \"" << readString(filename, field_name, matfp, matvar) << "\"";
206
207
			if (field->class_type==MAT_C_DOUBLE) {
				std::cout << " [";
Mathias Bavay's avatar
Mathias Bavay committed
208
				size_t count=1;
209
210
211
				for (int jj=0; jj<field->rank; jj++) {
					std::cout << field->dims[jj];
					if (jj<(field->rank-1)) std::cout << "x";
Mathias Bavay's avatar
Mathias Bavay committed
212
					count *= field->dims[jj];
213
214
				}
				std::cout << "]";
Mathias Bavay's avatar
Mathias Bavay committed
215
216
217
218
219
220
221
222
223
224
225
				if (count==1) {
					if (Mat_VarReadDataAll(matfp, field))
						throw InvalidFormatException("could not read field '"+field_name+"' in file '"+filename+"'", AT);
					const double val = static_cast<double*>(field->data)[0];
					if (field_name=="time") {
						Date timestep;
						timestep.setMatlabDate( val, OshdIO::in_dflt_TZ );
						std::cout << " = " << timestep.toString(Date::ISO_TZ);
					} else
						std::cout << " = " << val;
				}
226
227
228
229
230
231
232
233
234
235
236
237
238
			}
			if (field->class_type==MAT_C_CELL) {
				std::cout << " [";
				for (int jj=0; jj<field->rank; jj++) {
					std::cout << field->dims[jj];
					if (jj<(field->rank-1)) std::cout << "x";
				}
				std::cout << "]";
			}

			std::cout << "\n";
		}
	}
Mathias Bavay's avatar
Mathias Bavay committed
239
	std::cout << "</" << FileUtils::getFilename( filename ) << ">\n\n";
240
241
242
243
244
	Mat_VarFree(matvar);
	matvar = NULL;
	Mat_Close(matfp);
}

245
246
247
248
249
250
/**********************************************************************************
 * Now really implementing the OshdIO class
 **********************************************************************************/
const char* OshdIO::meteo_ext = ".mat";
const double OshdIO::in_dflt_TZ = 0.; //COSMO data is always GMT

251
OshdIO::OshdIO(const std::string& configfile) : cfg(configfile), cache_meteo_files(), vecMeta(), vecIDs(), params_map(), vecIdx(), 
252
               in_meteopath(), in_metafile(), debug(false)
253
254
255
256
{
	parseInputOutputSection();
}

257
OshdIO::OshdIO(const Config& cfgreader) : cfg(cfgreader), cache_meteo_files(), vecMeta(), vecIDs(), params_map(), vecIdx(), 
258
               in_meteopath(), in_metafile(), debug(false)
259
260
261
262
263
264
{
	parseInputOutputSection();
}

void OshdIO::parseInputOutputSection()
{
265
	cfg.getValue("OSHD_DEBUG", "INPUT", debug, IOUtils::nothrow);
266
	//IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
267
	
268
269
	cfg.getValues("STATION", "INPUT", vecIDs);
	cfg.getValue("METEOPATH", "Input", in_meteopath);
270
271
272
	bool is_recursive = false;
	cfg.getValue("METEOPATH_RECURSIVE", "Input", is_recursive, IOUtils::nothrow);
	scanMeteoPath(in_meteopath, is_recursive, cache_meteo_files);
273
	
274
	cfg.getValue("METAFILE", "INPUT", in_metafile);
275
	if (FileUtils::getFilename(in_metafile) == in_metafile) { //ie there is no path in the provided filename
276
277
278
		in_metafile = in_meteopath + "/" + in_metafile;
	}
	
279
	//fill the params mapping vector
280
	params_map.push_back( std::make_pair(MeteoData::ILWR, "ilwr") );
281
282
	params_map.push_back( std::make_pair(MeteoData::P, "pair") );
	params_map.push_back( std::make_pair(MeteoData::PSUM, "prec") );
283
284
	params_map.push_back( std::make_pair(MeteoData::RH, "rhum") );
	params_map.push_back( std::make_pair(MeteoData::TA, "tcor") ); //old:tair
285
286
287
288
	params_map.push_back( std::make_pair(MeteoData::VW, "wcor") );
	params_map.push_back( std::make_pair(MeteoData::DW, "wdir") );
}

289
void OshdIO::scanMeteoPath(const std::string& meteopath_in, const bool& is_recursive,  std::vector< std::pair<mio::Date,std::string> > &meteo_files)
290
291
{
	meteo_files.clear();
292
	std::list<std::string> dirlist( FileUtils::readDirectory(meteopath_in, meteo_ext, is_recursive) );
293
	std::list<std::string> prefix_list;
294

295
296
	//make sure each timestamp only appears once, ie remove duplicates
	std::list<std::string>::iterator it = dirlist.begin();
297
298
299
	while ((it != dirlist.end())) {
		const std::string filename( *it );
		const std::string::size_type date_pos = filename.find_first_of("0123456789", 0);
300
301
		if (date_pos!=string::npos) 
			prefix_list.push_back( filename.substr(date_pos) );
302
303
		it++;
	}
304
305
306
307
308
309
310
311
312
313
314
315
316
	
	prefix_list.sort();
	prefix_list.unique();
	
	//Convert date in every filename and cache it
	std::list<std::string>::const_iterator it_prefix = prefix_list.begin();
	while ((it_prefix != prefix_list.end())) {
		const std::string filename( *it_prefix );
		Date date;
		IOUtils::convertString(date, filename.substr(0,12), in_dflt_TZ);
		meteo_files.push_back( std::make_pair(date, filename) );
		it_prefix++;
	}
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
}

size_t OshdIO::getFileIdx(const Date& start_date) const
{
	if (cache_meteo_files.empty())
		throw InvalidArgumentException("No input files found or configured!", AT);

	//find which file we should open
	if (cache_meteo_files.size()==1) {
		return 0;
	} else {
		for (size_t idx=1; idx<cache_meteo_files.size(); idx++) {
			if (start_date>=cache_meteo_files[idx-1].first && start_date<cache_meteo_files[idx].first) {
				return --idx;
			}
		}

		//not found, we take the closest timestamp we have
		if (start_date<cache_meteo_files.front().first)
			return 0;
		else
			return cache_meteo_files.size()-1;
	}
}

342
void OshdIO::readStationData(const Date& /*date*/, std::vector<StationData>& vecStation)
343
{
344
345
346
	vecStation.clear();
	if (vecMeta.empty()) fillStationMeta();
	vecStation = vecMeta;
347
348
}

349
void OshdIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
350
                             std::vector< std::vector<MeteoData> >& vecMeteo)
351
{
352
353
354
355
	vecMeteo.clear();
	size_t file_idx = getFileIdx( dateStart );
	Date station_date( cache_meteo_files[file_idx].first );
	if (station_date<dateStart || station_date>dateEnd) return; //the requested period is NOT in the available files
356
357
358

	const size_t nr_files = cache_meteo_files.size();
	const size_t nrIDs = vecIDs.size();
359
	
360
	if (vecMeta.empty()) fillStationMeta(); //this also fills vecIdx
361
362
363
364
	vecMeteo.resize( nrIDs );
	do {
		//create empty MeteoData for the current timestep
		for (size_t jj=0; jj<nrIDs; jj++) {
365
			const MeteoData md( station_date, vecMeta[jj] );
366
367
368
369
370
			vecMeteo[jj].push_back( md );
		}
		
		//read the data and fill vecMeteo
		const std::string file_suffix( cache_meteo_files[ file_idx ].second );
371
		std::vector<double> vecData;
372
373
		for (size_t ii=0; ii<params_map.size(); ii++) {
			const MeteoData::Parameters param( params_map[ii].first );
374
			const std::string prefix( params_map[ii].second );
375
			const std::string filename( in_meteopath + "/" + prefix + "_" + file_suffix );
376
377
			vecData.resize( nrIDs, IOUtils::nodata );
			readFromFile(filename, param, station_date, vecData);
378

379
380
381
382
383
			for (size_t jj=0; jj<nrIDs; jj++)
				vecMeteo[jj].back()( param ) =  vecData[jj];
		}
		
		readSWRad(station_date, file_suffix, nrIDs, vecMeteo); //the short wave radiation is a little different...
384
385
		readPPhase(station_date, file_suffix, nrIDs, vecMeteo); //the precipitation phase is a little different...

386
387
388
389
390
		file_idx++;
		station_date = ((file_idx)<nr_files)? cache_meteo_files[file_idx].first : dateEnd+1.;
	} while (file_idx<nr_files && station_date<=dateEnd);
}

391
void OshdIO::readSWRad(const Date& station_date, const std::string& file_suffix, const size_t& nrIDs, std::vector< std::vector<MeteoData> >& vecMeteo) const
392
393
394
{
	std::vector<double> vecDir;
	vecDir.resize( nrIDs, IOUtils::nodata );
395
	const std::string filename_dir( in_meteopath + "/" + "idir" + "_" + file_suffix );
396
397
398
399
	readFromFile(filename_dir, MeteoData::ISWR, station_date, vecDir);
	
	std::vector<double> vecDiff;
	vecDiff.resize( nrIDs, IOUtils::nodata );
400
	const std::string filename_diff( in_meteopath + "/" + "idif" + "_" + file_suffix );
401
	readFromFile(filename_diff, MeteoData::ISWR, station_date, vecDiff);
402
403
404
405

	std::vector<double> vecAlbd;
	vecAlbd.resize( nrIDs, IOUtils::nodata );
	const std::string filename_albd( in_meteopath + "/" + "albd" + "_" + file_suffix );
406
	readFromFile(filename_albd, MeteoData::RSWR, station_date, vecAlbd); //We read ALBD and use it to build RSWR
407
	
408
	for (size_t jj=0; jj<nrIDs; jj++) {
409
		vecMeteo[jj].back()( MeteoData::ISWR ) =  vecDir[jj]+vecDiff[jj];
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
		vecMeteo[jj].back()( MeteoData::RSWR ) =  (vecDir[jj]+vecDiff[jj])*vecAlbd[jj];
	}
}

void OshdIO::readPPhase(const Date& station_date, const std::string& file_suffix, const size_t& nrIDs, std::vector< std::vector<MeteoData> >& vecMeteo) const
{
	std::vector<double> vecSnowLine;
	vecSnowLine.resize( nrIDs, IOUtils::nodata );
	const std::string filename_dir( in_meteopath + "/" + "snfl" + "_" + file_suffix );
	readFromFile(filename_dir, MeteoData::PSUM_PH, station_date, vecSnowLine);

	for (size_t jj=0; jj<nrIDs; jj++) {
		const double altitude = vecMeteo[jj].front().meta.getAltitude();
		if (altitude>(vecSnowLine[jj]+50.))
			vecMeteo[jj].back()( MeteoData::PSUM_PH ) = 0.;
		else if (altitude<(vecSnowLine[jj]-50.))
			vecMeteo[jj].back()( MeteoData::PSUM_PH ) = 1.;
		else
			vecMeteo[jj].back()( MeteoData::PSUM_PH ) = .5;
	}
430
431
}

432
void OshdIO::readFromFile(const std::string& filename, const MeteoData::Parameters& param, const Date& in_timestep, std::vector<double> &vecData) const
433
{
434
	if (debug) printFileStructure(filename);
435
436
437
438
439
	mat_t *matfp = Mat_Open(filename.c_str(), MAT_ACC_RDONLY);
	if ( NULL == matfp ) throw AccessException(filename, AT);

	//open the file and read some metadata
	matvar_t *matvar = Mat_VarReadInfo(matfp, "stat");
440
	if (matvar==NULL) throw NotFoundException("structure 'stat' not found in file '"+filename+"'", AT);
441
	if (matvar->class_type!=MAT_C_STRUCT) throw InvalidFormatException("The matlab file should contain 1 structure", AT);
442
	
443
444
445
446
447
448
449
450
451
452
453
454
455
456
	const std::string type( readString(filename, "type", matfp, matvar) );
	checkFieldType(param, type);

	//check that the timestep is as expected
	std::vector<double> vecTime;
	readDoubleVector(filename, "time", matfp, matvar, vecTime);
	if (vecTime.size()!=1) throw InvalidFormatException("one and only one time step must be present in the 'time' vector", AT);
	Date timestep;
	timestep.setMatlabDate( vecTime[0], in_dflt_TZ );
	if (in_timestep!=timestep) throw InvalidArgumentException("the in-file timestep and the filename time step don't match for for '"+filename+"'", AT);
	
	//check that each station is still at the same index, build the index cache if necessary
	std::vector<std::string> vecAcro;
	readStringVector(filename, "acro", matfp, matvar, vecAcro);
457
	const size_t nrIDs = vecIDs.size();
458
459
460
461
462
463
	for (size_t ii=0; ii<nrIDs; ii++) { //check that the IDs still match
		if (vecIDs[ii] != vecAcro[ vecIdx[ii] ])
			throw InvalidFormatException("station '"+vecIDs[ii]+"' is not listed in the same position as previously in file '"+filename+"'", AT);
	}
	
	//extract the data for the selected stations
464
	const std::string units( readString(filename, "dunit", matfp, matvar) );
465
466
467
468
	std::vector<double> vecRaw;
	readDoubleVector(filename, "data", matfp, matvar, vecRaw);
	if (vecAcro.size() != vecRaw.size()) throw InvalidFormatException("'acro' and 'data' arrays don't match in file '"+filename+"'", AT);
	for (size_t ii=0; ii<nrIDs; ii++)
469
		vecData[ii] = convertUnits( vecRaw[ vecIdx[ii] ], units, param);
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
	
	Mat_VarFree(matvar);
	Mat_Close(matfp);
}

void OshdIO::checkFieldType(const MeteoData::Parameters& param, const std::string& type)
{
	if (param==MeteoData::TA && type=="TA") return;
	if (param==MeteoData::RH && type=="RH") return;
	if (param==MeteoData::PSUM && type=="PREC") return;
	if (param==MeteoData::VW && type=="WS") return;
	if (param==MeteoData::DW && type=="WD") return;
	if (param==MeteoData::ILWR && type=="LWR") return;
	if (param==MeteoData::ISWR && type=="SWR") return;
	if (param==MeteoData::P && type=="other") return;
485
	if (param==MeteoData::PSUM_PH && type=="other") return;
486
	if (param==MeteoData::RSWR && type=="other") return; //this is in fact ALBD
487
488
489
490
	
	throw InvalidArgumentException("trying to read "+MeteoData::getParameterName(param)+" but found '"+type+"'", AT);
}

491
double OshdIO::convertUnits(const double& val, const std::string& units, const MeteoData::Parameters& param)
492
493
494
{
	if (units=="%") return val/100.;
	if (units=="cm") return val/100.;
495
496
497
498
	if (units=="mm") {
		if (param==MeteoData::PSUM) return val;
		else return val/1000.;
	}
499
500
501
502
503
504
505
506
507
	if (units=="\xB0\x43") return val+Cst::t_water_freezing_pt; //ISO-8859-1 hex for '°C'
// 	//usually skip these extra tests
// 	if (units=="\xB0") return val; //ISO-8859-1 hex for '°'
// 	if (units.empty()) return val;
// 	if (units=="Pa") return val;
// 	if (units=="W/m2") return val;
// 	if (units=="m/s") return val;
// 	else 
// 		throw IOException("Unknown units '"+units+"'", AT);
508
509
510
511
	
	return val;
}

512
513
void OshdIO::fillStationMeta()
{
514
	vecMeta.resize( vecIDs.size(), StationData() );
515
516
	mat_t *matfp = Mat_Open(in_metafile.c_str(), MAT_ACC_RDONLY);
	if ( NULL == matfp ) throw AccessException(in_metafile, AT);
517
518

	matvar_t *matvar = Mat_VarReadInfo(matfp, "statlist");
519
	if (matvar==NULL) throw NotFoundException("structure 'statlist' not found in file '"+in_metafile+"'", AT);
520
521
	
	std::vector<std::string> vecAcro;
522
	readStringVector(in_metafile, "acro", matfp, matvar, vecAcro);
523
524
	
	std::vector<std::string> vecNames;
525
	readStringVector(in_metafile, "name", matfp, matvar, vecNames);
526
527
	
	std::vector<double> easting, northing, altitude;
528
529
530
	readDoubleVector(in_metafile, "x", matfp, matvar, easting);
	readDoubleVector(in_metafile, "y", matfp, matvar, northing);
	readDoubleVector(in_metafile, "z", matfp, matvar, altitude);
531
532
533
534
	
	Mat_VarFree(matvar);
	Mat_Close(matfp);
	
535
536
537
538
539
540
	if (debug) {
		for (size_t ii=0; ii<vecAcro.size(); ii++) 
			std::cout << std::setw(8) << vecAcro[ii] << std::setw(40) << vecNames[ii] << std::setw(8) << easting[ii] << std::setw(8) << northing[ii] << std::setw(8) << altitude[ii] << "\n";
		std::cout << endl;
	}
	
541
542
543
544
545
546
547
548
549
550
	buildVecIdx(vecAcro);
	for (size_t ii=0; ii<vecIdx.size(); ii++) {
		Coords position("CH1903", "");
		const size_t idx = vecIdx[ii];
		position.setXY(easting[idx], northing[idx], altitude[idx]);
		const StationData sd(position, vecAcro[idx], vecNames[idx]);
		vecMeta[ii] = sd;
	}
}

551
552
553
554
555
void OshdIO::buildVecIdx(const std::vector<std::string>& vecAcro)
{
	const size_t nrIDs = vecIDs.size();
	if (nrIDs==0)
		throw InvalidArgumentException("Please provide at least one station ID to read!", AT);
556
	vecIdx.resize( nrIDs, 0 );
557
558
559
560
	
	for (size_t ii=0; ii<nrIDs; ii++) {
		bool found = false;
		for (size_t jj=0; jj<vecAcro.size(); jj++) {
561
			if (vecIDs[ii]==vecAcro[jj]) {
562
563
564
565
566
567
				vecIdx[ii] = jj;
				found = true;
				break;
			}
		}
		if (!found) 
568
			throw NotFoundException("station ID '"+vecIDs[ii]+"' could not be found in the provided data", AT);
569
570
571
572
	}
}

} //namespace