WSL/SLF GitLab Repository

Commit 1ecb5269 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

Finally, adding PSUM_PH to the MeteoData enum. And implementing an...

Finally, adding PSUM_PH to the MeteoData enum. And implementing an interpolation algorithm that generates distributed PSUM_PH by calling the parametrizations on each cell.
parent c0b10142
......@@ -63,6 +63,8 @@ InterpolationAlgorithm* AlgorithmFactory::getAlgorithm(const std::string& i_algo
return new LapseOrdinaryKrigingAlgorithm(i_mi, i_vecArgs, i_algoname, tsm, gdm);
} else if (algoname == "USER"){// read user provided grid
return new USERInterpolation(i_mi, i_vecArgs, i_algoname, tsm, gdm);
} else if (algoname == "PPHASE"){// precipitation phase parametrization
return new PPHASEInterpolation(i_mi, i_vecArgs, i_algoname, tsm, gdm);
} else if (algoname == "PSUM_SNOW"){// precipitation interpolation according to (Magnusson, 2010)
return new SnowPSUMInterpolation(i_mi, i_vecArgs, i_algoname, tsm, gdm);
} else {
......@@ -614,14 +616,15 @@ double ListonWindAlgorithm::getQualityRating(const Date& i_date, const MeteoData
}
}
if (nrOfMeasurments==0)
return 0.0;
if ( (param==MeteoData::VW && Interpol2D::allZeroes(vecDataVW)) ||
(param==MeteoData::DW && Interpol2D::allZeroes(vecDataDW)) ) {
inputIsAllZeroes = true;
return 0.9;
}
if (nrOfMeasurments==0)
return 0.0;
if( (nrOfMeasurments<vecDataVW.size()/2) || ( nrOfMeasurments<2 ) )
return 0.6;
......@@ -673,14 +676,15 @@ double RyanAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Para
}
}
if (nrOfMeasurments==0)
return 0.0;
if ( (param==MeteoData::VW && Interpol2D::allZeroes(vecDataVW)) ||
(param==MeteoData::DW && Interpol2D::allZeroes(vecDataDW)) ) {
inputIsAllZeroes = true;
return 0.9;
}
if (nrOfMeasurments==0)
return 0.0;
if (nrOfMeasurments<2)
return 0.6;
......@@ -972,6 +976,68 @@ void USERInterpolation::calculate(const DEMObject& dem, Grid2DObject& grid)
}
double PPHASEInterpolation::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
{
date = i_date;
param = in_param;
const size_t nArgs = vecArgs.size();
if (nArgs<1 || IOUtils::isNumeric(vecArgs[0]))
throw InvalidArgumentException("Wrong arguments supplied to the "+algo+" interpolation. Please provide the method to use and its arguments!", AT);
const std::string user_algo = IOUtils::strToUpper(vecArgs[0]);
if (user_algo=="THRESH") {
if (nArgs!=2)
throw InvalidArgumentException("Wrong number of arguments supplied to the "+algo+" interpolation for the "+user_algo+" method", AT);
IOUtils::convertString(fixed_thresh, vecArgs[1]);
model = THRESH;
} else if (user_algo=="RANGE") {
if (nArgs!=3)
throw InvalidArgumentException("Wrong number of arguments supplied to the "+algo+" interpolation for the "+user_algo+" method", AT);
double range_thresh1, range_thresh2;
IOUtils::convertString(range_thresh1, vecArgs[1]);
IOUtils::convertString(range_thresh2, vecArgs[2]);
if (range_thresh1==range_thresh2)
throw InvalidArgumentException(algo+" interpolation, "+user_algo+" method: the two provided threshold must be different", AT);
if (range_thresh1>range_thresh2)
std::swap(range_thresh1, range_thresh2);
range_start = range_thresh1;
range_norm = 1. / (range_thresh2-range_thresh1);
model = RANGE;
} else
throw InvalidArgumentException("Unknown parametrization \""+user_algo+"\" supplied to the "+algo+" interpolation", AT);
return 0.1;
}
void PPHASEInterpolation::calculate(const DEMObject& dem, Grid2DObject& grid)
{
info.clear(); info.str("");
Grid2DObject ta;
mi.interpolate(date, dem, MeteoData::TA, ta); //get TA interpolation from call back to Meteo2DInterpolator
grid.set(dem, IOUtils::nodata);
const size_t nrCells = dem.getNx() * dem.getNy();
if (model==THRESH) {
for (size_t ii=0; ii<nrCells; ii++) {
const double TA=ta(ii);
if (TA==IOUtils::nodata) continue;
grid(ii) = (TA>=fixed_thresh)? 1. : 0.;
}
} else if (model==RANGE) {
for (size_t ii=0; ii<nrCells; ii++) {
const double TA=ta(ii);
if (TA==IOUtils::nodata) continue;
const double tmp_rainfraction = range_norm * (TA - range_start);
grid(ii) = (tmp_rainfraction>1)? 1. : (tmp_rainfraction<0.)? 0. : tmp_rainfraction;
}
}
}
double SnowPSUMInterpolation::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
{
date = i_date;
......
......@@ -89,6 +89,7 @@ class Meteo2DInterpolator; // forward declaration, cyclic header include
* - RYAN: the wind direction is interpolated using IDW and then altered depending on the local slope (see RyanAlgorithm)
* - WINSTRAL: the solid precipitation is redistributed by wind according to (Winstral, 2002) (see WinstralAlgorithm)
* - PSUM_SNOW: precipitation interpolation according to (Magnusson, 2011) (see SnowPSUMInterpolation)
* - PPHASE: precipitation phase parametrization performed at each cell (see PPHASEInterpolation)
* - ODKRIG: ordinary kriging (see OrdinaryKrigingAlgorithm)
* - ODKRIG_LAPSE: ordinary kriging with lapse rate (see LapseOrdinaryKrigingAlgorithm)
* - USER: user provided grids to be read from disk (if available, see USERInterpolation)
......@@ -500,6 +501,39 @@ class USERInterpolation : public InterpolationAlgorithm {
std::string filename;
};
/**
* @class PPHASEInterpolation
* @brief Precipitation phase splitting generation
* This does not interpolate any measured precipitation phase but generates it for each point based on parametrizations, similarly to the PPHASE generator
* (see PPhaseGenerator).
*
* The methods that are offered are currently the following:
* - THRESH: a provided fixed air temperature threshold splits precipitation as either fully solid or fully liquid
* - RANGE: two air temperature thresholds provide the lower and upper range for fully solid / fully liquid precipitation.
* Within the provided range, a linear transition is assumed.
* @code
* PSUM::algorithms = PPHASE
* PSUM::pphase = THRESH 274.35
* @endcode
*/
class PPHASEInterpolation : public InterpolationAlgorithm {
public:
PPHASEInterpolation(Meteo2DInterpolator& i_mi,
const std::vector<std::string>& i_vecArgs,
const std::string& i_algo, TimeSeriesManager& i_tsmanager, GridsManager& i_gridsmanager)
: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, i_tsmanager, i_gridsmanager),
model(THRESH), fixed_thresh(IOUtils::nodata), range_start(IOUtils::nodata), range_norm(IOUtils::nodata) {}
virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
private:
typedef enum PARAMETRIZATION {
THRESH,
RANGE
} parametrization;
parametrization model;
double fixed_thresh, range_start, range_norm;
};
/**
* @class SnowPSUMInterpolation
* @brief Precipitation distribution according to the local slope and curvature.
......
......@@ -104,6 +104,7 @@ bool MeteoData::initStaticData()
static_meteoparamname[ILWR] = "ILWR";
static_meteoparamname[TAU_CLD]= "TAU_CLD";
static_meteoparamname[PSUM] = "PSUM";
static_meteoparamname[PSUM_PH] = "PSUM_PH";
s_default_paramname.push_back("P");
s_default_paramname.push_back("TA");
......@@ -119,6 +120,7 @@ bool MeteoData::initStaticData()
s_default_paramname.push_back("ILWR");
s_default_paramname.push_back("TAU_CLD");
s_default_paramname.push_back("PSUM");
s_default_paramname.push_back("PSUM_PH");
return true;
}
......
......@@ -117,7 +117,8 @@ class MeteoData {
ILWR, ///< Incoming long wave radiation (downwelling)
TAU_CLD, ///< Cloud transmissivity or ISWR/ISWR_clear_sky
PSUM, ///< Water equivalent of precipitations, either solid or liquid
lastparam=PSUM};
PSUM_PH, ///< Precipitation phase: between 0 (fully solid) and 1(fully liquid)
lastparam=PSUM_PH};
static const std::string& getParameterName(const size_t& parindex);
......
......@@ -231,9 +231,7 @@ void SMETIO::identify_fields(const std::vector<std::string>& fields, std::vector
}
//specific key mapping
if (key == "PSUM") {
indexes.push_back(md.getParameterIndex("PSUM"));
} else if (key == "OSWR") {
if (key == "OSWR") {
indexes.push_back(md.getParameterIndex("RSWR"));
} else if (key == "OLWR") {
md.addParameter("OLWR");
......@@ -610,6 +608,9 @@ void SMETIO::getFormatting(const size_t& param, int& prec, int& width)
} else if (param == MeteoData::PSUM){
prec = 3;
width = 6;
} else if (param == MeteoData::PSUM_PH){
prec = 3;
width = 4;
} else if (param == MeteoData::HS){
prec = 3;
width = 8;
......
Markdown is supported
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