Skip to content

Commit

Permalink
Merge commit '3bf95c7f71988d65b3e3e2921ee777f9710a1ede'
Browse files Browse the repository at this point in the history
  • Loading branch information
@ftomei
ftomei committed Nov 20, 2024
2 parents e1ad81e + 3bf95c7 commit 485dd19
Show file tree
Hide file tree
Showing 27 changed files with 2,410 additions and 312 deletions.
31 changes: 31 additions & 0 deletions agrolib/gis/gis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -593,6 +593,37 @@ namespace gis
return sqrtf((dx * dx) + (dy * dy));
}

std::vector<float> computeEuclideanDistanceStation2Area(std::vector<std::vector<int>>& cells,std::vector<std::vector<int>>& stations)
{
// è possibile sapere in quale cella (row,col) si trova la stazione?
std::vector<float> distance(stations.size());
for (int i=0;i<stations.size();i++)
{
distance[i]= (stations[i][0] - cells[0][0])*(stations[i][0] - cells[0][0])+(stations[i][1] - cells[0][1])*(stations[i][1] - cells[0][1]);
for (int j=1;j<cells[i].size();j++)
{
distance[i] = MINVALUE(distance[i],(stations[i][0] - cells[j][0])*(stations[i][0] - cells[j][0])+(stations[i][1] - cells[j][1])*(stations[i][1] - cells[j][1]));
}
distance[i] = sqrt(1.*distance[i]);
}
return distance;
}

std::vector<int> computeMetropolisDistanceStation2Area(std::vector<std::vector<int>>& cells,std::vector<std::vector<int>>& stations)
{
// è possibile sapere in quale cella (row,col) si trova la stazione?
std::vector<int> distance(stations.size());
for (int i=0;i<stations.size();i++)
{
distance[i]=abs(stations[i][0] - cells[0][0])+abs(stations[i][1] - cells[0][1]);
for (int j=1;j<cells[i].size();j++)
{
distance[i] = MINVALUE(distance[i],abs(stations[i][0] - cells[j][0])+abs(stations[i][1] - cells[j][1]));
}
}
return distance;
}

