WSL/SLF GitLab Repository

SMETIO.cc 29.5 KB
Newer Older
1
/***********************************************************************************/
2
/*  Copyright 2010 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/***********************************************************************************/
/* 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 "SMETIO.h"
19
#include <meteoio/IOUtils.h>
20
21
22
23
24

using namespace std;

namespace mio {
/**
25
 * @page smetio SMET
26
 * @section template_format Format
27
 * The Station METeo data files is a station centered, ascii file format that has been designed with flexibility and ease of use in mind. Please refer to its <a href="../SMET_specifications.pdf">official format specification</a> for more information.
28
29
 *
 * @section template_units Units
30
 * All units are MKSA, the only exception being the precipitations that are in mm/h. It is however possible to use  multipliers and offsets (but they must be specified in the file header).
31
32
33
 *
 * @section template_keywords Keywords
 * This plugin uses the following keywords:
34
35
 * - STATION#: input filename (in METEOPATH). As many meteofiles as needed may be specified
 * - METEOPATH: meteo files directory where to read/write the meteofiles; [Input] and [Output] sections
36
37
38
39
40
 * - METEOPARAM: output file format options (ASCII or BINARY that might be followed by GZIP)
 *
 * Example:
 * @code
 * [Input]
41
42
43
44
 * METEOPATH = ./input
 * STATION1 = uppper_station.smet
 * STATION2 = lower_station.smet
 * STATION3 = outlet_station.smet
45
46
47
48
 * [Output]
 * METEOPATH = ./output
 * METEOPARAM = ASCII GZIP
 * @endcode
49
50
 */

51
52
//TODO: keep a pointer to last read position, so we can fseek for the next read?

53
const std::string SMETIO::smet_version = "1.0";
54
const unsigned int SMETIO::buffer_reserve = 23*24*2; //kind of average size of a buffer for optimizing vectors
55
56
map<string, MeteoData::Parameters> SMETIO::mapParameterByName;
const bool SMETIO::__init = SMETIO::initStaticData();
57

58
bool SMETIO::initStaticData()
59
60
61
62
63
64
65
66
{
	for (unsigned int ii=0; ii<MeteoData::nrOfParameters; ii++){
		mapParameterByName[MeteoData::getParameterName(ii)] = MeteoData::Parameters(ii);
	}

	return true;
}

67
double& SMETIO::getParameter(const std::string& columnName, MeteoData& md)
68
69
70
71
72
73
{//HACK: the whole name mapping is a big hack. Replace with proper mapping!!!
	if(columnName=="OSWR") {
		MeteoData::Parameters paramindex = mapParameterByName["RSWR"];
		return md.param(paramindex);
	}
	if(columnName=="PSUM") {
74
		MeteoData::Parameters paramindex = mapParameterByName["HNW"];
75
76
77
		return md.param(paramindex);
	}

78
79
80
81
	MeteoData::Parameters paramindex = mapParameterByName[columnName];
	return md.param(paramindex);
}

82
void SMETIO::checkColumnNames(const std::vector<std::string>& vecColumns, const bool& locationInHeader)
83
{
84
85
86
87
	/*
	 * This function checks whether the sequence of keywords specified in the 
	 * [HEADER] section (key 'fields') is valid
	 */
88
89
90
	vector<unsigned int> paramcounter = vector<unsigned int>(MeteoData::nrOfParameters, 0);

	for (unsigned int ii=0; ii<vecColumns.size(); ii++){
91
92
93
94
95
96
97
98
		std::string column = vecColumns[ii];

		//column names mapping
		if(column=="OSWR") column="RSWR";
		if(column=="PSUM") column="HNW";

		if ((column == "timestamp") || (column == "longitude")
		    || (column == "latitude") || (column == "altitude")){
99
100
			//everything ok
		} else {
101
			map<string, MeteoData::Parameters>::iterator it = mapParameterByName.find(column);
102
103
104
			if (it == mapParameterByName.end())
				throw InvalidFormatException("Key 'fields' specified in [HEADER] section contains invalid names", AT);
			
105
			paramcounter.at(mapParameterByName[column])++;
106
107
		}
	}
108
109
	
	//Check for multiple usages of parameters
110
111
112
113
114
	for (unsigned int ii=0; ii<paramcounter.size(); ii++){
		if (paramcounter[ii] > 1)
			throw InvalidFormatException("In 'fields': Multiple use of " + MeteoData::getParameterName(ii), AT);
	}
	
115
116
	//If there is no location information in the [HEADER] section, then
	//location information must be part of fields
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
	if (!locationInHeader){
		unsigned int latcounter = 0, loncounter=0, altcounter=0;
		for (unsigned int ii=0; ii<vecColumns.size(); ii++){
			if (vecColumns[ii] == "longitude")
				loncounter++;
			else if (vecColumns[ii] == "latitude")
				latcounter++;
			else if (vecColumns[ii] == "altitude")
				altcounter++;
		}

		if ((latcounter != loncounter) && (loncounter != altcounter) && (altcounter != 1))
			throw InvalidFormatException("Key 'fields' must contain 'latitude', 'longitude' and 'altitude'", AT);
	}
}

