WSL/SLF GitLab Repository

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

A (more) intelligent handling of buffered grids has been implemented: the...

A (more) intelligent handling of buffered grids has been implemented: the grids are stored in the equivalent of a circular buffer of a given size (user defined, or 10 by default). This prevents running out of memory when processing large numbers of grids!

A method for computing quantiles has been implemented in libinterpol1D. A vector of data has to be provided as well as a vector of the desired quantiles (ex: 0.25, 0.5, 0.75) and a vector of values for these quantiles will be returned. So far, the algorithm that is used produces the same results as R in default operations. The possibility of choosing which quantiles algorithms might be added in the future.

The getJulianDayNumber method was just totally wrong... This has been fixed. One can also force a GMT day of year with a flag.

Finally, some documentation has been written/updated.
parent bc30c9c2
......@@ -39,6 +39,16 @@ BufferedIOHandler::~BufferedIOHandler() throw()
void BufferedIOHandler::bufferGrid(const Grid2DObject& in_grid2Dobj, const std::string& in_filename)
if(IndexBufferedGrids.size() >= max_grids) { //we need to remove the oldest grid
mapBufferedGrids.erase( mapBufferedGrids.find( IndexBufferedGrids[0] ) );
IndexBufferedGrids.erase( IndexBufferedGrids.begin() );
mapBufferedGrids[in_filename] = in_grid2Dobj;
IndexBufferedGrids.push_back( in_filename );
void BufferedIOHandler::read2DGrid(Grid2DObject& in_grid2Dobj, const std::string& in_filename)
std::map<std::string, Grid2DObject>::iterator it = mapBufferedGrids.find(in_filename);
......@@ -49,7 +59,7 @@ void BufferedIOHandler::read2DGrid(Grid2DObject& in_grid2Dobj, const std::string
Grid2DObject tmpgrid2D;
iohandler.read2DGrid(tmpgrid2D, in_filename);
mapBufferedGrids[in_filename] = tmpgrid2D;
bufferGrid(tmpgrid2D, in_filename);
in_grid2Dobj = tmpgrid2D;
......@@ -91,6 +101,7 @@ void BufferedIOHandler::readLanduse(Grid2DObject& in_grid2Dobj)
in_grid2Dobj = tmpgrid2D;
//HACK: manage buffering of assimilation grids! Why not considering them normal grids?
void BufferedIOHandler::readAssimilationData(const Date& in_date, Grid2DObject& in_grid2Dobj)
std::map<std::string, Grid2DObject>::iterator it = mapBufferedGrids.find("/:ASSIMILATIONDATA" + in_date.toString(Date::FULL));
......@@ -149,6 +160,9 @@ void BufferedIOHandler::setDfltBufferProperties()
buff_before = chunk_size * 0.1; //10% centering by default
max_grids = 10; //default number of grids to keep in buffer
cfg.getValue("BUFF_GRIDS", "General", max_grids, Config::nothrow);
double BufferedIOHandler::getAvgSamplingRate()
......@@ -42,6 +42,8 @@ namespace mio {
* - BUFF_BEFORE: alternate way of buffer centering: When rebuffering, the new date will be located BUFF_BEFORE days from the
* begining of the buffer (therefore, it takes a value in days); [General] section, optional. Only one of
* two centering option can be used.
* - BUFF_GRIDS: how many grids to keep in the buffer. If more grids have to be read, the oldest ones will be removed from
* the buffer. (10 by default)
* @author Thomas Egger
* @date 2009-07-25
......@@ -125,6 +127,7 @@ class BufferedIOHandler : public IOInterface {
void setDfltBufferProperties();
void bufferData(const Date& date_start, const Date& date_end, std::vector< METEO_TIMESERIE >& vecvecMeteo);
void bufferGrid(const Grid2DObject& in_grid2Dobj, const std::string& in_filename);
//private members
IOHandler& iohandler;
......@@ -132,11 +135,13 @@ class BufferedIOHandler : public IOInterface {
std::vector< METEO_TIMESERIE > vec_buffer_meteo;
std::map<std::string, Grid2DObject> mapBufferedGrids;
std::vector<std::string> IndexBufferedGrids;
Date buffer_start, buffer_end;
Duration chunk_size; ///< How much data to read at once
Duration buff_before; ///< How much data to read before the requested date in buffer
unsigned int chunks; ///< How many chuncks to buffer
unsigned int max_grids; ///< How many grids to buffer (grids, dems, landuse and assimilation grids together)
} //end namespace
......@@ -523,19 +523,26 @@ void Date::getDate(int& year_out, int& month_out, int& day_out, int& hour_out, i
* @brief Return year, month, day.
* @brief Return the julian day for the current date.
* Return the day of the year index for the current Date object
* @param gmt convert returned value to GMT? (default: false)
* @return julian day number
int Date::getJulianDayNumber() const {
int Date::getJulianDayNumber(const bool& gmt) const {
throw UnknownValueException("Date object is undefined!", AT);
//this is quite inefficient... we might want to deal with leap years with their rule + days arrays instead
int local_year, local_month, local_day, local_hour, local_minute;
getDate(local_year, local_month, local_day, local_hour, local_minute);
return (getJulianDayNumber(local_year, local_month, local_day));
if(gmt) {
const double first_day_of_year = getJulianDayNumber(gmt_year, 1, 1);
return (gmt_julian - first_day_of_year + 1);
} else {
const double local_julian = GMTToLocal(gmt_julian);
int local_year, local_month, local_day, local_hour, local_minute;
calculateValues(local_julian, local_year, local_month, local_day, local_hour, local_minute);
const double in_day_offset = 1./24.*((double)local_hour+1./60.*(double)local_minute) - 0.5;
const double first_day_of_year = static_cast<double>(getJulianDayNumber(local_year, 1, 1)) + in_day_offset;
return (local_julian - first_day_of_year + 1);
......@@ -125,7 +125,7 @@ class Date {
void getDate(int& year, int& month, int& day, int& hour, int& minute, const bool& gmt=false) const;
int getYear(const bool& gmt=false) const;
int getJulianDayNumber() const;
int getJulianDayNumber(const bool& gmt=false) const;
bool isLeapYear() const;
void rnd(const double& precision, const RND& type=CLOSEST);
......@@ -479,7 +479,7 @@ bool IOUtils::convertString(Date& t, const std::string& str, const double& time_
const std::string datum = str.substr(beg, end-beg);
const size_t d_len = datum.length();
if(d_len<8 || d_len>14) return false;
if(d_len<10 || d_len>14) return false;
if( convertString(year,datum.substr(0,4))==false ) return false;
if( convertString(month,datum.substr(4,2))==false ) return false;
if( convertString(day,datum.substr(6,2))==false ) return false;
......@@ -109,6 +109,14 @@ namespace IOUtils {
double angle_to_bearing(const double& angle);
* @brief Build a list of file in a given directory.
* The matching is very primitive: it only looks for the substring "pattern" in the file names.
* If this substrings exists, the file matches.
* @param path directory containing the files
* @param dirlist list of mathcing file names
* @param pattern optional pattern that must be part of the file names
void readDirectory(const std::string& path, std::list<std::string>& dirlist, const std::string& pattern = "");
bool validFileName(const std::string& filename);
......@@ -513,7 +513,7 @@ WARN_LOGFILE =
# directories like "/usr/src/myproject". Separate the files or directories
# with spaces.
INPUT = ./meteoio ./meteoio/plugins ./meteoio/meteolaws ./meteoio/meteofilters
INPUT = ./meteoio ./meteoio/plugins ./meteoio/meteolaws ./meteoio/meteofilters ./meteoio/meteostats
# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding, which is
......@@ -81,11 +81,24 @@ class Quadratic : public FitLeastSquare {
double f(const double& x);
//HACK: write copy constructor for Fit1D!!! needed to pass a Fit1D object by copy
* @class Fit1D
* @brief A class to perform 1D regressions
* It works on a time serie and uses matrix arithmetic to perform an arbitrary fit.
* It works on a time serie and uses either ad-hoc methods or matrix arithmetic to perform an arbitrary fit.
* Currently, the following models are supported:
* - SimpleLinear
* - NoisyLinear
* - SphericVario
* - LinVario
* - ExpVario
* - RatQuadVario
* - LinearLS
* - Quadratic
* However, the current way of specifying which model to use by providing the constructor with a string WILL change to
* an enum.
* @ingroup stats
* @author Mathias Bavay
......@@ -23,6 +23,62 @@ using namespace std;
namespace mio {
* @brief This function returns a vector of quantiles.
* The vector does not have to be sorted. See for more.
* This code is heavily inspired by Ken Wilder,
* (quantile method, replacing the nth-element call by direct access to a sorted vector).
* @param X vector to classify
* @param quartiles vector of quartiles, between 0 and 1
* @return vector of ordinates of the quantiles
std::vector<double> Interpol1D::quantiles(const std::vector<double>& X, const std::vector<double>& quartiles)
const size_t Xsize = X.size();
const size_t Qsize = quartiles.size();
if (Xsize == 0)
throw NoAvailableDataException("Trying to calculate quantiles with no data points", AT);
if (Qsize == 0)
throw NoAvailableDataException("No quantiles specified", AT);
//in order to properly escape nodata points, we need to copy in a temporary vector
vector<double> vecTemp;
for(size_t i=0; i<Xsize; i++) {
const double& value=X[i];
std::sort( vecTemp.begin(), vecTemp.end()); //since we will process several values, we sort the vector
//we will store results in a new vector
std::vector<double> vecResults(Qsize, IOUtils::nodata);
const size_t vecSize = vecTemp.size();
if (vecSize == 0) {
return vecResults; //ie: nodata values
if(vecSize == 1) {
std::fill(vecResults.begin(), vecResults.end(), vecTemp[0]);
return vecResults;
//compute quantiles
for(size_t ii=0; ii<Qsize; ii++) {
const double q = quartiles[ii];
if(q<=0.) vecResults[ii] = vecTemp[0];
else if(q>=1.) vecResults[ii] = vecTemp[vecSize-1];
else {
const double pos = static_cast<double>(vecSize - 1) * q;
const unsigned int ind = static_cast<unsigned int>(pos);
const double delta = pos - static_cast<double>(ind);
const double i1 = vecTemp[ind];
const double i2 = vecTemp[ind+1];
vecResults[ii] = i1 * (1.0 - delta) + i2 * delta;
return vecResults;
* @brief This function returns the vector of local derivatives, given a vector of abscissae and ordinates.
* The vectors must be sorted by ascending x. The derivatives will be centered if possible, left or right otherwise or nodata
......@@ -118,7 +174,7 @@ double Interpol1D::weightedMean(const double& d1, const double& d2, const double
* @brief This function returns the weighted aritmetic mean of two a vector.
* @brief This function returns the weighted aritmetic mean of a vector.
* See for more...
* @param vecData vector of values
* @param weight weights to apply to the mean
......@@ -37,6 +37,7 @@ namespace mio {
class Interpol1D {
static std::vector<double> quantiles(const std::vector<double>& X, const std::vector<double>& quartiles);
static std::vector<double> derivative(const std::vector<double>& X, const std::vector<double>& Y);
static void sort(std::vector<double>& X, std::vector<double>& Y);
static double weightedMean(const double& d1, const double& d2, const double& weight=1.);
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