WSL/SLF GitLab Repository

SMETIO.cc 31.7 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.1";
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
	for (size_t ii=0; ii<MeteoData::nrOfParameters; ii++){
61
62
63
64
65
66
		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
	 * This function checks whether the sequence of keywords specified in the
86
87
	 * [HEADER] section (key 'fields') is valid
	 */
88
	vector<size_t> paramcounter = vector<size_t>(MeteoData::nrOfParameters, 0);
89

90
	for (size_t 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
			if (it == mapParameterByName.end())
				throw InvalidFormatException("Key 'fields' specified in [HEADER] section contains invalid names", AT);
104

105
			paramcounter.at(mapParameterByName[column])++;
106
107
		}
	}
108

109
	//Check for multiple usages of parameters
110
	for (size_t ii=0; ii<paramcounter.size(); ii++){
111
112
113
		if (paramcounter[ii] > 1)
			throw InvalidFormatException("In 'fields': Multiple use of " + MeteoData::getParameterName(ii), AT);
	}
114

115
116
	//If there is no location information in the [HEADER] section, then
	//location information must be part of fields
117
	if (!locationInHeader){
118
119
		size_t latcounter = 0, loncounter=0, altcounter=0;
		for (size_t ii=0; ii<vecColumns.size(); ii++){
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
			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
{
	parseInputOutputSection();
139
	plugin_nodata = IOUtils::nodata;
140
141
}

142
SMETIO::SMETIO(const std::string& configfile) : IOInterface(NULL), cfg(configfile)
143
144
{
	parseInputOutputSection();
145
	plugin_nodata = IOUtils::nodata;
146
147
}

148
SMETIO::SMETIO(const Config& cfgreader) : IOInterface(NULL), cfg(cfgreader)
149
150
{
	parseInputOutputSection();
151
	plugin_nodata = IOUtils::nodata;
152
153
}

154
SMETIO::~SMETIO() throw()
155
156
157
158
{
	cleanup();
}

159
void SMETIO::cleanup() throw()
160
{
161
162
163
	//clear ios flags
	fout << resetiosflags(ios_base::fixed | ios_base::left);

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
	}
170
171
}

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

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

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

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

196
197
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
198
//ie: it should use the given date! (and TZ)
199
	vecStation.clear();
200
	vecStation.reserve(nr_stations);
201
202

	//Now loop through all requested stations, open the respective files and parse them
203
	for (size_t ii=0; ii<nr_stations; ii++){
204
205
		bool isAscii = true;
		string filename = vecFiles.at(ii); //filename of current station
206

207
208
209
210
211
212
213
214
215
216
217
218
219
		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;
220
		std::vector<double> vecUnitsOffset, vecUnitsMultiplier;
221
222
223
224
225
226
227
228
229
230
231
232
		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
233
			readHeader(eoln, filename, locationInHeader, timezone, sd, vecDataSequence, vecUnitsOffset, vecUnitsMultiplier);
234
235
236
237
238
239
240
			cleanup();
		} catch(std::exception& e) {
			cleanup();
			throw;
		}
		vecStation.push_back(sd);
	}
241
242
}

243
void SMETIO::parseInputOutputSection()
244
{
245
	//default timezones
246
	in_dflt_TZ = out_dflt_TZ = IOUtils::nodata;
247
248
	cfg.getValue("TIME_ZONE","Input",in_dflt_TZ,Config::nothrow);
	cfg.getValue("TIME_ZONE","Output",out_dflt_TZ,Config::nothrow);
249
250

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

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

261
262
263
		do {
			stringstream ss;
			filename = "";
264

265
266
			ss << "STATION" << counter;
			cfg.getValue(ss.str(), "Input", filename, Config::nothrow);
267

268
269
270
271
272
			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);
273

274
275
276
277
				vecFiles.push_back(file_and_path.str());
			}
			counter++;
		} while (filename != "");
278

279
280
		nr_stations = counter - 1;
	}
281

282
283
284
	if (vec_streampos.size() == 0) //the vec_streampos save file pointers for certain dates
		vec_streampos = vector< map<Date, std::streampos> >(nr_stations);

285
	//Parse output section: extract info on whether to write ASCII or BINARY format, gzipped or not
286
	outpath = "";
287
	outputIsAscii = true;
288
289
290
	outputIsGzipped = false;

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

	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;
307
	else
308
309
310
311
312
313
314
315
		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;
	}
316

317
318
}