//END STATIC SECTION


136
SMETIO::SMETIO(void (*delObj)(void*), const Config& i_cfg) : IOInterface(delObj), cfg(i_cfg)
137
138
139
140
{
	parseInputOutputSection();
}

141
SMETIO::SMETIO(const std::string& configfile) : IOInterface(NULL), cfg(configfile)
142
143
144
145
{
	parseInputOutputSection();
}

146
SMETIO::SMETIO(const Config& cfgreader) : IOInterface(NULL), cfg(cfgreader)
147
148
149
150
{
	parseInputOutputSection();
}

151
SMETIO::~SMETIO() throw()
152
153
154
155
{
	cleanup();
}

156
void SMETIO::cleanup() throw()
157
{
158
159
160
	//clear ios flags
	fout << resetiosflags(ios_base::fixed | ios_base::left);

161
162
163
164
165
166
167
168
	if (fin.is_open()) {//close fin if open
		fin.close();
	}
	if (fout.is_open()) {//close fout if open
		fout.close();
	}		
}

169
void SMETIO::read2DGrid(Grid2DObject& /*grid_out*/, const std::string& /*_name*/)
170
171
172
173
174
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

175
void SMETIO::readDEM(DEMObject& /*dem_out*/)
176
177
178
179
180
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

181
void SMETIO::readLanduse(Grid2DObject& /*landuse_out*/)
182
183
184
185
186
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

187
void SMETIO::readAssimilationData(const Date& /*date_in*/, Grid2DObject& /*da_out*/)
188
189
190
191
192
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

193
194
void SMETIO::readStationData(const Date&, std::vector<StationData>& vecStation)
{//big HACK: this is a barbaric code duplication!! Plus it should support coordinates in the data
195
//ie: it should use the given date! (and TZ)
196
	vecStation.clear();
197
	vecStation.reserve(nr_stations);
198
199

	//Now loop through all requested stations, open the respective files and parse them
200
	for (unsigned int ii=0; ii<nr_stations; ii++){
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
		bool isAscii = true;
		string filename = vecFiles.at(ii); //filename of current station
		
		if (!IOUtils::fileExists(filename))
			throw FileNotFoundException(filename, AT);

		fin.clear();
		fin.open (filename.c_str(), ios::in);
		if (fin.fail())
			throw FileAccessException(filename, AT);

		char eoln = IOUtils::getEoln(fin); //get the end of line character for the file

		//Go through file, save key value pairs
		string line="";
		std::vector<std::string> tmpvec, vecDataSequence;
217
		std::vector<double> vecUnitsOffset, vecUnitsMultiplier;
218
219
220
221
222
223
224
225
226
227
228
229
		double timezone = 0.0;
		bool locationInHeader = false;
		StationData sd;

		try {
			//1. Read signature
			getline(fin, line, eoln); //read complete signature line
			IOUtils::stripComments(line);
			IOUtils::readLineToVec(line, tmpvec);
			checkSignature(tmpvec, filename, isAscii);

			//2. Read Header
230
			readHeader(eoln, filename, locationInHeader, timezone, sd, vecDataSequence, vecUnitsOffset, vecUnitsMultiplier);
231
232
233
234
235
236
237
			cleanup();
		} catch(std::exception& e) {
			cleanup();
			throw;
		}
		vecStation.push_back(sd);
	}
238
239
}

