[out3Plot]: add functionality for assets

.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:

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"]))

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

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

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