WSL/SLF GitLab Repository

SMETIO.cc 27.6 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 (including the list of standard parameters: TA, TSS, TSG, RH, VW, DW, ISWR, OSWR, ILWR, OLWR, PINT, PSUM, HS).
28
 * This plugin can also provide special points, given as a SMET file containing either latitude/longitude/altitude or easting/northing/altitude. For the latter, the header must contain the epsg code (see example below).
29
 *
30
31
 * Non-standard parameters can also be given, such as extra snow temperatures. These parameters will then take the name that has been given in "fields", converted to uppercase. It is usually a good idea to number these parameters, such as TS1, TS2, TS3 for a serie of temperatures at various positions.
 *
32
 * @section template_units Units
33
 * 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).
34
35
36
 *
 * @section template_keywords Keywords
 * This plugin uses the following keywords:
37
38
 * - 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
39
 * - METEOPARAM: output file format options (ASCII or BINARY that might be followed by GZIP)
40
 * - SPECIALPTSFILE: a path+file name to the a file containing grid coordinates of special points of interest (for special outputs)
41
42
43
44
 *
 * Example:
 * @code
 * [Input]
45
 * METEO = SMET
46
47
48
49
 * METEOPATH = ./input
 * STATION1 = uppper_station.smet
 * STATION2 = lower_station.smet
 * STATION3 = outlet_station.smet
50
51
52
53
 * [Output]
 * METEOPATH = ./output
 * METEOPARAM = ASCII GZIP
 * @endcode
54
55
56
57
58
59
60
61
62
63
64
65
66
67
 *
 * Below is an example of special points input:
 * @code
 * SMET 1.1 ASCII
 * [HEADER]
 * station_id	= my_pts
 * epsg	= 21781
 * nodata	= -999
 * fields = easting northing altitude
 * [DATA]
 * 832781 187588 2115
 * 635954 80358 2428
 * @endcode
 *
68
69
 */

70
71
const std::string SMETIO::dflt_extension = ".smet";

72
SMETIO::SMETIO(const std::string& configfile)
73
        : cfg(configfile),
74
75
76
          coordin(), coordinparam(), coordout(), coordoutparam(),
          vec_smet_reader(), vecFiles(), outpath(), in_dflt_TZ(0.), out_dflt_TZ(0.),
          plugin_nodata(IOUtils::nodata), nr_stations(0), outputIsAscii(true), outputIsGzipped(false)
77
78
79
80
{
	parseInputOutputSection();
}

81
SMETIO::SMETIO(const Config& cfgreader)
82
        : cfg(cfgreader),
83
84
85
          coordin(), coordinparam(), coordout(), coordoutparam(),
          vec_smet_reader(), vecFiles(), outpath(), in_dflt_TZ(0.), out_dflt_TZ(0.),
          plugin_nodata(IOUtils::nodata), nr_stations(0), outputIsAscii(true), outputIsGzipped(false)
86
87
88
89
{
	parseInputOutputSection();
}

90
SMETIO::~SMETIO() throw()
91
92
93
94
{

}

95
void SMETIO::read2DGrid(Grid2DObject& /*grid_out*/, const std::string& /*_name*/)
96
97
98
99
100
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

101
102
103
104
105
106
void SMETIO::read2DGrid(Grid2DObject&, const MeteoGrids::Parameters&, const Date&)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

107
void SMETIO::readDEM(DEMObject& /*dem_out*/)
108
109
110
111
112
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

113
void SMETIO::readLanduse(Grid2DObject& /*landuse_out*/)
114
115
116
117
118
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

119
void SMETIO::readAssimilationData(const Date& /*date_in*/, Grid2DObject& /*da_out*/)
120
121
122
123
124
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

125
void SMETIO::readStationData(const Date&, std::vector<StationData>& vecStation)
126
{//HACK: It should support coordinates in the data, ie: it should use the given date! (and TZ)
127
	vecStation.clear();
128
	vecStation.reserve(nr_stations);
129
130

	//Now loop through all requested stations, open the respective files and parse them
131
	for (size_t ii=0; ii<vec_smet_reader.size(); ii++){
132
		StationData sd;
133
		smet::SMETReader& myreader = vec_smet_reader[ii];
134

135
		read_meta_data(myreader, sd);
136
137
		vecStation.push_back(sd);
	}
138
139
}