240
void SMETIO::parseInputOutputSection()
241
{
242
	//default timezones
243
	in_dflt_TZ = out_dflt_TZ = IOUtils::nodata;
244
245
246
247
	cfg.getValue("TZ","Input",in_dflt_TZ,Config::nothrow);
	cfg.getValue("TZ","Output",out_dflt_TZ,Config::nothrow);

	// Parse the [Input] and [Output] sections within Config object cfg
248
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
249
250

	//Parse input section: extract number of files to read and store filenames in vecFiles
251
	std::string inpath="", in_meteo="";
252
253
	cfg.getValue("METEO", "Input", in_meteo, Config::nothrow);
	if (in_meteo == "SMET") { //keep it synchronized with IOHandler.cc for plugin mapping!!
254
		cfg.getValue("METEOPATH", "Input", inpath);
255
256
		unsigned int counter = 1;
		string filename = "";
257

258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
		do {
			stringstream ss;
			filename = "";
			
			ss << "STATION" << counter;
			cfg.getValue(ss.str(), "Input", filename, Config::nothrow);
			
			if (filename != ""){
				stringstream file_and_path;
				file_and_path << inpath << "/" << filename;
				if (!IOUtils::validFileName(file_and_path.str())) //Check whether filename is valid
					throw InvalidFileNameException(file_and_path.str(), AT);
				
				vecFiles.push_back(file_and_path.str());
			}
			counter++;
		} while (filename != "");
		
		nr_stations = counter - 1;
	}
278

279
	//Parse output section: extract info on whether to write ASCII or BINARY format, gzipped or not
280
281
282
283
284
	outpath = "";
	outputIsAscii = true; 
	outputIsGzipped = false;

	vector<string> vecArgs;
285
286
	cfg.getValue("METEOPATH", "Output", outpath, Config::nothrow);
	cfg.getValue("METEOPARAM", "Output", vecArgs, Config::nothrow); //"ASCII|BINARY GZIP"
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309

	if (outpath == "")
		return;

	if (vecArgs.size() == 0)
		vecArgs.push_back("ASCII");

	if (vecArgs.size() > 2)
		throw InvalidFormatException("Too many values for key METEOPARAM", AT);

	if (vecArgs[0] == "BINARY")
		outputIsAscii = false;
	else if (vecArgs[0] == "ASCII")
		outputIsAscii = true;
	else 
		throw InvalidFormatException("The first value for key METEOPARAM may only be ASCII or BINARY", AT);

	if (vecArgs.size() == 2){
		if (vecArgs[1] != "GZIP")
			throw InvalidFormatException("The second value for key METEOPARAM may only be GZIP", AT);

		outputIsGzipped = true;
	}
310

311
312
}

313
void SMETIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
314
315
                           std::vector< std::vector<MeteoData> >& vecMeteo,
                           const unsigned int& stationindex)
316
{
317
	//Make sure that vecMeteo have the correct dimension and stationindex is valid
318
319
	unsigned int startindex=0, endindex=vecFiles.size();	
	if (stationindex != IOUtils::npos){
320
		if ((stationindex < vecFiles.size()) || (stationindex < vecMeteo.size())){
321
322
323
324
325
326
327
328
			startindex = stationindex;
			endindex = stationindex+1;
		} else {
			throw IndexOutOfBoundsException("Invalid stationindex", AT);
		} 

		vecMeteo[stationindex].clear();
	} else {
329
		vecMeteo.clear();
330
		vecMeteo = vector< vector<MeteoData> >(vecFiles.size());
331
		vecMeteo.reserve(nr_stations);
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
	}

	//Now loop through all requested stations, open the respective files and parse them
	for (unsigned int ii=startindex; ii<endindex; ii++){
		bool isAscii = true;
		string filename = vecFiles.at(ii); //filename of current station
		
		if (!IOUtils::fileExists(filename))
			throw FileNotFoundException(filename, AT);

		fin.clear();
		fin.open (filename.c_str(), ios::in);
		if (fin.fail())
			throw FileAccessException(filename, AT);

		char eoln = IOUtils::getEoln(fin); //get the end of line character for the file

		//Go through file, save key value pairs
		string line="";
		std::vector<std::string> tmpvec, vecDataSequence;
352
		std::vector<double> vecUnitsOffset, vecUnitsMultiplier;
353
354
355
356
357
358
359
360
361
362
363
364
		double timezone = 0.0;
		bool locationInHeader = false;
		StationData sd;

		try {
			//1. Read signature
			getline(fin, line, eoln); //read complete signature line
			IOUtils::stripComments(line);
			IOUtils::readLineToVec(line, tmpvec);
			checkSignature(tmpvec, filename, isAscii);

			//2. Read Header
365
			readHeader(eoln, filename, locationInHeader, timezone, sd, vecDataSequence, vecUnitsOffset, vecUnitsMultiplier);
366
			SMETIO::checkColumnNames(vecDataSequence, locationInHeader);
367

368
369
370
371
372
373
374
375
			//now, timezone MUST contain something valid
			if(timezone==IOUtils::nodata) {
				stringstream ss;
				ss << "No timezone information available for file " << filename;
				ss << " either in file header or io.ini";
				throw InvalidFormatException(ss.str(), AT);
			}

376
377
			//3. Read DATA
			if (isAscii){
378
				readDataAscii(eoln, filename, timezone, sd, vecDataSequence, vecUnitsOffset,
379
			                      vecUnitsMultiplier, dateStart, dateEnd, vecMeteo[ii]);
380
381
382
383
384
385
386
			} else {
				streampos currpos = fin.tellg();
				fin.close();
				fin.open (filename.c_str(), ios::in|ios::binary);				
				if (fin.fail()) throw FileAccessException(filename, AT);
				fin.seekg(currpos); //jump to binary data section

387
				readDataBinary(eoln, filename, timezone, sd, vecDataSequence, vecUnitsOffset,
388
				               vecUnitsMultiplier, dateStart, dateEnd, vecMeteo[ii]);
389
390
391
392
393
394
395
396
397
			}
			cleanup();
		} catch(std::exception& e) {
			cleanup();
			throw;
		}
	}
}

