Skip to content

Commit

Permalink
Merge commit 'b4188c5033669bf55546ce55b47c16220f974eb3'
Browse files Browse the repository at this point in the history
  • Loading branch information
@ftomei
ftomei committed Jan 22, 2025
2 parents df3af33 + b4188c5 commit 0513169
Show file tree
Hide file tree
Showing 6 changed files with 274 additions and 123 deletions.
63 changes: 63 additions & 0 deletions agrolib/interpolation/interpolation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1754,6 +1754,69 @@ std::vector <double> getfittingParameters(Crit3DProxyCombination myCombination,
return myParam;
}

void macroAreaDetrending(Crit3DMacroArea myArea, meteoVariable myVar, Crit3DInterpolationSettings &mySettings, Crit3DMeteoSettings* meteoSettings, Crit3DMeteoPoint* meteoPoints, std::vector <Crit3DInterpolationDataPoint> interpolationPoints, std::vector <Crit3DInterpolationDataPoint> &subsetInterpolationPoints, int elevationPos)
{
//take the parameters+combination for that area
mySettings.setFittingParameters(myArea.getParameters());
mySettings.setCurrentCombination(myArea.getCombination());

//find the fitting functions vector based on the length of the parameters vector for every proxy
std::vector<std::function<double (double, std::vector<double> &)> > fittingFunction;
for (int l = 0; l < myArea.getParameters().size(); l++)
{
if (myArea.getParameters()[l].size() == 2)
fittingFunction.push_back(functionLinear_intercept);
else if (myArea.getParameters()[l].size() == 4)
fittingFunction.push_back(lapseRatePiecewise_two);
else if (myArea.getParameters()[l].size() == 5)
fittingFunction.push_back(lapseRatePiecewise_three);
else if (myArea.getParameters()[l].size() == 6)
fittingFunction.push_back(lapseRatePiecewise_three_free);
}

mySettings.setFittingFunction(fittingFunction);

// create vector of macro area interpolation points
std::vector<int> temp = myArea.getMeteoPoints();
for (int l = 0; l < temp.size(); l++)
{
for (int k = 0; k < interpolationPoints.size(); k++)
{
if (interpolationPoints[k].index == temp[l])
{
subsetInterpolationPoints.push_back(interpolationPoints[k]);
}
}
}

//detrending
if (elevationPos != NODATA && myArea.getCombination().isProxyActive(elevationPos) && myArea.getCombination().isProxySignificant(elevationPos))
{
detrendingElevation(elevationPos, subsetInterpolationPoints, &mySettings);
}

detrendingOtherProxies(elevationPos, subsetInterpolationPoints, &mySettings);

Crit3DMeteoPoint* myMeteoPoints = new Crit3DMeteoPoint[unsigned(myArea.getMeteoPointsNr())];
std::vector<int> meteoPointsList = myArea.getMeteoPoints();

for (unsigned i = 0; i < meteoPointsList.size(); i++)
{
myMeteoPoints[i] = meteoPoints[meteoPointsList[i]];
}

if (mySettings.getUseTD() && getUseTdVar(myVar))
{
topographicDistanceOptimize(myVar, myMeteoPoints, myArea.getMeteoPointsNr(),
subsetInterpolationPoints, &mySettings, meteoSettings);
}

myMeteoPoints->clear();
delete[] myMeteoPoints;

return;
}

bool multipleDetrendingMain(std::vector <Crit3DInterpolationDataPoint> &myPoints,
Crit3DInterpolationSettings* mySettings, meteoVariable myVar, std::string &errorStr)
{
Expand Down
4 changes: 4 additions & 0 deletions agrolib/interpolation/interpolation.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,10 @@
Crit3DInterpolationSettings *mySettings, Crit3DClimateParameters *myClimate,
meteoVariable myVar, Crit3DTime myTime);

