Merge commit '83a9273a8b5d95e36d673bb5645471510100390c'

agrolib/criteriaModel/criteria1DCase.cpp

@@ -361,7 +361,7 @@ bool Crit1DCase::computeNumericalFluxes(const Crit3DDate &myDate, std::string &e
 }
 
 
-void Crit1DCase::saveWaterContent()
+void Crit1DCase::storeWaterContent()
 {
     prevWaterContent.clear();
     prevWaterContent.resize(soilLayers.size());
@@ -537,7 +537,7 @@ bool Crit1DCase::computeDailyModel(Crit3DDate &myDate, std::string &error)
 
     output.dailyIrrigation = 0;
     // Water fluxes (first computation)
-    saveWaterContent();
+    storeWaterContent();
     if (! computeWaterFluxes(myDate, error)) return false;
     // Irrigation
     double irrigation = checkIrrigationDemand(doy, prec, precTomorrow, output.dailyMaxTranspiration);
@@ -653,11 +653,11 @@ double Crit1DCase::getTotalWaterContent()
 
 
 /*!
- * \brief getWaterContent
+ * \brief getVolumetricWaterContent
  * \param computationDepth = computation soil depth  [cm]
- * \return volumetric water content at specific depth [-]
+ * \return volumetric water content (soil moisture) at specific depth [m3 m-3]
  */
-double Crit1DCase::getWaterContent(double computationDepth)
+double Crit1DCase::getVolumetricWaterContent(double computationDepth)
 {
     computationDepth /= 100;        // [cm] --> [m]
 
@@ -679,6 +679,34 @@ double Crit1DCase::getWaterContent(double computationDepth)
 }
 
 
+/*!
+ * \brief getDegreeOfSaturation
+ * \param computationDepth = computation soil depth  [cm]
+ * \return degree of saturation at specific depth [-]
+ */
+double Crit1DCase::getDegreeOfSaturation(double computationDepth)
+{
+    computationDepth /= 100;        // [cm] --> [m]
+    if (computationDepth <= 0 || computationDepth > mySoil.totalDepth)
+    {
+        return NODATA;
+    }
+
+    double upperDepth, lowerDepth;
+    for (unsigned int i = 1; i < soilLayers.size(); i++)
+    {
+        upperDepth = soilLayers[i].depth - soilLayers[i].thickness * 0.5;
+        lowerDepth = soilLayers[i].depth + soilLayers[i].thickness * 0.5;
+        if (computationDepth >= upperDepth && computationDepth <= lowerDepth)
+        {
+            return soilLayers[i].waterContent / soilLayers[i].SAT;
+        }
+    }
+
+    return NODATA;
+}
+
+
 /*!
  * \brief getWaterPotential
  * \param computationDepth = computation soil depth  [cm]
@@ -778,7 +806,7 @@ double Crit1DCase::getWaterDeficitSum(double computationDepth)
 
 /*!
  * \brief getWaterCapacity
- * \param computationDepth = computation soil depth  [cm]
+ * \param computationDepth = computation soil depth [cm]
  * \return sum of available water capacity (FC-WP) from zero to computationDepth (mm)
  */
 double Crit1DCase::getWaterCapacitySum(double computationDepth)

agrolib/criteriaModel/criteria1DCase.h

@@ -78,7 +78,8 @@
         bool initializeWaterContent(Crit3DDate myDate);
         bool computeDailyModel(Crit3DDate &myDate, std::string &error);
 
-        double getWaterContent(double computationDepth);
+        double getVolumetricWaterContent(double computationDepth);
+        double getDegreeOfSaturation(double computationDepth);
         double getWaterPotential(double computationDepth);
         double getFractionAW(double computationDepth);
         double getSlopeStability(double computationDepth);
@@ -99,7 +100,7 @@
         bool computeNumericalFluxes(const Crit3DDate &myDate, std::string &error);
         bool computeWaterFluxes(const Crit3DDate &myDate, std::string &error);
         double checkIrrigationDemand(int doy, double currentPrec, double nextPrec, double maxTranspiration);
-        void saveWaterContent();
+        void storeWaterContent();
         void restoreWaterContent();
         double getTotalWaterContent();
 

agrolib/criteriaModel/criteria1DProject.cpp

@@ -63,9 +63,11 @@ void Crit1DProject::initialize()
     lastSimulationDate = QDate(1800,1,1);
 
     outputString = "";
+
     // specific outputs
     waterDeficitDepth.clear();
     waterContentDepth.clear();
+    degreeOfSaturationDepth.clear();
     waterPotentialDepth.clear();
     availableWaterDepth.clear();
     fractionAvailableWaterDepth.clear();
@@ -257,6 +259,12 @@ bool Crit1DProject::readSettings()
             projectError = "Wrong water content depth in " + configFileName;
             return false;
         }
