Merge commit 'fd0c80bef62b515a2e1305baed8589ad59e14b96'

agrolib/climate/climate.cpp

@@ -4991,7 +4991,7 @@ bool monthlyAggregateDataGrid(Crit3DMeteoGridDbHandler* meteoGridDbHandler, QDat
 
 int computeAnnualSeriesOnPointFromDaily(QString *myError, Crit3DMeteoPointsDbHandler* meteoPointsDbHandler, Crit3DMeteoGridDbHandler* meteoGridDbHandler,
                                          Crit3DMeteoPoint* meteoPointTemp, Crit3DClimate* clima, bool isMeteoGrid, bool isAnomaly, Crit3DMeteoSettings* meteoSettings,
-                                        std::vector<float> &outputValues, bool dataAlreadyLoaded)
+                                        std::vector<float> &outputValues, std::vector<int> &vectorYears, bool dataAlreadyLoaded)
 {
     int validYears = 0;
     if (clima->param1IsClimate())
@@ -5046,6 +5046,7 @@ int computeAnnualSeriesOnPointFromDaily(QString *myError, Crit3DMeteoPointsDbHan
             if ( elaborationOnPoint(myError, meteoPointsDbHandler, nullptr, meteoPointTemp, clima, isMeteoGrid, startDate, endDate, isAnomaly, meteoSettings, dataAlreadyLoaded))
             {
                 validYears = validYears + 1;
+                vectorYears.push_back(myYear);
                 if(isAnomaly)
                 {
                     outputValues.push_back(meteoPointTemp->anomaly);

agrolib/climate/climate.h

@@ -150,8 +150,8 @@
                     Crit3DQuality *qualityCheck, Crit3DClimateParameters *climateParam, QString &myError);
 
     int computeAnnualSeriesOnPointFromDaily(QString *myError, Crit3DMeteoPointsDbHandler* meteoPointsDbHandler, Crit3DMeteoGridDbHandler* meteoGridDbHandler,
-                    Crit3DMeteoPoint* meteoPointTemp, Crit3DClimate* clima, bool isMeteoGrid, bool isAnomaly, 
-					Crit3DMeteoSettings* meteoSettings, std::vector<float> &outputValues, bool dataAlreadyLoaded);
+                    Crit3DMeteoPoint* meteoPointTemp, Crit3DClimate* clima, bool isMeteoGrid, bool isAnomaly,
+                    Crit3DMeteoSettings* meteoSettings, std::vector<float> &outputValues, std::vector<int> &vectorYears, bool dataAlreadyLoaded);
 
 	void computeClimateOnDailyData(Crit3DMeteoPoint meteoPoint, meteoVariable var, QDate firstDate, QDate lastDate,
                     int smooth, float* dataPresence, Crit3DQuality* qualityCheck, Crit3DClimateParameters* climateParam,

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/dbMeteoGrid/dbMeteoGrid.cpp

@@ -2233,6 +2233,96 @@ bool Crit3DMeteoGridDbHandler::loadGridMonthlyData(QString &myError, QString met
     return true;
 }
 
+bool Crit3DMeteoGridDbHandler::loadGridAllMonthlyData(QString &myError, QDate firstDate, QDate lastDate)
+{
+    myError = "";
+    QString table = "MonthlyData";
+
+    // set day to 1 to better comparison
+    firstDate.setDate(firstDate.year(), firstDate.month(), 1);
+    lastDate.setDate(lastDate.year(), lastDate.month(), 1);
+    int numberOfMonths = (lastDate.year()-firstDate.year())*12 + lastDate.month() - (firstDate.month()-1);
+
+    if (firstDate > _lastMonthlyDate || lastDate < _firstMonthlyDate)
+    {
+        return false;
+    }
+
+    QSqlQuery qry(_db);
+    QDate date;
+    unsigned row, col;
+    int year, month, varCode;
+    float value;
+    QString pointCode;
+    QString prevCode = "NODATA";
+    QString statement = QString("SELECT * FROM `%1` WHERE `PragaYear` BETWEEN %2 AND %3 ORDER BY `PragaYear`").arg(table).arg(firstDate.year()).arg(lastDate.year());
+    if( !qry.exec(statement) )
+    {
+        myError = qry.lastError().text();
+        return false;
+    }
+    else
+    {
+        while (qry.next())
+        {
+            if (!getValue(qry.value("PragaYear"), &year))
+            {
+                myError = "Missing PragaYear";
+                return false;
+            }
+
+            if (!getValue(qry.value("PragaMonth"), &month))
+            {
+                myError = "Missing PragaMonth";
+                return false;
+            }
+
+            date.setDate(year,month, 1);
+            if (date < firstDate || date > lastDate)
+            {
+                continue;
+            }
+
+            if (!getValue(qry.value("PointCode"), &pointCode))
+            {
+                myError = "Missing PointCode";
+                return false;
+            }
+
+            if (!getValue(qry.value("VariableCode"), &varCode))
+            {
+                myError = "Missing VariableCode";
+                return false;
+            }
+
+            if (!getValue(qry.value("Value"), &value))
+            {
+                myError = "Missing Value";
+            }
+
+            meteoVariable variable = getMonthlyVarEnum(varCode);
+            if (prevCode != pointCode) // new point
+            {
+                if (!_meteoGrid->findMeteoPointFromId(&row, &col, pointCode.toStdString()) )
+                {
+                    myError = "Missing MeteoPoint id";
+                    return false;
+                }
+
+                if (_meteoGrid->meteoPointPointer(row,col)->nrObsDataDaysM == 0)
+                {
+                    _meteoGrid->meteoPointPointer(row,col)->initializeObsDataM(numberOfMonths, firstDate.month(), firstDate.year());
+                }
+            }
+            if (! _meteoGrid->meteoPointPointer(row,col)->setMeteoPointValueM(getCrit3DDate(date), variable, value))
+                return false;
+            prevCode = pointCode;
+        }
+    }
+
+    return true;
+}
+
 
 std::vector<float> Crit3DMeteoGridDbHandler::loadGridDailyVar(QString *myError, QString meteoPoint,
                                     meteoVariable variable, QDate first, QDate last, QDate* firstDateDB)

agrolib/dbMeteoGrid/dbMeteoGrid.h

@@ -119,6 +119,7 @@
         bool loadGridHourlyDataFixedFields(QString &myError, QString meteoPoint, QDateTime first, QDateTime last);
         bool loadGridHourlyDataEnsemble(QString &myError, QString meteoPoint, int memberNr, QDateTime first, QDateTime last);
         bool loadGridMonthlyData(QString &myError, QString meteoPoint, QDate firstDate, QDate lastDate);
+        bool loadGridAllMonthlyData(QString &myError, QDate firstDate, QDate lastDate);
         bool loadGridMonthlySingleDate(QString &myError, const QString &meteoPoint, const QDate &myDate);
 
         std::vector<float> loadGridDailyVar(QString *myError, QString meteoPoint, meteoVariable variable, QDate first, QDate last, QDate *firstDateDB);

agrolib/drought/drought.cpp

@@ -224,6 +224,7 @@ bool Drought::computeSpiParameters()
     int nTot = 0;
     std::vector<float> mySums;
     std::vector<float> monthSeries;
+    float minPerc = meteoSettings->getMinimumPercentage();
 
     for (int j = indexStart; j<meteoPoint->nrObsDataDaysM; j++)
     {
@@ -249,7 +250,7 @@ bool Drought::computeSpiParameters()
                     break;
             }
         }
-        if ( (count / nTot) < (meteoSettings->getMinimumPercentage() / 100) )
+        if ( (float)count / nTot < (minPerc / 100) )
         {
             mySums[n] = NODATA;
         }
@@ -271,7 +272,7 @@ bool Drought::computeSpiParameters()
             }
         }
 
