WSL/SLF GitLab Repository

SMETIO.cc 27 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
20
21
22
23

using namespace std;

namespace mio {
/**
24
 * @page smetio SMET
25
 * @section template_format Format
26
 * 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.
27
28
 *
 * @section template_units Units
29
 * 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).
30
31
32
 *
 * @section template_keywords Keywords
 * This plugin uses the following keywords:
33
34
35
36
37
38
39
40
41
42
43
44
45
46
 * - METEOFILE#: input filename and path. As many meteofiles as needed may be specified
 * - METEOPATH: output directory where to write the output meteofiles
 * - METEOPARAM: output file format options (ASCII or BINARY that might be followed by GZIP)
 *
 * Example:
 * @code
 * [Input]
 * METEOFILE1 = ./input/uppper_station.smet
 * METEOFILE2 = ./input/lower_station.smet
 * METEOFILE3 = ./input/outlet_station.smet
 * [Output]
 * METEOPATH = ./output
 * METEOPARAM = ASCII GZIP
 * @endcode
47
48
 */

49
const std::string SMETIO::smet_version = "0.99";
50
const unsigned int SMETIO::buffer_reserve = 23*24*2; //kind of average size of a buffer for optimizing vectors
51
52
map<string, MeteoData::Parameters> SMETIO::mapParameterByName;
const bool SMETIO::__init = SMETIO::initStaticData();
53

54
bool SMETIO::initStaticData()
55
56
57
58
59
60
61
62
{
	for (unsigned int ii=0; ii<MeteoData::nrOfParameters; ii++){
		mapParameterByName[MeteoData::getParameterName(ii)] = MeteoData::Parameters(ii);
	}

	return true;
}

63
double& SMETIO::getParameter(const std::string& columnName, MeteoData& md)
64
65
66
67
68
69
{//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") {
70
		MeteoData::Parameters paramindex = mapParameterByName["HNW"];
71
72
73
		return md.param(paramindex);
	}

74
75
76
77
	MeteoData::Parameters paramindex = mapParameterByName[columnName];
	return md.param(paramindex);
}

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

	for (unsigned int ii=0; ii<vecColumns.size(); ii++){
87
88
89
90
91
92
93
94
		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")){
95
96
			//everything ok
		} else {
97
			map<string, MeteoData::Parameters>::iterator it = mapParameterByName.find(column);
98
99
100
			if (it == mapParameterByName.end())
				throw InvalidFormatException("Key 'fields' specified in [HEADER] section contains invalid names", AT);
			
101
			paramcounter.at(mapParameterByName[column])++;
102
103
		}
	}
104
105
	
	//Check for multiple usages of parameters
106
107
108
109
110
	for (unsigned int ii=0; ii<paramcounter.size(); ii++){
		if (paramcounter[ii] > 1)
			throw InvalidFormatException("In 'fields': Multiple use of " + MeteoData::getParameterName(ii), AT);
	}
	
111
112
	//If there is no location information in the [HEADER] section, then
	//location information must be part of fields
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
	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


132
SMETIO::SMETIO(void (*delObj)(void*), const Config& i_cfg) : IOInterface(delObj), cfg(i_cfg)
133
134
135
136
{
	parseInputOutputSection();
}

137
SMETIO::SMETIO(const std::string& configfile) : IOInterface(NULL), cfg(configfile)
138
139
140
141
{
	parseInputOutputSection();
}

142
SMETIO::SMETIO(const Config& cfgreader) : IOInterface(NULL), cfg(cfgreader)
143
144
145
146
{
	parseInputOutputSection();
}

147
SMETIO::~SMETIO() throw()
148
149
150
151
{
	cleanup();
}

152
void SMETIO::cleanup() throw()
153
{
154
155
156
	//clear ios flags
	fout << resetiosflags(ios_base::fixed | ios_base::left);

157
158
159
160
161
162
163
164
	if (fin.is_open()) {//close fin if open
		fin.close();
	}
	if (fout.is_open()) {//close fout if open
		fout.close();
	}		
}

165
void SMETIO::read2DGrid(Grid2DObject& /*grid_out*/, const std::string& /*_name*/)
166
167
168
169
170
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

