Merge commit '18cfb88ea417fd506456bdef346b6d48104ee9ef'

ARPA-SIMC · Sep 8, 2023 · 0536c22 · 0536c22
2 parents e50991c + 18cfb88
commit 0536c22
Show file tree

Hide file tree

Showing 2 changed files with 56 additions and 27 deletions.
diff --git a/agrolib/interpolation/interpolation.cpp b/agrolib/interpolation/interpolation.cpp
@@ -23,6 +23,7 @@
     [email protected]
 */
 
+#include <ostream>
 #include <stdlib.h>
 #include <math.h>
 #include <vector>
@@ -37,6 +38,7 @@
 #include "spatialControl.h"
 #include "interpolation.h"
 
+#include <iostream>
 
 using namespace std;
 
@@ -1320,7 +1322,10 @@ bool proxyValidity(std::vector <Crit3DInterpolationDataPoint> &myPoints, int pro
 
     *stdDev = float(sqrt(sum / (proxyValues.size() - 1)));
 
-    return (*stdDev > stdDevThreshold);
+    if (stdDevThreshold != NODATA)
+        return (*stdDev > stdDevThreshold);
+    else
+        return true;
 }
 
 Crit3DProxyCombination multipleDetrending(std::vector <Crit3DInterpolationDataPoint> &myPoints,
@@ -1340,23 +1345,19 @@ Crit3DProxyCombination multipleDetrending(std::vector <Crit3DInterpolationDataPo
         }
     }
 
-    if (nrPredictors == 0) return myCombination;
+    if (nrPredictors == 0) return myCombination;    
 
-    // proxy dispersion
+    // proxy spatial variability (1st step)
     float avg, stdDev;
-    std::vector <float> avgs;
-    std::vector <float> stdDevs;
     unsigned validNr;
     Crit3DProxyCombination outCombination = myCombination;
     unsigned pos;
-
     validNr = 0;
+
     for (pos=0; pos < int(mySettings->getProxyNr()); pos++)
     {
         if (myCombination.getValue(pos) && proxyValidity(myPoints, pos, mySettings->getProxy(pos)->getStdDevThreshold(), &avg, &stdDev))
         {
-            avgs.push_back(avg);
-            stdDevs.push_back(stdDev);
             outCombination.setValue(pos, true);
             validNr++;
         }
@@ -1366,22 +1367,15 @@ Crit3DProxyCombination multipleDetrending(std::vector <Crit3DInterpolationDataPo
 
     if (validNr == 0) return outCombination;
 
-    // z-score normalization
-    std::vector <float> rowPredictors;
-    std::vector <std::vector <float>> predictorsNorm;
-    std::vector <float> predictands;
-    std::vector <float> weights;
+    // exclude points with incomplete proxies
+    unsigned i;
+    std::vector <Crit3DInterpolationDataPoint> finalPoints;
     bool isValid;
     float proxyValue;
-    unsigned index = 0;
-    const int MIN_NR = 10;
-    unsigned i;
 
     for (i=0; i < myPoints.size(); i++)
     {
         isValid = true;
-        rowPredictors.clear();
-        index = 0;
         for (pos=0; pos < mySettings->getProxyNr(); pos++)
             if (outCombination.getValue(pos))
             {
@@ -1391,18 +1385,53 @@ Crit3DProxyCombination multipleDetrending(std::vector <Crit3DInterpolationDataPo
                     isValid = false;
                     break;
                 }
-                else
-                    rowPredictors.push_back((proxyValue - avgs[index]) / stdDevs[index]);
-
-                index++;
             }
 
-        if (isValid)
+        if (isValid) finalPoints.push_back(myPoints[i]);
+    }
+
+    // proxy spatial variability (2nd step)
+    std::vector <float> avgs;
+    std::vector <float> stdDevs;
+
+    for (pos=0; pos < int(mySettings->getProxyNr()); pos++)
+    {
+        if (myCombination.getValue(pos) && proxyValidity(finalPoints, pos, mySettings->getProxy(pos)->getStdDevThreshold(), &avg, &stdDev))
         {
-            predictorsNorm.push_back(rowPredictors);
-            predictands.push_back(myPoints[i].value);
-            weights.push_back(myPoints[i].regressionWeight);
+            avgs.push_back(avg);
+            stdDevs.push_back(stdDev);
+            outCombination.setValue(pos, true);
+            validNr++;
         }
+        else
+            outCombination.setValue(pos, false);
+    }
+
+    if (validNr == 0) return outCombination;
+
+    // z-score normalization
+    std::vector <float> rowPredictors;
+    std::vector <std::vector <float>> predictorsNorm;
+    std::vector <float> predictands;
+    std::vector <float> weights;
+    unsigned index = 0;
+    const int MIN_NR = 10;
+
+    for (i=0; i < finalPoints.size(); i++)
+    {
+        rowPredictors.clear();
+        index = 0;
+        for (pos=0; pos < mySettings->getProxyNr(); pos++)
+            if (outCombination.getValue(pos))
+            {
+                proxyValue = finalPoints[i].getProxyValue(pos);
+                rowPredictors.push_back((proxyValue - avgs[index]) / stdDevs[index]);
+                index++;
+            }
+
+        predictorsNorm.push_back(rowPredictors);
+        predictands.push_back(finalPoints[i].value);
+        weights.push_back(finalPoints[i].regressionWeight);
     }
 
     if (predictorsNorm.size() < MIN_NR)

diff --git a/agrolib/project/project.cpp b/agrolib/project/project.cpp
@@ -2225,7 +2225,6 @@ bool Project::interpolationDemLocalDetrending(meteoVariable myVar, const Crit3DT
                 myRaster->getRowCol(x, y, row, col);
                 if (! myRaster->isOutOfGrid(row, col))
                 {
-
                     std::vector <Crit3DInterpolationDataPoint> subsetInterpolationPoints;
                     localSelection(interpolationPoints, subsetInterpolationPoints, x, y, interpolationSettings);
                     preInterpolation(subsetInterpolationPoints, &interpolationSettings, meteoSettings, &climateParameters, meteoPoints, nrMeteoPoints, myVar, myTime);
@@ -2848,6 +2847,7 @@ void Project::saveProxies()
             if (myProxy->getProxyTable() != "") parameters->setValue("table", QString::fromStdString(myProxy->getProxyTable()));
             if (myProxy->getProxyField() != "") parameters->setValue("field", QString::fromStdString(myProxy->getProxyField()));
             if (myProxy->getGridName() != "") parameters->setValue("raster", getRelativePath(QString::fromStdString(myProxy->getGridName())));
+            if (myProxy->getStdDevThreshold() != NODATA) parameters->setValue("stddev_threshold", myProxy->getStdDevThreshold());
         parameters->endGroup();
     }
 }