diff --git a/.github/workflows/require-checklist.yml b/.github/workflows/require-checklist.yml
index 32e86c948..2ece0fc48 100644
--- a/.github/workflows/require-checklist.yml
+++ b/.github/workflows/require-checklist.yml
@@ -3,6 +3,8 @@ name: Require PR Checklist
 on:
   pull_request:
     types: [opened, edited, synchronize, reopened]
+  workflow_dispatch:
+
 
 jobs:
   require-checklist:
diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py
index e539aa29a..f02b5a9bf 100644
--- a/avaframe/com1DFA/com1DFA.py
+++ b/avaframe/com1DFA/com1DFA.py
@@ -2123,6 +2123,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
@@ -2366,6 +2367,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"]))
diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini
index 884f826a9..4b608f955 100644
--- a/avaframe/com1DFA/com1DFACfg.ini
+++ b/avaframe/com1DFA/com1DFACfg.ini
@@ -246,6 +246,8 @@ meshCellSizeThreshold = 0.001
 cleanRemeshedRasters = True
 # remesh if extent matches computational DEM within resizeThreshold x meshCellSize
 resizeThreshold = 3
+# remeshing option for raster input data
+remeshInterpMethod = default
 
 # Normal computation on rectangular grid
 # 4 triangles method        6 triangles method         8 triangles method
