root density nel 3D

ARPA-SIMC · Dec 7, 2023 · 3d40beb · 3d40beb
1 parent 63efa74
commit 3d40beb
Show file tree

Hide file tree

Showing 5 changed files with 138 additions and 2 deletions.
diff --git a/agrolib/crop/root.cpp b/agrolib/crop/root.cpp
@@ -492,5 +492,132 @@ namespace root
         return true;
     }
 
+
+    bool computeRootDensity3D(Crit3DCrop* myCrop, const soil::Crit3DSoil &currentSoil, int nrLayers,
+                              const std::vector<double> &layerDepth, const std::vector<double> &layerThickness)
+    {
+        // check soil
+        if (nrLayers == 0)
+        {
+            myCrop->roots.firstRootLayer = NODATA;
+            myCrop->roots.lastRootLayer = NODATA;
+            return false;
+        }
+
+        // Initialize
+        myCrop->roots.rootDensity.clear();
+         myCrop->roots.rootDensity.resize(nrLayers);
+        for (unsigned int i = 0; i < nrLayers; i++)
+        {
+            myCrop->roots.rootDensity[i] = 0.0;
+        }
+
+        if (myCrop->roots.currentRootLength <= 0 )
+            return true;
+
+        if (myCrop->roots.rootShape == GAMMA_DISTRIBUTION)
+            myCrop->roots.rootShape = CARDIOID_DISTRIBUTION;
+
+        int nrAtoms = int(currentSoil.totalDepth * 100) + 1;
+        double minimumThickness = 0.01;                                    // [m]
+
+        int numberOfRootedLayers, numberOfTopUnrootedLayers;
+        numberOfTopUnrootedLayers = int(round(myCrop->roots.rootDepthMin / minimumThickness));
+        numberOfRootedLayers = int(round(MINVALUE(myCrop->roots.currentRootLength, currentSoil.totalDepth) / minimumThickness));
+
+        // roots are still too short
+        if (numberOfRootedLayers == 0)
+            return true;
+
+        // check nr of thin layers
+        if ((numberOfTopUnrootedLayers + numberOfRootedLayers) > nrAtoms)
+        {
+            numberOfRootedLayers = nrAtoms - numberOfTopUnrootedLayers;
+        }
+
+        // initialize thin layers density
+        std::vector<double> densityThinLayers;
+        densityThinLayers.resize(nrAtoms);
+        for (int i=0; i < nrAtoms; i++)
+        {
+            densityThinLayers[i] = 0.;
+        }
+
+        if (myCrop->roots.rootShape == CARDIOID_DISTRIBUTION)
+        {
+            cardioidDistribution(myCrop->roots.shapeDeformation, numberOfRootedLayers,
+                                 numberOfTopUnrootedLayers, signed(nrAtoms), densityThinLayers);
+        }
+        else if (myCrop->roots.rootShape == CYLINDRICAL_DISTRIBUTION)
+        {
+            cylindricalDistribution(myCrop->roots.shapeDeformation, numberOfRootedLayers,
+                                    numberOfTopUnrootedLayers, signed(nrAtoms), densityThinLayers);
+        }
+
+
+        double maxLayerDepth = layerDepth[nrLayers-1] + layerThickness[nrLayers-1] * 0.5;
+        int atom = 0;
+        double currentDepth = atom / 100.;                                    // [m]
+        double rootDensitySum = 0.;
+        while (currentDepth <= maxLayerDepth && atom < nrAtoms)
+        {
+            for (int l = 0; l < nrLayers; l++)
+            {
+                double upperDepth = layerDepth[l] - layerThickness[l] * 0.5;
+                double lowerDepth = layerDepth[l] + layerThickness[l] * 0.5;
+                if (currentDepth >= upperDepth && currentDepth <= lowerDepth)
+                {
+                    myCrop->roots.rootDensity[l] += densityThinLayers[atom];
+                    rootDensitySum += densityThinLayers[atom];
+                    break;
+                }
+            }
+
+            atom++;
+            currentDepth = atom / 100.;                                    // [m]
+        }
+
+        if (rootDensitySum <= EPSILON)
+            return true;
+
+        double rootDensitySumSubset = 0.;
+        for (unsigned int l=0 ; l < nrLayers; l++)
+        {
+            int horIndex = currentSoil.getHorizonIndex(layerDepth[l]);
+            if (horIndex != int(NODATA))
+            {
+                double soilFraction = (1 - currentSoil.horizon[horIndex].coarseFragments);
+                myCrop->roots.rootDensity[l] *= soilFraction;
+                rootDensitySumSubset += myCrop->roots.rootDensity[l];
+            }
+        }
+
+        if (rootDensitySumSubset != rootDensitySum)
+        {
+            double ratio = rootDensitySum / rootDensitySumSubset;
+            for (unsigned int l=0 ; l < nrLayers; l++)
+            {
+                myCrop->roots.rootDensity[l] *= ratio;
+            }
+        }
+
+        myCrop->roots.firstRootLayer = 0;
+        unsigned int layer = 0;
+        while (layer < nrLayers && myCrop->roots.rootDensity[layer] == 0.0)
+        {
+            layer++;
+            (myCrop->roots.firstRootLayer)++;
+        }
+
+        myCrop->roots.lastRootLayer = myCrop->roots.firstRootLayer;
+        while (layer < nrLayers && myCrop->roots.rootDensity[layer] != 0.0)
+        {
+            myCrop->roots.lastRootLayer = signed(layer);
+            layer++;
+        }
+
+        return true;
+    }
+
 }
 
