[com1]: time dependent release file with x and y cooridnates and velocity components

avaframe/com1DFA/com1DFA.py

@@ -620,7 +620,17 @@ def prepareInputData(inputSimFiles, cfg):
     # TODO: remove if not required anymore
     # relThFieldData, _ = gI.initializeRelTh(cfg, dOHeader)
 
-    if cfg["INPUT"]["relThFile"] == "":
+    if relFile.suffix.lower() == ".csv":
+        if not cfg["GENERAL"].getboolean("timeDependentRelease"):
+            message = "A CSV release geometry requires timeDependentRelease=True"
+            log.error(message)
+            raise ValueError(message)
+        timeDepRelValues, _ = gI.getTimeDepRelCsv(cfg["INPUT"]["timeDepRelCsv"])
+        if "x" not in timeDepRelValues:
+            message = "A CSV release geometry requires x and y columns in the time dependent release CSV"
+            log.error(message)
+            raise ValueError(message)
+    elif cfg["INPUT"]["relThFile"] == "":
         # get line from release area polygon
         releaseLine = shpConv.readLine(relFile, "release1", demOri)
         releaseLine["file"] = relFile
@@ -646,13 +656,24 @@ def prepareInputData(inputSimFiles, cfg):
         log.info("Set %s for relThField" % relRasterPath)
     # get line from release area polygon
     if cfg["GENERAL"].getboolean("timeDependentRelease"):
-        releaseLine["type"] = "time dependent Release"
         timeDepRelValues, _ = gI.getTimeDepRelCsv(cfg["INPUT"]["timeDepRelCsv"])
-        releaseLine["thickness"] = [
-            timeDepRelValues["thickness"][timeDepRelValues["timeStep"] == 0].item()
-        ] * len(releaseLine["Name"])
+
+        if "x" in timeDepRelValues:
+            releaseLine, relThFieldData = debF.defineReleaseLineFromCoordinates(relFile, timeDepRelValues, dOHeader)
+        else:
+            releaseLine["thickness"] = [
+                                           timeDepRelValues["thickness"][timeDepRelValues["timeStep"] == 0].item()
+                                       ] * len(releaseLine["Name"])
+
+            if "velocity" in timeDepRelValues.keys():
+                releaseLine["velocity"] = timeDepRelValues["velocity"][timeDepRelValues["timeStep"] == 0]
+            else:
+                releaseLine["velocityX"] = timeDepRelValues["velocityX"][timeDepRelValues["timeStep"] == 0]
+                releaseLine["velocityY"] = timeDepRelValues["velocityY"][timeDepRelValues["timeStep"] == 0]
+                releaseLine["velocityZ"] = timeDepRelValues["velocityZ"][timeDepRelValues["timeStep"] == 0]
+
+        releaseLine["type"] = "time dependent Release"
         releaseLine["thicknessSource"] = ["csv file"] * len(releaseLine["Name"])
-        releaseLine["velocity"] = timeDepRelValues["velocity"][timeDepRelValues["timeStep"] == 0]
         releaseLine["timeDepRelValues"] = timeDepRelValues
 
     # get line from secondary release area polygon
@@ -1258,8 +1279,9 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
         relThField=relThField,
     )
 
-    if cfgGen.getboolean("timeDependentRelease") and releaseLine["velocity"] != 0:
-        particles = DFAfunC.updateInitialVelocity(cfgGen, particles, dem, releaseLine["velocity"])
+    if cfgGen.getboolean("timeDependentRelease") and "velocity" in releaseLine.keys():
+        if np.any(releaseLine["velocity"]) != 0:
+            particles = DFAfunC.updateInitialVelocity(cfgGen, particles, dem, releaseLine["velocity"])
     particles, fields = initializeFields(cfg, dem, particles, releaseLine)
 
     reportAreaInfo["Release area info"]["Model release volume [m3]"] = "%.0f" % (
@@ -1449,7 +1471,7 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName):
     return particles, fields, dem, reportAreaInfo
 
 