319
void SMETIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
320
321
                           std::vector< std::vector<MeteoData> >& vecMeteo,
                           const unsigned int& stationindex)
322
{
323
	//Make sure that vecMeteo have the correct dimension and stationindex is valid
324
	size_t startindex=0, endindex=vecFiles.size();
325
	if (stationindex != (unsigned int)IOUtils::npos){ //HACK do we really still need stationindex??
326
		if ((stationindex < vecFiles.size()) || (stationindex < vecMeteo.size())){
327
328
329
330
			startindex = stationindex;
			endindex = stationindex+1;
		} else {
			throw IndexOutOfBoundsException("Invalid stationindex", AT);
331
		}
332
333
334

		vecMeteo[stationindex].clear();
	} else {
335
		vecMeteo.clear();
336
		vecMeteo = vector< vector<MeteoData> >(vecFiles.size());
337
		vecMeteo.reserve(nr_stations);
338
339
340
	}

	//Now loop through all requested stations, open the respective files and parse them
341
	for (size_t ii=startindex; ii<endindex; ii++){
342
343
		bool isAscii = true;
		string filename = vecFiles.at(ii); //filename of current station
344

345
346
347
348
349
350
351
352
		if (!IOUtils::fileExists(filename))
			throw FileNotFoundException(filename, AT);

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

353
		const char eoln = IOUtils::getEoln(fin); //get the end of line character for the file
354
355
356
357

		//Go through file, save key value pairs
		string line="";
		std::vector<std::string> tmpvec, vecDataSequence;
358
		std::vector<double> vecUnitsOffset, vecUnitsMultiplier;
359
360
361
362
363
364
365
366
367
368
369
370
		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
371
			readHeader(eoln, filename, locationInHeader, timezone, sd, vecDataSequence, vecUnitsOffset, vecUnitsMultiplier);
372
			SMETIO::checkColumnNames(vecDataSequence, locationInHeader);
373

374
375
376
377
378
379
380
381
			//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);
			}

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

393
				readDataBinary(eoln, filename, timezone, sd, vecDataSequence, vecUnitsOffset,
394
				               vecUnitsMultiplier, dateStart, dateEnd, vecMeteo[ii]);
395
396
397
398
399
400
401
402
403
			}
			cleanup();
		} catch(std::exception& e) {
			cleanup();
			throw;
		}
	}
}

404
void SMETIO::readDataBinary(const char&, const std::string&, const double& timezone,
405
                            const StationData& sd, const std::vector<std::string>& vecDataSequence,
406
                            const std::vector<double>& vecUnitsOffset, std::vector<double>& vecUnitsMultiplier,
407
                            const Date& dateStart, const Date& dateEnd, std::vector<MeteoData>& vecMeteo)
408
{
409
	const size_t nrOfColumns = vecDataSequence.size();
410

411
412
	vecMeteo.reserve(buffer_reserve);

413
414
415
416
	while (!fin.eof()){
		MeteoData md;
		StationData tmpsd = sd;
		double lat=IOUtils::nodata, lon=IOUtils::nodata, alt=IOUtils::nodata;
417
		size_t poscounter = 0;
418

419
		for (size_t ii=0; ii<nrOfColumns; ii++){
420
			if (vecDataSequence[ii] == "timestamp"){
421
422
423
				double tmpval;
				fin.read(reinterpret_cast < char * > (&tmpval), sizeof(double));

424
				md.date.setDate(tmpval, timezone);
425
426
427
428

				if (md.date < dateStart)
					continue;
				if (md.date > dateEnd)
429
					return;
430
			} else {
431
				float tmpval;
432
				fin.read(reinterpret_cast < char * > (&tmpval), sizeof(float));
433
434
435
436
437
438
439
440
441
442
443
444
				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 {
445
446
447
					if(val==plugin_nodata)
						SMETIO::getParameter(vecDataSequence[ii], md) = IOUtils::nodata;
					else
448
						SMETIO::getParameter(vecDataSequence[ii], md) = (val * vecUnitsMultiplier[ii]) + vecUnitsOffset[ii];
449
				}
450
451
452
453
454
455
456
457
			}
		}

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

458
459
460
461
		if (poscounter == 3) {
			lat=IOUtils::standardizeNodata(lat, plugin_nodata);
			lon=IOUtils::standardizeNodata(lon, plugin_nodata);
			alt=IOUtils::standardizeNodata(alt, plugin_nodata);
462
			tmpsd.position.setLatLon(lat, lon, alt);
463
		}
464

465
		if (md.date >= dateStart){
466
			md.meta = tmpsd;
467
468
			vecMeteo.push_back(md);
		}
469
	}
470
471
}

