PINT is now properly supported in SMET. The OLWR/TSS conversion must now be...

PINT is now properly supported in SMET. The OLWR/TSS conversion must now be handled by a data generator (implemented as "TS_OLWR"). A bug in the Coords documentation has been fixed.

meteoio/GeneratorAlgorithms.cc

@@ -38,6 +38,8 @@ GeneratorAlgorithm* GeneratorAlgorithmFactory::getAlgorithm(const std::string& i
 		return new StandardPressureGenerator(vecArgs, i_algoname);
 	} else if (algoname == "RELHUM"){
 		return new RhGenerator(vecArgs, i_algoname);
+	} else if (algoname == "TS_OLWR"){
+		return new TsGenerator(vecArgs, i_algoname);
 	} else if (algoname == "CLEARSKY_LW"){
 		return new ClearSkyLWGenerator(vecArgs, i_algoname);
 	} else if (algoname == "ALLSKY_LW"){
@@ -246,6 +248,42 @@ bool RhGenerator::generate(const size_t& param, std::vector<MeteoData>& vecMeteo
 }
 
 
+const double TsGenerator::e_snow = .983; //snow emissivity (0.969 - 0.997)
+const double TsGenerator::e_soil = .9805; //grass emissivity (0.975 - 0.986)
+const double TsGenerator::snow_thresh = .1; //if snow height greater than this threshold -> snow albedo
+
+bool TsGenerator::generate(const size_t& param, MeteoData& md)
+{
+	double &value = md(param);
+	if(value == IOUtils::nodata) {
+		const double olwr = md("OLWR");
+		if (olwr==IOUtils::nodata) //nothing else we can do here
+			return false;
+
+		const double hs = md(MeteoData::HS);
+		const double ea = (hs==IOUtils::nodata)? .5*(e_snow+e_soil) : (hs>snow_thresh)? e_snow : e_soil;
+
+		//value = pow( olwr / ( ea * Cst::stefan_boltzmann ), 0.25);
+		value = Optim::invSqrt( Optim::invSqrt(olwr / ( ea * Cst::stefan_boltzmann )) );
+
+		if (value==IOUtils::nodata) return false;
+	}
+
+	return true; //all missing values could be filled
+}
+
+bool TsGenerator::generate(const size_t& param, std::vector<MeteoData>& vecMeteo)
+{
+	if(vecMeteo.empty()) return true;
+
+	for(size_t ii=0; ii<vecMeteo.size(); ii++) {
+		generate(param, vecMeteo[ii]);
+	}
+
+	return true; //all missing values could be filled
+}
+
+
 void ClearSkyLWGenerator::parse_args(const std::vector<std::string>& vecArgs)
 {
 	//Get the optional arguments for the algorithm: constant value to use

meteoio/GeneratorAlgorithms.h

@@ -71,6 +71,7 @@ namespace mio {
  * The keywords defining the algorithms are the following:
  * - STD_PRESS: standard atmospheric pressure as a function of the elevation of each station (see StandardPressureGenerator)
  * - RELHUM: relative humidity from other humidity measurements (see RhGenerator)
+ * - TS_OLWR: surface temperature from Outgoing Long Wave Radiation (see TssGenerator)
  * - CST: constant value as provided in argument (see ConstGenerator)
  * - SIN: sinusoidal variation (see SinGenerator)
  * - CLEARSKY_LW: use a clear sky model to generate ILWR from TA, RH (see ClearSkyLWGenerator)
@@ -209,6 +210,23 @@ class RhGenerator : public GeneratorAlgorithm {
 		bool generate(const size_t& param, std::vector<MeteoData>& vecMeteo);
 };
 
+/**
+ * @class TsGenerator
+ * @brief Surface temperature generator.
+ * Generate the surface temperature from the outgoing long wave (OLWR).
+ * @code
+ * TSS::generators = TS_OLWR
+ * @endcode
+ */
+class TsGenerator : public GeneratorAlgorithm {
+	public:
+		TsGenerator(const std::vector<std::string>& vecArgs, const std::string& i_algo)
+			: GeneratorAlgorithm(vecArgs, i_algo) { parse_args(vecArgs); }
+		bool generate(const size_t& param, MeteoData& md);
+		bool generate(const size_t& param, std::vector<MeteoData>& vecMeteo);
+	private:
+		static const double e_snow, e_soil, snow_thresh;
+};
 
 /**
  * @class ClearSkyLWGenerator

meteoio/dataClasses/Coords.cc

@@ -1607,8 +1607,8 @@ void Coords::WGS84_to_local(double lat_in, double long_in, double& east_out, dou
 
 /**
 * @brief Coordinate conversion: from local grid as given in coordparam to WGS84 Lat/Long
-* @param east_in easting coordinate (Swiss system)
-* @param north_in northing coordinate (Swiss system)
+* @param east_in easting coordinate (in local system)
+* @param north_in northing coordinate (in local system)
 * @param lat_out Decimal Latitude
 * @param long_out Decimal Longitude
 */

meteoio/plugins/SMETIO.cc

@@ -232,7 +232,7 @@ void SMETIO::identify_fields(const std::vector<std::string>& fields, std::vector
 		} else if (key == "OLWR") {
 			md.addParameter("OLWR");
 			indexes.push_back(md.getParameterIndex("OLWR"));
-		} else if (key == "PINT") {
+		} else if (key == "PINT") { //in mm/h
 			md.addParameter("PINT");
 			indexes.push_back(md.getParameterIndex("PINT"));
 		} else if (key == "julian") {
@@ -315,7 +315,7 @@ void SMETIO::copy_data(const smet::SMETReader& myreader,
 
 	if ((timestamps.empty()) && (!julian_present)) return; //nothing to do
 
-	const bool olwr_present = md.param_exists("OLWR");
+	const bool pint_present = md.param_exists("PINT");
 	const bool data_wgs84 = myreader.location_in_data(smet::WGS84);
 	const bool data_epsg = myreader.location_in_data(smet::EPSG);
 
@@ -332,6 +332,7 @@ void SMETIO::copy_data(const smet::SMETReader& myreader,
 
 	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
+	double previous_ts = IOUtils::nodata;
 	for (size_t ii = 0; ii<nr_of_lines; ii++){
 		vecMeteo.push_back(md);
 		MeteoData& tmp_md = vecMeteo.back();
@@ -375,18 +376,18 @@ void SMETIO::copy_data(const smet::SMETReader& myreader,
 			current_index++;
 		}
 
-		if ((olwr_present) && (tmp_md(MeteoData::TSS) == IOUtils::nodata)) {//HACK
-			tmp_md(MeteoData::TSS) = olwr_to_tss(tmp_md("OLWR"));
+		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;
+			}
 		}
-	}
-}
 
-double SMETIO::olwr_to_tss(const double& olwr) {
-	const double ea = 1.;
-	if(olwr==IOUtils::nodata) return IOUtils::nodata;
-	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
-	return pow( olwr / ( ea * Cst::stefan_boltzmann ), 0.25);
+		previous_ts = tmp_md.date.getJulian();
+	}
 }
 
 void SMETIO::readMeteoData(const Date& dateStart, const Date& dateEnd,