398
void SMETIO::readDataBinary(const char&, const std::string&, const double& timezone,
399
                            const StationData& sd, const std::vector<std::string>& vecDataSequence,
400
                            const std::vector<double>& vecUnitsOffset, std::vector<double>& vecUnitsMultiplier,
401
                            const Date& dateStart, const Date& dateEnd, std::vector<MeteoData>& vecMeteo)
402
403
404
{
	const unsigned int nrOfColumns = vecDataSequence.size();

405
406
	vecMeteo.reserve(buffer_reserve);

407
408
409
410
411
412
413
414
	while (!fin.eof()){
		MeteoData md;
		StationData tmpsd = sd;
		double lat=IOUtils::nodata, lon=IOUtils::nodata, alt=IOUtils::nodata;
		unsigned int poscounter = 0;

		for (unsigned int ii=0; ii<nrOfColumns; ii++){
			if (vecDataSequence[ii] == "timestamp"){
415
416
417
				double tmpval;
				fin.read(reinterpret_cast < char * > (&tmpval), sizeof(double));

418
				md.date.setDate(tmpval, timezone);
419
420
421
422

				if (md.date < dateStart)
					continue;
				if (md.date > dateEnd)
423
					return;
424
			} else {
425
426
427
428
429
430
431
432
433
434
435
436
437
438
				float tmpval;
				fin.read(reinterpret_cast < char * > (&tmpval), sizeof(float));			
				double val = (double)tmpval;

				if (vecDataSequence[ii] == "latitude"){
					lat = val;
					poscounter++;
				} else if (vecDataSequence[ii] == "longitude"){
					lon = val;
					poscounter++;
				} else if (vecDataSequence[ii] == "altitude"){
					alt = val;
					poscounter++;
				} else {
439
					SMETIO::getParameter(vecDataSequence[ii], md) = (val + vecUnitsOffset[ii]) * vecUnitsMultiplier[ii];
440
				}
441
442
443
444
445
446
447
448
449
450
451
			}
		}

		char c;
		fin.read(&c, sizeof(char));
		if (c != '\n')
			throw InvalidFormatException("Corrupted data in section [DATA]", AT);

		if (poscounter == 3)
			tmpsd.position.setLatLon(lat, lon, alt);

452
		if (md.date >= dateStart){
453
			md.meta = tmpsd;
454
455
			vecMeteo.push_back(md);
		}
456
457
458
	}	
}

459
void SMETIO::readDataAscii(const char& eoln, const std::string& filename, const double& timezone,
460
                           const StationData& sd, const std::vector<std::string>& vecDataSequence,
461
                           const std::vector<double>& vecUnitsOffset, std::vector<double>& vecUnitsMultiplier,
462
                           const Date& dateStart, const Date& dateEnd, std::vector<MeteoData>& vecMeteo)
