[out3Plot]: add functionality for assets

avaframe/com1DFA/com1DFA.py

@@ -2073,6 +2073,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
@@ -2316,6 +2317,8 @@ def DFAIterate(cfg, particles, fields, dem, inputSimLines, outDir, cuSimName, si
         tCPUtimeLoop = time.time() - startTime
         tCPU["timeLoop"] = tCPU["timeLoop"] + tCPUtimeLoop
     tCPU["nIter"] = nIter
+
+    tCPU["simDate"] = 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

@@ -996,7 +996,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 (
@@ -1087,7 +1087,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
@@ -1098,8 +1098,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,51 @@ 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
+    """
+
+    particlesTimeArrays["assetsValue"] = (
+        np.zeros((particlesTimeArrays["ID"].shape[0], particlesTimeArrays["ID"].shape[1])) * np.nan
+    )
+    particleAssets = np.zeros((dem["header"]["nrows"], dem["header"]["ncols"])) * np.nan
+    # loop over each particle
+    for pId in range(len(particlesTimeArrays["ID"][0, :])):
+        # loop over each infrastructure class (sorted from low to high)
+        for l in uniqueAssets:
+            # check if particle meets infra at any time step to preselect particles
+            if len(np.intersect1d(particlesTimeArrays["inCellDEM"][:, pId], assetsValues["value_%d" % l])):
+                # for all time steps check if particle location is within infra
+                # if yes mark all prior particle locations with this class
+                # overridden by higher classes
+                for m in reversed(range(len(particlesTimeArrays["inCellDEM"][:, 0]))):
+                    if particlesTimeArrays["inCellDEM"][m, pId] in assetsValues["value_%d" % l]:
+                        particlesTimeArrays["assetsValue"][0:m, pId] = l
+                        # in infra raster mark cells for this particle with infra class
+                        indX = np.asarray(particlesTimeArrays["indXDEM"][0 : m + 1, pId], dtype=int)
+                        indY = np.asarray(particlesTimeArrays["indYDEM"][0 : m + 1, pId], dtype=int)
+                        particleAssets[indY, indX] = np.where(
+                            particleAssets[indY, indX] >= l, particleAssets[indY, indX], l
+                        )
+
+    return particleAssets, particlesTimeArrays

avaframe/in1Data/getInput.py

@@ -1178,3 +1178,73 @@ def getTimeDepRelCsv(timeDepRelCsv):
         "velocity": timeDepRelDF["velocity"].to_numpy(dtype=np.float64),
     }
     return timeDepRelValues, timeDepRelDF
+
+
+def preprocessAssets(avalancheDir, dem, sizeThreshold):
+    """fetch an assets raster with same extent as simulation DEM and create assets class info
+
+    Paratmeters
+    ------------
+    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]
+
+    # 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 infra 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,69 @@
+"""
+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)
+
+# load DEM
+dem = gI.readDEM(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)

avaframe/tests/test_particleTools.py

@@ -293,7 +293,6 @@ def test_reshapeParticlesDicts():
     assert np.array_equal(particlesTimeArrays['velMag'], test)
     assert np.array_equal(particlesTimeArrays['uAcc'], np.asarray([[4., 5., 6., 7.], [4., 5., 7., 7.], [40., 50., 60., 70.]]))
 
-
     # setup required input
     partDict1 = {'velMag': np.asarray([1.,2., 4., 5.]),'uX': np.asarray([10.,20., 40., 50.]),
         'uAcc': np.asarray([4., 5., 6., 7.]), 'ID': np.asarray([0, 1, 2, 3]), 't': 0., 'nPart': 4}
@@ -307,3 +306,92 @@ def test_reshapeParticlesDicts():
         # call function to be tested
         particlesTimeArrays = particleTools.reshapeParticlesDicts(particlesList, ['velMag', 'uAcc', 't', 'ID'])
     assert str(e.value) == ("Number of particles changed throughout simulation")
+
+
+def test_createAssetsRasterFromParticleLocations():
+    """test creating a raster indicating particle trajectories colorcoded with infra classes"""
+
+    # setup required input data
+    dem = {
+        "header": {"nrows": 8, "ncols": 10, "cellsize": 1, "xllcenter": 0, "yllcenter": 0},
+        "rasterData": np.ones((8, 10)),
+    }
+    uniqueAssets = np.asarray([2, 3])
+    assets = {"header": {"nrows": 8, "ncols": 10, "cellsize": 1, "xllcenter": 0, "yllcenter": 0}}
+    assetsRaster = np.zeros((8, 10)) * np.nan
+    assetsRaster[1, 6] = 2.0
+    assetsRaster[3, 4] = 3.0
+    assetsRaster[5, 6] = 3.0
+    assetsRaster[6, 5] = 2.0
+    assets["rasterData"] = assetsRaster
+    # 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 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]
+
+    particlesTimeArrays = {
+        "ID": np.asarray([[1, 2], [1, 2], [1, 2], [1, 2], [1, 2], [1, 2], [1, 2]]),
+        "indXDEM": np.asarray(
+            [
+                [0, 4],
+                [1, 5],
+                [2, 6],
+                [3, 7],
+                [4, 8],
+                [5, 9],
+                [6, 9],
+            ]
+        ),
+        "indYDEM": np.asarray(
+            [
+                [7, 7],
+                [6, 6],
+                [5, 5],
+                [4, 4],
+                [3, 3],
+                [2, 2],
+                [1, 2],
+            ]
+        ),
+        "inCellDEM": np.asarray(
+            [
+                [cellNo[7, 0], cellNo[7, 4]],
+                [cellNo[6, 1], cellNo[6, 5]],
+                [cellNo[5, 2], cellNo[5, 6]],
+                [cellNo[4, 3], cellNo[4, 7]],
+                [cellNo[3, 4], cellNo[3, 8]],
+                [cellNo[2, 5], cellNo[2, 9]],
+                [cellNo[1, 6], cellNo[2, 9]],
+            ]
+        ),
+    }
+
+    # call function to be tested
+    particleAssets, particleTimeArrays = particleTools.createAssetsRasterFromParticleLocations(
+        particlesTimeArrays, dem, uniqueAssets, assetsValues
+    )
+
+    particleAssetsTest = np.zeros((8, 10)) * np.nan
+    particleAssetsTest[7, 0] = 3.0
+    particleAssetsTest[6, 1] = 3.0
+    particleAssetsTest[5, 2] = 3.0
+    particleAssetsTest[4, 3] = 3.0
+    particleAssetsTest[3, 4] = 3.0
+    particleAssetsTest[2, 5] = 2.0
+    particleAssetsTest[1, 6] = 2.0
+
+    particleAssetsTest[7, 4] = 3.0
+    particleAssetsTest[6, 5] = 3.0
+    particleAssetsTest[5, 6] = 3.0
+
+    print(particleAssets)
+    print(particleAssetsTest)
+    assert particleAssets.shape == (8, 10)
+    assert np.array_equal(particleAssets, particleAssetsTest, equal_nan=True)