add raster merge option: sum

  - check input parameters for sensible values
  - changes to checkInputParameterValues() and one-line in flowCore.py

avaframe/com4FlowPy/com4FlowPy.py

@@ -167,6 +167,9 @@ def com4FlowPyMain(cfgPath, cfgSetup):
     # check if input layers have same x,y dimensions
     checkInputLayerDimensions(modelParameters, modelPaths)
 
+    # check if input parameters are within physically sensible ranges
+    checkInputParameterValues(modelParameters, modelPaths)
+
     # get information on cellsize and nodata value from demHeader
     rasterAttributes = {}
 
@@ -213,7 +216,7 @@ def startLogging(modelParameters, forestParams, modelPaths, MPOptions):
     forestParams: dict
         input parameters regarding forest interaction (from .ini - file)
     modelPaths: dict
-        contains paths to input files
+        paths to input files, workDir, resDir, outputFileFormat, etc. (from .ini - file)
     MPOptions: dict
         contains parameters for multiprocessing (from .ini - file)
     """
@@ -256,6 +259,9 @@ def startLogging(modelParameters, forestParams, modelPaths, MPOptions):
         log.info("calculation with Infrastructure")
         log.info(f"{'INFRA LAYER:' : <14}{'%s'%modelPaths['infraPath'] : <5}")
         log.info("------------------------")
+    if modelParameters["fluxDistOldVersionBool"]:
+        log.info("Calculation using old (BUGGY!!) version of flux distribution!")
+        log.info("------------------------")
     for param, value in MPOptions.items():
         log.info(f"{'%s:'%param : <20}{value : <5}")
         # log.info("{}:\t{}".format(param,value))
@@ -349,9 +355,71 @@ def checkInputLayerDimensions(modelParameters, modelPaths):
         sys.exit(1)
 
 
+def checkInputParameterValues(modelParameters, modelPaths):
+    """check if the input parameters alpha, uMaxLimit/ zDeltaMaxLimit, exponent
+    are within a physically sensible range
+
+    Parameters
+    -----------
+    modelParameters: dict
+        model input parameters (from .ini - file)
+    modelPaths: dict
+        contains paths to input files
+    """
+    alpha = modelParameters['alpha']
+    if (alpha < 0 or alpha > 90):
+        log.error("Error: Alpha value is not within a physically sensible range ([0,90]).")
+        sys.exit(1)
+
+    zDelta = modelParameters['max_z']
+    if (zDelta < 0 or zDelta > 8848):
+        log.error("Error: zDeltaMaxLimit value is not within a physically sensible range ([0,8848]).")
+        sys.exit(1)
+
+    exp = modelParameters['exp']
+    if exp < 0:
+        log.error("Error: Exponent value is not within a physically sensible range (> 0).")
+        sys.exit(1)
+
+    _checkVarParams = True
+
+    if modelParameters["varAlphaBool"]:
+        rasterValues, header = io.read_raster(modelPaths["varAlphaPath"])
+        rasterValues[rasterValues < 0] = np.nan  # handle different noData values
+        if np.any(rasterValues > 90, where=~np.isnan(rasterValues)):
+            log.error("Error: Not all Alpha-raster values are within a physically sensible range ([0,90]),\
+                 in respective startcells the general alpha angle is used.")
+            _checkVarParams = False
+
+    if modelParameters["varUmaxBool"]:
+        rasterValues, header = io.read_raster(modelPaths["varUmaxPath"])
+        rasterValues[rasterValues < 0] = np.nan
+        if modelParameters["varUmaxType"].lower() == 'umax':
+            _maxVal = 1500  # ~sqrt(8848*2*9.81)
+        else:
+            _maxVal = 8848
+        if np.any(rasterValues > _maxVal, where=~np.isnan(rasterValues)):
+            log.error("Error: Not all zDeltaMaxLimit-raster values are within a physically sensible range \
+                ([0, 8848 m] or [0, 1500 m/s]), in respective startcells the general zDeltaMax value is used.")
+            _checkVarParams = False
+
+    if _checkVarParams:
+        log.info("All input parameters are within a physically sensible range.")
+        log.info("========================")
+    else:
+        log.info("NOT ALL variable input parameter rasters are within physically sensible ranges.")
+        log.info("!!PLEASE RE-CHECK Input Rasters and handle Results with Care!!")
+        log.info("========================")
+
+
 def tileInputLayers(modelParameters, modelPaths, rasterAttributes, tilingParameters):
-    """divides all used input layers into tiles
-    and saves the tiles in the temp folder
+    """computes the number of tiles (_tileCols, _tileRows) and tileOverlap (_U)
+    based on input layer dimensions and tilingParameters,
+    divides all used input layers into tiles
+    and saves the tiles to the temp folder
+
+    the function is a wrapper around the code in splitAndMerge.py,
+    where the actual tiling is handled.
 
     Parameters
     -----------
@@ -367,9 +435,10 @@ def tileInputLayers(modelParameters, modelPaths, rasterAttributes, tilingParamet
     Returns
     -----------
     nTiles: tuple
-        number of tiles
+        nTiles[0]: maximum index of tiles along rows
+        nTiles[1]: maximum index of tiles along columns
+        actual number of tiles = (nTiles[0] + 1) * (nTiles[1] + 1)
     """