472
void SMETIO::readDataAscii(const char& eoln, const std::string& filename, const double& timezone,
473
                           const StationData& sd, const std::vector<std::string>& vecDataSequence,
474
                           const std::vector<double>& vecUnitsOffset, std::vector<double>& vecUnitsMultiplier,
475
                           const Date& dateStart, const Date& dateEnd, std::vector<MeteoData>& vecMeteo, const size_t& stat_idx)
476
477
478
{
	string line = "";
	vector<string> tmpvec;
479
	const size_t nrOfColumns = vecDataSequence.size();
480

481
482
483
	tmpvec.reserve(nrOfColumns);
	vecMeteo.reserve(buffer_reserve);

484
485
486
487
488
489
	//The following 4 lines are an optimization to jump to the correct position in the file
	streampos current_fpointer = -1;  //the filepointer for the current valid date
	map<Date,streampos>::const_iterator it = vec_streampos.at(stat_idx).find(dateStart);
	if (it != vec_streampos.at(stat_idx).end())
		fin.seekg(it->second); //jump to position in the file

490
	while (!fin.eof()){
491
		streampos tmp_fpointer = fin.tellg();
492
493
494
495
496
		getline(fin, line, eoln);
		IOUtils::stripComments(line);
		IOUtils::trim(line);
		if (line == "") continue; //Pure comment lines and empty lines are ignored

497
		const size_t ncols = IOUtils::readLineToVec(line, tmpvec);
498
499

		if (ncols != nrOfColumns)
500
			throw InvalidFormatException("In "+ filename + ": Invalid amount of data in data line "+line, AT);
501
502
503

		MeteoData md;
		StationData tmpsd = sd;
504
		double lat=IOUtils::nodata, lon=IOUtils::nodata, alt=IOUtils::nodata;
505
506
		size_t poscounter = 0;
		for (size_t ii=0; ii<nrOfColumns; ii++){
507
			if (vecDataSequence[ii] == "timestamp"){
508
				if (!IOUtils::convertString(md.date, tmpvec[ii], timezone, std::dec))
509
					throw InvalidFormatException("In "+filename+": Timestamp "+tmpvec[ii]+" invalid in data line", AT);
510
511
				if (md.date < dateStart)
					continue;
512
513
514
				if (md.date > dateEnd) {
					//save stream position and the corresponding end date
					if (current_fpointer != ((ifstream::pos_type)-1)) vec_streampos.at(stat_idx)[dateEnd] = current_fpointer;
515
					return;
516
				}
517
518

			} else if (vecDataSequence[ii] == "latitude"){
519
				if (!IOUtils::convertString(lat, tmpvec[ii]))
520
521
522
					throw InvalidFormatException("In "+filename+": Latitude invalid", AT);
				poscounter++;
			} else if (vecDataSequence[ii] == "longitude"){
523
				if (!IOUtils::convertString(lon, tmpvec[ii]))
524
525
526
					throw InvalidFormatException("In "+filename+": Longitude invalid", AT);
				poscounter++;
			} else if (vecDataSequence[ii] == "altitude"){
527
				if (!IOUtils::convertString(alt, tmpvec[ii]))
528
529
530
					throw InvalidFormatException("In "+filename+": Altitude invalid", AT);
				poscounter++;
			} else {
531
532
				double val;
				if (!IOUtils::convertString(val, tmpvec[ii]))
533
					throw InvalidFormatException("In "+filename+": Invalid value for param", AT);
534
535
536
				if(val==plugin_nodata)
					SMETIO::getParameter(vecDataSequence[ii], md) = IOUtils::nodata;
				else
537
					SMETIO::getParameter(vecDataSequence[ii], md) = (val * vecUnitsMultiplier[ii]) + vecUnitsOffset[ii];
538
539
540
			}
		}

541
542
543
544
		if (poscounter == 3) {
			lat=IOUtils::standardizeNodata(lat, plugin_nodata);
			lon=IOUtils::standardizeNodata(lon, plugin_nodata);
			alt=IOUtils::standardizeNodata(alt, plugin_nodata);
545
			tmpsd.position.setLatLon(lat, lon, alt);
546
		}
547

548
		if (md.date >= dateStart){
549
			md.meta = tmpsd;
550
			vecMeteo.push_back(md);
551
			current_fpointer = tmp_fpointer; //save this file pointer, it's a valid one for sure
552
		}
553
554
555
	}
}

