libsmet.cc

/***********************************************************************************/
/*  Copyright 2009 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* 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/>.
*/
#include <meteoio/plugins/libsmet.h>
#include <errno.h>
#include <string.h>
#include <limits>

#if defined _WIN32 || defined __MINGW32__
	#include <windows.h>
#else
	#include <sys/stat.h>
#endif

using namespace std;

namespace smet {

const char* SMETCommon::smet_version = "1.1";
set<string> SMETCommon::all_mandatory_header_keys = set<std::string>();
set<string> SMETCommon::all_optional_header_keys  = set<std::string>();
set<string> SMETCommon::all_decimal_header_values = set<std::string>();

SMETException::SMETException(const std::string& message, const std::string& position)
              : msg( (position.empty())? "At unknown position: " + message : position + ": " + message) {}

SMETException::~SMETException() throw() {}

const char* SMETException::what() const throw()
{
	return msg.c_str();
}

const bool SMETCommon::__init = SMETCommon::initStaticData();
bool SMETCommon::initStaticData()
{
	all_mandatory_header_keys.insert("station_id");
	all_mandatory_header_keys.insert("nodata");
	all_mandatory_header_keys.insert("fields");

	all_optional_header_keys.insert("latitude");
	all_optional_header_keys.insert("longitude");
	all_optional_header_keys.insert("altitude");
	all_optional_header_keys.insert("easting");
	all_optional_header_keys.insert("northing");
	all_optional_header_keys.insert("epsg");

	all_optional_header_keys.insert("station_name");
	all_optional_header_keys.insert("tz");
	all_optional_header_keys.insert("creation");
	all_optional_header_keys.insert("source");
	all_optional_header_keys.insert("units_offset");
	all_optional_header_keys.insert("units_multiplier");
	all_optional_header_keys.insert("comment");

	all_decimal_header_values.insert("nodata");
	all_decimal_header_values.insert("latitude");
	all_decimal_header_values.insert("longitude");
	all_decimal_header_values.insert("altitude");
	all_decimal_header_values.insert("easting");
	all_decimal_header_values.insert("northing");
	all_decimal_header_values.insert("epsg");
	all_decimal_header_values.insert("tz");

	return true;
}

#if defined _WIN32 || defined __MINGW32__
bool SMETCommon::fileExists(const std::string& filename)
{
	const DWORD attributes = GetFileAttributes( filename.c_str() );

	if (attributes==INVALID_FILE_ATTRIBUTES || attributes==FILE_ATTRIBUTE_VIRTUAL
	     || attributes==FILE_ATTRIBUTE_DIRECTORY || attributes==FILE_ATTRIBUTE_DEVICE)
		return false;

	return true;
}
#else
bool SMETCommon::fileExists(const std::string& filename)
{
	struct stat buffer ;

	if ((stat( filename.c_str(), &buffer))!=0) {//File exists if stat returns 0
		return false;
	}

	if (S_ISREG(buffer.st_mode) || S_ISFIFO(buffer.st_mode) || S_ISLNK(buffer.st_mode))
		return true;
	else
		return false; //exclude char device, block device, sockets, etc
}
#endif

bool SMETCommon::validFileAndPath(const std::string& filename)
{
#if defined _WIN32 || defined __MINGW32__ || defined __CYGWIN__
	const size_t startpos = filename.find_first_not_of(" \t\n"); // Find the first character position after excluding leading blank spaces
	const size_t invalid_char = filename.find_first_of("\000*:<>?|"); //find possible invalid characters
#else
	const size_t startpos = filename.find_first_not_of(" \t\n"); // Find the first character position after excluding leading blank spaces
	const size_t invalid_char = filename.find_first_of("\000"); //find possible invalid characters
#endif

	if ((startpos!=0) || (invalid_char!=std::string::npos) || (filename==".") || (filename=="..")) {
		return false;
	}
	return true;
}

void SMETCommon::copy_file(const std::string& src, const std::string& dest)
{
	if (src == dest) return; //copying to the same file doesn't make sense, but is no crime either

	if (!fileExists(src)) throw SMETException("File '"+src+"' does not exist", AT);
	std::ifstream fin(src.c_str(), std::ios::binary);
	if (fin.fail()) throw SMETException("Failed to open file '"+src+"'", AT);

	if (!validFileAndPath(dest)) throw SMETException("Destination file name '"+dest+"' is invalid", AT);
	std::ofstream fout(dest.c_str(), std::ios::binary);
	if (fout.fail()) {
		fin.close();
		throw SMETException("Failed to open destination file '"+dest+"'", AT);
	}

	fout << fin.rdbuf();

	fin.close();
	fout.close();
}

double SMETCommon::convert_to_double(const std::string& in_string)
{
	char* conversion_end = NULL;
	const double conversion_value = strtod(in_string.c_str(), &conversion_end);

	if (*conversion_end == '\0') {
		return conversion_value;
	} else {
		throw SMETException("Value \"" + in_string + "\" cannot be converted to double", SMET_AT);
	}
}

int SMETCommon::convert_to_int(const std::string& in_string)
{
	istringstream ss(in_string);
	int value;
	if (!(ss >> value)) throw SMETException("Value \"" + in_string + "\" cannot be converted to int", SMET_AT);

	return value;
}

void SMETCommon::trim(std::string& str)
{
	const std::string whitespaces(" \t\f\v\n\r");
	const size_t startpos = str.find_first_not_of(whitespaces); // Find the first character position after excluding leading blank spaces
	const size_t endpos = str.find_last_not_of(whitespaces); // Find the first character position from reverse af

	// if all spaces or empty return an empty string
	if (( std::string::npos == startpos ) || ( std::string::npos == endpos)) {
		str.clear();
	} else {
		str.erase(endpos+1); //right trim
		str.erase(0, startpos); //left trim
	}
}

bool SMETCommon::readKeyValuePair(const std::string& in_line, const std::string& delimiter,
                                  std::map<std::string,std::string>& out_map)
{
	//size_t pos = in_line.find(delimiter); //first occurence of '='

	size_t pos = std::string::npos;
	if ((delimiter==" ") || (delimiter=="\t")) {
		pos = in_line.find_first_of(" \t", 0);
	} else {
		pos = in_line.find(delimiter); //first occurence of '='
	}


	if (pos != std::string::npos) { //ignore in_lines that are empty or without '='
		std::string key = in_line.substr(0, pos);
		std::string value = in_line.substr(pos + 1);

		SMETCommon::trim(key);
		SMETCommon::trim(value);

		if (key.empty() || value.empty()) {
			return false;
		}

		out_map[key] = value;
	} else {
		return false;
	}

	return true;
}

void SMETCommon::stripComments(std::string& str)
{
	const size_t found = str.find_first_of("#;");
	if (found != std::string::npos){
		str.erase(found); //rest of line disregarded
	}
}

void SMETCommon::toUpper(std::string& str)
{
	std::transform(str.begin(), str.end(),str.begin(), ::toupper);
}

char SMETCommon::getEoln(std::istream& fin)
{
	std::streambuf* pbuf;
	char tmp = '0';
	int chars = 0;

	const std::streampos position = fin.tellg();

	do {
		fin.get(tmp);
		chars++;

		if ((tmp == '\r') || (tmp == '\n')) {
			char peekc = tmp;
			while ((!fin.eof() && ((peekc=='\r') || (peekc=='\n')))) {
				tmp = peekc;
				fin.get(peekc);
				chars++;
			}
			pbuf = fin.rdbuf();
			pbuf->pubseekpos(position); //rewind
			fin.clear(); //reset eof flag, etc
			return tmp;
		}
	} while ((chars < 3000) && (!fin.eof()));

	pbuf = fin.rdbuf();
	pbuf->pubseekpos(position); //rewind
	fin.clear(); //reset eof flag, etc

	return '\n';
}

bool SMETCommon::is_decimal(const std::string& value)
{
	istringstream ss(value);
	double doublevalue;
	if (!(ss >> doublevalue)) return false;

	return true;
}

size_t SMETCommon::readLineToVec(const std::string& line_in, std::vector<std::string>& vec_string)
{
	vec_string.clear();
	std::istringstream iss(line_in); //construct inputstream with the string line as input

	std::string tmp_string;
	while (!iss.eof()) {
		iss >> std::skipws >> tmp_string;

		if (!tmp_string.empty()) {
			vec_string.push_back( tmp_string );
		}
		tmp_string.clear();
	}

	return vec_string.size();
}


////////////////////////////////////////////////////////////
//// SMETWriter class
SMETWriter::SMETWriter(const std::string& in_filename, const SMETType& in_type)
           : other_header_keys(), ascii_precision(), ascii_width(), header(), mandatory_header_keys(),
             filename(in_filename), nodata_string(), smet_type(in_type), nodata_value(-999.), nr_of_fields(0),
             julian_field(0), timestamp_field(0), location_wgs84(0), location_epsg(0),
             location_in_header(false), location_in_data_wgs84(false), location_in_data_epsg(false),
             timestamp_present(false), julian_present(false), file_is_binary(false), append_mode(false), append_possible(false) {}

//what the caller MUST provide: 
//	* the fields header (for checks)
//	* the ascii width and precision
//	* nodata_value provided through set_header_value() calls
//HACK: the multipliers and offsets are not handled yet...
//for now, we do the most basic append: all the data after the first insertion point is deleted
SMETWriter::SMETWriter(const std::string& in_filename, const std::string& in_fields, const double& in_nodata)
           : other_header_keys(), ascii_precision(), ascii_width(), header(), mandatory_header_keys(),
             filename(in_filename), nodata_string(), smet_type(ASCII), nodata_value(in_nodata), nr_of_fields(0),
             julian_field(0), timestamp_field(0), location_wgs84(0), location_epsg(0),
             location_in_header(false), location_in_data_wgs84(false), location_in_data_epsg(false),
             timestamp_present(false), julian_present(false), file_is_binary(false), append_mode(true), append_possible(false)
{
	SMETReader reader(filename);
	smet_type = (reader.isAscii)? ASCII : BINARY;
	nodata_string = reader.get_header_value("nodata");
	//nodata_value = reader.nodata_value; //we trust the value provided to the constructor
	
	nr_of_fields = reader.nr_of_fields; //in SMETReader, this does NOT contain the timestamp
	julian_field = reader.julian_field;
	timestamp_field = reader.timestamp_field;
	timestamp_present = reader.timestamp_present;
	julian_present = reader.julian_present;
	
	location_wgs84 = reader.location_wgs84;
	location_epsg = reader.location_epsg;
	location_in_header = reader.location_in_header(WGS84) || reader.location_in_header(EPSG);
	location_in_data_wgs84 = reader.location_in_data(WGS84);
	location_in_data_epsg = reader.location_in_data(EPSG);
	
	//check that the fields match
	std::vector<std::string> vecFields;
	SMETCommon::readLineToVec(in_fields, vecFields);
	if (timestamp_present) {
		if (vecFields[timestamp_field]!="timestamp") {
			std::ostringstream ss;
			ss << "Timestamp should be in field at position " << timestamp_field << "when  appending data to file '" << filename << "' but instead '" << reader.get_field_name(timestamp_field) << "' was found";
			throw SMETException(ss.str(), AT);
		}
		vecFields.erase( vecFields.begin()+timestamp_field );
	}
	if (julian_present) {
		if (vecFields[julian_field]!="julian") {
			std::ostringstream ss;
			ss << "Julian should be in field at position " << julian_field << "when  appending data to file '" << filename << "' but instead '" << reader.get_field_name(julian_field) << "' was found";
			throw SMETException(ss.str(), AT);
		}
		vecFields.erase( vecFields.begin()+julian_field );
	}
	
	if (nr_of_fields!=vecFields.size()) {
		std::ostringstream ss;
		ss << "Trying to write " << vecFields.size() << " fields in file '" << filename << "' that has " << nr_of_fields << " fields";
		throw SMETException(ss.str(), AT);
	}
	for (size_t ii=0; ii<nr_of_fields; ii++) {
		if (vecFields[ii] != reader.get_field_name(ii)) {
			std::ostringstream ss;
			ss << "Trying to write field '" << vecFields[ii] << "' at position " << ii << " in file '" << filename << "' when it should be '" << reader.get_field_name(ii);
			throw SMETException(ss.str(), AT);
		}
	}
	
	//adjust the number of fields to reflect the potential presence of timestamp and julian
	if (timestamp_present) nr_of_fields++;
	if (julian_present) nr_of_fields++;
}

void SMETWriter::set_header_value(const std::string& key, const double& value)
{
	//check if key is decimal, transform to string and add to header
	if (SMETCommon::all_decimal_header_values.find(key) != SMETCommon::all_decimal_header_values.end()){
		ostringstream ss;
		if ((key == "latitude") || (key == "longitude") || (key == "easting") || (key == "northing")){
			ss << fixed << setprecision(6) << value;
		} else if (key == "altitude"){
			ss << fixed << setprecision(1) << value;
		} else if ((key == "epsg") || (key == "tz")){
			ss << fixed << setprecision(0) << value;
		} else {
			ss << value; //for nodata
		}

		set_header_value(key, ss.str());
	} else { //It's a non-standard header value
		ostringstream ss;
		ss << value; //for nodata
		set_header_value(key, ss.str());
	}
}

void SMETWriter::set_header_value(const std::string& key, const std::string& value)
{
	//check if header key/value pair is valid
	if (valid_header_pair(key, value)){
		header[key] = value;

		if ((SMETCommon::all_optional_header_keys.find(key) == SMETCommon::all_optional_header_keys.end())
		    && (SMETCommon::all_mandatory_header_keys.find(key) == SMETCommon::all_mandatory_header_keys.end()))
			other_header_keys.push_back(key);
	} else {
		throw SMETException("Invalid, inconsistent or unknown key/value pair: " + key + " = " + value, SMET_AT);
	}
}

bool SMETWriter::valid_header_pair(const std::string& key, const std::string& value)
{
	bool key_ok = true;

	if (SMETCommon::all_mandatory_header_keys.find(key) != SMETCommon::all_mandatory_header_keys.end()){
		mandatory_header_keys.insert(key);
		key_ok = true;
	}

	//nodata value needs extra treatment
	if (key == "nodata"){
		istringstream ss(value);
		if (!(ss >> nodata_value)) return false;
		nodata_string = value;
	}

	//Using WGS_84 coordinate system
	if (key == "latitude")  location_wgs84 |= 1;
	else if (key == "longitude") location_wgs84 |= 2;
	else if (key == "altitude")  location_wgs84 |= 4;

	//Using an EPSG coordinate system
	if (key == "easting")  location_epsg |= 1;
	else if (key == "northing") location_epsg |= 2;
	else if (key == "altitude") location_epsg |= 4;
	else if (key == "epsg")     location_epsg |= 8;

	//Now do some value checks
	if (key == "epsg"){
		istringstream ss(value);
		int intvalue;
		if (!(ss >> intvalue)) return false;
	}

	if (SMETCommon::all_decimal_header_values.find(key) != SMETCommon::all_decimal_header_values.end()){
		//check if value is a double value
		if (!SMETCommon::is_decimal(value)) return false;
	}

	if ((location_epsg == 15) || (location_wgs84 == 7))
		location_in_header = true;

	//Rudimentary check on keys: fields, units_offset, units_multiplier
	if ((key == "fields") || (key == "units_offset") || (key == "units_multiplier"))
		key_ok = check_fields(key, value);

	return key_ok;
}

bool SMETWriter::check_fields(const std::string& key, const std::string& value)
{
	vector<string> tmp_vec;
	const size_t counter = SMETCommon::readLineToVec(value, tmp_vec);

	//Firstly: Check if number of fields is consistent
	if (nr_of_fields != 0){
		if (counter != nr_of_fields) return false;
	} else {
		nr_of_fields = counter;
	}

	if (key == "fields"){
        size_t count_wgs84 = 0, count_epsg = 0;

		//check if location is in data and if timestamp is present
		for (size_t ii = 0; ii<tmp_vec.size(); ii++){
			if (tmp_vec[ii] == "latitude") count_wgs84++;
			else if (tmp_vec[ii] == "longitude") count_wgs84++;
			else if (tmp_vec[ii] == "easting") count_epsg++;
			else if (tmp_vec[ii] == "northing") count_epsg++;

			if (tmp_vec[ii] == "altitude") {
				count_wgs84++;
				count_epsg++;
			}

			if (tmp_vec[ii] == "timestamp"){
				if (timestamp_present) return false; //no duplicate timestamp field allowed
				timestamp_present = true;
				timestamp_field = ii;
			}

			if (tmp_vec[ii] == "julian"){
				if (julian_present) return false; //no duplicate julian field allowed
				julian_present = true;
				julian_field = ii;
			}
		}

		if (count_wgs84 == 3)
			location_in_data_wgs84 = true;

		if (count_epsg == 3)
			location_in_data_epsg = true;
	} else {
		//every value in units_offset and units_multiplier must be a decimal
		for (vector<string>::const_iterator it = tmp_vec.begin(); it != tmp_vec.end(); ++it)
			if (!SMETCommon::is_decimal(*it)) return false;
	}

	return true;
}

bool SMETWriter::valid_header()
{
	if (mandatory_header_keys.size() != 3)
		return false;

	if (location_in_data_epsg){ //EPSG code must be in the header anyway
		const map<string,string>::const_iterator it = header.find("epsg");
		if (it == header.end()){
			return false;
		} else {
			return true;
		}
	}

	if (location_in_header || location_in_data_wgs84)
		return true;

	return false;
}

void SMETWriter::write(const std::vector<std::string>& vec_timestamp, const std::vector<double>& data)
{
	if (vec_timestamp.empty() || data.empty()) return;
	if (!SMETCommon::validFileAndPath(filename)) throw SMETException("Invalid file name \""+filename+"\"", AT);
	errno = 0;
	
	ofstream fout;
	if (append_mode) {
		if (!append_possible) {
			//check where to insert the new data
			SMETReader reader(filename);
			const std::string last_timestamp = reader.getLastTimestamp();
			if (last_timestamp.empty() || last_timestamp>=vec_timestamp[0])
				reader.truncate_file(vec_timestamp[0]);
			
			append_possible = true;
		}
		fout.open(filename.c_str(), ios::binary | ofstream::app);
	} else { //normal mode
		fout.open(filename.c_str(), ios::binary);
	}
	
	if (fout.fail()) {
		ostringstream ss;
		ss << "Error opening file \"" << filename << "\" for writing, possible reason: " << strerror(errno);
		throw SMETException(ss.str(), SMET_AT);
	}

	if (!append_mode) write_header(fout); //Write the header info, always in ASCII format

	if (nr_of_fields == 0) {
		fout.close();
		return;
	}

	if (!timestamp_present) {
		fout.close();
		throw SMETException("No timestamp present when writing file \""+filename+"\", use write(const vector<double>& data)", SMET_AT);
	}

	const size_t nr_of_data_fields = nr_of_fields - 1;
	if (nr_of_fields<=1) { //avoid division by zero in the next block
		cerr << "[W] Attempting to write a dataset that contains no fields to file \"" + filename + "\"!\n";
		fout.close();
		return;
	}

	const size_t nr_of_lines = data.size() / (nr_of_fields-1);
	if ((nr_of_lines != vec_timestamp.size()) || ((data.size() % (nr_of_fields-1)) != 0)) {
		fout.close();
		ostringstream os;
		os << "Inconsistency between the number of timestamp (" << vec_timestamp.size() << ")  and number of data points (" << data.size() << "/" << nr_of_fields-1 << ") detected for file \""+filename+"\", recheck your data";
		throw SMETException(os.str(), SMET_AT);
	}

	std::vector<double> current_data(nr_of_fields-1);
	check_formatting();

	if (smet_type == ASCII) {
		for (size_t ii=0; ii<nr_of_lines; ii++) {
			const size_t offset = ii*(nr_of_fields-1);
			if (!data.empty())
				copy(data.begin()+offset, data.begin()+offset+nr_of_data_fields, current_data.begin());
			write_data_line_ascii(vec_timestamp[ii], current_data, fout);
		}
	} else {
		fout.close();
		throw SMETException("Cannot write binary file \""+filename+"\" with a timestamp, use julian instead", SMET_AT);
	}

	fout.close();
}

void SMETWriter::write(const std::vector<double>& data)
{
	if (!SMETCommon::validFileAndPath(filename)) throw SMETException("Invalid file name \""+filename+"\"", AT);
	errno = 0;
	std::ofstream fout;
	fout.open(filename.c_str(), ios::binary);
	if (fout.fail()) {
		ostringstream ss;
		ss << "Error opening file \"" << filename << "\" for writing, possible reason: " << strerror(errno);
		throw SMETException(ss.str(), SMET_AT);
	}

	write_header(fout); //Write the header info, always in ASCII format

	if (nr_of_fields == 0){
		fout.close();
		return;
	}

	const size_t nr_of_lines = data.size() / nr_of_fields;
	if ((data.size() % nr_of_fields) != 0) {
		fout.close();
		throw SMETException("Inconsistency between data and header fields detected in file \""+filename+"\", recheck your data", SMET_AT);
	}

	std::vector<double> current_data(nr_of_fields);
	check_formatting();

	if (smet_type == ASCII){
		for (size_t ii=0; ii<nr_of_lines; ii++){
			if (!data.empty())
				copy(data.begin()+ii*nr_of_fields, data.begin()+ii*nr_of_fields+nr_of_fields, current_data.begin());
			write_data_line_ascii("0000-01-01T00:00", current_data, fout); //dummy time
		}
	} else {
		for (size_t ii=0; ii<nr_of_lines; ii++){
			if (!data.empty())
				copy(data.begin()+ii*nr_of_fields, data.begin()+ii*nr_of_fields+nr_of_fields, current_data.begin());
			write_data_line_binary(current_data, fout);
		}

		file_is_binary = false;
	}

	fout.close();
}

void SMETWriter::write_header(std::ofstream& fout)
{
	if (!valid_header()) {
		fout.close();
		throw SMETException("The header data you supplied is not valid, file \""+filename+"\" cannot be written", SMET_AT);
	}

	//write signature
	fout << "SMET " << SMETCommon::smet_version << " ";
	if (smet_type == ASCII) fout << "ASCII" << "\n";
	else fout << "BINARY" << "\n";

	fout << "[HEADER]" << "\n";
	fout << "station_id       = " << header["station_id"] << "\n";

	map<string,string>::const_iterator it = header.find("station_name");
	if (it != header.end())
		fout << "station_name     = " << it->second << "\n";

	if (location_in_header){
		if (location_wgs84 == 7){
			fout << "latitude         = " << header["latitude"]  << "\n";
			fout << "longitude        = " << header["longitude"] << "\n";
			fout << "altitude         = " << header["altitude"]  << "\n";
		}

		if (location_epsg == 15){
			fout << "easting          = " << header["easting"]   << "\n";
			fout << "northing         = " << header["northing"]  << "\n";
			if (location_wgs84 != 7)
				fout << "altitude        = " << header["altitude"]  << "\n";
			fout << "epsg             = " << header["epsg"]  << "\n";
		}
	} else {
		if (location_in_data_epsg)
			fout << "epsg             = " << header["epsg"]      << "\n";
	}

	fout << "nodata           = " << header["nodata"] << "\n";

	//Optional header keys:
	it = header.find("tz");
	if (it != header.end())
		fout << "tz               = " << it->second << "\n";

	it = header.find("creation");
	if (it != header.end())
		fout << "creation         = " << it->second << "\n";

	it = header.find("source");
	if (it != header.end())
		fout << "source           = " << it->second << "\n";

	it = header.find("units_offset");
	if (it != header.end())
		fout << "units_offset     = " << it->second << "\n";

	it = header.find("units_multiplier");
	if (it != header.end())
		fout << "units_multiplier = " << it->second << "\n";

	it = header.find("comment");
	if (it != header.end())
		fout << "comment          = " << it->second << "\n";

	for (size_t ii=0; ii<other_header_keys.size(); ii++){
		fout << other_header_keys[ii] << " = " << header[other_header_keys[ii]] << "\n";
	}

	fout << "fields           = " << header["fields"] << "\n";
	fout << "[DATA]" << endl;
}

void SMETWriter::write_data_line_binary(const std::vector<double>& data, std::ofstream& fout)
{
	const char eoln = '\n';

	if (!file_is_binary){ //open it as binary file
		fout.close();
		errno = 0;
		fout.open(filename.c_str(), ios::out | ios::app | ios::binary); //reopen as binary file
		if (fout.fail()) {
			ostringstream ss;
			ss << "Error writing to file \"" << filename << "\", possible reason: " << strerror(errno);
			throw SMETException(ss.str(), SMET_AT);
		}

		file_is_binary = true;
	}

	for (size_t ii = 0; ii < data.size(); ii++){
		if (julian_present && (julian_field == ii)){
			const double julian = data[ii];
			fout.write((char*)&julian, sizeof(double)); //the julian date is written in 64bit IEEE754 precision
		} else {
			const float val = (float)data[ii];
			fout.write((char*)&val, sizeof(float)); //normal data fields are written in 32bit IEEE754 precision
		}
	}

	fout.write((char*)&eoln, sizeof(char));
}

void SMETWriter::write_data_line_ascii(const std::string& timestamp, const std::vector<double>& data, std::ofstream& fout)
{
	fout.fill(' ');
	fout << right;
	fout << fixed;

	if ((data.empty()) && timestamp_present) fout << timestamp;

	for (size_t ii = 0; ii < data.size(); ii++){
		if (ii > 0) fout << " ";
		if (timestamp_present && (timestamp_field == ii))	fout << timestamp << " ";

		fout << setw(ascii_width[ii]) << setprecision(ascii_precision[ii]);
		if (data[ii] == nodata_value) fout << nodata_string; //to have a nicer representation
		else fout << data[ii];
	}
	fout << "\n";
}

void SMETWriter::check_formatting()
{
	const size_t nr_of_data_fields = (timestamp_present)? nr_of_fields-1 : nr_of_fields;

	if ((ascii_precision.size() != nr_of_data_fields) || (ascii_width.size() != nr_of_data_fields)){
		ascii_precision.resize(nr_of_data_fields, 3);
		ascii_width.resize(nr_of_data_fields, 8);
	}

	if (julian_present){
		ascii_precision[julian_field] = 8;
		ascii_width[julian_field] = 16;
	}
}

void SMETWriter::set_width(const std::vector<int>& vec_width)
{
	ascii_width = vec_width;
}

void SMETWriter::set_precision(const std::vector<int>& vec_precision)
{
	ascii_precision = vec_precision;
}


////////////////////////////////////////////////////////////
//// SMETReader class
const size_t SMETReader::streampos_every_n_lines = 2000; //save streampos every 2000 lines of data

SMETReader::SMETReader(const std::string& in_fname)
            : data_start_fpointer(), vec_offset(), vec_multiplier(), vec_fieldnames(),
              header(), indexer(),
              filename(in_fname), timestamp_start("-4714-11-24T00:00"),
              timestamp_end("9999-12-31T00:00"), nodata_value(-999.),
              julian_start(0.), julian_end(5373483.5),
              nr_of_fields(0), timestamp_field(0), julian_field(0),
              location_wgs84(0), location_epsg(0), location_data_wgs84(0), location_data_epsg(0),
              eoln('\n'),
              timestamp_present(false), julian_present(false), isAscii(true), mksa(true),
              timestamp_interval(false), julian_interval(false)
{
	std::ifstream fin; //Input file streams
	if (!SMETCommon::fileExists(filename)) throw SMETException("File '"+filename+"' does not exists", AT); //prevent invalid filenames
	errno = 0;
	fin.open (filename.c_str(), ios::in|ios::binary); //ascii does end of line translation, which messes up the pointer code
	if (fin.fail()) {
		ostringstream ss;
		ss << "Error opening file \"" << filename << "\" for reading, possible reason: " << strerror(errno);
		ss << " Please check file existence and permissions!";
		throw SMETException(ss.str(), SMET_AT);
	}

	try {
		eoln = SMETCommon::getEoln(fin); //get the end of line character for the file
		read_header(fin);
		process_header();
	} catch(...){
		cleanup(fin); //closes file
		throw;
	}
	cleanup(fin); //closes file
}

SMETReader::~SMETReader() {}

void SMETReader::cleanup(std::ifstream& fin) throw()
{
	if (fin.is_open()) //close fin if open
		fin.close();
}

void SMETReader::convert_to_MKSA(const bool& in_mksa)
{
	mksa = in_mksa;
}

std::string SMETReader::get_field_name(const size_t& nr_of_field)
{
	if (nr_of_field < nr_of_fields){
		return vec_fieldnames[nr_of_field];
	} else {
		ostringstream ss;
		ss << "Trying to access field #" << nr_of_field << " (starting from 0) of " << nr_of_fields << " fields in file \"" << filename << "\". ";
		ss << "This is out of bounds!";
		throw SMETException(ss.str(), SMET_AT);
	}
}

size_t SMETReader::get_nr_of_fields() const
{
	return nr_of_fields;
}

void SMETReader::get_units_conversion(std::vector<double>& offset, std::vector<double>& multiplier) const
{
	multiplier = vec_multiplier;
	offset = vec_offset;
}

bool SMETReader::location_in_header(const LocationType& type) const
{
	if ((location_epsg == 15) && (type == EPSG)){
		return true;
	} else if ((location_wgs84 == 7) && (type == WGS84)){
		return true;
	}

	return false;
}

bool SMETReader::location_in_data(const LocationType& type) const
{
	if ((location_epsg == 8) && (location_data_epsg == 7) && (type == EPSG)){
		return true;
	} else if ((location_data_wgs84 == 7) && (type == WGS84)){
		return true;
	}

	return false;
}

void SMETReader::process_header()
{
	vector<string> tmp_vec;
	set<string> obligatory_keys;
	for (map<string,string>::iterator it = header.begin(); it != header.end(); ++it){
		if (SMETCommon::all_mandatory_header_keys.find(it->first) != SMETCommon::all_mandatory_header_keys.end())
			obligatory_keys.insert(it->first);

		if (it->first == "fields"){
			SMETCommon::readLineToVec(it->second, tmp_vec);
			string newfields;
			if (!tmp_vec.empty()){
				for (size_t ii=0; ii<tmp_vec.size(); ii++){
					if (tmp_vec[ii] == "timestamp"){
						timestamp_present = true;
						timestamp_field = ii;
					} else {
						if (nr_of_fields > 0) newfields += " ";
						vec_fieldnames.push_back(tmp_vec[ii]);
						newfields += tmp_vec[ii];
						nr_of_fields++;
					}

					if (tmp_vec[ii] == "julian"){
						julian_present = true;
						julian_field = ii;
					}

					//Using WGS_84 coordinate system
					if (tmp_vec[ii] == "latitude")  location_data_wgs84 |= 1;
					else if (tmp_vec[ii] == "longitude") location_data_wgs84 |= 2;
					else if (tmp_vec[ii] == "altitude")  location_data_wgs84 |= 4;

					//Using an EPSG coordinate system
					if (tmp_vec[ii] == "easting")  location_data_epsg |= 1;
					else if (tmp_vec[ii] == "northing") location_data_epsg |= 2;
					else if (tmp_vec[ii] == "altitude") location_data_epsg |= 4;
				}
				it->second = newfields;
			}
		}

		if (it->first == "nodata")
			nodata_value = SMETCommon::convert_to_double(it->second);

		//Using WGS_84 coordinate system
		if (it->first == "latitude")  location_wgs84 |= 1;
		else if (it->first == "longitude") location_wgs84 |= 2;
		else if (it->first == "altitude")  location_wgs84 |= 4;

		//Using an EPSG coordinate system
		if (it->first == "easting")  location_epsg |= 1;
		else if (it->first == "northing") location_epsg |= 2;
		else if (it->first == "altitude") location_epsg |= 4;
		else if (it->first == "epsg")     location_epsg |= 8;

		//Now do a value check on EPSG
		if (it->first == "epsg")
			SMETCommon::convert_to_int(it->second);
	}

	if (get_header_value("units_offset") != ""){
		vector<string> tmp_offset;
		string offsetstring;
		SMETCommon::readLineToVec(get_header_value("units_offset"), tmp_offset);

		for (size_t ii=0; ii<tmp_offset.size(); ii++){
			if (!timestamp_present || (ii != timestamp_field)){
				if (ii>0) offsetstring += " ";
				offsetstring += tmp_offset[ii];
				vec_offset.push_back(SMETCommon::convert_to_double(tmp_offset[ii]));
			}
		}
		header["units_offset"] = offsetstring;
	} else {
		vec_offset.resize(nr_of_fields, 0.0);
	}

	if (get_header_value("units_multiplier") != ""){
		vector<string> tmp_multiplier;
		string multiplierstring;
		SMETCommon::readLineToVec(get_header_value("units_multiplier"), tmp_multiplier);

		for (size_t ii=0; ii<tmp_multiplier.size(); ii++){
			if (!timestamp_present || (ii != timestamp_field)){
				if (ii>0) multiplierstring += " ";
				multiplierstring += tmp_multiplier[ii];
				vec_multiplier.push_back(SMETCommon::convert_to_double(tmp_multiplier[ii]));
			}
		}
		header["units_multiplier"] = multiplierstring;
	} else {
		vec_multiplier.resize(nr_of_fields, 1.0);
	}

	if (obligatory_keys.size() != SMETCommon::all_mandatory_header_keys.size())
		throw SMETException("\""+filename+"\" is not a valid SMET file, mandatory information in header is missing", SMET_AT);

	if ((location_wgs84 == 7) || (location_epsg == 15) || (location_data_wgs84 == 7)
	    || ((location_epsg == 8) && (location_data_epsg == 7))){
		//location info present somewhere
	} else {
		throw SMETException("\""+filename+"\" is not a valid SMET file, mandatory location info is missing (header and data)", SMET_AT);
	}

	const size_t nr_offset = vec_offset.size();
	const size_t nr_multiplier = vec_multiplier.size();
	if (nr_offset!=nr_of_fields || nr_multiplier!=nr_of_fields) {
		const size_t Nfields = (timestamp_present)? nr_of_fields+1 : nr_of_fields;
		const size_t Noffset = (timestamp_present)? nr_offset+1 : nr_offset;
		const size_t Nmultiplier =(timestamp_present)? nr_multiplier+1 : nr_multiplier;
		std::ostringstream ss;