Merge commit '866fa25539f662168042b38177cb217674f513ce'

agrolib/interpolation/interpolation.cpp

@@ -1872,6 +1872,7 @@ bool multipleDetrendingOtherProxiesFitting(int elevationPos, std::vector <Crit3D
     }
 
     // proxy spatial variability (1st step)
+    // this is done before check of incomplete to keep as many points as possible
     unsigned validNr;
     validNr = 0;
 
@@ -1947,6 +1948,7 @@ bool multipleDetrendingOtherProxiesFitting(int elevationPos, std::vector <Crit3D
     }
 
     // proxy spatial variability (2nd step)
+    // to be done because we might have excluded some points
     validNr = 0;
     for (int pos=0; pos < proxyNr; pos++)
     {
@@ -1972,20 +1974,20 @@ bool multipleDetrendingOtherProxiesFitting(int elevationPos, std::vector <Crit3D
     std::vector <double> predictands;
     std::vector <double> weights;
 
-    for (i=0; i < myPoints.size(); i++)
+    for (i=0; i < othersPoints.size(); i++)
     {
         rowPredictors.clear();
         for (int pos=0; pos < proxyNr; pos++)
-            if (pos != elevationPos && myCombination.isProxyActive(pos) && mySettings->getCurrentCombination().isProxySignificant(pos) && checkLapseRateCode(myPoints[i].lapseRateCode, mySettings->getUseLapseRateCode(), false))
+            if (pos != elevationPos && myCombination.isProxyActive(pos) && mySettings->getCurrentCombination().isProxySignificant(pos))
             {
-                proxyValue = myPoints[i].getProxyValue(pos);
+                proxyValue = othersPoints[i].getProxyValue(pos);
                 rowPredictors.push_back(proxyValue);
             }
 
         predictors.push_back(rowPredictors);
-        predictands.push_back(myPoints[i].value);
-        if (!isEqual(myPoints[i].regressionWeight, NODATA))
-            weights.push_back(myPoints[i].regressionWeight);
+        predictands.push_back(othersPoints[i].value);
+        if (!isEqual(othersPoints[i].regressionWeight, NODATA))
+            weights.push_back(othersPoints[i].regressionWeight);
         else
             weights.push_back(1);
     }
@@ -2252,8 +2254,7 @@ void optimalDetrending(meteoVariable myVar, Crit3DMeteoPoint* &myMeteoPoints, in
             if (mySettings->getUseTD() && getUseTdVar(myVar))
                 topographicDistanceOptimize(myVar, myMeteoPoints, nrMeteoPoints, interpolationPoints, mySettings, meteoSettings);
 
-            // TODO CATE : serve un check nei settings per excludeOutsideDem
-            if (computeResiduals(myVar, myMeteoPoints, nrMeteoPoints, interpolationPoints, mySettings, meteoSettings, false, true))
+            if (computeResiduals(myVar, myMeteoPoints, nrMeteoPoints, interpolationPoints, mySettings, meteoSettings, mySettings->getUseExcludeStationsOutsideDEM(), true))
             {
                 avgError = computeErrorCrossValidation(myMeteoPoints, nrMeteoPoints);
                 if (! isEqual(avgError, NODATA) && (isEqual(minError, NODATA) || avgError < minError))

agrolib/interpolation/interpolationSettings.cpp

@@ -455,6 +455,7 @@ void Crit3DInterpolationSettings::initialize()
     useTD = false;
     useLocalDetrending = false;
 	useGlocalDetrending = false;
+    useExcludeStationsOutsideDEM = false;
     topoDist_maxKh = 128;
     useDewPoint = true;
     useInterpolatedTForRH = true;
@@ -568,6 +569,9 @@ void Crit3DInterpolationSettings::setUseThermalInversion(bool myValue)
     selectedCombination.setUseThermalInversion(myValue);
 }
 
+void Crit3DInterpolationSettings::setUseExcludeStationsOutsideDEM(bool myValue)
+{ useExcludeStationsOutsideDEM = myValue; }
+
 void Crit3DInterpolationSettings::setUseTD(bool myValue)
 { useTD = myValue;}
 
@@ -589,6 +593,9 @@ void Crit3DInterpolationSettings::setUseDewPoint(bool myValue)
 bool Crit3DInterpolationSettings::getUseThermalInversion()
 { return (useThermalInversion);}
 
+bool Crit3DInterpolationSettings::getUseExcludeStationsOutsideDEM()
+{ return (useExcludeStationsOutsideDEM); }
+
 bool Crit3DInterpolationSettings::getUseDewPoint()
 { return (useDewPoint);}
 

agrolib/interpolation/interpolationSettings.h

@@ -173,6 +173,7 @@
 
         float minRegressionR2;
         bool useThermalInversion;
+        bool useExcludeStationsOutsideDEM;
         bool useTD;
         bool useLocalDetrending;
 		bool useGlocalDetrending;
@@ -231,6 +232,9 @@
         void setUseThermalInversion(bool myValue);
         bool getUseThermalInversion();
 
+        bool getUseExcludeStationsOutsideDEM();
+        void setUseExcludeStationsOutsideDEM(bool myValue);
+
         void setUseTD(bool myValue);
         bool getUseTD();
 

agrolib/project/dialogInterpolation.cpp

@@ -111,6 +111,10 @@ DialogInterpolation::DialogInterpolation(Project *myProject)
     thermalInversionEdit->setChecked(_interpolationSettings->getUseThermalInversion());
     layoutDetrending->addWidget(thermalInversionEdit);
 
+    excludeStationsOutsideDEM = new QCheckBox(tr("exclude meteo stations outside DEM"));
+    excludeStationsOutsideDEM->setChecked(_interpolationSettings->getUseExcludeStationsOutsideDEM());
+    layoutDetrending->addWidget(excludeStationsOutsideDEM);
+
     optimalDetrendingEdit = new QCheckBox(tr("optimal detrending"));
     optimalDetrendingEdit->setChecked(_interpolationSettings->getUseBestDetrending());
     if (_interpolationSettings->getUseBestDetrending())
@@ -345,6 +349,7 @@ void DialogInterpolation::accept()
     _interpolationSettings->setUseDoNotRetrend(doNotRetrendEdit->isChecked());
     _interpolationSettings->setUseRetrendOnly(retrendOnlyEdit->isChecked());
     _interpolationSettings->setUseThermalInversion(thermalInversionEdit->isChecked());
+    _interpolationSettings->setUseExcludeStationsOutsideDEM(excludeStationsOutsideDEM->isChecked());
     _interpolationSettings->setUseDewPoint(useDewPointEdit->isChecked());
     _interpolationSettings->setUseInterpolatedTForRH((useInterpolTForRH->isChecked()));
     _interpolationSettings->setMinRegressionR2(QLocale().toFloat(minRegressionR2Edit.text()));

agrolib/project/dialogInterpolation.h

@@ -24,6 +24,7 @@
             QLineEdit maxTdMultiplierEdit;
             QLineEdit minPointsLocalDetrendingEdit;
             QCheckBox* lapseRateCodeEdit;
+            QCheckBox* excludeStationsOutsideDEM;
             QCheckBox* thermalInversionEdit;
             QCheckBox* optimalDetrendingEdit;
             QCheckBox* multipleDetrendingEdit;

agrolib/project/project.cpp

@@ -652,6 +652,9 @@ bool Project::loadParameters(QString parametersFileName)
                 qualityInterpolationSettings.setUseThermalInversion(parametersSettings->value("thermalInversion").toBool());
             }
 
+            if (parametersSettings->contains("excludeStationsOutsideDEM"))
+                interpolationSettings.setUseExcludeStationsOutsideDEM(parametersSettings->value("excludeStationsOutsideDEM").toBool());
+
             if (parametersSettings->contains("topographicDistance"))
                 interpolationSettings.setUseTD(parametersSettings->value("topographicDistance").toBool());
 
@@ -2182,7 +2185,7 @@ bool Project::loadTopographicDistanceMaps(bool onlyWithData, bool showInfo)
     return true;
 }
 
-bool Project::glocalWeightsMaps(float windowWidth)
+bool Project::writeGlocalWeightsMaps(float windowWidth)
 {
     //controlla se c'è mappa aree
 
@@ -2283,7 +2286,7 @@ bool Project::loadGlocalAreasMap()
 
     if (!QFile::exists(fileNameMap + ".flt"))
     {
-        errorString = "Could not find " + fileNameMap + " file";
+        errorString = "Could not find file:\n" + fileNameMap;
         return false;
     }
 
@@ -2335,10 +2338,7 @@ bool Project::loadGlocalStationsAndCells(bool isGrid)
 
     //assegnazione pesi a ogni cella di ogni area
     if (!loadGlocalWeightMaps(myAreas, isGrid))
-    {
-        errorString = "Error loading glocal weight maps";
         return false;
-    }
 
     interpolationSettings.setMacroAreas(myAreas);
 
@@ -2349,13 +2349,21 @@ bool Project::loadGlocalWeightMaps(std::vector<Crit3DMacroArea> &myAreas, bool i
 {
     QString mapsFolder = projectPath + PATH_GLOCAL;
     if (! QDir(mapsFolder).exists())
-        glocalWeightsMaps(10000);
+    {
+        errorString = "Glocal dir not found. Please create glocal weight maps";
+        return false;
+    }
+
+    if (QDir(mapsFolder).isEmpty(QDir::Files))
+    {
+        errorString = "Glocal dir empty. Please create glocal weight maps";
+        return false;
+    }
 
     std::string myError;
     gis::Crit3DRasterGrid* macroAreasGrid = new gis::Crit3DRasterGrid();
     std::string fileName = mapsFolder.toStdString() + "glocalWeight_";
 
-
     std::vector<float> areaCells;
     int nrCols, nrRows;
     double myX, myY;
@@ -2372,44 +2380,55 @@ bool Project::loadGlocalWeightMaps(std::vector<Crit3DMacroArea> &myAreas, bool i
         nrRows = meteoGridDbHandler->gridStructure().header().nrRows;
     }
 
+    unsigned nrAreasWithCells = 0;
+
     for (int i = 0; i < myAreas.size(); i++)
     {
         if (QFile::exists(QString::fromStdString(fileName) + QString::number(i) + ".flt"))
         {
-        if (gis::readEsriGrid(fileName + std::to_string(i), macroAreasGrid, myError))
-        {
-            for (int row = 0; row < nrRows; row++)
+            if (readEsriGrid(fileName + std::to_string(i), macroAreasGrid, myError))
             {
-                for (int col = 0; col < nrCols; col++)
+                for (int row = 0; row < nrRows; row++)
                 {
-                    if (isGrid)
+                    for (int col = 0; col < nrCols; col++)
                     {
-                        myX = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.utm.x;
-                        myY = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.utm.y;
-                    }
-                    else
-                        DEM.getXY(row, col, myX, myY);
+                        if (isGrid)
+                        {
+                            myX = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.utm.x;
+                            myY = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.utm.y;
+                        }
+                        else
+                            DEM.getXY(row, col, myX, myY);
 
-                    myValue = macroAreasGrid->getValueFromXY(myX, myY);
+                        myValue = macroAreasGrid->getValueFromXY(myX, myY);
 
-                    if (!isEqual(myValue, NODATA) && !isEqual(myValue, 0))
-                    {
-                        areaCells.push_back(row*nrCols+col);
-                        areaCells.push_back(myValue);
+                        if (!isEqual(myValue, NODATA) && !isEqual(myValue, 0))
+                        {
+                            areaCells.push_back(row*nrCols+col);
+                            areaCells.push_back(myValue);
+                        }
                     }
                 }
+                if (areaCells.size() > 0) nrAreasWithCells++;
+                myAreas[i].setAreaCells(areaCells);
+                areaCells.clear();
+            }
+            else
+            {
+                macroAreasGrid->clear();
+                errorString = "macroareas not found";
+                return false;
             }
-            myAreas[i].setAreaCells(areaCells);
-            areaCells.clear();
-        }
-        else {
-            macroAreasGrid->clear();
-            return false;
-        }
         }
     }
     macroAreasGrid->clear();
 
+    if (nrAreasWithCells == 0)
+    {
+        errorString = "No valid files found in glocal dir. Please create glocal weight maps";
+        return false;
+    }
+
     return true;
 }
 
@@ -2713,7 +2732,7 @@ bool Project::interpolationCv(meteoVariable myVar, const Crit3DTime& myTime)
 
     if (! interpolationSettings.getUseLocalDetrending())
     {
-        if (! computeResiduals(myVar, meteoPoints, nrMeteoPoints, interpolationPoints, &interpolationSettings, meteoSettings, true, true))
+        if (! computeResiduals(myVar, meteoPoints, nrMeteoPoints, interpolationPoints, &interpolationSettings, meteoSettings, interpolationSettings.getUseExcludeStationsOutsideDEM(), true))
             return false;
     }
     else
@@ -3744,6 +3763,7 @@ void Project::saveInterpolationParameters()
         parametersSettings->setValue("lapseRateCode", interpolationSettings.getUseLapseRateCode());
         parametersSettings->setValue("meteogrid_upscalefromdem", interpolationSettings.getMeteoGridUpscaleFromDem());
         parametersSettings->setValue("thermalInversion", interpolationSettings.getUseThermalInversion());
+        parametersSettings->setValue("excludeStationsOutsideDEM", interpolationSettings.getUseExcludeStationsOutsideDEM());
         parametersSettings->setValue("topographicDistance", interpolationSettings.getUseTD());
         parametersSettings->setValue("localDetrending", interpolationSettings.getUseLocalDetrending());
         parametersSettings->setValue("topographicDistanceMaxMultiplier", QString::number(interpolationSettings.getTopoDist_maxKh()));
@@ -4426,8 +4446,11 @@ void Project::showLocalProxyGraph(gis::Crit3DGeoPoint myPoint)
 
     if (interpolationSettings.getUseGlocalDetrending() && ! interpolationSettings.isGlocalReady())
     {
-        if (loadGlocalAreasMap())
-            loadGlocalStationsAndCells(false);
+        if (! loadGlocalAreasMap() || ! loadGlocalStationsAndCells(false))
+        {
+            logError();
+            return;
+        }
     }
 
     localProxyWidget = new Crit3DLocalProxyWidget(myUtm.x, myUtm.y, myZDEM, myZGrid, this->gisSettings, &interpolationSettings, meteoPoints, nrMeteoPoints, currentVariable, currentFrequency, currentDate, currentHour, quality, &qualityInterpolationSettings, meteoSettings, &climateParameters, checkSpatialQuality);

agrolib/project/project.h

@@ -268,7 +268,7 @@
         bool loadGlocalWeightMaps(std::vector<Crit3DMacroArea> &myAreas, bool isGrid);
         bool loadGlocalStationsCsv(QString fileName, std::vector<std::vector<std::string> > &areaPoints);
         bool groupCellsInArea(std::vector<int> &areaPoints, unsigned index, bool isGrid);
-        bool glocalWeightsMaps(float windowWidth);
+        bool writeGlocalWeightsMaps(float windowWidth);
 		void macroAreaDetrending(Crit3DMacroArea myArea, meteoVariable myVar, std::vector<Crit3DInterpolationDataPoint> interpolationPoints, std::vector<Crit3DInterpolationDataPoint> &subsetInterpolationPoints, int elevationPos);
 
 

agrolib/proxyWidget/localProxyWidget.cpp

@@ -357,44 +357,44 @@ void Crit3DLocalProxyWidget::plot()
         areaCode = gis::getValueFromXY(*interpolationSettings->getMacroAreasMap(), x, y);
         if (areaCode < interpolationSettings->getMacroAreas().size())
         {
-        Crit3DMacroArea myArea = interpolationSettings->getMacroAreas()[areaCode];
-        std::vector<int> stations = myArea.getMeteoPoints();
-        if (detrended.isChecked())
-        {
-            outInterpolationPoints.clear();
-
-            checkAndPassDataToInterpolation(quality, myVar, meteoPoints, nrMeteoPoints, getCurrentTime(), SQinterpolationSettings,
-                                            interpolationSettings, meteoSettings, climateParam,
-                                            outInterpolationPoints, checkSpatialQuality, errorStdStr);
-
-            for (int k = 0; k < stations.size(); k++)
+            Crit3DMacroArea myArea = interpolationSettings->getMacroAreas()[areaCode];
+            std::vector<int> stations = myArea.getMeteoPoints();
+            if (detrended.isChecked())
             {
-                for (int j = 0; j < outInterpolationPoints.size(); j++)
-                    if (outInterpolationPoints[j].index == stations[k])
-                    {
-                        subsetInterpolationPoints.push_back(outInterpolationPoints[j]);
-                    }
-            }
+                outInterpolationPoints.clear();
 
-            detrending(subsetInterpolationPoints, interpolationSettings->getSelectedCombination(), interpolationSettings, climateParam, myVar, getCurrentTime());
-        }
-        else
-        {
-            outInterpolationPoints.clear();
-            checkAndPassDataToInterpolation(quality, myVar, meteoPoints, nrMeteoPoints, getCurrentTime(), SQinterpolationSettings,
-                                            interpolationSettings, meteoSettings, climateParam,
-                                            outInterpolationPoints, checkSpatialQuality, errorStdStr);
+                checkAndPassDataToInterpolation(quality, myVar, meteoPoints, nrMeteoPoints, getCurrentTime(), SQinterpolationSettings,
+                                                interpolationSettings, meteoSettings, climateParam,
+                                                outInterpolationPoints, checkSpatialQuality, errorStdStr);
 
-            for (int k = 0; k < stations.size(); k++)
+                for (int k = 0; k < stations.size(); k++)
+                {
+                    for (int j = 0; j < outInterpolationPoints.size(); j++)
+                        if (outInterpolationPoints[j].index == stations[k])
+                        {
+                            subsetInterpolationPoints.push_back(outInterpolationPoints[j]);
+                        }
+                }
+
+                detrending(subsetInterpolationPoints, interpolationSettings->getSelectedCombination(), interpolationSettings, climateParam, myVar, getCurrentTime());
+            }
+            else
             {
-                for (int j = 0; j < outInterpolationPoints.size(); j++)
-                    if (outInterpolationPoints[j].index == stations[k])
-                    {
-                        subsetInterpolationPoints.push_back(outInterpolationPoints[j]);
-                    }
+                outInterpolationPoints.clear();
+                checkAndPassDataToInterpolation(quality, myVar, meteoPoints, nrMeteoPoints, getCurrentTime(), SQinterpolationSettings,
+                                                interpolationSettings, meteoSettings, climateParam,
+                                                outInterpolationPoints, checkSpatialQuality, errorStdStr);
+
+                for (int k = 0; k < stations.size(); k++)
+                {
+                    for (int j = 0; j < outInterpolationPoints.size(); j++)
+                        if (outInterpolationPoints[j].index == stations[k])
+                        {
+                            subsetInterpolationPoints.push_back(outInterpolationPoints[j]);
+                        }
+                }
             }
         }
-        }
     }
     QList<QPointF> pointListPrimary, pointListSecondary, pointListSupplemental, pointListMarked;
     QMap< QString, QPointF > idPointMap1;