+        depthList = projectSettings->value("degreeOfSaturation").toStringList();
+        if (! setVariableDepth(depthList, degreeOfSaturationDepth))
+        {
+            projectError = "Wrong degree of saturation depth in " + configFileName;
+            return false;
+        }
         depthList = projectSettings->value("waterPotential").toStringList();
         if (! setVariableDepth(depthList, waterPotentialDepth))
         {
@@ -1544,7 +1552,12 @@ bool Crit1DProject::createOutputTable(QString &myError)
     // specific depth variables
     for (unsigned int i = 0; i < waterContentDepth.size(); i++)
     {
-        QString fieldName = "SWC_" + QString::number(waterContentDepth[i]);
+        QString fieldName = "VWC_" + QString::number(waterContentDepth[i]);
+        queryString += ", " + fieldName + " REAL";
+    }
+    for (unsigned int i = 0; i < degreeOfSaturationDepth.size(); i++)
+    {
+        QString fieldName = "DEGSAT_" + QString::number(degreeOfSaturationDepth[i]);
         queryString += ", " + fieldName + " REAL";
     }
     for (unsigned int i = 0; i < waterPotentialDepth.size(); i++)
@@ -1604,7 +1617,12 @@ void Crit1DProject::updateOutput(Crit3DDate myDate, bool isFirst)
         // specific depth variables
         for (unsigned int i = 0; i < waterContentDepth.size(); i++)
         {
-            QString fieldName = "SWC_" + QString::number(waterContentDepth[i]);
+            QString fieldName = "VWC_" + QString::number(waterContentDepth[i]);
+            outputString += ", " + fieldName;
+        }
+        for (unsigned int i = 0; i < degreeOfSaturationDepth.size(); i++)
+        {
+            QString fieldName = "DEGSAT_" + QString::number(degreeOfSaturationDepth[i]);
             outputString += ", " + fieldName;
         }
         for (unsigned int i = 0; i < waterPotentialDepth.size(); i++)
@@ -1669,7 +1687,11 @@ void Crit1DProject::updateOutput(Crit3DDate myDate, bool isFirst)
     // specific depth variables
     for (unsigned int i = 0; i < waterContentDepth.size(); i++)
     {
-        outputString += "," + QString::number(myCase.getWaterContent(waterContentDepth[i]), 'g', 4);
+        outputString += "," + QString::number(myCase.getVolumetricWaterContent(waterContentDepth[i]), 'g', 4);
+    }
+    for (unsigned int i = 0; i < degreeOfSaturationDepth.size(); i++)
+    {
+        outputString += "," + QString::number(myCase.getDegreeOfSaturation(degreeOfSaturationDepth[i]), 'g', 4);
     }
     for (unsigned int i = 0; i < waterPotentialDepth.size(); i++)
     {

agrolib/criteriaModel/criteria1DProject.h

@@ -91,6 +91,7 @@
 
         // specific output
         std::vector<int> waterContentDepth;
+        std::vector<int> degreeOfSaturationDepth;
         std::vector<int> waterPotentialDepth;
         std::vector<int> waterDeficitDepth;
         std::vector<int> awcDepth;

agrolib/crop/root.cpp

@@ -591,6 +591,7 @@ namespace root
             }
         }
 