-def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", relThField="", thName="rel"):
+def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", relThField="", thName="rel", timestep=0.0):
     """Initialize DFA simulation
 
     Create particles and fields dictionary according to config parameters
@@ -1469,6 +1491,8 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel
         if the release thickness is not uniform, give here the releaseRaster
     thName: str
         name rel, secondaryRel
+    timestep: float
+        for log: timestep at which particles are released
 
     Returns
     -------
@@ -1631,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 releaseLine["initializedFrom"] != "csvfile"
     ):
         if debugPlot:
             xyParticlesAll = {"x": particles["x"], "y": particles["y"]}
@@ -1667,8 +1692,8 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel
     particles["nPPK"] = nPPK
 
     log.info(
-        "Initialized particles. MTot = %.2f kg, %s particles in %.2f cells."
-        % (particles["mTot"], particles["nPart"], relCells)
+        "Initialized particles in t = %.2f s. MTot = %.2f kg, %s particles in %.2f cells."
+        % (timestep, particles["mTot"], particles["nPart"], relCells)
     )
     log.info(
         "Mass per particle = %.2f kg and particles per cell = %.2f."
@@ -2860,7 +2885,8 @@ def releaseSecRelArea(cfg, particles, fields, dem, zPartArray0, reportAreaInfo):
             if secondaryReleaseInfo["initializedFrom"] == "shapefile":
                 secRelInfo = shpConv.extractFeature(secondaryReleaseInfo, count)
                 secRelInfo["rasterData"] = secRelRaster
-                secRelParticles = initializeParticles(cfg, secRelInfo, dem, thName="secondaryRel")
+                secRelParticles = initializeParticles(cfg, secRelInfo, dem, thName="secondaryRel",
+                                                      timestep=particles["t"])
             else:
                 secondaryReleaseInfo["rasterData"] = secRelRaster
                 secRelParticles = initializeParticles(
@@ -2869,6 +2895,7 @@ def releaseSecRelArea(cfg, particles, fields, dem, zPartArray0, reportAreaInfo):
                     dem,
                     relThField=secRelRaster,
                     thName="secondaryRel",
+                    timestep=particles["t"]
                 )
             # release secondary release area by just appending the particles
             log.info(
@@ -3436,7 +3463,12 @@ def prepareVarSimDict(standardCfg, inputSimFiles, variationDict, simNameExisting
             timeDepRelValues, _ = gI.getTimeDepRelCsv(cfgSim["INPUT"]["timeDepRelCsv"])
             cfgSim["INPUT"]["timeDepRelTimeStep"] = str(timeDepRelValues["timeStep"])
             cfgSim["INPUT"]["timeDepRelThickness"] = str(timeDepRelValues["thickness"])
-            cfgSim["INPUT"]["timeDepRelVelocity"] = str(timeDepRelValues["velocity"])
+            if "velocity" in timeDepRelValues:
+                cfgSim["INPUT"]["timeDepRelVelocity"] = str(timeDepRelValues["velocity"])
+            else:
+                cfgSim["INPUT"]["timeDepRelVelocityX"] = str(timeDepRelValues["velocityX"])
+                cfgSim["INPUT"]["timeDepRelVelocityY"] = str(timeDepRelValues["velocityY"])
+                cfgSim["INPUT"]["timeDepRelVelocityZ"] = str(timeDepRelValues["velocityZ"])
         else:
             cfgSim["INPUT"]["timeDepRelCsv"] = ""
 
@@ -3851,10 +3883,20 @@ def initializeRelVol(cfg, demVol, releaseFile, radius, releaseType="primary", ti
 
     # check if release thickness provided as field or constant value
     if cfg["INPUT"][(typeTh + "File")] != "":
-        # read relThField from file
-        relThFilePath = pathlib.Path(cfg["GENERAL"]["avalancheDir"], "Inputs", cfg["INPUT"][typeTh + "File"])
-        relThFieldFull = IOf.readRaster(relThFilePath)
-        relThField = relThFieldFull["rasterData"]
+        if releaseType == "timeDepRel":
+            # compute total initialized thickness
+            timeDepRelValues, _ = gI.getTimeDepRelCsv(timeDepRelFile)
+            # for time dependent release use the release volume summed up over all timesteps
+            relThField = np.zeros((demVol["header"]["nrows"], demVol["header"]["ncols"]))
+            for ts in np.unique(timeDepRelValues["timeStep"]):
+                idx = np.where(timeDepRelValues["timeStep"] == ts)[0]
+                thRaster = gI.timeDepRelCoordsToRaster(timeDepRelValues, idx, demVol["header"], parameter="thickness")
+                relThField += thRaster
+        else:
+            # read relThField from file
+            relThFilePath = pathlib.Path(cfg["GENERAL"]["avalancheDir"], "Inputs", cfg["INPUT"][typeTh + "File"])
+            relThFieldFull = IOf.readRaster(relThFilePath)
+            relThField = relThFieldFull["rasterData"]
 
         # mask the relThField with raster from polygon
         releaseLineMask = np.ma.masked_where(relThField == 0.0, relThField)

avaframe/com1DFA/com1DFACfg.ini

@@ -133,8 +133,9 @@ 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 timestep, thickness and velocity
+# in a csv file in Inputs/REL. Optional x and y columns define a one-cell release location
+# for each timestep; when x/y are provided, no release polygon is required.
 timeDependentRelease = False
 # specify one or multiple particular time dependent release files,
 # provide name of csv file with or without extension .csv

avaframe/com1DFA/debrisFunctions.py

@@ -56,22 +56,20 @@ def initializeTimeDepRelease(cfg, inputSimLines, particles, fields, dem, zPartAr
     if np.isclose(t, timeDepRelValues["timeStep"], atol=atol, rtol=0).any():
         iTup = np.where(np.isclose(t, timeDepRelValues["timeStep"], atol=atol, rtol=0))
         # iTup is a tuple containing an array with one value in the first position, so we can extract the index:
-        i = iTup[0].item()
-        log.info(
-            "add release at timestep: %.2f s with thickness %s m and velocity %s m/s"
-            % (t, timeDepRelValues["thickness"][i], timeDepRelValues["velocity"][i])
-        )
+        i = iTup[0]
+
         # similar workflow to secondary release!
         particles, zPartArray0 = addReleaseParticles(
             cfg,
             particles,
             inputSimLines,
-            timeDepRelValues["thickness"][i],
-            timeDepRelValues["velocity"][i],
+            timeDepRelValues,
             dem,
             zPartArray0,
+            timeDepRelIndex=i,
         )
         particles = DFAfunC.getNeighborsC(particles, dem)
+
         # update fields (compute grid values)
         if fields["computeTA"]:
             particles = DFAfunC.computeTrajectoryAngleC(particles, zPartArray0)
@@ -80,7 +78,9 @@ def initializeTimeDepRelease(cfg, inputSimLines, particles, fields, dem, zPartAr
     return particles, fields, zPartArray0
 
 
-def addReleaseParticles(cfg, particles, inputSimLines, thickness, velocityMag, dem, zPartArray0):
+def addReleaseParticles(
+        cfg, particles, inputSimLines, timeDepRelValues, dem, zPartArray0, timeDepRelIndex
+):
     """
     add new particles initialized by a time dependent release to particles that are in the flow already
 
