WSL/SLF GitLab Repository

OshdIO.cc 26 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, "ilwc") );
281
	params_map.push_back( std::make_pair(MeteoData::P, "pair") );
282
	params_map.push_back( std::make_pair(MeteoData::PSUM, "prec") ); //in mm/ts
283
	params_map.push_back( std::make_pair(MeteoData::RH, "rcor") ); //old: rhum
284
	params_map.push_back( std::make_pair(MeteoData::TA, "tcor") ); //old:tair
285
	params_map.push_back( std::make_pair(MeteoData::VW, "wcor") ); //old: wind
286
287
288
	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< struct file_index > &meteo_files)
290
291
{
	meteo_files.clear();
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324

	std::list<std::string> dirlist( FileUtils::readDirectory(meteopath_in, "prec", is_recursive) ); //we consider that if we have found one parameter, the others are also there

	std::map<std::string, size_t> mapIdx; //make sure each timestamp only appears once, ie remove duplicates
	for (std::list<std::string>::iterator it = dirlist.begin(); it != dirlist.end(); ++it) {
		const std::string file_and_path( *it );
		const std::string filename( FileUtils::getFilename(file_and_path) );

		//we need to split the file name into its components: parameter, date, run_date
		const std::string::size_type pos_param = filename.find('_');
		if (pos_param==string::npos) continue;
		const std::string::size_type date_start = filename.find_first_of("0123456789");
		if (date_start==string::npos) continue;
		const std::string::size_type date_end = filename.find('_', date_start);
		if (date_end==string::npos) continue;
		const std::string::size_type rundate_start = filename.rfind('_');
		if (rundate_start==string::npos) continue;
		const std::string::size_type rundate_end = filename.rfind('.');
		if (rundate_end==string::npos) continue;

		const std::string date_str( filename.substr(date_start, date_end-date_start) );
		const std::string run_date( filename.substr(rundate_start+1, rundate_end-rundate_start-1) );

		//do we already have an entry for this date?
		size_t idx = IOUtils::npos;
		const std::map<std::string, size_t>::const_iterator it_map = mapIdx.find( date_str );
		if (it_map!=mapIdx.end()) {
			idx = it_map->second;
			if (meteo_files[idx].run_date>run_date) continue;
		}

		//we don't have an entry or it is too old -> create new entry / replace existing one
		const std::string path( FileUtils::getPath(file_and_path) );
325
		Date date;
326
327
328
329
330
331
332
333
334
		IOUtils::convertString(date, date_str, in_dflt_TZ);
		const file_index elem(date, path, filename.substr(pos_param), run_date);
		if (idx==IOUtils::npos) {
			meteo_files.push_back( elem );
			mapIdx[ date_str ] = meteo_files.size()-1;
		} else {
			meteo_files[ idx] = elem;
			mapIdx[ date_str ] = idx;
		}
335
	}
336
337

	std::sort(meteo_files.begin(), meteo_files.end());
338
339
340
341
342
343
344
345
346
347
348
349
}

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++) {
350
			if (start_date>=cache_meteo_files[idx-1].date && start_date<cache_meteo_files[idx].date) {
351
352
353
354
355
				return --idx;
			}
		}

		//not found, we take the closest timestamp we have
356
		if (start_date<cache_meteo_files.front().date)
357
358
359
360
361
362
			return 0;
		else
			return cache_meteo_files.size()-1;
	}
}

363
void OshdIO::readStationData(const Date& /*date*/, std::vector<StationData>& vecStation)
364
{
365
366
367
	vecStation.clear();
	if (vecMeta.empty()) fillStationMeta();
	vecStation = vecMeta;
368
369
}

370
void OshdIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
371
                             std::vector< std::vector<MeteoData> >& vecMeteo)
