WSL/SLF GitLab Repository

Commit 284f1213 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

The filters have been slightly simplified, some IOUtils functions optimized...

The filters have been slightly simplified, some IOUtils functions optimized (after reading an interesting post about the handling of rvalues), a bug fixed in PNGIO on Windows, the macros MAX and MIN are not used anymore (relying on std::max and std::min instead). A new data generator is on its way to generate HNW from HS differences (but this will require architectural changes to the DataGenerator).
parent 065e9ead
......@@ -481,7 +481,7 @@ void DEMObject::getPointsBetween(Coords point1, Coords point2, std::vector<GRID_
if(ix1==ix2) {
//special case of vertical alignement
for(int iy=MIN(iy1,iy2); iy<=MAX(iy1,iy2); iy++) {
for(int iy=min(iy1,iy2); iy<=max(iy1,iy2); iy++) {
GRID_POINT_2D pts;
pts.ix = ix1;
pts.iy = iy;
......@@ -496,7 +496,7 @@ void DEMObject::getPointsBetween(Coords point1, Coords point2, std::vector<GRID_
for(int ix=ix1; ix<=ix2; ix++) {
//extension of the line segment (ix, ix+1) along the Y axis
int y1 = (int)floor( a*(double)ix+b );
//const int y2 = MIN( (int)floor( a*((double)ix+1)+b ) , iy2);
//const int y2 = min( (int)floor( a*((double)ix+1)+b ) , iy2);
int y2 = (int)floor( a*((double)ix+1)+b );
if(ix==ix2 && y1==iy2) {
//we don't want to overshoot when reaching the target cell
......@@ -869,21 +869,21 @@ double DEMObject::steepestGradient(double A[4][4]) {
if(A[2][2]!=IOUtils::nodata) {
if(A[1][1]!=IOUtils::nodata)
smax = MAX( smax, fabs(A[2][2] - A[1][1])/(cellsize*sqrt2) );
smax = max( smax, fabs(A[2][2] - A[1][1])/(cellsize*sqrt2) );
if(A[1][2]!=IOUtils::nodata)
smax = MAX( smax, fabs(A[2][2] - A[1][2])/(cellsize) );
smax = max( smax, fabs(A[2][2] - A[1][2])/(cellsize) );
if(A[1][3]!=IOUtils::nodata)
smax = MAX( smax, fabs(A[2][2] - A[1][3])/(cellsize*sqrt2) );
smax = max( smax, fabs(A[2][2] - A[1][3])/(cellsize*sqrt2) );
if(A[2][1]!=IOUtils::nodata)
smax = MAX( smax, fabs(A[2][2] - A[2][1])/(cellsize) );
smax = max( smax, fabs(A[2][2] - A[2][1])/(cellsize) );
if(A[2][3]!=IOUtils::nodata)
smax = MAX( smax, fabs(A[2][2] - A[2][3])/(cellsize) );
smax = max( smax, fabs(A[2][2] - A[2][3])/(cellsize) );
if(A[3][1]!=IOUtils::nodata)
smax = MAX( smax, fabs(A[2][2] - A[3][1])/(cellsize*sqrt2) );
smax = max( smax, fabs(A[2][2] - A[3][1])/(cellsize*sqrt2) );
if(A[3][2]!=IOUtils::nodata)
smax = MAX( smax, fabs(A[2][2] - A[3][2])/(cellsize) );
smax = max( smax, fabs(A[2][2] - A[3][2])/(cellsize) );
if(A[3][3]!=IOUtils::nodata)
smax = MAX( smax, fabs(A[2][2] - A[3][3])/(cellsize*sqrt2) );
smax = max( smax, fabs(A[2][2] - A[3][3])/(cellsize*sqrt2) );
}
if(smax<0.)
......
......@@ -39,6 +39,8 @@ GeneratorAlgorithm* GeneratorAlgorithmFactory::getAlgorithm(const std::string& i
return new UnsworthGenerator(vecArgs, i_algoname);
} else if (algoname == "POT_RADIATION"){
return new PotRadGenerator(vecArgs, i_algoname);
} else if (algoname == "HS_SWE"){
return new HSSweGenerator(vecArgs, i_algoname);
} else {
throw IOException("The generator algorithm '"+algoname+"' is not implemented" , AT);
}
......@@ -430,5 +432,159 @@ double PotRadGenerator::getSolarIndex(const double& ta, const double& rh, const
return karsten_Si;
}
const double HSSweGenerator::soil_albedo = .23; //grass
const double HSSweGenerator::snow_albedo = .85; //snow
const double HSSweGenerator::snow_thresh = .1; //if snow height greater than this threshold -> snow albedo
const double HSSweGenerator::thresh_rh = .7;
const double HSSweGenerator::thresh_Dt = 3.;
const double HSSweGenerator::thresh_iswr = 30.; //minimum radiation to consider it is daylight
const double HSSweGenerator::thresh_solarIndex = .4; //threshold for clear sky
const bool HSSweGenerator::soft = true;
void HSSweGenerator::parse_args(const std::vector<std::string>& vecArgs)
{
if(vecArgs.size()>0) { //incorrect arguments, throw an exception
throw InvalidArgumentException("Wrong number of arguments supplied for the "+algo+" generator", AT);
}
}
bool HSSweGenerator::generate(const size_t& /*param*/, MeteoData& /*md*/)
{//HACK: modify prototype so we can get the full vector + the index of the replacement
return false; //all missing values could be filled
}
bool HSSweGenerator::generate(const size_t& param, std::vector<MeteoData>& vecMeteo)
{
if(param!=MeteoData::HNW)
throw InvalidArgumentException("Trying to use "+algo+" generator on " + MeteoData::getParameterName(param) + " but it can only be applied to HNW!!", AT);
if(vecMeteo.empty()) return true;
//Find first point that is not IOUtils::nodata
size_t last_good = IOUtils::npos;
for (size_t ii=0; ii<vecMeteo.size(); ii++){
if (vecMeteo[ii](MeteoData::HS) != IOUtils::nodata){
last_good = ii;
break;
}
}
if (last_good == IOUtils::npos) //can not find a good point to start
return false;
const double lat = vecMeteo.front().meta.position.getLat();
const double lon = vecMeteo.front().meta.position.getLon();
const double alt = vecMeteo.front().meta.position.getAltitude();
if(lat!=IOUtils::nodata && lon!=IOUtils::nodata && alt!=IOUtils::nodata) {
sun.setLatLon(lat, lon, alt);
sun_ok = true;
} else {
sun_ok = false;
}
bool all_filled = (last_good>0)? false : true;
for(size_t ii=last_good+1; ii<vecMeteo.size(); ii++) {
const double HS_curr = vecMeteo[ii](MeteoData::HS);
if(HS_curr==IOUtils::nodata) continue;
const size_t start_idx = last_good+1;
const double HS_prev = vecMeteo[last_good](MeteoData::HS);
const double HS_delta = HS_curr - HS_prev;
//if HS_delta<0 or ==0, we don't do anything,
//we let the user define a Cst generator to set the remaining nodata to 0
if(HS_delta>0.) {
const double rho = newSnowDensity(vecMeteo[ii]);
const double precip = HS_delta * rho; //in kg/m2 or mm
distributeHNW(precip, start_idx, ii, vecMeteo);
} else
all_filled = false;
last_good=ii;
}
return all_filled;
}
double HSSweGenerator::newSnowDensity(const MeteoData& md) const
{ //Zwart parametrization
const double vw = max(2., md(MeteoData::VW));
const double rh = md(MeteoData::RH);
const double ta = md(MeteoData::TA) - Cst::t_water_triple_pt;
const double beta01=3.28, beta1=0.03, beta02=-0.36, beta2=-0.75, beta3=0.3;
double arg = beta01 + beta1*ta + beta2*asin(sqrt(rh)) + beta3*log10(vw);
if(ta>=-14.)
arg += beta02; // += beta2*ta;
return min( pow(10., arg), 250. ); //limit the density to 250 kg/m3
}
void HSSweGenerator::distributeHNW(const double& precip, const size_t& start_idx, const size_t& end_idx, std::vector<MeteoData>& vecMeteo)
{
const size_t nr_elems = end_idx-start_idx+1;
std::vector<unsigned char> vecScores(end_idx-start_idx+1);
//assign a score for each timestep according to the possibility of precipitation
size_t nr_score1 = 0, nr_score2 = 0, nr_score3 = 0; //count how many in each possible score categories
for(size_t ii=0; ii<nr_elems; ii++) {
unsigned char score = 0;
const double rh = vecMeteo[ii+start_idx](MeteoData::RH);
const double ta = vecMeteo[ii+start_idx](MeteoData::TA);
const double tss = vecMeteo[ii+start_idx](MeteoData::TSS);
const double hs = vecMeteo[ii+start_idx](MeteoData::HS);
//trying to get some radiation information
double iswr = vecMeteo[ii+start_idx](MeteoData::ISWR);
const double rswr = vecMeteo[ii+start_idx](MeteoData::RSWR);
double albedo = .5;
if(iswr!=IOUtils::nodata && rswr!=IOUtils::nodata)
albedo = rswr / iswr;
else if(hs!=IOUtils::nodata)
albedo = (hs>=snow_thresh)? snow_albedo : soil_albedo;
if(iswr==IOUtils::nodata && rswr!=IOUtils::nodata)
iswr = rswr / albedo;
//assign the scores
if(rh!=IOUtils::nodata && rh>=thresh_rh) //enough moisture for precip
score++;
if (ta!=IOUtils::nodata && tss!=IOUtils::nodata && (ta-tss)<=thresh_Dt ) //cloudy sky condition
score++;
if(iswr!=IOUtils::nodata && iswr>thresh_iswr && sun_ok) { //low radiation compared to clear sky
sun.setDate(vecMeteo[ii+start_idx].date.getJulian(true), 0.);
const double p=vecMeteo[ii+start_idx](MeteoData::P);
if(p==IOUtils::nodata)
sun.calculateRadiation(ta, rh, albedo);
else
sun.calculateRadiation(ta, rh, p, albedo);
double toa, direct, diffuse;
sun.getHorizontalRadiation(toa, direct, diffuse);
const double solarIndex = iswr / (direct+diffuse);
if(solarIndex<=thresh_solarIndex)
score++;
}
if(score==1) nr_score1++;
if(score==2) nr_score2++;
if(score==3) nr_score3++;
vecScores[ii] = score;
}
//distribute the precipitation on the time steps that have the highest scores
const unsigned char winning_scores = (nr_score3>0)? 3 : (nr_score2>0)? 2 : (nr_score1>0)? 1 : 0;
const size_t nr_winning_scores = (nr_score3>0)? nr_score3 : (nr_score2>0)? nr_score2 : (nr_score1>0)? nr_score1 : nr_elems;
const double precip_increment = precip / double(nr_winning_scores);
for(size_t ii=0; ii<nr_elems; ii++) {
if(winning_scores>0 && vecScores[ii]==winning_scores) {
vecMeteo[ii+start_idx](MeteoData::HNW) = precip_increment;
} else
vecMeteo[ii+start_idx](MeteoData::HNW) = 0.; //all other time steps set to 0
}
}
} //namespace
......@@ -225,6 +225,24 @@ class PotRadGenerator : public GeneratorAlgorithm {
static const double soil_albedo, snow_albedo, snow_thresh; //to try using rswr if not iswr is given
};
class HSSweGenerator : public GeneratorAlgorithm {
public:
HSSweGenerator(const std::vector<std::string>& vecArgs, const std::string& i_algo)
: GeneratorAlgorithm(vecArgs, i_algo), sun(), sun_ok(false) { parse_args(vecArgs); }
bool generate(const size_t& param, MeteoData& md);
bool generate(const size_t& param, std::vector<MeteoData>& vecMeteo);
private:
void parse_args(const std::vector<std::string>& vecArgs);
double newSnowDensity(const MeteoData& md) const;
void distributeHNW(const double& precip, const size_t& start_idx, const size_t& end_idx, std::vector<MeteoData>& vecMeteo);
SunObject sun;
static const double soil_albedo, snow_albedo, snow_thresh; //to try using rswr if not iswr is given
static const double thresh_rh, thresh_Dt, thresh_iswr, thresh_solarIndex;
static const bool soft;
bool sun_ok;
};
} //end namespace mio
#endif
......@@ -123,34 +123,32 @@ void copy_file(const std::string& src, const std::string& dest)
fout.close();
}
std::string cleanPath(const std::string& in_path, const bool& resolve)
std::string cleanPath(std::string in_path, const bool& resolve)
{
if(!resolve) { //do not resolve links, relative paths, etc
std::string out_path(in_path);
std::replace(out_path.begin(), out_path.end(), '\\', '/');
return out_path;
std::replace(in_path.begin(), in_path.end(), '\\', '/');
return in_path;
} else {
#ifdef WIN32
//if this would not suffice, see http://pdh11.blogspot.ch/2009/05/pathcanonicalize-versus-what-it-says-on.html
char **ptr = NULL;
char *out_buff = (char*)calloc(MAX_PATH, sizeof(char));
const DWORD status = GetFullPathName(in_path.c_str(), MAX_PATH, out_buff, ptr);
std::string out_path = (status!=0 && status<=MAX_PATH)? out_buff : in_path;
if(status!=0 && status<=MAX_PATH) in_path = out_buff;
free(out_buff);
std::replace(out_path.begin(), out_path.end(), '\\', '/');
return out_path;
std::replace(in_path.begin(), in_path.end(), '\\', '/');
return in_path;
#else //POSIX
std::string out_path(in_path);
std::replace(out_path.begin(), out_path.end(), '\\', '/');
std::replace(in_path.begin(), in_path.end(), '\\', '/');
char *real_path = realpath(out_path.c_str(), NULL); //POSIX
char *real_path = realpath(in_path.c_str(), NULL); //POSIX
if(real_path!=NULL) {
const std::string tmp(real_path);
free(real_path);
return tmp;
} else
return out_path; //something failed in realpath, keep it as it is
return in_path; //something failed in realpath, keep it as it is
#endif
}
}
......@@ -221,16 +219,16 @@ void toLower(std::string& str) {
std::transform(str.begin(), str.end(), str.begin(), ::tolower);
}
std::string strToUpper(const std::string &str) {
std::string strcopy(str.size(), 0);
std::transform(str.begin(), str.end(), strcopy.begin(), ::toupper);
return strcopy;
std::string strToUpper(std::string str) {
//based on http://cpp-next.com/archive/2009/08/want-speed-pass-by-value/
//it is better to let the compiler copy (or not!) the parameters
std::transform(str.begin(), str.end(), str.begin(), ::toupper);
return str;
}
std::string strToLower(const std::string &str) {
std::string strcopy(str.size(), 0);
std::transform(str.begin(), str.end(), strcopy.begin(), ::tolower);
return strcopy;
std::string strToLower(std::string str) {
std::transform(str.begin(), str.end(), str.begin(), ::tolower);
return str;
}
bool isNumeric(std::string str, const unsigned int& nBase)
......@@ -734,7 +732,7 @@ size_t seek(const Date& soughtdate, const std::vector<MeteoData>& vecM, const bo
const double curr_date = soughtdate.getJulian(true);
const double raw_pos = (curr_date-start_date) / (end_date-start_date) * static_cast<double>(max_idx); //always >=0
const size_t start_idx = (size_t)floor(raw_pos*.9);
const size_t end_idx = MIN( (size_t)ceil(raw_pos*1.1), max_idx);
const size_t end_idx = std::min( (size_t)ceil(raw_pos*1.1), max_idx);
//first and last index of the search interval, either using our initial guess or the full vector
size_t first = (curr_date >= vecM[start_idx].date.getJulian(true))? start_idx : 0;
......
......@@ -36,13 +36,6 @@
#include <limits>
#include <cmath>
#ifndef MAX
#define MAX(x,y) (((x) < (y)) ? (y) : (x))
#endif
#ifndef MIN
#define MIN(x,y) (((x) < (y)) ? (x) : (y))
#endif
#ifndef C_TO_K
#define C_TO_K( T ) ( T + Cst::t_water_freezing_pt ) // degree Celsius to kelvin
#endif
......@@ -144,7 +137,7 @@ namespace IOUtils {
* @param in_path the path string to cleanup
* @param resolve resolve links, convert relative paths, etc? (default=false)
*/
std::string cleanPath(const std::string& in_path, const bool& resolve=false);
std::string cleanPath(std::string in_path, const bool& resolve=false);
/**
* @brief returns the extension part of a given filename.
......@@ -204,9 +197,9 @@ namespace IOUtils {
std::string &key, std::string &value, const bool& setToUpperCase=false);
void toUpper(std::string& str);
std::string strToUpper(const std::string &str);
std::string strToUpper(std::string str);
void toLower(std::string& str);
std::string strToLower(const std::string &str);
std::string strToLower(std::string str);
bool isNumeric(std::string input, const unsigned int& nBase=10);
size_t readLineToVec(const std::string& line_in, std::vector<double>& vec_data);
size_t readLineToVec(const std::string& line_in, std::vector<std::string>& vecString);
......
......@@ -70,8 +70,8 @@ void MeteoProcessor::getWindowSize(ProcessingProperties& o_properties) const
void MeteoProcessor::compareProperties(const ProcessingProperties& newprop, ProcessingProperties& current)
{
current.points_before = MAX(current.points_before, newprop.points_before);
current.points_after = MAX(current.points_after, newprop.points_after);
current.points_before = max(current.points_before, newprop.points_before);
current.points_after = max(current.points_after, newprop.points_after);
if (newprop.time_before > current.time_before)
current.time_before = newprop.time_before;
......
......@@ -23,7 +23,7 @@ using namespace std;
namespace mio {
FilterMAD::FilterMAD(const std::vector<std::string>& vec_args) : WindowedFilter("MAD")
FilterMAD::FilterMAD(const std::vector<std::string>& vec_args, const std::string& name) : WindowedFilter(name)
{
parse_args(vec_args);
......
......@@ -43,7 +43,7 @@ namespace mio {
class FilterMAD : public WindowedFilter {
public:
FilterMAD(const std::vector<std::string>& vec_args);
FilterMAD(const std::vector<std::string>& vec_args, const std::string& name);
virtual void process(const unsigned int& param, const std::vector<MeteoData>& ivec,
std::vector<MeteoData>& ovec);
......
......@@ -21,8 +21,8 @@ using namespace std;
namespace mio {
FilterMax::FilterMax(const std::vector<std::string>& vec_args)
: FilterBlock("MAX"), max_val(0.), max_soft(0.), is_soft(true)
FilterMax::FilterMax(const std::vector<std::string>& vec_args, const std::string& name)
: FilterBlock(name), max_val(0.), max_soft(0.), is_soft(true)
{
parse_args(vec_args);
properties.stage = ProcessingProperties::both; //for the rest: default values
......
......@@ -42,7 +42,7 @@ namespace mio {
class FilterMax : public FilterBlock {
public:
FilterMax(const std::vector<std::string>& vec_args);
FilterMax(const std::vector<std::string>& vec_args, const std::string& name);
virtual void process(const unsigned int& param, const std::vector<MeteoData>& ivec,
std::vector<MeteoData>& ovec);
......
......@@ -22,7 +22,7 @@ using namespace std;
namespace mio {
FilterMeanAvg::FilterMeanAvg(const std::vector<std::string>& vec_args) : WindowedFilter("MEAN_AVG")
FilterMeanAvg::FilterMeanAvg(const std::vector<std::string>& vec_args, const std::string& name) : WindowedFilter(name)
{
parse_args(vec_args);
......
......@@ -49,7 +49,7 @@ namespace mio {
class FilterMeanAvg : public WindowedFilter {
public:
FilterMeanAvg(const std::vector<std::string>& vec_args);
FilterMeanAvg(const std::vector<std::string>& vec_args, const std::string& name);
virtual void process(const unsigned int& param, const std::vector<MeteoData>& ivec,
std::vector<MeteoData>& ovec);
......
......@@ -22,7 +22,7 @@ using namespace std;
namespace mio {
FilterMedianAvg::FilterMedianAvg(const std::vector<std::string>& vec_args) : WindowedFilter("MEAN_AVG")
FilterMedianAvg::FilterMedianAvg(const std::vector<std::string>& vec_args, const std::string& name) : WindowedFilter(name)
{
parse_args(vec_args);
......
......@@ -52,7 +52,7 @@ namespace mio {
class FilterMedianAvg : public WindowedFilter {
public:
FilterMedianAvg(const std::vector<std::string>& vec_args);
FilterMedianAvg(const std::vector<std::string>& vec_args, const std::string& name);
virtual void process(const unsigned int& param, const std::vector<MeteoData>& ivec,
std::vector<MeteoData>& ovec);
......
......@@ -21,8 +21,8 @@ using namespace std;
namespace mio {
FilterMin::FilterMin(const std::vector<std::string>& vec_args)
: FilterBlock("MIN"), min_val(0.), min_soft(0.), is_soft(true)
FilterMin::FilterMin(const std::vector<std::string>& vec_args, const std::string& name)
: FilterBlock(name), min_val(0.), min_soft(0.), is_soft(true)
{
parse_args(vec_args);
properties.stage = ProcessingProperties::both; //for the rest: default values
......
......@@ -44,7 +44,7 @@ namespace mio {
class FilterMin : public FilterBlock {
public:
FilterMin(const std::vector<std::string>& vec_args);
FilterMin(const std::vector<std::string>& vec_args, const std::string& name);
virtual void process(const unsigned int& param, const std::vector<MeteoData>& ivec,
std::vector<MeteoData>& ovec);
......
......@@ -21,8 +21,8 @@ using namespace std;
namespace mio {
FilterMinMax::FilterMinMax(const std::vector<std::string>& vec_args)
: FilterBlock("MIN_MAX"), min_val(0.), max_val(0.), min_soft(0.), max_soft(0.), is_soft(true)
FilterMinMax::FilterMinMax(const std::vector<std::string>& vec_args, const std::string& name)
: FilterBlock(name), min_val(0.), max_val(0.), min_soft(0.), max_soft(0.), is_soft(true)
{
parse_args(vec_args);
......
......@@ -48,7 +48,7 @@ namespace mio {
class FilterMinMax : public FilterBlock {
public:
FilterMinMax(const std::vector<std::string>& vec_args);
FilterMinMax(const std::vector<std::string>& vec_args, const std::string& name);
virtual void process(const unsigned int& param, const std::vector<MeteoData>& ivec,
std::vector<MeteoData>& ovec);
......
......@@ -22,7 +22,8 @@ using namespace std;
namespace mio {
FilterRate::FilterRate(const std::vector<std::string>& vec_args) : FilterBlock("RATE"), min_rate_of_change(0.), max_rate_of_change(0.)
FilterRate::FilterRate(const std::vector<std::string>& vec_args, const std::string& name)
: FilterBlock(name), min_rate_of_change(0.), max_rate_of_change(0.)
{
parse_args(vec_args);
properties.stage = ProcessingProperties::both; //for the rest: default values
......
......@@ -43,7 +43,7 @@ namespace mio {
*/
class FilterRate : public FilterBlock {
public:
FilterRate(const std::vector<std::string>& vec_args);
FilterRate(const std::vector<std::string>& vec_args, const std::string& name);
virtual void process(const unsigned int& param, const std::vector<MeteoData>& ivec,
std::vector<MeteoData>& ovec);
......
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