[com1]: Add source line as start condition (initial)

avaframe/com1DFA/com1DFA.py

@@ -484,6 +484,8 @@ def prepareReleaseEntrainment(cfg, rel, inputSimLines):
         inputSimLines["releaseLine"]["thicknessSource"] = ["csv file"] * len(
             inputSimLines["releaseLine"]["Name"]
         )
+    elif cfg["GENERAL"]["constMassFlow"] != "":
+        inputSimLines["releaseLine"]["massFlowTot"] = cfg["GENERAL"].getfloat("constMassFlow")
     elif cfg["INPUT"]["relThFile"] == "":
         # otherwise release thickness is read from ini or shape file
         releaseLine = setThickness(cfg, inputSimLines["releaseLine"], "relTh")
@@ -642,6 +644,7 @@ def prepareInputData(inputSimFiles, cfg):
         releaseLine["thickness"] = "from raster"
         log.info("Set %s for relThField" % relRasterPath)
     # get line from release area polygon
+    # TODO:  use this for sourceline?
     if cfg["GENERAL"].getboolean("timeDependentRelease"):
         releaseLine["type"] = "time dependent Release"
         timeDepRelValues, _ = gI.getTimeDepRelCsv(cfg["INPUT"]["timeDepRelCsv"])
@@ -1188,19 +1191,43 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
         releaseLine = inputSimLines["releaseLine"]
         # create release area raster if not read from file
         if inputSimLines["releaseLine"]["initializedFrom"] == "shapefile":
-            # check if release features overlap between features
-            geoTrans.prepareArea(releaseLine, dem, thresholdPointInPoly, combine=True, checkOverlap=True)
+            # if release shp file is a line, find cells taht are crossed by the line
+            if inputSimLines["releaseLine"]["shapeTypeName"] in ["POLYLINE", "POLYLINEZ"]:
+                releaseLine = geoTrans.getCellsAlongLine(demOri["header"], releaseLine, addBuffer=False)
+                if cfg["GENERAL"]["constMassFlow"] != "":
+                    # compute release mass per raster cell and timestep
+                    releaseLine["massFlowCell"] = (
+                            releaseLine["massFlowTot"]
+                            / np.nansum(releaseLine["cellsCrossed"])
+                            * cfg["GENERAL"].getfloat("dt")
+                    )
+                    # compute thickness per release cell: mass / area / density
+                    releaseLine["rasterData"] = (
+                            releaseLine["cellsCrossed"].reshape(dem["header"]["nrows"], dem["header"]["ncols"])
+                            * releaseLine["massFlowCell"]
+                            / dem["areaRaster"]
+                            / cfg["GENERAL"].getfloat("rho")
+                    )
+                else:
+                    releaseLine["rasterData"] = (
+                            releaseLine["cellsCrossed"].reshape(dem["header"]["nrows"], dem["header"]["ncols"])
+                            * releaseLine["thickness"]
+                    )
 
-            # if no release thickness field or function - set release according to shapefile or ini file
-            # this is a list of release rasters that we want to combine
-            releaseLine = geoTrans.prepareArea(
-                releaseLine,
-                dem,
-                np.sqrt(2),
-                thList=releaseLine["thickness"],
-                combine=True,
-                checkOverlap=False,
-            )
+            else:
+                # check if release features overlap between features
+                geoTrans.prepareArea(releaseLine, dem, thresholdPointInPoly, combine=True, checkOverlap=True)
+
+                # if no release thickness field or function - set release according to shapefile or ini file
+                # this is a list of release rasters that we want to combine
+                releaseLine = geoTrans.prepareArea(
+                    releaseLine,
+                    dem,
+                    np.sqrt(2),
+                    thList=releaseLine["thickness"],
+                    combine=True,
+                    checkOverlap=False,
+                )
 
     # set relRaster
     relRaster = releaseLine["rasterData"]
@@ -1628,6 +1655,7 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel
         not cfg.getboolean("iniStep")
         and not cfg.getboolean("initialiseParticlesFromFile")
         and len(relThField) == 0
