[com4FlowPy] minor Bugfix, documentation update and added feature for "forestFriciton", "forestDetrainment" and "forestFrictionLayer" model options

avaframe/com4FlowPy/com4FlowPy.py

@@ -130,7 +130,8 @@ def com4FlowPyMain(cfgPath, cfgSetup):
         forestParams["velThForDetrain"] = cfgSetup.getfloat("velThForDetrain")  # float(cfgSetup["velThForDetrain"])
         # 'forestFrictionLayer' parameter
         forestParams["fFrLayerType"] = cfgSetup.get("forestFrictionLayerType")
-        forestParams["nSkipForest"] = cfgSetup.getint("skipForestCells")
+        # skipForestDist - no forest friciton effect assumed while distance along path <= skipForestDist
+        forestParams["skipForestDist"] = cfgSetup.getfloat("skipForestDist")
 
     else:
         modelPaths["forestPath"] = ""
@@ -612,8 +613,6 @@ def checkConvertReleaseShp2Tif(modelPaths):
             dem = IOf.readRaster(modelPaths["demPath"])
             demHeader = IOf.readRasterHeader(modelPaths["demPath"])
         elif ext in ['.tif', '.tiff', '.TIF', '.TIFF']:
-            #dem = io.read_raster(modelPaths["demPath"])
-            #demHeader = io.read_header(modelPaths["demPath"])
             _errMsg = "using release area in '.shp' format currently only supported in combination with '.asc' DEMs"
             log.error(_errMsg)
             raise ValueError(_errMsg)

avaframe/com4FlowPy/com4FlowPyCfg.ini

@@ -136,13 +136,12 @@ velThForDetrain   = 0
 # ['absolute', 'relative']
 forestFrictionLayerType = absolute
 
-# ['1 cell', '2 cells']
-# if value == 1 forest effect is not considered at the start cell (default behaviour)
-# if value == 2 also first cell after start cell has no forest effect
-# background --> for some processes (e.g.) rockfall it might make sense, if a first
-# acceleration phase is allowed before accounting for forest effects
-# NOTE: this currently only works with 'forestFrictionLayer' module!!
-skipForestCells = 1 
+# skip Forest Effect (added forest friction) for first x meters (calculated in 3D - XYZ)
+# should allow an initial acceleration phase of processes starting in or directly above
+# dense forests (these would in many cases stop otherwise)
+# if e.g. skipForestDist = 40, no added forestFriction will be assumed until 40 m 3D-distance
+# along the path from the startCell.
+skipForestDist = 0 
 
 #++++++++++++ Method to calculate flux distribution
 # We fixed a bug in flowClass.py, which affects the distribution of the remaining flux, 

avaframe/com4FlowPy/flowClass.py

