[out3Plot]: add functionality for assets

.github/workflows/runTestSinglePython.yml

@@ -23,7 +23,7 @@ jobs:
       - uses: actions/checkout@v5
         with:
           fetch-depth: 0
-      - uses: prefix-dev/setup-pixi@v0.9.1
+      - uses: prefix-dev/setup-pixi@v0.9.6
         with:
           cache: false
       - name: Pixi run pytest

avaframe/com1DFA/com1DFA.py

@@ -2075,6 +2075,7 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
         "timePos": 0.0,
         "timeNeigh": 0.0,
         "timeField": 0.0,
+        "simTimestamp": 0.0,
     }
 
     # Load configuration settings
@@ -2318,6 +2319,8 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
         tCPUtimeLoop = time.time() - startTime
         tCPU["timeLoop"] = tCPU["timeLoop"] + tCPUtimeLoop
     tCPU["nIter"] = nIter
+
+    tCPU["simTimestamp"] = datetime.now().strftime("%Y%m%d_%Hh%Mm%Ss")
     log.info("Ending computation at time t = %f s", t - dt)
     log.debug("Saving results for time step t = %f s", t - dt)
     log.info("MTot = %f kg, %s particles" % (particles["mTot"], particles["nPart"]))

avaframe/com1DFA/deriveParameterSet.py

@@ -1003,7 +1003,7 @@ def checkExtentAndCellSize(cfg, inputFile, dem, fileType):
     rT = float(cfg["GENERAL"]["resizeThreshold"])
     cT = float(cfg["GENERAL"]["meshCellSizeThreshold"])
 
