WSL/SLF GitLab Repository

Commit 04b213bf authored by Mathias Bavay's avatar Mathias Bavay
Browse files

A bug/inconsistency has been fixed in the Winstral algorithm: when pixels are...

A bug/inconsistency has been fixed in the Winstral algorithm: when pixels are excluded because their air temperature is too high for snow, they should not contribute to the min/max erosion/deposition! A new version is now in dveelopment that takes into account the local wind field as well as a threshold on the wind speed. This is not documented as it is still highly experiemental.
parent 42ab416c
......@@ -57,6 +57,8 @@ InterpolationAlgorithm* AlgorithmFactory::getAlgorithm(const std::string& i_algo
return new RyanAlgorithm(i_mi, i_vecArgs, i_algoname, tsm, gdm);
} else if (algoname == "WINSTRAL"){// Winstral wind exposure factor
return new WinstralAlgorithm(i_mi, i_vecArgs, i_algoname, tsm, gdm);
} else if (algoname == "WINSTRAL++"){// Winstral/Liston wind exposure factor
return new WinstralListonAlgorithm(i_mi, i_vecArgs, i_algoname, tsm, gdm);
} else if (algoname == "ODKRIG"){// ordinary kriging
return new OrdinaryKrigingAlgorithm(i_mi, i_vecArgs, i_algoname, tsm, gdm);
} else if (algoname == "ODKRIG_LAPSE"){// ordinary kriging with lapse rate
......@@ -931,6 +933,115 @@ void WinstralAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
}
const double WinstralListonAlgorithm::dmax = 300.;
WinstralListonAlgorithm::WinstralListonAlgorithm(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), base_algo("IDW_LAPSE"), ref_station(),
inputIsAllZeroes(false)
{
const size_t nr_args = vecArgs.size();
if (nr_args==2) {
base_algo = IOUtils::strToUpper( vecArgs[0] );
ref_station = vecArgs[1];
return;
} else if (nr_args!=2)
throw InvalidArgumentException("Wrong number of arguments supplied for the "+algo+" algorithm", AT);
}
double WinstralListonAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
{
//This algorithm is only valid for PSUM (we could add HS later)
if (in_param!=MeteoData::PSUM)
return 0.0;
date = i_date;
param = in_param;
nrOfMeasurments = getData(date, param, vecData, vecMeta);
inputIsAllZeroes = Interpol2D::allZeroes(vecData);
if (inputIsAllZeroes) return 0.99;
if (nrOfMeasurments==0) return 0.0;
if (nrOfMeasurments==1 && ref_station.empty()) { //ie: still using default base_algo
base_algo = "CST";
}
//check that the necessary wind data is available
if (!ref_station.empty()) {
if (!windIsAvailable(vecMeteo, ref_station))
return 0.0;
}
return 0.99;
}
void WinstralListonAlgorithm::initGrid(const DEMObject& dem, Grid2DObject& grid)
{
//initialize precipitation grid with user supplied algorithm (IDW_LAPSE by default)
vector<string> vecArgs2;
mi.getArgumentsForAlgorithm(MeteoData::getParameterName(param), base_algo, vecArgs2);
auto_ptr<InterpolationAlgorithm> algorithm(AlgorithmFactory::getAlgorithm(base_algo, mi, vecArgs2, tsmanager, gridsmanager));
algorithm->getQualityRating(date, param);
algorithm->calculate(dem, grid);
info << algorithm->getInfo();
}
bool WinstralListonAlgorithm::windIsAvailable(const std::vector<MeteoData>& vecMeteo, const std::string& ref_station)
{
if (ref_station.empty()) {
for (size_t ii=0; ii<vecMeteo.size(); ii++) {
const double VW = vecMeteo[ii](MeteoData::VW);
const double DW = vecMeteo[ii](MeteoData::DW);
if (VW!=IOUtils::nodata && DW!=IOUtils::nodata)
return true; //at least one station is enough
}
} else {
double VW, DW;
getSynopticWind(vecMeteo, ref_station, VW, DW);
if (VW!=IOUtils::nodata && DW!=IOUtils::nodata)
return true;
}
return false;
}
void WinstralListonAlgorithm::getSynopticWind(const std::vector<MeteoData>& vecMeteo, const std::string& ref_station, double& VW, double& DW)
{
for (size_t ii=0; ii<vecMeteo.size(); ++ii) {
if (vecMeteo[ii].meta.stationID==ref_station) {
VW = vecMeteo[ii](MeteoData::VW);
DW = vecMeteo[ii](MeteoData::DW);
return;
}
}
}
void WinstralListonAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
{
info.clear(); info.str("");
//if all data points are zero, simply fill the grid with zeroes
if (inputIsAllZeroes) {
Interpol2D::constant(0., dem, grid);
return;
}
//get initial PSUM grid
initGrid(dem, grid);
//get meteo fields interpolation from call back to Meteo2DInterpolator
Grid2DObject ta, dw, vw;
mi.interpolate(date, dem, MeteoData::TA, ta);
mi.interpolate(date, dem, MeteoData::DW, dw);
mi.interpolate(date, dem, MeteoData::VW, vw);
//alter the field with Winstral and the chosen wind direction
Interpol2D::Winstral(dem, ta, dw, vw, dmax, grid);
}
std::string USERInterpolation::getGridFileName() const
{
//HACK: use read2DGrid(grid, MeteoGrid::Parameters, Date) instead?
......
......@@ -473,6 +473,23 @@ class WinstralAlgorithm : public InterpolationAlgorithm {
static const double dmax;
};
class WinstralListonAlgorithm : public InterpolationAlgorithm {
public:
WinstralListonAlgorithm(Meteo2DInterpolator& i_mi,
const std::vector<std::string>& i_vecArgs,
const std::string& i_algo, TimeSeriesManager& i_tsmanager, GridsManager& i_gridsmanager);
virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
private:
void initGrid(const DEMObject& dem, Grid2DObject& grid);
static bool windIsAvailable(const std::vector<MeteoData>& vecMeteo, const std::string& ref_station);
static void getSynopticWind(const std::vector<MeteoData>& vecMeteo, const std::string& ref_station, double& VW, double& DW);
std::string base_algo, ref_station;
bool inputIsAllZeroes;
static const double dmax;
};
/**
* @class USERInterpolation
* @brief Reads user provided gridded data on the disk.
......
......@@ -749,6 +749,38 @@ void Interpol2D::WinstralSX(const DEMObject& dem, const double& dmax, const doub
}
}
void Interpol2D::WinstralSX(const DEMObject& dem, const double& dmax, const Grid2DObject& DW, Grid2DObject& grid)
{
if (!DW.isSameGeolocalization(dem)){
throw IOException("Requested grid DW doesn't match the geolocalization of the DEM", AT);
}
grid.set(dem, IOUtils::nodata);
const double dmin = 0.;
const double bearing_inc = 5.;
const double bearing_width = 30.;
const size_t ncols = dem.getNx(), nrows = dem.getNy();
for (size_t jj = 0; jj<nrows; jj++) {
for (size_t ii = 0; ii<ncols; ii++) {
if (dem(ii,jj)==IOUtils::nodata) continue;
const double in_bearing = DW(ii,jj);
double bearing1 = fmod( in_bearing - bearing_width/2., 360. );
double bearing2 = fmod( in_bearing + bearing_width/2., 360. );
if (bearing1>bearing2) std::swap(bearing1, bearing2);
double sum = 0.;
unsigned short count=0;
for (double bearing=bearing1; bearing<=bearing2; bearing += bearing_inc) {
sum += atan( getTanMaxSlope(dem, dmin, dmax, bearing, ii, jj) );
count++;
}
grid(ii,jj) = (count>0)? sum/(double)count : IOUtils::nodata;
}
}
}
/**
* @brief Alter a precipitation field with the Winstral Sx exposure coefficient
* This implements the wind exposure coefficient (Sx) for one bearing as in
......@@ -772,6 +804,11 @@ void Interpol2D::Winstral(const DEMObject& dem, const Grid2DObject& TA, const do
//compute wind exposure factor
Grid2DObject Sx;
WinstralSX(dem, dmax, in_bearing, Sx);
//don't change liquid precipitation
for (size_t ii=0; ii<Sx.getNx()*Sx.getNy(); ii++) {
if (TA(ii)>Cst::t_water_freezing_pt) Sx(ii)=IOUtils::nodata;
}
//get the scaling parameters
const double min_sx = Sx.grid2D.getMin(); //negative
......@@ -780,7 +817,6 @@ void Interpol2D::Winstral(const DEMObject& dem, const Grid2DObject& TA, const do
//erosion: fully eroded at min_sx
for (size_t ii=0; ii<Sx.getNx()*Sx.getNy(); ii++) {
if (TA(ii)>Cst::t_water_freezing_pt) continue; //don't change liquid precipitation
const double sx = Sx(ii);
if (sx==IOUtils::nodata) continue;
double &val = grid(ii);
......@@ -802,7 +838,56 @@ void Interpol2D::Winstral(const DEMObject& dem, const Grid2DObject& TA, const do
//-> we now have the proper scaling factor so we can deposit in individual cells
const double ratio = sum_erosion/sum_deposition;
for (size_t ii=0; ii<Sx.getNx()*Sx.getNy(); ii++) {
if (TA(ii)>Cst::t_water_freezing_pt) continue; //don't change liquid precipitation
const double sx = Sx(ii);
if (sx==IOUtils::nodata) continue;
double &val = grid(ii);
if (sx>0.) {
const double deposited = ratio * sx/max_sx;
val += deposited;
}
}
}
void Interpol2D::Winstral(const DEMObject& dem, const Grid2DObject& TA, const Grid2DObject& DW, const Grid2DObject& VW, const double& dmax, Grid2DObject& grid)
{
const double vw_thresh = 5.; //m/s
//compute wind exposure factor
Grid2DObject Sx;
WinstralSX(dem, dmax, DW, Sx);
//don't change liquid precipitation
for (size_t ii=0; ii<Sx.getNx()*Sx.getNy(); ii++) {
if (TA(ii)>Cst::t_water_freezing_pt) Sx(ii)=IOUtils::nodata;
}
//get the scaling parameters
const double min_sx = Sx.grid2D.getMin(); //negative
const double max_sx = Sx.grid2D.getMax(); //positive
double sum_erosion=0., sum_deposition=0.;
//erosion: fully eroded at min_sx
for (size_t ii=0; ii<Sx.getNx()*Sx.getNy(); ii++) {
const double sx = Sx(ii);
if (sx==IOUtils::nodata || VW(ii)<vw_thresh) continue; //low wind speed pixels don't contribute to erosion
double &val = grid(ii);
if (sx<0.) {
const double eroded = val * sx/min_sx;
sum_erosion += eroded;
val -= eroded;
}
else { //at this point, we can only compute the sum of deposition
const double deposited = sx/max_sx;
sum_deposition += deposited;
}
}
//no cells can take the eroded mass or no cells even got freezing temperatures
if (sum_deposition==0 || sum_erosion==0) return;
//deposition: garantee mass balance conservation
//-> we now have the proper scaling factor so we can deposit in individual cells
const double ratio = sum_erosion/sum_deposition;
for (size_t ii=0; ii<Sx.getNx()*Sx.getNy(); ii++) {
const double sx = Sx(ii);
if (sx==IOUtils::nodata) continue;
double &val = grid(ii);
......
......@@ -67,6 +67,7 @@ class Interpol2D {
static void RyanWind(const DEMObject& dem, Grid2DObject& VW, Grid2DObject& DW);
static void Winstral(const DEMObject& dem, const Grid2DObject& TA, const double& dmax, const double& in_bearing, Grid2DObject& grid);
static void Winstral(const DEMObject& dem, const Grid2DObject& TA, const Grid2DObject& DW, const Grid2DObject& VW, const double& dmax, Grid2DObject& grid);
static bool allZeroes(const std::vector<double>& vecData);
......@@ -97,6 +98,7 @@ class Interpol2D {
static double depositAroundCell(const DEMObject& dem, const size_t& ii, const size_t& jj, const double& precip, Grid2DObject &grid);
static void WinstralSX(const DEMObject& dem, const double& dmax, const double& in_bearing, Grid2DObject& grid);
static void WinstralSX(const DEMObject& dem, const double& dmax, const Grid2DObject& DW, Grid2DObject& grid);
//weighting methods
static double weightInvDist(const double& d2);
......
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