WSL/SLF GitLab Repository

SMETIO.cc 27.3 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 smetio_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 Points Of Interest, 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 smetio_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
 * @section smetio_keywords Keywords
36
 * 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
 * - POIFILE: a path+file name to the a file containing grid coordinates of 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
 * Below is an example of Points Of Interest input:
56
57
58
59
60
61
62
63
64
65
66
67
 * @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
156
157
158
159
160
161
162
163
164
165
166
		std::vector<std::string> vecFilenames;
		cfg.getValues("STATION", "INPUT", vecFilenames);

		for (size_t ii=0; ii<vecFilenames.size(); ii++) {
			const string filename = vecFilenames[ii];
			const string extension = IOUtils::getExtension(filename);
			const std::string file_and_path = (extension!="")? inpath+"/"+filename : inpath+"/"+filename+dflt_extension;

			if (!IOUtils::validFileName(file_and_path)) //Check whether filename is valid
				throw InvalidFileNameException(file_and_path, AT);
			vecFiles.push_back(file_and_path);
			vec_smet_reader.push_back(smet::SMETReader(file_and_path));
167
		}
168
	}
169

170
	//Parse output section: extract info on whether to write ASCII or BINARY format, gzipped or not
171
	outpath.clear();
172
	outputIsAscii = true;
173
174
175
	outputIsGzipped = false;

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

179
	if (outpath.empty())
180
181
		return;

182
	if (vecArgs.empty())
183
184
185
186
187
188
189
190
191
		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;
192
	else
193
194
195
196
197
198
199
200
		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;
	}
201

202
203
}

204
void SMETIO::identify_fields(const std::vector<std::string>& fields, std::vector<size_t>& indexes,
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
                             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++){
220
221
222
223
224
225
226
227
228
229
		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"));
230
231
		} else if (key == "OSWR") {
			indexes.push_back(md.getParameterIndex("RSWR"));
232
		} else if (key == "OLWR") {
233
234
			md.addParameter("OLWR");
			indexes.push_back(md.getParameterIndex("OLWR"));
235
		} else if (key == "PINT") { //in mm/h
236
237
			md.addParameter("PINT");
			indexes.push_back(md.getParameterIndex("PINT"));
238
		} else if (key == "julian") {
239
240
			julian_present = true;
			indexes.push_back(IOUtils::npos);
241
		} else if (key == "latitude") {
242
			indexes.push_back(IOUtils::npos-1);
243
		} else if (key == "longitude") {
244
			indexes.push_back(IOUtils::npos-2);
245
		} else if (key == "easting") {
246
			indexes.push_back(IOUtils::npos-3);
247
		} else if (key == "northing") {
248
			indexes.push_back(IOUtils::npos-4);
249
		} else if (key == "altitude") {
250
251
			indexes.push_back(IOUtils::npos-5);
		} else {
252
			//this is an extra parameter, we convert to uppercase
253
			std::string extra_param = key;
254
255
256
			IOUtils::toUpper(extra_param);
			md.addParameter(extra_param);
			indexes.push_back(md.getParameterIndex(extra_param));
257
258
259
260
261
262
263
264
265
266
267
		}
	}
}

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
	 */
268
	const double nodata_value = myreader.get_header_doublevalue("nodata");
269
270
271

	meta.position.setProj(coordin, coordinparam); //set the default projection from config file
	if (myreader.location_in_header(smet::WGS84)){
272
273
274
		const double lat = myreader.get_header_doublevalue("latitude");
		const double lon = myreader.get_header_doublevalue("longitude");
		const double alt = myreader.get_header_doublevalue("altitude");
275
		meta.position.setLatLon(lat, lon, alt);
276
	}
277
278

	if (myreader.location_in_header(smet::EPSG)){
279
280
281
282
		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));
283
		meta.position.setEPSG(epsg); //this needs to be set before calling setXY(...)
284
285
		meta.position.setXY(east, north, alt);
	}
286

287
288
289
	meta.stationID = myreader.get_header_value("station_id");
	meta.stationName = myreader.get_header_value("station_name");

290
	const bool data_epsg = myreader.location_in_data(smet::EPSG);
291
	if (data_epsg){
292
293
		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;
294
295
296
297
		meta.position.setEPSG(epsg);
	}
}