diff --git a/agrolib/crop/root.h b/agrolib/crop/root.h
@@ -54,6 +54,8 @@
 
         double getRootLengthDD(const Crit3DRoot &myRoot, double currentDD, double emergenceDD);
         bool computeRootDensity(Crit3DCrop* myCrop, const std::vector<soil::Crit3DLayer> &soilLayers);
+        bool computeRootDensity3D(Crit3DCrop* myCrop, const soil::Crit3DSoil &currentSoil, int nrLayers,
+                                  const std::vector<double> &layerDepth, const std::vector<double> &layerThickness);
 
         int getNrAtoms(const std::vector<soil::Crit3DLayer> &soilLayers, double &minThickness, std::vector<int> &atoms);
     }

diff --git a/agrolib/soil/soil.cpp b/agrolib/soil/soil.cpp
@@ -188,7 +188,7 @@ namespace soil
         nrHorizons = nrHorizons - 1;
     }
 
-    int Crit3DSoil::getHorizonIndex(double depth)
+    int Crit3DSoil::getHorizonIndex(double depth) const
     {
        for (unsigned int index = 0; index < nrHorizons; index++)
        {

diff --git a/agrolib/soil/soil.h b/agrolib/soil/soil.h
@@ -201,7 +201,7 @@
             void cleanSoil();
             void addHorizon(int nHorizon, const Crit3DHorizon &newHorizon);
             void deleteHorizon(int nHorizon);
-            int getHorizonIndex(double depth);
+            int getHorizonIndex(double depth) const;
 
             bool setSoilLayers(double layerThicknessMin, double geometricFactor,
                                std::vector<Crit3DLayer> &soilLayers, std::string &myError);

diff --git a/bin/CRITERIA3D/shared/project3D.cpp b/bin/CRITERIA3D/shared/project3D.cpp
@@ -32,6 +32,7 @@
 
 #include "math.h"
 #include "utilities.h"
+#include "root.h"
 
 #include <QUuid>
 #include <QApplication>
@@ -1388,6 +1389,12 @@ double Project3D::assignTranspiration(int row, int col, double lai, double degre
     if (cropList[cropIndex].roots.currentRootLength <= 0)
         return actualTranspiration;
 
+    // update root density
+    if (! root::computeRootDensity3D(&(cropList[cropIndex]), soilList[soilIndex], nrLayers, layerDepth, layerThickness))
+    {
+        return actualTranspiration;
+    }
+
     // TODO assign root density from degreedays
     /*
     if (roots.firstRootLayer == NODATA) return 0;