Merge branch 'master' of https://github.com/ARPA-SIMC/PRAGA

agrolib/interpolation/interpolation.cpp

@@ -1021,7 +1021,7 @@ float gaussWeighted(vector <Crit3DInterpolationDataPoint> &myPointList)
 }
 */
 
-void localSelection_new(std::vector<Crit3DInterpolationDataPoint> &inputPoints, std::vector<Crit3DInterpolationDataPoint> &selectedPoints,
+/*void localSelection_new(std::vector<Crit3DInterpolationDataPoint> &inputPoints, std::vector<Crit3DInterpolationDataPoint> &selectedPoints,
                     float x, float y, float z, Crit3DInterpolationSettings& mySettings)
 {
     std::vector<Crit3DInterpolationDataPoint> tempPoints;
@@ -1078,7 +1078,7 @@ void localSelection_new(std::vector<Crit3DInterpolationDataPoint> &inputPoints,
 
     return;
 }
-
+*/
 
 
 
@@ -1432,13 +1432,12 @@ bool proxyValidity(std::vector <Crit3DInterpolationDataPoint> &myPoints, int pro
 bool proxyValidityWeighted(std::vector <Crit3DInterpolationDataPoint> &myPoints, int proxyPos, float stdDevThreshold, double &avg, double &stdDev)
 {
     std::vector<float> proxyValues;
-    const int MIN_NR = 10;
 
     std::vector<double> data, weights;
     data.resize(myPoints.size());
     weights.resize(myPoints.size());
 
-    for (int i = 0; i < myPoints.size(); i++)
+    for (unsigned i = 0; i < myPoints.size(); i++)
     {
         data[i] = myPoints[i].getProxyValue(proxyPos);
         weights[i] = myPoints[i].regressionWeight;
@@ -1780,7 +1779,6 @@ bool multipleDetrendingElevation(Crit3DProxyCombination elevationCombination, st
 
     // proxy spatial variability (1st step)
     double avg, stdDev;
-    unsigned validNr;
     validNr = 0;
 
     if (proxyValidityWeighted(elevationPoints, elevationPos, elevationProxy->getStdDevThreshold(), avg, stdDev))
@@ -1874,7 +1872,7 @@ bool multipleDetrendingElevation(Crit3DProxyCombination elevationCombination, st
         }
     }
 
-    if (mySettings->getUseLocalDetrending() && elevationPoints.size() < mySettings->getMinPointsLocalDetrending())
+    if (mySettings->getUseLocalDetrending() && elevationPoints.size() < unsigned(mySettings->getMinPointsLocalDetrending()))
     {
         elevationProxy->setIsSignificant(false);
         Crit3DProxyCombination myCombination = mySettings->getSelectedCombination();

agrolib/interpolation/interpolationSettings.cpp

@@ -374,7 +374,7 @@ void Crit3DInterpolationSettings::setChosenElevationFunction(TFittingFunction ch
     {
         if (chosenFunction == getProxy(elPos)->getFittingFunctionName() && !getProxy(elPos)->getFittingParametersRange().empty())
         {
-            std::vector tempParam = getProxy(elPos)->getFittingParametersRange();
+            std::vector <double> tempParam = getProxy(elPos)->getFittingParametersRange();
 
             if (chosenFunction == piecewiseTwo)
             {

agrolib/interpolation/spatialControl.cpp

@@ -133,6 +133,71 @@ bool computeResiduals(meteoVariable myVar, Crit3DMeteoPoint* meteoPoints, int nr
     return true;
 }
 
+bool computeResidualsLocalDetrending(meteoVariable myVar, Crit3DTime myTime, Crit3DMeteoPoint* meteoPoints, int nrMeteoPoints,
+                                              std::vector <Crit3DInterpolationDataPoint> &interpolationPoints, Crit3DInterpolationSettings* settings,
+                                              Crit3DMeteoSettings* meteoSettings, Crit3DClimateParameters* climateParameters,
+                                     bool excludeOutsideDem, bool excludeSupplemental)
+{
+
+    if (myVar == noMeteoVar) return false;
+
+    std::vector <double> myProxyValues;
+    bool isValid;
+    std::string errorStdString;
+
+    for (int i = 0; i < nrMeteoPoints; i++)
+    {
+        myProxyValues = meteoPoints[i].getProxyValues();
+
+        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)
+        {
+            float myValue = meteoPoints[i].currentValue;
+
+            std::vector <Crit3DInterpolationDataPoint> subsetInterpolationPoints;
+            localSelection(interpolationPoints, subsetInterpolationPoints, meteoPoints->point.utm.x, meteoPoints->point.utm.y, meteoPoints->point.z, *settings);
+            if (! preInterpolation(subsetInterpolationPoints, settings, meteoSettings,
+                                  climateParameters, meteoPoints, nrMeteoPoints, myVar, myTime, errorStdString))
+            {
+                return false;
+            }
+
+            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 (myValue != NODATA)
+                {
+                    if (myValue < meteoSettings->getRainfallThreshold())
+                        myValue=0.;
+                }
+
+                if (interpolatedValue != NODATA)
+                {
+                    if (interpolatedValue < meteoSettings->getRainfallThreshold())
+                        interpolatedValue=0.;
+                }
+            }
+
+            // TODO derived var
+
+            if ((interpolatedValue != NODATA) && (myValue != NODATA))
+            {
+                meteoPoints[i].residual = interpolatedValue - myValue;
+            }
+        }
+    }
+
+    return true;
+}
 
 float computeErrorCrossValidation(Crit3DMeteoPoint* myPoints, int nrMeteoPoints)
 {

agrolib/interpolation/spatialControl.h

@@ -30,6 +30,9 @@
 
     bool computeResiduals(meteoVariable myVar, Crit3DMeteoPoint* meteoPoints, int nrMeteoPoints,
                           std::vector <Crit3DInterpolationDataPoint> &interpolationPoints, Crit3DInterpolationSettings* settings, Crit3DMeteoSettings* meteoSettings, bool excludeOutsideDem, bool excludeSupplemental);
+    bool computeResidualsLocalDetrending(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);
 

agrolib/project/project.cpp

@@ -2248,7 +2248,6 @@ bool Project::computeStatisticsCrossValidation(crossValidationStatistics* myStat
     return true;
 }
 
-
 bool Project::interpolationCv(meteoVariable myVar, const Crit3DTime& myTime, crossValidationStatistics *myStats)
 {
     // TODO local detrending
@@ -2263,7 +2262,7 @@ bool Project::interpolationCv(meteoVariable myVar, const Crit3DTime& myTime, cro
         myVar == airRelHumidity))
     {
         logError("Cross validation is not available for " + QString::fromStdString(getVariableString(myVar))
-                 + "\n Deactive 'Interpolate relative humidity using dew point' option.");
+                 + "\n Deactivate option 'Interpolate relative humidity using dew point'");
         return false;
     }
 
@@ -2294,15 +2293,24 @@ bool Project::interpolationCv(meteoVariable myVar, const Crit3DTime& myTime, cro
     if (interpolationSettings.getUseMultipleDetrending())
         interpolationSettings.clearFitting();
 
-    if (! preInterpolation(interpolationPoints, &interpolationSettings, meteoSettings, &climateParameters,
+    if (! interpolationSettings.getUseLocalDetrending() && ! preInterpolation(interpolationPoints, &interpolationSettings, meteoSettings, &climateParameters,
                           meteoPoints, nrMeteoPoints, myVar, myTime, errorStdStr))
     {
         logError("Error in function preInterpolation:\n" + QString::fromStdString(errorStdStr));
         return false;
     }
 
-    if (! computeResiduals(myVar, meteoPoints, nrMeteoPoints, interpolationPoints, &interpolationSettings, meteoSettings, true, true))
-        return false;
+    if (! interpolationSettings.getUseLocalDetrending())
+    {
+        if (! computeResiduals(myVar, meteoPoints, nrMeteoPoints, interpolationPoints, &interpolationSettings, meteoSettings, true, true))
+            return false;
+    }
+    else
+    {
+        if (! computeResidualsLocalDetrending(myVar, myTime, meteoPoints, nrMeteoPoints, interpolationPoints,
+                                             &interpolationSettings, meteoSettings, &climateParameters, true, true))
+            return false;
+    }
 
     if (! computeStatisticsCrossValidation(myStats))
         return false;
@@ -2311,6 +2319,7 @@ bool Project::interpolationCv(meteoVariable myVar, const Crit3DTime& myTime, cro
 }
 
 
+
 bool Project::interpolationDem(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster)
 {
     std::vector <Crit3DInterpolationDataPoint> interpolationPoints;

src/mainWindow.cpp

@@ -3114,36 +3114,20 @@ void MainWindow::on_actionInterpolationCrossValidation_triggered()
 
                         if  (myProxy->getIsSignificant())
                         {
-                            if (getProxyPragaName(myProxy->getName()) == proxyHeight)
+                            cvOutput << myProxy->getName() << std::endl;
+
+                            if (! myProject.interpolationSettings.getUseLocalDetrending())
                             {
-                                if (myProxy->getFittingFunctionName() == piecewiseTwo)
+                                if (getProxyPragaName(myProxy->getName()) == proxyHeight)
                                 {
-                                    cvOutput << "h0: " << par[i][0] << std::endl;
-                                    cvOutput << "t0: " << par[i][1] << std::endl;
-                                    cvOutput << "slope0: " << par[i][2] << std::endl;
-                                    cvOutput << "slope1: " << par[i][3] << std::endl;
+                                    cvOutput << "function: " << getKeyStringElevationFunction(myProxy->getFittingFunctionName()) << std::endl;
+                                    for (int j=0; j < par[i].size(); j++)
+                                        cvOutput << "par" << j << ": " << par[i][j] << std::endl;
+
                                 }
                             }
                         }
                     }
-                            /*
-
-
-                                if (interpolationSettings->getProxy(proxyPos)->getFittingFunctionName() == piecewiseThree)
-                                {
-                                    std::vector <double> xVector;
-                                    for (int m = xMin; m < xMax; m += 5)
-                                        xVector.push_back(m);
-
-                                    for (int p = 0; p < int(xVector.size()); p++)
-                                    {
-                                        point.setX(xVector[p]);
-                                        point.setY(lapseRatePiecewise_three(xVector[p], parameters[proxyPos]));
-                                        point_vector.append(point);
-                                    }
-                                }
-
-                    }*/
                 }
             }
         }