WSL/SLF GitLab Repository

Commit 5b5cd758 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

Several 1D statistical methods has been implemented while the existing ones...

Several 1D statistical methods has been implemented while the existing ones have been made nodata safe.

Added a redefinition of the << operator for IOPlugin and IOHandler. This should be convenient for debugging!

Some small optimizations have been brought into the spatial interpolations following some tests with cachegrind. Nothing huge, but small improvements.

parent 9cc19d6d
......@@ -8,7 +8,8 @@ IF(PROJ4)
ADD_DEFINITIONS(-DPROJ4)
ENDIF(PROJ4)
SET(meteoio_sources
SET(meteoio_sources
IOPlugin.cc
MeteoProcessor.cc
ResamplingAlgorithms.cc
Meteo1DInterpolator.cc
......
......@@ -223,7 +223,7 @@ void Grid2DObject::setValues(const unsigned int& _ncols, const unsigned int& _nr
llcorner = _llcorner;
}
bool Grid2DObject::isSameGeolocalization(const Grid2DObject& target)
bool Grid2DObject::isSameGeolocalization(const Grid2DObject& target) const
{
if( ncols==target.ncols && nrows==target.nrows &&
llcorner==target.llcorner &&
......@@ -251,16 +251,13 @@ std::ostream& operator<<(std::ostream& os, const Grid2DObject& grid)
using namespace mio; //HACK for POPC
void Grid2DObject::Serialize(POPBuffer &buf, bool pack)
{
if (pack)
{
if (pack) {
buf.Pack(&ncols,1);
buf.Pack(&nrows,1);
buf.Pack(&cellsize,1);
marshal_Coords(buf, llcorner, 0, FLAG_MARSHAL, NULL);
marshal_DOUBLE2D(buf, grid2D, 0, FLAG_MARSHAL, NULL);
}
else
{
} else {
buf.UnPack(&ncols,1);
buf.UnPack(&nrows,1);
buf.UnPack(&cellsize,1);
......
......@@ -133,7 +133,7 @@ class Grid2DObject {
* @param target (Grid2DObject) grid to compare to
* @return (bool) true if same geolocalization
*/
bool isSameGeolocalization(const Grid2DObject& target);
bool isSameGeolocalization(const Grid2DObject& target) const;
Array2D<double> grid2D; ///<the grid itself (simple 2D table containing the values for each point)
unsigned int ncols; ///<number of columns in the grid
......
......@@ -90,7 +90,7 @@ IOHandler::IOHandler(const std::string& configfile) : IOInterface(NULL), cfg(con
//Copy constructor
#ifdef _POPC_
//IOHandler::IOHandler(const IOHandler& aio) : cfg(aio.cfg), fileio(aio.cfg), bormaio(aio.cfg), imisio(aio.cfg){
//Nothing else so far
//Nothing else so far //HACK for POPC
//}
#else
IOHandler::IOHandler(const IOHandler& aio) : IOInterface(NULL), cfg(aio.cfg), fileio(aio.cfg)
......@@ -196,10 +196,10 @@ IOInterface* IOHandler::getPlugin(const std::string& cfgkey, const std::string&
return (mapit->second).io;
}
void IOHandler::read2DGrid(Grid2DObject& _grid, const std::string& _filename)
void IOHandler::read2DGrid(Grid2DObject& out_grid, const std::string& _filename)
{
IOInterface *plugin = getPlugin("GRID2D", "Input");
plugin->read2DGrid(_grid, _filename);
plugin->read2DGrid(out_grid, _filename);
}
void IOHandler::readDEM(DEMObject& dem_out)
......@@ -285,4 +285,21 @@ void IOHandler::write2DGrid(const Grid2DObject& grid_in, const std::string& name
plugin->write2DGrid(grid_in, name);
}
std::ostream& operator<<(std::ostream& os, const IOHandler& data)
{
os << "<IOHandler>\n";
os << data.cfg;
os << "<mapPlugins>\n";
os << setw(10) << "Keyword" << " = " << IOPlugin::header << "\n";
std::map<std::string, IOPlugin::IOPlugin>::const_iterator it1;
for (it1=data.mapPlugins.begin(); it1 != data.mapPlugins.end(); it1++){
os << setw(10) << it1->first << " = " << it1->second;
}
os << "</mapPlugins>\n";
os << "</IOHandler>\n";
return os;
}
} //end namespace
......@@ -47,7 +47,7 @@ class IOHandler : public IOInterface {
~IOHandler() throw();
//methods defined in the IOInterface class
virtual void read2DGrid(Grid2DObject& dem_out, const std::string& parameter="");
virtual void read2DGrid(Grid2DObject& out_grid, const std::string& parameter="");
virtual void readDEM(DEMObject& dem_out);
virtual void readLanduse(Grid2DObject& landuse_out);
virtual void readStationData(const Date& date,
......@@ -65,6 +65,8 @@ class IOHandler : public IOInterface {
virtual void readSpecialPoints(std::vector<Coords>& pts);
virtual void write2DGrid(const Grid2DObject& grid_in, const std::string& name);
friend std::ostream& operator<<(std::ostream& os, const IOHandler& data);
private:
void loadDynamicPlugins();
void loadPlugin(const std::string& libname, const std::string& classname,
......
......@@ -49,7 +49,7 @@ parclass IOHandler {
~IOHandler();
//methods defined in the IOInterface class
virtual void read2DGrid([out]Grid2DObject& dem_out, const std::string& parameter="");
virtual void read2DGrid([out]Grid2DObject& out_grid, const std::string& parameter="");
virtual void readDEM([out]DEMObject& dem_out);
virtual void readLanduse([out]Grid2DObject& landuse_out);
virtual void readStationData([in]const Date& date,
......@@ -67,6 +67,8 @@ parclass IOHandler {
virtual void readSpecialPoints([out,proc=marshal_vec_coords]std::vector<Coords>& pts);
virtual void write2DGrid([in]const Grid2DObject& grid_in, [in]const std::string& name);
friend std::ostream& operator<<(std::ostream& os, const IOHandler& data);
private:
void loadDynamicPlugins();
void loadPlugin(const std::string& libname, const std::string& classname,
......
/***********************************************************************************/
/* Copyright 2010 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 <meteoio/IOPlugin.h>
namespace mio {
const std::string IOPlugin::header="<IOPlugin> libname, classname, &class, &lib</IOPlugin>";
std::ostream& operator<<(std::ostream& os, const IOPlugin& data) {
const unsigned int pt_w=8;
os << "<IOPlugin>" << std::setw(21) << data.libname << "," << std::setw(10) << data.classname;
os << "," << std::showbase << std::setw(pt_w+2) << data.io;
os << "," << std::showbase << std::setw(pt_w+2) << data.dynLibrary << "</IOPlugin>\n";
return os;
}
} //end namespace
......@@ -51,7 +51,12 @@ class IOPlugin {
*/
IOPlugin(std::string _s1, std::string _s2, IOInterface *p1, DynamicLibrary *p2) : libname(_s1), classname(_s2), io(p1), dynLibrary(p2){}
IOPlugin() : libname(""), classname(""), io(NULL), dynLibrary(NULL){}
friend std::ostream& operator<<(std::ostream& os, const IOPlugin& data);
static const std::string header; //to contain a helpful header for understanding the output of <<
};
} //end namespace
#endif
......@@ -74,7 +74,9 @@ InterpolationAlgorithm* AlgorithmFactory::getAlgorithm(const std::string& _algon
return new RHAlgorithm(_mi, _dem, _vecMeteo, _vecStation, _vecArgs, _algoname);
} else if (algoname == "WIND_CURV"){
return new SimpleWindInterpolationAlgorithm(_mi, _dem, _vecMeteo, _vecStation, _vecArgs, _algoname);
} else if (algoname == "USER"){
/*} else if (algoname == "KRIG"){
return new SimpleKrigingAlgorithm(_mi, _dem, _vecMeteo, _vecStation, _vecArgs, _algoname);
*/} else if (algoname == "USER"){
return new USERinterpolation(_mi, _dem, _vecMeteo, _vecStation, _vecArgs, _algoname);
} else {
throw IOException("The interpolation algorithm '"+algoname+"' is not implemented" , AT);
......
......@@ -149,10 +149,11 @@ void Meteo2DInterpolator::checkMinMax(const double& minval, const double& maxval
for (unsigned int ii=0; ii<nx; ii++){
double& value = gridobj.grid2D(ii,jj);
if (value == IOUtils::nodata){
//Do nothing
} else if (value < minval){
continue;
}
if (value < minval) {
value = minval;
} else if (value > maxval){
} else if (value > maxval) {
value = maxval;
}
}
......
......@@ -57,51 +57,68 @@ double Interpol1D::linearInterpolation(const double& d1, const double& d2, const
double Interpol1D::arithmeticMean(const std::vector<double>& vecData)
{//This method is not safe against vector containing nodata
{
if (vecData.size() == 0)
throw NoAvailableDataException("Trying to calculate an arithmetic mean with no data points", AT);
unsigned int count=0;
double sum = 0.0;
for (unsigned int ii=0; ii<vecData.size(); ii++){
sum += vecData[ii];
const double value=vecData[ii];
if(value!=IOUtils::nodata) {
sum += vecData[ii];
count++;
}
}
return (sum/(double)vecData.size());
if(count>0)
return (sum/(double)count);
else
return IOUtils::nodata;
}
double Interpol1D::getMedian(const std::vector<double>& vecData)
{//This method is not safe against vector containing nodata
{
if (vecData.size() == 0)
throw NoAvailableDataException("Trying to calculate a median with no data points", AT);
vector<double> vecTemp;
for(unsigned int i=0; i<vecData.size(); i++) {
const double& value=vecData[i];
if(value!=IOUtils::nodata)
vecTemp.push_back(value);
}
vector<double> vecTemp(vecData); //copy by value of vecData
if (vecTemp.size() == 0)
return IOUtils::nodata;
double median = 0.0;
unsigned int vecSize = vecTemp.size();
unsigned int middle = (unsigned int)(vecSize/2);
const unsigned int vecSize = vecTemp.size();
const unsigned int middle = (unsigned int)(vecSize/2);
sort(vecTemp.begin(), vecTemp.end());
if ((vecSize % 2) == 1){ //uneven
median = vecTemp.at(middle);
return vecTemp.at(middle);
} else { //use arithmetic mean of element n/2 and n/2-1
median = Interpol1D::linearInterpolation(vecTemp.at(middle-1), vecTemp.at(middle), 0.5);
return Interpol1D::linearInterpolation(vecTemp.at(middle-1), vecTemp.at(middle), 0.5);
}
return median;
}
double Interpol1D::getMedianAverageDeviation(const std::vector<double>& vecData)
{//This method is not safe against vector containing nodata
{
if (vecData.size() == 0)
throw NoAvailableDataException("Trying to calculate MAD with no data points", AT);
vector<double> vecWindow(vecData);
double median = Interpol1D::getMedian(vecData);
double median = Interpol1D::getMedian(vecWindow);
if(median==IOUtils::nodata)
return IOUtils::nodata;
//Calculate vector of deviations and write each value back into the vecWindow
for(unsigned int ii=0; ii<vecWindow.size(); ii++){
vecWindow[ii] = std::abs(vecWindow[ii] - median);
double& value=vecWindow[ii];
if(value!=IOUtils::nodata)
value = std::abs(value - median);
}
//Calculate the median of the deviations
......@@ -110,4 +127,119 @@ double Interpol1D::getMedianAverageDeviation(const std::vector<double>& vecData)
return mad;
}
double Interpol1D::variance(const std::vector<double>& X)
{
const unsigned int n = X.size();
unsigned int count=0;
double sum=0., sum_sq=0.;
for(unsigned int i=0; i<n; i++) {
const double value=X[i];
if(value!=IOUtils::nodata) {
sum += value;
sum_sq += value*value;
count++;
}
}
if(count<=1) return IOUtils::nodata;
const double mean = sum/(double)count;
return (sum_sq - sum*mean)/((double)count-1.);
}
double Interpol1D::std_dev(const std::vector<double>& X)
{
return sqrt(variance(X));
}
double Interpol1D::covariance(const std::vector<double>& X, const std::vector<double>& Y)
{
if(X.size()!=Y.size())
throw IOException("Vectors should have the same size for covariance!", AT);
const unsigned int n = X.size();
if(n==0) return IOUtils::nodata;
const double X_mean = Interpol1D::arithmeticMean(X);
const double Y_mean = Interpol1D::arithmeticMean(Y);
if(X_mean==IOUtils::nodata || Y_mean==IOUtils::nodata)
return IOUtils::nodata;
unsigned int count=0;
double sum=0.;
for(unsigned int i=0; i<n; i++) {
if(X[i]!=IOUtils::nodata && Y[i]!=IOUtils::nodata) {
sum += (X[i] - X_mean) * (Y[i] - Y_mean);
count++;
}
}
if(count<=1) return IOUtils::nodata;
return sum/((double)count-1.);
}
/**
* @brief Computes the linear regression coefficients fitting the points given as X and Y in two vectors
* the linear regression has the form Y = aX + b with a regression coefficient r (it is nodata safe)
* @param X vector of X coordinates
* @param Y vector of Y coordinates (same order as X)
* @param a slope of the linear regression
* @param b origin of the linear regression
* @param r linear regression coefficient
*/
void Interpol1D::LinRegression(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r)
{ //check arguments
const unsigned int n=X.size();
if(n==0)
throw NoAvailableDataException("Trying to calculate linear regression with no data points", AT);
if(n!=Y.size())
throw IOException("Vectors should have the same size for linear regression!", AT);
//computing x_avg and y_avg
double x_avg=0., y_avg=0.;
unsigned int count=0;
for (unsigned int i=0; i<n; i++) {
if(X[i]!=IOUtils::nodata && Y[i]!=IOUtils::nodata) {
x_avg += X[i];
y_avg += Y[i];
count++;
}
}
if(count==0)
throw NoAvailableDataException("Trying to calculate linear regression with no valid data points", AT);
x_avg /= (double)count;
y_avg /= (double)count;
//computing sx, sy, sxy
double sx=0., sy=0., sxy=0.;
for (unsigned int i=0; i<n; i++) {
if(X[i]!=IOUtils::nodata && Y[i]!=IOUtils::nodata) {
sx += (X[i]-x_avg) * (X[i]-x_avg);
sy += (Y[i]-y_avg) * (Y[i]-y_avg);
sxy += (X[i]-x_avg) * (Y[i]-y_avg);
}
}
//computing the regression line
const double epsilon = 1e-6;
if(sx <= abs(x_avg)*epsilon) { //sx and sy are always positive
//all points have same X -> we return a constant value that is the average
a = 0.;
b = y_avg;
r = 1.;
std::cerr << "[W] Computing linear regression on data at identical X\n";
return;
}
a = sxy / sx;
b = y_avg - a*x_avg;
if(sy==0) {
//horizontal line: all y's are equals
r = 1.;
} else {
//any other line
r = sxy / sqrt(sx*sy);
}
}
//TODO: with variable changes, compute log and exp regression
} //namespace
......@@ -19,6 +19,7 @@
#define __LIBINTERPOL1D_H__
#include <meteoio/IOExceptions.h>
#include <meteoio/IOUtils.h>
#include <cmath>
#include <vector>
......@@ -34,6 +35,11 @@ class Interpol1D {
static double arithmeticMean(const std::vector<double>& vecData);
static double getMedian(const std::vector<double>& vecData);
static double getMedianAverageDeviation(const std::vector<double>& vecData);
static double variance(const std::vector<double>& X);
static double std_dev(const std::vector<double>& X);
static double covariance(const std::vector<double>& z1, const std::vector<double>& z2);
static void LinRegression(const std::vector<double>& X, const std::vector<double>& Y, double& a, double& b, double& r);
};
} //end namespace
......
......@@ -125,7 +125,7 @@ double Interpol2D::weightInvDist(const double& d2)
}
double Interpol2D::weightInvDistSqrt(const double& d2)
{
return invSqrt( invSqrt(d2) ); //we use the optimized approximation for 1/sqrt
return fastSqrt_Q3( invSqrt(d2) ); //we use the optimized approximation for 1/sqrt
}
double Interpol2D::weightInvDist2(const double& d2)
{
......@@ -137,56 +137,6 @@ double Interpol2D::weightInvDistN(const double& d2)
}
//Data regression models
/**
* @brief Computes the linear regression coefficients fitting the points given as X and Y in two vectors
* the linear regression has the form Y = aX + b with a regression coefficient r
* @param X vector of X coordinates
* @param Y vector of Y coordinates (same order as X)
* @param a slope of the linear regression
* @param b origin of the linear regression
* @param r linear regression coefficient
*/
void Interpol2D::LinRegressionCore(const std::vector<double>& X, const std::vector<double>& Y, const unsigned int imax, double& a, double& b, double& r)
{//finds the linear regression for points (x,y,z,Value)
double x_avg=0., y_avg=0.;
double sx=0., sy=0., sxy=0.;
//computing x_avg and y_avg
for (unsigned int i=0; i<imax; i++) {
x_avg += X[i];
y_avg += Y[i];
}
x_avg /= (double)imax;
y_avg /= (double)imax;
//computing sx, sy, sxy
for (unsigned int i=0; i<imax; i++) {
sx += (X[i]-x_avg) * (X[i]-x_avg);
sy += (Y[i]-y_avg) * (Y[i]-y_avg);
sxy += (X[i]-x_avg) * (Y[i]-y_avg);
}
//computing the regression line
const double epsilon = 1e-6;
if(sx <= x_avg*epsilon) {
//all points have same X -> we return a constant value that is the average
a = 0.;
b = y_avg;
r = 1.;
std::cerr << "[W] Computing linear regression on data at identical X\n";
return;
}
a = sxy / sx;
b = y_avg - a*x_avg;
if(sy==0) {
//horizontal line: all y's are equals
r = 1.;
} else {
//any other line
r = sxy / sqrt(sx*sy);
}
}
/**
* @brief Computes the linear regression coefficients fitting the points given as X and Y in two vectors
* the linear regression has the form Y = aX + b with a regression coefficient r. If the regression coefficient is not good enough, a bad point is looked removed.
......@@ -200,9 +150,9 @@ int Interpol2D::LinRegression(const std::vector<double>& in_X, const std::vector
//finds the linear regression for points (x,y,z,Value)
const double r_thres=0.7;
//we want at least 4 points AND 85% of the initial data set kept in the regression
const double min_dataset=floor(0.85*(double)in_X.size());
const unsigned int min_dataset=(unsigned int)floor(0.85*(double)in_X.size());
const unsigned int min_pts=(min_dataset>4)?min_dataset:4;
unsigned int nb_pts = in_X.size();
const unsigned int nb_pts = in_X.size();
double a,b,r;
if (nb_pts==2) {
......@@ -216,23 +166,21 @@ int Interpol2D::LinRegression(const std::vector<double>& in_X, const std::vector
return EXIT_FAILURE;
}
LinRegressionCore(in_X, in_Y, nb_pts, coeffs[1], coeffs[2], coeffs[3]);
Interpol1D::LinRegression(in_X, in_Y, coeffs[1], coeffs[2], coeffs[3]);
if(fabs(coeffs[3])>=r_thres)
return EXIT_SUCCESS;
std::vector<double> X(in_X), Y(in_Y);
unsigned int nb_valid_pts=nb_pts;
while(fabs(coeffs[3])<r_thres && nb_pts>min_pts) {
while(fabs(coeffs[3])<r_thres && nb_valid_pts>min_pts) {
//we try to remove the one point in the data set that is the worst
coeffs[3]=0.;
unsigned int index_bad=0;
for (unsigned int i=0; i<nb_pts; i++) {
//removing alternatively each point
//index nb_pts is used as temporary storage
const double X_tmp=X[i]; X[i]=X[nb_pts-1];
const double Y_tmp=Y[i]; Y[i]=Y[nb_pts-1];
LinRegressionCore(X, Y, nb_pts-1, a, b, r);
X[i]=X_tmp;
//invalidating alternatively each point
const double Y_tmp=Y[i]; Y[i]=IOUtils::nodata;
Interpol1D::LinRegression(X, Y, a, b, r);
Y[i]=Y_tmp;
if (fabs(r)>fabs(coeffs[3])) {
......@@ -243,9 +191,8 @@ int Interpol2D::LinRegression(const std::vector<double>& in_X, const std::vector
}
}
//the worst point has been found, we overwrite it
X[index_bad]=X[nb_pts-1];
Y[index_bad]=Y[nb_pts-1];
nb_pts--;
Y[index_bad]=IOUtils::nodata;
nb_valid_pts--;
}
//check if r is reasonnable
......@@ -301,12 +248,12 @@ double Interpol2D::LinProject(const double& value, const double& altitude, const
* @param grid 2D array to fill
*/
void Interpol2D::stdPressureGrid2DFill(const DEMObject& dem, Grid2DObject& grid) {
grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);
grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);
//provide each point with an altitude dependant pressure... it is worth what it is...
for (unsigned int j=0; j<grid.nrows; j++) {
for (unsigned int i=0; i<grid.ncols; i++) {
const double cell_altitude=dem.grid2D(i,j);
const double& cell_altitude=dem.grid2D(i,j);
if (cell_altitude!=IOUtils::nodata) {
grid.grid2D(i,j) = lw_AirPressure(cell_altitude);
} else {
......@@ -374,14 +321,18 @@ double Interpol2D::IDWCore(const double& x, const double& y, const std::vector<d
const std::vector<StationData>& vecStations_in)
{
//The value at any given cell is the sum of the weighted contribution from each source
const unsigned int n_stations=vecStations_in.size();
double parameter=0., norm=0.;
const double scale = 1.e6;
for (unsigned int i=0; i<(unsigned int)vecStations_in.size(); i++) {
for (unsigned int i=0; i<n_stations; i++) {
/*const double weight=1./(HorizontalDistance(x, y, vecStations_in[i].position.getEasting(),
vecStations_in[i].position.getNorthing()) + 1e-6);*/
const double DX=x-vecStations_in[i].position.getEasting();
const double DY=y-vecStations_in[i].position.getNorthing();
const double weight = invSqrt( DX*DX + DY*DY ); //use the optimized 1/sqrt approximation
const Coords& position = vecStations_in[i].position;
const double DX=x-position.getEasting();
const double DY=y-position.getNorthing();
const double weight = invSqrt( DX*DX + DY*DY + scale ); //use the optimized 1/sqrt approximation
//const double weight = weightInvDistSqrt( (DX*DX + DY*DY) );
parameter += weight*vecData_in[i];
norm += weight;
}
......@@ -407,11 +358,12 @@ void Interpol2D::LapseIDW(const std::vector<double>& vecData_in, const std::vect
Grid2DObject& grid)
{ //multiple source stations: lapse rate projection, IDW Krieging, re-projection
const double ref_altitude = getReferenceAltitude(dem);
const unsigned int n_stations=vecStations_in.size();
grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);
std::vector<double> vecTref(vecStations_in.size(), 0.0); // init to 0.0
for (unsigned int i=0; i<(unsigned int)vecStations_in.size(); i++) {
for (unsigned int i=0; i<n_stations; i++) {
vecTref[i] = funcptr(vecData_in[i], vecStations_in[i].position.getAltitude(),
ref_altitude, vecCoefficients);
}
......