463
464
465
466
467
{
	string line = "";
	vector<string> tmpvec;
	const unsigned int nrOfColumns = vecDataSequence.size();

468
469
470
	tmpvec.reserve(nrOfColumns);
	vecMeteo.reserve(buffer_reserve);

471
	while (!fin.eof()){
472
		//HACK nodata mapping is NOT done!!!!!!
473
		//something like lat = IOUtils::standardizeNodata(lat, plugin_nodata); should be done
474
475
476
477
478
479
480
481
		getline(fin, line, eoln);
		IOUtils::stripComments(line);
		IOUtils::trim(line);
		if (line == "") continue; //Pure comment lines and empty lines are ignored

		unsigned int ncols = IOUtils::readLineToVec(line, tmpvec);

		if (ncols != nrOfColumns)
482
			throw InvalidFormatException("In "+ filename + ": Invalid amount of data in data line "+line, AT);
483
484
485
486
487
488
489

		MeteoData md;
		StationData tmpsd = sd;
		double lat, lon, alt;
		unsigned int poscounter = 0;
		for (unsigned int ii=0; ii<nrOfColumns; ii++){
			if (vecDataSequence[ii] == "timestamp"){
490
				if (!IOUtils::convertString(md.date, tmpvec[ii], timezone, std::dec))
491
					throw InvalidFormatException("In "+filename+": Timestamp "+tmpvec[ii]+" invalid in data line", AT);
492
493
494
				if (md.date < dateStart)
					continue;
				if (md.date > dateEnd)
495
					return;
496
497
498
499
500
501
502
503
504
505
506
507
508
509

			} else if (vecDataSequence[ii] == "latitude"){
				if (!IOUtils::convertString(lat, tmpvec[ii])) 
					throw InvalidFormatException("In "+filename+": Latitude invalid", AT);
				poscounter++;
			} else if (vecDataSequence[ii] == "longitude"){
				if (!IOUtils::convertString(lon, tmpvec[ii])) 
					throw InvalidFormatException("In "+filename+": Longitude invalid", AT);
				poscounter++;
			} else if (vecDataSequence[ii] == "altitude"){
				if (!IOUtils::convertString(alt, tmpvec[ii])) 
					throw InvalidFormatException("In "+filename+": Altitude invalid", AT);
				poscounter++;
			} else {
510
511
				double val;
				if (!IOUtils::convertString(val, tmpvec[ii]))
512
					throw InvalidFormatException("In "+filename+": Invalid value for param", AT);
513
				SMETIO::getParameter(vecDataSequence[ii], md) = (val + vecUnitsOffset[ii]) * vecUnitsMultiplier[ii];
514
515
516
517
518
519
			}
		}

		if (poscounter == 3)
			tmpsd.position.setLatLon(lat, lon, alt);

520
		if (md.date >= dateStart){
521
			md.meta = tmpsd;
522
523
			vecMeteo.push_back(md);
		}
524
525
526
	}
}

527
void SMETIO::readHeader(const char& eoln, const std::string& filename, bool& locationInHeader,
528
529
                        double& timezone, StationData& sd, std::vector<std::string>& vecDataSequence,
                        std::vector<double>& vecUnitsOffset, std::vector<double>& vecUnitsMultiplier)
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
{
	string line="";
	while (!fin.eof() && (fin.peek() != '['))
		getline(fin, line, eoln);
	
	getline(fin, line, eoln);
	IOUtils::stripComments(line);
	IOUtils::trim(line);
	IOUtils::toUpper(line);

	if (line != "[HEADER]")
		throw InvalidFormatException("Section " + line + " in "+ filename + " invalid", AT);

	std::map<std::string,std::string> mapHeader;
	while (!fin.eof() && (fin.peek() != '[')){
		getline(fin, line, eoln);

		IOUtils::stripComments(line);
		IOUtils::trim(line);

550
		if (line != "") {
551
552
			if (!IOUtils::readKeyValuePair(line, "=", mapHeader))
				throw InvalidFormatException("Invalid key value pair in section [Header]", AT);
553
		}
554
555
556
557
558
	}

	//Now extract info from mapHeader
	IOUtils::getValueForKey(mapHeader, "station_id", sd.stationID);
	IOUtils::getValueForKey(mapHeader, "station_name", sd.stationName, IOUtils::nothrow);
559
	timezone = in_dflt_TZ;
560
	IOUtils::getValueForKey(mapHeader, "tz", timezone, IOUtils::nothrow);
561
	sd.position.setProj(coordin, coordinparam); //set the default projection from config file
562

563
564
565
	//trying to read easting/northing
	double easting=IOUtils::nodata, northing=IOUtils::nodata, alt=IOUtils::nodata;
	short int epsg=IOUtils::snodata;
566
567
568
569
	IOUtils::getValueForKey(mapHeader, "epsg", epsg, IOUtils::nothrow);
	if(epsg!=IOUtils::snodata) {
		sd.position.setEPSG(epsg);
	}
570
	
571
	IOUtils::getValueForKey(mapHeader, "easting", easting, IOUtils::nothrow);
572
	if (easting != IOUtils::nodata){
573
574
575
576
577
578
579
580
581
582
		IOUtils::getValueForKey(mapHeader, "northing", northing);
		IOUtils::getValueForKey(mapHeader, "altitude", alt);
		sd.position.setXY(easting, northing, alt);
		locationInHeader = true;
	} else {
		locationInHeader = false;
	}

	//now trying to read lat/long (second, so that it has precedence over east/north coordinates)
	double lat=IOUtils::nodata, lon=IOUtils::nodata;
583
	IOUtils::getValueForKey(mapHeader, "latitude", lat, IOUtils::nothrow);
584
	if (lat != IOUtils::nodata){ 
585
586
587
588
589
590
591
		IOUtils::getValueForKey(mapHeader, "longitude", lon);
		IOUtils::getValueForKey(mapHeader, "altitude", alt);
		sd.position.setLatLon(lat, lon, alt);
		locationInHeader = true;
	}

	IOUtils::getValueForKey(mapHeader, "fields", vecDataSequence);
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625

	//trying to read units offsets
	std::vector<std::string> vecTmp;
	IOUtils::getValueForKey(mapHeader, "units_offset", vecTmp, IOUtils::nothrow); //TODO: implement these!!
	vecUnitsOffset.clear();
	if(vecTmp.size()==0) {
		for(unsigned int i=0; i<vecDataSequence.size(); i++)
			vecUnitsOffset.push_back(0.);
	} else if(vecTmp.size()==vecDataSequence.size()) {
		for(unsigned int i=0; i<vecDataSequence.size(); i++) {
			double tmp;
			IOUtils::convertString(tmp, vecTmp[i]);
			vecUnitsOffset.push_back(tmp);
		}
	} else {
		throw InvalidFormatException("header lines \"units_offset\" and \"fields\" do not match in "+ filename, AT);
	}

	//trying to read units multipliers
	vecTmp.clear();
	IOUtils::getValueForKey(mapHeader, "units_multiplier", vecTmp, IOUtils::nothrow);
	vecUnitsMultiplier.clear();
	if(vecTmp.size()==0) {
		for(unsigned int i=0; i<vecDataSequence.size(); i++)
			vecUnitsMultiplier.push_back(1.);
	} else if(vecTmp.size()==vecDataSequence.size()) {
		for(unsigned int i=0; i<vecDataSequence.size(); i++) {
			double tmp;
			IOUtils::convertString(tmp, vecTmp[i]);
			vecUnitsMultiplier.push_back(tmp);
		}
	} else {
		throw InvalidFormatException("header lines \"units_multipliers\" and \"fields\" do not match in "+ filename, AT);
	}
626
627
628
629
630
631
632
633
634
635
636

	//Read [DATA] section tag
	getline(fin, line, eoln);
	IOUtils::stripComments(line);
	IOUtils::trim(line);
	IOUtils::toUpper(line);

	if (line != "[DATA]")
		throw InvalidFormatException("Section " + line + " in "+ filename + " invalid, expected [DATA]", AT);
}