void macroAreaDetrending(Crit3DMacroArea myArea, meteoVariable myVar, Crit3DInterpolationSettings &mySettings, Crit3DMeteoSettings* meteoSettings,
Crit3DMeteoPoint* meteoPoints, std::vector <Crit3DInterpolationDataPoint> interpolationPoints,
std::vector <Crit3DInterpolationDataPoint> &subsetInterpolationPoints, int elevationPos);

bool multipleDetrendingMain(std::vector <Crit3DInterpolationDataPoint> &myPoints,
Crit3DInterpolationSettings* mySettings, meteoVariable myVar, std::string &errorStr);

Expand Down
232 changes: 179 additions & 53 deletions agrolib/interpolation/spatialControl.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,39 +94,42 @@ bool computeResiduals(meteoVariable myVar, Crit3DMeteoPoint* meteoPoints, int nr

meteoPoints[i].residual = NODATA;

isValid = (! excludeSupplemental || checkLapseRateCode(meteoPoints[i].lapseRateCode, settings->getUseLapseRateCode(), false));
isValid = (isValid && (! excludeOutsideDem || meteoPoints[i].isInsideDem));

if (isValid && meteoPoints[i].quality == quality::accepted)
if (meteoPoints[i].active)
{
float myValue = meteoPoints[i].currentValue;

float interpolatedValue = interpolate(interpolationPoints, settings, meteoSettings, myVar,
float(meteoPoints[i].point.utm.x),
float(meteoPoints[i].point.utm.y),
float(meteoPoints[i].point.z),
myProxyValues, false);
isValid = (! excludeSupplemental || checkLapseRateCode(meteoPoints[i].lapseRateCode, settings->getUseLapseRateCode(), false));
isValid = (isValid && (! excludeOutsideDem || meteoPoints[i].isInsideDem));

if ( myVar == precipitation || myVar == dailyPrecipitation)
if (isValid && meteoPoints[i].quality == quality::accepted)
{
if (myValue != NODATA)
{
if (myValue < meteoSettings->getRainfallThreshold())
myValue=0.;
}
float myValue = meteoPoints[i].currentValue;

if (interpolatedValue != NODATA)
float interpolatedValue = interpolate(interpolationPoints, settings, meteoSettings, myVar,
float(meteoPoints[i].point.utm.x),
float(meteoPoints[i].point.utm.y),
float(meteoPoints[i].point.z),
myProxyValues, false);

if ( myVar == precipitation || myVar == dailyPrecipitation)
{
if (interpolatedValue < meteoSettings->getRainfallThreshold())
interpolatedValue=0.;
if (myValue != NODATA)
{
if (myValue < meteoSettings->getRainfallThreshold())
myValue=0.;
}

if (interpolatedValue != NODATA)
{
if (interpolatedValue < meteoSettings->getRainfallThreshold())
interpolatedValue=0.;
}
}
}

// TODO derived var
// TODO derived var

if ((interpolatedValue != NODATA) && (myValue != NODATA))
{
meteoPoints[i].residual = myValue - interpolatedValue;
if ((interpolatedValue != NODATA) && (myValue != NODATA))
{
meteoPoints[i].residual = myValue - interpolatedValue;
}
}
}
}
Expand All @@ -153,48 +156,171 @@ bool computeResidualsLocalDetrending(meteoVariable myVar, Crit3DTime myTime, Cri

meteoPoints[i].residual = NODATA;

isValid = (! excludeSupplemental || checkLapseRateCode(meteoPoints[i].lapseRateCode, settings->getUseLapseRateCode(), false));
isValid = (isValid && (! excludeOutsideDem || meteoPoints[i].isInsideDem));

if (isValid && meteoPoints[i].quality == quality::accepted)
if (meteoPoints[i].active)
{
float myValue = meteoPoints[i].currentValue;

std::vector <Crit3DInterpolationDataPoint> subsetInterpolationPoints;
localSelection(interpolationPoints, subsetInterpolationPoints, float(meteoPoints[i].point.utm.x),
float(meteoPoints[i].point.utm.y), *settings, false);
if (! preInterpolation(subsetInterpolationPoints, settings, meteoSettings,
climateParameters, meteoPoints, nrMeteoPoints, myVar, myTime, errorStdString))
isValid = (! excludeSupplemental || checkLapseRateCode(meteoPoints[i].lapseRateCode, settings->getUseLapseRateCode(), false));
isValid = (isValid && (! excludeOutsideDem || meteoPoints[i].isInsideDem));

if (isValid && meteoPoints[i].quality == quality::accepted)
{
return false;
}
float myValue = meteoPoints[i].currentValue;

float interpolatedValue = interpolate(subsetInterpolationPoints, settings, meteoSettings, myVar,
float(meteoPoints[i].point.utm.x),
float(meteoPoints[i].point.utm.y),
float(meteoPoints[i].point.z),
myProxyValues, false);
std::vector <Crit3DInterpolationDataPoint> subsetInterpolationPoints;
localSelection(interpolationPoints, subsetInterpolationPoints, float(meteoPoints[i].point.utm.x),
float(meteoPoints[i].point.utm.y), *settings, false);
if (! preInterpolation(subsetInterpolationPoints, settings, meteoSettings,
climateParameters, meteoPoints, nrMeteoPoints, myVar, myTime, errorStdString))
{
return false;
}

if ( myVar == precipitation || myVar == dailyPrecipitation)
{
if (myValue != NODATA)
float interpolatedValue = interpolate(subsetInterpolationPoints, settings, meteoSettings, myVar,
float(meteoPoints[i].point.utm.x),
float(meteoPoints[i].point.utm.y),
float(meteoPoints[i].point.z),
myProxyValues, false);

if ( myVar == precipitation || myVar == dailyPrecipitation)
{
if (myValue < meteoSettings->getRainfallThreshold())
myValue=0.;
if (myValue != NODATA)
{
if (myValue < meteoSettings->getRainfallThreshold())
myValue=0.;
}

if (interpolatedValue != NODATA)
{
if (interpolatedValue < meteoSettings->getRainfallThreshold())
interpolatedValue=0.;
}
}

if (interpolatedValue != NODATA)
// TODO derived var

if ((interpolatedValue != NODATA) && (myValue != NODATA))
{
if (interpolatedValue < meteoSettings->getRainfallThreshold())
interpolatedValue=0.;
meteoPoints[i].residual = myValue - interpolatedValue;
}
}
}
}

