diff --git a/agrolib/criteriaModel/criteria1DCase.cpp b/agrolib/criteriaModel/criteria1DCase.cpp index 4f0db7f2..b11707e1 100644 --- a/agrolib/criteriaModel/criteria1DCase.cpp +++ b/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) diff --git a/agrolib/criteriaModel/criteria1DCase.h b/agrolib/criteriaModel/criteria1DCase.h index 2a2bd033..b561ec6e 100644 --- a/agrolib/criteriaModel/criteria1DCase.h +++ b/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(); diff --git a/agrolib/criteriaModel/criteria1DProject.cpp b/agrolib/criteriaModel/criteria1DProject.cpp index d76d51a0..8148c6ba 100644 --- a/agrolib/criteriaModel/criteria1DProject.cpp +++ b/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++) { diff --git a/agrolib/criteriaModel/criteria1DProject.h b/agrolib/criteriaModel/criteria1DProject.h index 426cb494..29c5775e 100644 --- a/agrolib/criteriaModel/criteria1DProject.h +++ b/agrolib/criteriaModel/criteria1DProject.h @@ -91,6 +91,7 @@ // specific output std::vector waterContentDepth; + std::vector degreeOfSaturationDepth; std::vector waterPotentialDepth; std::vector waterDeficitDepth; std::vector awcDepth; diff --git a/agrolib/crop/root.cpp b/agrolib/crop/root.cpp index 02ddc3f8..7001da48 100644 --- a/agrolib/crop/root.cpp +++ b/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) diff --git a/agrolib/gis/gis.cpp b/agrolib/gis/gis.cpp index b0beb0d3..b229f805 100644 --- a/agrolib/gis/gis.cpp +++ b/agrolib/gis/gis.cpp @@ -459,7 +459,7 @@ namespace gis } - void convertNodataRasterGrid(Crit3DRasterGrid& myGrid) + void convertFlagToNodata(Crit3DRasterGrid& myGrid) { if (myGrid.header->flag == NODATA) return; diff --git a/agrolib/gis/gis.h b/agrolib/gis/gis.h index 6ee692f9..94caea8d 100644 --- a/agrolib/gis/gis.h +++ b/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); diff --git a/agrolib/pragaProject/pragaProject.cpp b/agrolib/pragaProject/pragaProject.cpp index a11a76f0..de64abe5 100644 --- a/agrolib/pragaProject/pragaProject.cpp +++ b/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 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() { diff --git a/agrolib/pragaProject/pragaProject.h b/agrolib/pragaProject/pragaProject.h index 4d5db5cb..21cfffb0 100644 --- a/agrolib/pragaProject/pragaProject.h +++ b/agrolib/pragaProject/pragaProject.h @@ -142,6 +142,7 @@ bool monthlyAggregateVariablesGrid(const QDate &firstDate, const QDate &lastDate, QList &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); diff --git a/agrolib/project/project.cpp b/agrolib/project/project.cpp index f64e6f54..f9adbdab 100644 --- a/agrolib/project/project.cpp +++ b/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); diff --git a/agrolib/soil/soil.h b/agrolib/soil/soil.h index 8d8b66fb..24aaa2c1 100644 --- a/agrolib/soil/soil.h +++ b/agrolib/soil/soil.h @@ -152,6 +152,9 @@ Crit3DDriessen Driessen; Crit3DHorizon(); + + double getSoilFraction() + { return (1.0 - coarseFragments); } }; diff --git a/agrolib/soil/soilDbTools.cpp b/agrolib/soil/soilDbTools.cpp index ebddbcf3..7b72d803 100644 --- a/agrolib/soil/soilDbTools.cpp +++ b/agrolib/soil/soilDbTools.cpp @@ -10,15 +10,21 @@ #include #include #include +#include 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;