140
void SMETIO::parseInputOutputSection()
141
{
142
	//default timezones
143
	in_dflt_TZ = out_dflt_TZ = IOUtils::nodata;
144
145
	cfg.getValue("TIME_ZONE","Input",in_dflt_TZ,IOUtils::nothrow);
	cfg.getValue("TIME_ZONE","Output",out_dflt_TZ,IOUtils::nothrow);
146
147

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

	//Parse input section: extract number of files to read and store filenames in vecFiles
151
	std::string inpath, in_meteo;
152
	cfg.getValue("METEO", "Input", in_meteo, IOUtils::nothrow);
153
	if (in_meteo == "SMET") { //keep it synchronized with IOHandler.cc for plugin mapping!!
154
		cfg.getValue("METEOPATH", "Input", inpath);
155
		size_t counter = 1;
156
		string filename;
157

158
159
		do {
			stringstream ss;
160
			filename.clear();
161

162
			ss << "STATION" << counter;
163
			cfg.getValue(ss.str(), "Input", filename, IOUtils::nothrow);
164

165
			if (!filename.empty()){
166
				if(IOUtils::getExtension(filename)=="") filename += dflt_extension; //default extension
167
				ostringstream file_and_path;
168
169
170
				file_and_path << inpath << "/" << filename;
				if (!IOUtils::validFileName(file_and_path.str())) //Check whether filename is valid
					throw InvalidFileNameException(file_and_path.str(), AT);
171

172
173
174
				vecFiles.push_back(file_and_path.str());
			}
			counter++;
175
		} while (!filename.empty());
176

177
		nr_stations = counter - 1;
178

179
180
181
		for (size_t ii=0; ii<vecFiles.size(); ii++){
			vec_smet_reader.push_back(smet::SMETReader(vecFiles[ii]));
		}
182
	}
183

184
	//Parse output section: extract info on whether to write ASCII or BINARY format, gzipped or not
185
	outpath.clear();
186
	outputIsAscii = true;
187
188
189
	outputIsGzipped = false;

	vector<string> vecArgs;
190
191
	cfg.getValue("METEOPATH", "Output", outpath, IOUtils::nothrow);
	cfg.getValue("METEOPARAM", "Output", vecArgs, IOUtils::nothrow); //"ASCII|BINARY GZIP"
192

193
	if (outpath.empty())
194
195
		return;

196
	if (vecArgs.empty())
197
198
199
200
201
202
203
204
205
		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;
206
	else
207
208
209
210
211
212
213
214
		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;
	}
215

216
217
}