return true;
}

bool computeResidualsGlocalDetrending(meteoVariable myVar, Crit3DTime myTime, Crit3DMeteoPoint* meteoPoints, int nrMeteoPoints,
std::vector <Crit3DInterpolationDataPoint> &interpolationPoints, Crit3DInterpolationSettings* settings,
Crit3DMeteoSettings* meteoSettings, Crit3DClimateParameters* climateParameters,
bool excludeOutsideDem, bool excludeSupplemental)
{

// TODO derived var
//TODO: glocal cv with grid ONLY (no DEM)

if ((interpolatedValue != NODATA) && (myValue != NODATA))
if (myVar == noMeteoVar) return false;

std::vector <double> myProxyValues;
bool isValid;
std::string errorStdString;
std::vector <Crit3DMacroArea> macroAreas = settings->getMacroAreas();

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

for (int j = 0; j < nrMeteoPoints; j++)
{
meteoPoints[j].residual = NODATA;
}

//ciclo sulle aree
for (int k = 0; k < macroAreas.size(); k++)
{
Crit3DMacroArea myArea = macroAreas[k];
std::vector<Crit3DInterpolationDataPoint> areaInterpolationPoints;
std::vector<int> meteoPointsList = myArea.getMeteoPoints();
std::vector<float> areaCells;

//if (! myArea.getAreaCellsGrid().empty() || ! myArea.getAreaCellsDEM().empty() )
if (! myArea.getAreaCellsDEM().empty())
{
if (! meteoPointsList.empty())
{
meteoPoints[i].residual = myValue - interpolatedValue;
//un solo detrending per ogni area
macroAreaDetrending(myArea, myVar, *settings, meteoSettings, meteoPoints, interpolationPoints, areaInterpolationPoints, elevationPos);

//ciclo sui meteopoint dell'area
for (int i = 0; i < meteoPointsList.size(); i++)
{
myProxyValues = meteoPoints[meteoPointsList[i]].getProxyValues();

if (meteoPoints[meteoPointsList[i]].active)
{
//peso della stazione nell'area attuale in base alla sua posizione
int row, col;
float weight = NODATA;
int temp = NODATA;

//valido solo per DEM
areaCells = myArea.getAreaCellsDEM();
gis::Crit3DPoint point = meteoPoints[meteoPointsList[i]].point;
std::string name = meteoPoints[meteoPointsList[i]].name;
std::string id = meteoPoints[meteoPointsList[i]].id;

gis::getRowColFromXY(*settings->getCurrentDEM()->header, meteoPoints[meteoPointsList[i]].point.utm, &row, &col);
temp = settings->getCurrentDEM()->header->nrCols*row + col;

for (int k = 0; k < areaCells.size(); k = k + 2)
{
if (areaCells[k] == temp)
weight = areaCells[k+1];
}

isValid = (! excludeSupplemental || checkLapseRateCode(meteoPoints[meteoPointsList[i]].lapseRateCode, settings->getUseLapseRateCode(), false));
isValid = (isValid && (! excludeOutsideDem || meteoPoints[meteoPointsList[i]].isInsideDem));

if (isValid && meteoPoints[meteoPointsList[i]].quality == quality::accepted)
{
float myValue = meteoPoints[meteoPointsList[i]].currentValue;

float interpolatedValue = interpolate(areaInterpolationPoints, settings, meteoSettings, myVar,
float(meteoPoints[meteoPointsList[i]].point.utm.x),
float(meteoPoints[meteoPointsList[i]].point.utm.y),
float(meteoPoints[meteoPointsList[i]].point.z),
myProxyValues, false);

if ( myVar == precipitation || myVar == dailyPrecipitation)
{
if (myValue != NODATA)
{
if (myValue < meteoSettings->getRainfallThreshold())
myValue=0.;
}

if (interpolatedValue != NODATA)
{
if (interpolatedValue < meteoSettings->getRainfallThreshold())
interpolatedValue=0.;
}
}

// TODO derived var

if (!isEqual(interpolatedValue, NODATA) && !isEqual(myValue, NODATA) && !isEqual(weight, NODATA))
{
if (isEqual(meteoPoints[meteoPointsList[i]].residual, NODATA))
meteoPoints[meteoPointsList[i]].residual = (myValue - interpolatedValue)*weight;
else
{
meteoPoints[meteoPointsList[i]].residual += (myValue - interpolatedValue)*weight;
}
}
}
}
}
}
}
}
Expand Down
5 changes: 5 additions & 0 deletions agrolib/interpolation/spatialControl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,11 @@
std::vector <Crit3DInterpolationDataPoint> &interpolationPoints, Crit3DInterpolationSettings* settings,
Crit3DMeteoSettings* meteoSettings, Crit3DClimateParameters *climateParameters, bool excludeOutsideDem, bool excludeSupplemental);

bool computeResidualsGlocalDetrending(meteoVariable myVar, Crit3DTime myTime, Crit3DMeteoPoint* meteoPoints, int nrMeteoPoints,
std::vector <Crit3DInterpolationDataPoint> &interpolationPoints, Crit3DInterpolationSettings* settings,
Crit3DMeteoSettings* meteoSettings, Crit3DClimateParameters* climateParameters,
bool excludeOutsideDem, bool excludeSupplemental);

float computeErrorCrossValidation(Crit3DMeteoPoint *myPoints, int nrMeteoPoints);

bool spatialQualityControl(meteoVariable myVar, Crit3DMeteoPoint* meteoPoints, int nrMeteoPoints,
Expand Down
Loading

0 comments on commit 0513169

Please sign in to comment.