171
void SMETIO::readDEM(DEMObject& /*dem_out*/)
172
173
174
175
176
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

177
void SMETIO::readLanduse(Grid2DObject& /*landuse_out*/)
178
179
180
181
182
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

183
void SMETIO::readAssimilationData(const Date& /*date_in*/, Grid2DObject& /*da_out*/)
184
185
186
187
188
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

189
190
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
191
//ie: it should use the given date!
192
	vecStation.clear();
193
	vecStation.reserve(nr_stations);
194
195

	//Now loop through all requested stations, open the respective files and parse them
196
	for (unsigned int ii=0; ii<nr_stations; ii++){
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
		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;
		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
			readHeader(eoln, filename, locationInHeader, timezone, sd, vecDataSequence);
			cleanup();
		} catch(std::exception& e) {
			cleanup();
			throw;
		}
		vecStation.push_back(sd);
	}
233
234
}

235
void SMETIO::parseInputOutputSection()
236
{
237
	/*
238
	 * Parse the [Input] and [Output] sections within Config object cfg
239
	 */
240
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
241
242

	//Parse input section: extract number of files to read and store filenames in vecFiles
243
244
245
246
247
248
249
	unsigned int counter = 1;
	string filename = "";
	do {
		stringstream ss;
		filename = "";
		
		ss << "METEOFILE" << counter;
250
		cfg.getValue(ss.str(), "Input", filename, Config::nothrow);
251
252
253
254
255
256
257
258
259
260

		if (filename != ""){
			if (!IOUtils::validFileName(filename)) //Check whether filename is valid
				throw InvalidFileNameException(filename, AT);

			vecFiles.push_back(filename);
		}
		counter++;
	} while (filename != "");

261
262
	nr_stations = counter - 1;

263
	//Parse output section: extract info on whether to write ASCII or BINARY format, gzipped or not
264
265
266
267
268
	outpath = "";
	outputIsAscii = true; 
	outputIsGzipped = false;

	vector<string> vecArgs;
269
270
	cfg.getValue("METEOPATH", "Output", outpath, Config::nothrow);
	cfg.getValue("METEOPARAM", "Output", vecArgs, Config::nothrow); //"ASCII|BINARY GZIP"
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293

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

295
296
}

297
void SMETIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
                            std::vector< std::vector<MeteoData> >& vecMeteo, 
                            std::vector< std::vector<StationData> >& vecStation,
                            const unsigned int& stationindex)
{
	//Make sure that vecMeteo/vecStation have the correct dimension and stationindex is valid
	unsigned int startindex=0, endindex=vecFiles.size();	
	if (stationindex != IOUtils::npos){
		if ((stationindex < vecFiles.size()) || (stationindex < vecMeteo.size()) || (stationindex < vecStation.size())){
			startindex = stationindex;
			endindex = stationindex+1;
		} else {
			throw IndexOutOfBoundsException("Invalid stationindex", AT);
		} 

		vecMeteo[stationindex].clear();
		vecStation[stationindex].clear();
	} else {
		vecMeteo.clear();
		vecStation.clear();
		
		vecMeteo = vector< vector<MeteoData> >(vecFiles.size());
		vecStation = vector< vector<StationData> >(vecFiles.size());
320
321
		vecMeteo.reserve(nr_stations);
		vecStation.reserve(nr_stations);
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
	}

	//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;
		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
			readHeader(eoln, filename, locationInHeader, timezone, sd, vecDataSequence);
355
			SMETIO::checkColumnNames(vecDataSequence, locationInHeader);
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378

			//3. Read DATA
			if (isAscii){
				readDataAscii(eoln, filename, timezone, sd, vecDataSequence, dateStart, dateEnd,
						    vecMeteo[ii], vecStation[ii]);
			} 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

				readDataBinary(eoln, filename, timezone, sd, vecDataSequence, dateStart, dateEnd,
							vecMeteo[ii], vecStation[ii]);
			}
			cleanup();
		} catch(std::exception& e) {
			cleanup();
			throw;
		}
	}
}