298
void SMETIO::copy_data(const smet::SMETReader& myreader,
299
300
                       const std::vector<std::string>& timestamps,
                       const std::vector<double>& mydata, std::vector<MeteoData>& vecMeteo)
301
302
303
304
305
306
{
	/*
	 * 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
	 */
307
	const string myfields = myreader.get_header_value("fields");
308
309
	vector<string> fields;
	IOUtils::readLineToVec(myfields, fields);
310
311
312
313

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

316
	if ((timestamps.empty()) && (!julian_present)) return; //nothing to do
317

318
	const bool pint_present = md.param_exists("PINT");
319
320
	const bool data_wgs84 = myreader.location_in_data(smet::WGS84);
	const bool data_epsg = myreader.location_in_data(smet::EPSG);
321
322
323

	read_meta_data(myreader, md.meta);

324
	const double nodata_value = myreader.get_header_doublevalue("nodata");
325
	double current_timezone = myreader.get_header_doublevalue("tz");
326
327
	if (current_timezone == nodata_value)
		current_timezone = in_dflt_TZ;
328
	const bool timestamp_present = myreader.contains_timestamp();
329

330
331
	const size_t nr_of_fields = indexes.size();
	const size_t nr_of_lines = mydata.size() / nr_of_fields;
332
333
334

	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
335
	double previous_ts = IOUtils::nodata;
336
337
338
	for (size_t ii = 0; ii<nr_of_lines; ii++){
		vecMeteo.push_back(md);
		MeteoData& tmp_md = vecMeteo.back();
339

340
341
342
		if (timestamp_present)
			IOUtils::convertString(tmp_md.date, timestamps[ii], current_timezone);

343
		//Copy data points
344
		for (size_t jj=0; jj<nr_of_fields; jj++){
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
			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)
365
					tmp_md(indexes[jj]) = IOUtils::nodata;
366
				else
367
					tmp_md(indexes[jj]) = current_data;
368
369
370
371
372
373
374
375
376
377
			}

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

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

			current_index++;
		}
378

379
380
381
382
383
384
385
386
		if ((pint_present) && (tmp_md(MeteoData::HNW) == IOUtils::nodata)) {
			const double pint = tmp_md("PINT");
			if (pint==0.) {
				tmp_md(MeteoData::HNW) = 0.;
			} else if (previous_ts!=IOUtils::nodata) {
				const double acc_period = (tmp_md.date.getJulian() - previous_ts) * 24.; //in hours
				tmp_md(MeteoData::HNW) = pint * acc_period;
			}
387
		}
388

389
390
		previous_ts = tmp_md.date.getJulian();
	}
391
392
}

393
void SMETIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
394
                           std::vector< std::vector<MeteoData> >& vecMeteo,
395
                           const size_t& stationindex)
396
{
397
	//Make sure that vecMeteo have the correct dimension and stationindex is valid
398
	size_t startindex=0, endindex=vecFiles.size();
399
	if (stationindex != (size_t)IOUtils::npos){ //HACK do we really still need stationindex??
400
		if ((stationindex < vecFiles.size()) || (stationindex < vecMeteo.size())){
401
402
403
404
			startindex = stationindex;
			endindex = stationindex+1;
		} else {
			throw IndexOutOfBoundsException("Invalid stationindex", AT);
405
		}
406
407
408

		vecMeteo[stationindex].clear();
	} else {
409
		vecMeteo.clear();
410
		vecMeteo = vector< vector<MeteoData> >(vecFiles.size());
411
		vecMeteo.reserve(nr_stations);
412
413
414
	}

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

418
419
420
		if (!IOUtils::fileExists(filename))
			throw FileNotFoundException(filename, AT);

421
422
		smet::SMETReader& myreader = vec_smet_reader.at(ii);
		myreader.convert_to_MKSA(true); // we want converted values for MeteoIO
423

424
425
		vector<double> mydata; //sequentially store all data in the smet file
		vector<string> mytimestamps;
426

427
428
429
		if (myreader.contains_timestamp()){
			myreader.read(dateStart.toString(Date::ISO), dateEnd.toString(Date::ISO), mytimestamps, mydata);
		} else {
430
			myreader.read(dateStart.getJulian(), dateEnd.getJulian(), mydata);
431
		}
432

433
		copy_data(myreader, mytimestamps, mydata, vecMeteo[ii]);
434
	}