-
     _tileCOLS = int(tilingParameters["tileSize"] / rasterAttributes["cellsize"])
     _tileROWS = int(tilingParameters["tileSize"] / rasterAttributes["cellsize"])
     _U = int(tilingParameters["tileOverlap"] / rasterAttributes["cellsize"])
@@ -458,10 +527,10 @@ def mergeAndWriteResults(modelPaths, modelOptions):
     # Merge calculated tiles
     zDelta = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta")
     flux = SPAM.mergeRaster(modelPaths["tempDir"], "res_flux")
-    cellCounts = SPAM.mergeRaster(modelPaths["tempDir"], "res_count")
-    zDeltaSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta_sum")
-    routFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_rout_flux_sum")
-    depFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_dep_flux_sum")
+    cellCounts = SPAM.mergeRaster(modelPaths["tempDir"], "res_count", method='sum')
+    zDeltaSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_z_delta_sum", method='sum')
+    routFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_rout_flux_sum", method='sum')
+    depFluxSum = SPAM.mergeRaster(modelPaths["tempDir"], "res_dep_flux_sum", method='sum')
     fpTa = SPAM.mergeRaster(modelPaths["tempDir"], "res_fp")
     slTa = SPAM.mergeRaster(modelPaths["tempDir"], "res_sl")
     travelLength = SPAM.mergeRaster(modelPaths["tempDir"], "res_travel_length")
@@ -507,9 +576,6 @@ def mergeAndWriteResults(modelPaths, modelOptions):
         io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / "com4_{}_{}_travelLength{}".format(_uid,
         _ts, _oF), travelLength)
 
-    # TODO: List of result files, which are produced should be specified also in the .ini file!!!!
-    # NOTE: Probably good to have "default" output files (z_delta,FP_travel_angle,cell_counts)
-    #      and only write other output files if set accordingly
     # NOTE:
     # if not modelOptions["infraBool"]:  # if no infra
         # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("cell_counts%s" %(output_format)),cell_counts)

avaframe/com4FlowPy/com4FlowPyCfg.ini

@@ -144,11 +144,14 @@ forestFrictionLayerType = absolute
 # NOTE: this currently only works with 'forestFrictionLayer' module!!
 skipForestCells = 1 
 
-#++++++++++++ Calculate flux distribution
-# We found a bug in the distribution of the remaining flux, if a cell receives flux
-# smaller than the provided threshold (in flowClass.py). The default now is a 
-# calculation with a fixed bug. The old version (with minor bug) can be switched on 
-# with the flag.
+#++++++++++++ Method to calculate flux distribution
+# We fixed a bug in flowClass.py, which affects the distribution of the remaining flux, 
+# if a cell receives flux smaller than the provided flux_threshold. 
+#
+# The default now (post Jan. 2025) is a calculation with the fixed bug!
+# 
+# For backward compatibility the old version (prior to Jan. 2025 - with minor bug) can 
+# be switched on by setting "fluxDistOldVersion = True".
 
 fluxDistOldVersion = False
 

avaframe/com4FlowPy/flowClass.py

@@ -167,7 +167,7 @@ def add_parent(self, parent):
         Parameters
         -----------
         parent: class Cell