void getRowColFromXY(const Crit3DRasterHeader& myHeader, double myX, double myY, int *row, int *col)
{
*row = (myHeader.nrRows - 1) - int(floor((myY - myHeader.llCorner.y) / myHeader.cellSize));
Expand Down
2 changes: 2 additions & 0 deletions agrolib/gis/gis.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -195,6 +195,8 @@

float computeDistance(float x1, float y1, float x2, float y2);
double computeDistancePoint(Crit3DUtmPoint *p0, Crit3DUtmPoint *p1);
std::vector<float> computeEuclideanDistanceStation2Area(std::vector<std::vector<int>>& cells,std::vector<std::vector<int>>& stations);
std::vector<int> computeMetropolisDistanceStation2Area(std::vector<std::vector<int>>& cells,std::vector<std::vector<int>>& stations);
bool updateMinMaxRasterGrid(Crit3DRasterGrid *rasterGrid);
void convertFlagToNodata(Crit3DRasterGrid& myGrid);
bool updateColorScale(Crit3DRasterGrid* rasterGrid, int row0, int col0, int row1, int col1);
Expand Down
189 changes: 162 additions & 27 deletions agrolib/interpolation/interpolation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@
#include "meteo.h"

#include <functional>

#include <cmath>

using namespace std;

Expand Down Expand Up @@ -1093,12 +1093,12 @@ void localSelection(vector <Crit3DInterpolationDataPoint> &inputPoints, vector <
for (unsigned long i = 0; i < inputPoints.size() ; i++)
inputPoints[i].distance = gis::computeDistance(x, y, float((inputPoints[i]).point->utm.x), float((inputPoints[i]).point->utm.y));

unsigned int nrValid = 0;
int nrValid = 0;
float stepRadius = 7500; // [m]
float r0 = 0; // [m]
float r1 = stepRadius; // [m]
unsigned int i;
unsigned int nrPrimaries = 0;
int nrPrimaries = 0;

int maxDistance = 0;
while ((!mySettings.getUseLapseRateCode() && nrValid < minPoints) || (mySettings.getUseLapseRateCode() && nrPrimaries < minPoints))
Expand All @@ -1122,15 +1122,9 @@ void localSelection(vector <Crit3DInterpolationDataPoint> &inputPoints, vector <
r1 += stepRadius;
}

double temp = 0;

if (maxDistance != 0)
for (i=0; i< selectedPoints.size(); i++)
{
/*temp = (MAXVALUE(1 - selectedPoints[i].distance / maxDistance, EPSILON))*1000;
temp = round(temp);
selectedPoints[i].regressionWeight = temp/1000;*/

selectedPoints[i].regressionWeight = MAXVALUE(1 - selectedPoints[i].distance / maxDistance, EPSILON);

//selectedPoints[i].regressionWeight = MAXVALUE(std::exp(-selectedPoints[i].distance*selectedPoints[i].distance/((0.5*maxDistance)*(0.5*maxDistance))),EPSILON);
Expand Down Expand Up @@ -1518,9 +1512,9 @@ void calculateFirstGuessCombinations(Crit3DProxy* myProxy)
std::vector<double> tempParam = myProxy->getFittingParametersRange();
std::vector <int> firstGuessPosition = myProxy->getFittingFirstGuess();
std::vector <double> tempFirstGuess;
int numSteps = 13;
int numSteps = 15;
std::vector <double> stepSize;
int nrParam = int(tempParam.size()/2);
unsigned nrParam = int(tempParam.size()/2);

double min_,max_;
for (unsigned j=0; j < nrParam; j++)
Expand Down Expand Up @@ -1727,9 +1721,9 @@ bool multipleDetrendingMain(std::vector <Crit3DInterpolationDataPoint> &myPoints
elevationPos = pos;
}

if (mySettings->getCurrentCombination().isProxyActive(elevationPos))
if (elevationPos != NODATA && mySettings->getCurrentCombination().isProxyActive(elevationPos))
{
if (!multipleDetrendingElevationFitting(elevationPos, myPoints, mySettings, myVar, errorStr))
if (!multipleDetrendingElevationFitting(elevationPos, myPoints, mySettings, myVar, errorStr, true))
return false;

detrendingElevation(elevationPos, myPoints, mySettings);
Expand All @@ -1745,7 +1739,7 @@ bool multipleDetrendingMain(std::vector <Crit3DInterpolationDataPoint> &myPoints
}

bool multipleDetrendingElevationFitting(int elevationPos, std::vector <Crit3DInterpolationDataPoint> &myPoints,
Crit3DInterpolationSettings* mySettings, meteoVariable myVar, std::string &errorStr)
Crit3DInterpolationSettings* mySettings, meteoVariable myVar, std::string &errorStr, bool isWeighted)
{
mySettings->getProxy(elevationPos)->setRegressionR2(NODATA);
if (! getUseDetrendingVar(myVar)) return true;
Expand Down Expand Up @@ -1792,7 +1786,7 @@ bool multipleDetrendingElevationFitting(int elevationPos, std::vector <Crit3DInt
std::vector<double> parametersDelta;
std::vector<double> parameters;
std::vector<double> stepSize;
int numSteps = 10;
int numSteps = 15;
std::function<double(double, std::vector<double>&)> myFunc;


Expand All @@ -1801,14 +1795,19 @@ bool multipleDetrendingElevationFitting(int elevationPos, std::vector <Crit3DInt
{
errorStr = "couldn't prepare the fitting parameters for proxy: elevation.";
return false;
}
}

std::vector<std::vector<double>> firstGuessCombinations = mySettings->getProxy(elevationPos)->getFirstGuessCombinations();
// multiple non linear fitting
double R2 = interpolation::bestFittingMarquardt_nDimension_singleFunction(*(myFunc.target<double(*)(double, std::vector<double>&)>()), 400, 4, parametersMin, parametersMax, parameters, parametersDelta,
stepSize, numSteps, 1000, 0.002, 0.005, predictors, predictands, weights,firstGuessCombinations);
double R2 = NODATA;
if (isWeighted)
R2 = interpolation::bestFittingMarquardt_nDimension_singleFunction(*(myFunc.target<double(*)(double, std::vector<double>&)>()), 4, parametersMin, parametersMax, parameters, parametersDelta,
1000, 0.002, 0.005, predictors, predictands, weights,firstGuessCombinations);
else
R2 = interpolation::bestFittingMarquardt_nDimension_singleFunction(*(myFunc.target<double(*)(double, std::vector<double>&)>()), 4, parametersMin, parametersMax, parameters, parametersDelta,
1000, 0.002, 0.005, predictors, predictands,firstGuessCombinations);

mySettings->getProxy(elevationPos)->setRegressionR2(R2);
mySettings->getProxy(elevationPos)->setRegressionR2(float(R2));

std::vector<std::vector<double>> newParameters;
newParameters.push_back(parameters);
Expand Down Expand Up @@ -1985,7 +1984,10 @@ bool multipleDetrendingOtherProxiesFitting(int elevationPos, std::vector <Crit3D

predictors.push_back(rowPredictors);
predictands.push_back(myPoints[i].value);
weights.push_back(myPoints[i].regressionWeight);
if (!isEqual(myPoints[i].regressionWeight, NODATA))
weights.push_back(myPoints[i].regressionWeight);
else
weights.push_back(1);
}

if (mySettings->getUseLocalDetrending() && othersPoints.size() < mySettings->getMinPointsLocalDetrending())
Expand Down Expand Up @@ -2038,8 +2040,12 @@ void detrendingOtherProxies(int elevationPos, std::vector<Crit3DInterpolationDat
myFunc.erase(myFunc.begin());
}

//if no other proxies are active/significant, return
if (parameters.empty()) return;

std::vector<std::vector<double>> fullParameters = parameters;
std::vector<std::function<double(double, std::vector<double>&)>> fullFunc = myFunc;

// detrending
float detrendValue;
for (int i = 0; i < myPoints.size(); i++)
Expand All @@ -2051,37 +2057,158 @@ void detrendingOtherProxies(int elevationPos, std::vector<Crit3DInterpolationDat
if (pos != elevationPos && mySettings->getCurrentCombination().isProxyActive(pos) && mySettings->getCurrentCombination().isProxySignificant(pos))
{
proxyValue = myPoints[i].getProxyValue(pos);
proxyValues.push_back(double(proxyValue));
if (!isEqual(proxyValue, NODATA))
proxyValues.push_back(double(proxyValue));
else {
parameters.erase(parameters.begin()+proxyValues.size());
myFunc.erase(myFunc.begin()+proxyValues.size());
}
}
}

detrendValue = float(functionSum(myFunc, proxyValues, parameters));
myPoints[i].value -= detrendValue;

myFunc = fullFunc;
parameters = fullParameters;
}
}

bool glocalDetrendingFitting(std::vector <Crit3DInterpolationDataPoint> &myPoints, Crit3DInterpolationSettings* mySettings, meteoVariable myVar, std::string& errorStr)
{
std::vector<Crit3DInterpolationDataPoint> subsetPoints;
std::vector<Crit3DMacroArea> macroAreas = mySettings->getMacroAreas();
int areasNr = macroAreas.size();
int i = 0;

int elevationPos = NODATA;
for (unsigned int pos=0; pos < mySettings->getSelectedCombination().getProxySize(); pos++)
{
if (getProxyPragaName(mySettings->getProxy(pos)->getName()) == proxyHeight)
elevationPos = pos;
}

//create the subset of points starting from the meteopoints vector saved in every macro area
for (i = 0; i < areasNr; i++) //ATTENZIONE INDICI ?
{
std::vector<int> temp = macroAreas[i].getMeteoPoints();
if (! temp.empty())
{
std::vector<Crit3DInterpolationDataPoint> subsetPoints;
for (int l = 0; l < temp.size(); l++)
{

for (int k = 0; k < myPoints.size(); k++)
if (myPoints[k].index == temp[l])
{
subsetPoints.push_back(myPoints[k]);
break;
}
}

mySettings->setCurrentCombination(mySettings->getSelectedCombination());
mySettings->clearFitting();

//fitting and elevation detrending (necessary for the fitting of other proxies)
if (elevationPos != NODATA && mySettings->getSelectedCombination().isProxyActive(elevationPos))
{
if (!multipleDetrendingElevationFitting(elevationPos, subsetPoints, mySettings, myVar, errorStr, false)) return false;
detrendingElevation(elevationPos, subsetPoints, mySettings);
}
if (!multipleDetrendingOtherProxiesFitting(elevationPos, subsetPoints, mySettings, myVar, errorStr))
return false;

//save parameters and combination in the macro area
macroAreas[i].setParameters(mySettings->getFittingParameters());
macroAreas[i].setCombination(mySettings->getCurrentCombination());
}
subsetPoints.clear();
}
//update macro areas in settings
mySettings->setMacroAreas(macroAreas);
return true;
}

double goldenSectionSearch(meteoVariable myVar,Crit3DMeteoPoint* &myMeteoPoints,int nrMeteoPoints,
std::vector <Crit3DInterpolationDataPoint> &interpolationPoints,Crit3DInterpolationSettings* mySettings,
Crit3DMeteoSettings* meteoSettings, double a, double b)
{
// this function finds the minimum by means of golden section method
double tol = 1;
//const double phi = (1 + std::sqrt(5)) / 2; // golden section
double x1 = b - (b - a) / GOLDEN_SECTION;
double x2 = a + (b - a) / GOLDEN_SECTION;
int counter=0;
while (std::abs(b - a) > tol && counter<100)
{
counter++;
if (topographicDistanceInternalFunction(myVar,myMeteoPoints,nrMeteoPoints,interpolationPoints, mySettings,meteoSettings, x1) <
topographicDistanceInternalFunction(myVar,myMeteoPoints,nrMeteoPoints,interpolationPoints, mySettings,meteoSettings,x2))
{
b = x2;
x2 = x1;
x1 = b - (b - a) / GOLDEN_SECTION;
mySettings->addToKhSeries(float(x1), topographicDistanceInternalFunction(myVar,myMeteoPoints,nrMeteoPoints,interpolationPoints, mySettings,meteoSettings, x1));
}
else
{
a = x1;
x1 = x2;
x2 = a + (b - a) / GOLDEN_SECTION;
mySettings->addToKhSeries(float(x2), topographicDistanceInternalFunction(myVar,myMeteoPoints,nrMeteoPoints,interpolationPoints, mySettings,meteoSettings, x2));
}
}
mySettings->addToKhSeries(float((a + b) / 2), topographicDistanceInternalFunction(myVar,myMeteoPoints,nrMeteoPoints,interpolationPoints, mySettings,meteoSettings, (a + b) / 2));
return (a + b) / 2; // approximated minimum
}


double topographicDistanceInternalFunction(meteoVariable myVar,
Crit3DMeteoPoint* &myMeteoPoints,
int nrMeteoPoints,
std::vector <Crit3DInterpolationDataPoint> &interpolationPoints,
Crit3DInterpolationSettings* mySettings,
Crit3DMeteoSettings* meteoSettings, double khFloat)
{
float avgError = 0;
int kh = int(khFloat);
mySettings->setTopoDist_Kh(kh);
if (computeResiduals(myVar, myMeteoPoints, nrMeteoPoints, interpolationPoints, mySettings, meteoSettings, true, true))
{
avgError = computeErrorCrossValidation(myMeteoPoints, nrMeteoPoints);
}
return avgError;
}

void topographicDistanceOptimize(meteoVariable myVar,
Crit3DMeteoPoint* &myMeteoPoints,
int nrMeteoPoints,
std::vector <Crit3DInterpolationDataPoint> &interpolationPoints,
Crit3DInterpolationSettings* mySettings,
Crit3DMeteoSettings* meteoSettings)
{
float avgError;
//float avgError;

mySettings->initializeKhSeries();

int kh = 0;
int bestKh = kh;
float bestError = NODATA;
//float kh = 0;
double bestKh = 0;
//float bestError = NODATA;
std::vector<float> avgErrorVec;
float khMin = 0;
float khMax = mySettings->getTopoDist_maxKh();
//kh = 0.5*(khMin+khMax);
bestKh = goldenSectionSearch(myVar, myMeteoPoints, nrMeteoPoints,interpolationPoints,mySettings,meteoSettings, khMin,khMax);
mySettings->setTopoDist_Kh(int(bestKh));
/*
while (kh <= mySettings->getTopoDist_maxKh())
{
mySettings->setTopoDist_Kh(kh);
if (computeResiduals(myVar, myMeteoPoints, nrMeteoPoints, interpolationPoints, mySettings, meteoSettings, true, true))
{
avgError = computeErrorCrossValidation(myMeteoPoints, nrMeteoPoints);
avgErrorVec.push_back(avgError);
if (isEqual(bestError, NODATA) || avgError < bestError)
{
bestError = avgError;
Expand All @@ -2094,6 +2221,7 @@ void topographicDistanceOptimize(meteoVariable myVar,
}
mySettings->setTopoDist_Kh(bestKh);
*/
}


Expand Down Expand Up @@ -2170,7 +2298,14 @@ bool preInterpolation(std::vector <Crit3DInterpolationDataPoint> &myPoints, Crit
if (!mySettings->getUseLocalDetrending())
setMultipleDetrendingHeightTemperatureRange(mySettings);

if (! multipleDetrendingMain(myPoints, mySettings, myVar, errorStr)) return false;
if (mySettings->getUseGlocalDetrending())
{
glocalDetrendingFitting(myPoints, mySettings, myVar, errorStr);
}
else
{
if (! multipleDetrendingMain(myPoints, mySettings, myVar, errorStr)) return false;
}
}
else
{
Expand All @@ -2181,7 +2316,7 @@ bool preInterpolation(std::vector <Crit3DInterpolationDataPoint> &myPoints, Crit
}
}

if (mySettings->getUseTD() && getUseTdVar(myVar))
if (mySettings->getUseTD() && getUseTdVar(myVar) && ! mySettings->getUseGlocalDetrending())
{
topographicDistanceOptimize(myVar, myMeteoPoints, nrMeteoPoints, myPoints, mySettings, meteoSettings);
}
Expand Down
Loading

0 comments on commit 485dd19

Please sign in to comment.