218
void SMETIO::identify_fields(const std::vector<std::string>& fields, std::vector<size_t>& indexes,
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
                             bool& julian_present, MeteoData& md)
{
	/*
	 * This function associates a parameter index for MeteoData objects with the
	 * lineup of field types in a SMET header. The following SMET fields are treated
	 * exceptionally:
	 * - julian, associated with IOUtils::npos
	 * - latitude, associated with IOUtils::npos-1
	 * - longitude, associated with IOUtils::npos-2
	 * - easting, associated with IOUtils::npos-3
	 * - norhting, associated with IOUtils::npos-4
	 * - altitude, associated with IOUtils::npos-5
	 * If a paramter is unknown in the fields section, then it is added as separate field to MeteoData
	 */
	for (size_t ii=0; ii<fields.size(); ii++){
234
235
236
237
238
239
240
241
242
243
		const string& key = fields[ii];

		if(md.param_exists(key)) {
			indexes.push_back(md.getParameterIndex(key));
			continue;
		}

		//specific key mapping
		if (key == "PSUM") {
			indexes.push_back(md.getParameterIndex("HNW"));
244
245
		} else if (key == "OSWR") {
			indexes.push_back(md.getParameterIndex("RSWR"));
246
		} else if (key == "OLWR") {
247
248
			md.addParameter("OLWR");
			indexes.push_back(md.getParameterIndex("OLWR"));
249
		} else if (key == "PINT") {
250
251
			md.addParameter("PINT");
			indexes.push_back(md.getParameterIndex("PINT"));
252
		} else if (key == "julian") {
253
254
			julian_present = true;
			indexes.push_back(IOUtils::npos);
255
		} else if (key == "latitude") {
256
			indexes.push_back(IOUtils::npos-1);
257
		} else if (key == "longitude") {
258
			indexes.push_back(IOUtils::npos-2);
259
		} else if (key == "easting") {
260
			indexes.push_back(IOUtils::npos-3);
261
		} else if (key == "northing") {
262
			indexes.push_back(IOUtils::npos-4);
263
		} else if (key == "altitude") {
264
265
			indexes.push_back(IOUtils::npos-5);
		} else {
266
			//this is an extra parameter, we convert to uppercase
267
			std::string extra_param = key;
268
269
270
			IOUtils::toUpper(extra_param);
			md.addParameter(extra_param);
			indexes.push_back(md.getParameterIndex(extra_param));
271
272
273
274
275
276
277
278
279
280
281
		}
	}
}

void SMETIO::read_meta_data(const smet::SMETReader& myreader, StationData& meta)
{
	/*
	 * This function reads in the header data provided by a SMETReader object.
	 * SMETReader objects read all the header info upon construction and can subsequently
	 * be queried for that info
	 */
282
	const double nodata_value = myreader.get_header_doublevalue("nodata");
283
284
285

	meta.position.setProj(coordin, coordinparam); //set the default projection from config file
	if (myreader.location_in_header(smet::WGS84)){
286
287
288
		const double lat = myreader.get_header_doublevalue("latitude");
		const double lon = myreader.get_header_doublevalue("longitude");
		const double alt = myreader.get_header_doublevalue("altitude");
289
		meta.position.setLatLon(lat, lon, alt);
290
	}
291
292

	if (myreader.location_in_header(smet::EPSG)){
293
294
295
296
		const double east  = myreader.get_header_doublevalue("easting");
		const double north = myreader.get_header_doublevalue("northing");
		const double alt   = myreader.get_header_doublevalue("altitude");
		const short int epsg  = (short int)(floor(myreader.get_header_doublevalue("epsg") + 0.1));
297
		meta.position.setEPSG(epsg); //this needs to be set before calling setXY(...)
298
299
		meta.position.setXY(east, north, alt);
	}
300

301
302
303
	meta.stationID = myreader.get_header_value("station_id");
	meta.stationName = myreader.get_header_value("station_name");

304
	const bool data_epsg = myreader.location_in_data(smet::EPSG);
305
	if (data_epsg){
306
307
		const double d_epsg = myreader.get_header_doublevalue("epsg");
		const short int epsg = (d_epsg != nodata_value)? (short int)(floor(d_epsg + 0.1)): IOUtils::snodata;
308
309
310
311
		meta.position.setEPSG(epsg);
	}
}

312
void SMETIO::copy_data(const smet::SMETReader& myreader,
313
314
                       const std::vector<std::string>& timestamps,
                       const std::vector<double>& mydata, std::vector<MeteoData>& vecMeteo)