435
436
}

437
void SMETIO::writeMeteoData(const std::vector< std::vector<MeteoData> >& vecMeteo, const std::string&)
438
439
{
	//Loop through all stations
440
	for (size_t ii=0; ii<vecMeteo.size(); ii++){
441
		//1. check consistency of station data position -> write location in header or data section
442
		StationData sd;
443
		sd.position.setProj(coordout, coordoutparam);
444
		const bool isConsistent = checkConsistency(vecMeteo.at(ii), sd);
445

446
		if (sd.stationID.empty()){
447
			ostringstream ss;
448
449
450
451
			ss << "Station" << ii+1;
			sd.stationID = ss.str();
		}

452
		const string filename = outpath + "/" + sd.stationID + ".smet";
453
454
455
		if (!IOUtils::validFileName(filename)) //Check whether filename is valid
			throw InvalidFileNameException(filename, AT);

456
		//2. check which meteo parameter fields are actually in use
457
		const size_t nr_of_parameters = getNrOfParameters(sd.stationID, vecMeteo[ii]);
458
459
		vector<bool> vecParamInUse = vector<bool>(nr_of_parameters, false);
		vector<string> vecColumnName = vector<string>(nr_of_parameters, "NULL");
460
		double timezone = IOUtils::nodata; //time zone of the data
461
		checkForUsedParameters(vecMeteo[ii], nr_of_parameters, timezone, vecParamInUse, vecColumnName);
462
		if(out_dflt_TZ != IOUtils::nodata) timezone=out_dflt_TZ; //if the user set an output time zone, all will be converted to it
463

464
		try {
465
			const smet::SMETType type = (outputIsAscii)? smet::ASCII : smet::BINARY;
466

467
			smet::SMETWriter mywriter(filename, type, outputIsGzipped);
468
			generateHeaderInfo(sd, outputIsAscii, isConsistent, timezone,
469
                               nr_of_parameters, vecParamInUse, vecColumnName, mywriter);
470
471
472

			vector<string> vec_timestamp;
			vector<double> vec_data;
473
			for (size_t jj=0; jj<vecMeteo[ii].size(); jj++) {
474
				if (outputIsAscii){
475
					if (out_dflt_TZ != IOUtils::nodata) { //user-specified time zone
476
477
478
479
480
						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));
481
					}
482
483
484
485
486
				} else {
					double julian;
					if(out_dflt_TZ!=IOUtils::nodata) {
						Date tmp_date(vecMeteo[ii][jj].date);
						tmp_date.setTimeZone(out_dflt_TZ);
487
						julian = tmp_date.getJulian();
488
					} else {
489
						julian = vecMeteo[ii][jj].date.getJulian();
490
491
492
					}
					vec_data.push_back(julian);
				}
493

494
				if (!isConsistent) { //Meta data changes
495
496
497
498
					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());
				}
499

500
				for (size_t kk=0; kk<nr_of_parameters; kk++) {
501
					if (vecParamInUse[kk])
502
						vec_data.push_back(vecMeteo[ii][jj](kk)); //add data value
503
				}
504
			}
505

506
507
508
			if (outputIsAscii) mywriter.write(vec_timestamp, vec_data);
			else mywriter.write(vec_data);

509
		} catch(exception&) {
510
			throw;
511
512
513
514
		}
	}
}

515
void SMETIO::generateHeaderInfo(const StationData& sd, const bool& i_outputIsAscii, const bool& isConsistent,
516
517
518
                                const double& timezone, const size_t& nr_of_parameters,
                                const std::vector<bool>& vecParamInUse, const std::vector<std::string>& vecColumnName,
                                smet::SMETWriter& mywriter)
519
520
521
522
523
524
525
526
527
528
{
	/**
	 * 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)
	 */
529
	ostringstream ss;
530
531

	mywriter.set_header_value("station_id", sd.stationID);
532
	if (!sd.stationName.empty())
533
534
535
536
537
		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

538
	if (i_outputIsAscii) {
539
540
541
542
543
544
545
546
547
548
549
550
551
552
		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());
553

554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
		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;
574

575
576
577
578
579
580
581
582
583
584
585
			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);
}