+        // normalize root density
         if (rootDensitySumSubset != rootDensitySum)
         {
             double ratio = rootDensitySum / rootDensitySumSubset;
@@ -600,6 +601,7 @@ namespace root
             }
         }
 
+        // find first and last root layers
         myCrop->roots.firstRootLayer = 0;
         unsigned int layer = 0;
         while (layer < nrLayers && myCrop->roots.rootDensity[layer] == 0.0)

agrolib/gis/gis.cpp

@@ -459,7 +459,7 @@ namespace gis
     }
 
 
-    void convertNodataRasterGrid(Crit3DRasterGrid& myGrid)
+    void convertFlagToNodata(Crit3DRasterGrid& myGrid)
     {
         if (myGrid.header->flag == NODATA)
             return;

agrolib/gis/gis.h

@@ -197,7 +197,7 @@
         float computeDistance(float x1, float y1, float x2, float y2);
         double computeDistancePoint(Crit3DUtmPoint *p0, Crit3DUtmPoint *p1);
         bool updateMinMaxRasterGrid(Crit3DRasterGrid *rasterGrid);
-        void convertNodataRasterGrid(Crit3DRasterGrid& myGrid);
+        void convertFlagToNodata(Crit3DRasterGrid& myGrid);
         bool updateColorScale(Crit3DRasterGrid* rasterGrid, int row0, int col0, int row1, int col1);
 
         void getRowColFromXY(const Crit3DRasterHeader& myHeader, double myX, double myY, int *row, int *col);

agrolib/pragaProject/pragaProject.cpp

@@ -3460,6 +3460,115 @@ bool PragaProject::computeDroughtIndexPoint(droughtIndex index, int timescale, i
     return true;
 }
 
+bool PragaProject::computeDroughtIndexPointGUI(droughtIndex index, int timescale, int refYearStart, int refYearEnd, QDate myDate)
+{
+    // check meteo point
+    if (! meteoPointsLoaded)
+    {
+        logError("No meteo point");
+        return false;
+    }
+
+    // check ref years
+    if (refYearStart > refYearEnd)
+    {
+        logError("Wrong reference years");
+        return false;
+    }
+
+    QDate firstDate(refYearStart,1,1);
+    QDate lastDate;
+    int maxYear = std::max(refYearEnd,myDate.year());
+    if (maxYear == QDate::currentDate().year())
+    {
+        lastDate.setDate(maxYear, QDate::currentDate().month(),1);
+    }
+    else
+    {
+        lastDate.setDate(maxYear,12,1);
+    }
+
+    bool loadHourly = false;
+    bool loadDaily = true;
+    bool showInfo = true;
+    float value = NODATA;
+    QString indexStr;
+
+    if (index == INDEX_SPI)
+    {
+        indexStr = "SPI";
+    }
+    else if (index == INDEX_SPEI)
+    {
+        indexStr = "SPEI";
+    }
+    else if (index == INDEX_DECILES)
+    {
+        indexStr = "DECILES";
+    }
+    else
+    {
+        logError("Unknown index");
+        return false;
+    }
+
+    if (!loadMeteoPointsData(firstDate, lastDate, loadHourly, loadDaily, showInfo))
+    {
+        logError("There are no data");
+        return false;
+    }
+
+    int step = 0;
+    if (showInfo)
+    {
+        QString infoStr = "Compute drought - Meteo Points";
+        step = setProgressBar(infoStr, nrMeteoPoints);
+    }
+
+    std::vector<meteoVariable> dailyMeteoVar;
+    dailyMeteoVar.push_back(dailyPrecipitation);
+    dailyMeteoVar.push_back(dailyReferenceEvapotranspirationHS);
+    bool res = false;
+
+    int nrMonths = (lastDate.year()-firstDate.year())*12+lastDate.month()-(firstDate.month()-1);
+    for (int i=0; i < nrMeteoPoints; i++)
+    {
+        if (showInfo && (i % step) == 0)
+        {
+            updateProgressBar(i);
+        }
+
+        // compute monthly data
+        meteoPoints[i].initializeObsDataM(nrMonths, firstDate.month(), firstDate.year());
+        for(int j = 0; j < dailyMeteoVar.size(); j++)
+        {
+            meteoPoints[i].computeMonthlyAggregate(getCrit3DDate(firstDate), getCrit3DDate(lastDate), dailyMeteoVar[j], meteoSettings, quality, &climateParameters);
+        }
+        Drought mydrought(index, refYearStart, refYearEnd, getCrit3DDate(myDate), &(meteoPoints[i]), meteoSettings);
+        if (timescale > 0)
+        {
+            mydrought.setTimeScale(timescale);
+        }
+        if (index == INDEX_DECILES)
+        {
+            if (mydrought.computePercentileValuesCurrentDay())
+            {
+                value = mydrought.getCurrentPercentileValue();
+            }
+        }
+        else if (index == INDEX_SPI || index == INDEX_SPEI)
+        {
+            value = mydrought.computeDroughtIndex();
+        }
+        meteoPoints[i].elaboration = value;
+        if (value != NODATA)
+        {
+            res = true;
+        }
+    }
+    return res;
+}
+
 bool PragaProject::activeMeteoGridCellsWithDEM()
 {
 

agrolib/pragaProject/pragaProject.h

@@ -142,6 +142,7 @@
         bool monthlyAggregateVariablesGrid(const QDate &firstDate, const QDate &lastDate, QList <meteoVariable> &variablesList);
         bool computeDroughtIndexAll(droughtIndex index, int firstYear, int lastYear, QDate date, int timescale, meteoVariable myVar);
         bool computeDroughtIndexPoint(droughtIndex index, int timescale, int refYearStart, int refYearEnd);
+        bool computeDroughtIndexPointGUI(droughtIndex index, int timescale, int refYearStart, int refYearEnd, QDate myDate);
         void showPointStatisticsWidgetPoint(std::string idMeteoPoint);
         void showHomogeneityTestWidgetPoint(std::string idMeteoPoint);
         void showSynchronicityTestWidgetPoint(std::string idMeteoPoint);

agrolib/project/project.cpp

@@ -996,7 +996,7 @@ bool Project::loadDEM(QString myFileName)
         QString infoStr = "WARNING: " + QString::number(DEM.header->flag) + " is not a valid NODATA value for DEM!";
         infoStr += " It will be converted in: " + QString::number(NODATA);
         logInfo(infoStr);
-        gis::convertNodataRasterGrid(DEM);
+        gis::convertFlagToNodata(DEM);
     }
 
     setColorScale(noMeteoTerrain, DEM.colorScale);

agrolib/soil/soil.h

@@ -152,6 +152,9 @@
             Crit3DDriessen Driessen;
 
             Crit3DHorizon();
+
+            double getSoilFraction()
+                { return (1.0 - coarseFragments); }
         };
 
 

agrolib/soil/soilDbTools.cpp

@@ -10,15 +10,21 @@
 #include <QSqlError>
 #include <QUuid>
 #include <QVariant>
+#include <QFile>
 
 
 bool openDbSoil(const QString &dbSoilName, QSqlDatabase &dbSoil, QString &errorStr)
 {
+    if (! QFile(dbSoilName).exists())
+    {
+        errorStr = "Soil database doesn't exist:\n" + dbSoilName;
+        return false;
+    }
 
     dbSoil = QSqlDatabase::addDatabase("QSQLITE", QUuid::createUuid().toString());
     dbSoil.setDatabaseName(dbSoilName);
 
-    if (!dbSoil.open())
+    if (! dbSoil.open())
     {
        errorStr = "Connection with database fail";
        return false;