556
void SMETIO::readHeader(const char& eoln, const std::string& filename, bool& locationInHeader,
557
558
                        double& timezone, StationData& sd, std::vector<std::string>& vecDataSequence,
                        std::vector<double>& vecUnitsOffset, std::vector<double>& vecUnitsMultiplier)
559
560
561
562
{
	string line="";
	while (!fin.eof() && (fin.peek() != '['))
		getline(fin, line, eoln);
563

564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
	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);

579
		if (line != "") {
580
581
			if (!IOUtils::readKeyValuePair(line, "=", mapHeader))
				throw InvalidFormatException("Invalid key value pair in section [Header]", AT);
582
		}
583
584
585
	}

	//Now extract info from mapHeader
586
	IOUtils::getValueForKey(mapHeader, "nodata", plugin_nodata);
587
588
	IOUtils::getValueForKey(mapHeader, "station_id", sd.stationID);
	IOUtils::getValueForKey(mapHeader, "station_name", sd.stationName, IOUtils::nothrow);
589
	timezone = in_dflt_TZ;
590
	IOUtils::getValueForKey(mapHeader, "tz", timezone, IOUtils::nothrow);
591
592
	if(timezone==plugin_nodata) //if a nodata was given in the header, we replace it with the io.ini timezone
		timezone = in_dflt_TZ;
593
	sd.position.setProj(coordin, coordinparam); //set the default projection from config file
594

595
596
597
	//trying to read easting/northing
	double easting=IOUtils::nodata, northing=IOUtils::nodata, alt=IOUtils::nodata;
	short int epsg=IOUtils::snodata;
598
599
600
601
	IOUtils::getValueForKey(mapHeader, "epsg", epsg, IOUtils::nothrow);
	if(epsg!=IOUtils::snodata) {
		sd.position.setEPSG(epsg);
	}
602

603
	IOUtils::getValueForKey(mapHeader, "easting", easting, IOUtils::nothrow);
604
	easting=IOUtils::standardizeNodata(easting, plugin_nodata);
605
	if (easting != IOUtils::nodata){
606
607
		IOUtils::getValueForKey(mapHeader, "northing", northing);
		IOUtils::getValueForKey(mapHeader, "altitude", alt);
608
609
		northing=IOUtils::standardizeNodata(northing, plugin_nodata);
		alt=IOUtils::standardizeNodata(alt, plugin_nodata);
610
611
612
613
614
615
616
617
		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;
618
	IOUtils::getValueForKey(mapHeader, "latitude", lat, IOUtils::nothrow);
619
	lat=IOUtils::standardizeNodata(lat, plugin_nodata);
620
	if (lat != IOUtils::nodata){
621
622
		IOUtils::getValueForKey(mapHeader, "longitude", lon);
		IOUtils::getValueForKey(mapHeader, "altitude", alt);
623
624
		lon=IOUtils::standardizeNodata(lon, plugin_nodata);
		alt=IOUtils::standardizeNodata(alt, plugin_nodata);
625
626
627
628
629
		sd.position.setLatLon(lat, lon, alt);
		locationInHeader = true;
	}

	IOUtils::getValueForKey(mapHeader, "fields", vecDataSequence);
630
631
632
633
634
635

	//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) {
636
		for(size_t i=0; i<vecDataSequence.size(); i++)
637
638
			vecUnitsOffset.push_back(0.);
	} else if(vecTmp.size()==vecDataSequence.size()) {
639
		for(size_t i=0; i<vecDataSequence.size(); i++) {
640
641
642
643
644
645
646
647
648
649
650
651
652
			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) {
653
		for(size_t i=0; i<vecDataSequence.size(); i++)
654
655
			vecUnitsMultiplier.push_back(1.);
	} else if(vecTmp.size()==vecDataSequence.size()) {
656
		for(size_t i=0; i<vecDataSequence.size(); i++) {
657
658
659
660
661
662
663
			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);
	}
664
665
666
667
668
669
670
671
672
673
674

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

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

680
	std::string version = vecSignature[1];
681
	if ((version != "0.9") && (version != "0.95") && (version != "0.99") && (version != "1.0") && (version != smet_version))
682
683
		throw InvalidFormatException("Unsupported file format version for file " + filename, AT);

684
685
686
687
	if(version=="0.9" || version=="0.95" || version=="0.99" || version=="1.0") {
		std::cout << "[W] SMET specification 1.1 changes the priorities of units_multiplier and units_offset. Please check/update your files and bring them to 1.1!!\n";
	}

688
689
	const std::string type = vecSignature[2];
	if (type == "ASCII")
690
		isAscii = true;
691
	else if (type == "BINARY")
692
		isAscii = false;
693
	else
694
		throw InvalidFormatException("The 3rd column in the signature of file " + filename + " must be either ASCII or BINARY", AT);
695
696
}