-        if (n / (mySums.size()/12) >= meteoSettings->getMinimumPercentage() / 100)
+        if ((float)n / (mySums.size()/12) >= minPerc / 100)
         {
             generalizedGammaFitting(monthSeries, n, &(currentGamma[myMonth-1].beta), &(currentGamma[myMonth-1].gamma),  &(currentGamma[myMonth-1].pzero));
         }
@@ -309,6 +310,7 @@ bool Drought::computeSpeiParameters()
     std::vector<float> mySums;
     std::vector<float> monthSeries;
     std::vector<float> pwm(3);
+    float minPerc = meteoSettings->getMinimumPercentage();
 
     for (int j = indexStart; j<meteoPoint->nrObsDataDaysM; j++)
     {
@@ -334,7 +336,7 @@ bool Drought::computeSpeiParameters()
                     break;
             }
         }
-        if ( (count / nTot) < (meteoSettings->getMinimumPercentage() / 100))
+        if ( (float)count / nTot < (minPerc / 100))
         {
             mySums[n] = NODATA;
         }
@@ -356,7 +358,7 @@ bool Drought::computeSpeiParameters()
             }
         }
 
-        if (n / (mySums.size()/12) >= meteoSettings->getMinimumPercentage() / 100)
+        if ((float)n / (mySums.size()/12) >= minPerc / 100)
         {
             // Sort values
             sorting::quicksortAscendingFloat(monthSeries, 0, monthSeries.size()-1);
@@ -389,6 +391,7 @@ bool Drought::computePercentileValuesCurrentDay()
     std::vector<float> myValues;
     int nValid = 0;
     int nTot = 0;
+    float minPerc = meteoSettings->getMinimumPercentage();
 
     for (int j = indexStart+myMonth-1; j < meteoPoint->nrObsDataDaysM ; j=j+12)
     {
@@ -411,7 +414,7 @@ bool Drought::computePercentileValuesCurrentDay()
     }
     if (nTot > 0)
     {
-        if (nValid/nTot >= meteoSettings->getMinimumPercentage() / 100)
+        if ((float)nValid/nTot >= minPerc / 100)
         {
             int index = (date.year - meteoPoint->obsDataM[0]._year)*12 + date.month -meteoPoint->obsDataM[0]._month; // starts from 0
             if (index < meteoPoint->nrObsDataDaysM)