637
void SMETIO::checkSignature(const std::vector<std::string>& vecSignature, const std::string& filename, bool& isAscii)
638
{
639
	if ((vecSignature.size() != 3) || (vecSignature[0] != "SMET"))
640
641
		throw InvalidFormatException("The signature of file " + filename + " is invalid", AT);

642
	std::string version = vecSignature[1];
643
	if ((version != "0.9") && (version != "0.95") && (version != "0.99") && (version != smet_version))
644
645
		throw InvalidFormatException("Unsupported file format version for file " + filename, AT);

646
647
	const std::string type = vecSignature[2];
	if (type == "ASCII")
648
		isAscii = true;
649
	else if (type == "BINARY")
650
651
		isAscii = false;
	else 
652
		throw InvalidFormatException("The 3rd column in the signature of file " + filename + " must be either ASCII or BINARY", AT);
653
654
}

655
void SMETIO::writeMeteoData(const std::vector< std::vector<MeteoData> >& vecMeteo, const std::string&)
656
657
658
659
660
{
	//Loop through all stations
	for (unsigned int ii=0; ii<vecMeteo.size(); ii++){
		//1. check consitency of station data position -> write location in header or data section
		StationData sd;
661
		sd.position.setProj(coordout, coordoutparam);
662
		bool isConsistent = checkConsistency(vecMeteo.at(ii), sd);
663
664
665
666
667
668
669
670
		
		if (sd.stationID == ""){
			stringstream ss;
			ss << "Station" << ii+1;

			sd.stationID = ss.str();
		}

671
		string filename = outpath + "/" + sd.stationID + ".smet";
672
673
674
675
676
677
678
		if (!IOUtils::validFileName(filename)) //Check whether filename is valid
			throw InvalidFileNameException(filename, AT);

		try {
			fout.open(filename.c_str());
			if (fout.fail()) throw FileAccessException(filename.c_str(), AT);

679
			fout << "SMET " << smet_version << " ";
680
681
682
683
684
685
686
687
688
			if (outputIsAscii)
				fout << "ASCII" << endl;
			else 
				fout << "BINARY" << endl;

			//2. write header, but first check which meteo parameter fields are actually in use
			vector<bool> vecParamInUse = vector<bool>(MeteoData::nrOfParameters, false); 
			double timezone = IOUtils::nodata;
			checkForUsedParameters(vecMeteo[ii], timezone, vecParamInUse);
689
690
691
692
693
			if(out_dflt_TZ!=IOUtils::nodata) {
				writeHeaderSection(isConsistent, sd, out_dflt_TZ, vecParamInUse);
			} else {
				writeHeaderSection(isConsistent, sd, timezone, vecParamInUse);
			}
694
695
696
			
			//3. write data depending on ASCII/BINARY
			if (outputIsAscii) {
697
				writeDataAscii(isConsistent, vecMeteo[ii], vecParamInUse);
698
699
700
701
702
			} else {
				fout << "[DATA]" << endl;
				fout.close();
				fout.open(filename.c_str(), ios::out | ios::app | ios::binary); //reopen as binary file
				if (fout.fail()) throw FileAccessException(filename.c_str(), AT);
703
				writeDataBinary(isConsistent, vecMeteo[ii], vecParamInUse);
704
705
706
707
708
709
710
711
712
713
714
715
			}
			
			cleanup();

			//4. gzip file or not
		} catch (exception& e){
			cleanup();
			throw;
		}
	}
}

