WSL/SLF GitLab Repository

SMETIO.cc 27.2 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! (and TZ)
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
239
240
241
242
	//default timezones
	in_dflt_TZ = out_dflt_TZ = 0.;
	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
243
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
244
245

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

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

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

264
265
	nr_stations = counter - 1;

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

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

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

298
299
}

300
void SMETIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
                            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());
323
324
		vecMeteo.reserve(nr_stations);
		vecStation.reserve(nr_stations);
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
355
356
357
	}

	//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);
358
			SMETIO::checkColumnNames(vecDataSequence, locationInHeader);
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381

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

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

389
390
391
	vecMeteo.reserve(buffer_reserve);
	vecStation.reserve(buffer_reserve);

392
393
	while (!fin.eof()){
		MeteoData md;
394
395
		if ((timezone != IOUtils::nodata) && (timezone != 0.0))
			md.date.setTimeZone(timezone);
396
397
398
399
400
401
		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"){
402
403
404
405
				double tmpval;
				fin.read(reinterpret_cast < char * > (&tmpval), sizeof(double));

				md.date.setDate(tmpval);
406
407
408
409

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

		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;
442
443
444
445
		if (md.date >= dateStart){
			vecMeteo.push_back(md);
			vecStation.push_back(tmpsd);
		}
446
447
448
	}	
}

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

458
459
460
461
	tmpvec.reserve(nrOfColumns);
	vecMeteo.reserve(buffer_reserve);
	vecStation.reserve(buffer_reserve);

462
463
464
465
466
467
468
469
470
	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)
471
			throw InvalidFormatException("In "+ filename + ": Invalid amount of data in data line "+line, AT);
472
473

		MeteoData md;
474
475
		if ((timezone != IOUtils::nodata) && (timezone != 0.0))
			md.date.setTimeZone(timezone);
476
477
478
479
480
481
		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]))
482
					throw InvalidFormatException("In "+filename+": Timestamp "+tmpvec[ii]+" invalid in data line", AT);
483
484
485
				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
	}

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

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

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

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

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

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

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

611
void SMETIO::writeMeteoData(const std::vector< std::vector<MeteoData> >& vecMeteo,
612
613
614
615
616
617
618
619
					    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;
620
		sd.position.setProj(coordout, coordoutparam);
621
622
623
624
625
626
627
628
629
		bool isConsistent = checkConsistency(vecStation.at(ii), sd);
		
		if (sd.stationID == ""){
			stringstream ss;
			ss << "Station" << ii+1;

			sd.stationID = ss.str();
		}

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

638
			fout << "SMET " << smet_version << " ";
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
670
			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;
		}
	}
}

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

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

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

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

		if (!writeLocationInHeader){ //Meta data changes
711
712
713
			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();
714
715
716
		}

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

729
void SMETIO::setFormatting(const MeteoData::Parameters& paramindex)
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
{
	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);
}

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

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

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

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

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

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


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

803
bool SMETIO::checkConsistency(const std::vector<StationData>& vecStation, StationData& sd)
804
805
806
807
808
809
810
811
812
813
814
815
{
	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;
}

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

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


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

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

} //namespace