@@ -53,6 +53,7 @@ def __init__(
         # every iteration of calc_z_delta(self)
 
         self.min_distance = 0  # minimal distance to start-cell (i.e. along shortest path) min_distance >=
+        self.minDistXYZ = 0  # minimal distance to start-cell (Actual 3D lenght, not projected!!)
         self.max_distance = 0  # NOTE: self.max_distance is never used - maybe remove!?
         self.min_gamma = 0  # NOTE: self.min_gamma (assumingly minimal travel angle to cell) never used - maybe remove!?
         self.max_gamma = 0
@@ -68,6 +69,7 @@ def __init__(
 
             self.forestBool = True
             self.forestModule = forestParams["forestModule"]
+            self.skipForestDist = forestParams["skipForestDist"]
 
             # forestInteraction:
             self.forestInteraction = forestParams["forestInteraction"]
@@ -91,7 +93,7 @@ def __init__(
                 _vThFr = self.noFrictionEffectV
                 _vThDe = self.noDetrainmentEffectV
                 _sqrt2xG = self._SQRT2 * 9.81
-                self.noFricitonEffectZdelta = (_vThFr * _vThFr) / _sqrt2xG
+                self.noFrictionEffectZDelta = (_vThFr * _vThFr) / _sqrt2xG
                 self.noDetrainmentEffectZdelta = (_vThDe * _vThDe) / _sqrt2xG
 
             elif self.forestModule == "forestFrictionLayer":
@@ -103,7 +105,6 @@ def __init__(
 
                 self.AlphaFor = max(self.AlphaFor, self.alpha)  # Friction in Forest can't be lower than without forest
                 self.tanAlphaFor = np.tan(np.deg2rad(self.AlphaFor))
-                self.nSkipForestCells = forestParams["nSkipForest"]
 
             # NOTE: This is a quick hack to check if all values for Detrainment are set to 0 (as provided in the
             #      .ini file)
@@ -132,7 +133,7 @@ def __init__(
             self.maxAddedDetrainmentForest = 0
             self.minAddedDetrainmentForest = 0
             self.noDetrainmentEffectV = 0
-            self.noFricitonEffectZdelta = 0
+            self.noFrictionEffectZDelta = 0
             self.noDetrainmentEffectZdelta = 0
 
         if type(startcell) == bool:  # if a boolean variable (i.e.'True') is passed to the constructor
@@ -176,18 +177,41 @@ def add_parent(self, parent):
             if parent.forestIntCount < (self.forestIntCount - self.isForest):
                 self.forestIntCount = parent.forestIntCount + self.isForest
 
+    def calcDistMin(self, calc3D=False):
+        """
+        function calculates the projected horizontal (self.min_distance) and 3D (self.minDistXYZ) length
+        of the shortest flow path from the start-cell to the current cell.
+
+        Parameters
+        -----------
+        calc3D: bool
+            if True, the 3D distance is computed additionally (default: False)
+        """
+        if calc3D:
+            _ldistMin = []
+            _lDistMinXYZ = []
+            for parent in self.lOfParents:
+                _dx = abs(parent.colindex - self.colindex) * self.cellsize
+                _dy = abs(parent.rowindex - self.rowindex) * self.cellsize
+                _dz = abs(parent.altitude - self.altitude)
+                _ldistMin.append(math.sqrt(_dx * _dx + _dy * _dy) + parent.min_distance)
+                _lDistMinXYZ.append(math.sqrt(_dy*_dy + _dy*_dy + _dz*_dz) + parent.minDistXYZ)
+            self.min_distance = np.amin(_ldistMin)
+            self.minDistXYZ = np.amin(_lDistMinXYZ)
+        else:
+            _ldistMin = []
+            for parent in self.lOfParents:
+                _dx = abs(parent.colindex - self.colindex) * self.cellsize
+                _dy = abs(parent.rowindex - self.rowindex) * self.cellsize
+                _ldistMin.append(math.sqrt(_dx * _dx + _dy * _dy) + parent.min_distance)
+            self.min_distance = np.amin(_ldistMin)
+
     def calc_fp_travelangle(self):
         """function calculates the travel-angle along the shortest flow-path from the start-cell
         to the current cell. The travel-angle along the shortest flow-path is equivalent to the
         maximum travel angle along all paths from the startcell to this cell.
         """
-        _ldistMin = []  #
         _dh = self.startcell.altitude - self.altitude  # elevation difference from cell to start-cell
-        for parent in self.lOfParents:
-            _dx = abs(parent.colindex - self.colindex)
-            _dy = abs(parent.rowindex - self.rowindex)
-            _ldistMin.append(math.sqrt(_dx * _dx + _dy * _dy) * self.cellsize + parent.min_distance)
-        self.min_distance = np.amin(_ldistMin)
         self.max_gamma = np.rad2deg(np.arctan(_dh / self.min_distance))
 
     def calc_sl_travelangle(self):
@@ -213,41 +237,40 @@ def calc_z_delta(self):
         self.z_gamma = self.altitude - self.dem_ng
         ds = np.array([[self._SQRT2, 1, self._SQRT2], [1, 0, 1], [self._SQRT2, 1, self._SQRT2]])
 
+        if (not self.is_start):
+            if (not self.forestBool):
+                self.calcDistMin()
+            else:
+                self.calcDistMin(calc3D=True)
+
         if self.forestBool:
+
             if self.forestModule == "forestFrictionLayer":
-                # default behavior - forest effect only neglected for start-cells
-                if (self.nSkipForestCells == 1) and (not self.is_start):
-                    _tanAlpha = self.tanAlphaFor
-                # forest effect also neglected for direct successors to the start-cell if nSkipForestCells==2
-                elif ((self.nSkipForestCells == 2) and (not self.is_start) and
-                      (True not in [x.is_start for x in self.lOfParents])):
+                if (not self.is_start) and (self.skipForestDist < self.minDistXYZ):
                     _tanAlpha = self.tanAlphaFor
                 else:
                     _tanAlpha = self.tanAlpha
 
-            elif (self.forestModule == "forestFriction") or (self.forestModule == "forestDetrainment"):
-
-                if (self.forestBool) and (self.FSI > 0.0) and (not self.is_start):
+            if self.forestModule in ["forestFriction", "forestDetrainment"]:
+                if (not self.is_start) and (self.FSI > 0.) and (self.skipForestDist < self.minDistXYZ):
                     # if forestBool, we assume that forestFriciton is activated
                     # and if FSI > 0 then we also calculate _tanAlpha with forestEffect
                     # NOTE: We also don't assume a forest Effect on potential Start Zells, since this should
                     #      ideally be handled by a separate release-area algorithm in the pre-processing
                     # NOTE-TODO: The rest of this implementation is also just copy+pasted from 'foreste_detraiment'
                     #      branch and not yet fully tested!!
-                    if self.z_delta < self.noFricitonEffectZdelta:
+                    if self.z_delta < self.noFrictionEffectZDelta:
                         # friction at rest v=0 would be applied to start cells
                         _rest = self.maxAddedFrictionForest * self.FSI
                         # rise over run
-                        _slope = (_rest - self.minAddedFrictionForest) / (0 - self.noFricitonEffectZdelta)
+                        _slope = (_rest - self.minAddedFrictionForest) / (0 - self.noFrictionEffectZDelta)
                         # y = mx + b, shere z_delta is the x
                         friction = max(self.minAddedFrictionForest, _slope * self.z_delta + _rest)
-
                         _alpha_calc = self.alpha + max(0, friction)  # NOTE: not sure what this does, seems redundant!
                     else:
                         _alpha_calc = self.alpha + self.minAddedFrictionForest
 
                     _tanAlpha = np.tan(np.deg2rad(_alpha_calc))
-
                 else:
                     _tanAlpha = self.tanAlpha
 

avaframe/com4FlowPy/flowCore.py

@@ -383,8 +383,8 @@ def calculation(args):
         the maximum of flux in every cell
     countArray: numpy array
         the number of hits (GMF paths) in every cell
-    zDeltaPathList: list
-        containing the max zDelta Arrays of all paths
+    zDeltaSumArray: numpy array
+        the maximum of zDelta in every cell per path and the sum over the paths
     backcalc: numpy array
           Array with back calculation, still TODO!!!
     fpTravelAngleArray: numpy array

avaframe/runCom4FlowPy.py

@@ -60,7 +60,7 @@ def main():
         initProj.cleanModuleFiles(avalancheDir, com4FlowPy, deleteOutput=False)
 
         # Start logging
-        log = logUtils.initiateLogger(avalancheDir, logName)
+        log = logUtils.initiateLogger(avalancheDir, logName+'_'+uid)
         log.info("==================================")
         log.info("MAIN SCRIPT")
         log.info("Current avalanche: %s", avalancheDir)
@@ -109,19 +109,26 @@ def main():
     # if customPaths == True --> check
     elif cfgCustomPaths["useCustomPaths"] == "True":
         # if "useCustomPaths" == True, we don't need the AvaDir Info for the
-        # creation of the simulaiton uid
+        # creation of the simulation uid
         uid = cfgUtils.cfgHash(cfg)
         cfgPath = {}
 
         # Handling Custom directory creation
         workDir = pathlib.Path(cfgCustomPaths["workDir"])
 
+        if not os.path.isdir(workDir):
+            try:
+                os.makedirs(workDir)
+            except Exception as e:
+                return e
+
+        log = logUtils.initiateLogger(workDir, logName+'_'+uid)
+
         timeString = datetime.now().strftime("%Y%m%d_%H%M%S")
         try:
             os.makedirs(workDir / "res_{}".format(uid))  # (time_string))
             res_dir = workDir / "res_{}".format(uid)   # (time_string)
         except FileExistsError:
-            log.info("folder with same name already exists - aborting")
             log.info("simulation results folder with same .ini parameters already exists: simulation {}".format(uid))
             sys.exit(1)
         try:
@@ -130,7 +137,6 @@ def main():
         except FileExistsError:
             log.info("temp folder for simualtion {} already exists - aborting".format(uid))
             sys.exit(1)
-        log = logUtils.initiateLogger(res_dir, logName)
 
         # writing config to .json file
         successToJSON = writeCfgJSON(cfg, uid, workDir)

docs/moduleCom4FlowPy.rst

@@ -233,22 +233,27 @@ deleted after completion of the model run (can be useful for calculation of larg
 Output
 -------
 
-All outputs are in the .tif or in .asc raster format in the same resolution and extent as the input raster layers.
+All outputs are written in *'.tif'* or in *'.asc'* raster format (controlable via the ``outputFileFormat`` option in ``(local_)com4FlowPyCfg.ini``, default is *'.tif'*) in the same resolution and extent as the input raster layers.
+You can customize which output rasters are written at the end of the model run by selecting the desired output files through the ``outputFiles`` option in ``(local_)com4FlowPyCfg.ini``.
+
+By default the following four output layers are written to disk at the end of the model run:
 
 - ``zdelta``: the maximum z_delta of all paths for every raster cell (geometric measure of process magnitude, can be associated to kinetic energy/velocity)
-- ``flux``: The maximum routing flux of all paths for every raster cell
-- ``zDeltaSum``: z_delta summed up over all paths on every raster cell
 - ``cellCounts``: number of paths/release cells that route flux through a raster cell
+- ``travelLength``: the travel length along the flow path
 - ``fpTravelAngle``: the gamma angle along the flow path
+
+In addition these output layers are also available:
+
+- ``flux``: The maximum routing flux of all paths for every raster cell
+- ``zDeltaSum``: z_delta summed up over all paths on every raster cell
 - ``slTravelAngle``: gamma angle calculated along a straight-line between release cell and current cell
-- ``travelLength``: the travel length along the flow path
 - ``routFluxSum``: routing flux summed up over all paths
 - ``depFluxSum``: deposited flux summed up over all paths
-- ``forestInteraction``: only if ``forestInteraction = True``: minimum number of forested raster cells a path runs through
 
-.. Note::
-  * **please interpret** ``cell_counts.tif`` **with caution, since absolute cell_count values do currently not reflect the number of release-cells which route flux through a cell - we are currently fixing the implementation of this feature**
-  * we are also working on making the output files configurable via the ``com4FlowPyCfg.ini`` file for improved flexibility (different output files might be desirable for different applications)
+If ``forestInteraction = True`` this layer will be written automatically (no need to separately define in ``outputFiles``):
+
+- ``forestInteraction``: minimum number of forested raster cells a path runs through
 
  .. Model Parameterisation
  .. ------------------------

docs/theoryCom4FlowPy.rst

@@ -12,6 +12,7 @@ com4FlowPy theory
 
   Initial information pertaining to the application of :py:mod:`com4FlowPy` with consideration of forest-effects can be found in:
      - D'Amboise, C.J.L; Teich, M.; Hormes, A.; Steger, S. and Berger, F. (2021). Modeling Protective Forests for Gravitational Natural Hazards and How It Relates to Risk-Based Decision Support Tools. in: *Protective forests as Ecosystem-based solution for Disaster Risk Reduction (ECO-DRR), Teich et al., 2021*. (https://www.intechopen.com/chapters/78979)
+     - Huber, A.; Saxer, L.; Hesselbach, C.; Neuhauser, M.; D'Amboise, C.J.L. and Teich, M. (2024): Regional-scale avalanche modeling with com4FlowPy - potential and limitations for considering avalanche-forest interaction along the avalanche track. *Proceedings, International Snow Science Workshop, Tromso, Norway* (https://arc.lib.montana.edu/snow-science/item.php?id=3194)
 
 
 Background and motivation
@@ -288,13 +289,16 @@ autoATES
 ~~~~~~~~~~~~~~~~~~~
 - Toft, H.B.; Sykes, J.; Schauer, A.; Hendrikx, J. and Hetland, A. (2024). AutoATES v2.0: Automated Avalanche Terrain Exposure Scale mapping. *Nat. Hazards Earth Syst. Sci., 24*, 1779–1793 (https://doi.org/10.5194/nhess-24-1779-2024)
 - Sykes, J.; Toft, H.B.; Haegeli, P. and Statham, G. (2024). Automated Avalanche Terrain Exposure Scale (ATES) mapping – local validation and optimization in western Canada. *Nat. Hazards Earth Syst. Sci., 24*, 947–971 (https://doi.org/10.5194/nhess-24-947-2024)
+- Panayotov, M.; Markov, K.; Tsvetanov, N.; Huber, A.; Hesselbach, C. and Teich, M. (2024). Avalanche Hazard Mapping using the Avalanche Terrain Exposure Scale (ATES) in the high mountain ranges of Bulgaria. *Proceedings, International Snow Science Workshop, Tromso, Norway* (https://arc.lib.montana.edu/snow-science/item.php?id=3366)
 - Spannring, P.; Hesselbach, C.; Mitterer, C. and Fischer, J.-T. (2024). Classification of avalanche terrain: a open-source model chain for the Avalanche Terrain Exposure Scale. *Proceedings Interpraevent 2024, Vienna, Austria* (https://www.interpraevent.at/en/proceeding/proceedings-ip-2024)
+- Spannring, P. (2024). Comparison of two avalanche terrain classification approaches : Avalanche Terrain Exposure scale - Classified Avalanche Terrain. *Masters' Thesis*. University of Innsbruck. (https://ulb-dok.uibk.ac.at/urn/urn:nbn:at:at-ubi:1-155858)
 - von Avis, C.D.; Sykes, J. and Tutt, B. (2023). Development of large scale automated Avalanche Terrain Exposure Scale (ATES) ratings in collaboration with local avalanche experts. *Proceedings, International Snow Science Workshop, Bend, OR, USA* (http://arc.lib.montana.edu/snow-science/item/2998)
 - Huber, A.; Hesselbach, C.; Oesterle, F.; Neuhauser, M.; Adams, M.; Plörer, M.; Stephan, L.; Toft, H.B.; Sykes, J.; Mitterer, C. and Fischer, J.-T. (2023). AutoATES Austria - Testing and application of an automated model-chain for avalanche terrain classification in the Austrian Alps. *Proceedings, International Snow Science Workshop, Bend, OR, USA* (http://arc.lib.montana.edu/snow-science/item/2989)
 - Hesselbach, C. (2023). Adaptaion and application of an automated Avalanche Terrain Classification in Austria. *Masters' Thesis*. University of Life Sciences (BOKU), Vienna (https://forschung.boku.ac.at/de/publications/175549)
 - Schumacher, J.; Toft, H.B.; McLean, J.P.; Hauglin, M.; Astrup, R. and Breidenbach, J. (2022). The utility of forest attribute maps for automated Avalanche Terrain Exposure Scale (ATES) modelling. Scandinavian Journal of Forest Research, 37:4, 264-275 (https://doi.org/10.1080/02827581.2022.2096921)
 
 other
 ~~~~~~~~~~~~~~~~~~~
+- Huber, A.; Saxer, L.; Hesselbach, C.; Neuhauser, M.; D'Amboise, C.J.L. and Teich, M. (2024): Regional-scale avalanche modeling with com4FlowPy - potential and limitations for considering avalanche-forest interaction along the avalanche track. *Proceedings, International Snow Science Workshop, Tromso, Norway* (https://arc.lib.montana.edu/snow-science/item.php?id=3194)
 - Perzl, F.; Huber, A.; Fromm, R. and Teich, M. (2024). Estimation of potential snow avalanche hazard probability in areas below protective forests in Austria. *Proceedings Interpraevent 2024, Vienna, Austria* (https://www.interpraevent.at/en/proceeding/proceedings-ip-2024)
 - D'Amboise, C.J.L; Teich, M.; Hormes, A.; Steger, S. and Berger, F. (2021). Modeling Protective Forests for Gravitational Natural Hazards and How It Relates to Risk-Based Decision Support Tools. in: *Protective forests as Ecosystem-based solution for Disaster Risk Reduction (ECO-DRR), Teich et al., 2021*. (https://www.intechopen.com/chapters/78979)