WSL/SLF GitLab Repository

Commit 9d2dae3e authored by Mathias Bavay's avatar Mathias Bavay
Browse files

A tricky bug has been found (by Leonardo, once again!) when accumulating...

A tricky bug has been found (by Leonardo, once again!) when accumulating precipitation beyond the original data set: the accumulation was creating "random" precipitation by accumulating data that did not exist. This has been found and fixed (a few conditions were missing to ensure that no etrapolations would be performed when summing partial contributions).
parent c035c2f0
......@@ -66,7 +66,7 @@ std::string Accumulate::toString() const
size_t Accumulate::findStartOfPeriod(const std::vector<MeteoData>& vecM, const size_t& index, const Date& dateStart)
{
size_t start_idx = IOUtils::npos;
for (size_t idx=index; idx--> 0; ) {
for (size_t idx=index+1; idx--> 0; ) { //because idx gets decremented right away
if ( vecM[idx].date <= dateStart) {
start_idx = idx;
break;
......@@ -92,9 +92,9 @@ double Accumulate::complexSampling(const std::vector<MeteoData>& vecM, const siz
{//to keep in mind: start_idx is last index <= dateStart and index is first index >= resampling_date
double sum = IOUtils::nodata;
//resample begining point, in the [start_idx ; start_idx+1] interval
const double start_value = partialAccumulateAtRight(vecM, paramindex, start_idx, dateStart);
if (start_value!=IOUtils::nodata)
sum=start_value;
const double start_val = partialAccumulateAtRight(vecM, paramindex, start_idx, dateStart);
if (start_val!=IOUtils::nodata)
sum = start_val;
else if (strict) return IOUtils::nodata;
//sum all whole periods AFTER the begining point, in the [start_idx+2 ; index-1] interval
......@@ -107,11 +107,20 @@ double Accumulate::complexSampling(const std::vector<MeteoData>& vecM, const siz
}
//resample end point, in the [index-1 ; index] interval
const double end_val = partialAccumulateAtLeft(vecM, paramindex, index-1, resampling_date);
if (end_val!=IOUtils::nodata) {
if (sum!=IOUtils::nodata) sum += end_val;
else sum = end_val;
} else if (strict) return IOUtils::nodata;
if (vecM[index].date>=resampling_date) { //we have enough data points to cover the whole accumulating period
const double end_val = partialAccumulateAtLeft(vecM, paramindex, index-1, resampling_date);
if (end_val!=IOUtils::nodata) {
if (sum!=IOUtils::nodata) sum += end_val;
else sum = end_val;
} else if (strict) return IOUtils::nodata;
} else { //some data points are missing to cover the whole accumulating period
if (strict) return IOUtils::nodata;
const double end_val = vecM[index](paramindex);
if (end_val!=IOUtils::nodata) {
if (sum!=IOUtils::nodata) sum += end_val;
else sum = end_val;
}
}
return sum;
}
......@@ -129,8 +138,8 @@ void Accumulate::resample(const size_t& index, const ResamplingPosition& positio
md(paramindex) = IOUtils::nodata;
const Date resampling_date = md.date;
const Date dateStart(resampling_date.getJulian() - accumulate_period, resampling_date.getTimeZone());
const Date resampling_date( md.date );
const Date dateStart( resampling_date.getJulian() - accumulate_period, resampling_date.getTimeZone() );
const size_t start_idx = findStartOfPeriod(vecM, index, dateStart);
if (start_idx==IOUtils::npos) {//No acceptable starting point found
std::cerr << "[W] Could not accumulate " << vecM.at(0).getNameForParameter(paramindex) << ": ";
......@@ -138,6 +147,12 @@ void Accumulate::resample(const size_t& index, const ResamplingPosition& positio
return;
}
if (position==ResamplingAlgorithms::end && start_idx>=(vecM.size()-1)) { //the first element before "position" is the last of the vector
//the period that should be accumulated is not even in vecM
if (!strict) md(paramindex) = 0.;
return;
}
if ((index - start_idx) <= 1) {//easy upsampling when start & stop are in the same input time step
//upsampling (for example, generate 15min values from hourly data)
const double sum = easySampling(vecM, paramindex, index, start_idx, dateStart, resampling_date); //if resampling was unsuccesful, sum==IOUtils::nodata
......
......@@ -66,9 +66,10 @@ double ResamplingAlgorithms::partialAccumulateAtLeft(const std::vector<MeteoData
const double jul1 = vecM[pos].date.getJulian(true);
const double jul2 = vecM[end].date.getJulian(true);
const double jul_curr = curr_date.getJulian(true);
const double left_accumulation = linearInterpolation(jul1, 0., jul2, valend, curr_date.getJulian(true));
if (jul_curr>jul2) return IOUtils::nodata; //jul1 and jul2 should frame jul_curr, if not some data points are missing
const double left_accumulation = linearInterpolation(jul1, 0., jul2, valend, jul_curr);
return left_accumulation;
}
......@@ -84,7 +85,9 @@ double ResamplingAlgorithms::partialAccumulateAtRight(const std::vector<MeteoDat
const double jul1 = vecM[pos].date.getJulian(true);
const double jul2 = vecM[end].date.getJulian(true);
const double jul_curr = curr_date.getJulian(true);
if (jul_curr<jul1) return IOUtils::nodata; //jul1 and jul2 should frame jul_curr, if not some data points are missing
const double left_accumulation = linearInterpolation(jul1, 0., jul2, valend, curr_date.getJulian(true));
return valend - left_accumulation;
......@@ -107,7 +110,7 @@ void ResamplingAlgorithms::getNearestValidPts(const size_t& pos, const size_t& p
indexP2=IOUtils::npos;
const Date dateStart = resampling_date - window_size;
for (size_t ii=pos; ii-- >0; ) {
for (size_t ii=pos; ii-- >0; ) { //because idx gets decremented right away
if (vecM[ii].date < dateStart) break;
if (vecM[ii](paramindex) != IOUtils::nodata) {
indexP1 = ii;
......
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