diff --git a/avaframe/com1DFA/deriveParameterSet.py b/avaframe/com1DFA/deriveParameterSet.py
index 2e0b2e512..cfa199629 100644
--- a/avaframe/com1DFA/deriveParameterSet.py
+++ b/avaframe/com1DFA/deriveParameterSet.py
@@ -965,7 +965,7 @@ def checkRasterMeshSize(cfgSim, rasterFile, typeIndicator="DEM", onlySearch=Fals
     return pathToRaster
 
 
-def checkExtentAndCellSize(cfg, inputFile, dem, fileType):
+def checkExtentAndCellSize(cfg, inputFile, dem, fileType, nanInsideDEMCheck=True):
     """check if extent of inputFile is within resizeThreshold of dem, if so resize and save to remeshedRasters
 
     Parameters
@@ -978,12 +978,14 @@ def checkExtentAndCellSize(cfg, inputFile, dem, fileType):
         dictionary with info on DEM
     fileType: str
         name of fileType
+    nanInsideDEMCheck: bool
+        if True check if in remeshed file only nans where also DEM has nans
     """
 
     inputField = IOf.readRaster(inputFile)
 
     # Check for NaN values before remeshing - not allowed in non-DEM input rasters
-    if fileType.upper() != "DEM":
+    if fileType.upper() != "DEM" and fileType.upper() != "ASSETS":
         if np.isnan(inputField["rasterData"]).any():
             message = "In %s file (%s) nan values found - this is not allowed" % (
                 fileType,
@@ -996,18 +998,12 @@ def checkExtentAndCellSize(cfg, inputFile, dem, fileType):
     demHeader = dem["header"]
 
     # check if negative values in raster file
-    if np.any(inputField["rasterData"] < 0):
-        message = "In %s file (%s) negative values found - this is not allowed" % (
-            fileType,
-            inputFile.name,
-        )
-        log.error(message)
-        raise AssertionError(message)
+    checkNegativeInRaster(inputField, fileType, inputFile)
 
     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 (
@@ -1020,7 +1016,9 @@ def checkExtentAndCellSize(cfg, inputFile, dem, fileType):
         remeshedFlag = "No"
     else:
         # resize data, project data from inputFile onto computational domain
-        inputField["rasterData"], _ = geoTrans.resizeData(inputField, dem)
+        inputField["rasterData"], _ = geoTrans.resizeData(
+            inputField, dem, interp=cfg["GENERAL"]["remeshInterpMethod"]
+        )
 
         # add warning
         log.warning(
@@ -1065,40 +1063,65 @@ def checkExtentAndCellSize(cfg, inputFile, dem, fileType):
         returnStr = str(pathlib.Path("remeshedRasters", outFile.name))
         remeshedFlag = "Yes"
 
-    # check if no data values only where DEM also has no data values
-    # if remeshed - potentially nans at edges, if inputField has a different origin/extent
-    # for example if extent of DEM is larger -> nans in this region in remeshed input file
-    # first maks nans values of DEM as there, also inputField is allowed to have nans
-    nanDEMMasked = np.where(np.isnan(dem["rasterData"]), -9999, inputField["rasterData"])
-    # search for indices where nans come from remeshing use difference in extent prior to remeshing
-    # if diff is negative on Left side DEM is smaller
-    # if diff is negative on right side DEM is larger
-    nNeglectRowsMin = int(np.ceil(abs(min(0, diffY0 / dem["header"]["cellsize"]))))
-    nNeglectRowsMax = int(np.ceil(abs(max(0, diffY1 / dem["header"]["cellsize"]))))
-    # if diff is negative on lower side DEM is smaller
-    # if diff is negative on right side DEM is larger
-    nNeglectColsMin = int(np.ceil(abs(min(0, diffX0 / dem["header"]["cellsize"]))))
-    nNeglectColsMax = int(np.ceil(abs(max(0, diffX1 / dem["header"]["cellsize"]))))
-    maxRows = dem["header"]["nrows"] - nNeglectRowsMin
-    maxCols = dem["header"]["ncols"] - nNeglectColsMin
-
-    # add mask where nans come from remeshing
-    # change order because diff is with origin lower!!
-    nanDEMMaskedLimited = nanDEMMasked[nNeglectRowsMax:maxRows, nNeglectColsMax:maxCols]
-
-    # check if nan values in raster file data (excluding nans from remeshing and where also DEM has nans)
-    if np.any(np.isnan(nanDEMMaskedLimited)):
-        message = "In %s file (%s) nan values found inside DEM extent - this is not allowed" % (
+    if nanInsideDEMCheck:
+        # check if no data values only where DEM also has no data values
+        # if remeshed - potentially nans at edges, if inputField has a different origin/extent
+        # for example if extent of DEM is larger -> nans in this region in remeshed input file
+        # first maks nans values of DEM as there, also inputField is allowed to have nans
+        nanDEMMasked = np.where(np.isnan(dem["rasterData"]), -9999, inputField["rasterData"])
+        # search for indices where nans come from remeshing use difference in extent prior to remeshing
+        # if diff is negative on Left side DEM is smaller
+        # if diff is negative on right side DEM is larger
+        nNeglectRowsMin = int(np.ceil(abs(min(0, diffY0 / dem["header"]["cellsize"]))))
+        nNeglectRowsMax = int(np.ceil(abs(max(0, diffY1 / dem["header"]["cellsize"]))))
+        # if diff is negative on lower side DEM is smaller
+        # if diff is negative on right side DEM is larger
+        nNeglectColsMin = int(np.ceil(abs(min(0, diffX0 / dem["header"]["cellsize"]))))
+        nNeglectColsMax = int(np.ceil(abs(max(0, diffX1 / dem["header"]["cellsize"]))))
+        maxRows = dem["header"]["nrows"] - nNeglectRowsMin
+        maxCols = dem["header"]["ncols"] - nNeglectColsMin
+
+        # add mask where nans come from remeshing
+        # change order because diff is with origin lower!!
+        nanDEMMaskedLimited = nanDEMMasked[nNeglectRowsMax:maxRows, nNeglectColsMax:maxCols]
+
+        # check if nan values in raster file data (excluding nans from remeshing and where also DEM has nans)
+        if np.any(np.isnan(nanDEMMaskedLimited)):
+            message = "In %s file (%s) nan values found inside DEM extent - this is not allowed" % (
+                fileType,
+                inputFile.name,
+            )
+            log.error(message)
+            raise AssertionError(message)
+
+    return returnStr, outFile, remeshedFlag
+
+
+def checkNegativeInRaster(inputField, fileType, inputFile):
+    """check if inputField['rasterData'] contains negative values, if True error
+
+    Parameters
+    -------------
+    inputField: dict
+        dictionary with header and rasterData to be checked
+    fileType: str
+        file type used in error message
+    inputFile: pathlib.Path
+        path to input file - used in error message
+
+    """
+
+    # check if negative values in raster file
+    if np.any(inputField["rasterData"] < 0):
+        message = "In %s file (%s) negative values found - this is not allowed" % (
             fileType,
             inputFile.name,
         )
         log.error(message)
         raise AssertionError(message)
 
-    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
@@ -1109,8 +1132,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..3a4e05df3 100644
--- a/avaframe/com1DFA/particleTools.py
+++ b/avaframe/com1DFA/particleTools.py
@@ -16,7 +16,7 @@
 import avaframe.com1DFA.DFAtools as DFAtls
 import avaframe.in3Utils.geoTrans as geoTrans
 import avaframe.com1DFA.DFAfunctionsCython as DFAfunC
-
+import avaframe.out3Plot.outDebugPlots as outDebug
 
 # create local logger
 # change log level in calling module to DEBUG to see log messages
@@ -1009,3 +1009,178 @@ 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
+            )
+
+    # find indices of all cells that were affected by particles
+    xyIndAllUnique = findUniqueCellIndices(particlesTimeArrays)
+
+    # set all locations where particles were but not affecting assets to class -1.0
+    testArray = np.full((dem["header"]["nrows"], dem["header"]["ncols"]), np.nan)
+    testArray[xyIndAllUnique[:, 1], xyIndAllUnique[:, 0]] = -1.0
+    particleAssets = np.where(~np.isnan(particleAssets), particleAssets, testArray)
+
+    return particleAssets, particlesTimeArrays
+
+
+def findUniqueCellIndices(particlesTimeArrays):
+    """find indices of all cells that are saved for each particle as the particle location was within that cell for the
+    respective time step
+
+    Parameters
+    --------------
+    particlesTimeArrays: dict
+        dictionary with time series of properties of particles, keys: property each with a timeSteps x number of
+        particles array
+
+    Returns
+    ----------
+    xyIndAllUnique: np.ndarray
+        array with 2 columns: column 1: X indices of all affected cells column 2: Y indices of all affected cells
+
+    """
+
+    xyIndAll = np.column_stack(
+        (particlesTimeArrays["indXDEM"].flatten(), particlesTimeArrays["indYDEM"].flatten())
+    )
+    xyIndAllUnique = np.unique(xyIndAll, axis=0)
+    xyIndAllUnique = xyIndAllUnique[~np.isnan(xyIndAllUnique).any(axis=1)]
+    xyIndAllUnique = np.asarray(xyIndAllUnique, dtype=int)
+
+    return xyIndAllUnique
+
+
+def interpolateParticlesTrajectories(dem, particlesTimeArrays, cellSizeFactor, debugPlot=False):
+    """interpolate particle trajectories to have more closely spaced values important if affected grid cells
+    should be identified without gaps
+
+       Parameters
+       --------------
+       dem: dict
+           dictionary with info on dem header and rasterData
+       particlesTimeArrays: dict
+           dictionary with time series of properties of particles
+        cellSizeFactor: float
+            dem mesh cellsize x cellSizeFactor will give desired distance of interpolated points
+
+       Returns
+       ----------
+       pLong: dict
+            updated dictionary with time series of particles location, affected cells
+
+    """
+
+    # fetch initial information
+    ncols = dem["header"]["ncols"]
+    cellSize = dem["header"]["cellsize"]
+    x = particlesTimeArrays["x"]
+    y = particlesTimeArrays["y"]
+    z = particlesTimeArrays["z"]
+    nTime, nPart = particlesTimeArrays["ID"].shape
+    # computational mesh has origin 0,0 in com1DFA required for projectOnRaster for z coordinates in prepareLine function
+    demComputation = {
+        "header": {"xllcenter": 0, "yllcenter": 0, "cellsize": dem["header"]["cellsize"]},
+        "rasterData": dem["rasterData"],
+    }
+
+    # initialize dict for interpolated particle trajectories
+    pLong = {}
+    indLongest = 0
+    for k in range(nPart):
+        pTraj = {"x": x[:, k], "y": y[:, k], "z": z[:, k]}
+        # interpolate particle trajectories
+        pTraj, _ = geoTrans.prepareLineStrict(
+            demComputation,
+            pTraj,
+            np.floor(dem["header"]["cellsize"] * cellSizeFactor),
+            Point=None,
+        )
+        # plot particle trajectories original and interpolated result on grid
+        if debugPlot:
+            outDebug.plotParticleTrajOnGrid(x[:, k], y[:, k], pTraj["x"], pTraj["y"], dem)
+
+        pLong["%s_pTraj" % particlesTimeArrays["ID"][0, k]] = pTraj
+        if (pTraj["x"].shape[0]) > indLongest:
+            indLongest = pTraj["x"].shape[0]
+
+    # initialize arrays
+    pLong["indXDEM"] = np.zeros((indLongest, nPart))
+    pLong["indYDEM"] = np.zeros((indLongest, nPart))
+    pLong["inCellDEM"] = np.zeros((indLongest, nPart))
+    pLong["ID"] = np.zeros((indLongest, nPart))
+    pLong["x"] = np.zeros((indLongest, nPart))
+    pLong["y"] = np.zeros((indLongest, nPart))
+    # pLong["z"] = np.zeros((indLongest, nPart))
+
+    # find grid cells for each particle location
+    for k in range(nPart):
+        x1 = np.full(indLongest, np.nan)
+        lenX = pLong["%s_pTraj" % particlesTimeArrays["ID"][0, k]]["x"].shape[0]
+        x1[0:lenX] = pLong["%s_pTraj" % particlesTimeArrays["ID"][0, k]]["x"]
+        y1 = np.full(indLongest, np.nan)
+        y1[0:lenX] = pLong["%s_pTraj" % particlesTimeArrays["ID"][0, k]]["y"]
+        # z1 = np.full(indLongest, np.nan)
+        # z1[0:lenX] = pLong["%s_pTraj" % particlesTimeArrays["ID"][0, k]]["z"]
+        # find cell indices
+        pLong["indXDEM"][:, k] = np.round(x1[:] / cellSize)
+        pLong["indYDEM"][:, k] = np.round(y1[:] / cellSize)
+        # get index of cell containing the particle
+        pLong["inCellDEM"][:, k] = pLong["indXDEM"][:, k] + ncols * pLong["indYDEM"][:, k]
+        pLong["ID"][:, k] = particlesTimeArrays["ID"][0, k]
+        pLong["x"][:, k] = x1
+        pLong["y"][:, k] = y1
+        # pLong["z"][:, k] = z1
+
+    return pLong
diff --git a/avaframe/in1Data/getInput.py b/avaframe/in1Data/getInput.py
index 4f9c451db..31e6f6f62 100644
--- a/avaframe/in1Data/getInput.py
+++ b/avaframe/in1Data/getInput.py
@@ -1266,3 +1266,66 @@ 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, cfg):
+    """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
+    cfg: configparser object
+        configuration settings for remeshing
+
+    Returns
+    ---------
+    uniqueAssets: list
+        list of asset class values ordered from low to high
+    assets: dict
+        dictionary with assets info
+    assetsValues: dict
+        dictionary with affected cell numbers for each asset class value
+
+
+    """
+
+    # load infrastructure data
+    assetsPath = pathlib.Path(avalancheDir, "Inputs", "INFRA")
+    assetsFileList = list(assetsPath.glob("*.asc")) + list(assetsPath.glob("*.tif"))
+    if len(assetsFileList) > 1:
+        message = "More than one assets class file found in %s, this is not allowed" % (str(assetsPath))
+        log.error(message)
+        raise ValueError(message)
+    elif len(assetsFileList) == 0:
+        message = "No assets class file found in %s" % (assetsPath)
+        log.error(message)
+        raise FileNotFoundError(message)
+    else:
+        assetsFile = assetsFileList[0]
+
+    # check extent and cell size of assets raster - if not aligned with computational mesh - remesh
+    pathToAssets, pathToAssetsFull, remeshedAssets = dP.checkExtentAndCellSize(
+        cfg, assetsFile, dem, "ASSETS", nanInsideDEMCheck=False
+    )
+    assets = IOf.readRaster(pathToAssetsFull, noDataToNan=True)
+
+    # 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
diff --git a/avaframe/in3Utils/cfgUtils.py b/avaframe/in3Utils/cfgUtils.py
index 2c341b02f..ed791b129 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/in3Utils/geoTrans.py b/avaframe/in3Utils/geoTrans.py
index 12e6984b8..f3c0428fe 100644
--- a/avaframe/in3Utils/geoTrans.py
+++ b/avaframe/in3Utils/geoTrans.py
@@ -168,7 +168,7 @@ def projectOnGrid(x, y, Z, csz=1, xllc=0, yllc=0, interp="bilinear", getXYField=
         return z, ioob
 
 
-def resizeData(raster, rasterRef):
+def resizeData(raster, rasterRef, interp="default"):
     """
     Reproject raster on a grid of shape rasterRef, raster adapts cellsize and extend of rasterRef
 
@@ -178,6 +178,8 @@ def resizeData(raster, rasterRef):
         raster dictionary
     rasterRef : dict
         reference raster dictionary
+    interp: str
+        default - refers to bilinear, other option is nearest
 
     Returns
     -------
@@ -189,6 +191,18 @@ def resizeData(raster, rasterRef):
     if rU.isEqualASCheader(raster["header"], rasterRef["header"]):
         return raster["rasterData"], rasterRef["rasterData"]
     else:
+        if interp == "default":
+            interpMethod = "bilinear"
+        elif interp == "nearest":
+            interpMethod = "nearest"
+        else:
+            message = (
+                'Interpolation method "%s" not recognized (valid options: "default" (bilinear), "nearest" (nearest-neighbor))'
+                % (interp)
+            )
+            log.error(message)
+            raise NameError(message)
+
         headerRef = rasterRef["header"]
         ncols = headerRef["ncols"]
         nrows = headerRef["nrows"]
@@ -199,7 +213,7 @@ def resizeData(raster, rasterRef):
         ygrid = np.linspace(yllc, yllc + (nrows - 1) * csz, nrows)
         X, Y = np.meshgrid(xgrid, ygrid)
         Points = {"x": X, "y": Y}
-        Points, _ = projectOnRaster(raster, Points, interp="bilinear")
+        Points, _ = projectOnRaster(raster, Points, interp=interpMethod)
         bilinearData = Points["z"]
 
         if np.isnan(bilinearData).any():
@@ -231,7 +245,8 @@ def remeshData(rasterDict, cellSizeNew, remeshOption="griddata", interpMethod="c
         Check the scipy documentation for more details
         default is 'griddata'
     interpMethod: str
-        interpolation order to use for the interpolation ('linear', 'cubic' or 'quintic')
+        interpolation order to use for the interpolation ('default', 'linear', 'cubic' or 'quintic')
+        default refers to cubic
     larger: Boolean
         if true (default) output grid is at least as big as the input
 
@@ -241,8 +256,21 @@ def remeshData(rasterDict, cellSizeNew, remeshOption="griddata", interpMethod="c
         remeshed data dict with data as numpy array and header info
 
     """
+
     header = rasterDict["header"]
 
+    # TODO: shall we allow other options?
+    if interpMethod == "default":
+        interpMethod = "cubic"
+    elif interpMethod in ["nearest", "linear", "cubic"] and remeshOption == "griddata":
+        interpMethod = interpMethod
+    elif interpMethod in ["linear", "quintic", "cubic"] and remeshOption == "RectBivariateSpline":
+        interpMethod = interpMethod
+    else:
+        message = 'Interpolation method "%s" not recognized' % (interpMethod)
+        log.error(message)
+        raise NameError(message)
+
     # fetch shape info and get new mesh info
     xGrid, yGrid, _, _ = makeCoordGridFromHeader(header)
     xGridNew, yGridNew, ncolsNew, nrowsNew = makeCoordGridFromHeader(
@@ -306,7 +334,7 @@ def remeshData(rasterDict, cellSizeNew, remeshOption="griddata", interpMethod="c
     return remeshedRaster
 
 
-def remeshDataRio(rasterFile, cellSizeNew, larger=True):
+def remeshDataRio(rasterFile, cellSizeNew, resamplingM="default", larger=True):
     """resample raster data using rasterio to change effective cell size to cellSizeNew by specifying an output array
     of specified size, default resampling option is set to cubic
 
@@ -316,6 +344,8 @@ def remeshDataRio(rasterFile, cellSizeNew, larger=True):
         path to file with raster data options asci or tif
     cellSizeNew: float
         desired spatial resolution of new raster dataset, cellsize
+    resamplingM: str
+        options are: default: cubic resampling, nearest
     larger: Boolean
         if true (default) output grid is at least as big as the input
 
@@ -347,6 +377,16 @@ def remeshDataRio(rasterFile, cellSizeNew, larger=True):
     else:
         srcCrs = src.crs
 
+    # set resampling method
+    if resamplingM == "default":
+        resamplingMethod = Resampling.cubic
+    elif resamplingM == "nearest":
+        resamplingMethod = Resampling.nearest
+    else:
+        message = 'Interpolation method "%s" not recognized' % (resamplingM)
+        log.error(message)
+        raise NameError(message)
+
     data, transform = rasterio.warp.reproject(
         source=src.read(),
         destination=np.empty((src.count, height, width)) * np.nan,
@@ -356,7 +396,7 @@ def remeshDataRio(rasterFile, cellSizeNew, larger=True):
         dst_crs=srcCrs,
         src_nodata=src.nodata,
         dst_nodata=src.nodata,
-        resampling=Resampling.cubic,
+        resampling=resamplingMethod,
     )
 
     data = np.where(data == src.nodata, np.nan, data)
@@ -373,6 +413,7 @@ def remeshDataRio(rasterFile, cellSizeNew, larger=True):
 
     # create remeshed raster dictionary
     remeshedRaster = {"rasterData": data[0], "header": headerRemeshed}
+    log.info('Remeshing of %s complete using resampling method "%s"' % (rasterFile.name, resamplingM))
 
     return remeshedRaster
 
@@ -430,13 +471,19 @@ def remeshRaster(rasterFile, cfgSim, typeIndicator="DEM", onlySearch=False, lega
     )
     if legacy:
         remeshedRaster = remeshData(
-            raster, cszRasterNew, remeshOption="griddata", interpMethod="cubic", larger=False
+            raster,
+            cszRasterNew,
+            remeshOption="griddata",
+            interpMethod=cfgSim["GENERAL"]["remeshInterpMethod"],
+            larger=False,
         )
         log.info("Legacy option used for remeshing")
         flipArg = True
     else:
         log.info("Using rasterio resampling")
-        remeshedRaster = remeshDataRio(rasterFile, cszRasterNew, larger=False)
+        remeshedRaster = remeshDataRio(
+            rasterFile, cszRasterNew, cfgSim["GENERAL"]["remeshInterpMethod"], larger=False
+        )
         flipArg = False
 
     # save remeshed raster
@@ -2129,3 +2176,51 @@ def checkDBOverlap(DBXl, DBXr, DBYl, DBYr):
     if not DBrLine.is_simple or not DBlLine.is_simple:
         message = "Domain transformation for given path_aimec - curvature of provided line leads to folding"
         log.warning(message)
+
+
+def interpolateLineLinear(lineDict, distance):
+    """interpolate a lines x, y coordinates using numpy.interp performing a one-dimensional piecewise linear
+    interpolation
+
+    Parameters
+    ------------
+    lineDict: dict
+        dictionary with x, y coordinates
+    distance: float
+        distance between new points
+
+    Returns
+    ---------
+    lineDict: dict
+        updated dictionary with x, y coordinates
+
+    """
+
+    # fetch x, y coors from lineDict
+    x = lineDict["x"]
+    y = lineDict["y"]
+
+    # check if duplicate points in lineDict coordinates
+    indexNonDup = np.where(np.abs(np.diff(x)) + np.abs(np.diff(y)) > 0)
+    xcoor = x[indexNonDup]
+    xNew = np.append(xcoor, x[-1])
+    ycoor = y[indexNonDup]
+    yNew = np.append(ycoor, y[-1])
+
+    s = computeLengthOfLine2D(xNew, yNew)
+    s = np.append(0, s)
+    nPoints = int(np.ceil(s[-1] / distance))
+
+    # create a count for the number of given coordinate points in lineDict
+    xPoints = np.linspace(0, s[-1], nPoints)
+    yPoints = np.linspace(0, s[-1], nPoints)
+
+    # perform interpolation on each coordinate
+    X1 = np.interp(xPoints, s, xNew)
+    Y1 = np.interp(yPoints, s, yNew)
+
+    # set new x, y coordinates
+    lineDict["x"] = X1
+    lineDict["y"] = Y1
+
+    return lineDict
diff --git a/avaframe/out3Plot/outDebugPlots.py b/avaframe/out3Plot/outDebugPlots.py
index 00e103508..e9cc68b76 100644
--- a/avaframe/out3Plot/outDebugPlots.py
+++ b/avaframe/out3Plot/outDebugPlots.py
@@ -422,13 +422,59 @@ def plotParticlesRelease(particles, relRaster, releaseLine, dem, cfg, xyParticle
         xyParticlesAll["x"] + dem["originalHeader"]["xllcenter"],
         xyParticlesAll["y"] + dem["originalHeader"]["yllcenter"],
         "+g",
-    )
+        )
     # only particles that have not been removed
     ax.plot(
         particles["x"] + dem["originalHeader"]["xllcenter"],
         particles["y"] + dem["originalHeader"]["yllcenter"],
         "*r",
-    )
+        )
     ax.plot(releaseLine["x"], releaseLine["y"], "-b")
     ax.set_title("mass/rho: %.2fm3" % (volParticles))
     plt.show()
+
+
+def plotParticleTrajOnGrid(x, y, xNew, yNew, dem):
+    """plot particle trajectory old and new (interpolated) on grid
+        current use: interpolateParticlesTrajectories in com1DFA/particleTools.py
+
+    Parameters
+    -----------
+    x, y, xNew, yNew : array
+        old and new (interpolated) coordinates of particle trajectories
+    dem: dict
+        dictionary with header and rasterData of computational DEM 
+    """
+    extentCellCenters, extentCellCorners = pU.createExtentMinMax(
+        dem["rasterData"], dem["header"], originLLCenter=True
+    )
+
+    # figure
+    fig, ax = plt.subplots(nrows=1, ncols=1)
+    # Minor ticks
+    ax.set_xticks(
+        np.arange(extentCellCorners[0], extentCellCorners[1], dem["header"]["cellsize"]), minor=True
+    )
+    ax.set_yticks(
+        np.arange(extentCellCorners[2], extentCellCorners[3], dem["header"]["cellsize"]), minor=True
+    )
+    # Gridlines based on minor ticks
+    ax.grid(which="minor", color="gray", linestyle="-", linewidth=2)
+    ax.plot(
+        x + dem["header"]["xllcenter"],
+        y + dem["header"]["yllcenter"],
+        "-",
+        color="blue",
+        linewidth=5,
+        alpha=0.25,
+    )
+    ax.plot(
+        x + dem["header"]["xllcenter"],
+        y + dem["header"]["yllcenter"],
+        "+",
+        color="blue",
+        markersize=15,
+    )
+    ax.plot(xNew + dem["header"]["xllcenter"], yNew + dem["header"]["yllcenter"], "-*r")
+    ax.axis("equal")
+    plt.show()
diff --git a/avaframe/out3Plot/outParticlesAnalysis.py b/avaframe/out3Plot/outParticlesAnalysis.py
index 2ad479a10..4e106220d 100644
--- a/avaframe/out3Plot/outParticlesAnalysis.py
+++ b/avaframe/out3Plot/outParticlesAnalysis.py
@@ -838,3 +838,101 @@ def readMeasuredParticleData(avalancheDir, demHeader, pData=""):
     mParticles["y"] = mParticles["y"] - demHeader["yllcenter"]
 
     return mParticles
+
+
+def plotParticlesAssets(dem, assets, particleAssets, outDir, plotName, title):
+    """Plot locations of particles over all timesteps that reached assets color coded with highes assets class
+
+    Parameters
+    -----------
+    dem: dict
+        dictionary with dem info
+    assets: dict
+        dictionary with assets info
+    particleAssets: numpy ndarray
+        array with particle trajectories colorcoded with assets classes
+    outDir: pathlib.Path
+        path to output directory
+    plotName: str
+        name of plot file
+    title: str
+        plot title
+    """
+    extentCellCenters, extentCellCorners = pU.createExtentMinMax(
+        dem["rasterData"], dem["header"], originLLCenter=True
+    )
+
+    fig, ax = plt.subplots(ncols=1)
+    ax.set_title(title)
+    # add DEM hillshade with contour lines
+    # set extent in meters using cellSize and llcenter location
+    _, _ = pU.addHillShadeContours(ax, dem["rasterData"], dem["header"]["cellsize"], extentCellCenters)
+
+    # create common colormap for assets info and particles assets info
+    vMin = np.nanmin(assets["rasterData"])
+    vMax = np.nanmax(assets["rasterData"])
+    cmap1, colorsNew, levelsNew, norm = pU.makeColorMap(
+        {"cmap": cm.hawaii.reversed()}, vMin, vMax, continuous=True
+    )
+    # set all cells that have been affected by particles but don't belong to the identified trajectories affecting assets
+    # to white
+    cmap1.set_under("white")
+    ax.imshow(
+        np.where(assets["rasterData"] != 0, assets["rasterData"], np.nan),
+        alpha=1.0,
+        extent=extentCellCenters,
+        origin="lower",
+        cmap=cmap1,
+        zorder=11,
+        vmin=vMin,
+        vmax=vMax,
+    )
+
+    im1 = ax.imshow(
+        np.where(particleAssets != 0, particleAssets, np.nan),
+        extent=extentCellCenters,
+        origin="lower",
+        alpha=0.6,
+        cmap=cmap1,
+        zorder=10,
+        vmin=vMin,
+        vmax=vMax,
+    )
+    # separate plot so that alpha=1
+    ax.imshow(
+        np.where(particleAssets == -1.0, particleAssets, np.nan),
+        extent=extentCellCenters,
+        origin="lower",
+        alpha=0.2,
+        cmap=cmap1,
+        zorder=10,
+        vmin=vMin,
+        vmax=vMax,
+    )
+    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)))
+
+
+def checkSavingTimeStepParticles(timeStepInfo, limitValue=2.0):
+    """check if saved particles are of high enough temporal resolution
+
+    Parameters
+    -----------
+    timeStepInfo: list, np array
+        list of saved time steps
+    limitValue: float
+        allowed difference between consecutive saved time steps
+    """
+    deltaT = np.diff(timeStepInfo)
+    if np.any(deltaT > limitValue):
+        message = (
+            "Saving time step of simulation particle Info exceeds two seconds - this can lead to errors in analysis"
+            "set saving time step to <= %.2f" % limitValue
+        )
+        log.error(message)
+        raise AssertionError(message)
diff --git a/avaframe/runScripts/runParticlesAssetsInfo.py b/avaframe/runScripts/runParticlesAssetsInfo.py
new file mode 100644
index 000000000..163953fc8
--- /dev/null
+++ b/avaframe/runScripts/runParticlesAssetsInfo.py
@@ -0,0 +1,112 @@
+"""
+Run creating a raster with numerical particle trajectories colorcoded with assets classes, highest overrides lower classes
+"""
+
+import pathlib
+import time
+import configparser
+
+# 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
+
+# +++++++++REQUIRED+++++++++++++
+# if particle locations are saved e.g. every second, the resulting assets raster might
+# show gaps as particles travelled further within this time step, to avoid these gaps
+# option to perform interpolation - can lead to errors if particle locations too spaced out
+# recommendation is particle saving time step every second
+interpolateParticlesTrajectoriesFlag = True
+# if interpolateParticlesTrajectoriesFlag the mesh cellsize x cellSizeFactor will give desired
+# distance of interpolated points
+cellSizeFactor = 0.5
+resizeThreshold = 3
+meshCellSizeThreshold = 0.001
+useCompression = True
+remeshInterpMethod = "nearest"
+# ++++++++++++++++++++++++++++++
+
+# load avalanche directory
+cfgMain = cfgUtils.getGeneralConfig()
+avalancheDir = cfgMain["MAIN"]["avalancheDir"]
+outDir = pathlib.Path(avalancheDir, "Outputs", "out3Plot", "particleAnalysis")
+fU.makeADir(outDir)
+cfg = configparser.ConfigParser()
+cfg["GENERAL"] = {
+    "avalancheDir": avalancheDir,
+    "meshCellSizeThreshold": meshCellSizeThreshold,
+    "resizeThreshold": resizeThreshold,
+    "remeshInterpMethod": remeshInterpMethod,
+}
+cfg["EXPORTS"] = {"useCompression": useCompression}
+
+# 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():
+    startTime = time.time()
+    # fetch simName
+    simName = row["simName"]
+    dem = rU.readRaster(pathlib.Path(avalancheDir, "Inputs", row["DEM"]))
+    log.info("Find particle trajectories for simulation: %s" % (simName))
+    # add info on meshCellSize
+    cfg["GENERAL"]["meshCellSize"] = str(dem["header"]["cellsize"])
+
+    # fetch info on infrastructure
+    uniqueAssets, assets, assetsValues = gI.preprocessAssets(avalancheDir, dem, cfg)
+
+    # read particles saved from com1DFA simulation
+    Particles, timeStepInfo = pT.readPartFromPickle(
+        pathlib.Path(avalancheDir), simName=simName, flagAvaDir=True, comModule="com1DFA"
+    )
+
+    # Error if time step is larger than 2 seconds
+    oP.checkSavingTimeStepParticles(timeStepInfo)
+
+    # create time series of particles arrays
+    particlesTimeArrays = pT.reshapeParticlesDicts(
+        Particles, ["ID", "indXDEM", "indYDEM", "x", "y", "z", "inCellDEM"]
+    )
+
+    if interpolateParticlesTrajectoriesFlag:
+        # interpolate particle trajectories +
+        pLong = pT.interpolateParticlesTrajectories(dem, particlesTimeArrays, cellSizeFactor)
+        particleTimeInfo = pLong.copy()
+    else:
+        particleTimeInfo = particlesTimeArrays.copy()
+
+    # derive info on which particles interacted with infrastructure
+    particleAssets, particleTimeInfo = pT.createAssetsRasterFromParticleLocations(
+        particleTimeInfo, dem, uniqueAssets, assetsValues
+    )
+    # export raster
+    rU.writeResultToRaster(
+        dem["header"], particleAssets, (outDir / ("particleAssetsInfo_%s" % simName)), flip=True
+    )
+
+    # create plot
+    plotName = "particleAssetsInfo_%s" % simName
+    if interpolateParticlesTrajectoriesFlag:
+        plotTitle = "Particle trajectories (interpolated) color-coded with asset classes"
+    else:
+        plotTitle = "Particle trajectories color-coded with asset classes"
+    _ = oP.plotParticlesAssets(dem, assets, particleAssets, outDir, plotName, plotTitle)
+
+    timeNeeded = "%.2f" % (time.time() - startTime)
+    log.info("computation took: %s s " % timeNeeded)
diff --git a/avaframe/tests/test_deriveParameterSet.py b/avaframe/tests/test_deriveParameterSet.py
index 3c1315c88..0898c95c6 100644
--- a/avaframe/tests/test_deriveParameterSet.py
+++ b/avaframe/tests/test_deriveParameterSet.py
@@ -609,7 +609,12 @@ def test_checkExtentAndCellSize(tmp_path):
     # setup required inputs
     testDir = pathlib.Path(tmp_path, "test")
     cfg = configparser.ConfigParser()
-    cfg["GENERAL"] = {"resizeThreshold": 3.0, "meshCellSize": 1.0, "meshCellSizeThreshold": 0.0001}
+    cfg["GENERAL"] = {
+        "resizeThreshold": 3.0,
+        "meshCellSize": 1.0,
+        "meshCellSizeThreshold": 0.0001,
+        "remeshInterpMethod": "default",
+    }
     cfg["GENERAL"]["avalancheDir"] = str(testDir)
     cfg["EXPORTS"] = {"useCompression": "True"}
     inDir = testDir / "Inputs"
@@ -662,6 +667,14 @@ def test_checkExtentAndCellSize(tmp_path):
     assert remeshedFlag == "Yes"
     assert outFile.name == testFile.split("/")[1]
 
+    cfg["GENERAL"]["remeshInterpMethod"] = "bilinear"
+    with pytest.raises(NameError) as e:
+        assert dP.checkExtentAndCellSize(cfg, inputFile, dem, "mu")
+    assert 'Interpolation method "%s" not recognized' % (cfg["GENERAL"]["remeshInterpMethod"]) in str(
+        e.value
+    )
+
+    cfg["GENERAL"]["remeshInterpMethod"] = "default"
     inputFile2 = inDirR / "inputFile1.asc"
     headerInput2 = {
         "nrows": 4,
@@ -798,6 +811,70 @@ def test_checkExtentAndCellSize(tmp_path):
     testFile5, outFile5, remeshedFlag5 = dP.checkExtentAndCellSize(cfg, inputFile, dem, "DEM")
     assert remeshedFlag5 == "No"
 
+    # test nearest meshing
+    inDirR = pathlib.Path(tmp_path, "avaTestRemeshing")
+    fU.makeADir(inDirR)
+
+    # configuration settings
+    cfg = configparser.ConfigParser()
+    cfg["GENERAL"] = {
+        "resizeThreshold": 3.0,
+        "meshCellSize": 2.0,
+        "meshCellSizeThreshold": 0.0001,
+        "remeshInterpMethod": "default",
+    }
+    cfg["GENERAL"]["avalancheDir"] = str(inDirR)
+    cfg["EXPORTS"] = {"useCompression": "True"}
+    cfg["GENERAL"]["remeshInterpMethod"] = "nearest"
+
+    # create DEM
+    demField = np.ones((10, 12))
+    dem = {
+        "header": {
+            "nrows": 10,
+            "ncols": 12,
+            "xllcenter": 1,
+            "yllcenter": 5,
+            "cellsize": 2,
+            "nodata_value": -9999,
+            "driver": "AAIGrid",
+        },
+        "rasterData": demField,
+    }
+    dem["header"]["transform"] = IOf.transformFromASCHeader(dem["header"])
+    dem["header"]["crs"] = rasterio.crs.CRS()
+    IOf.writeResultToRaster(dem["header"], demField, inDirR / "demTest", flip=False)
+
+    # create inputField
+    inputFile = inDirR / "inputFile.asc"
+    headerInput = {
+        "nrows": 5,
+        "ncols": 6,
+        "xllcenter": 1.1,
+        "yllcenter": 5.1,
+        "cellsize": 4,
+        "nodata_value": -9999,
+        "driver": "AAIGrid",
+    }
+
+    headerInput["transform"] = IOf.transformFromASCHeader(headerInput)
+    headerInput["crs"] = rasterio.crs.CRS()
+    inField = np.ones((5, 6))
+    inField[2, 4] = 10.0
+    IOf.writeResultToRaster(headerInput, inField, inputFile.parent / inputFile.stem, flip=False)
+
+    # remesh inputFile to match DEM extent and cellSize
+    outFilePath, _, _ = dP.checkExtentAndCellSize(cfg, inputFile, dem, "mu", nanInsideDEMCheck=False)
+    outFile = IOf.readRaster((inDirR / "Inputs" / outFilePath))
+    testRaster = np.ones((10, 12))
+    testRaster[4:6, 8:10] = 10.0
+
+    assert np.array_equal(testRaster, outFile["rasterData"])
+    assert outFile["header"]["nrows"] == 10
+    assert outFile["rasterData"].shape[1] == 12
+    assert outFile["header"]["xllcenter"] == 1.0
+    assert outFile["header"]["yllcenter"] == 5.0
+
 
 # Produced by AI (test):
 
diff --git a/avaframe/tests/test_geoTrans.py b/avaframe/tests/test_geoTrans.py
index 95e9b578b..ea81ced81 100644
--- a/avaframe/tests/test_geoTrans.py
+++ b/avaframe/tests/test_geoTrans.py
@@ -603,6 +603,7 @@ def test_remeshDEM(tmp_path):
         "meshCellSizeThreshold": "0.0001",
         "meshCellSize": "2.",
         "avalancheDir": str(avaDir),
+        "remeshInterpMethod": "default",
     }
 
     # call function
@@ -667,6 +668,47 @@ def test_remeshDEM(tmp_path):
     assert dataNew2["rasterData"].shape[1] == dataSol["header"]["ncols"]
     assert testRes2
 
+    cfg["GENERAL"]["remeshInterpMethod"] = "bilinear"
+    cfg["GENERAL"]["meshCellSize"] = "7."
+    with pytest.raises(NameError) as e:
+        assert geoTrans.remeshRaster(avaDEM1, cfg, legacy=True)
+    assert 'Interpolation method "%s" not recognized' % (cfg["GENERAL"]["remeshInterpMethod"]) in str(
+        e.value
+    )
+
+    # copy input data for remeshing
+    avaDir4 = pathlib.Path(tmp_path, "avaTestMesh")
+    inputDir1 = dirPath / ".." / "data" / "avaParabola"
+    inputDEM1 = inputDir1 / "Inputs" / "DEM_PF_Topo.asc"
+    fU.makeADir((avaDir4 / "Inputs"))
+    avaDEM4 = avaDir4 / "Inputs" / "DEM_PF_Topo.asc"
+    shutil.copy(inputDEM1, avaDEM4)
+
+    cfg["GENERAL"]["remeshInterpMethod"] = "cubic"
+    cfg["GENERAL"]["meshCellSize"] = "8."
+    cfg["GENERAL"]["avalancheDir"] = str(avaDir4)
+
+    testRes4 = geoTrans.remeshRaster(avaDEM4, cfg, legacy=True)
+    fullP2 = avaDir1 / "Inputs" / testRes4
+
+    dataNew4 = IOf.readRaster(fullP2)
+    dataSol = IOf.readRaster(inputDEM)
+
+    # compare solution to result from function
+    testRes44 = np.allclose(dataNew4["rasterData"], dataSol["rasterData"], atol=1.0e-6)
+
+    assert dataNew4["rasterData"].shape[0] == dataSol["header"]["nrows"]
+    assert dataNew4["rasterData"].shape[1] == dataSol["header"]["ncols"]
+    assert testRes44
+
+    cfg["GENERAL"]["remeshInterpMethod"] = "bilinear"
+    cfg["GENERAL"]["meshCellSize"] = "4.5"
+    with pytest.raises(NameError) as e:
+        assert geoTrans.remeshRaster(avaDEM1, cfg, legacy=False)
+    assert 'Interpolation method "%s" not recognized' % (cfg["GENERAL"]["remeshInterpMethod"]) in str(
+        e.value
+    )
+
 
 def test_isCounterClockWise():
     """test isCounterClockWise"""
diff --git a/avaframe/tests/test_getInput.py b/avaframe/tests/test_getInput.py
index 553556fec..2f504c3c8 100644
--- a/avaframe/tests/test_getInput.py
+++ b/avaframe/tests/test_getInput.py
@@ -19,6 +19,7 @@
 import avaframe.in3Utils.geoTrans as geoTrans
 import avaframe.com1DFA.DFAtools as DFAtls
 from avaframe.com1DFA import com1DFA
+import avaframe.in2Trans.rasterUtils as rU
 import logging
 
 
@@ -1650,3 +1651,74 @@ def test_checkTimeDepRelease():
     with pytest.raises(ValueError) as e:
         getInput.checkTimeDepRelease(timeDepRelValues, timeDepRelCsv)
     assert ("The initial velocity provided in %s can not be negative." % (timeDepRelCsv)) in str(e.value)
+
+
+def test_preprocessAssets(tmp_path):
+    """test creating asset info"""
+
+    # setup required inputs
+
+    dem = {
+        "header": {
+            "xllcenter": 0,
+            "yllcenter": 0,
+            "cellsize": 2,
+            "nrows": 10,
+            "ncols": 11,
+            "driver": "AAIGrid",
+            "nodata_value": np.nan,
+            "crs": None,
+        },
+        "rasterData": np.array(
+            [
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+            ]
+        ),
+    }
+    transform = rU.transformFromASCHeader(dem["header"])
+    dem["header"]["transform"] = transform
+
+    avaTestDir = pathlib.Path(tmp_path, "avaTest")
+    assetsDir = avaTestDir / "Inputs" / "INFRA"
+    fU.makeADir(assetsDir)
+    assetsFile = assetsDir / "assets_ASSETS"
+    assetsArray = np.full((dem["header"]["nrows"], dem["header"]["ncols"]), np.nan)
+    assetsArray[0, 0] = 1.0
+    assetsArray[1, 0:2] = 2.0
+    assetsArray[2, 0] = 3.0
+    assetsArray[4, 0] = 0.0
+    rU.writeResultToRaster(dem["header"], assetsArray, assetsFile, flip=True)
+
+    cfg = configparser.ConfigParser()
+    cfg["GENERAL"] = {
+        "avalancheDir": avaTestDir,
+        "meshCellSizeThreshold": "0.0001",
+        "resizeThreshold": "3.",
+        "remeshInterpMethod": "nearest",
+    }
+    cfg["EXPORTS"] = {"useCompression": "True"}
+    # call function to be tested
+    uniqueAssets, assets, assetsValues = getInput.preprocessAssets(avaTestDir, dem, cfg)
+
+    assert uniqueAssets == [1.0, 2.0, 3.0]
+    assert assets["header"]["xllcenter"] == dem["header"]["xllcenter"]
+    assert np.array_equal(assets["rasterData"], assetsArray, equal_nan=True)
+    assert assetsValues["value_1"] == [0]
+    assert assetsValues["value_2"] == [11, 12]
+    assert assetsValues["value_3"] == [22]
+
+    assetsArray[5, 0] = -1.0
+    rU.writeResultToRaster(dem["header"], assetsArray, assetsFile, flip=True)
+
+    with pytest.raises(AssertionError) as e:
+        getInput.preprocessAssets(avaTestDir, dem, cfg)
+    assert "In ASSETS file (assets_ASSETS.asc) negative values found - this is not allowed" in str(e.value)
diff --git a/avaframe/tests/test_outParticlesAnalysis.py b/avaframe/tests/test_outParticlesAnalysis.py
index c63ca27c1..27effce7c 100644
--- a/avaframe/tests/test_outParticlesAnalysis.py
+++ b/avaframe/tests/test_outParticlesAnalysis.py
@@ -6,6 +6,7 @@
 import pathlib
 import configparser
 import matplotlib.pyplot as plt
+import pytest
 
 # Local imports
 import avaframe.out3Plot.outParticlesAnalysis as oA
@@ -95,3 +96,22 @@ def test_velocityEnvelopeThalweg():
     assert np.array_equal(dictVelAltThalweg["medianVelocity"], velMagMean)
     assert np.array_equal(dictVelAltThalweg["maxSxyz"], sxyzMax)
     assert np.array_equal(dictVelAltThalweg["minSxyz"], sxyzMin)
+
+
+def test_checkSavingTimeStepParticles():
+    """check if time step exceeds limitValue"""
+
+    # setup input
+    timeStepInfo = [0.0, 1.999, 2.0, 4.0, 4.5]
+    limitValue = 2.0
+
+    oA.checkSavingTimeStepParticles(timeStepInfo, limitValue=limitValue)
+
+    timeStepInfo = [0.0, 2.0001, 2.4, 4.0, 4.5]
+
+    with pytest.raises(AssertionError) as e:
+        oA.checkSavingTimeStepParticles(timeStepInfo, limitValue=limitValue)
+    assert (
+        "aving time step of simulation particle Info exceeds two seconds - this can lead to errors in analysis"
+        in str(e.value)
+    )
diff --git a/avaframe/tests/test_particleTools.py b/avaframe/tests/test_particleTools.py
index 3cb5820b7..d4138f043 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,155 @@ 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
+    particleAssetsTest[4, 7] = -1.0
+    particleAssetsTest[3, 8] = -1.0
+    particleAssetsTest[2, 9] = -1.0
+
+    assert particleAssets.shape == (8, 10)
+    assert np.array_equal(particleAssets, particleAssetsTest, equal_nan=True)
+
+
+def findUniqueCellIndices():
+    """test finding indices that are unique in particleTimeArrays indXDEM and indYDEM"""
+
+    # setup required input
+    indXDEM = np.asarray([[1, 2, 3], [4, 5, 6], [1, 2, 10], [4, 8, 9]])
+    indYDEM = np.asarray([[10, 20, 30], [40, 50, 60], [10, 20, 50], [40, 80, 90]])
+    particleTimeArrays = {"indXDEM": indXDEM, "indYDEM": indYDEM}
+
+    xyIndAllUnique = particleTools.findUniqueCellIndices(particleTimeArrays)
+
+    testArray1 = np.asarray(
+        [[1, 10], [2, 20], [3, 30], [4, 40], [5, 50], [6, 60], [10, 50], [8, 80], [9, 90]]
+    )
+    testArray = np.asarray(testArray1, dtype=int)
+
+    assert np.array_equal(xyIndAllUnique, testArray, equal_nan=True)
+
+
+def test_interpolateParticlesTrajectories():
+    """test if interpolation is done at desired distances and that start and end are kept"""
+
+    dem = {
+        "header": {"xllcenter": 0, "yllcenter": 0, "cellsize": 2, "nrows": 10, "ncols": 11},
+        "rasterData": np.array(
+            [
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+                [50, 40, 30, 20, 10, 0, 0, 0, 0, 0, 0],
+            ]
+        ),
+    }
+    particlesTimeArrays = {"x": np.asarray([[1, 2], [2, 3], [3, 4], [8, 9]])}
+    particlesTimeArrays["y"] = np.asarray([[1, 2], [2, 3], [3.7, 4], [8, 9]])
+    particlesTimeArrays["z"] = np.asarray([[40, 40], [30, 30], [20, 20], [0, 0]])
+    particlesTimeArrays["ID"] = np.asarray([[1, 2], [1, 2], [1, 2], [1, 2]])
+
+    pLong = particleTools.interpolateParticlesTrajectories(dem, particlesTimeArrays, 0.5, debugPlot=False)
+
+    assert np.array_equal(pLong["2_pTraj"]["x"], pLong["2_pTraj"]["y"])
+    assert (np.diff(pLong["2_pTraj"]["s"]) < (0.5 * dem["header"]["cellsize"])).any()
+    assert pLong["1_pTraj"]["x"][0] == particlesTimeArrays["x"][0, 0]
+    assert pLong["1_pTraj"]["x"][-1] == particlesTimeArrays["x"][-1, 0]
+    assert pLong["2_pTraj"]["y"][0] == particlesTimeArrays["y"][0, 1]
+    assert pLong["2_pTraj"]["y"][-1] == particlesTimeArrays["y"][-1, 1]
+    assert pLong["1_pTraj"]["y"][0] == particlesTimeArrays["y"][0, 0]
+    assert pLong["1_pTraj"]["y"][-1] == particlesTimeArrays["y"][-1, 0]
+    assert pLong["2_pTraj"]["x"][0] == particlesTimeArrays["x"][0, 1]
+    assert pLong["2_pTraj"]["x"][-1] == particlesTimeArrays["x"][-1, 1]
+
+    pLong2 = particleTools.interpolateParticlesTrajectories(dem, particlesTimeArrays, 1, debugPlot=False)
+
+    assert np.array_equal(pLong2["indXDEM"][:, 1], np.asarray([1, 2, 2, 3, 4, 4]))
+    assert np.array_equal(pLong2["indYDEM"][:, 1], np.asarray([1, 2, 2, 3, 4, 4]))
diff --git a/docs/moduleOut3Plot.rst b/docs/moduleOut3Plot.rst
index 7bd6e5ff3..ae12156ca 100644
--- a/docs/moduleOut3Plot.rst
+++ b/docs/moduleOut3Plot.rst
@@ -245,3 +245,33 @@ 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`` (see Note) 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|z|inCellDEM|nPart*.
+In addition, an assets raster file in ``avalancheDir/Inputs/INFRA`` is required. Classes have to be > 0,
+negative and zero values are treated as no data values. Preferably, the extent and resolution of the provided
+assets raster should match the extent and resolution of the simulation DEM. If extents or resolution do not match,
+remeshing will be performed. However, this can potentially introduce geometrical artefacts in the assets layer
+(a corresponding warning will be written to the log-file). In order to avoid introducing new classes as a result
+of interpolation, the default setting in the corresponding run script (parameter ``remeshInterpMethod``)
+is using a nearest-neighbor based interpolation. To perform the analysis (requires a prior :py:mod:`com1DFA.com1DFA`
+simulation run using the settings described above):
+run:
+  ::
+
+    pixi run python runScripts/runParticlesAssetsInfo.py
+
+.. Note::
+    The setting of the saving time step ``tSteps`` has a strong effect on the results of the assets analysis.
+    Choosing a saving time step close to the computational time step dt (default value is 0.1s), will result in
+    most accurate results. When choosing a significantly larger time step, derived particle trajectories will lead
+    to gaps in the analysis, hence cells will not be attributed and color-coded correctly. Interpolation of particle
+    trajectories (setting: ``interpolateParticlesTrajectoriesFlag = True``), can help to determine which cells are
+    affected but also with this option, if a too large saving time step (> 1 second) is chosen, errors have to be
+    expected. Hence, the default setting is ``interpolateParticlesTrajectoriesFlag = True`` and if saving time steps
+    exceed 2 seconds an error is raised.
\ No newline at end of file