379
void SMETIO::readDataBinary(const char&, const std::string&, const double& timezone,
380
381
382
                            const StationData& sd, const std::vector<std::string>& vecDataSequence,
                            const Date& dateStart, const Date& dateEnd,
                            std::vector<MeteoData>& vecMeteo, std::vector<StationData>& vecStation)
383
384
385
{
	const unsigned int nrOfColumns = vecDataSequence.size();

386
387
388
	vecMeteo.reserve(buffer_reserve);
	vecStation.reserve(buffer_reserve);

389
390
391
392
393
394
395
396
	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"){
397
398
399
400
				double tmpval;
				fin.read(reinterpret_cast < char * > (&tmpval), sizeof(double));

				md.date.setDate(tmpval);
401
402
403
404
405
406
407

				if ((timezone != IOUtils::nodata) && (timezone != 0.0))
					md.date.setTimeZone(timezone);

				if (md.date < dateStart)
					continue;
				if (md.date > dateEnd)
408
					return;
409
			} else {
410
411
412
413
414
415
416
417
418
419
420
421
422
423
				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 {
424
					SMETIO::getParameter(vecDataSequence[ii], md) = val;
425
				}
426
427
428
429
430
431
432
433
434
435
436
437
438
439
			}
		}

		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);

		//cout << "===" << endl;
		//cout << sd << endl;
		//cout << md.date.toString(Date::ISO) << endl;
440
441
442
443
		if (md.date >= dateStart){
			vecMeteo.push_back(md);
			vecStation.push_back(tmpsd);
		}
444
445
446
	}	
}

447
void SMETIO::readDataAscii(const char& eoln, const std::string& filename, const double& timezone,
448
449
450
                           const StationData& sd, const std::vector<std::string>& vecDataSequence,
                           const Date& dateStart, const Date& dateEnd,
                           std::vector<MeteoData>& vecMeteo, std::vector<StationData>& vecStation)
451
452
453
454
455
{
	string line = "";
	vector<string> tmpvec;
	const unsigned int nrOfColumns = vecDataSequence.size();

456
457
458
459
	tmpvec.reserve(nrOfColumns);
	vecMeteo.reserve(buffer_reserve);
	vecStation.reserve(buffer_reserve);

460
461
462
463
464
465
466
467
468
	while (!fin.eof()){
		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)
469
			throw InvalidFormatException("In "+ filename + ": Invalid amount of data in data line "+line, AT);
470
471
472
473
474
475
476
477

		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"){
				if (!IOUtils::convertString(md.date, tmpvec[ii]))
478
					throw InvalidFormatException("In "+filename+": Timestamp "+tmpvec[ii]+" invalid in data line", AT);
479
480
481
482
483
484
485

				if ((timezone != IOUtils::nodata) && (timezone != 0.0))
					md.date.setTimeZone(timezone);

				if (md.date < dateStart)
					continue;
				if (md.date > dateEnd)
486
					return;
487
488
489
490
491
492
493
494
495
496
497
498
499
500

			} 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 {
501
				if (!IOUtils::convertString(SMETIO::getParameter(vecDataSequence[ii], md), tmpvec[ii]))
502
503
504
505
506
507
508
					throw InvalidFormatException("In "+filename+": Invalid value for param", AT);
			}
		}

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

509
510
511
512
		if (md.date >= dateStart){
			vecMeteo.push_back(md);
			vecStation.push_back(tmpsd);
		}
513
514
515
	}
}

516
void SMETIO::readHeader(const char& eoln, const std::string& filename, bool& locationInHeader,
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
                         double& timezone, StationData& sd, std::vector<std::string>& vecDataSequence)
{
	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);

538
		if (line != "") {
539
540
			if (!IOUtils::readKeyValuePair(line, "=", mapHeader))
				throw InvalidFormatException("Invalid key value pair in section [Header]", AT);
541
		}
542
543
544
545
546
547
	}

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

550
551
552
	//trying to read easting/northing
	double easting=IOUtils::nodata, northing=IOUtils::nodata, alt=IOUtils::nodata;
	short int epsg=IOUtils::snodata;
553
554
555
556
	IOUtils::getValueForKey(mapHeader, "epsg", epsg, IOUtils::nothrow);
	if(epsg!=IOUtils::snodata) {
		sd.position.setEPSG(epsg);
	}
557
	
558
559
560
561
562
563
564
565
566
567
568
569
	IOUtils::getValueForKey(mapHeader, "easting", easting, IOUtils::nothrow);
	if (easting != IOUtils::nodata){ //HACK
		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;
570
571
572
573
574
575
576
577
578
	IOUtils::getValueForKey(mapHeader, "latitude", lat, IOUtils::nothrow);
	if (lat != IOUtils::nodata){ //HACK
		IOUtils::getValueForKey(mapHeader, "longitude", lon);
		IOUtils::getValueForKey(mapHeader, "altitude", alt);
		sd.position.setLatLon(lat, lon, alt);
		locationInHeader = true;
	}

	IOUtils::getValueForKey(mapHeader, "fields", vecDataSequence);
579
	//IOUtils::getValueForKey(mapHeader, "units_offset", vecUnitsOffset, IOUtils::nothrow); //TODO: implement these!!
580
	//IOUtils::getValueForKey(mapHeader, "units_multiplier", vecUnitsMultiplier, IOUtils::nothrow);
581
582
583
584
585
586
587
588
589
590
591

	//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);
}

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

597
	std::string version = vecSignature[1];
598
	if ((version != "0.9") && (version != "0.95") && (version != smet_version))
599
600
		throw InvalidFormatException("Unsupported file format version for file " + filename, AT);

601
602
	const std::string type = vecSignature[2];
	if (type == "ASCII")
603
		isAscii = true;
604
	else if (type == "BINARY")
605
606
		isAscii = false;
	else 
607
		throw InvalidFormatException("The 3rd column in the signature of file " + filename + " must be either ASCII or BINARY", AT);
608
609
}

610
void SMETIO::writeMeteoData(const std::vector< std::vector<MeteoData> >& vecMeteo,
611
612
613
614
615
616
617
618
					    const std::vector< std::vector<StationData> >& vecStation,
					    const std::string&)
{

	//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;
619
		sd.position.setProj(coordout, coordoutparam);
620
621
622
623
624
625
626
627
628
		bool isConsistent = checkConsistency(vecStation.at(ii), sd);
		
		if (sd.stationID == ""){
			stringstream ss;
			ss << "Station" << ii+1;

			sd.stationID = ss.str();
		}

629
		string filename = outpath + "/" + sd.stationID + ".smet";
630
631
632
633
634
635
636
		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);

637
			fout << "SMET " << smet_version << " ";
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
			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);
			writeHeaderSection(isConsistent, sd, timezone, vecParamInUse);
			
			//3. write data depending on ASCII/BINARY
			if (outputIsAscii) {
				writeDataAscii(isConsistent, vecMeteo[ii], vecStation[ii], vecParamInUse);
			} 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);
				writeDataBinary(isConsistent, vecMeteo[ii], vecStation[ii], vecParamInUse);
			}
			
			cleanup();

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

670
void SMETIO::writeDataBinary(const bool& writeLocationInHeader, const std::vector<MeteoData>& vecMeteo,
671
672
673
674
675
                              const std::vector<StationData>& vecStation, const std::vector<bool>& vecParamInUse)
{
	char eoln = '\n';

	for (unsigned int ii=0; ii<vecMeteo.size(); ii++){
676
677
678
679
		float val = 0;
		double julian = vecMeteo[ii].date.getJulianDate();
		
		fout.write((char*)&julian, sizeof(double));
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699

		if (!writeLocationInHeader){ //Meta data changes
			val = (float)vecStation[ii].position.getLat();
			fout.write((char*)&val, sizeof(float));
			val = (float)vecStation[ii].position.getLon();
			fout.write((char*)&val, sizeof(float));
			val = (float)vecStation[ii].position.getAltitude();
			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));
	}
}

