Merging from https://gitlabext.wsl.ch/snow-models/ up to MeteoIO: ff0…

…46fdd and SNOWPACK: 428a31ef.

Source/meteoio/CMakeLists.txt

@@ -90,8 +90,8 @@ SET(PLUGIN_ALPUG ON CACHE BOOL "Compilation ALPUG ON or OFF")
 SET(PLUGIN_ARCIO ON CACHE BOOL "Compilation ARCIO ON or OFF")
 SET(PLUGIN_ARGOSIO OFF CACHE BOOL "Compilation ARGOSIO ON or OFF")
 SET(PLUGIN_ARPSIO ON CACHE BOOL "Compilation ARPSIO ON or OFF")
-SET(PLUGIN_CSVIO ON CACHE BOOL "Compilation CsvIO ON or OFF")
 SET(PLUGIN_COSMOXMLIO OFF CACHE BOOL "Compilation COSMOXMLIO ON or OFF")
+SET(PLUGIN_CSVIO ON CACHE BOOL "Compilation CsvIO ON or OFF")
 SET(PLUGIN_DBO OFF CACHE BOOL "Compilation DBO ON or OFF")
 SET(PLUGIN_GEOTOPIO ON CACHE BOOL "Compilation GEOTOPIO ON or OFF")
 SET(PLUGIN_GOESIO ON CACHE BOOL "Compilation GOESIO ON or OFF")
@@ -105,9 +105,10 @@ SET(PLUGIN_PGMIO ON CACHE BOOL "Compilation PGMIO ON or OFF")
 SET(PLUGIN_PMODIO OFF CACHE BOOL "Compilation PMODIO ON or OFF")
 SET(PLUGIN_PNGIO OFF CACHE BOOL "Compilation PNGIO ON or OFF")
 SET(PLUGIN_PSQLIO OFF CACHE BOOL "Compilation PSQLIO ON or OFF")
+SET(PLUGIN_SASEIO OFF CACHE BOOL "Compilation SASEIO ON or OFF")
 SET(PLUGIN_SMETIO ON CACHE BOOL "Compilation SMETIO ON or OFF")
 SET(PLUGIN_SNIO ON CACHE BOOL "Compilation SNIO ON or OFF")
-SET(PLUGIN_SASEIO OFF CACHE BOOL "Compilation SASEIO ON or OFF")
+SET(PLUGIN_WWCSIO OFF CACHE BOOL "Compilation WWCSIO ON or OFF")
 SET(PLUGIN_ZRXPIO OFF CACHE BOOL "Compilation ZRXPIO ON or OFF")
 SET(PROJ OFF CACHE BOOL "Use PROJ for the class MapProj ON or OFF")
 

Source/meteoio/meteoio/Config.cc

