From a31d16a5b0ccc82941e3510bd4e5bfd3c135ac97 Mon Sep 17 00:00:00 2001 From: Fausto Tomei Date: Sun, 26 Nov 2023 19:22:11 +0100 Subject: [PATCH] update 3D --- agrolib/crop/crop.cpp | 60 +- agrolib/crop/crop.h | 2 +- agrolib/crop/cropDbTools.cpp | 6 +- agrolib/interpolation/interpolation.cpp | 9 +- agrolib/snow/snow.cpp | 3 +- agrolib/utilities/utilities.cpp | 7 + agrolib/utilities/utilities.h | 1 + bin/CRITERIA3D/criteria3DProject.cpp | 109 ++- bin/CRITERIA3D/criteria3DProject.h | 17 +- bin/CRITERIA3D/mainwindow.cpp | 19 +- bin/CRITERIA3D/shared/project3D.cpp | 1160 +++++++++++------------ bin/CRITERIA3D/shared/project3D.h | 18 +- 12 files changed, 693 insertions(+), 718 deletions(-) diff --git a/agrolib/crop/crop.cpp b/agrolib/crop/crop.cpp index 8b030d153..385e3c19a 100644 --- a/agrolib/crop/crop.cpp +++ b/agrolib/crop/crop.cpp @@ -51,14 +51,9 @@ void Crit3DCrop::clear() idCrop = ""; type = HERBACEOUS_ANNUAL; - /*! - * \brief roots - */ roots.clear(); - /*! - * \brief crop cycle - */ + // crop cycle sowingDoy = NODATA; currentSowingDoy = NODATA; doyStartSenescence = NODATA; @@ -74,17 +69,13 @@ void Crit3DCrop::clear() degreeDaysDecrease = NODATA; degreeDaysEmergence = NODATA; - /*! - * \brief water need - */ + // water need kcMax = NODATA; psiLeaf = NODATA; stressTolerance = NODATA; fRAW = NODATA; - /*! - * \brief irrigation - */ + // irrigation irrigationShift = NODATA; irrigationVolume = NODATA; degreeDaysStartIrrigation = NODATA; @@ -93,9 +84,7 @@ void Crit3DCrop::clear() doyEndIrrigation = NODATA; maxSurfacePuddle = NODATA; - /*! - * \brief variables - */ + // variables isLiving = false; isEmerged = false; LAIstartSenescence = NODATA; @@ -330,7 +319,7 @@ bool Crit3DCrop::isRootStatic() const /*! * \brief getSurfaceWaterPonding - * \return maximum height of water ponding [mm] + * \return maximum height of water pond [mm] */ double Crit3DCrop::getSurfaceWaterPonding() const { @@ -406,7 +395,7 @@ bool Crit3DCrop::needReset(Crit3DDate myDate, double latitude, double waterTable // reset of (already initialized) crop -// TODO: smart start (using sowing doy and cycle) +// TODO: smart start (using meteo settings) void Crit3DCrop::resetCrop(unsigned int nrLayers) { // roots @@ -452,7 +441,7 @@ bool Crit3DCrop::dailyUpdate(const Crit3DDate &myDate, double latitude, const st unsigned int nrLayers = unsigned(soilLayers.size()); - // check start/end crop cycle (update isLiving) + // check start/end crop cycle if (needReset(myDate, latitude, waterTableDepth)) { resetCrop(nrLayers); @@ -520,29 +509,26 @@ bool Crit3DCrop::restore(const Crit3DDate &myDate, double latitude, const std::v } -// Liangxia Zhang, Zhongmin Hu, Jiangwen Fan, Decheng Zhou & Fengpei Tang, 2014 -// A meta-analysis of the canopy light extinction coefficient in terrestrial ecosystems -// "Cropland had the highest value of K (0.62), followed by broadleaf forest (0.59), -// shrubland (0.56), grassland (0.50), and needleleaf forest (0.45)" -double Crit3DCrop::getSurfaceCoverFraction() +/*! \brief getCoveredSurfaceFraction + * \ref Liangxia Zhang, Zhongmin Hu, Jiangwen Fan, Decheng Zhou & Fengpei Tang, 2014 + * A meta-analysis of the canopy light extinction coefficient in terrestrial ecosystems + * "Cropland had the highest value of K (0.62), followed by broadleaf forest (0.59) + * shrubland (0.56), grassland (0.50), and needleleaf forest (0.45)" + * \return covered surface fraction [-] + */ +double Crit3DCrop::getCoveredSurfaceFraction() { - double k = 0.6; // [-] light extinction coefficient + if (idCrop == "" || ! isLiving || LAI < EPSILON) return 0; - if (idCrop == "" || ! isLiving || LAI < EPSILON) - { - return 0; - } - else - { - return 1 - exp(-k * LAI); - } + double k = 0.6; // [-] light extinction coefficient + return 1 - exp(-k * LAI); } double Crit3DCrop::getMaxEvaporation(double ET0) { - double evapMax = ET0 * (1.0 - getSurfaceCoverFraction()); - // TODO check + double evapMax = ET0 * (1.0 - getCoveredSurfaceFraction()); + // TODO check evaporation on free water return evapMax * 0.66; } @@ -552,9 +538,9 @@ double Crit3DCrop::getMaxTranspiration(double ET0) if (idCrop == "" || ! isLiving || LAI < EPSILON) return 0; - double SCF = this->getSurfaceCoverFraction(); - double kcmaxFactor = 1 + (kcMax - 1) * SCF; - return ET0 * (SCF * kcmaxFactor); + double coverSurfFraction = getCoveredSurfaceFraction(); + double kcFactor = 1 + (kcMax - 1) * coverSurfFraction; + return ET0 * coverSurfFraction * kcFactor; } diff --git a/agrolib/crop/crop.h b/agrolib/crop/crop.h index 35766f176..55a9b4aff 100644 --- a/agrolib/crop/crop.h +++ b/agrolib/crop/crop.h @@ -95,7 +95,7 @@ bool restore(const Crit3DDate &myDate, double latitude, const std::vector &soilLayers, double currentWaterTable, std::string &myError); - double getSurfaceCoverFraction(); + double getCoveredSurfaceFraction(); double getMaxEvaporation(double ET0); double getMaxTranspiration(double ET0); double getSurfaceWaterPonding() const; diff --git a/agrolib/crop/cropDbTools.cpp b/agrolib/crop/cropDbTools.cpp index c7358a32e..6904172b7 100644 --- a/agrolib/crop/cropDbTools.cpp +++ b/agrolib/crop/cropDbTools.cpp @@ -67,7 +67,11 @@ bool loadCropParameters(const QSqlDatabase &dbCrop, QString idCrop, Crit3DCrop & myCrop.roots.rootDepthMax = query.value("root_depth_max").toDouble(); myCrop.roots.actualRootDepthMax = myCrop.roots.rootDepthMax; - if (! getValue(query.value("roots_additional_cohesion"), &(myCrop.roots.rootsAdditionalCohesion))) + if (fieldExists(query, "roots_additional_cohesion")) + { + getValue(query.value("roots_additional_cohesion"), &(myCrop.roots.rootsAdditionalCohesion)); + } + else { // default: no mechanical effect of roots myCrop.roots.rootsAdditionalCohesion = 0; diff --git a/agrolib/interpolation/interpolation.cpp b/agrolib/interpolation/interpolation.cpp index 74262544e..d76526227 100644 --- a/agrolib/interpolation/interpolation.cpp +++ b/agrolib/interpolation/interpolation.cpp @@ -23,7 +23,6 @@ ftomei@arpae.it */ -#include #include #include #include @@ -1172,7 +1171,7 @@ float retrend(meteoVariable myVar, vector myProxyValues, Crit3DInterpola if (! getUseDetrendingVar(myVar)) return 0.; - float retrendValue = 0.; + double retrendValue = 0.; double myProxyValue; Crit3DProxy* myProxy = nullptr; Crit3DProxyCombination myCombination = mySettings->getCurrentCombination(); @@ -1223,7 +1222,7 @@ float retrend(meteoVariable myVar, vector myProxyValues, Crit3DInterpola } } - return retrendValue; + return float(retrendValue); } @@ -1395,6 +1394,7 @@ std::vector getfittingParameters(Crit3DProxyCombination myCombination, return myParam; } + bool multipleDetrending(std::vector &myPoints, Crit3DInterpolationSettings* mySettings, meteoVariable myVar) { if (! getUseDetrendingVar(myVar)) return false; @@ -1519,7 +1519,7 @@ bool multipleDetrending(std::vector &myPoints, Cr setFittingParameters(myCombination, mySettings, myFunc, parametersMin, parametersMax, parametersDelta, parameters); // multiple non linear fitting - int nSteps = interpolation::bestFittingMarquardt_nDimension(&functionSum, myFunc, 10000, 5, parametersMin, parametersMax, parameters, parametersDelta, + interpolation::bestFittingMarquardt_nDimension(&functionSum, myFunc, 10000, 5, parametersMin, parametersMax, parameters, parametersDelta, 100, EPSILON, 0.01, predictors, predictands, false, weights); mySettings->setFittingFunction(myFunc); @@ -1547,6 +1547,7 @@ bool multipleDetrending(std::vector &myPoints, Cr return true; } + void topographicDistanceOptimize(meteoVariable myVar, Crit3DMeteoPoint* &myMeteoPoints, int nrMeteoPoints, diff --git a/agrolib/snow/snow.cpp b/agrolib/snow/snow.cpp index eb5d46301..85e050705 100644 --- a/agrolib/snow/snow.cpp +++ b/agrolib/snow/snow.cpp @@ -392,7 +392,8 @@ void Crit3DSnow::computeSnowBrooksModel() double freeze_melt = 0; // [mm] freeze (positive) or melt (negative) // pag.53 (3.25) - /*! net amount of liquid water that freezes (heat added to the snow pack) and ice that melts (heat removed from the snow pack) */ + /*! net amount of liquid water that freezes (heat added to the snow pack) + * and ice that melts (heat removed from the snow pack) */ double w = (prevInternalEnergy + QTotal) / (LATENT_HEAT_FUSION * WATER_DENSITY); // [m] if (w < 0) { diff --git a/agrolib/utilities/utilities.cpp b/agrolib/utilities/utilities.cpp index 3f100d208..78b82876b 100644 --- a/agrolib/utilities/utilities.cpp +++ b/agrolib/utilities/utilities.cpp @@ -46,6 +46,13 @@ QList getFieldsUpperCase(const QSqlQuery& query) } +bool fieldExists(const QSqlQuery &query, const QString fieldName) +{ + QList fieldList = getFieldsUpperCase(query); + return fieldList.contains(fieldName.toUpper()); +} + + // return boolean (false if recordset is not valid) bool getValue(QVariant myRs) { diff --git a/agrolib/utilities/utilities.h b/agrolib/utilities/utilities.h index 134fa0042..47a5c88f0 100644 --- a/agrolib/utilities/utilities.h +++ b/agrolib/utilities/utilities.h @@ -30,6 +30,7 @@ QList getFields(QSqlDatabase* db_, QString tableName); QList getFields(const QSqlQuery& query); QList getFieldsUpperCase(const QSqlQuery &query); + bool fieldExists(const QSqlQuery &query, const QString fieldName); bool getValue(QVariant myRs); bool getValue(QVariant myRs, int* myValue); diff --git a/bin/CRITERIA3D/criteria3DProject.cpp b/bin/CRITERIA3D/criteria3DProject.cpp index 97dd3e8f3..0eb187a30 100644 --- a/bin/CRITERIA3D/criteria3DProject.cpp +++ b/bin/CRITERIA3D/criteria3DProject.cpp @@ -39,27 +39,6 @@ #include -Crit3DProcesses::Crit3DProcesses() -{ - initialize(); -} - - -void Crit3DProcesses::initialize() -{ - computeMeteo = false; - computeRadiation = false; - computeWater = false; - computeEvaporation = false; - computeCrop = false; - computeSnow = false; - computeSolutes = false; - computeHeat = false; - computeAdvectiveHeat = false; - computeLatentHeat = false; -} - - Crit3DProject::Crit3DProject() : Project3D() { _saveOutputRaster = false; @@ -90,14 +69,11 @@ bool Crit3DProject::initializeCriteria3DModel() if (! initializeWaterBalance3D()) { clearWaterBalance3D(); - logError("Criteria3D model NOT initialized."); + logError("Criteria3D model is NOT initialized."); return false; } - if (processes.computeCrop) - { - initializeCrop(); - } + initializeCrop(); isCriteria3DInitialized = true; logInfoGUI("Criteria3D model initialized"); @@ -108,12 +84,19 @@ bool Crit3DProject::initializeCriteria3DModel() void Crit3DProject::initializeCrop() { + // initialize LAI and degree days map to NODATA + laiMap.initializeGrid(*(DEM.header)); + degreeDaysMap.initializeGrid(*(DEM.header)); + + if (! processes.computeCrop) + { + // nothing to do + return; + } + dailyTminMap.initializeGrid(*(DEM.header)); dailyTmaxMap.initializeGrid(*(DEM.header)); - degreeDaysMap.initializeGrid(*(DEM.header)); - laiMap.initializeGrid(*(DEM.header)); - for (int row = 0; row < DEM.header->nrRows; row++) { for (int col = 0; col < DEM.header->nrCols; col++) @@ -229,10 +212,10 @@ void Crit3DProject::dailyUpdateCrop() /*! - * \brief computeRealET + * \brief assignETreal * assign soil evaporation and crop transpiration for the whole domain */ -void Crit3DProject::computeRealET() +void Crit3DProject::assignETreal() { totalEvaporation = 0; totalTranspiration = 0; @@ -254,13 +237,45 @@ void Crit3DProject::computeRealET() // assign real evaporation double realEvap = assignEvaporation(row, col, lai); // [mm] - double flow = area * (realEvap / 1000.); // [m3 h-1] - totalEvaporation += flow; + double evapFlow = area * (realEvap / 1000.); // [m3 h-1] + totalEvaporation += evapFlow; // assign real transpiration - double realTransp = assignTranspiration(row, col, lai); // [mm] - flow = area * (realTransp / 1000.); // [m3 h-1] - totalTranspiration += flow; + double realTransp = assignTranspiration(row, col, lai); // [mm] + double traspFlow = area * (realTransp / 1000.); // [m3 h-1] + totalTranspiration += traspFlow; + } + } + } +} + + +void Crit3DProject::assignPrecipitation() +{ + // initialize + totalPrecipitation = 0; + + double area = DEM.header->cellSize * DEM.header->cellSize; + + // precipitation + for (long row = 0; row < indexMap.at(0).header->nrRows; row++) + { + for (long col = 0; col < indexMap.at(0).header->nrCols; col++) + { + int surfaceIndex = indexMap.at(0).value[row][col]; + if (surfaceIndex != indexMap.at(0).header->flag) + { + double waterSource = 0; + float prec = hourlyMeteoMaps->mapHourlyPrec->value[row][col]; + if (! isEqual(prec, hourlyMeteoMaps->mapHourlyPrec->header->flag)) + waterSource += prec; + + if (waterSource > 0) + { + double flow = area * (waterSource / 1000.); // [m3 h-1] + waterSinkSource[unsigned(surfaceIndex)] += flow / 3600.; // [m3 s-1] + totalPrecipitation += flow; + } } } } @@ -874,7 +889,7 @@ bool Crit3DProject::computeSnowModel() if (! snowMaps.isInitialized) { - if (! this->initializeSnowModel()) + if (! initializeSnowModel()) return false; } @@ -891,20 +906,20 @@ bool Crit3DProject::computeSnowModel() DEM.getRowCol(x, y, row, col); if (! gis::isOutOfGridRowCol(row, col, DEM)) { - this->computeSnowPoint(row, col); + computeSnowPoint(row, col); } } } } else { - for (long row = 0; row < DEM.header->nrRows; row++) + for (int row = 0; row < DEM.header->nrRows; row++) { - for (long col = 0; col < DEM.header->nrCols; col++) + for (int col = 0; col < DEM.header->nrCols; col++) { if (! isEqual(DEM.value[row][col], DEM.header->flag)) { - this->computeSnowPoint(row, col); + computeSnowPoint(row, col); } } } @@ -989,7 +1004,7 @@ bool Crit3DProject::modelHourlyCycle(QDateTime myTime, const QString& hourlyOutp if (processes.computeSnow) { - // check conflitti con evaporation + // check evaporation // check snowmelt -> surface H0 if (! computeSnowModel()) { @@ -1014,19 +1029,15 @@ bool Crit3DProject::modelHourlyCycle(QDateTime myTime, const QString& hourlyOutp { saveHourlyMeteoOutput(referenceEvapotranspiration, hourlyOutputPath, myTime); } - } - if (processes.computeEvaporation && (! processes.computeCrop)) - { - computeBareSoilEvaporation(); + assignETreal(); + + qApp->processEvents(); } if (processes.computeCrop) { updateDailyTemperatures(); - - computeRealET(); - qApp->processEvents(); } diff --git a/bin/CRITERIA3D/criteria3DProject.h b/bin/CRITERIA3D/criteria3DProject.h index 6bafbf728..bd14b488c 100644 --- a/bin/CRITERIA3D/criteria3DProject.h +++ b/bin/CRITERIA3D/criteria3DProject.h @@ -23,19 +23,6 @@ #include - class Crit3DProcesses - { - public: - - bool computeMeteo, computeRadiation, computeWater; - bool computeEvaporation, computeCrop, computeSnow, computeSolutes; - bool computeHeat, computeAdvectiveHeat, computeLatentHeat; - - Crit3DProcesses(); - void initialize(); - }; - - class Crit3DProject : public Project3D { @@ -59,7 +46,6 @@ gis::Crit3DRasterGrid dailyTmaxMap; Crit3DSnow snowModel; - Crit3DProcesses processes; bool modelPause, modelStop; @@ -71,7 +57,8 @@ bool initializeCriteria3DModel(); void initializeCrop(); void dailyUpdateCrop(); - void computeRealET(); + void assignETreal(); + void assignPrecipitation(); bool runModels(QDateTime firstTime, QDateTime lastTime); diff --git a/bin/CRITERIA3D/mainwindow.cpp b/bin/CRITERIA3D/mainwindow.cpp index 87d0e4afb..f6c565c7a 100644 --- a/bin/CRITERIA3D/mainwindow.cpp +++ b/bin/CRITERIA3D/mainwindow.cpp @@ -2175,6 +2175,11 @@ void MainWindow::on_actionWaterFluxes_settings_triggered() void MainWindow::on_actionCriteria3D_Initialize_triggered() { + myProject.processes.initialize(); + myProject.processes.computeMeteo = true; + myProject.processes.computeRadiation = true; + myProject.processes.computeWater = true; + myProject.processes.computeEvaporation = true; myProject.processes.computeCrop = true; if (myProject.initializeCriteria3DModel()) @@ -2214,13 +2219,6 @@ void MainWindow::on_actionCriteria3D_compute_current_hour_triggered() QDateTime currentTime = myProject.getCurrentTime(); - myProject.processes.initialize(); - myProject.processes.computeMeteo = true; - myProject.processes.computeRadiation = true; - myProject.processes.computeWater = true; - myProject.processes.computeEvaporation = true; - myProject.processes.computeCrop = true; - startModels(currentTime, currentTime); } @@ -2237,13 +2235,6 @@ void MainWindow::on_actionCriteria3D_run_models_triggered() if (! selectDates(firstTime, lastTime)) return; - myProject.processes.initialize(); - myProject.processes.computeMeteo = true; - myProject.processes.computeRadiation = true; - myProject.processes.computeWater = false; - myProject.processes.computeEvaporation = true; - myProject.processes.computeCrop = true; - startModels(firstTime, lastTime); } diff --git a/bin/CRITERIA3D/shared/project3D.cpp b/bin/CRITERIA3D/shared/project3D.cpp index 163ab2380..f20eeb487 100644 --- a/bin/CRITERIA3D/shared/project3D.cpp +++ b/bin/CRITERIA3D/shared/project3D.cpp @@ -24,7 +24,6 @@ */ #include "commonConstants.h" -#include "basicMath.h" #include "cropDbTools.h" #include "project3D.h" #include "soilFluxes3D.h" @@ -57,6 +56,27 @@ void WaterFluxesParameters::initialize() } +Crit3DProcesses::Crit3DProcesses() +{ + initialize(); +} + + +void Crit3DProcesses::initialize() +{ + computeMeteo = false; + computeRadiation = false; + computeWater = false; + computeEvaporation = false; + computeCrop = false; + computeSnow = false; + computeSolutes = false; + computeHeat = false; + computeAdvectiveHeat = false; + computeLatentHeat = false; +} + + Project3D::Project3D() : Project() { initializeProject3D(); @@ -174,21 +194,28 @@ bool Project3D::loadProject3DSettings() bool Project3D::initializeWaterBalance3D() { - logInfo("\nInitialize Waterbalance..."); + logInfo("\nInitialize 3D water balance..."); + // check soil if (!soilMap.isLoaded || soilList.size() == 0) { logInfo("WARNING: soil map or soil db is missing: only surface fluxes will be computed."); waterFluxesParameters.computeOnlySurface = true; } - // TODO considerare i vari casi - if (! landUseMap.isLoaded || landUnitList.empty()) + // check crop + if (processes.computeCrop) { - logInfo("WARNING: land use map or crop db is missing: default properties will be used (FALLOW)."); - landUnitList.clear(); - Crit3DLandUnit deafultLandUnit; - landUnitList.push_back(deafultLandUnit); + if (! landUseMap.isLoaded || landUnitList.empty()) + { + logInfo("WARNING: land use map or crop db is missing: crop computation will be deactivated."); + processes.computeCrop = false; + + // use default crop per surface properties + landUnitList.clear(); + Crit3DLandUnit deafultLandUnit; + landUnitList.push_back(deafultLandUnit); + } } // set computation depth @@ -212,12 +239,12 @@ bool Project3D::initializeWaterBalance3D() } logInfo("Computation depth: " + QString::number(computationSoilDepth) + " m"); - // Layers depth + // set layers depth setSoilLayers(); setLayersDepth(); logInfo("Nr of layers: " + QString::number(nrLayers)); - // nr of nodes + // set nr of nodes setIndexMaps(); logInfo("Nr of nodes: " + QString::number(nrNodes)); if (nrNodes == 0) @@ -228,11 +255,11 @@ bool Project3D::initializeWaterBalance3D() waterSinkSource.resize(nrNodes); - // Boundary + // set boundary if (!setLateralBoundary()) return false; logInfo("Lateral boundary computed"); - // Initiale soil fluxes + // initialize soil fluxes int myResult = soilFluxes3D::initialize(long(nrNodes), int(nrLayers), nrLateralLink, true, false, false); if (isCrit3dError(myResult, errorString)) { @@ -241,7 +268,7 @@ bool Project3D::initializeWaterBalance3D() } logInfo("Memory initialized"); - // Set properties for soil surface (roughness, pond) + // set properties for soil surface (roughness, pond) if (! setCrit3DSurfaces()) { logError(); @@ -264,144 +291,32 @@ bool Project3D::initializeWaterBalance3D() logInfo("Topology initialized"); if (! setCrit3DNodeSoil()) + { + logError(); return false; + } logInfo("Node properties initialized"); soilFluxes3D::setHydraulicProperties(MODIFIEDVANGENUCHTEN, MEAN_LOGARITHMIC, waterFluxesParameters.horizVertRatioConductivity); - double vmax = 10.0; // [m s-1] - double minimumDeltaT = DEM.header->cellSize / vmax; // [m] + double vmax = 10.0; // [m s-1] + double minimumDeltaT = DEM.header->cellSize / vmax; // [m] soilFluxes3D::setNumericalParameters(minimumDeltaT, 3600, 100, 10, 12, 3); // precision //soilFluxes3D::setNumericalParameters(minimumDeltaT, 3600, 100, 10, 12, 2); // speedy //soilFluxes3D::setNumericalParameters(minimumDeltaT, 3600, 100, 10, 12, 1); // very speedy (high error) - if (!initializeMatricPotential(waterFluxesParameters.initialWaterPotential)) // [m] - return false; - - logInfo("3D water balance initialized"); - return true; -} - - -bool Project3D::loadSoilDatabase(QString fileName) -{ - if (fileName == "") - { - logError("Missing Soil DB filename"); - return false; - } - - soilDbFileName = fileName; - fileName = getCompleteFileName(fileName, PATH_SOIL); - - if (! loadAllSoils(fileName, soilList, texturalClassList, geotechnicsClassList, fittingOptions, errorString)) + if (! initializeMatricPotential(waterFluxesParameters.initialWaterPotential)) // [m] { logError(); return false; } - nrSoils = unsigned(soilList.size()); - - logInfo("Soil database = " + fileName); - return true; -} - - -bool Project3D::loadCropDatabase(QString fileName) -{ - if (fileName == "") - { - logError("Missing Crop DB filename"); - return false; - } - - cropDbFileName = fileName; - fileName = getCompleteFileName(fileName, PATH_SOIL); - - QSqlDatabase dbCrop; - dbCrop = QSqlDatabase::addDatabase("QSQLITE", QUuid::createUuid().toString()); - dbCrop.setDatabaseName(fileName); - - if (!dbCrop.open()) - { - logError("Connection with crop database fail"); - return false; - } - - // land unit list - if (! loadLandUnitList(dbCrop, landUnitList, errorString)) - { - logError("Error in reading land_units table\n" + errorString); - return false; - } - - // crop list (same index of landUnitsList) - cropList.resize(landUnitList.size()); - for (int i = 0; i < landUnitList.size(); i++) - { - if (landUnitList[i].idCrop == "") continue; - - if (! loadCropParameters(dbCrop, landUnitList[i].idCrop, cropList[i], errorString)) - { - QString infoStr = "Error in reading crop data: " + landUnitList[i].idCrop; - logError(infoStr + "\n" + errorString); - return false; - } - } - logInfo("Crop/landUse database = " + fileName); + logInfo("3D water balance initialized"); return true; } -void Project3D::setSoilLayers() - { - double nextThickness; - double prevThickness = minThickness; - double depth = minThickness * 0.5; - - nrLayers = 1; - while (depth < computationSoilDepth) - { - nextThickness = MINVALUE(maxThickness, prevThickness * thickFactor); - depth = depth + (prevThickness + nextThickness) * 0.5; - prevThickness = nextThickness; - nrLayers++; - } -} - - -// set thickness and depth (center) of layers [m] -void Project3D::setLayersDepth() -{ - unsigned int lastLayer = nrLayers-1; - layerDepth.resize(nrLayers); - layerThickness.resize(nrLayers); - - layerDepth[0] = 0.0; - layerThickness[0] = 0.0; - - if (nrLayers <= 1) return; - - layerThickness[1] = minThickness; - layerDepth[1] = minThickness * 0.5; - - for (unsigned int i = 2; i < nrLayers; i++) - { - if (i == lastLayer) - { - layerThickness[i] = computationSoilDepth - (layerDepth[i-1] + layerThickness[i-1] / 2.0); - } - else - { - layerThickness[i] = MINVALUE(maxThickness, layerThickness[i-1] * thickFactor); - } - - layerDepth[i] = layerDepth[i-1] + (layerThickness[i-1] + layerThickness[i]) * 0.5; - } -} - - void Project3D::setIndexMaps() { indexMap.resize(nrLayers); @@ -471,6 +386,12 @@ bool Project3D::setLateralBoundary() bool Project3D::setCrit3DSurfaces() { + if (landUnitList.empty()) + { + errorString = "Error in setCrit3DSurfaces: missing land use data"; + return false; + } + for (int i = 0; i < landUnitList.size(); i++) { int result = soilFluxes3D::setSurfaceProperties(i, landUnitList[i].roughness, landUnitList[i].pond); @@ -699,14 +620,15 @@ bool Project3D::initializeMatricPotential(float psi) if (isCrit3dError(myResult, error)) { - logError("setCrit3DMatricPotential: " + error + " in row:" - + QString::number(row) + " col:" + QString::number(col)); + errorString = "setCrit3DMatricPotential: " + error + " in row:" + + QString::number(row) + " col:" + QString::number(col); return false; } } } } } + return true; } @@ -732,8 +654,11 @@ bool Project3D::setCrit3DNodeSoil() if (landUnitIndex != NODATA) soilFluxes3D::setNodeSurface(index, landUnitIndex); else - logError("Wrong surface definition in row, col: " - + QString::number(row) + "," + QString::number(col)); + { + errorString = "Wrong surface definition in row, col: " + + QString::number(row) + "," + QString::number(col); + return false; + } } else { @@ -744,10 +669,9 @@ bool Project3D::setCrit3DNodeSoil() if (horizonIndex == NODATA) { errorString = "function setCrit3DNodeSoil:\n No horizon definition in soil " - + QString::fromStdString(soilList[unsigned(soilIndex)].code) - + " depth: " + QString::number(layerDepth[layer]) - +"\nCheck soil totalDepth"; - logError(); + + QString::fromStdString(soilList[unsigned(soilIndex)].code) + + " depth: " + QString::number(layerDepth[layer]) + + "\nCheck soil totalDepth"; return false; } @@ -756,8 +680,8 @@ bool Project3D::setCrit3DNodeSoil() // check error if (isCrit3dError(myResult, errorString)) { - logError("setCrit3DNodeSoil:" + errorString + " in soil nr: " + QString::number(soilIndex) - + " horizon nr:" + QString::number(horizonIndex)); + errorString = "setCrit3DNodeSoil:" + errorString + " in soil nr: " + QString::number(soilIndex) + + " horizon nr:" + QString::number(horizonIndex); return false; } } @@ -830,419 +754,282 @@ bool Project3D::initializeSoilMoisture(int month) } -int Project3D::getSoilIndex(long row, long col) +/*! \brief computeWaterBalance3D + * \param totalTimeStep [s] + */ +void Project3D::computeWaterBalance3D(double totalTimeStep) { - if ( !soilIndexMap.isLoaded) - return NODATA; + double previousWaterContent = soilFluxes3D::getTotalWaterContent(); - int soilIndex = int(soilIndexMap.getValueFromRowCol(row, col)); + logInfo("total water [m^3]: " + QString::number(previousWaterContent)); + logInfo("precipitation [m^3]: " + QString::number(totalPrecipitation)); + logInfo("evaporation [m^3]: " + QString::number(-totalEvaporation)); + logInfo("transpiration [m^3]: " + QString::number(-totalTranspiration)); + logInfo("Compute water flow..."); - if (soilIndex == int(soilIndexMap.header->flag)) - return NODATA; + soilFluxes3D::initializeBalance(); - return soilIndex; -} + currentSeconds = 0; // [s] + double showTime = 60; // [s] + int currentStep = 0; + while (currentSeconds < totalTimeStep) + { + currentSeconds += soilFluxes3D::computeStep(totalTimeStep - currentSeconds); + + if (showEachTimeStep) + { + if (currentSeconds < totalTimeStep && int(currentSeconds / showTime) > currentStep) + { + currentStep = int(currentSeconds / showTime); + emit updateOutputSignal(); + } + } + } + double runoff = soilFluxes3D::getBoundaryWaterSumFlow(BOUNDARY_RUNOFF); + logInfo("runoff [m^3]: " + QString::number(runoff)); -bool Project3D::isWithinSoil(int soilIndex, double depth) -{ - if (soilIndex == int(NODATA) || soilIndex >= int(soilList.size())) return false; + double freeDrainage = soilFluxes3D::getBoundaryWaterSumFlow(BOUNDARY_FREEDRAINAGE); + logInfo("free drainage [m^3]: " + QString::number(freeDrainage)); - // check if depth is lower than lowerDepth of last horizon - unsigned int lastHorizon = soilList[unsigned(soilIndex)].nrHorizons -1; - double lowerDepth = soilList[unsigned(soilIndex)].horizon[lastHorizon].lowerDepth; + double lateralDrainage = soilFluxes3D::getBoundaryWaterSumFlow(BOUNDARY_FREELATERALDRAINAGE); + logInfo("lateral drainage [m^3]: " + QString::number(lateralDrainage)); - return (depth <= lowerDepth); + double currentWaterContent = soilFluxes3D::getTotalWaterContent(); + double forecastWaterContent = previousWaterContent + runoff + freeDrainage + lateralDrainage + + totalPrecipitation - totalEvaporation - totalTranspiration; + double massBalanceError = currentWaterContent - forecastWaterContent; + logInfo("Mass balance error [m^3]: " + QString::number(massBalanceError)); } -// upper depth of soil layer [m] +// ----------------------------------------- CROP and LAND USE ----------------------------------- -double Project3D::getSoilLayerTop(unsigned int i) +bool Project3D::loadCropDatabase(QString fileName) { - return layerDepth[i] - layerThickness[i] / 2.0; -} + if (fileName == "") + { + logError("Missing Crop DB filename"); + return false; + } -// lower depth of soil layer [m] -double Project3D::getSoilLayerBottom(unsigned int i) -{ - return layerDepth[i] + layerThickness[i] / 2.0; -} + cropDbFileName = fileName; + fileName = getCompleteFileName(fileName, PATH_SOIL); -// soil layer index from soil depth -int Project3D::getSoilLayerIndex(double depth) -{ - unsigned int i= 0; - while (depth > getSoilLayerBottom(i)) + QSqlDatabase dbCrop; + dbCrop = QSqlDatabase::addDatabase("QSQLITE", QUuid::createUuid().toString()); + dbCrop.setDatabaseName(fileName); + + if (!dbCrop.open()) { - if (i == nrLayers-1) + logError("Connection with crop database fail"); + return false; + } + + // land unit list + if (! loadLandUnitList(dbCrop, landUnitList, errorString)) + { + logError("Error in reading land_units table\n" + errorString); + return false; + } + + // crop list (same index of landUnitsList) + cropList.resize(landUnitList.size()); + for (int i = 0; i < landUnitList.size(); i++) + { + if (landUnitList[i].idCrop == "") continue; + + if (! loadCropParameters(dbCrop, landUnitList[i].idCrop, cropList[i], errorString)) { - logError("getSoilLayerIndex: wrong soil depth."); - return INDEX_ERROR; + QString infoStr = "Error in reading crop data: " + landUnitList[i].idCrop; + logError(infoStr + "\n" + errorString); + return false; } - i++; } - return signed(i); + logInfo("Crop/landUse database = " + fileName); + return true; } -bool Project3D::saveHourlyMeteoOutput(meteoVariable myVar, const QString& myPath, QDateTime myTime) +int Project3D::getLandUnitIdUTM(double x, double y) { - gis::Crit3DRasterGrid* myRaster = getHourlyMeteoRaster(myVar); - if (myRaster == nullptr) return false; + if (! landUseMap.isLoaded) + return NODATA; - QString fileName = getOutputNameHourly(myVar, myTime); - QString outputFileName = myPath + fileName; + int id = int(gis::getValueFromXY(landUseMap, x, y)); - std::string errStr; - if (! gis::writeEsriGrid(outputFileName.toStdString(), myRaster, errStr)) + if (id == int(landUseMap.header->flag)) { - logError(QString::fromStdString(errStr)); - return false; + return NODATA; } else - return true; + { + return id; + } } -bool Project3D::aggregateAndSaveDailyMap(meteoVariable myVar, aggregationMethod myAggregation, const Crit3DDate& myDate, - const QString& dailyPath, const QString& hourlyPath) +int Project3D::getLandUnitIdGeo(double lat, double lon) { - std::string myError; - int myTimeStep = int(3600. / meteoSettings->getHourlyIntervals()); - Crit3DTime myTimeIni(myDate, myTimeStep); - Crit3DTime myTimeFin(myDate.addDays(1), 0.); - - gis::Crit3DRasterGrid* myMap = new gis::Crit3DRasterGrid(); - myMap->initializeGrid(DEM); - gis::Crit3DRasterGrid* myAggrMap = new gis::Crit3DRasterGrid(); - myAggrMap->initializeGrid(DEM); + double x, y; + gis::latLonToUtmForceZone(gisSettings.utmZone, lat, lon, &x, &y); - long myRow, myCol; - int nrAggrMap = 0; + return getLandUnitIdUTM(x, y); +} - for (Crit3DTime myTime = myTimeIni; myTime<=myTimeFin; myTime=myTime.addSeconds(myTimeStep)) - { - QString hourlyFileName = getOutputNameHourly(myVar, getQDateTime(myTime)); - if (gis::readEsriGrid((hourlyPath + hourlyFileName).toStdString(), myMap, myError)) - { - if (myTime == myTimeIni) - { - for (myRow = 0; myRow < myAggrMap->header->nrRows; myRow++) - for (myCol = 0; myCol < myAggrMap->header->nrCols; myCol++) - myAggrMap->value[myRow][myCol] = myMap->value[myRow][myCol]; - nrAggrMap++; - } - else - { - if (myAggregation == aggrMin) - gis::mapAlgebra(myAggrMap, myMap, myAggrMap, operationMin); - else if (myAggregation == aggrMax) - gis::mapAlgebra(myAggrMap, myMap, myAggrMap, operationMax); - else if (myAggregation == aggrSum || myAggregation == aggrAverage || myAggregation == aggrIntegral) - gis::mapAlgebra(myAggrMap, myMap, myAggrMap, operationSum); - else - { - logError("wrong aggregation type in function 'aggregateAndSaveDailyMap'"); - return(false); - } - nrAggrMap++; - } - } - } - - if (myAggregation == aggrAverage) - gis::mapAlgebra(myAggrMap, nrAggrMap, myAggrMap, operationDivide); - else if (myAggregation == aggrSum) +int Project3D::getLandUnitIndexRowCol(int row, int col) +{ + if (! landUseMap.isLoaded || landUnitList.empty()) { - if (myVar == globalIrradiance || myVar == directIrradiance || myVar == diffuseIrradiance || myVar == reflectedIrradiance) - gis::mapAlgebra(myAggrMap, float(myTimeStep / 1000000.0), myAggrMap, operationProduct); + return 0; // default } - meteoVariable dailyVar = getDailyMeteoVarFromHourly(myVar, myAggregation); - QString varName = QString::fromStdString(MapDailyMeteoVarToString.at(dailyVar)); - - QString filename = getOutputNameDaily(varName , "", getQDate(myDate)); - - QString outputFileName = dailyPath + filename; - bool isOk = gis::writeEsriGrid(outputFileName.toStdString(), myAggrMap, myError); - - myMap->clear(); - myAggrMap->clear(); + double x, y; + DEM.getXY(row, col, x, y); - if (! isOk) + int id = getLandUnitIdUTM(x, y); + if (id == NODATA) { - logError("aggregateMapToDaily: " + QString::fromStdString(myError)); - return false; + return NODATA; } - return true; + return getLandUnitIndex(landUnitList, id); } -// compute evaporation [mm] -// TODO check pond -double Project3D::assignEvaporation(int row, int col, double lai) -{ - double depthCoeff, thickCoeff, layerCoeff; - double availableWater; // [mm] +// ------------------------------------ SOIL -------------------------------------- - double const MAX_EVAPORATION_DEPTH = 0.2; // [m] - int lastEvapLayer; // [-] index - if (computationSoilDepth >= MAX_EVAPORATION_DEPTH) - { - lastEvapLayer = getSoilLayerIndex(MAX_EVAPORATION_DEPTH); - } - else +bool Project3D::loadSoilDatabase(QString fileName) +{ + if (fileName == "") { - lastEvapLayer = getSoilLayerIndex(computationSoilDepth); + logError("Missing Soil DB filename"); + return false; } - double area = DEM.header->cellSize * DEM.header->cellSize; // [m2] - - double et0 = double(hourlyMeteoMaps->mapHourlyET0->value[row][col]); // [mm] - double potentialEvaporation = getMaxEvaporation(et0, lai); // [mm] - double realEvaporation = 0; // [mm] + soilDbFileName = fileName; + fileName = getCompleteFileName(fileName, PATH_SOIL); - for (unsigned int layer=0; layer <= unsigned(lastEvapLayer); layer++) + if (! loadAllSoils(fileName, soilList, texturalClassList, geotechnicsClassList, fittingOptions, errorString)) { - long nodeIndex = long(indexMap.at(layer).value[row][col]); - - // layer coefficient - if (layer == 0) - { - // surface: water level [m] -> [mm] - availableWater = getCriteria3DVar(availableWaterContent, nodeIndex) * 1000; - layerCoeff = 1; - } - else - { - // sub-surface: volumetric awc [m^3 m^-3] - availableWater = getCriteria3DVar(availableWaterContent, nodeIndex); - // [m] -> [mm] - availableWater *= layerThickness[layer] * 1000; - - depthCoeff = layerDepth[layer] / MAX_EVAPORATION_DEPTH; - thickCoeff = layerThickness[layer] / 0.04; - layerCoeff = exp(-EULER * depthCoeff) * thickCoeff; - } - - double residualEvap = potentialEvaporation - realEvaporation; // [mm] - double layerEvap = MINVALUE(potentialEvaporation * layerCoeff, residualEvap); // [mm] - layerEvap = MINVALUE(layerEvap, availableWater); - - if (layerEvap > 0) - { - realEvaporation += layerEvap; - double flow = area * (layerEvap / 1000); // [m3 h-1] - waterSinkSource.at(unsigned(nodeIndex)) -= (flow / 3600); // [m3 s-1] - } + logError(); + return false; } + nrSoils = unsigned(soilList.size()); - return realEvaporation; + logInfo("Soil database = " + fileName); + return true; } -// assign real crop transpiration -// return sum of crop transpiration over the soil column -double Project3D::assignTranspiration(int row, int col, double lai) -{ - double soilColumnTranspirationSum = 0; - - // check - /*if (idCrop == "" || ! isLiving) return 0; - if (roots.rootDepth <= roots.rootDepthMin) return 0; - if (roots.firstRootLayer == NODATA) return 0; - if (maxTranspiration < EPSILON) return 0; - - double thetaWP; // [m3 m-3] volumetric water content at Wilting Point - double cropWP; // [mm] wilting point specific for crop - double waterSurplusThreshold; // [mm] water surplus stress threshold - double waterScarcityThreshold; // [mm] water scarcity stress threshold - double WSS; // [] water surplus stress - - double TRs=0.0; // [mm] actual transpiration with only water scarsity stress - double TRe=0.0; // [mm] actual transpiration with only water surplus stress - double totRootDensityWithoutStress = 0.0; // [-] - double redistribution = 0.0; // [mm] +void Project3D::setSoilLayers() + { + double nextThickness; + double prevThickness = minThickness; + double depth = minThickness * 0.5; - // initialize - unsigned int nrLayers = unsigned(soilLayers.size()); - bool* isLayerStressed = new bool[nrLayers]; - for (unsigned int i = 0; i < nrLayers; i++) + nrLayers = 1; + while (depth < computationSoilDepth) { - isLayerStressed[i] = false; - layerTranspiration[i] = 0; + nextThickness = MINVALUE(maxThickness, prevThickness * thickFactor); + depth = depth + (prevThickness + nextThickness) * 0.5; + prevThickness = nextThickness; + nrLayers++; } +} - // water surplus - if (isWaterSurplusResistant()) - WSS = 0.0; - else - WSS = 0.5; - for (unsigned int i = unsigned(roots.firstRootLayer); i <= unsigned(roots.lastRootLayer); i++) - { - // [mm] - waterSurplusThreshold = soilLayers[i].SAT - (WSS * (soilLayers[i].SAT - soilLayers[i].FC)); +// set thickness and depth (center) of layers [m] +void Project3D::setLayersDepth() +{ + unsigned int lastLayer = nrLayers-1; + layerDepth.resize(nrLayers); + layerThickness.resize(nrLayers); - thetaWP = soil::thetaFromSignPsi(-soil::cmTokPa(psiLeaf), *(soilLayers[i].horizonPtr)); - // [mm] - cropWP = thetaWP * soilLayers[i].thickness * soilLayers[i].soilFraction * 1000.0; + layerDepth[0] = 0.0; + layerThickness[0] = 0.0; - // [mm] - waterScarcityThreshold = soilLayers[i].FC - fRAW * (soilLayers[i].FC - cropWP); + if (nrLayers <= 1) return; - if ((soilLayers[i].waterContent - waterSurplusThreshold) > EPSILON) - { - // WATER SURPLUS - layerTranspiration[i] = maxTranspiration * roots.rootDensity[i] * - ((soilLayers[i].SAT - soilLayers[i].waterContent) - / (soilLayers[i].SAT - waterSurplusThreshold)); + layerThickness[1] = minThickness; + layerDepth[1] = minThickness * 0.5; - TRe += layerTranspiration[i]; - TRs += maxTranspiration * roots.rootDensity[i]; - isLayerStressed[i] = true; - } - else if (soilLayers[i].waterContent < waterScarcityThreshold) + for (unsigned int i = 2; i < nrLayers; i++) + { + if (i == lastLayer) { - // WATER SCARSITY - if (soilLayers[i].waterContent <= cropWP) - { - layerTranspiration[i] = 0; - } - else - { - layerTranspiration[i] = maxTranspiration * roots.rootDensity[i] * - ((soilLayers[i].waterContent - cropWP) / (waterScarcityThreshold - cropWP)); - } - - TRs += layerTranspiration[i]; - TRe += maxTranspiration * roots.rootDensity[i]; - isLayerStressed[i] = true; + layerThickness[i] = computationSoilDepth - (layerDepth[i-1] + layerThickness[i-1] / 2.0); } else { - // normal conditions - layerTranspiration[i] = maxTranspiration * roots.rootDensity[i]; - - TRs += layerTranspiration[i]; - TRe += layerTranspiration[i]; - - if ((soilLayers[i].waterContent - layerTranspiration[i]) > waterScarcityThreshold) - { - isLayerStressed[i] = false; - totRootDensityWithoutStress += roots.rootDensity[i]; - } - else - { - isLayerStressed[i] = true; - } - } - } - - // WATER STRESS [-] - double firstWaterStress = 1 - (TRs / maxTranspiration); - - // Hydraulic redistribution - // the movement of water from moist to dry soil through plant roots - // TODO add numerical process - if (firstWaterStress > EPSILON && totRootDensityWithoutStress > EPSILON) - { - // redistribution acts on not stressed roots - redistribution = MINVALUE(firstWaterStress, totRootDensityWithoutStress) * maxTranspiration; - - for (int i = roots.firstRootLayer; i <= roots.lastRootLayer; i++) - { - if (! isLayerStressed[i]) - { - double addLayerTransp = redistribution * (roots.rootDensity[unsigned(i)] / totRootDensityWithoutStress); - layerTranspiration[unsigned(i)] += addLayerTransp; - TRs += addLayerTransp; - } + layerThickness[i] = MINVALUE(maxThickness, layerThickness[i-1] * thickFactor); } - } - - waterStress = 1 - (TRs / maxTranspiration); - for (int i = roots.firstRootLayer; i <= roots.lastRootLayer; i++) - { - totalTranspiration += layerTranspiration[unsigned(i)]; + layerDepth[i] = layerDepth[i-1] + (layerThickness[i-1] + layerThickness[i]) * 0.5; } - */ - - return soilColumnTranspirationSum; } -// input: timeStep [s] -void Project3D::computeWaterBalance3D(double totalTimeStep) +int Project3D::getSoilIndex(long row, long col) { - double previousWaterContent = soilFluxes3D::getTotalWaterContent(); + if ( !soilIndexMap.isLoaded) + return NODATA; - this->logInfo("total water [m^3]: " + QString::number(previousWaterContent)); - this->logInfo("precipitation [m^3]: " + QString::number(totalPrecipitation)); - this->logInfo("evaporation [m^3]: " + QString::number(-totalEvaporation)); - this->logInfo("transpiration [m^3]: " + QString::number(-totalTranspiration)); - this->logInfo("Compute water flow..."); + int soilIndex = int(soilIndexMap.getValueFromRowCol(row, col)); - soilFluxes3D::initializeBalance(); + if (soilIndex == int(soilIndexMap.header->flag)) + return NODATA; - currentSeconds = 0; // [s] - double showTime = 60; // [s] - int currentStep = 0; - while (currentSeconds < totalTimeStep) - { - currentSeconds += soilFluxes3D::computeStep(totalTimeStep - currentSeconds); + return soilIndex; +} - if (showEachTimeStep) - { - if (currentSeconds < totalTimeStep && int(currentSeconds / showTime) > currentStep) - { - currentStep = int(currentSeconds / showTime); - emit updateOutputSignal(); - } - } - } - double runoff = soilFluxes3D::getBoundaryWaterSumFlow(BOUNDARY_RUNOFF); - this->logInfo("runoff [m^3]: " + QString::number(runoff)); +bool Project3D::isWithinSoil(int soilIndex, double depth) +{ + if (soilIndex == int(NODATA) || soilIndex >= int(soilList.size())) return false; - double freeDrainage = soilFluxes3D::getBoundaryWaterSumFlow(BOUNDARY_FREEDRAINAGE); - this->logInfo("free drainage [m^3]: " + QString::number(freeDrainage)); + // check if depth is lower than lowerDepth of last horizon + unsigned int lastHorizon = soilList[unsigned(soilIndex)].nrHorizons -1; + double lowerDepth = soilList[unsigned(soilIndex)].horizon[lastHorizon].lowerDepth; - double lateralDrainage = soilFluxes3D::getBoundaryWaterSumFlow(BOUNDARY_FREELATERALDRAINAGE); - this->logInfo("lateral drainage [m^3]: " + QString::number(lateralDrainage)); + return (depth <= lowerDepth); +} - double currentWaterContent = soilFluxes3D::getTotalWaterContent(); - double forecastWaterContent = previousWaterContent + runoff + freeDrainage + lateralDrainage - + totalPrecipitation - totalEvaporation - totalTranspiration; - double massBalanceError = currentWaterContent - forecastWaterContent; - this->logInfo("Mass balance error [m^3]: " + QString::number(massBalanceError)); + +// upper depth of soil layer [m] +double Project3D::getSoilLayerTop(unsigned int i) +{ + return layerDepth[i] - layerThickness[i] / 2.0; } +// lower depth of soil layer [m] +double Project3D::getSoilLayerBottom(unsigned int i) +{ + return layerDepth[i] + layerThickness[i] / 2.0; +} -bool Project3D::updateCrop(QDateTime myTime) +// soil layer index from soil depth +int Project3D::getSoilLayerIndex(double depth) { - for (long row = 0; row < DEM.header->nrRows ; row++) + unsigned int i= 0; + while (depth > getSoilLayerBottom(i)) { - for (long col = 0; col < DEM.header->nrCols; col++) + if (i == nrLayers-1) { - if (int(DEM.value[row][col]) != int(DEM.header->flag)) - { - // TODO read crop - // compute LAI and kc - // state variables - } + logError("getSoilLayerIndex: wrong soil depth."); + return INDEX_ERROR; } + i++; } - return true; + return signed(i); } +// ------------------------------------ INPUT MAP ------------------------------------------ + bool Project3D::interpolateHourlyMeteoVar(meteoVariable myVar, const QDateTime& myTime) { if (myVar == airRelHumidity && interpolationSettings.getUseDewPoint()) @@ -1288,38 +1075,151 @@ bool Project3D::interpolateAndSaveHourlyMeteo(meteoVariable myVar, const QDateTi } -void Project3D::assignPrecipitation() +// ----------------------------------------- OUTPUT MAP ------------------------------------------ + +bool Project3D::setCriteria3DMap(criteria3DVariable var, int layerIndex) { - // initialize - totalPrecipitation = 0; + if (layerIndex >= indexMap.size()) + { + errorString = "Layer is not defined: " + QString::number(layerIndex); + return false; + } - double area = DEM.header->cellSize * DEM.header->cellSize; + criteria3DMap.initializeGrid(indexMap.at(layerIndex)); + + for (int row = 0; row < indexMap.at(layerIndex).header->nrRows; row++) + { + for (int col = 0; col < indexMap.at(layerIndex).header->nrCols; col++) + { + long nodeIndex = indexMap.at(layerIndex).value[row][col]; + if (nodeIndex != indexMap.at(layerIndex).header->flag) + { + double value = getCriteria3DVar(var, nodeIndex); + + if (value == NODATA) + { + criteria3DMap.value[row][col] = criteria3DMap.header->flag; + } + else + { + if (var == waterContent && layerIndex == 0) + { + // surface + value *= 1000; // [m] -> [mm] + } + criteria3DMap.value[row][col] = value; + } + } + else + { + criteria3DMap.value[row][col] = criteria3DMap.header->flag; + } + } + } + + gis::updateMinMaxRasterGrid(&criteria3DMap); + return true; +} + + +bool Project3D::saveHourlyMeteoOutput(meteoVariable myVar, const QString& myPath, QDateTime myTime) +{ + gis::Crit3DRasterGrid* myRaster = getHourlyMeteoRaster(myVar); + if (myRaster == nullptr) return false; + + QString fileName = getOutputNameHourly(myVar, myTime); + QString outputFileName = myPath + fileName; + + std::string errStr; + if (! gis::writeEsriGrid(outputFileName.toStdString(), myRaster, errStr)) + { + logError(QString::fromStdString(errStr)); + return false; + } + else + return true; +} + + +bool Project3D::aggregateAndSaveDailyMap(meteoVariable myVar, aggregationMethod myAggregation, const Crit3DDate& myDate, + const QString& dailyPath, const QString& hourlyPath) +{ + std::string myError; + int myTimeStep = int(3600. / meteoSettings->getHourlyIntervals()); + Crit3DTime myTimeIni(myDate, myTimeStep); + Crit3DTime myTimeFin(myDate.addDays(1), 0.); + + gis::Crit3DRasterGrid* myMap = new gis::Crit3DRasterGrid(); + myMap->initializeGrid(DEM); + gis::Crit3DRasterGrid* myAggrMap = new gis::Crit3DRasterGrid(); + myAggrMap->initializeGrid(DEM); - // precipitation - for (long row = 0; row < indexMap.at(0).header->nrRows; row++) + long myRow, myCol; + int nrAggrMap = 0; + + for (Crit3DTime myTime = myTimeIni; myTime<=myTimeFin; myTime=myTime.addSeconds(myTimeStep)) { - for (long col = 0; col < indexMap.at(0).header->nrCols; col++) + QString hourlyFileName = getOutputNameHourly(myVar, getQDateTime(myTime)); + if (gis::readEsriGrid((hourlyPath + hourlyFileName).toStdString(), myMap, myError)) { - int surfaceIndex = indexMap.at(0).value[row][col]; - if (surfaceIndex != indexMap.at(0).header->flag) + if (myTime == myTimeIni) { - double waterSource = 0; - float prec = hourlyMeteoMaps->mapHourlyPrec->value[row][col]; - if (! isEqual(prec, hourlyMeteoMaps->mapHourlyPrec->header->flag)) - waterSource += prec; + for (myRow = 0; myRow < myAggrMap->header->nrRows; myRow++) + for (myCol = 0; myCol < myAggrMap->header->nrCols; myCol++) + myAggrMap->value[myRow][myCol] = myMap->value[myRow][myCol]; - if (waterSource > 0) + nrAggrMap++; + } + else + { + if (myAggregation == aggrMin) + gis::mapAlgebra(myAggrMap, myMap, myAggrMap, operationMin); + else if (myAggregation == aggrMax) + gis::mapAlgebra(myAggrMap, myMap, myAggrMap, operationMax); + else if (myAggregation == aggrSum || myAggregation == aggrAverage || myAggregation == aggrIntegral) + gis::mapAlgebra(myAggrMap, myMap, myAggrMap, operationSum); + else { - double flow = area * (waterSource / 1000.); // [m3 h-1] - totalPrecipitation += flow; - waterSinkSource[unsigned(surfaceIndex)] += flow / 3600.; // [m3 s-1] + logError("wrong aggregation type in function 'aggregateAndSaveDailyMap'"); + return(false); } + nrAggrMap++; } } } + + if (myAggregation == aggrAverage) + gis::mapAlgebra(myAggrMap, nrAggrMap, myAggrMap, operationDivide); + else if (myAggregation == aggrSum) + { + if (myVar == globalIrradiance || myVar == directIrradiance || myVar == diffuseIrradiance || myVar == reflectedIrradiance) + gis::mapAlgebra(myAggrMap, float(myTimeStep / 1000000.0), myAggrMap, operationProduct); + } + + meteoVariable dailyVar = getDailyMeteoVarFromHourly(myVar, myAggregation); + QString varName = QString::fromStdString(MapDailyMeteoVarToString.at(dailyVar)); + + QString filename = getOutputNameDaily(varName , "", getQDate(myDate)); + + QString outputFileName = dailyPath + filename; + bool isOk = gis::writeEsriGrid(outputFileName.toStdString(), myAggrMap, myError); + + myMap->clear(); + myAggrMap->clear(); + + if (! isOk) + { + logError("aggregateMapToDaily: " + QString::fromStdString(myError)); + return false; + } + + return true; } + +// ---------------------------------- SINK / SOURCE ----------------------------------- + bool Project3D::setSinkSource() { for (unsigned long i = 0; i < nrNodes; i++) @@ -1335,157 +1235,246 @@ bool Project3D::setSinkSource() return true; } +/*! \brief getCoveredSurfaceFraction + * \param leafAreaIndex [m2 m-2] + * \ref Liangxia Zhang, Zhongmin Hu, Jiangwen Fan, Decheng Zhou & Fengpei Tang, 2014 + * A meta-analysis of the canopy light extinction coefficient in terrestrial ecosystems + * "Cropland had the highest value of K (0.62), followed by broadleaf forest (0.59) + * shrubland (0.56), grassland (0.50), and needleleaf forest (0.45)" + * \return covered surface fraction [-] + */ +double getCoveredSurfaceFraction(double leafAreaIndex) +{ + if (leafAreaIndex < EPSILON) return 0; + + double k = 0.6; // [-] light extinction coefficient + return 1 - exp(-k * leafAreaIndex); +} + -void Project3D::computeBareSoilEvaporation() +/*! \brief getMaxEvaporation + * \param ET0: potential evapo-transpiration [mm] + * \param leafAreaIndex [m2 m-2] + * \return maximum soil evaporation [mm] + */ +double getMaxSoilEvaporation(double ET0, double leafAreaIndex) { - // initialize - totalEvaporation = 0; + double evapMax = ET0 * (1.0 - getCoveredSurfaceFraction(leafAreaIndex)); + // TODO check evaporation on free water + return evapMax * 0.66; +} - double area = DEM.header->cellSize * DEM.header->cellSize; - for (int row = 0; row < indexMap.at(0).header->nrRows; row++) +/*! \brief getMaxCropTranspiration + * \param ET0: potential evapo-transpiration [mm] + * \param leafAreaIndex [m2 m-2] + * \param maxKc: maximum crop coefficient [-] + * \return maximum crop transpiration [mm] + */ +double getMaxCropTranspiration(double ET0, double leafAreaIndex, double maxKc) +{ + double covSurfFraction = getCoveredSurfaceFraction(leafAreaIndex); + double kcFactor = 1 + (maxKc - 1) * covSurfFraction; + return ET0 * covSurfFraction * kcFactor; +} + + +// compute evaporation [mm] +// TODO check pond +double Project3D::assignEvaporation(int row, int col, double lai) +{ + double depthCoeff, thickCoeff, layerCoeff; + double availableWater; // [mm] + + double const MAX_EVAPORATION_DEPTH = 0.2; // [m] + int lastEvapLayer; // [-] index + if (computationSoilDepth >= MAX_EVAPORATION_DEPTH) + { + lastEvapLayer = getSoilLayerIndex(MAX_EVAPORATION_DEPTH); + } + else + { + lastEvapLayer = getSoilLayerIndex(computationSoilDepth); + } + + double area = DEM.header->cellSize * DEM.header->cellSize; // [m2] + + double et0 = double(hourlyMeteoMaps->mapHourlyET0->value[row][col]); // [mm] + double potentialEvaporation = getMaxSoilEvaporation(et0, lai); // [mm] + double realEvaporation = 0; // [mm] + + for (unsigned int layer=0; layer <= unsigned(lastEvapLayer); layer++) { - for (int col = 0; col < indexMap.at(0).header->nrCols; col++) + long nodeIndex = long(indexMap.at(layer).value[row][col]); + + // layer coefficient + if (layer == 0) { - int surfaceIndex = indexMap.at(0).value[row][col]; - if (surfaceIndex != indexMap.at(0).header->flag) - { - double lai = 0; + // surface: water level [m] -> [mm] + availableWater = getCriteria3DVar(availableWaterContent, nodeIndex) * 1000; + layerCoeff = 1; + } + else + { + // sub-surface: volumetric awc [m^3 m^-3] + availableWater = getCriteria3DVar(availableWaterContent, nodeIndex); + // [m] -> [mm] + availableWater *= layerThickness[layer] * 1000; - double realEvap = assignEvaporation(row, col, lai); // [mm] - double flow = area * (realEvap / 1000.); // [m3 h-1] - totalEvaporation += flow; - } + depthCoeff = layerDepth[layer] / MAX_EVAPORATION_DEPTH; + thickCoeff = layerThickness[layer] / 0.04; + layerCoeff = exp(-EULER * depthCoeff) * thickCoeff; } - } -} + double residualEvap = potentialEvaporation - realEvaporation; // [mm] + double layerEvap = MINVALUE(potentialEvaporation * layerCoeff, residualEvap); // [mm] + layerEvap = MINVALUE(layerEvap, availableWater); -// crop transpiration -// TODO -/* - totalTranspiration = 0; - for (unsigned int layerIndex=1; layerIndex < nrLayers; layerIndex++) - { - for (long row = 0; row < indexMap.at(size_t(layerIndex)).header->nrRows; row++) + if (layerEvap > 0) { - for (long col = 0; col < indexMap.at(size_t(layerIndex)).header->nrCols; col++) - { - int nodeIndex = indexMap.at(size_t(layerIndex)).value[row][col]; - if (nodeIndex != indexMap.at(size_t(layerIndex)).header->flag) - { - float transp = outputPlantMaps->transpirationLayerMaps[layerIndex]->value[row][col]; - if (int(transp) != int(outputPlantMaps->transpirationLayerMaps[layerIndex]->header->flag)) - { - flow = area * (transp / 1000.0); //[m^3/h] - totalTranspiration += flow; - waterSinkSource.at(unsigned(nodeIndex)) -= flow / 3600.0; //[m^3/s] - } - } - } + realEvaporation += layerEvap; + double flow = area * (layerEvap / 1000); // [m3 h-1] + waterSinkSource.at(unsigned(nodeIndex)) -= (flow / 3600); // [m3 s-1] } } - */ + return realEvaporation; +} +// assign real crop transpiration +// return sum of crop transpiration over the soil column +double Project3D::assignTranspiration(int row, int col, double lai) +{ + double transpirationSum = 0; + // check + /*if (idCrop == "" || ! isLiving) return 0; + if (roots.rootDepth <= roots.rootDepthMin) return 0; + if (roots.firstRootLayer == NODATA) return 0; + if (maxTranspiration < EPSILON) return 0; + double thetaWP; // [m3 m-3] volumetric water content at Wilting Point + double cropWP; // [mm] wilting point specific for crop + double waterSurplusThreshold; // [mm] water surplus stress threshold + double waterScarcityThreshold; // [mm] water scarcity stress threshold + double WSS; // [] water surplus stress + double TRs=0.0; // [mm] actual transpiration with only water scarsity stress + double TRe=0.0; // [mm] actual transpiration with only water surplus stress + double totRootDensityWithoutStress = 0.0; // [-] + double redistribution = 0.0; // [mm] -bool Project3D::setCriteria3DMap(criteria3DVariable var, int layerIndex) -{ - if (layerIndex >= indexMap.size()) + // initialize + unsigned int nrLayers = unsigned(soilLayers.size()); + bool* isLayerStressed = new bool[nrLayers]; + for (unsigned int i = 0; i < nrLayers; i++) { - errorString = "Layer is not defined: " + QString::number(layerIndex); - return false; + isLayerStressed[i] = false; + layerTranspiration[i] = 0; } - criteria3DMap.initializeGrid(indexMap.at(layerIndex)); + // water surplus + if (isWaterSurplusResistant()) + WSS = 0.0; + else + WSS = 0.5; - for (int row = 0; row < indexMap.at(layerIndex).header->nrRows; row++) + for (unsigned int i = unsigned(roots.firstRootLayer); i <= unsigned(roots.lastRootLayer); i++) { - for (int col = 0; col < indexMap.at(layerIndex).header->nrCols; col++) + // [mm] + waterSurplusThreshold = soilLayers[i].SAT - (WSS * (soilLayers[i].SAT - soilLayers[i].FC)); + + thetaWP = soil::thetaFromSignPsi(-soil::cmTokPa(psiLeaf), *(soilLayers[i].horizonPtr)); + // [mm] + cropWP = thetaWP * soilLayers[i].thickness * soilLayers[i].soilFraction * 1000.0; + + // [mm] + waterScarcityThreshold = soilLayers[i].FC - fRAW * (soilLayers[i].FC - cropWP); + + if ((soilLayers[i].waterContent - waterSurplusThreshold) > EPSILON) { - long nodeIndex = indexMap.at(layerIndex).value[row][col]; - if (nodeIndex != indexMap.at(layerIndex).header->flag) - { - double value = getCriteria3DVar(var, nodeIndex); + // WATER SURPLUS + layerTranspiration[i] = maxTranspiration * roots.rootDensity[i] * + ((soilLayers[i].SAT - soilLayers[i].waterContent) + / (soilLayers[i].SAT - waterSurplusThreshold)); - if (value == NODATA) - { - criteria3DMap.value[row][col] = criteria3DMap.header->flag; - } - else - { - if (var == waterContent && layerIndex == 0) - { - // surface - value *= 1000; // [m] -> [mm] - } - criteria3DMap.value[row][col] = value; - } + TRe += layerTranspiration[i]; + TRs += maxTranspiration * roots.rootDensity[i]; + isLayerStressed[i] = true; + } + else if (soilLayers[i].waterContent < waterScarcityThreshold) + { + // WATER SCARSITY + if (soilLayers[i].waterContent <= cropWP) + { + layerTranspiration[i] = 0; } else { - criteria3DMap.value[row][col] = criteria3DMap.header->flag; + layerTranspiration[i] = maxTranspiration * roots.rootDensity[i] * + ((soilLayers[i].waterContent - cropWP) / (waterScarcityThreshold - cropWP)); } - } - } - - gis::updateMinMaxRasterGrid(&criteria3DMap); - return true; -} - -int Project3D::getLandUnitIdUTM(double x, double y) -{ - if (! landUseMap.isLoaded) - return NODATA; + TRs += layerTranspiration[i]; + TRe += maxTranspiration * roots.rootDensity[i]; + isLayerStressed[i] = true; + } + else + { + // normal conditions + layerTranspiration[i] = maxTranspiration * roots.rootDensity[i]; - int id = int(gis::getValueFromXY(landUseMap, x, y)); + TRs += layerTranspiration[i]; + TRe += layerTranspiration[i]; - if (id == int(landUseMap.header->flag)) - { - return NODATA; - } - else - { - return id; + if ((soilLayers[i].waterContent - layerTranspiration[i]) > waterScarcityThreshold) + { + isLayerStressed[i] = false; + totRootDensityWithoutStress += roots.rootDensity[i]; + } + else + { + isLayerStressed[i] = true; + } + } } -} - - -int Project3D::getLandUnitIdGeo(double lat, double lon) -{ - double x, y; - gis::latLonToUtmForceZone(gisSettings.utmZone, lat, lon, &x, &y); - - return getLandUnitIdUTM(x, y); -} + // WATER STRESS [-] + double firstWaterStress = 1 - (TRs / maxTranspiration); -int Project3D::getLandUnitIndexRowCol(int row, int col) -{ - if (! landUseMap.isLoaded || landUnitList.empty()) + // Hydraulic redistribution + // the movement of water from moist to dry soil through plant roots + // TODO add numerical process + if (firstWaterStress > EPSILON && totRootDensityWithoutStress > EPSILON) { - return 0; // default + // redistribution acts on not stressed roots + redistribution = MINVALUE(firstWaterStress, totRootDensityWithoutStress) * maxTranspiration; + + for (int i = roots.firstRootLayer; i <= roots.lastRootLayer; i++) + { + if (! isLayerStressed[i]) + { + double addLayerTransp = redistribution * (roots.rootDensity[unsigned(i)] / totRootDensityWithoutStress); + layerTranspiration[unsigned(i)] += addLayerTransp; + TRs += addLayerTransp; + } + } } - double x, y; - DEM.getXY(row, col, x, y); + waterStress = 1 - (TRs / maxTranspiration); - int id = getLandUnitIdUTM(x, y); - if (id == NODATA) + for (int i = roots.firstRootLayer; i <= roots.lastRootLayer; i++) { - return NODATA; + totalTranspiration += layerTranspiration[unsigned(i)]; } + */ - return getLandUnitIndex(landUnitList, id); + return transpirationSum; } -// ------------------------- other functions ------------------------- +// ------------------------------ other functions ---------------------------------- bool isCrit3dError(int result, QString& error) { @@ -1652,16 +1641,3 @@ QString getDailyPrefixFromVar(QDate myDate, criteria3DVariable myVar) return fileName; } - -double getMaxEvaporation(double ET0, double LAI) -{ - const double ke = 0.6; //[-] light extinction factor - const double maxEvaporationRatio = 0.66; - - double Kc = exp(-ke * LAI); - return(ET0 * Kc * maxEvaporationRatio); -} - - - - diff --git a/bin/CRITERIA3D/shared/project3D.h b/bin/CRITERIA3D/shared/project3D.h index abc594589..ba07f257a 100644 --- a/bin/CRITERIA3D/shared/project3D.h +++ b/bin/CRITERIA3D/shared/project3D.h @@ -40,6 +40,19 @@ }; + class Crit3DProcesses + { + public: + + bool computeMeteo, computeRadiation, computeWater; + bool computeEvaporation, computeCrop, computeSnow, computeSolutes; + bool computeHeat, computeAdvectiveHeat, computeLatentHeat; + + Crit3DProcesses(); + void initialize(); + }; + + class Project3D : public Project { Q_OBJECT @@ -62,6 +75,7 @@ bool isCriteria3DInitialized; bool showEachTimeStep; + Crit3DProcesses processes; WaterFluxesParameters waterFluxesParameters; QString soilDbFileName; @@ -144,8 +158,6 @@ bool interpolateHourlyMeteoVar(meteoVariable myVar, const QDateTime& myTime); double assignEvaporation(int row, int col, double lai); double assignTranspiration(int row, int col, double lai); - void computeBareSoilEvaporation(); - void assignPrecipitation(); bool setSinkSource(); void computeWaterBalance3D(double totalTimeStep); bool updateCrop(QDateTime myTime); @@ -165,7 +177,5 @@ float readDataHourly(meteoVariable myVar, QString hourlyPath, QDateTime myTime, QString myArea, int row, int col); bool readHourlyMap(meteoVariable myVar, QString hourlyPath, QDateTime myTime, QString myArea, gis::Crit3DRasterGrid* myGrid); - double getMaxEvaporation(double ET0, double LAI); - #endif // PROJECT3D_H