697
void SMETIO::writeMeteoData(const std::vector< std::vector<MeteoData> >& vecMeteo, const std::string&)
698
699
{
	//Loop through all stations
700
	for (size_t ii=0; ii<vecMeteo.size(); ii++){
701
702
		//1. check consitency of station data position -> write location in header or data section
		StationData sd;
703
		sd.position.setProj(coordout, coordoutparam);
704
		const bool isConsistent = checkConsistency(vecMeteo.at(ii), sd);
705

706
707
708
709
710
711
		if (sd.stationID == ""){
			stringstream ss;
			ss << "Station" << ii+1;
			sd.stationID = ss.str();
		}

712
		const string filename = outpath + "/" + sd.stationID + ".smet";
713
714
715
716
717
718
719
		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);

720
			fout << "SMET " << smet_version << " ";
721
722
			if (outputIsAscii)
				fout << "ASCII" << endl;
723
			else
724
725
726
				fout << "BINARY" << endl;

			//2. write header, but first check which meteo parameter fields are actually in use
727
			vector<bool> vecParamInUse = vector<bool>(MeteoData::nrOfParameters, false);
728
729
			double timezone = IOUtils::nodata;
			checkForUsedParameters(vecMeteo[ii], timezone, vecParamInUse);
730
731
732
733
734
			if(out_dflt_TZ!=IOUtils::nodata) {
				writeHeaderSection(isConsistent, sd, out_dflt_TZ, vecParamInUse);
			} else {
				writeHeaderSection(isConsistent, sd, timezone, vecParamInUse);
			}
735

736
737
			//3. write data depending on ASCII/BINARY
			if (outputIsAscii) {
738
				writeDataAscii(isConsistent, vecMeteo[ii], vecParamInUse);
739
740
741
742
743
			} 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);
744
				writeDataBinary(isConsistent, vecMeteo[ii], vecParamInUse);
745
			}
746

747
748
749
750
751
752
753
754
755
756
			cleanup();

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

757
void SMETIO::writeDataBinary(const bool& writeLocationInHeader, const std::vector<MeteoData>& vecMeteo,
758
                             const std::vector<bool>& vecParamInUse)
759
{
760
	const char eoln = '\n';
761

762
	for (size_t ii=0; ii<vecMeteo.size(); ii++){
763
		float val = 0;
764
765
766
767
768
769
770
771

		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();
772
		}
773
		fout.write((char*)&julian, sizeof(double));
774
775

		if (!writeLocationInHeader){ //Meta data changes
776
			val = (float)vecMeteo[ii].meta.position.getLat();
777
			fout.write((char*)&val, sizeof(float));
778
			val = (float)vecMeteo[ii].meta.position.getLon();
779
			fout.write((char*)&val, sizeof(float));
780
			val = (float)vecMeteo[ii].meta.position.getAltitude();
781
782
783
			fout.write((char*)&val, sizeof(float));
		}

784
		for (size_t jj=0; jj<MeteoData::nrOfParameters; jj++){
785
786
787
788
789
790
791
792
793
			if (vecParamInUse[jj]){
				val = (float)vecMeteo[ii].param(jj);
				fout.write((char*)&val, sizeof(float));
			}
		}
		fout.write((char*)&eoln, sizeof(char));
	}
}

794
void SMETIO::writeDataAscii(const bool& writeLocationInHeader, const std::vector<MeteoData>& vecMeteo,
795
                            const std::vector<bool>& vecParamInUse)
796
797
{
	fout << "[DATA]" << endl;
798
799
	fout.fill(' ');
	fout << right;
800
	fout << fixed;
801
	for (size_t ii=0; ii<vecMeteo.size(); ii++){
802
803
804
805
806
807
808
		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);
		}
809
810

		if (!writeLocationInHeader){ //Meta data changes
811
812
813
			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();
814
815
		}