315
316
317
318
319
320
{
	/*
	 * This function parses the data read from a SMETReader object, a vector<double>,
	 * and copies the values into their respective places in the MeteoData structure
	 * Meta data, whether in header or in data is also handled
	 */
321
	const string myfields = myreader.get_header_value("fields");
322
323
	vector<string> fields;
	IOUtils::readLineToVec(myfields, fields);
324
325
326
327

	bool julian_present = false;
	MeteoData md;
	vector<size_t> indexes;
328
329
	identify_fields(fields, indexes, julian_present, md);

330
	if ((timestamps.empty()) && (!julian_present)) return; //nothing to do
331

332
333
334
	const bool olwr_present = md.param_exists("OLWR");
	const bool data_wgs84 = myreader.location_in_data(smet::WGS84);
	const bool data_epsg = myreader.location_in_data(smet::EPSG);
335
336
337

	read_meta_data(myreader, md.meta);

338
	const double nodata_value = myreader.get_header_doublevalue("nodata");
339
	double current_timezone = myreader.get_header_doublevalue("tz");
340
341
	if (current_timezone == nodata_value)
		current_timezone = in_dflt_TZ;
342
	const bool timestamp_present = myreader.contains_timestamp();
343

344
345
	const size_t nr_of_fields = indexes.size();
	const size_t nr_of_lines = mydata.size() / nr_of_fields;
346
347
348
349
350
351

	double lat=IOUtils::nodata, lon=IOUtils::nodata, east=IOUtils::nodata, north=IOUtils::nodata, alt=IOUtils::nodata;
	size_t current_index = 0; //index to vec_data
	for (size_t ii = 0; ii<nr_of_lines; ii++){
		vecMeteo.push_back(md);
		MeteoData& tmp_md = vecMeteo.back();
352

353
354
355
		if (timestamp_present)
			IOUtils::convertString(tmp_md.date, timestamps[ii], current_timezone);

356
		//Copy data points
357
		for (size_t jj=0; jj<nr_of_fields; jj++){
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
			const double& current_data = mydata[current_index];
			if (indexes[jj] >= IOUtils::npos-5){ //the special fields have high indexes
				if (indexes[jj] == IOUtils::npos){
					if (!timestamp_present){
						if (current_data != nodata_value)
							tmp_md.date.setDate(current_data, current_timezone);
					}
				} else if (indexes[jj] == IOUtils::npos-1){
					lat = current_data;
				} else if (indexes[jj] == IOUtils::npos-2){
					lon = current_data;
				} else if (indexes[jj] == IOUtils::npos-3){
					east = current_data;
				} else if (indexes[jj] == IOUtils::npos-4){
					north = current_data;
				} else if (indexes[jj] == IOUtils::npos-5){
					alt = current_data;
				}
			} else {
				if (current_data == nodata_value)
378
					tmp_md(indexes[jj]) = IOUtils::nodata;
379
				else
380
					tmp_md(indexes[jj]) = current_data;
381
382
383
384
385
386
387
388
389
390
			}

			if (data_epsg)
				tmp_md.meta.position.setXY(east, north, alt);

			if (data_wgs84)
				tmp_md.meta.position.setXY(lat, lon, alt);

			current_index++;
		}
391

392
393
394
		if ((olwr_present) && (tmp_md(MeteoData::TSS) == IOUtils::nodata)) {//HACK
			tmp_md(MeteoData::TSS) = olwr_to_tss(tmp_md("OLWR"));
		}
395
396
397
	}
}

398
399
400
double SMETIO::olwr_to_tss(const double& olwr) {
	const double ea = 1.;
	if(olwr==IOUtils::nodata) return IOUtils::nodata;
401
402
	if(olwr<0.) return IOUtils::nodata; //since olwr is NOT filtered, making sure no arithmetic exception would happen
	if(olwr>1e4) return IOUtils::nodata; //since olwr is NOT filtered, making sure no arithmetic exception would happen
403
	return pow( olwr / ( ea * Cst::stefan_boltzmann ), 0.25);
404
405
}

406
void SMETIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
407
                           std::vector< std::vector<MeteoData> >& vecMeteo,
408
                           const size_t& stationindex)