372
{
373
374
	vecMeteo.clear();
	size_t file_idx = getFileIdx( dateStart );
375
	Date station_date( cache_meteo_files[file_idx].date );
376
	if (station_date<dateStart || station_date>dateEnd) return; //the requested period is NOT in the available files
377
378
379

	const size_t nr_files = cache_meteo_files.size();
	const size_t nrIDs = vecIDs.size();
380
	
381
	if (vecMeta.empty()) fillStationMeta(); //this also fills vecIdx
382
383
384
385
	vecMeteo.resize( nrIDs );
	do {
		//create empty MeteoData for the current timestep
		for (size_t jj=0; jj<nrIDs; jj++) {
386
			const MeteoData md( station_date, vecMeta[jj] );
387
388
389
390
			vecMeteo[jj].push_back( md );
		}
		
		//read the data and fill vecMeteo
391
392
		const std::string file_suffix( cache_meteo_files[ file_idx ].file_suffix );
		const std::string path( in_meteopath + "/" + cache_meteo_files[ file_idx ].path );
393
		std::vector<double> vecData;
394
395
		for (size_t ii=0; ii<params_map.size(); ii++) {
			const MeteoData::Parameters param( params_map[ii].first );
396
			const std::string filename( path + "/" + params_map[ii].second + file_suffix );
397
398
			vecData.resize( nrIDs, IOUtils::nodata );
			readFromFile(filename, param, station_date, vecData);
399

400
401
402
403
			for (size_t jj=0; jj<nrIDs; jj++)
				vecMeteo[jj].back()( param ) =  vecData[jj];
		}
		
404
405
		readSWRad(station_date, path, file_suffix, nrIDs, vecMeteo); //the short wave radiation is a little different...
		readPPhase(station_date, path, file_suffix, nrIDs, vecMeteo); //the precipitation phase is a little different...
406

407
		file_idx++;
408
		station_date = ((file_idx)<nr_files)? cache_meteo_files[file_idx].date : dateEnd+1.;
409
410
411
	} while (file_idx<nr_files && station_date<=dateEnd);
}

412
void OshdIO::readSWRad(const Date& station_date, const std::string& path, const std::string& file_suffix, const size_t& nrIDs, std::vector< std::vector<MeteoData> >& vecMeteo) const
413
414
415
{
	std::vector<double> vecDir;
	vecDir.resize( nrIDs, IOUtils::nodata );
416
	const std::string filename_dir( path + "/" + "idrc" + file_suffix );
417
418
419
420
	readFromFile(filename_dir, MeteoData::ISWR, station_date, vecDir);
	
	std::vector<double> vecDiff;
	vecDiff.resize( nrIDs, IOUtils::nodata );
421
	const std::string filename_diff( path + "/" + "idfc" + file_suffix );
422
	readFromFile(filename_diff, MeteoData::ISWR, station_date, vecDiff);
423
424

	std::vector<double> vecAlbd;
425
	const std::string filename_albd( path + "/" + "albd" + file_suffix );
426
427
428
429
	if (FileUtils::fileExists(filename_albd)) {
		vecAlbd.resize( nrIDs, IOUtils::nodata );
		readFromFile(filename_albd, MeteoData::RSWR, station_date, vecAlbd); //We read ALBD and use it to build RSWR
	}
430
	
431
	const double albedo = !vecAlbd.empty();
432
	for (size_t jj=0; jj<nrIDs; jj++) {
433
		vecMeteo[jj].back()( MeteoData::ISWR ) =  vecDir[jj]+vecDiff[jj];
434
		if (albedo) vecMeteo[jj].back()( MeteoData::RSWR ) =  (vecDir[jj]+vecDiff[jj])*vecAlbd[jj];
435
436
437
	}
}

438
void OshdIO::readPPhase(const Date& station_date, const std::string& path, const std::string& file_suffix, const size_t& nrIDs, std::vector< std::vector<MeteoData> >& vecMeteo) const
439
{
440
	const std::string filename( path + "/" + "snfl" + file_suffix );
441

442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
	if (FileUtils::fileExists(filename)) {
		const double half_elevation_band = 50.;  //we consider that there are mixed precip in the elevation range snow_line ± half_elevation_band
		std::vector<double> vecSnowLine;
		vecSnowLine.resize( nrIDs, IOUtils::nodata );
		readFromFile(filename, 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]+half_elevation_band))
				vecMeteo[jj].back()( MeteoData::PSUM_PH ) = 0.;
			else if (altitude<(vecSnowLine[jj]-half_elevation_band))
				vecMeteo[jj].back()( MeteoData::PSUM_PH ) = 1.;
			else
				vecMeteo[jj].back()( MeteoData::PSUM_PH ) = .5;
		}
457
	}
458
459
}