816
		for (size_t jj=0; jj<MeteoData::nrOfParameters; jj++){
817
818
819
820
821
822
823
			fout << " ";
			if (vecParamInUse[jj]){
				setFormatting(MeteoData::Parameters(jj));
				if (vecMeteo[ii].param(jj) == IOUtils::nodata)
					fout << setprecision(0);
				fout << vecMeteo[ii].param(jj);
			}
824
825
826
827
828
		}
		fout << endl;
	}
}

829
void SMETIO::setFormatting(const MeteoData::Parameters& paramindex)
830
831
832
{
	if ((paramindex == MeteoData::TA) || (paramindex == MeteoData::TSS) || (paramindex == MeteoData::TSG))
		fout << setw(8) << setprecision(2);
833
	else if ((paramindex == MeteoData::VW) || (paramindex == MeteoData::VW_MAX))
834
835
836
837
838
839
840
841
842
843
844
845
846
		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);
}

847
void SMETIO::writeHeaderSection(const bool& writeLocationInHeader, const StationData& sd,
848
                                const double& timezone, const std::vector<bool>& vecParamInUse)
849
850
{
	fout << "[HEADER]" << endl;
851
	fout << "station_id   = " << sd.getStationID() << endl;
852
853
	if (sd.getStationName() != "")
		fout << "station_name = " << sd.getStationName() << endl;
854

855
	fout << fixed;
856
	if (writeLocationInHeader){ //TODO: only write if != nodata
857
858
859
860
861
862
		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";
863
864
	}

865
	fout << "nodata       = " << setw(7) << setprecision(0) << IOUtils::nodata << "\n";
866

867
	if ((timezone != IOUtils::nodata) && (timezone != 0.0))
868
		fout << "tz           = " << setw(7)  << setprecision(0) << timezone << "\n";
869

870
	fout << "fields       = timestamp";
871
872
873
874
875

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

876
	for (size_t ii=0; ii<MeteoData::nrOfParameters; ii++){
877
878
879
880
881
882
		if (vecParamInUse[ii]) {
			std::string column=MeteoData::getParameterName(ii);
			if(column=="RSWR") column="OSWR";
			if(column=="HNW") column="PSUM";
			fout << " " << column;
		}
883
884
885
886
887
	}
	fout << endl;
}


888
void SMETIO::checkForUsedParameters(const std::vector<MeteoData>& vecMeteo, double& timezone,
889
                                    std::vector<bool>& vecParamInUse)
890
{
891
892
	for (size_t ii=0; ii<vecMeteo.size(); ii++){
		for (size_t jj=0; jj<MeteoData::nrOfParameters; jj++){
893
894
895
896
			if (!vecParamInUse[jj])
				if (vecMeteo[ii].param(jj) != IOUtils::nodata)
					vecParamInUse[jj] = true;
		}
897
	}
898
899
900
901
902

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

903
bool SMETIO::checkConsistency(const std::vector<MeteoData>& vecMeteo, StationData& sd)
904
{
905
906
	if (vecMeteo.size() > 0) //HACK to get the station data even when encoutering bug 87
		sd = vecMeteo[0].meta;
907

908
	for (size_t ii=1; ii<vecMeteo.size(); ii++){
909
910
911
912
913
914
		const Coords& p1 = vecMeteo[ii-1].meta.position;
		const Coords& p2 = vecMeteo[ii].meta.position;
		if (p1 != p2) {
			//we don't mind if p1==nodata or p2==nodata
			if(p1.isNodata()==false && p2.isNodata()==false) return false;
		}
915
916
917
918
919
	}

	return true;
}

920
void SMETIO::readSpecialPoints(std::vector<Coords>&)
921
922
923
924
925
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

926
void SMETIO::write2DGrid(const Grid2DObject& /*grid_in*/, const std::string& /*name*/)
927
928
929
930
931
932
933
934
935
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}


#ifndef _METEOIO_JNI
extern "C"
{
936
937
938
#define COMPILE_PLUGIN
#include "exports.h"

939
	METEOIO_EXPORT void deleteObject(void* obj) {
940
941
942
		delete reinterpret_cast<PluginObject*>(obj);
	}

943
	METEOIO_EXPORT void* loadObject(const string& classname, const Config& cfg) {
944
		if(classname == "SMETIO") {
945
			//cerr << "Creating dynamic handle for " << classname << endl;
946
			return new SMETIO(deleteObject, cfg);
947
948
949
950
951
952
953
954
		}
		//cerr << "Could not load " << classname << endl;
		return NULL;
	}
}
#endif

} //namespace