716
void SMETIO::writeDataBinary(const bool& writeLocationInHeader, const std::vector<MeteoData>& vecMeteo,
717
                             const std::vector<bool>& vecParamInUse)
718
719
720
721
{
	char eoln = '\n';

	for (unsigned int ii=0; ii<vecMeteo.size(); ii++){
722
		float val = 0;
723
724
725
726
727
728
729
730
731

		double julian;
		if(out_dflt_TZ!=IOUtils::nodata) {
			Date tmp_date(vecMeteo[ii].date);
			tmp_date.setTimeZone(out_dflt_TZ);
			julian = tmp_date.getJulianDate();
		} else {
			julian = vecMeteo[ii].date.getJulianDate();
		}		
732
		fout.write((char*)&julian, sizeof(double));
733
734

		if (!writeLocationInHeader){ //Meta data changes
735
			val = (float)vecMeteo[ii].meta.position.getLat();
736
			fout.write((char*)&val, sizeof(float));
737
			val = (float)vecMeteo[ii].meta.position.getLon();
738
			fout.write((char*)&val, sizeof(float));
739
			val = (float)vecMeteo[ii].meta.position.getAltitude();
740
741
742
743
744
745
746
747
748
749
750
751
752
			fout.write((char*)&val, sizeof(float));
		}

		for (unsigned int jj=0; jj<MeteoData::nrOfParameters; jj++){
			if (vecParamInUse[jj]){
				val = (float)vecMeteo[ii].param(jj);
				fout.write((char*)&val, sizeof(float));
			}
		}
		fout.write((char*)&eoln, sizeof(char));
	}
}

753
void SMETIO::writeDataAscii(const bool& writeLocationInHeader, const std::vector<MeteoData>& vecMeteo,
754
                            const std::vector<bool>& vecParamInUse)
755
756
{
	fout << "[DATA]" << endl;
757
758
	fout.fill(' ');
	fout << right;
759
	for (unsigned int ii=0; ii<vecMeteo.size(); ii++){
760
761
762
763
764
765
766
		if(out_dflt_TZ!=IOUtils::nodata) {
			Date tmp_date(vecMeteo[ii].date);
			tmp_date.setTimeZone(out_dflt_TZ);
			fout << tmp_date.toString(Date::ISO);
		} else {
			fout << vecMeteo[ii].date.toString(Date::ISO);
		}
767
768

		if (!writeLocationInHeader){ //Meta data changes
769
770
771
			fout << " " << setw(12) << setprecision(6) << vecMeteo[ii].meta.position.getLat();
			fout << " " << setw(12) << setprecision(6) << vecMeteo[ii].meta.position.getLon();
			fout << " " << setw(8)  << setprecision(2) << vecMeteo[ii].meta.position.getAltitude();
772
773
774
		}

		for (unsigned int jj=0; jj<MeteoData::nrOfParameters; jj++){
775
776
777
778
779
780
781
			fout << " ";
			if (vecParamInUse[jj]){
				setFormatting(MeteoData::Parameters(jj));
				if (vecMeteo[ii].param(jj) == IOUtils::nodata)
					fout << setprecision(0);
				fout << vecMeteo[ii].param(jj);
			}
782
783
784
785
786
		}
		fout << endl;
	}
}

787
void SMETIO::setFormatting(const MeteoData::Parameters& paramindex)
788
789
790
{
	if ((paramindex == MeteoData::TA) || (paramindex == MeteoData::TSS) || (paramindex == MeteoData::TSG))
		fout << setw(8) << setprecision(2);
791
	else if ((paramindex == MeteoData::VW) || (paramindex == MeteoData::VW_MAX))
792
793
794
795
796
797
798
799
800
801
802
803
804
		fout << setw(6) << setprecision(1);
	else if (paramindex == MeteoData::DW)
		fout << setw(5) << setprecision(0);
	else if ((paramindex == MeteoData::ISWR) || (paramindex == MeteoData::RSWR) || (paramindex == MeteoData::ILWR))
		fout << setw(6) << setprecision(0);
	else if (paramindex == MeteoData::HNW)
		fout << setw(6) << setprecision(3);
	else if (paramindex == MeteoData::HS)
		fout << setw(8) << setprecision(3);
	else if (paramindex == MeteoData::RH)
		fout << setw(7) << setprecision(3);
}