409
{
410
	//Make sure that vecMeteo have the correct dimension and stationindex is valid
411
	size_t startindex=0, endindex=vecFiles.size();
412
	if (stationindex != (size_t)IOUtils::npos){ //HACK do we really still need stationindex??
413
		if ((stationindex < vecFiles.size()) || (stationindex < vecMeteo.size())){
414
415
416
417
			startindex = stationindex;
			endindex = stationindex+1;
		} else {
			throw IndexOutOfBoundsException("Invalid stationindex", AT);
418
		}
419
420
421

		vecMeteo[stationindex].clear();
	} else {
422
		vecMeteo.clear();
423
		vecMeteo = vector< vector<MeteoData> >(vecFiles.size());
424
		vecMeteo.reserve(nr_stations);
425
426
427
	}

	//Now loop through all requested stations, open the respective files and parse them
428
	for (size_t ii=startindex; ii<endindex; ii++){
429
		const string& filename = vecFiles.at(ii); //filename of current station
430

431
432
433
		if (!IOUtils::fileExists(filename))
			throw FileNotFoundException(filename, AT);

434
435
		smet::SMETReader& myreader = vec_smet_reader.at(ii);
		myreader.convert_to_MKSA(true); // we want converted values for MeteoIO
436

437
438
		vector<double> mydata; //sequentially store all data in the smet file
		vector<string> mytimestamps;
439

440
441
442
		if (myreader.contains_timestamp()){
			myreader.read(dateStart.toString(Date::ISO), dateEnd.toString(Date::ISO), mytimestamps, mydata);
		} else {
443
			myreader.read(dateStart.getJulian(), dateEnd.getJulian(), mydata);
444
		}
445

446
		copy_data(myreader, mytimestamps, mydata, vecMeteo[ii]);
447
	}
448
449
}

450
void SMETIO::writeMeteoData(const std::vector< std::vector<MeteoData> >& vecMeteo, const std::string&)
451
452
{
	//Loop through all stations
453
	for (size_t ii=0; ii<vecMeteo.size(); ii++){
454
		//1. check consistency of station data position -> write location in header or data section
455
		StationData sd;
456
		sd.position.setProj(coordout, coordoutparam);
457
		const bool isConsistent = checkConsistency(vecMeteo.at(ii), sd);
458

459
		if (sd.stationID.empty()){
460
461
462
463
464
			stringstream ss;
			ss << "Station" << ii+1;
			sd.stationID = ss.str();
		}

465
		const string filename = outpath + "/" + sd.stationID + ".smet";
466
467
468
		if (!IOUtils::validFileName(filename)) //Check whether filename is valid
			throw InvalidFileNameException(filename, AT);

469
		//2. check which meteo parameter fields are actually in use
470
		const size_t nr_of_parameters = getNrOfParameters(sd.stationID, vecMeteo[ii]);
471
472
		vector<bool> vecParamInUse = vector<bool>(nr_of_parameters, false);
		vector<string> vecColumnName = vector<string>(nr_of_parameters, "NULL");
473
		double timezone = IOUtils::nodata; //time zone of the data
474
		checkForUsedParameters(vecMeteo[ii], nr_of_parameters, timezone, vecParamInUse, vecColumnName);
475
		if(out_dflt_TZ != IOUtils::nodata) timezone=out_dflt_TZ; //if the user set an output time zone, all will be converted to it
476

477
		try {
478
			const smet::SMETType type = (outputIsAscii)? smet::ASCII : smet::BINARY;
479

480
			smet::SMETWriter mywriter(filename, type, outputIsGzipped);
481
			generateHeaderInfo(sd, outputIsAscii, isConsistent, timezone,
482
                               nr_of_parameters, vecParamInUse, vecColumnName, mywriter);
483
484
485

			vector<string> vec_timestamp;
			vector<double> vec_data;
486
			for (size_t jj=0; jj<vecMeteo[ii].size(); jj++) {
487
				if (outputIsAscii){
488
					if (out_dflt_TZ != IOUtils::nodata) { //user-specified time zone
489
490
491
492
493
						Date tmp_date(vecMeteo[ii][jj].date);
						tmp_date.setTimeZone(out_dflt_TZ);
						vec_timestamp.push_back(tmp_date.toString(Date::ISO));
					} else {
						vec_timestamp.push_back(vecMeteo[ii][jj].date.toString(Date::ISO));
494
					}
495
496
497
498
499
				} else {
					double julian;
					if(out_dflt_TZ!=IOUtils::nodata) {
						Date tmp_date(vecMeteo[ii][jj].date);
						tmp_date.setTimeZone(out_dflt_TZ);
500
						julian = tmp_date.getJulian();
501
					} else {
502
						julian = vecMeteo[ii][jj].date.getJulian();
503
504
505
					}
					vec_data.push_back(julian);
				}
506

507
				if (!isConsistent) { //Meta data changes
508
509
510
511
					vec_data.push_back(vecMeteo[ii][jj].meta.position.getLat());
					vec_data.push_back(vecMeteo[ii][jj].meta.position.getLon());
					vec_data.push_back(vecMeteo[ii][jj].meta.position.getAltitude());
				}
512

513
				for (size_t kk=0; kk<nr_of_parameters; kk++) {
514
					if (vecParamInUse[kk])
515
						vec_data.push_back(vecMeteo[ii][jj](kk)); //add data value
516
				}
517
			}
518

519
520
521
			if (outputIsAscii) mywriter.write(vec_timestamp, vec_data);
			else mywriter.write(vec_data);

522
		} catch(exception&) {
523
			throw;
524
525
526
527
		}
	}
}

