diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index 85a6cf81d..79004ba5a 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/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"])) diff --git a/avaframe/com1DFA/deriveParameterSet.py b/avaframe/com1DFA/deriveParameterSet.py index c80bcc8ca..92ff2f981 100644 --- a/avaframe/com1DFA/deriveParameterSet.py +++ b/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 diff --git a/avaframe/com1DFA/particleTools.py b/avaframe/com1DFA/particleTools.py index 6676407f0..0b3968dbf 100644 --- a/avaframe/com1DFA/particleTools.py +++ b/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 diff --git a/avaframe/in1Data/getInput.py b/avaframe/in1Data/getInput.py index acf4e7ac2..a27dcd7e0 100644 --- a/avaframe/in1Data/getInput.py +++ b/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 diff --git a/avaframe/in3Utils/cfgUtils.py b/avaframe/in3Utils/cfgUtils.py index 36d476e33..e947cc036 100644 --- a/avaframe/in3Utils/cfgUtils.py +++ b/avaframe/in3Utils/cfgUtils.py @@ -1036,6 +1036,7 @@ def readConfigurationInfoFromDone(avaDir, specDir="", latest=False): "nSave", "nIter", "simName", + "simTimestamp", ] ], how="left", diff --git a/avaframe/out3Plot/outParticlesAnalysis.py b/avaframe/out3Plot/outParticlesAnalysis.py index 256a63d97..4ac448d3b 100644 --- a/avaframe/out3Plot/outParticlesAnalysis.py +++ b/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))) diff --git a/avaframe/runScripts/runParticlesAssetsInfo.py b/avaframe/runScripts/runParticlesAssetsInfo.py new file mode 100644 index 000000000..292d2f868 --- /dev/null +++ b/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) diff --git a/avaframe/tests/test_particleTools.py b/avaframe/tests/test_particleTools.py index 3cb5820b7..d020d9998 100644 --- a/avaframe/tests/test_particleTools.py +++ b/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) diff --git a/docs/moduleOut3Plot.rst b/docs/moduleOut3Plot.rst index 7bd6e5ff3..cbe34e09f 100644 --- a/docs/moduleOut3Plot.rst +++ b/docs/moduleOut3Plot.rst @@ -245,3 +245,15 @@ To run +particle assets information +============================= +:py:mod:`out3Plot.particleAnalysisPlots` can also be used to create a plot that shows the cells affected by the particle +trajectories color-coded according to different assets classes (from high to low). This functionality is implemented only +for :py:mod:`com1DFA.com1DFA` and relies on particle dictionaries that are saved for each simulation run. For this, +adding *particles* to the ``resType`` and adjusting the desired saving time step in ``tSteps`` in your local copy +of ```com1DFACfg.ini`` is required. To reduce the amount of data that is saved, consider only exporting the required +particle properties by setting ``exportParticlePorperties`` to: *ID|indXDEM|indYDEM|x|y|inCellDEM|nPart*. In addition, +an assets raster file (with the same extent and cell size as the simulation DEM) in ``avalancheDir/Inputs/REFDATA`` with +the suffix *_ASSETS* is required. + +