-            added parent
+            Cell object to add to the list of parent cells 'lOfParents'
         """
         self.lOfParents.append(parent)
         if self.forestInteraction:
@@ -191,7 +191,11 @@ def calc_fp_travelangle(self):
         self.max_gamma = np.rad2deg(np.arctan(_dh / self.min_distance))
 
     def calc_sl_travelangle(self):
-        """ function calculates the travel-angle between the start cell and the current cell.
+        """function calculates the travel-angle between the start cell and the current cell
+        using the shortest connection or straight line (sl) between the two cells.
+        The travel angle calculated with this shortest horizontal distance 'sl_gamma' is always
+        larger or equal to the travel angle  'max_gamma', which is calculated along the
+        (potentially curved) flow path (fp). (vgl. Meißl, 1998)
         """
         _dx = abs(self.startcell.colindex - self.colindex)
         _dy = abs(self.startcell.rowindex - self.rowindex)
@@ -202,7 +206,7 @@ def calc_sl_travelangle(self):
 
     def calc_z_delta(self):
         """
-        function calculates zDelta (velocity line height) to the eligible neighbours
+        function calculates zDelta (velocity or energy line height) to the eligible neighbours
         NOTE: forestFriction related mechanics are implemented here!
         """
         self.z_delta_neighbour = np.zeros((3, 3))
@@ -367,20 +371,32 @@ def calc_distribution(self):
         # "0<n<8" AND "m=0" is not handled!!! (in this case flux is "lost")
 
         if self.fluxDistOldVersionBool:
+            """
+            legacy version of code (pre 01-2025) WITH BUG (is used if self.fluxDistOldVersionBool)
+            here 'count' returns the number of neighbor/child cells that receive flux > 0, but
+            below the flux_threshold.
+            """
             count = ((0 < self.dist) & (self.dist < threshold)).sum()
         else:
-            count = (self.dist >= threshold).sum()  # this is the correct way to calculate count
-        # TODO: make this the default, but keep option to use "old" version with minor Bug for backward compatibility of
-        # model results
+            """
+            default/correct version with FIXED BUG (used if self.fluxDistOldVersionBool == False)
+            her 'count' returns the number of neighbor/child cells which receive
+            flux >= the flux_threshold.
+            """
+            count = (self.dist >= threshold).sum()
+        # massToDistribute ~ sum of flux of child cells below the flux_threshold
         mass_to_distribute = np.sum(self.dist[self.dist < threshold])
-        """Checking if flux is distributed to a field that isn't taking in account, when then distribute it equally to
-         the other fields"""
+
         if mass_to_distribute > 0 and count > 0:
+            # local flux redistribution to eligible child cells
             self.dist[self.dist > threshold] += mass_to_distribute / count
             self.dist[self.dist < threshold] = 0
         if np.sum(self.dist) < self.flux and count > 0:
+            # correction/flux conservation for potential rounding losses
             self.dist[self.dist > threshold] += (self.flux - np.sum(self.dist)) / count
-        if count == 0:  # TODO: not 100% sure if this is right
+        if count == 0:
+            # if all child cells are below flux_threshold, the flux is deposited
+            # TODO: check alternatives (e.g. 'global' redistribution or within generations)
             self.fluxDep = self.flux
         row_local, col_local = np.where(self.dist > threshold)
 

avaframe/com4FlowPy/flowCore.py

@@ -444,7 +444,7 @@ def calculation(args):
     countArray = np.zeros_like(dem, dtype=np.int32)
 
     fpTravelAngleArray = np.zeros_like(dem, dtype=np.float32)  # fp = Flow Path
-    slTravelAngleArray = np.zeros_like(dem, dtype=np.float32) * 90  # sl = Straight Line
+    slTravelAngleArray = np.zeros_like(dem, dtype=np.float32)  # sl = Straight Line
 
     travelLengthArray = np.zeros_like(dem, dtype=np.float32)
 

avaframe/com4FlowPy/rasterIo.py

@@ -14,17 +14,18 @@
 
 def read_header(input_file):
     """
-    Reads in the header of the raster file
+    Reads the header of the raster file
+    raster file should be readable by rasterio (e.g. .tif, .asc)
 
     Parameters
     -----------
     input_file: str
-        directory to raster file
+        path to raster file
 
     Returns
     -----------
     header: dict
-        header of raster
+        header of raster file in style of ASCII-Rasters
     """
 
     raster = rasterio.open(input_file)