+            and not releaseLine["shapeTypeName"] in ["POLYLINE", "POLYLINEZ"]
     ):
         if debugPlot:
             xyParticlesAll = {"x": particles["x"], "y": particles["y"]}
@@ -2263,6 +2291,20 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
             particles, fields, zPartArray0 = debF.initializeTimeDepRelease(
                 cfg, inputSimLines, particles, fields, dem, zPartArray0, t
             )
+        elif inputSimLines["releaseLine"]["shapeTypeName"] in ["POLYLINE", "POLYLINEZ"]:
+            particlesRelease = com1DFA.initializeParticles(
+                cfgGen,
+                inputSimLines["releaseLine"],
+                dem,
+            )
+            particles = particleTools.mergeParticleDict(particles, particlesRelease)
+            zPartArray0 = np.append(zPartArray0, copy.deepcopy(particlesRelease["z"]))
+            particles = DFAfunC.getNeighborsC(particles, dem)
+            # update fields (compute grid values)
+            if fields["computeTA"]:
+                particles = DFAfunC.computeTrajectoryAngleC(particles, zPartArray0)
+            particles, fields = DFAfunC.updateFieldsC(cfg["GENERAL"], particles, dem, fields)
+
         # Perform computations
         particles, fields, zPartArray0, tCPU, dem = computeEulerTimeStep(
             cfgGen,
@@ -2746,6 +2788,7 @@ def computeEulerTimeStep(
     # loop version of the compute force
     log.debug("Compute Force C")
     particles, force, fields = DFAfunC.computeForceC(cfg, particles, fields, dem, frictType, resistanceType)
+
     tCPUForce = time.time() - startTime
     tCPU["timeForce"] = tCPU["timeForce"] + tCPUForce
     # compute lateral force (SPH component of the calculation)
@@ -2797,7 +2840,6 @@ def computeEulerTimeStep(
         particles, zPartArray0, reportAreaInfo = releaseSecRelArea(
             cfg, particles, fields, dem, zPartArray0, reportAreaInfo
         )
-
     # get particles location (neighbours for sph)
     startTime = time.time()
     log.debug("get Neighbours C")

avaframe/com1DFA/com1DFACfg.ini

@@ -133,8 +133,8 @@ entThDistVariation =
 entTh =
 
 #+++++++++++++General start conditions: time dependent release
-# if timeDependentRelease is True (and relThFromFile is True), provide the the timesteps, thickness and velocity
-# for a releases in a csv-file in the REL folder
+# if timeDependentRelease is True (and relThFromFile is True), provide the timesteps, thickness and velocity
+# for releases in a csv-file in the REL folder
 timeDependentRelease = False
 # specify one or multiple particular time dependent release files,
 # provide name of csv file with or without extension .csv
@@ -153,6 +153,9 @@ thresholdPointInRel = 0.01
 # (distance = (timestep[i] - timestep[i-1]) * velocity)
 timeStepDistance = 5
 
+# test for mass flow (kg/s)
+constMassFlow =
+
 #++++++++++++Time stepping parameters
 # fixed time step (also used as first time step when using CFL) [s]
 dt = 0.1

avaframe/com1DFA/com1DFATools.py

@@ -19,7 +19,7 @@
 from avaframe.com1DFA import com1DFA
 from avaframe.in1Data import getInput as gI
 from avaframe.in3Utils import cfgUtils
-from avaframe.in2Trans import rasterUtils as IOf
+from avaframe.in2Trans import shpConversion as shpConv
 
 # create local logger
 # change log level in calling module to DEBUG to see log messages
@@ -351,6 +351,8 @@ def initializeInputs(avalancheDir, cleanRemeshedRasters, module=com1DFA):
     # fetch input data - dem, release-, entrainment- and resistance areas (and secondary release areas)
     inputSimFilesAll = gI.getInputDataCom1DFA(avalancheDir)
 
+    shpConv.checkShpType(inputSimFilesAll)
+
     # get thickness of release and entrainment areas (and secondary release areas) -if thFromShp = True
     inputSimFilesAll = gI.getThicknessInputSimFiles(inputSimFilesAll)
 

avaframe/com1DFA/debrisFunctions.py

@@ -110,28 +110,35 @@ def addReleaseParticles(cfg, particles, inputSimLines, thickness, velocityMag, d
     """
     relLine = inputSimLines["releaseLine"]
     relLine["header"] = dem["originalHeader"].copy()
-    relLine = geoTrans.prepareArea(
-        relLine,
-        dem,
-        np.sqrt(2),
-        thList=[thickness] * len(relLine["Name"]),
-        combine=True,
-        checkOverlap=False,
-    )
+    if relLine["shapeTypeName"] in ["POLYLINE", "POLYLINEZ"]:
+        relLine["rasterData"] = (
+                relLine["cellsCrossed"].reshape(dem["header"]["nrows"], dem["header"]["ncols"]) * thickness
+        )
+    else:
+        relLine = geoTrans.prepareArea(
+            relLine,
+            dem,
+            np.sqrt(2),
+            thList=[thickness] * len(relLine["Name"]),
+            combine=True,
+            checkOverlap=False,
+        )
 
-    # check if already existing particles are within the release polygon
-    # it's possible that there are still a few particles in the polygon with low velocities
-    # TODO: could think of a threshold of number of particles that are still allowed in the polygons?
-    mask = geoTrans.getParticlesInPolygon(particles, relLine, cfg["GENERAL"].getfloat("thresholdPointInRel"))
-    if np.sum(mask) > 0:
-        # if there is at least one particle within the polygon (including the buffer):
-        message = (
-            "Already existing particles are within the release polygon, which can cause numerical instabilities (at timestep: %02f s)"
-            % (particles["t"] + particles["dt"])
+        # check if already existing particles are within the release polygon
+        # it's possible that there are still a few particles in the polygon with low velocities
+        # TODO: could think of a threshold of number of particles that are still allowed in the polygons?
+        mask = geoTrans.getParticlesInPolygon(
+            particles, relLine, cfg["GENERAL"].getfloat("thresholdPointInRel")
         )
-        # timestep in particles is not updated yet
-        log.error(message)
-        raise ValueError(message)
+        if np.sum(mask) > 0:
+            # if there is at least one particle within the polygon (including the buffer):
+            message = (
+                    "Already existing particles are within the release polygon, which can cause numerical instabilities (at timestep: %02f s)"
+                    % (particles["t"] + particles["dt"])
+            )
+            # timestep in particles is not updated yet
+            log.error(message)
+            raise ValueError(message)
 
     particlesRelease = com1DFA.initializeParticles(
         cfg["GENERAL"],

avaframe/in2Trans/shpConversion.py

@@ -99,6 +99,7 @@ def SHP2Array(infile, defname=None):
     # Start reading the shapefile
     records = sf.shapeRecords()
     shps = sf.shapes()
+    shapeTypeName = sf.shapeTypeName
 
     SHPdata = {}
     SHPdata["sks"] = sks
@@ -243,6 +244,7 @@ def SHP2Array(infile, defname=None):
     SHPdata["rotAngle"] = rotAngleList
     SHPdata["direc"] = direcList
     SHPdata["offset"] = offsetList
+    SHPdata["shapeTypeName"] = shapeTypeName
 
     sf.close()
 
@@ -643,3 +645,29 @@ def readShapefile(inputShp):
                 srs = f.read().strip()
 
     return fields, fieldNames, properties, geometries, srs
+
+
+def checkShpType(inputSimFiles):
+    """
+    checks if the shp type is same in all release files,
+    raise error if not
+
+    Parameters
+    ----------
+    inputSimFiles: dict
+        relFiles: contains paths to release files
+        entResInfo: contains information on file types
+    """
+    shpTypes = []
+    for releaseA in inputSimFiles["relFiles"]:
+        # fetch thickness and id info from input data
+        if inputSimFiles["entResInfo"]["relThFileType"] == ".shp":
+            sf = shapefile.Reader(str(releaseA))
+            shpTypes.append(sf.shapeType)
+
+            if len(set(shpTypes)) > 1:
+                message = (
+                    "Release shapefiles have inconsistent shape types, provide either lines or polygons."
+                )
+                log.error(message)
+                raise AssertionError(message)

avaframe/tests/test_geoTrans.py

@@ -21,7 +21,6 @@
 import avaframe.in3Utils.geoTrans as geoTrans
 import avaframe.com1DFA.DFAtools as DFAtls
 
-
 log = logging.getLogger(__name__)
 
 
@@ -1174,20 +1173,20 @@ def test_getNormalMesh(capfd):
 
         atol = 1e-10
         TestNX = np.allclose(
-            Nx[1: n - 1, 1: m - 1],
-            (-a * np.ones(np.shape(Y)) / np.sqrt(1 + a * a + b * b))[1: n - 1, 1: m - 1],
+            Nx[1 : n - 1, 1 : m - 1],
+            (-a * np.ones(np.shape(Y)) / np.sqrt(1 + a * a + b * b))[1 : n - 1, 1 : m - 1],
             atol=atol,
         )
         assert TestNX
         TestNY = np.allclose(
-            Ny[1: n - 1, 1: m - 1],
-            (-b * np.ones(np.shape(Y)) / np.sqrt(1 + a * a + b * b))[1: n - 1, 1: m - 1],
+            Ny[1 : n - 1, 1 : m - 1],
+            (-b * np.ones(np.shape(Y)) / np.sqrt(1 + a * a + b * b))[1 : n - 1, 1 : m - 1],
             atol=atol,
         )
         assert TestNY
         TestNZ = np.allclose(
-            Nz[1: n - 1, 1: m - 1],
-            (np.ones(np.shape(Y)) / np.sqrt(1 + a * a + b * b))[1: n - 1, 1: m - 1],
+            Nz[1 : n - 1, 1 : m - 1],
+            (np.ones(np.shape(Y)) / np.sqrt(1 + a * a + b * b))[1 : n - 1, 1 : m - 1],
             atol=atol,
         )
         assert TestNZ
@@ -1207,20 +1206,58 @@ def test_getNormalMesh(capfd):
         #        print((1 / np.sqrt(1 + 4*a*a*X*X + 4*b*b*Y*Y))[1:n-1, 1:m-1])
         atol = 1e-10
         TestNX = np.allclose(
-            Nx[1: n - 1, 1: m - 1],
-            (-2 * a * X / np.sqrt(1 + 4 * a * a * X * X + 4 * b * b * Y * Y))[1: n - 1, 1: m - 1],
+            Nx[1 : n - 1, 1 : m - 1],
+            (-2 * a * X / np.sqrt(1 + 4 * a * a * X * X + 4 * b * b * Y * Y))[1 : n - 1, 1 : m - 1],
             atol=atol,
         )
         assert TestNX
         TestNY = np.allclose(
-            Ny[1: n - 1, 1: m - 1],
-            (-2 * b * Y / np.sqrt(1 + 4 * a * a * X * X + 4 * b * b * Y * Y))[1: n - 1, 1: m - 1],
+            Ny[1 : n - 1, 1 : m - 1],
+            (-2 * b * Y / np.sqrt(1 + 4 * a * a * X * X + 4 * b * b * Y * Y))[1 : n - 1, 1 : m - 1],
             atol=atol,
         )
         assert TestNY
         TestNZ = np.allclose(
-            Nz[1: n - 1, 1: m - 1],
-            (1 / np.sqrt(1 + 4 * a * a * X * X + 4 * b * b * Y * Y))[1: n - 1, 1: m - 1],
+            Nz[1 : n - 1, 1 : m - 1],
+            (1 / np.sqrt(1 + 4 * a * a * X * X + 4 * b * b * Y * Y))[1 : n - 1, 1 : m - 1],
             atol=atol,
         )
         assert TestNZ
+
+
+def test_getCellsAlongLine():
+    """test for a straight line"""
+    header = {}
+    header["ncols"] = 10
+    header["nrows"] = 10
+    header["cellsize"] = 10
+    header["xllcenter"] = 5
+    header["yllcenter"] = 5
+
+    lineDict = {}
+    lineDict["x"] = np.array([5, 5])
+    lineDict["y"] = np.array([14, 45])
+
+    testCrossedCells = np.zeros((header["nrows"], header["ncols"])).astype(int)
+    testCrossedCells[[1, 2, 3, 4], [0, 0, 0, 0]] = 1
+
+    testLine = geoTrans.getCellsAlongLine(header, lineDict, addBuffer=False)
+    crossedCells = testLine["cellsCrossed"].reshape(header["nrows"], header["ncols"])
+
+    assert np.all(testCrossedCells == crossedCells)
+    assert np.all(lineDict["x"] == testLine["x"])
+    assert np.all(lineDict["y"] == testLine["y"])
+
+    lineDict = {}
+    lineDict["x"] = np.array([15, 45])
+    lineDict["y"] = np.array([15, 35])
+
+    testCrossedCells = np.zeros((header["nrows"], header["ncols"])).astype(int)
+    testCrossedCells[[1, 1, 2, 2, 3, 3], [1, 2, 2, 3, 3, 4]] = 1
+
+    testLine = geoTrans.getCellsAlongLine(header, lineDict, addBuffer=False)
+    crossedCells = testLine["cellsCrossed"].reshape(header["nrows"], header["ncols"])
+
+    assert np.all(testCrossedCells == crossedCells)
+    assert np.all(lineDict["x"] == testLine["x"])
+    assert np.all(lineDict["y"] == testLine["y"])

docs/com1DFAAlgorithm.rst

@@ -34,7 +34,8 @@ Initialize release, entrainment and resistance areas
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 Read and check shapefiles according to the configuration (check consistency between
 what is required by the configuration file and what is available in the ``Inputs`` folder).
-Convert shapefile features (polygons) to rasters (:py:func:`in3Utils.geoTrans.prepareArea`).
+Convert shapefile features to rasters (for polygons: :py:func:`in3Utils.geoTrans.prepareArea`,
+for lines: :py:func:`in3Utils.geoTrans.getCellsAlongLine`).
 Check consistency of rasters according to the following rules:
 
   - multiple release features in the release and secondary release shapefiles
@@ -144,6 +145,9 @@ If the release is time dependent, particles are initialized in provided timestep
 the particles have this initial velocity (magnitude) in direction of the steepest descent,
 it is computed in the following function: :py:func:`com1DFA.DFAfunctionsCython.updateInitialVelocity`.
 
+If the release features are lines and the release is **not** time dependent, in every time step
+particles are initialized.
+
 
 Particle properties
 ^^^^^^^^^^^^^^^^^^^^

docs/moduleCom1DFA.rst

@@ -55,12 +55,12 @@ or GeoTIFF format, or shape files and specified for the respective input type be
 
 * **digital elevation model as raster file. The format of the DEM determines the format of the output files.**
 
-* release area scenario as (multi-) polygon shapefile OR raster file (in Inputs/REL; only shapefiles OR raster files)
-    - either polygon shapefile(s):
+* release area scenario as (multi-) polygon or line shapefile OR raster file (in Inputs/REL; only shapefiles OR raster files)
+    - either polygon or line shapefile(s):
         - the release area polygon must not contain any "holes" or inner rings
         - multiple features are allowed
         - recommended attributes are *name*, *thickness* (see :ref:`moduleCom1DFA:Release-, entrainment thickness settings`) and *ci95* (see :ref:`moduleAna4Stats:probAna - Probability maps`)
-    - or raster file(s):
+- or raster file(s):
         - cells with non-zero values define the release area
         - cell value can be read as thickness (measured normal to the slope) (see :ref:`moduleCom1DFA:Release-, entrainment thickness settings`)
         - negative values and no-data values are not allowed