528
void SMETIO::generateHeaderInfo(const StationData& sd, const bool& i_outputIsAscii, const bool& isConsistent,
529
530
531
                                const double& timezone, const size_t& nr_of_parameters,
                                const std::vector<bool>& vecParamInUse, const std::vector<std::string>& vecColumnName,
                                smet::SMETWriter& mywriter)
532
533
534
535
536
537
538
539
540
541
542
543
544
{
	/**
	 * This procedure sets all relevant information for the header in the SMETWriter object mywriter
	 * The following key/value pairs are set for the header:
	 * - station_id, station_name (if present)
	 * - nodata (set to IOUtils::nodata)
	 * - fields (depending on ASCII/BINARY format and whether the meta data is part of the header or data)
	 * - timezone
	 * - meta data (lat/lon/alt or east/north/alt/epsg if not part of data section)
	 */
	stringstream ss;

	mywriter.set_header_value("station_id", sd.stationID);
545
	if (!sd.stationName.empty())
546
547
548
549
550
		mywriter.set_header_value("station_name", sd.stationName);
	mywriter.set_header_value("nodata", IOUtils::nodata);

	vector<int> myprecision, mywidth; //set meaningful precision/width for each column

551
	if (i_outputIsAscii) {
552
553
554
555
556
557
558
559
560
561
562
563
564
565
		ss << "timestamp";
	} else {
		ss << "julian";
		myprecision.push_back(8);
		mywidth.push_back(16);
	}

	if (isConsistent) {
		mywriter.set_header_value("latitude", sd.position.getLat());
		mywriter.set_header_value("longitude", sd.position.getLon());
		mywriter.set_header_value("easting", sd.position.getEasting());
		mywriter.set_header_value("northing", sd.position.getNorthing());
		mywriter.set_header_value("altitude", sd.position.getAltitude());
		mywriter.set_header_value("epsg", (double)sd.position.getEPSG());
566

567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
		if ((timezone != IOUtils::nodata) && (timezone != 0.0))
			mywriter.set_header_value("tz", timezone);
	} else {
		ss << " latitude longitude altitude";
		myprecision.push_back(8); //for latitude
		mywidth.push_back(11);    //for latitude
		myprecision.push_back(8); //for longitude
		mywidth.push_back(11);    //for longitude
		myprecision.push_back(1); //for altitude
		mywidth.push_back(7);     //for altitude
	}

	//Add all other used parameters
	int tmpwidth, tmpprecision;
	for (size_t ll=0; ll<nr_of_parameters; ll++) {
		if (vecParamInUse[ll]) {
			string column = vecColumnName.at(ll);
			if (column == "RSWR") column = "OSWR";
			if (column == "HNW")  column = "PSUM";
			ss << " " << column;
587

588
589
590
591
592
593
594
595
596
597
598
			getFormatting(ll, tmpprecision, tmpwidth);
			myprecision.push_back(tmpprecision);
			mywidth.push_back(tmpwidth);
		}
	}

	mywriter.set_header_value("fields", ss.str());
	mywriter.set_width(mywidth);
	mywriter.set_precision(myprecision);
}