700
void SMETIO::writeDataAscii(const bool& writeLocationInHeader, const std::vector<MeteoData>& vecMeteo,
701
702
703
                             const std::vector<StationData>& vecStation, const std::vector<bool>& vecParamInUse)
{
	fout << "[DATA]" << endl;
704
705
	fout.fill(' ');
	fout << right;
706
707
708
709
	for (unsigned int ii=0; ii<vecMeteo.size(); ii++){
		fout << vecMeteo[ii].date.toString(Date::ISO);

		if (!writeLocationInHeader){ //Meta data changes
710
711
712
			fout << " " << setw(12) << setprecision(6) << vecStation[ii].position.getLat();
			fout << " " << setw(12) << setprecision(6) << vecStation[ii].position.getLon();
			fout << " " << setw(8)  << setprecision(2) << vecStation[ii].position.getAltitude();
713
714
715
		}

		for (unsigned int jj=0; jj<MeteoData::nrOfParameters; jj++){
716
717
718
719
720
721
722
			fout << " ";
			if (vecParamInUse[jj]){
				setFormatting(MeteoData::Parameters(jj));
				if (vecMeteo[ii].param(jj) == IOUtils::nodata)
					fout << setprecision(0);
				fout << vecMeteo[ii].param(jj);
			}
723
724
725
726
727
		}
		fout << endl;
	}
}

728
void SMETIO::setFormatting(const MeteoData::Parameters& paramindex)
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
{
	if ((paramindex == MeteoData::TA) || (paramindex == MeteoData::TSS) || (paramindex == MeteoData::TSG))
		fout << setw(8) << setprecision(2);
	else if (paramindex == MeteoData::VW)
		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);
}

746
void SMETIO::writeHeaderSection(const bool& writeLocationInHeader, const StationData& sd,
747
748
749
                                 const double& timezone, const std::vector<bool>& vecParamInUse)
{
	fout << "[HEADER]" << endl;
750
	fout << "station_id   = " << sd.getStationID() << endl;
751
752
	if (sd.getStationName() != "")
		fout << "station_name = " << sd.getStationName() << endl;
753
754

	if (writeLocationInHeader){ //TODO: only write if != nodata
755
		fout << fixed;
756
757
758
759
760
761
		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";
762
763
	}

764
	fout << "nodata       = " << setw(7) << setprecision(0) << IOUtils::nodata << "\n";
765
766
	
	if ((timezone != IOUtils::nodata) && (timezone != 0.0))
767
		fout << "tz           = " << setw(7)  << setprecision(0) << timezone << "\n";
768

769
	fout << "fields       = timestamp";
770
771
772
773
774
775

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

	for (unsigned int ii=0; ii<MeteoData::nrOfParameters; ii++){
776
777
778
779
780
781
		if (vecParamInUse[ii]) {
			std::string column=MeteoData::getParameterName(ii);
			if(column=="RSWR") column="OSWR";
			if(column=="HNW") column="PSUM";
			fout << " " << column;
		}
782
783
784
785
786
	}
	fout << endl;
}


787
void SMETIO::checkForUsedParameters(const std::vector<MeteoData>& vecMeteo, double& timezone,
788
789
790
791
792
793
794
795
796
797
798
799
800
801
                                     std::vector<bool>& vecParamInUse)
{
	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();
}

802
bool SMETIO::checkConsistency(const std::vector<StationData>& vecStation, StationData& sd)
803
804
805
806
807
808
809
810
811
812
813
814
{
	for (unsigned int ii=1; ii<vecStation.size(); ii++){
		if (vecStation[ii].position != vecStation[ii-1].position)
			return false;
	}

	if (vecStation.size() > 0)
		sd = vecStation[0];

	return true;
}

815
void SMETIO::readSpecialPoints(std::vector<Coords>&)
816
817
818
819
820
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

821
void SMETIO::write2DGrid(const Grid2DObject& /*grid_in*/, const std::string& /*name*/)
822
823
824
825
826
827
828
829
830
831
832
833
834
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}


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

835
	void* loadObject(const string& classname, const Config& cfg) {
836
		if(classname == "SMETIO") {
837
			//cerr << "Creating dynamic handle for " << classname << endl;
838
			return new SMETIO(deleteObject, cfg);
839
840
841
842
843
844
845
846
		}
		//cerr << "Could not load " << classname << endl;
		return NULL;
	}
}
#endif

} //namespace