460
void OshdIO::readFromFile(const std::string& filename, const MeteoData::Parameters& param, const Date& in_timestep, std::vector<double> &vecData) const
461
{
462
	if (debug) printFileStructure(filename);
463
464
465
466
467
	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");
468
	if (matvar==NULL) throw NotFoundException("structure 'stat' not found in file '"+filename+"'", AT);
469
	if (matvar->class_type!=MAT_C_STRUCT) throw InvalidFormatException("The matlab file should contain 1 structure", AT);
470
	
471
472
473
474
475
476
477
478
479
480
481
482
483
484
	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);
485
	const size_t nrIDs = vecIDs.size();
486
487
488
489
490
491
	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
492
	const std::string units( readString(filename, "dunit", matfp, matvar) );
493
494
495
496
	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++)
497
		vecData[ii] = convertUnits( vecRaw[ vecIdx[ii] ], units, param);
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
	
	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;
513
	if (param==MeteoData::PSUM_PH && type=="other") return;
514
	if (param==MeteoData::RSWR && type=="other") return; //this is in fact ALBD
515
516
517
518
	
	throw InvalidArgumentException("trying to read "+MeteoData::getParameterName(param)+" but found '"+type+"'", AT);
}

519
double OshdIO::convertUnits(const double& val, const std::string& units, const MeteoData::Parameters& param)
520
521
522
{
	if (units=="%") return val/100.;
	if (units=="cm") return val/100.;
523
524
525
526
	if (units=="mm") {
		if (param==MeteoData::PSUM) return val;
		else return val/1000.;
	}
527
528
529
530
531
532
533
534
535
	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);
536
537
538
539
	
	return val;
}

540
541
void OshdIO::fillStationMeta()
{
542
	vecMeta.resize( vecIDs.size(), StationData() );
543
544
	mat_t *matfp = Mat_Open(in_metafile.c_str(), MAT_ACC_RDONLY);
	if ( NULL == matfp ) throw AccessException(in_metafile, AT);
545
546

	matvar_t *matvar = Mat_VarReadInfo(matfp, "statlist");
547
	if (matvar==NULL) throw NotFoundException("structure 'statlist' not found in file '"+in_metafile+"'", AT);
548
549
	
	std::vector<std::string> vecAcro;
550
	readStringVector(in_metafile, "acro", matfp, matvar, vecAcro);
551
552
	
	std::vector<std::string> vecNames;
553
	readStringVector(in_metafile, "name", matfp, matvar, vecNames);
554
555
	
	std::vector<double> easting, northing, altitude;
556
557
558
	readDoubleVector(in_metafile, "x", matfp, matvar, easting);
	readDoubleVector(in_metafile, "y", matfp, matvar, northing);
	readDoubleVector(in_metafile, "z", matfp, matvar, altitude);
559
560
561
562
	
	Mat_VarFree(matvar);
	Mat_Close(matfp);
	
563
564
565
566
567
568
	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;
	}
	
569
570
571
	buildVecIdx(vecAcro);
	for (size_t ii=0; ii<vecIdx.size(); ii++) {
		const size_t idx = vecIdx[ii];
572
		Coords position("CH1903", "");
573
		position.setXY(easting[idx], northing[idx], altitude[idx]);
574
575
576
577
578
579
580
581
		std::string name( vecNames[idx] );

		//if the network name has been appended, remove it. We also remove spaces, just in case
		const size_t netz_pos = name.find(" (");
		if (netz_pos!=std::string::npos) name.erase(netz_pos);
		std::replace( name.begin(), name.end(), ' ', '_');

		const StationData sd(position, vecAcro[idx], name);
582
583
584
585
		vecMeta[ii] = sd;
	}
}

586
587
588
589
590
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);
591
	vecIdx.resize( nrIDs, 0 );
592
593
594
595
	
	for (size_t ii=0; ii<nrIDs; ii++) {
		bool found = false;
		for (size_t jj=0; jj<vecAcro.size(); jj++) {
596
			if (vecIDs[ii]==vecAcro[jj]) {
597
598
599
600
601
602
				vecIdx[ii] = jj;
				found = true;
				break;
			}
		}
		if (!found) 
603
			throw NotFoundException("station ID '"+vecIDs[ii]+"' could not be found in the provided data", AT);
604
605
606
607
	}
}

} //namespace