-    diffX0, diffX1, diffY0, diffY1 = checkSizeExtent(inputField, demHeader, inputFile, fileType, rT)
+    diffX0, diffX1, diffY0, diffY1 = checkSizeExtent(inputField, demHeader, fileType, rT)
 
     # check if identical extent, if so use unchanged
     if (
@@ -1094,7 +1094,7 @@ def checkExtentAndCellSize(cfg, inputFile, dem, fileType):
     return returnStr, outFile, remeshedFlag
 
 
-def checkSizeExtent(inputField, demHeader, inputFile, fileType, rT):
+def checkSizeExtent(inputField, demHeader, fileType, rT):
     """check if extent of an inputfield matches the extent of the DEM
     and also cellSize in case of RELTH files
     optionally within a specified threshold
@@ -1105,8 +1105,6 @@ def checkSizeExtent(inputField, demHeader, inputFile, fileType, rT):
         dictionary with header and data of input field
     demHeader: dict
         header of DEM
-    inputFile: pathlib Path
-        path to input field
     fileType: str
         name of file type of input field
     rT: float

avaframe/com1DFA/particleTools.py

@@ -1009,3 +1009,59 @@ def savePartDictToPickle(partDict, fName):
     fi = open(fName, "wb")
     pickle.dump(partDict, fi)
     fi.close()
+
+
+def createAssetsRasterFromParticleLocations(particlesTimeArrays, dem, uniqueAssets, assetsValues):
+    """create a raster indicating particle trajectories colorcoded with assets classes, highest overrides lower classes
+
+     Parameters
+    -----------
+    particlesTimeArrays: dict
+        dictionary with time series of properties of particles
+    dem: dict
+        dictionary with dem information header nrows, ncols required
+    uniqueAssets: list
+        list of assets class values sorted from low to high
+    assetsValues: dict
+        dictionary with for each infrastructure class value the affected cell numbers
+
+    Returns
+    ---------
+    particleAssets: numpy ndarray
+        array with particle trajectories colorcoded with assets classes
+    particlesTimeArrays: dict
+        updated with assets value
+
+    """
+
+    # initialize particle assets arrays with nans
+    nTime, nPart = particlesTimeArrays["ID"].shape
+    particlesTimeArrays["assetsValue"] = np.full((nTime, nPart), np.nan)
+    particleAssets = np.full((dem["header"]["nrows"], dem["header"]["ncols"]), np.nan)
+
+    # process classes from low to high so higher classes naturally override lower ones
+    for assetClass in uniqueAssets:
+        # find for each particle if its trajectory has an overlap with an asset class
+        assetCells = np.asarray(assetsValues["value_%d" % assetClass])
+        inAsset = np.isin(particlesTimeArrays["inCellDEM"], assetCells)
+
+        # loop over all particles
+        for pId in range(nPart):
+            # if particle trajectory has overlap mark all time steps until the last time step it has overlap
+            # if not leave loop
+            hitTimes = np.where(inAsset[:, pId])[0]
+            if len(hitTimes) == 0:
+                continue
+
+            # mark particle trajectory up to the last time it has overlap with this asset class
+            mMax = hitTimes[-1]
+            particlesTimeArrays["assetsValue"][:mMax, pId] = assetClass
+            # find indices of respective cells and mark them in particleAssets array
+            indX = particlesTimeArrays["indXDEM"][: mMax + 1, pId].astype(int)
+            indY = particlesTimeArrays["indYDEM"][: mMax + 1, pId].astype(int)
+            # only overwrite cells not already set to a higher class
+            particleAssets[indY, indX] = np.where(
+                particleAssets[indY, indX] >= assetClass, particleAssets[indY, indX], assetClass
+            )
+
+    return particleAssets, particlesTimeArrays

avaframe/in1Data/getInput.py

@@ -1266,3 +1266,74 @@ def checkTimeDepRelease(timeDepRelValues, timeDepRelCsv):
             message = "The initial velocity provided in %s can not be negative." % (timeDepRelCsv)
             log.error(message)
             raise ValueError(message)
+
+
+def preprocessAssets(avalancheDir, dem, sizeThreshold):
+    """fetch an assets raster with same extent as simulation DEM and create assets class info
+       nan values and 0 are disregarded as asset class
+
+    Parameters
+    ------------
+    avalancheDir: str or pathlib.Path
+        path to avalanche directory
+    dem: dict
+        dictionary with info on DEM
+    sizeThreshold: float
+        threshold for determining if asset layer location of origin and extent are in line with
+        computational mesh of simulation
+
+    Returns
+    ---------
+    uniqueAssets: list
+        list of asset class values ordered from low to high
+    infra: dict
+        dictionary with assets info
+    infraValues: dict
+        dictionary with affected cell numbers for each asset class value
+
+
+    """
+
+    # load infrastructure data
+    infraPath = pathlib.Path(avalancheDir, "Inputs", "REFDATA")
+    assetsSuffix = "*_ASSETS.asc"
+    assetsFileList = list(infraPath.glob(assetsSuffix))
+    if len(assetsFileList) > 1:
+        message = "More than one assets class file (%s) found in %s, this is not allowed" % (
+            assetsSuffix,
+            str(infraPath),
+        )
+        log.error(message)
+        raise ValueError(message)
+    elif len(assetsFileList) == 0:
+        message = "No assets class file (%s) found in %s" % (assetsSuffix)
+        log.error(message)
+        raise FileNotFoundError(message)
+    else:
+        assetsFile = list(infraPath.glob(assetsSuffix))[0]
+    assets = IOf.readRaster(assetsFile, noDataToNan=True)
+
+    # check if location of origin, extent and cellsize are in line with computational mesh of simulation
+    # NOTE: no remeshing if cellsize doesn't match so that no smearing of classes or boundaries of assets due to interpolation
+    _, _, _, _ = dP.checkSizeExtent(assets, dem["header"], "assets", sizeThreshold)
+    if np.abs(dem["header"]["cellsize"] - assets["header"]["cellsize"]) > sizeThreshold:
+        message = "Assets raster (%s) cell size does not match simulation dem cell size" % assetsFile.name
+        log.error(message)
+        raise AssertionError(message)
+
+    # create cell number
+    cellNo = np.zeros((dem["header"]["nrows"], dem["header"]["ncols"]))
+    for m in range(dem["header"]["nrows"]):
+        for k in range(dem["header"]["ncols"]):
+            cellNo[m, k] = k + m * dem["header"]["ncols"]
+    # fetch available assets classes and sort low to high class
+    uniqueAssets = np.sort(np.unique(assets["rasterData"]))
+    uniqueAssets = [i for i in uniqueAssets if i != 0 and not np.isnan(i)]
+
+    # fetch corresponding cell numbers for all infrastructure classes
+    assetsValues = {}
+    for i in uniqueAssets:
+        assetsArray = np.where(assets["rasterData"] == i, cellNo, np.nan)
+        assetsValues["value_%d" % i] = [int(k) for k in assetsArray.flatten() if np.isnan(k) == False]
+
+    return uniqueAssets, assets, assetsValues

avaframe/in3Utils/cfgUtils.py

@@ -1036,6 +1036,7 @@ def readConfigurationInfoFromDone(avaDir, specDir="", latest=False):
                             "nSave",
                             "nIter",
                             "simName",
+                            "simTimestamp",
                         ]
                     ],
                     how="left",

avaframe/out3Plot/outParticlesAnalysis.py

@@ -829,3 +829,48 @@ def readMeasuredParticleData(avalancheDir, demHeader, pData=""):
     mParticles["y"] = mParticles["y"] - demHeader["yllcenter"]
 
     return mParticles
+
+
+def plotParticlesInfra(dem, infra, particleInfra, outDir, plotName):
+    """Plot locations of particles over all timesteps that reached infrastructure color coded with highes infra class
+    Parameters
+    -----------
+    dem: dict
+        dictionary with dem info
+    infra: dict
+        dictionary with infrastructure info
+    particleInfra: dict
+        dictionary with color coded particles trajectories
+    outDir: pathlib.Path
+        path to output directory
+    plotName: str
+        name of plot file
+    """
+    extentCellCenters, extentCellCorners = pU.createExtentMinMax(
+        dem["rasterData"], dem["header"], originLLCenter=True
+    )
+
+    fig, ax = plt.subplots(ncols=1)
+    ax.set_title("Particle trajectories color-coded with asset classes")
+    # add DEM hillshade with contour lines
+    # set extent in meters using cellSize and llcenter location
+    _, _ = pU.addHillShadeContours(ax, dem["rasterData"], dem["header"]["cellsize"], extentCellCenters)
+    ax.imshow(
+        np.where(infra["rasterData"] != 0, infra["rasterData"], np.nan),
+        alpha=1.0,
+        extent=extentCellCenters,
+        origin="lower",
+        cmap=cm.hawaii,
+        zorder=11,
+    )
+    im1 = ax.imshow(
+        particleInfra, extent=extentCellCenters, origin="lower", alpha=0.6, cmap=cm.hawaii, zorder=10
+    )
+
+    ax.set_xlabel("x [m]")
+    ax.set_ylabel("y [m]")
+    fig.colorbar(im1, ax=ax)
+
+    # save and or plot
+    plotPath = pU.saveAndOrPlot({"pathResult": outDir}, plotName, fig)
+    log.info("Plot for %s successfully saved at %s" % (plotName, str(plotPath)))

avaframe/runScripts/runParticlesAssetsInfo.py

@@ -0,0 +1,66 @@
+"""
+Run creating a raster with numerical particle trajectories colorcoded with assets classes, highest overrides lower classes
+"""
+
+import pathlib
+
+# Local imports
+import avaframe.in1Data.getInput as gI
+from avaframe.in3Utils import cfgUtils
+import avaframe.out3Plot.outParticlesAnalysis as oP
+import avaframe.com1DFA.particleTools as pT
+import avaframe.in2Trans.rasterUtils as rU
+import avaframe.in3Utils.fileHandlerUtils as fU
+from avaframe.in3Utils import logUtils
+
+# load avalanche directory
+cfgMain = cfgUtils.getGeneralConfig()
+avalancheDir = cfgMain["MAIN"]["avalancheDir"]
+outDir = pathlib.Path(avalancheDir, "Outputs", "out3Plot", "particleAnalysis")
+fU.makeADir(outDir)
+
+# log file name; leave empty to use default runLog.log
+logName = "runParticlesAssetsInfo"
+# Start logging
+log = logUtils.initiateLogger(avalancheDir, logName)
+log.info("MAIN SCRIPT")
+log.info("Current avalanche: %s", avalancheDir)
+
+# create data frame that lists all available simulations
+inputsDF, resTypeList = fU.makeSimFromResDF(avalancheDir, "com1DFA")
+# load dataFrame for all configurations
+configurationDF = cfgUtils.createConfigurationInfo(avalancheDir, comModule="com1DFA")
+# Merge inputsDF with the configurationDF. Make sure to keep the indexing from inputs and to merge on 'simName'
+inputsDF = inputsDF.reset_index().merge(configurationDF, on=["simName", "modelType"]).set_index("index")
+
+# loop over all sims found
+for index, row in inputsDF.iterrows():
+    # fetch simName
+    simName = row["simName"]
+    dem = rU.readRaster(pathlib.Path(avalancheDir, "Inputs", row["DEM"]))
+    log.info("Find particle trajectories for simulation: %s" % (simName))
+
+    # fetch info on infrastructure
+    uniqueAssets, assets, assetsValues = gI.preprocessAssets(avalancheDir, dem, sizeThreshold=1.0e-10)
+
+    # read particles saved from com1DFA simulation
+    Particles, timeStepInfo = pT.readPartFromPickle(
+        pathlib.Path(avalancheDir), simName=simName, flagAvaDir=True, comModule="com1DFA"
+    )
+    # create time series of particles arrays
+    particlesTimeArrays = pT.reshapeParticlesDicts(
+        Particles, ["ID", "indXDEM", "indYDEM", "x", "y", "inCellDEM"]
+    )
+
+    # derive info on which particles interacted with infrastructure
+    particleAssets, particlesTimeArrays = pT.createAssetsRasterFromParticleLocations(
+        particlesTimeArrays, dem, uniqueAssets, assetsValues
+    )
+    # export raster
+    rU.writeResultToRaster(
+        dem["header"], particleAssets, (outDir / ("particleAssetsInfo_%s" % simName)), flip=True
+    )
+
+    # create plot
+    plotName = "particleAssetsInfo_%s" % simName
+    _ = oP.plotParticlesInfra(dem, assets, particleAssets, outDir, plotName)