WSL/SLF GitLab Repository

Commit ce2301d8 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

The virtual stations have been re-enabled. A new generator has been...

The virtual stations have been re-enabled. A new generator has been implemented to generate ISWR from RSWR or the opposite. Some documentation has been improved. The Black Globe Temperature index should now work.
parent 3b884d89
......@@ -42,6 +42,8 @@ GeneratorAlgorithm* GeneratorAlgorithmFactory::getAlgorithm(const std::string& i
return new RhGenerator(vecArgs, i_algoname);
} else if (algoname == "TS_OLWR"){
return new TsGenerator(vecArgs, i_algoname);
} else if (algoname == "ISWR_ALBEDO"){
return new IswrAlbedoGenerator(vecArgs, i_algoname);
} else if (algoname == "CLEARSKY_LW"){
return new ClearSkyLWGenerator(vecArgs, i_algoname);
} else if (algoname == "ALLSKY_LW"){
......@@ -278,13 +280,59 @@ bool TsGenerator::generate(const size_t& param, std::vector<MeteoData>& vecMeteo
{
if(vecMeteo.empty()) return true;
bool status = true;
for(size_t ii=0; ii<vecMeteo.size(); ii++) {
generate(param, vecMeteo[ii]);
if (!generate(param, vecMeteo[ii]))
status = false;
}
return status;
}
const double IswrAlbedoGenerator::soil_albedo = .23; //grass
const double IswrAlbedoGenerator::snow_albedo = .85; //snow
const double IswrAlbedoGenerator::snow_thresh = .1; //if snow height greater than this threshold -> snow albedo
bool IswrAlbedoGenerator::generate(const size_t& param, MeteoData& md)
{
double &value = md(param);
if(value == IOUtils::nodata) {
const double HS=md(MeteoData::HS), RSWR=md(MeteoData::RSWR), ISWR=md(MeteoData::ISWR);
double albedo = .5;
if (HS!=IOUtils::nodata)
albedo = (HS>=snow_thresh)? snow_albedo : soil_albedo;
if (param==MeteoData::ISWR && (RSWR!=IOUtils::nodata && HS!=IOUtils::nodata)) {
value = RSWR / albedo;
return true;
}
if (param==MeteoData::RSWR && (ISWR!=IOUtils::nodata && HS!=IOUtils::nodata)) {
value = ISWR * albedo;
return true;
}
return false;
}
return true; //all missing values could be filled
}
bool IswrAlbedoGenerator::generate(const size_t& param, std::vector<MeteoData>& vecMeteo)
{
if(vecMeteo.empty()) return true;
bool status = true;
for(size_t ii=0; ii<vecMeteo.size(); ii++) {
if (!generate(param, vecMeteo[ii]))
status = false;
}
return status;
}
void ClearSkyLWGenerator::parse_args(const std::vector<std::string>& vecArgs)
{
......
......@@ -72,6 +72,7 @@ namespace mio {
* - 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)
* - ISWR_ALBEDO: ISWR from RSWR or RSWR from ISWR with either a snow or a soil albedo, depending on HS (see IswrAlbedoGenerator)
* - 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)
......@@ -228,6 +229,25 @@ class TsGenerator : public GeneratorAlgorithm {
static const double e_snow, e_soil, snow_thresh;
};
/**
* @class IswrAlbedoGenerator
* @brief Incoming or reflected short wave generator.
* Generate the incoming short wave radiation from the reflected short wave radiation or the opposite. The albedo
* ie either a grassy soil albedo or a snow albedo depending on the snow height.
* @code
* ISWR::generators = ISWR_ALBEDO
* @endcode
*/
class IswrAlbedoGenerator : public GeneratorAlgorithm {
public:
IswrAlbedoGenerator(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 soil_albedo, snow_albedo, snow_thresh;
};
/**
* @class ClearSkyLWGenerator
* @brief ILWR clear sky parametrization
......
......@@ -121,7 +121,7 @@ class Meteo2DInterpolator; // forward declaration, cyclic header include
* @section interpol2D_biblio Bibliography
* The interpolation algorithms have been inspired by the following papers:
* - <i>"A Meteorological Distribution System for High-Resolution Terrestrial Modeling (MicroMet)"</i>, Liston and Elder, Journal of Hydrometeorology <b>7</b> (2006), 217-234.
* - <i>"Simulating wind fields and snow redistribution using terrain-based parameters to model snow accumulation and melt over a semi-arid mountain catchment"</i>, Adam Winstral and Danny Marks, Hydrological Processes <b>16</b> (2002), 3585– 3603. DOI: 10.1002/hyp.1238 [NOT YET IMPLEMENTED]
* - <i>"Simulating wind fields and snow redistribution using terrain-based parameters to model snow accumulation and melt over a semi-arid mountain catchment"</i>, Adam Winstral and Danny Marks, Hydrological Processes <b>16</b> (2002), 3585– 3603. DOI: 10.1002/hyp.1238
* - <i>"Quantitative evaluation of different hydrological modelling approaches in a partly glacierized Swiss watershed"</i>, Jan Magnusson, Daniel Farinotti, Tobias Jonas and Mathias Bavay, Hydrological Processes, 2010, under review.
* - <i>"Modelling runoff from highly glacierized alpine catchments in a changing climate"</i>, Matthias Huss, Daniel Farinotti, Andreas Bauder and Martin Funk, Hydrological Processes, <b>22</b>, 3888-3902, 2008.
* - <i>"Geostatistics for Natural Resources Evaluation"</i>, Pierre Goovaerts, Oxford University Press, Applied Geostatistics Series, 1997, 483 p., ISBN 0-19-511538-4
......
......@@ -346,6 +346,7 @@ void Meteo2DInterpolator::initVirtualStations()
std::vector<std::string> vecStation;
cfg.getValues("Vstation", "INPUT", vecStation);
for(size_t ii=0; ii<vecStation.size(); ii++) {
//The coordinate specification is given as either: "easting northing epsg" or "lat lon"
Coords tmp(coordin, coordinparam, vecStation[ii]);
if(!tmp.isNodata())
v_coords.push_back( tmp );
......@@ -388,7 +389,7 @@ size_t Meteo2DInterpolator::getVirtualStationsData(const Date& i_date, METEO_SET
//get data from real input stations
METEO_SET vecTrueMeteo;
//getStationsMeteoData(i_date, vecTrueMeteo); //HACK
tsmanager.getMeteoData(i_date, vecTrueMeteo);
if (vecTrueMeteo.empty()) return 0;
if (v_params.empty()) {
......
......@@ -23,26 +23,87 @@
namespace mio {
//we start easy: this is NOT a true ring buffer, just the same as what we previously had (for now)
/**
* @class MeteoBuffer
* @brief A class to buffer meteorological data.
* This class buffers MeteoData objects. It is currently NOT a proper ring buffer, this should come
* in a later implementation.
*
* @ingroup data_str
* @author Mathias Bavay
* @date 2014-12-01
*/
class MeteoBuffer {
public:
MeteoBuffer() : ts_buffer(), ts_start(), ts_end() {};
/**
* @brief Get buffer data for a specific date
* @param date A Date object representing the date/time for the sought MeteoData objects
* @param vecMeteo A vector of MeteoData objects to be filled with data
* @return true if the data was in the buffer
*/
bool get(const Date& date, METEO_SET &vecMeteo) const;
/**
* @brief Get buffer data between two dates
* @param date_start A Date object representing the beginning of an interval (inclusive)
* @param date_end A Date object representing the end of an interval (inclusive)
* @param vecMeteo A vector of vector<MeteoData> objects to be filled with data
* @return true if the data was in the buffer
*/
bool get(const Date& date_start, const Date& date_end, std::vector< METEO_SET > &vecMeteo) const;
/**
* @brief Returns the average sampling rate in the data.
* This computes the average sampling rate of the data that is contained in the buffer. This is a quick
* estimate, centered on how often a station measures "something" (ie, how many timestamps do we have
* for this station in the buffer). if the station measures TA at h+0 and h+30 and
* RH at h+15 and h+45, it would return 4 measurements per hour. If the station measures TA and RH at h+0 and h+30,
* it would return 2 measurements per hour.
* @return average sampling rate in Hz, nodata if the buffer is empty
*/
double getAvgSamplingRate() const;
/**
* @brief Returns the begining of the buffer.
* This is the start date of the <b>request</b> that was given to the IOHandler. If there was no data
* at this date, then the date of the first data would be greater.
* @return start date of the buffer
*/
Date getBufferStart() const;
/**
* @brief Returns the end of the buffer.
* This is the end date of the <b>request</b> that was given to the IOHandler. If there was no data
* at this date, then the date of the last data would be less.
* @return end date of the buffer
*/
Date getBufferEnd() const;
/**
* @brief Check if the buffer is empty
* @return true if the buffer is empty
*/
bool empty() const;
/**
* @brief Clear the buffer; the data is deleted and the start and end dates reset to <i>undef</i>
*/
void clear();
/**
* @brief Add data representing the available data between two dates.
* @param date_start A Date object representing the beginning of an interval (inclusive)
* @param date_end A Date object representing the end of an interval (inclusive)
* @param vecMeteo A vector of vector<MeteoData> objects providing the data
*/
void push(const Date& date_start, const Date& date_end, const std::vector< METEO_SET >& vecMeteo);
//void push(const Date& date, const METEO_SET& vecMeteo);
const std::string toString() const;
//HACK: these should be removed in order to hide the internals!
//but this requires a re-write of MeteoProcessor
//HACK: these should be removed in order to hide the internals! But this requires a re-write of MeteoProcessor
std::vector< METEO_SET >& getBuffer();
void setBufferStart(const Date& date);
void setBufferEnd(const Date& date);
......
......@@ -121,15 +121,15 @@ double Atmosphere::wetBulbTemperature(const double& T, const double& RH, const d
*/
double Atmosphere::blackGlobeTemperature(const double& TA, const double& RH, const double& VW, const double& iswr_dir, const double& iswr_diff, const double& cos_Z)
{
const double a=1, b=1, c=1; //HACK: get real values!
const double S = iswr_dir + iswr_diff;
const double h = a * pow(S, b) * pow(cos_Z, c);
//const double a=1, b=1, c=1; //HACK: get real values!
//const double h = a * pow(S, b) * pow(cos_Z, c);
const double h = 0.315; //personnal communication from S. Amburn
const double emissivity = 0.575 * pow(RH*waterSaturationPressure(TA), 1./7.);
const double B = S * (iswr_dir/(4.*Cst::stefan_boltzmann*cos_Z) + 1.2/Cst::stefan_boltzmann*iswr_diff) + emissivity*Optim::pow4(TA);
const double C = h * pow(VW, 0.58) / 5.3865e-8;
const double Tg = (B + C*TA + 7680000) / (C + 256000);
const double Tg = (B + C*TA + 7680000.) / (C + 256000.);
return Tg;
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment