WSL/SLF GitLab Repository

Commit 6f86bda4 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

A new "NOISE" filter has been added, to add noise to the inputs. An improved...

A new "NOISE" filter has been added, to add noise to the inputs. An improved check and error message has been implemented in FilterSuppr. An updated reference added to the ESOLIP generator.
parent 58999ecb
......@@ -389,7 +389,10 @@ class AllSkySWGenerator : public GeneratorAlgorithm {
* @brief Generate precipitation from changes in snow height.
* This implements the approach laid out in
* Mair et al., <i>" ESOLIP–estimate of solid and liquid precipitation at sub-daily time resolution by combining snow height
* and rain gauge measurements"</i>, Hydrology and Earth System Sciences Discussions, <b>10(7)</b>, 8683-8714, 2013.
* and rain gauge measurements"</i>, Hydrology and Earth System Sciences Discussions, <b>10(7)</b>, 8683-8714, 2013. or
* Mair E., Leitinger G., Della Chiesa S., Niedrist G., Tappeiner U., Bertoldi G., <i>"A simple method to combine snow height and
* meteorological observations to estimate winter precipitation at sub-daily resolution"</i>, Journal of Hydrological Sciences,
* in revision, 2015.
* The snow density relies on Zwart, <i>"Significance of new-snow properties for snowcover development"</i>,master's thesis,
* Institute for Marine and Atmospheric Research, University of Utrecht, 78 pp, 2007.
*
......
......@@ -22,6 +22,7 @@ SET(meteoFilters_sources
meteoFilters/ProcUndercatch_Forland.cc
meteoFilters/ProcAdd.cc
meteoFilters/ProcMult.cc
meteoFilters/ProcNoise.cc
meteoFilters/ProcExpSmoothing.cc
meteoFilters/ProcWMASmoothing.cc
meteoFilters/FilterBlock.cc
......
......@@ -64,7 +64,8 @@ void FilterSuppr::parse_args(std::vector<std::string> vec_args)
throw InvalidArgumentException("Wrong number of arguments for filter " + getName(), AT);
if (nrArgs==1) {
IOUtils::convertString(range, vec_args[0]);
if (!IOUtils::convertString(range, vec_args[0]))
throw InvalidArgumentException("Invalid range \""+vec_args[0]+"\" specified for the "+getName()+" filter.", AT);
if (range<0. || range>1.)
throw InvalidArgumentException("Wrong range for filter " + getName() + ", it should be between 0 and 1", AT);
}
......
/***********************************************************************************/
/* Copyright 2014 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
/***********************************************************************************/
/* This file is part of MeteoIO.
MeteoIO is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
MeteoIO is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with MeteoIO. If not, see <http://www.gnu.org/licenses/>.
*/
#include <ctime>
#include <cstdlib>
#include <meteoio/meteoFilters/ProcNoise.h>
#include <meteoio/MathOptim.h>
#include <meteoio/meteoLaws/Meteoconst.h>
using namespace std;
namespace mio {
ProcNoise::ProcNoise(const std::vector<std::string>& vec_args, const std::string& name)
: ProcessingBlock(name), range(IOUtils::nodata), distribution(), type()
{
parse_args(vec_args);
properties.stage = ProcessingProperties::first;
}
void ProcNoise::process(const unsigned int& param, const std::vector<MeteoData>& ivec,
std::vector<MeteoData>& ovec)
{
srand( static_cast<unsigned int>(time(NULL)) );
ovec = ivec;
if (type=='a' && distribution=='u') uniform_add(param, ovec);
else if (type=='m' && distribution=='u') uniform_mult(param, ovec);
else if (type=='a' && distribution=='n') normal_add(param, ovec);
else if (type=='m' && distribution=='n') normal_mult(param, ovec);
}
//add a number between -range/2 and +range/2
void ProcNoise::uniform_add(const unsigned int& param, std::vector<MeteoData>& ovec) const
{
for (size_t ii=0; ii<ovec.size(); ii++){
double& tmp = ovec[ii](param);
if (tmp == IOUtils::nodata) continue; //preserve nodata values
tmp += (2.*static_cast<double>(rand())/(RAND_MAX) - 1.) * range;
}
}
//add a number between -range*val/2 and +range*val/2
void ProcNoise::uniform_mult(const unsigned int& param, std::vector<MeteoData>& ovec) const
{
for (size_t ii=0; ii<ovec.size(); ii++){
double& tmp = ovec[ii](param);
if (tmp == IOUtils::nodata) continue; //preserve nodata values
tmp *= (1. + (2.*static_cast<double>(rand())/(RAND_MAX)-1.) * range);
}
}
void ProcNoise::normal_add(const unsigned int& param, std::vector<MeteoData>& ovec) const
{
for (size_t ii=0; ii<ovec.size(); ii++){
double& tmp = ovec[ii](param);
if (tmp == IOUtils::nodata) continue; //preserve nodata values
tmp += getBoxMuller() * range;
}
}
void ProcNoise::normal_mult(const unsigned int& param, std::vector<MeteoData>& ovec) const
{
for (size_t ii=0; ii<ovec.size(); ii++){
double& tmp = ovec[ii](param);
if (tmp == IOUtils::nodata) continue; //preserve nodata values
tmp += (1. + getBoxMuller() * range);
}
}
// Box–Muller method for normally distributed random numbers
// see https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
// This generate a normally distributed signal of mean=0 and std_dev=1
// For numerical reasons, the extremes will always be less than 7 std_dev from the mean
double ProcNoise::getBoxMuller() const
{
const double U = static_cast<double>(rand())/(RAND_MAX);
const double V = static_cast<double>(rand())/(RAND_MAX);
return Optim::fastSqrt_Q3(-2.*log(U)) * cos(2.*Cst::PI*V);
}
void ProcNoise::parse_args(std::vector<std::string> vec_args)
{
if (vec_args.size()!=3){
throw InvalidArgumentException("Wrong number of arguments for filter " + getName() + ": it requires a type, a distribution and a range", AT);
}
const string type_str=IOUtils::strToUpper( vec_args[0] );
if (type_str=="ADD")
type='a';
else if (type_str=="MULT")
type='m';
else
throw InvalidArgumentException("Invalid type \""+type_str+"\" specified for the "+getName()+" filter", AT);
const string distribution_str=IOUtils::strToUpper( vec_args[1] );
if (distribution_str=="UNIFORM")
distribution='u';
else if (distribution_str=="NORMAL")
distribution='n';
else
throw InvalidArgumentException("Invalid distribution \""+distribution_str+"\" specified for the "+getName()+" filter", AT);
if (!IOUtils::convertString(range, vec_args[2]))
throw InvalidArgumentException("Invalid range specified for the "+getName()+" filter", AT);
}
} //end namespace
/***********************************************************************************/
/* Copyright 2014 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
/***********************************************************************************/
/* This file is part of MeteoIO.
MeteoIO is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
MeteoIO is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with MeteoIO. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef __NOISE_H__
#define __NOISE_H__
#include <meteoio/meteoFilters/ProcessingBlock.h>
#include <vector>
#include <string>
namespace mio {
/**
* @class ProcNoise
* @brief Generate a noise signal to modify the input.
* The noise signal is either added ("add") to the input or used as a fraction and multiplied by the input signal ("mult"). This filter
* always takes three arguments: a <i>type</i> specifying if the noise is added or multiplied, a <i>distribution</i> specifying
* the random numbers distribution and a <i>range</i> that is the scaling factor to apply.
* The following random distributions are supported:
* + "uniform": the range represents the maximum amplitude of the noise;
* + "normal": the range represents the standard deviation of the noise.
*
* @note When multiplying the input signal by the noise, the range is interpreted as a fraction. For example, to modify the input by
* +/-10%, select "mult" with "uniform" and "0.1" as a range.
* @ingroup processing
* @author Mathias Bavay
* @date 2015-06-12
* @code
* #add a normally distributed noise (mean=0) of standard deviation 5 to TA
* TA::filter1 = Noise
* TA::arg1 = add normal 5
*
* #add +/-10% of variation, uniformly distributed to HS
* HS::filter1 = Noise
* HS::arg1 = mult uniform 0.1
* @endcode
*/
class ProcNoise : public ProcessingBlock {
public:
ProcNoise(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);
private:
void uniform_add(const unsigned int& param, std::vector<MeteoData>& ovec) const;
void uniform_mult(const unsigned int& param, std::vector<MeteoData>& ovec) const;
void normal_add(const unsigned int& param, std::vector<MeteoData>& ovec) const;
void normal_mult(const unsigned int& param, std::vector<MeteoData>& ovec) const;
double getBoxMuller() const;
void parse_args(std::vector<std::string> vec_args);
double range;
char distribution, type;
};
} //end namespace
#endif
......@@ -43,6 +43,7 @@
#include <meteoio/meteoFilters/ProcUnshade.h>
#include <meteoio/meteoFilters/ProcAdd.h>
#include <meteoio/meteoFilters/ProcMult.h>
#include <meteoio/meteoFilters/ProcNoise.h>
#include <meteoio/meteoFilters/ProcExpSmoothing.h>
#include <meteoio/meteoFilters/ProcWMASmoothing.h>
......@@ -93,9 +94,10 @@ namespace mio {
* - UNHEATED_RAINGAUGE: detection of snow melting in a rain gauge, see FilterUnheatedHNW
*
* Some data transformations are also supported besides filtering, both very basic and generic data transformations:
* - SUPPR: delete all data, see FilterSuppr
* - SUPPR: delete data, see FilterSuppr
* - ADD: adds a given offset to the data, see ProcAdd
* - MULT: multiply the data by a given factor, see ProcMult
* - NOISE: add noise to the data, see ProcNoise
*
* As well as more specific data transformations:
* - EXP_SMOOTHING: exponential smoothing of data, see ProcExpSmoothing
......@@ -156,6 +158,8 @@ ProcessingBlock* BlockFactory::getBlock(const std::string& blockname, const std:
return new ProcMult(vec_args, blockname, root_path);
} else if (blockname == "ADD"){
return new ProcAdd(vec_args, blockname, root_path);
} else if (blockname == "NOISE"){
return new ProcNoise(vec_args, blockname);
} else if (blockname == "EXP_SMOOTHING"){
return new ProcExpSmoothing(vec_args, blockname);
} else if (blockname == "WMA_SMOOTHING"){
......
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