586
void SMETIO::getFormatting(const size_t& param, int& prec, int& width)
587
{
588
589
	/**
	 * When writing a SMET file, different meteo parameters require a different
590
591
	 * format with regard to precision and width when printing.
	 * This procedure sets the precision and width for each known parameter and
592
593
	 * defaults to a width of 8 and precision of 3 digits for each unknown parameter.
	 */
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
	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;
618
619
620
	}
}

621
622
623
size_t SMETIO::getNrOfParameters(const std::string& stationname, const std::vector<MeteoData>& vecMeteo)
{
	/**
624
	 * This function loops through all MeteoData objects present in vecMeteo and returns the
625
626
627
628
629
630
631
632
633
634
635
	 * 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;
636
		const size_t current_size = vecMeteo[ii].getNrOfParameters();
637
638
639
640
641

		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
642
			cerr << "[w] While writing SMET file: Inconsistency in number of meteo "
643
644
645
646
647
648
649
650
651
652
653
654
655
656
				<< "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)
657
{
658
659
660
661
662
663
	/**
	 * 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]
	 */
664
665
666
667
	for (size_t ii=0; ii<vecMeteo.size(); ii++) {
		for (size_t jj=0; jj<nr_parameters; jj++) {
			if (!vecParamInUse[jj]) {
				if (vecMeteo[ii](jj) != IOUtils::nodata) {
668
					vecParamInUse[jj] = true;
669
670
671
					vecColumnName.at(jj) = vecMeteo[ii].getNameForParameter(jj);
				}
			}
672
		}
673
	}
674

675
	if (!vecMeteo.empty())
676
677
678
		timezone = vecMeteo[0].date.getTimeZone();
}

679
bool SMETIO::checkConsistency(const std::vector<MeteoData>& vecMeteo, StationData& sd)
680
{
681
682
683
684
685
686
	/**
	 * 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
	 */

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

690
	for (size_t ii=1; ii<vecMeteo.size(); ii++){
691
692
693
694
695
696
		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;
		}
697
698
699
700
701
	}

	return true;
}

702
void SMETIO::readPOI(std::vector<Coords>& pts)
703
{
704
	std::string filename;
705
	cfg.getValue("POIFILE", "Input", filename);
706
707
708
709
710
711
712
713
714
715
716
	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");

717
	if (myreader.location_in_data(smet::WGS84)==true) {
718
719
		size_t lat_fd=IOUtils::unodata, lon_fd=IOUtils::unodata;
		size_t alt_fd=IOUtils::unodata;
720
		for (size_t ii=0; ii<nr_fields; ii++) {
721
			const string tmp = myreader.get_field_name(ii);
722
723
724
			if (tmp=="latitude") lat_fd=ii;
			if (tmp=="longitude") lon_fd=ii;
			if (tmp=="altitude") alt_fd=ii;
725
		}
726
		for (size_t ii=0; ii<vec_data.size(); ii+=nr_fields) {
727
728
729
730
			Coords point;
			point.setLatLon(vec_data[ii+lat_fd], vec_data[ii+lon_fd], vec_data[ii+alt_fd]);
			pts.push_back(point);
		}
731
732
	} else if (myreader.location_in_data(smet::EPSG)==true) {
		if (epsg==(int)floor(smet_nodata + 0.1))
733
734
735
736
			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;
737
		for (size_t ii=0; ii<nr_fields; ii++) {
738
			const string tmp = myreader.get_field_name(ii);
739
740
741
			if (tmp=="easting") east_fd=ii;
			if (tmp=="northing") north_fd=ii;
			if (tmp=="altitude") alt_fd=ii;
742
		}
743
744
745
746
		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) {
747
748
749
750
751
752
753
754
			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);
	}
755
756
}

757
void SMETIO::write2DGrid(const Grid2DObject& /*grid_in*/, const std::string& /*name*/)
758
759
760
761
762
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

763
764
765
766
767
768
void SMETIO::write2DGrid(const Grid2DObject&, const MeteoGrids::Parameters&, const Date&)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

769
} //namespace