805
void SMETIO::writeHeaderSection(const bool& writeLocationInHeader, const StationData& sd,
806
                                const double& timezone, const std::vector<bool>& vecParamInUse)
807
808
{
	fout << "[HEADER]" << endl;
809
	fout << "station_id   = " << sd.getStationID() << endl;
810
811
	if (sd.getStationName() != "")
		fout << "station_name = " << sd.getStationName() << endl;
812
813

	if (writeLocationInHeader){ //TODO: only write if != nodata
814
		fout << fixed;
815
816
817
818
819
820
		fout << "latitude     = " << setw(14) << setprecision(6) << sd.position.getLat() << "\n";
		fout << "longitude    = " << setw(14) << setprecision(6) << sd.position.getLon() << "\n";
		fout << "altitude     = " << setw(9)  << setprecision(1) << sd.position.getAltitude() << "\n";
		fout << "easting      = " << setw(14) << setprecision(6) << sd.position.getEasting() << "\n";
		fout << "northing     = " << setw(14) << setprecision(6) << sd.position.getNorthing() << "\n";
		fout << "epsg         = " << setw(7)  << setprecision(0) << sd.position.getEPSG() << "\n";
821
822
	}

823
	fout << "nodata       = " << setw(7) << setprecision(0) << IOUtils::nodata << "\n";
824
825
	
	if ((timezone != IOUtils::nodata) && (timezone != 0.0))
826
		fout << "tz           = " << setw(7)  << setprecision(0) << timezone << "\n";
827

828
	fout << "fields       = timestamp";
829
830
831
832
833
834

	if (!writeLocationInHeader){
		fout << " latitude longitude altitude";
	}

	for (unsigned int ii=0; ii<MeteoData::nrOfParameters; ii++){
835
836
837
838
839
840
		if (vecParamInUse[ii]) {
			std::string column=MeteoData::getParameterName(ii);
			if(column=="RSWR") column="OSWR";
			if(column=="HNW") column="PSUM";
			fout << " " << column;
		}
841
842
843
844
845
	}
	fout << endl;
}


846
void SMETIO::checkForUsedParameters(const std::vector<MeteoData>& vecMeteo, double& timezone,
847
                                    std::vector<bool>& vecParamInUse)
848
849
850
851
852
853
854
855
856
857
858
859
860
{
	for (unsigned int ii=0; ii<vecMeteo.size(); ii++){
		for (unsigned int jj=0; jj<MeteoData::nrOfParameters; jj++){
			if (!vecParamInUse[jj])
				if (vecMeteo[ii].param(jj) != IOUtils::nodata)
					vecParamInUse[jj] = true;
		}
	}	

	if (vecMeteo.size() > 0)
		timezone = vecMeteo[0].date.getTimeZone();
}

861
bool SMETIO::checkConsistency(const std::vector<MeteoData>& vecMeteo, StationData& sd)
862
{
863
864
	if (vecMeteo.size() > 0) //HACK to get the station data even when encoutering bug 87
		sd = vecMeteo[0].meta;
865

866
867
	for (unsigned int ii=1; ii<vecMeteo.size(); ii++){
		if (vecMeteo[ii].meta.position != vecMeteo[ii-1].meta.position)
868
			return false; //BUG: if one is nodata -> we consider that positions are not consistent
869
870
871
872
873
	}

	return true;
}

874
void SMETIO::readSpecialPoints(std::vector<Coords>&)
875
876
877
878
879
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

880
void SMETIO::write2DGrid(const Grid2DObject& /*grid_in*/, const std::string& /*name*/)
881
882
883
884
885
886
887
888
889
890
891
892
893
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}


#ifndef _METEOIO_JNI
extern "C"
{
	void deleteObject(void* obj) {
		delete reinterpret_cast<PluginObject*>(obj);
	}

894
	void* loadObject(const string& classname, const Config& cfg) {
895
		if(classname == "SMETIO") {
896
			//cerr << "Creating dynamic handle for " << classname << endl;
897
			return new SMETIO(deleteObject, cfg);
898
899
900
901
902
903
904
905
		}
		//cerr << "Could not load " << classname << endl;
		return NULL;
	}
}
#endif

} //namespace