@@ -269,11 +269,12 @@ std::vector<std::string> Config::getKeys(std::string keymatch,
 	} else {
 		keymatch = section + keymatch;
 
-		for (std::map<string,string>::const_iterator it=properties.begin(); it != properties.end(); ++it) {
-			const size_t found_pos = (it->first).find(keymatch, 0);
+		//for (std::map<string,string>::const_iterator it=properties.begin(); it != properties.end(); ++it) {
+		for (const auto& it : properties) {
+			const size_t found_pos = (it.first).find(keymatch, 0);
 			if (found_pos==0) { //found it!
 				const size_t section_len = section.length();
-				const std::string key( (it->first).substr(section_len) ); //from pos to the end
+				const std::string key( (it.first).substr(section_len) ); //from pos to the end
 				vecResult.push_back( key );
 			}
 		}
@@ -291,8 +292,8 @@ void Config::write(const std::string& filename) const
 	try {
 		std::string current_section;
 		unsigned int sectioncount = 0;
-		for (std::map<string,string>::const_iterator it=properties.begin(); it != properties.end(); ++it) {
-			const std::string key_full( it->first );
+		for (const auto& it : properties) {
+			const std::string key_full( it.first );
 			const std::string section( ConfigParser::extract_section(key_full) );
 
 			if (current_section != section) {
@@ -304,7 +305,7 @@ void Config::write(const std::string& filename) const
 			}
 
 			const size_t key_start = key_full.find_first_of(":");
-			const std::string value( it->second );
+			const std::string value( it.second );
 			if (value.empty()) continue;
 
 			if (key_start!=string::npos) //start after the "::" marking the section prefix
@@ -327,7 +328,7 @@ const std::string Config::toString() const {
 	std::ostringstream os;
 	os << "<Config>\n";
 	os << "Source: " << sourcename << "\n";
-	for (map<string,string>::const_iterator it = properties.begin(); it != properties.end(); ++it){
+	for (auto it = properties.begin(); it != properties.end(); ++it){
 		os << (*it).first << " -> " << (*it).second << "\n";
 	}
 	os << "</Config>\n";

Source/meteoio/meteoio/Timer.cc

@@ -23,7 +23,7 @@
 	#include <windows.h>
 	#undef max
 	#undef min
-	#include <cstdlib>
+	#include <threadpoollegacyapiset.h>
 #else
 	#include <sys/time.h>
 	#include <unistd.h>

Source/meteoio/meteoio/dataClasses/DEMAlgorithms.cc

@@ -19,11 +19,16 @@
 #include <cmath>
 #include <limits.h>
 #include <algorithm>
+#include <fstream>
+#include <sstream>
+//#include <cerrno>
+#include <cstring>
 
 #include <meteoio/dataClasses/DEMAlgorithms.h>
 #include <meteoio/dataClasses/DEMObject.h>
 #include <meteoio/MathOptim.h>
 #include <meteoio/IOUtils.h>
+#include <meteoio/FileUtils.h>
 #include <meteoio/meteoLaws/Meteoconst.h> //for math constants
 
 /**
@@ -90,7 +95,7 @@ double DEMAlgorithms::getSearchDistance(const DEMObject& dem)
 * @param[in] ix1 x index of the origin point
 * @param[in] iy1 y index of the origin point
 * @param[in] bearing direction given by a compass bearing
-* @return tangente of angle above the horizontal (in deg)
+* @return tangente of angle above the horizontal (in deg) or IOUtils::nodata if the point (ix1, iy1) does not fit within the provided DEM
 */
 double DEMAlgorithms::getHorizon(const DEMObject& dem, const size_t& ix1, const size_t& iy1, const double& bearing)
 {
@@ -100,6 +105,7 @@ double DEMAlgorithms::getHorizon(const DEMObject& dem, const size_t& ix1, const
 
 	const int dimx = (signed)dem.grid2D.getNx();
 	const int dimy = (signed)dem.grid2D.getNy();
+	if ((signed)ix1>dimx || (signed)iy1>dimy) return IOUtils::nodata; //in case the point does not fith within the provided DEM
 	if (ix1==0 || (signed)ix1==dimx-1 || iy1==0 || (signed)iy1==dimy-1) return 0.; //a border cell is not shadded
 
 	const double cell_alt = dem.grid2D(ix1, iy1);
@@ -135,29 +141,200 @@ double DEMAlgorithms::getHorizon(const DEMObject& dem, const size_t& ix1, const
 * @param[in] dem DEM to work with
 * @param[in] point the origin point
 * @param[in] bearing direction given by a compass bearing
-* @return tangente of angle above the horizontal (in deg)
+* @return tangente of angle above the horizontal (in deg) or IOUtils::nodata if the provided point does not fir within the 
+* provided DEM
 */
-double DEMAlgorithms::getHorizon(const DEMObject& dem, const Coords& point, const double& bearing)
+double DEMAlgorithms::getHorizon(const DEMObject& dem, Coords point, const double& bearing)
 {
-	const int ix1 = (int)point.getGridI();
-	const int iy1 = (int)point.getGridJ();
+	if (!point.indexIsValid()) {
+		if (!dem.gridify(point)) return IOUtils::nodata;
+	}
+
+	const size_t ix1 = static_cast<size_t>( point.getGridI() );
+	const size_t iy1 = static_cast<size_t>( point.getGridJ() );
 	return getHorizon(dem, ix1, iy1, bearing);
 }
 
 /**
-* @brief Returns the horizon from a given point looking 360 degrees around by increments
+* @brief Returns the horizon from a given point looking 360 degrees around by increments. If the provided point does not
+* fit within the provided DEM, an empty result set is returned.
 * @param[in] dem DEM to work with
 * @param[in] point the origin point
 * @param[in] increment to the bearing between two angles
-* @param[out] horizon vector of heights above a given angle
-*
+* @return horizon vector of heights as a function of azimuth
 */
-void DEMAlgorithms::getHorizon(const DEMObject& dem, const Coords& point, const double& increment, std::vector< std::pair<double,double> >& horizon)
+std::vector< std::pair<double,double> > DEMAlgorithms::getHorizonScan(const DEMObject& dem, Coords point, const double& increment)
 {
+	std::vector< std::pair<double,double> > horizon;
+	if (!point.indexIsValid()) {
+		if (!dem.gridify(point)) return horizon;
+	}
+
+	const size_t ix1 = static_cast<size_t>( point.getGridI() );
+	const size_t iy1 = static_cast<size_t>( point.getGridJ() );
 	for (double bearing=0.0; bearing <360.; bearing += increment) {
-		const double tan_alpha = getHorizon(dem, point, bearing);
-		horizon.push_back( make_pair(bearing, atan(tan_alpha)*Cst::to_deg) );
+		const double tan_alpha = getHorizon(dem, ix1, iy1, bearing);
+		if (tan_alpha!=IOUtils::nodata)
+			horizon.push_back( make_pair(bearing, atan(tan_alpha)*Cst::to_deg) );
+	}
+
+	return horizon;
+}
+
+//custom function for sorting the horizons
+struct sort_horizons {
+	bool operator()(const std::pair<double,double> &left, const std::pair<double,double> &right) {
+		if (left.first < right.first) return true; else return false;
+	}
+};
+
+/**
+* @brief Read the horizons from a given set of points looking 360 degrees around provided in a file
+* @details The file containing the horizons is made of any number of lines with the following structure: 
+* `{stationID} {azimuth} {elevation}` where the azimuth is given in degrees North (as read from a compass) and the
+* elevation as degrees above the horizontal. For example:
+* @code
+* STB2 0 2
+* STB2 210 18
+* WFJ 180 5
+* WFJ 130 10
+* STB2 180 25
+* @endcode
+* 
+* @param[in] where Description of the caller to be used in error messages, such as 'Filter::shade'
+* @param[in] filename the file and path containing the horizon
+* @return a map containing for each stationID the horizon vector of heights as a function of azimuth 
+*/
+std::map< std::string, std::vector< std::pair<double,double> > > DEMAlgorithms::readHorizonScan(const std::string& where, const std::string& filename)
+{
+	std::ifstream fin( filename.c_str() );
+	if (fin.fail()) {
+		std::ostringstream ss;
+		ss << where << ": " << "error opening file \"" << filename << "\", possible reason: " << std::strerror(errno);
+		throw AccessException(ss.str(), AT);
+	}
+
+	const char eoln = FileUtils::getEoln(fin); //get the end of line character for the file
+	std::map< std::string, std::vector< std::pair<double,double> > > horizon;
+
+	try {
+		size_t lcount = 0;
+		double azimuth, value;
+		std::string stationID, line;
+		do {
+			lcount++;
+			getline(fin, line, eoln); //read complete line
+			IOUtils::stripComments(line);
+			IOUtils::trim(line);
+			if (line.empty()) continue;
+
+			std::istringstream iss(line);
+			iss.setf(std::ios::fixed);
+			iss.precision(std::numeric_limits<double>::digits10);
+			iss >> std::skipws >> stationID;
+			iss >> std::skipws >> azimuth;
+			if ( !iss || azimuth<0. || azimuth>360.) {
+				std::ostringstream ss;
+				ss << "Invalid azimuth in file " << filename << " at line " << lcount;
+				throw InvalidArgumentException(ss.str(), AT);
+			}
+			iss >> std::skipws >> value;
+			if ( !iss ){
+				std::ostringstream ss;
+				ss << "Invalid value in file " << filename << " at line " << lcount;
+				throw InvalidArgumentException(ss.str(), AT);
+			}
+
+			horizon[stationID].push_back( make_pair(azimuth, value) );
+		} while (!fin.eof());
+		fin.close();
+	} catch (const std::exception&){
+		if (fin.is_open()) {//close fin if open
+			fin.close();
+		}
+		throw;
+	}
+
+	if (horizon.empty()) throw InvalidArgumentException(where+", no valid horizon found in file '"+filename+"'", AT);
+	for (auto& it : horizon) {
+		if (it.second.empty()) throw InvalidArgumentException(where+", no valid horizon for station '"+it.first+"' in file '"+filename+"'", AT);
+		std::sort(it.second.begin(), it.second.end(), sort_horizons());
+	}
+
+	return horizon;
+}
+
+/**
+* @brief Write to a file the horizons from a given set of points looking 360 degrees around
+* @details The file containing the horizons is made of any number of lines with the following structure: 
+* `{stationID} {azimuth} {elevation}` where the azimuth is given in degrees North (as read from a compass) and the
+* elevation as degrees above the horizontal. For example:
+* @code
+* STB2 0 2
+* STB2 210 18
+* WFJ 180 5
+* WFJ 130 10
+* STB2 180 25
+* @endcode
+* 
+* @param[in] horizon a map of vectors of (azimuth, elevation) coordinates defining the horizon, per stationID
+* @param[in] filename the file and path where to write the horizon
+*/
+void DEMAlgorithms::writeHorizons(const std::map< std::string, std::vector< std::pair<double,double> > >& horizon, const std::string& filename)
+{
+	if (!FileUtils::validFileAndPath(filename)) throw InvalidNameException(filename,AT);
+	std::ofstream fout(filename.c_str(), ios::out);
+	if (fout.fail()) {
+		std::ostringstream ss;
+		ss << "error opening file \"" << filename << "\" for writing, possible reason: " << std::strerror(errno);
+		throw AccessException(ss.str(), AT);
+	}
+
+	for (const auto& it : horizon) {
+		const std::string stationID( it.first );
+		for (size_t ii=0; ii<it.second.size(); ii++) 
+			fout << stationID << " " << it.second[ii].first << " " << it.second[ii].second << "\n";
+	}
+
+	if (fout.is_open()) //close fout if open
+		fout.close();
+}
+
+/**
+* @brief Linearly interpolate the horizon height for any given azimuth given a set of (azimuth, elevation) points
+* @param[in] horizon a set of (azimuth, elevation) coordinates defining the horizon
+* @param[in] azimuth the azimuth where to perform the interpolation
+* @return the interpolated elevation for the provided azimuth
+*/
+double DEMAlgorithms::getHorizon(const std::vector< std::pair<double,double> > &horizon, const double& azimuth)
+{
+	if (horizon.empty())
+		throw InvalidArgumentException("attempting to interpolate a horizon elevation from an empty set of horizon points", AT);
+
+	//special (sick) case: only one value given for the horizon -> we consider this is a constant for all azimuths
+	if (horizon.size()==1)
+		return horizon.front().second;
+
+	const std::vector< std::pair<double, double> >::const_iterator next = std::upper_bound(horizon.begin(), horizon.end(), make_pair(azimuth, 0.), sort_horizons()); //first element that is > azimuth
+
+	double x1, y1, x2, y2;
+	if (next!=horizon.begin() && next!=horizon.end()) { //normal case
+		const size_t ii = next - horizon.begin();
+		x1 = horizon[ii-1].first;
+		y1 = horizon[ii-1].second;
+		x2 = horizon[ii].first;
+		y2 = horizon[ii].second;
+	} else {
+		x1 = horizon.back().first - 360.;
+		y1 = horizon.back().second;
+		x2 = horizon.front().first;
+		y2 = horizon.front().second;
 	}
+
+	const double a = (y2 - y1) / (x2 - x1);
+	const double b = y2 - a * x2;
+
+	return a*azimuth + b;
 }
 
 /**

Source/meteoio/meteoio/dataClasses/DEMAlgorithms.h

@@ -33,8 +33,11 @@ class DEMAlgorithms {
 	public:
 		static Grid2DObject getHillshade(const DEMObject& dem, const double& elev, const double& azimuth);
 		static double getHorizon(const DEMObject& dem, const size_t& ix1, const size_t& iy1, const double& bearing);
-		static double getHorizon(const DEMObject& dem, const Coords& point, const double& bearing);
-		static void getHorizon(const DEMObject& dem, const Coords& point, const double& increment, std::vector< std::pair<double,double> >& horizon);
+		static double getHorizon(const DEMObject& dem, Coords point, const double& bearing);
+		static std::vector< std::pair<double,double> > getHorizonScan(const DEMObject& dem, Coords point, const double& increment);
+		static std::map< std::string, std::vector< std::pair<double,double> > > readHorizonScan(const std::string& where, const std::string& filename);
+		static double getHorizon(const std::vector< std::pair<double,double> > &horizon, const double& azimuth);
+		static void writeHorizons(const std::map< std::string, std::vector< std::pair<double,double> > >& horizon, const std::string& filename);
         static double getCellSkyViewFactor(const DEMObject& dem, const size_t& ii, const size_t& jj);
 
 	private:

Source/meteoio/meteoio/dataGenerators/AllSkyLWGenerator.cc

@@ -26,6 +26,13 @@
 
 namespace mio {
 
+AllSkyLWGenerator::AllSkyLWGenerator(const std::vector< std::pair<std::string, std::string> >& vecArgs, const std::string& i_algo, const std::string& i_section, const double& TZ, const Config &i_cfg)
+                   : TauCLDGenerator(vecArgs, i_algo, i_section, TZ, i_cfg), sun(), model(OMSTEDT)
+{ 
+	//TauCLDGenerator will do its own arguments parsing, then AllSkyLWGenerator
+	parse_args(vecArgs); 
+}
+
 void AllSkyLWGenerator::parse_args(const std::vector< std::pair<std::string, std::string> >& vecArgs)
 {
 	const std::string where( section+"::"+algo );
@@ -42,30 +49,18 @@ void AllSkyLWGenerator::parse_args(const std::vector< std::pair<std::string, std
 			else if (user_algo=="UNSWORTH") model = UNSWORTH;
 			else if (user_algo=="CRAWFORD") {
 				model = CRAWFORD;
-				cloudiness_model = TauCLDGenerator::CLF_CRAWFORD;
+				if (cloudiness_model==TauCLDGenerator::DEFAULT)
+					cloudiness_model = TauCLDGenerator::CLF_CRAWFORD;
 			} else if (user_algo=="LHOMME") {
 				model = LHOMME;
-				cloudiness_model = TauCLDGenerator::CLF_LHOMME;
+				if (cloudiness_model==TauCLDGenerator::DEFAULT)
+					cloudiness_model = TauCLDGenerator::CLF_LHOMME;
 			} else
 				throw InvalidArgumentException("Unknown parametrization \""+user_algo+"\" supplied for "+where, AT);
 
 			has_type = true;
-		} else if (vecArgs[ii].first=="CLOUDINESS_TYPE") {
-			user_cloudiness = IOUtils::strToUpper(vecArgs[ii].second);
-		} else if (vecArgs[ii].first=="USE_RSWR") {
-			IOUtils::parseArg(vecArgs[ii], where, use_rswr);
-		} else if (vecArgs[ii].first=="USE_RAD_THRESHOLD") {
-			IOUtils::parseArg(vecArgs[ii], where, use_rad_threshold);
 		}
 	}
-
-	if (!user_cloudiness.empty()) {
-		if (user_cloudiness=="LHOMME") cloudiness_model = TauCLDGenerator::CLF_LHOMME;
-		else if (user_cloudiness=="KASTEN") cloudiness_model = TauCLDGenerator::KASTEN;
-		else if (user_cloudiness=="CRAWFORD") cloudiness_model = TauCLDGenerator::CLF_CRAWFORD;
-		else
-			throw InvalidArgumentException("Unknown parametrization \""+user_cloudiness+"\" supplied for "+where, AT);
-	}
 
 	if (!has_type) throw InvalidArgumentException("Please provide a TYPE for "+where, AT);
 }
@@ -98,7 +93,7 @@ bool AllSkyLWGenerator::generate(const size_t& param, MeteoData& md)
 			sun.setDate(julian_gmt, 0.);
 
 			bool is_night;
-			cloudiness = TauCLDGenerator::getCloudiness(cloudiness_model, md, use_rswr, use_rad_threshold, sun, is_night);
+			cloudiness = TauCLDGenerator::getCloudiness(md, sun, is_night);
 			if (cloudiness==IOUtils::nodata && !is_night) return false;
 
 			if (is_night) { //interpolate the cloudiness over the night