@@ -49,14 +50,14 @@ def read_raster(input_file):
     Parameters
     -----------
     input_file: str
-        directory to raster file
+        path to raster file
 
     Returns
     -----------
     my_array: np.array
-        raster that is read in
+        numpy array with values read in from the raster file
     header: dict
-        header of raster
+        header of raster file in style of ASCII-Rasters
     """
 
     header = read_header(input_file)
@@ -66,21 +67,21 @@ def read_raster(input_file):
     return my_array, header
 
 
-def output_raster(file, file_out, raster):
+def output_raster(referenceFile, file_out, raster):
     """
     Saves raster
 
     Parameters
     -----------
-    file: str
-        directory to raster file to reference on, mostly DEM
+    referenceFile: str
+        path to raster file to reference on, mostly DEM
     file_out: str
         path for the outputfile, possible extends are .asc or .tif
     raster: np.array
         raster (array) that is saved
     """
 
-    raster_trans = rasterio.open(file)
+    raster_trans = rasterio.open(referenceFile)
     try:
         crs = rasterio.crs.CRS.from_dict(raster_trans.crs.data)
     except:

avaframe/com4FlowPy/splitAndMerge.py

@@ -3,6 +3,7 @@
 
 """
     Functions to handle the raster tiles.
+     Tiling is intended to manage processing of large(r) computational domains.
 """
 
 import logging
@@ -180,7 +181,7 @@ def mergeRaster(inDirPath, fName, method='max'):
     method provided through the function parameters
 
     Parameters
-    -----------
+    ----------
     inDirPath: str
         Path to the temporary files, that are results for each tile
     fName : str
@@ -189,9 +190,10 @@ def mergeRaster(inDirPath, fName, method='max'):
         method, how the tiles should be merged (default: max)
         method 'min' calculates the minimum of input raster tiles,
         if the minimum is < 0, then 0 is used
+        method 'sum' calculates the sum of the raster tiles
 
     Returns
-    -----------
+    -------
     mergedRas : numpy array
         merged raster
     """
@@ -204,7 +206,8 @@ def mergeRaster(inDirPath, fName, method='max'):
 
     mergedRas = np.zeros((extL[0], extL[1]))
     # create Raster with original size
-    mergedRas[:, :] = np.nan
+    if method != 'sum':
+        mergedRas[:, :] = np.nan
 
     for i in range(nTiles[0] + 1):
         for j in range(nTiles[1] + 1):
@@ -222,6 +225,10 @@ def mergeRaster(inDirPath, fName, method='max'):
                     np.where((mergedRas[pos[0][0]:pos[0][1], pos[1][0]:pos[1][1]] >= 0) & (smallRas >= 0),
                     np.fmin(mergedRas[pos[0][0]:pos[0][1], pos[1][0]:pos[1][1]], smallRas),
                     np.fmax(mergedRas[pos[0][0]:pos[0][1], pos[1][0]:pos[1][1]], smallRas))
+            if method == 'sum':
+                mergedRas[pos[0][0]: pos[0][1], pos[1][0]: pos[1][1]] = np.add(
+                    mergedRas[pos[0][0]: pos[0][1], pos[1][0]: pos[1][1]], smallRas
+                )
             del smallRas
             log.info("appended result %s_%i_%i", fName, i, j)
     return mergedRas

docs/moduleCom4FlowPy.rst

@@ -235,12 +235,12 @@ Output
 
 All outputs are in the .tif or in .asc raster format in the same resolution and extent as the input raster layers.
 
-- ``z_delta``: the maximum z_delta of all paths for every raster cell (geometric measure of process magnitude, can be associated to kinetic energy/velocity)
+- ``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
-- ``z_delta_sum``: z_delta summed up over all paths on every raster cell
-- ``cell_counts``: number of paths that route flux through a raster cell
-- ``FP_travel_angle``: the gamma angle along the flow path
-- ``SL_travel_angle``: Saves the gamma angle, while the distances are calculated via a straight line from the release cell to the current 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
+- ``fpTravelAngle``: the gamma angle along the flow path
+- ``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