599
void SMETIO::getFormatting(const size_t& param, int& prec, int& width)
600
{
601
602
	/**
	 * When writing a SMET file, different meteo parameters require a different
603
604
	 * format with regard to precision and width when printing.
	 * This procedure sets the precision and width for each known parameter and
605
606
	 * defaults to a width of 8 and precision of 3 digits for each unknown parameter.
	 */
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
	if ((param == MeteoData::TA) || (param == MeteoData::TSS) || (param == MeteoData::TSG)){
		prec = 2;
		width = 8;
	} else if ((param == MeteoData::VW) || (param == MeteoData::VW_MAX)){
		prec = 1;
		width = 6;
	} else if (param == MeteoData::DW){
		prec = 0;
		width = 5;
	} else if ((param == MeteoData::ISWR) || (param == MeteoData::RSWR) || (param == MeteoData::ILWR)){
		prec = 0;
		width = 6;
	} else if (param == MeteoData::HNW){
		prec = 3;
		width = 6;
	} else if (param == MeteoData::HS){
		prec = 3;
		width = 8;
	} else if (param == MeteoData::RH){
		prec = 3;
		width = 7;
	} else {
		prec = 3;
		width = 8;
631
632
633
	}
}

634
635
636
size_t SMETIO::getNrOfParameters(const std::string& stationname, const std::vector<MeteoData>& vecMeteo)
{
	/**
637
	 * This function loops through all MeteoData objects present in vecMeteo and returns the
638
639
640
641
642
643
644
645
646
647
648
	 * number of meteo parameters that the MeteoData objects have. If there is an inconsistency
	 * in the number of meteo parameters in use within the vector of MeteoData then a warning
	 * is printed and MeteoData::nrOfParameters is returned, thus all additional meteo parameters
	 * that might be in use are ignored.
	 */

	size_t actual_nr_of_parameters = IOUtils::npos;
	bool has_one_element = false;

	for (size_t ii=0; ii<vecMeteo.size(); ii++){
		has_one_element = true;
649
		const size_t current_size = vecMeteo[ii].getNrOfParameters();
650
651
652
653
654

		if (actual_nr_of_parameters == IOUtils::npos){
			actual_nr_of_parameters = current_size;
		} else if (actual_nr_of_parameters != current_size){
			//There is an inconsistency in the fields, print out a warning and proceed
655
			cerr << "[w] While writing SMET file: Inconsistency in number of meteo "
656
657
658
659
660
661
662
663
664
665
666
667
668
669
				<< "parameters for station " << stationname << endl;
			actual_nr_of_parameters = MeteoData::nrOfParameters;
			break;
		}
	}

	if (!has_one_element)
		return MeteoData::nrOfParameters;

	return actual_nr_of_parameters;
}

void SMETIO::checkForUsedParameters(const std::vector<MeteoData>& vecMeteo, const size_t& nr_parameters, double& timezone,
                                    std::vector<bool>& vecParamInUse, std::vector<std::string>& vecColumnName)
670
{
671
672
673
674
675
676
	/**
	 * This procedure loops through all MeteoData objects present in vecMeteo and finds out which
	 * meteo parameters are actually in use, i. e. have at least one value that differs from IOUtils::nodata.
	 * If a parameter is in use, then vecParamInUse[index_of_parameter] is set to true and the column
	 * name is set in vecColumnName[index_of_parameter]
	 */
677
	for (size_t ii=0; ii<vecMeteo.size(); ii++){
678
679
		for (size_t jj=0; jj<nr_parameters; jj++){
			if (!vecParamInUse[jj]){
680
				if (vecMeteo[ii](jj) != IOUtils::nodata){
681
					vecParamInUse[jj] = true;
682
683
684
					vecColumnName.at(jj) = vecMeteo[ii].getNameForParameter(jj);
				}
			}
685
		}
686
	}
687

688
	if (!vecMeteo.empty())
689
690
691
		timezone = vecMeteo[0].date.getTimeZone();
}