@@ -92,10 +92,8 @@ def addReleaseParticles(cfg, particles, inputSimLines, thickness, velocityMag, d
         particles dictionary at t that are in the flow already
     inputSimLines : dict
         dictionary with input data dictionaries (releaseLine,...)
-    thickness: float
-        thickness of current release
-    velocityMag: float
-        velocity of current release
+    timeDepRelValues: dict
+        time dependent release values
     dem: dict
         dictionary with info on DEM data
     zPartArray0: numpy array
@@ -108,43 +106,106 @@ def addReleaseParticles(cfg, particles, inputSimLines, thickness, velocityMag, d
     zPartArray0: dict
         dictionary containing z - value of particles at timestep 0
     """
-    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,
-    )
+    thickness = timeDepRelValues["thickness"][timeDepRelIndex]
+    if "velocity" in timeDepRelValues:
+        velocityMag = timeDepRelValues["velocity"][timeDepRelIndex]
 
-    # 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"])
+    relLine = copy.deepcopy(inputSimLines["releaseLine"])
+    relLine["header"] = dem["originalHeader"].copy()
+    if relLine["initializedFrom"] == "csvfile":
+        relLine["rasterData"] = gI.timeDepRelCoordsToRaster(
+            relLine["timeDepRelValues"], timeDepRelIndex, dem["originalHeader"], parameter="thickness"
         )
-        # timestep in particles is not updated yet
-        log.error(message)
-        raise ValueError(message)
+        relThField = relLine["rasterData"]
+    else:
+        relThField = ""
+        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:
+            message = (
+                    "Already existing particles are within the release polygon, which can cause numerical instabilities (at timestep: %02f s)"
+                    % (particles["t"] + particles["dt"])
+            )
+            log.error(message)
+            raise ValueError(message)
 
     particlesRelease = com1DFA.initializeParticles(
         cfg["GENERAL"],
         relLine,
         dem,
+        relThField=relThField,
+        timestep=particles["t"] + particles["dt"]
     )
-    particlesRelease = DFAfunC.updateInitialVelocity(cfg["GENERAL"], particlesRelease, dem, velocityMag)
+
+    if "velocity" in timeDepRelValues and "x" not in timeDepRelValues:
+        particlesRelease = DFAfunC.updateInitialVelocity(cfg["GENERAL"], particlesRelease, dem, velocityMag)
+
+
+
+    elif "velocityX" in timeDepRelValues and "x" in timeDepRelValues:
+        for uComp, timedepParameter in zip(["ux", "uy", "uz"], ["velocityX", "velocityY", "velocityZ"]):
+            raster = gI.timeDepRelCoordsToRaster(relLine["timeDepRelValues"], timeDepRelIndex,
+                                                 dem["originalHeader"], parameter=timedepParameter
+                                                 )
+            rasterDict = {"header": dem["header"], "rasterData": raster}
+            particlesRelease, _ = geoTrans.projectOnRaster(rasterDict,
+                                                           particlesRelease, outData=uComp)
+
+        particlesRelease["uMag"] = np.sqrt(
+            particlesRelease["ux"] ** 2 + particlesRelease["uy"] ** 2 + particlesRelease["uz"] ** 2)
+
     particles = particleTools.mergeParticleDict(particles, particlesRelease)
     # save initial z position for travel angle computation
     zPartArray0 = np.append(zPartArray0, copy.deepcopy(particlesRelease["z"]))
     return particles, zPartArray0
 
 
+def defineReleaseLineFromCoordinates(relFile, timeDepRelValues, demHeader):
+    """
+    define the release line and its thickness raster from coordinates read from time dependent release values
+
+    Parameters
+    ----------
+    relFile: pathlib.Path
+        directory to release file (csv file)
+    timeDepRelValues: dict
+        time dependent release values
+    demHeader: dict
+        header of DEM
+
+    Returns
+    -------
+    releaseLine: dict
+        dictionary for release line containing thickness raster data
+    relThFieldData: numpy array
+        release thickness raster data
+    """
+    releaseLine = {
+        "file": relFile,
+        "Name": [relFile.stem],
+        "initializedFrom": "csvfile"
+    }
+    initialIndex = np.where(timeDepRelValues["timeStep"] == 0)[0]
+    releaseLine["rasterData"] = gI.timeDepRelCoordsToRaster(timeDepRelValues, initialIndex, demHeader,
+                                                            parameter="thickness")
+    relThFieldData = releaseLine["rasterData"]
+    # TODO: define thickness for output report, now mean of thickness values
+    releaseLine["thickness"] = np.nanmean(np.where(relThFieldData == 0, np.nan, relThFieldData))
+
+    return releaseLine, relThFieldData
+
+
 def prepareTimeDepRelLine(releaseLine, cfg):
     """
     read time dependent release values and return them as a dictionary containing: