WSL/SLF GitLab Repository

Commit 84602c43 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

Following the latest changes, this is a general code cleanup in all the...

Following the latest changes, this is a general code cleanup in all the classes involved in spatial interpolations. The LIDW method has been commented out since this works very poorly and would need to be properly redone.
parent a7d259a3
......@@ -32,7 +32,6 @@ InterpolationAlgorithm* AlgorithmFactory::getAlgorithm(const std::string& i_algo
const std::vector<std::string>& i_vecArgs, IOManager& iom)
{
const std::string algoname( IOUtils::strToUpper(i_algoname) );
//IOUtils::toUpper(algoname);
if (algoname == "CST"){// constant fill
return new ConstAlgorithm(i_mi, i_date, i_dem, i_vecArgs, i_algoname, iom);
......@@ -44,9 +43,9 @@ InterpolationAlgorithm* AlgorithmFactory::getAlgorithm(const std::string& i_algo
return new IDWAlgorithm(i_mi, i_date, i_dem, i_vecArgs, i_algoname, iom);
} else if (algoname == "IDW_LAPSE"){// Inverse Distance Weighting with an elevation lapse rate fill
return new IDWLapseAlgorithm(i_mi, i_date, i_dem, i_vecArgs, i_algoname, iom);
} else if (algoname == "LIDW_LAPSE"){// Inverse Distance Weighting with an elevation lapse rate fill, restricted to a local scale
} /*else if (algoname == "LIDW_LAPSE"){// Inverse Distance Weighting with an elevation lapse rate fill, restricted to a local scale
return new LocalIDWLapseAlgorithm(i_mi, i_date, i_dem, i_vecArgs, i_algoname, iom);
} else if (algoname == "RH"){// relative humidity interpolation
}*/ else if (algoname == "RH"){// relative humidity interpolation
return new RHAlgorithm(i_mi, i_date, i_dem, i_vecArgs, i_algoname, iom);
} else if (algoname == "ILWR"){// long wave radiation interpolation
return new ILWRAlgorithm(i_mi, i_date, i_dem, i_vecArgs, i_algoname, iom);
......@@ -78,8 +77,15 @@ InterpolationAlgorithm::InterpolationAlgorithm(Meteo2DInterpolator& i_mi,
size_t InterpolationAlgorithm::getData(const MeteoData::Parameters& i_param,
std::vector<double>& o_vecData) const
{
vector<StationData> tmp_vecMeta;
return getData(i_param, o_vecData, tmp_vecMeta);
o_vecData.clear();
for (size_t ii=0; ii<vecMeteo.size(); ii++){
const double& val = vecMeteo[ii](i_param);
if (val != IOUtils::nodata) {
o_vecData.push_back(val);
}
}
return o_vecData.size();
}
size_t InterpolationAlgorithm::getData(const MeteoData::Parameters& i_param,
......@@ -87,7 +93,6 @@ size_t InterpolationAlgorithm::getData(const MeteoData::Parameters& i_param,
{
o_vecData.clear();
o_vecMeta.clear();
for (size_t ii=0; ii<vecMeteo.size(); ii++){
const double& val = vecMeteo[ii](i_param);
if (val != IOUtils::nodata){
......@@ -104,9 +109,10 @@ size_t InterpolationAlgorithm::getStationAltitudes(const std::vector<StationData
{
o_vecData.clear();
for (size_t ii=0; ii<i_vecMeta.size(); ii++){
const double& val = i_vecMeta[ii].position.getAltitude();
if (val != IOUtils::nodata)
o_vecData.push_back(val);
const double& alt = i_vecMeta[ii].position.getAltitude();
if (alt != IOUtils::nodata) {
o_vecData.push_back(alt);
}
}
return o_vecData.size();
......@@ -124,7 +130,7 @@ std::string InterpolationAlgorithm::getInfo() const
os << " station";
else
os << " stations";
std::string tmp( info.str() );
const std::string tmp( info.str() );
if(!tmp.empty()) {
os << ", " << tmp;
}
......@@ -235,7 +241,7 @@ double StandardPressureAlgorithm::getQualityRating() const
}
void StandardPressureAlgorithm::calculate(Grid2DObject& grid) {
Interpol2D::stdPressureGrid2DFill(dem, grid);
Interpol2D::stdPressure(dem, grid);
}
......@@ -262,7 +268,7 @@ void ConstAlgorithm::calculate(Grid2DObject& grid) {
if (nrOfMeasurments == 0)
throw IOException("Interpolation FAILED for parameter " + MeteoData::getParameterName(param), AT);
Interpol2D::constantGrid2DFill(Interpol1D::arithmeticMean(vecData), dem, grid);
Interpol2D::constant(Interpol1D::arithmeticMean(vecData), dem, grid);
}
......@@ -301,8 +307,7 @@ void ConstLapseRateAlgorithm::calculate(Grid2DObject& grid)
getTrend(vecAltitudes, vecData, trend);
info << trend.getInfo();
detrend(trend, vecAltitudes, vecData);
const double avgData = Interpol1D::arithmeticMean(vecData);
Interpol2D::constantGrid2DFill(avgData, dem, grid);
Interpol2D::constant(Interpol1D::arithmeticMean(vecData), dem, grid);
retrend(trend, grid);
}
......@@ -376,9 +381,9 @@ double LocalIDWLapseAlgorithm::getQualityRating() const
return 0.7;
}
void LocalIDWLapseAlgorithm::calculate(Grid2DObject& grid)
void LocalIDWLapseAlgorithm::calculate(Grid2DObject& /*grid*/)
{
if (nrOfMeasurments == 0)
/*if (nrOfMeasurments == 0)
throw IOException("Interpolation FAILED for parameter " + MeteoData::getParameterName(param), AT);
//Set regression coefficients
......@@ -394,7 +399,7 @@ void LocalIDWLapseAlgorithm::calculate(Grid2DObject& grid)
double r2=0.;
Interpol2D::LocalLapseIDW(vecData, vecMeta, dem, nrOfNeighbors, grid, r2); //HACK
info << "r^2=" << Optim::pow2( r2 );
info << "r^2=" << Optim::pow2( r2 );*/
}
......@@ -544,9 +549,7 @@ double SimpleWindInterpolationAlgorithm::getQualityRating() const
return 0.0;
//This algorithm requires the curvatures
unsigned int nx=0, ny=0;
dem.curvature.size(nx,ny);
if (nx==0 || ny==0) {
if (dem.curvature.isEmpty()) {
std::cerr << "[W] WIND_CURV spatial interpolations algorithm selected but no dem curvature available! Skipping algorithm...\n";
return 0.0;
}
......@@ -716,7 +719,7 @@ void OrdinaryKrigingAlgorithm::calculate(Grid2DObject& grid)
retrend(trend, grid);
}
double OrdinaryKrigingAlgorithm::computeVariogram()
bool OrdinaryKrigingAlgorithm::computeVariogram()
{//return variogram fit of covariance between stations i and j
std::string vario_model;
if (vecArgs.empty()){
......@@ -758,7 +761,7 @@ double OrdinaryKrigingAlgorithm::computeVariogram()
}
throw IOException("The variogram model "+variogram.getName()+" could not be fitted to the available data!", AT);
}
return 1.;
return true;
}
} //namespace
......
......@@ -77,13 +77,13 @@ class Meteo2DInterpolator; // forward declaration, cyclic header include
* - CST_LAPSE: constant value reprojected to the elevation of the cell (see ConstLapseRateAlgorithm)
* - IDW: Inverse Distance Weighting averaging (see IDWAlgorithm)
* - IDW_LAPSE: Inverse Distance Weighting averaging with reprojection to the elevation of the cell (see IDWLapseAlgorithm)
* - LIDW_LAPSE: IDW_LAPSE restricted to a local scale (n neighbor stations, see LocalIDWLapseAlgorithm)
* - RH: the dew point temperatures are interpolated using IDW_LAPSE, then reconverted locally to relative humidity (see RHAlgorithm)
* - ILWR: the incoming long wave radiation is converted to emissivity and then interpolated (see ILWRAlgorithm)
* - WIND_CURV: the wind field (VW and DW) is interpolated using IDW_LAPSE and then altered depending on the local curvature and slope (taken from the DEM, see SimpleWindInterpolationAlgorithm)
* - HNW_SNOW: precipitation interpolation according to (Magnusson, 2011) (see SnowHNWInterpolation)
* - ODKRIG: ordinary kriging (see OrdinaryKrigingAlgorithm)
* - USER: user provided grids to be read from disk (if available, see USERInterpolation)
* <!-- - LIDW_LAPSE: IDW_LAPSE restricted to a local scale (n neighbor stations, see LocalIDWLapseAlgorithm) -->
*
* @section interpol2D_lapse Lapse rates
* Several algorithms use elevation trends, currently modelled as a linear relation. The slope of this linear relation can
......@@ -479,7 +479,7 @@ class OrdinaryKrigingAlgorithm : public InterpolationAlgorithm {
virtual double getQualityRating() const;
virtual void calculate(Grid2DObject& grid);
private:
double computeVariogram();
bool computeVariogram();
Fit1D variogram;
};
......
......@@ -113,7 +113,7 @@ bool SimpleLinear::checkInputs()
return true;
}
double SimpleLinear::f(const double& x) {
double SimpleLinear::f(const double& x) const {
return Lambda.at(0)*x + Lambda.at(1);
}
......@@ -168,7 +168,7 @@ bool NoisyLinear::fit()
}
//regression models using the standard least square algorithm
double SphericVario::f(const double& x) {
double SphericVario::f(const double& x) const {
//c0>=0, cs>=0, as>=0
const double c0 = Lambda.at(0);
const double cs = Lambda.at(1);
......@@ -192,7 +192,7 @@ void SphericVario::setDefaultGuess() {
Lambda.push_back( *max_element(X.begin(), X.end()) );
}
double LinVario::f(const double& x) {
double LinVario::f(const double& x) const {
//c0>=0, b1>=0
const double c0 = Lambda.at(0);
const double bl = Lambda.at(1);
......@@ -216,7 +216,7 @@ void LinVario::setDefaultGuess() {
Lambda.push_back( slope );
}
double ExpVario::f(const double& x) {
double ExpVario::f(const double& x) const {
//c0>=0, ce>=0, ae>=0
const double c0 = Lambda.at(0);
const double ce = Lambda.at(1);
......@@ -241,7 +241,7 @@ void ExpVario::setDefaultGuess() {
Lambda.push_back( 1. );
}
double RatQuadVario::f(const double& x) {
double RatQuadVario::f(const double& x) const {
//c0>=0, cr>=0, ar>=0
const double c0 = Lambda.at(0);
const double cr = Lambda.at(1);
......@@ -266,7 +266,7 @@ void RatQuadVario::setDefaultGuess() {
Lambda.push_back( 1. );
}
double LinearLS::f(const double& x) {
double LinearLS::f(const double& x) const {
const double y = Lambda.at(0)*x + Lambda.at(1); //Lambda is a vector
return y;
}
......@@ -283,7 +283,7 @@ void LinearLS::setDefaultGuess() {
Lambda.push_back( Y[xzero_idx] );
}
double Quadratic::f(const double& x) {
double Quadratic::f(const double& x) const {
const double y = Lambda.at(0)*x*x + Lambda.at(1)*x + Lambda.at(2); //Lambda is a vector
return y;
}
......
......@@ -30,7 +30,7 @@ class Zero : public FitModel {
Zero() {fit_ready = true; nParam = 0; min_nb_pts = 0; regname = "Zero";};
void setData(const std::vector<double>& /*in_X*/, const std::vector<double>& /*in_Y*/) { };
bool fit() { return true;};
double f(const double& /*x*/) {return 0.;};
double f(const double& /*x*/) const {return 0.;};
};
class SimpleLinear : public FitModel {
......@@ -38,7 +38,7 @@ class SimpleLinear : public FitModel {
SimpleLinear() : fixed_lapse_rate(IOUtils::nodata) {fit_ready = false; nParam = 2; min_nb_pts = 3; regname = "SimpleLinear";};
void setData(const std::vector<double>& in_X, const std::vector<double>& in_Y);
bool fit();
double f(const double& x);
double f(const double& x) const;
void setLapseRate(const double& in_lapse_rate) {fixed_lapse_rate = in_lapse_rate; fit_ready = false; min_nb_pts=1;};
protected:
bool checkInputs();
......@@ -55,42 +55,42 @@ class SphericVario : public FitLeastSquare {
public:
SphericVario() {fit_ready = false; nParam = 3; min_nb_pts = 4; regname = "SphericVario";};
void setDefaultGuess();
double f(const double& x);
double f(const double& x) const;
};
class LinVario : public FitLeastSquare {
public:
LinVario() {fit_ready = false; nParam = 2; min_nb_pts = 3; regname = "LinVario";};
void setDefaultGuess();
double f(const double& x);
double f(const double& x) const;
};
class ExpVario : public FitLeastSquare {
public:
ExpVario() {fit_ready = false; nParam = 3; min_nb_pts = 4; regname = "ExpVario";};
void setDefaultGuess();
double f(const double& x);
double f(const double& x) const;
};
class RatQuadVario : public FitLeastSquare {
public:
RatQuadVario() {fit_ready = false; nParam = 3; min_nb_pts = 4; regname = "RatQuadVario";};
void setDefaultGuess();
double f(const double& x);
double f(const double& x) const;
};
class LinearLS : public FitLeastSquare {
public:
LinearLS() {fit_ready = false; nParam = 2; min_nb_pts = 3; regname = "LinearLS";};
void setDefaultGuess();
double f(const double& x);
double f(const double& x) const;
};
class Quadratic : public FitLeastSquare {
public:
Quadratic() {fit_ready = false; nParam = 2; min_nb_pts = 3; regname = "Quadratic";};
void setDefaultGuess();
double f(const double& x);
double f(const double& x) const;
};
/**
......
......@@ -24,14 +24,14 @@ using namespace std;
namespace mio {
void FitModel::getParams(std::vector<double>& o_coefficients) {
void FitModel::getParams(std::vector<double>& o_coefficients) const {
if(fit_ready!=true) {
throw InvalidArgumentException("The regression has not yet being computed!", AT);
}
o_coefficients = Lambda;
}
std::string FitModel::getInfo() {
std::string FitModel::getInfo() const {
return infoString;
}
......
......@@ -35,10 +35,10 @@ class FitModel {
void setGuess(const std::vector<double>& lambda_in);
virtual void setLapseRate(const double& /*lapse_rate*/) {throw InvalidArgumentException("Lapse rates can only be forced for linear regressions!", AT);};
virtual bool fit() = 0;
virtual double f(const double& x) = 0;
void getParams(std::vector<double>& o_coefficients);
std::string getName() {return regname;};
std::string getInfo();
virtual double f(const double& x) const = 0;
void getParams(std::vector<double>& o_coefficients) const;
std::string getName() const {return regname;};
std::string getInfo() const;
FitModel& operator =(const FitModel& source);
std::string toString() const;
protected:
......@@ -69,7 +69,7 @@ class FitLeastSquare : public FitModel {
FitLeastSquare() {};
void setData(const std::vector<double>& in_X, const std::vector<double>& in_Y);
bool fit();
virtual double f(const double& x) = 0;
virtual double f(const double& x) const = 0;
protected:
virtual void setDefaultGuess(); //set defaults guess values. Called by setData
......
......@@ -27,21 +27,6 @@ using namespace std;
namespace mio {
const double Interpol2D::wind_ys = 0.58;
const double Interpol2D::wind_yc = 0.42;
const double Interpol2D::bilin_inflection = 1200.;
double Interpol2D::getReferenceAltitude(const DEMObject& dem)
{
double ref_altitude = 1500.;
if(dem.min_altitude!=IOUtils::nodata && dem.max_altitude!=IOUtils::nodata) {
//we use the median elevation as the reference elevation for reprojections
ref_altitude = 0.5 * (dem.min_altitude+dem.max_altitude);
} else {
//since there is nothing else that we can do, we use an arbitrary elevation
ref_altitude = 1500.;
}
return ref_altitude;
}
//Usefull functions
/**
......@@ -119,6 +104,18 @@ void Interpol2D::getNeighbors(const double& x, const double& y,
sort (list.begin(), list.end());
}
//convert a vector of stations into two vectors of eastings and northings
void Interpol2D::buildPositionsVectors(const std::vector<StationData>& vecStations, std::vector<double>& vecEastings, std::vector<double>& vecNorthings)
{
vecEastings.resize( vecStations.size() );
vecNorthings.resize( vecStations.size() );
for (size_t i=0; i<vecStations.size(); i++) {
const Coords& position = vecStations[i].position;
vecEastings[i] = position.getEasting();
vecNorthings[i] = position.getNorthing();
}
}
//these weighting functions take the square of a distance as an argument and return a weight
inline double Interpol2D::weightInvDist(const double& d2)
{
......@@ -137,111 +134,6 @@ inline double Interpol2D::weightInvDistN(const double& d2)
return pow( Optim::invSqrt(d2) , dist_pow); //we use the optimized approximation for 1/sqrt
}
//Data regression models
/**
* @brief Computes the linear regression coefficients fitting the points given as X and Y in two vectors
* It relies on Interpol1D::NoisyLinRegression.
* @param in_X vector of X coordinates
* @param in_Y vector of Y coordinates (same order as X)
* @param coeffs a,b,r coefficients in a vector
* @return EXIT_SUCCESS or EXIT_FAILURE
*/
int Interpol2D::LinRegression(const std::vector<double>& in_X, const std::vector<double>& in_Y, std::vector<double>& coeffs)
{
std::stringstream mesg;
const int code = Interpol1D::NoisyLinRegression(in_X, in_Y, coeffs[1], coeffs[2], coeffs[3], mesg);
cerr << mesg.str();
return code;
}
//temporary solution while we migrate the regression classes
int Interpol2D::BiLinRegression(const std::vector<double>& in_X, const std::vector<double>& in_Y, std::vector<double>& coeffs)
{
std::vector<double> params;
const int code = Interpol1D::twoLinRegression(in_X, in_Y, Interpol2D::bilin_inflection, params);
coeffs[1] = params[1]; coeffs[2] = params[2]; coeffs[3] = 1.;
coeffs[4] = params[3]; coeffs[5] = params[4]; coeffs[6] = 1.;
return code;
}
//Now, the core interpolation functions: they project a given parameter to a reference altitude, given a constant lapse rate
//example: Ta projected to 1500m with a rate of -0.0065K/m
/**
* @brief Projects a given parameter to another elevation:
* This implementation keeps the value constant as a function of the elevation.
* This interface has to follow the interface of *LapseRateProjectPtr
* @param value original value
* @param altitude altitude of the original value
* @param new_altitude altitude of the reprojected value
* @param coeffs coefficients to use for the projection
* @return reprojected value
*/
double Interpol2D::ConstProject(const double& value, const double&, const double&, const std::vector<double>&)
{
return value;
}
/**
* @brief Projects a given parameter to another elevation:
* This implementation assumes a linear dependency of the value as a function of the elevation.
* This interface has to follow the interface of *LapseRateProjectPtr
* @param value original value
* @param altitude altitude of the original value
* @param new_altitude altitude of the reprojected value
* @param coeffs coefficients to use for the projection
* @return reprojected value
*/
double Interpol2D::LinProject(const double& value, const double& altitude, const double& new_altitude, const std::vector<double>& coeffs)
{
//linear lapse: coeffs must have been already computed
if (coeffs.empty()) {
throw IOException("Linear regression coefficients not initialized", AT);
}
return (value + coeffs[1] * (new_altitude - altitude));
}
/**
* @brief Projects a given parameter to another elevation:
* This implementation assumes a 2 segments linear dependency of the value as a function of the elevation.
* It uses Interpol2D::bilin_inflection as the inflection point altitude. This interface has to follow the interface of *LapseRateProjectPtr
* @param value original value
* @param altitude altitude of the original value
* @param new_altitude altitude of the reprojected value
* @param coeffs coefficients to use for the projection
* @return reprojected value
*/
double Interpol2D::BiLinProject(const double& value, const double& altitude, const double& new_altitude, const std::vector<double>& coeffs)
{
//linear lapse: coeffs must have been already computed
if (coeffs.empty()) {
throw IOException("Linear regression coefficients not initialized", AT);
}
if(altitude<=bilin_inflection)
return (value + coeffs[1] * (new_altitude - altitude));
else
return (value + coeffs[4] * (new_altitude - altitude));
}
/**
* @brief Projects a given parameter to another elevation:
* This implementation assumes that the value increases by a given fraction as a function of the elevation.
* This interface has to follow the interface of *LapseRateProjectPtr
* @param value original value
* @param altitude altitude of the original value
* @param new_altitude altitude of the reprojected value
* @param coeffs coefficients to use for the projection
* @return reprojected value
*/
double Interpol2D::FracProject(const double& value, const double& altitude, const double& new_altitude, const std::vector<double>& coeffs)
{
//linear lapse: coeffs must have been already computed
if (coeffs.empty()) {
throw IOException("Linear regression coefficients not initialized", AT);
}
return (value * (1. + coeffs[1] * (new_altitude - altitude)));
}
//Filling Functions
/**
* @brief Grid filling function:
......@@ -249,8 +141,8 @@ double Interpol2D::FracProject(const double& value, const double& altitude, cons
* @param dem array of elevations (dem)
* @param grid 2D array to fill
*/
void Interpol2D::stdPressureGrid2DFill(const DEMObject& dem, Grid2DObject& grid) {
grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);
void Interpol2D::stdPressure(const DEMObject& dem, Grid2DObject& grid) {
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++) {
......@@ -272,7 +164,7 @@ void Interpol2D::stdPressureGrid2DFill(const DEMObject& dem, Grid2DObject& grid)
* @param dem array of elevations (dem). This is needed in order to know if a point is "nodata"
* @param grid 2D array to fill
*/
void Interpol2D::constantGrid2DFill(const double& value, const DEMObject& dem, Grid2DObject& grid)
void Interpol2D::constant(const double& value, const DEMObject& dem, Grid2DObject& grid)
{
grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);
......@@ -288,49 +180,6 @@ void Interpol2D::constantGrid2DFill(const double& value, const DEMObject& dem, G
}
}
/**
* @brief Grid filling function:
* This implementation fills a flat grid with a constant value and then reprojects it to the terrain's elevation.
* for example, the air temperature measured at one point at 1500m would be given as value, the 1500m as altitude and the dem would allow to reproject this temperature on the full DEM using the detrending function provided as pointer (with its previously calculated coefficients).
* @param value value to put in the grid
* @param altitude altitude of the "value"
* @param dem array of elevations (dem)
* @param vecCoefficients vector of detrending coefficients
* @param funcptr detrending function pointer (that uses the detrending coefficients)
* @param grid 2D array to fill
*/
void Interpol2D::constantLapseGrid2DFill(const double& value, const double& altitude,
const DEMObject& dem, const std::vector<double>& vecCoefficients,
const LapseRateProjectPtr& funcptr, Grid2DObject& grid)
{
grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);
//fills a data table with constant values and then reprojects it to the DEM's elevation from a given altitude
//the laspe rate parameters must have been set before
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);
if (cell_altitude!=IOUtils::nodata) {
grid.grid2D(i,j) = funcptr(value, altitude, cell_altitude, vecCoefficients);
} else {
grid.grid2D(i,j) = IOUtils::nodata;
}
}
}
}
//convert a vector of stations into two vectors of eastings and northings
void Interpol2D::buildPositionsVectors(const std::vector<StationData>& vecStations, std::vector<double>& vecEastings, std::vector<double>& vecNorthings)
{
vecEastings.resize( vecStations.size() );
vecNorthings.resize( vecStations.size() );
for (size_t i=0; i<vecStations.size(); i++) {
const Coords& position = vecStations[i].position;
vecEastings[i] = position.getEasting();
vecNorthings[i] = position.getNorthing();
}
}
double Interpol2D::IDWCore(const double& x, const double& y, const std::vector<double>& vecData_in,
const std::vector<double>& vecEastings, const std::vector<double>& vecNorthings)
{
......@@ -340,69 +189,16 @@ double Interpol2D::IDWCore(const double& x, const double& y, const std::vector<d
const double scale = 1.e6;
for (size_t 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-vecEastings[i];
const double DY=y-vecNorthings[i];
const double weight = Optim::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;
}
return (parameter/norm); //normalization
}
/**
* @brief Grid filling function:
* This implementation fills a flat grid using Inverse Distance Weighting and then reproject it to the terrain's elevation.
* for example, the air temperatures measured at several stations would be given as values, the stations altitude and positions
* as positions and projected to a flat grid. Afterward, the grid would be reprojected to the correct elevation as given
* by the dem would using the detrending function provided as pointer (with its previously calculated coefficients).
* @param vecData_in input values to use for the IDW
* @param vecStations_in position of the "values" (altitude and coordinates)
* @param dem array of elevations (dem)
* @param vecCoefficients vector of detrending coefficients
* @param funcptr detrending function pointer (that uses the detrending coefficients)
* @param grid 2D array to fill
*/
void Interpol2D::LapseIDW(const std::vector<double>& vecData_in, const std::vector<StationData>& vecStations_in,
const DEMObject& dem, const std::vector<double>& vecCoefficients,
const LapseRateProjectPtr& funcptr,
Grid2DObject& grid)
{ //multiple source stations: lapse rate projection, IDW Krieging, re-projection
const double ref_altitude = getReferenceAltitude(dem);
const size_t 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 (size_t i=0; i<n_stations; i++) {
vecTref[i] = funcptr(vecData_in[i], vecStations_in[i].position.getAltitude(),
ref_altitude, vecCoefficients);
}
std::vector<double> vecEastings, vecNorthings;
buildPositionsVectors(vecStations_in, vecEastings, vecNorthings);