692
bool SMETIO::checkConsistency(const std::vector<MeteoData>& vecMeteo, StationData& sd)
693
{
694
695
696
697
698
699
	/**
	 * This function checks whether all the MeteoData elements in vecMeteo are consistent
	 * regarding their meta data (position information, station name). If they are consistent
	 * true is returned, otherwise false
	 */

700
	if (!vecMeteo.empty()) //to get the station data even when in bug 87 conditions
701
		sd = vecMeteo[0].meta;
702

703
	for (size_t ii=1; ii<vecMeteo.size(); ii++){
704
705
706
707
708
709
		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;
		}
710
711
712
713
714
	}

	return true;
}

715
void SMETIO::readSpecialPoints(std::vector<Coords>& pts)
716
{
717
	std::string filename;
718
719
720
721
722
723
724
725
726
727
728
729
	cfg.getValue("SPECIALPTSFILE", "Input", filename);
	if (!IOUtils::fileExists(filename)) {
		throw FileNotFoundException(filename, AT);
	}

	smet::SMETReader myreader(filename);
	vector<double> vec_data;
	myreader.read(vec_data);
	const size_t nr_fields = myreader.get_nr_of_fields();
	const int epsg = myreader.get_header_intvalue("epsg");
	const double smet_nodata = myreader.get_header_doublevalue("nodata");

730
	if (myreader.location_in_data(smet::WGS84)==true) {
731
732
		size_t lat_fd=IOUtils::unodata, lon_fd=IOUtils::unodata;
		size_t alt_fd=IOUtils::unodata;
733
		for (size_t ii=0; ii<nr_fields; ii++) {
734
			const string tmp = myreader.get_field_name(ii);
735
736
737
			if (tmp=="latitude") lat_fd=ii;
			if (tmp=="longitude") lon_fd=ii;
			if (tmp=="altitude") alt_fd=ii;
738
		}
739
		for (size_t ii=0; ii<vec_data.size(); ii+=nr_fields) {
740
741
742
743
			Coords point;
			point.setLatLon(vec_data[ii+lat_fd], vec_data[ii+lon_fd], vec_data[ii+alt_fd]);
			pts.push_back(point);
		}
744
745
	} else if (myreader.location_in_data(smet::EPSG)==true) {
		if (epsg==(int)floor(smet_nodata + 0.1))
746
747
748
749
			throw InvalidFormatException("In file \""+filename+"\", missing EPSG code in header!", AT);

		size_t east_fd=IOUtils::unodata, north_fd=IOUtils::unodata;
		size_t alt_fd=IOUtils::unodata;
750
		for (size_t ii=0; ii<nr_fields; ii++) {
751
			const string tmp = myreader.get_field_name(ii);
752
753
754
			if (tmp=="easting") east_fd=ii;
			if (tmp=="northing") north_fd=ii;
			if (tmp=="altitude") alt_fd=ii;
755
		}
756
757
758
759
		if ((east_fd == IOUtils::unodata) || (north_fd == IOUtils::unodata) || (alt_fd == IOUtils::unodata))
			throw InvalidFormatException("File \""+filename+"\" does not contain all data fields necessary for EPSG coordinates", AT);

		for (size_t ii=0; ii<vec_data.size(); ii+=nr_fields) {
760
761
762
763
764
765
766
767
			Coords point;
			point.setEPSG(epsg);
			point.setXY(vec_data[ii+east_fd], vec_data[ii+north_fd], vec_data[ii+alt_fd]);
			pts.push_back(point);
		}
	} else {
		throw InvalidFormatException("File \""+filename+"\" does not contain expected location information in DATA section!", AT);
	}
768
769
}

770
void SMETIO::write2DGrid(const Grid2DObject& /*grid_in*/, const std::string& /*name*/)
771
772
773
774
775
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

776
777
778
779
780
781
void SMETIO::write2DGrid(const Grid2DObject&, const MeteoGrids::Parameters&, const Date&)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

782
} //namespace