diff --git a/avaframe/com1DFA/com1DFA.py b/avaframe/com1DFA/com1DFA.py index c5f15700a..e928523d3 100644 --- a/avaframe/com1DFA/com1DFA.py +++ b/avaframe/com1DFA/com1DFA.py @@ -620,7 +620,17 @@ def prepareInputData(inputSimFiles, cfg): # TODO: remove if not required anymore # relThFieldData, _ = gI.initializeRelTh(cfg, dOHeader) - if cfg["INPUT"]["relThFile"] == "": + if relFile.suffix.lower() == ".csv": + if not cfg["GENERAL"].getboolean("timeDependentRelease"): + message = "A CSV release geometry requires timeDependentRelease=True" + log.error(message) + raise ValueError(message) + timeDepRelValues, _ = gI.getTimeDepRelCsv(cfg["INPUT"]["timeDepRelCsv"]) + if "x" not in timeDepRelValues: + message = "A CSV release geometry requires x and y columns in the time dependent release CSV" + log.error(message) + raise ValueError(message) + elif cfg["INPUT"]["relThFile"] == "": # get line from release area polygon releaseLine = shpConv.readLine(relFile, "release1", demOri) releaseLine["file"] = relFile @@ -646,13 +656,24 @@ def prepareInputData(inputSimFiles, cfg): log.info("Set %s for relThField" % relRasterPath) # get line from release area polygon if cfg["GENERAL"].getboolean("timeDependentRelease"): - releaseLine["type"] = "time dependent Release" timeDepRelValues, _ = gI.getTimeDepRelCsv(cfg["INPUT"]["timeDepRelCsv"]) - releaseLine["thickness"] = [ - timeDepRelValues["thickness"][timeDepRelValues["timeStep"] == 0].item() - ] * len(releaseLine["Name"]) + + if "x" in timeDepRelValues: + releaseLine, relThFieldData = debF.defineReleaseLineFromCoordinates(relFile, timeDepRelValues, dOHeader) + else: + releaseLine["thickness"] = [ + timeDepRelValues["thickness"][timeDepRelValues["timeStep"] == 0].item() + ] * len(releaseLine["Name"]) + + if "velocity" in timeDepRelValues.keys(): + releaseLine["velocity"] = timeDepRelValues["velocity"][timeDepRelValues["timeStep"] == 0] + else: + releaseLine["velocityX"] = timeDepRelValues["velocityX"][timeDepRelValues["timeStep"] == 0] + releaseLine["velocityY"] = timeDepRelValues["velocityY"][timeDepRelValues["timeStep"] == 0] + releaseLine["velocityZ"] = timeDepRelValues["velocityZ"][timeDepRelValues["timeStep"] == 0] + + releaseLine["type"] = "time dependent Release" releaseLine["thicknessSource"] = ["csv file"] * len(releaseLine["Name"]) - releaseLine["velocity"] = timeDepRelValues["velocity"][timeDepRelValues["timeStep"] == 0] releaseLine["timeDepRelValues"] = timeDepRelValues # get line from secondary release area polygon @@ -1258,8 +1279,9 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName): relThField=relThField, ) - if cfgGen.getboolean("timeDependentRelease") and releaseLine["velocity"] != 0: - particles = DFAfunC.updateInitialVelocity(cfgGen, particles, dem, releaseLine["velocity"]) + if cfgGen.getboolean("timeDependentRelease") and "velocity" in releaseLine.keys(): + if np.any(releaseLine["velocity"]) != 0: + particles = DFAfunC.updateInitialVelocity(cfgGen, particles, dem, releaseLine["velocity"]) particles, fields = initializeFields(cfg, dem, particles, releaseLine) reportAreaInfo["Release area info"]["Model release volume [m3]"] = "%.0f" % ( @@ -1449,7 +1471,7 @@ def initializeSimulation(cfg, outDir, demOri, inputSimLines, logName): return particles, fields, dem, reportAreaInfo -def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", relThField="", thName="rel"): +def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", relThField="", thName="rel", timestep=0.0): """Initialize DFA simulation Create particles and fields dictionary according to config parameters @@ -1469,6 +1491,8 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel if the release thickness is not uniform, give here the releaseRaster thName: str name rel, secondaryRel + timestep: float + for log: timestep at which particles are released Returns ------- @@ -1631,6 +1655,7 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel not cfg.getboolean("iniStep") and not cfg.getboolean("initialiseParticlesFromFile") and len(relThField) == 0 + and releaseLine["initializedFrom"] != "csvfile" ): if debugPlot: xyParticlesAll = {"x": particles["x"], "y": particles["y"]} @@ -1667,8 +1692,8 @@ def initializeParticles(cfg, releaseLine, dem, inputSimLines="", logName="", rel particles["nPPK"] = nPPK log.info( - "Initialized particles. MTot = %.2f kg, %s particles in %.2f cells." - % (particles["mTot"], particles["nPart"], relCells) + "Initialized particles in t = %.2f s. MTot = %.2f kg, %s particles in %.2f cells." + % (timestep, particles["mTot"], particles["nPart"], relCells) ) log.info( "Mass per particle = %.2f kg and particles per cell = %.2f." @@ -2860,7 +2885,8 @@ def releaseSecRelArea(cfg, particles, fields, dem, zPartArray0, reportAreaInfo): if secondaryReleaseInfo["initializedFrom"] == "shapefile": secRelInfo = shpConv.extractFeature(secondaryReleaseInfo, count) secRelInfo["rasterData"] = secRelRaster - secRelParticles = initializeParticles(cfg, secRelInfo, dem, thName="secondaryRel") + secRelParticles = initializeParticles(cfg, secRelInfo, dem, thName="secondaryRel", + timestep=particles["t"]) else: secondaryReleaseInfo["rasterData"] = secRelRaster secRelParticles = initializeParticles( @@ -2869,6 +2895,7 @@ def releaseSecRelArea(cfg, particles, fields, dem, zPartArray0, reportAreaInfo): dem, relThField=secRelRaster, thName="secondaryRel", + timestep=particles["t"] ) # release secondary release area by just appending the particles log.info( @@ -3436,7 +3463,12 @@ def prepareVarSimDict(standardCfg, inputSimFiles, variationDict, simNameExisting timeDepRelValues, _ = gI.getTimeDepRelCsv(cfgSim["INPUT"]["timeDepRelCsv"]) cfgSim["INPUT"]["timeDepRelTimeStep"] = str(timeDepRelValues["timeStep"]) cfgSim["INPUT"]["timeDepRelThickness"] = str(timeDepRelValues["thickness"]) - cfgSim["INPUT"]["timeDepRelVelocity"] = str(timeDepRelValues["velocity"]) + if "velocity" in timeDepRelValues: + cfgSim["INPUT"]["timeDepRelVelocity"] = str(timeDepRelValues["velocity"]) + else: + cfgSim["INPUT"]["timeDepRelVelocityX"] = str(timeDepRelValues["velocityX"]) + cfgSim["INPUT"]["timeDepRelVelocityY"] = str(timeDepRelValues["velocityY"]) + cfgSim["INPUT"]["timeDepRelVelocityZ"] = str(timeDepRelValues["velocityZ"]) else: cfgSim["INPUT"]["timeDepRelCsv"] = "" @@ -3851,10 +3883,20 @@ def initializeRelVol(cfg, demVol, releaseFile, radius, releaseType="primary", ti # check if release thickness provided as field or constant value if cfg["INPUT"][(typeTh + "File")] != "": - # read relThField from file - relThFilePath = pathlib.Path(cfg["GENERAL"]["avalancheDir"], "Inputs", cfg["INPUT"][typeTh + "File"]) - relThFieldFull = IOf.readRaster(relThFilePath) - relThField = relThFieldFull["rasterData"] + if releaseType == "timeDepRel": + # compute total initialized thickness + timeDepRelValues, _ = gI.getTimeDepRelCsv(timeDepRelFile) + # for time dependent release use the release volume summed up over all timesteps + relThField = np.zeros((demVol["header"]["nrows"], demVol["header"]["ncols"])) + for ts in np.unique(timeDepRelValues["timeStep"]): + idx = np.where(timeDepRelValues["timeStep"] == ts)[0] + thRaster = gI.timeDepRelCoordsToRaster(timeDepRelValues, idx, demVol["header"], parameter="thickness") + relThField += thRaster + else: + # read relThField from file + relThFilePath = pathlib.Path(cfg["GENERAL"]["avalancheDir"], "Inputs", cfg["INPUT"][typeTh + "File"]) + relThFieldFull = IOf.readRaster(relThFilePath) + relThField = relThFieldFull["rasterData"] # mask the relThField with raster from polygon releaseLineMask = np.ma.masked_where(relThField == 0.0, relThField) diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index 4b608f955..386282123 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/avaframe/com1DFA/com1DFACfg.ini @@ -133,8 +133,9 @@ entThDistVariation = entTh = #+++++++++++++General start conditions: time dependent release -# if timeDependentRelease is True (and relThFromFile is True), provide the the timesteps, thickness and velocity -# for a releases in a csv-file in the REL folder +# if timeDependentRelease is True (and relThFromFile is True), provide timestep, thickness and velocity +# in a csv file in Inputs/REL. Optional x and y columns define a one-cell release location +# for each timestep; when x/y are provided, no release polygon is required. timeDependentRelease = False # specify one or multiple particular time dependent release files, # provide name of csv file with or without extension .csv diff --git a/avaframe/com1DFA/debrisFunctions.py b/avaframe/com1DFA/debrisFunctions.py index 08410304d..f91f1c1e3 100644 --- a/avaframe/com1DFA/debrisFunctions.py +++ b/avaframe/com1DFA/debrisFunctions.py @@ -56,22 +56,20 @@ def initializeTimeDepRelease(cfg, inputSimLines, particles, fields, dem, zPartAr if np.isclose(t, timeDepRelValues["timeStep"], atol=atol, rtol=0).any(): iTup = np.where(np.isclose(t, timeDepRelValues["timeStep"], atol=atol, rtol=0)) # iTup is a tuple containing an array with one value in the first position, so we can extract the index: - i = iTup[0].item() - log.info( - "add release at timestep: %.2f s with thickness %s m and velocity %s m/s" - % (t, timeDepRelValues["thickness"][i], timeDepRelValues["velocity"][i]) - ) + i = iTup[0] + # similar workflow to secondary release! particles, zPartArray0 = addReleaseParticles( cfg, particles, inputSimLines, - timeDepRelValues["thickness"][i], - timeDepRelValues["velocity"][i], + timeDepRelValues, dem, zPartArray0, + timeDepRelIndex=i, ) particles = DFAfunC.getNeighborsC(particles, dem) + # update fields (compute grid values) if fields["computeTA"]: particles = DFAfunC.computeTrajectoryAngleC(particles, zPartArray0) @@ -80,7 +78,9 @@ def initializeTimeDepRelease(cfg, inputSimLines, particles, fields, dem, zPartAr return particles, fields, zPartArray0 -def addReleaseParticles(cfg, particles, inputSimLines, thickness, velocityMag, dem, zPartArray0): +def addReleaseParticles( + cfg, particles, inputSimLines, timeDepRelValues, dem, zPartArray0, timeDepRelIndex +): """ add new particles initialized by a time dependent release to particles that are in the flow already @@ -92,10 +92,8 @@ def addReleaseParticles(cfg, particles, inputSimLines, thickness, velocityMag, d particles dictionary at t that are in the flow already inputSimLines : dict dictionary with input data dictionaries (releaseLine,...) - thickness: float - thickness of current release - velocityMag: float - velocity of current release + timeDepRelValues: dict + time dependent release values dem: dict dictionary with info on DEM data zPartArray0: numpy array @@ -108,43 +106,106 @@ def addReleaseParticles(cfg, particles, inputSimLines, thickness, velocityMag, d zPartArray0: dict dictionary containing z - value of particles at timestep 0 """ - relLine = inputSimLines["releaseLine"] - relLine["header"] = dem["originalHeader"].copy() - relLine = geoTrans.prepareArea( - relLine, - dem, - np.sqrt(2), - thList=[thickness] * len(relLine["Name"]), - combine=True, - checkOverlap=False, - ) + thickness = timeDepRelValues["thickness"][timeDepRelIndex] + if "velocity" in timeDepRelValues: + velocityMag = timeDepRelValues["velocity"][timeDepRelIndex] - # check if already existing particles are within the release polygon - # it's possible that there are still a few particles in the polygon with low velocities - # TODO: could think of a threshold of number of particles that are still allowed in the polygons? - mask = geoTrans.getParticlesInPolygon(particles, relLine, cfg["GENERAL"].getfloat("thresholdPointInRel")) - if np.sum(mask) > 0: - # if there is at least one particle within the polygon (including the buffer): - message = ( - "Already existing particles are within the release polygon, which can cause numerical instabilities (at timestep: %02f s)" - % (particles["t"] + particles["dt"]) + relLine = copy.deepcopy(inputSimLines["releaseLine"]) + relLine["header"] = dem["originalHeader"].copy() + if relLine["initializedFrom"] == "csvfile": + relLine["rasterData"] = gI.timeDepRelCoordsToRaster( + relLine["timeDepRelValues"], timeDepRelIndex, dem["originalHeader"], parameter="thickness" ) - # timestep in particles is not updated yet - log.error(message) - raise ValueError(message) + relThField = relLine["rasterData"] + else: + relThField = "" + relLine = geoTrans.prepareArea( + relLine, + dem, + np.sqrt(2), + thList=[thickness] * len(relLine["Name"]), + combine=True, + checkOverlap=False, + ) + + # check if already existing particles are within the release polygon + # it's possible that there are still a few particles in the polygon with low velocities + # TODO: could think of a threshold of number of particles that are still allowed in the polygons? + mask = geoTrans.getParticlesInPolygon(particles, relLine, cfg["GENERAL"].getfloat("thresholdPointInRel")) + if np.sum(mask) > 0: + message = ( + "Already existing particles are within the release polygon, which can cause numerical instabilities (at timestep: %02f s)" + % (particles["t"] + particles["dt"]) + ) + log.error(message) + raise ValueError(message) particlesRelease = com1DFA.initializeParticles( cfg["GENERAL"], relLine, dem, + relThField=relThField, + timestep=particles["t"] + particles["dt"] ) - particlesRelease = DFAfunC.updateInitialVelocity(cfg["GENERAL"], particlesRelease, dem, velocityMag) + + if "velocity" in timeDepRelValues and "x" not in timeDepRelValues: + particlesRelease = DFAfunC.updateInitialVelocity(cfg["GENERAL"], particlesRelease, dem, velocityMag) + + + + elif "velocityX" in timeDepRelValues and "x" in timeDepRelValues: + for uComp, timedepParameter in zip(["ux", "uy", "uz"], ["velocityX", "velocityY", "velocityZ"]): + raster = gI.timeDepRelCoordsToRaster(relLine["timeDepRelValues"], timeDepRelIndex, + dem["originalHeader"], parameter=timedepParameter + ) + rasterDict = {"header": dem["header"], "rasterData": raster} + particlesRelease, _ = geoTrans.projectOnRaster(rasterDict, + particlesRelease, outData=uComp) + + particlesRelease["uMag"] = np.sqrt( + particlesRelease["ux"] ** 2 + particlesRelease["uy"] ** 2 + particlesRelease["uz"] ** 2) + particles = particleTools.mergeParticleDict(particles, particlesRelease) # save initial z position for travel angle computation zPartArray0 = np.append(zPartArray0, copy.deepcopy(particlesRelease["z"])) return particles, zPartArray0 +def defineReleaseLineFromCoordinates(relFile, timeDepRelValues, demHeader): + """ + define the release line and its thickness raster from coordinates read from time dependent release values + + Parameters + ---------- + relFile: pathlib.Path + directory to release file (csv file) + timeDepRelValues: dict + time dependent release values + demHeader: dict + header of DEM + + Returns + ------- + releaseLine: dict + dictionary for release line containing thickness raster data + relThFieldData: numpy array + release thickness raster data + """ + releaseLine = { + "file": relFile, + "Name": [relFile.stem], + "initializedFrom": "csvfile" + } + initialIndex = np.where(timeDepRelValues["timeStep"] == 0)[0] + releaseLine["rasterData"] = gI.timeDepRelCoordsToRaster(timeDepRelValues, initialIndex, demHeader, + parameter="thickness") + relThFieldData = releaseLine["rasterData"] + # TODO: define thickness for output report, now mean of thickness values + releaseLine["thickness"] = np.nanmean(np.where(relThFieldData == 0, np.nan, relThFieldData)) + + return releaseLine, relThFieldData + + def prepareTimeDepRelLine(releaseLine, cfg): """ read time dependent release values and return them as a dictionary containing: diff --git a/avaframe/in1Data/getInput.py b/avaframe/in1Data/getInput.py index bc1fa8d9d..5911c70ee 100644 --- a/avaframe/in1Data/getInput.py +++ b/avaframe/in1Data/getInput.py @@ -220,6 +220,13 @@ def getInputDataCom1DFA(avaDir): relFiles = sorted( list(releaseDir.glob("*.shp")) + list(releaseDir.glob("*.tif")) + list(releaseDir.glob("*.asc")) ) + timeDepRelFiles = sorted(releaseDir.glob("*.csv")) + # A time-dependent release csv file with x/y columns can be its own release + # geometry. Only use csv files as release scenarios when no conventional + # release geometry was supplied, preserving the established csv file+polygon + # workflow. + if len(relFiles) == 0 and len(timeDepRelFiles) > 0: + relFiles = timeDepRelFiles relSuffixList = [relF.suffix for relF in relFiles] if ".shp" in relSuffixList and (".asc" in relSuffixList or ".tif" in relSuffixList): @@ -227,7 +234,7 @@ def getInputDataCom1DFA(avaDir): log.error(message) raise AssertionError(message) if len(relFiles) == 0: - message = "No release area is found - provide a .shp or .asc or .tif file" + message = "No release area is found - provide a .shp, .asc, .tif, or time-dependent release .csv file" log.error(message) raise FileNotFoundError(message) else: @@ -301,7 +308,6 @@ def getInputDataCom1DFA(avaDir): entResInfo["resRemeshed"] = "No" entResInfo["bhdRemeshed"] = "No" - timeDepRelFiles = sorted(list(releaseDir.glob("*.csv"))) if len(timeDepRelFiles) > 0: entResInfo["timeDepRelCsvAvailable"] = "Yes" else: @@ -526,7 +532,7 @@ def updateThicknessCfg(inputSimFiles, cfgInitial): # update configuration with thickness value to be used for simulations cfgInitial = dP.getThicknessValue(cfgInitial, inputSimFiles, releaseA, "relTh") cfgInitial["INPUT"]["relThFile"] = "" - if inputSimFiles["entResInfo"]["relThFileType"] != ".shp": + if inputSimFiles["entResInfo"]["relThFileType"] in [".asc", ".tif"]: cfgInitial["INPUT"]["relThFile"] = str( pathlib.Path("REL", releaseA + inputSimFiles["entResInfo"]["relThFileType"]) ) @@ -735,7 +741,7 @@ def fetchReleaseFile(inputSimFiles, releaseScenario, cfgSim, releaseList): # update config entry for release scenario, thickness and id cfgSim["INPUT"]["releaseScenario"] = str(releaseScenario) # check if release thickness is read from shapefile or raster file - if releaseScenarioPath.suffix in [".asc", ".tif"]: + if releaseScenarioPath.suffix in [".asc", ".tif", ".csv"]: # raster file - set relThFile path cfgSim["INPUT"]["relThFile"] = str( releaseScenarioPath.parts[-2] + "/" + releaseScenarioPath.parts[-1] @@ -1202,7 +1208,6 @@ def getTimeDepRelCsv(timeDepRelCsv): log.error(message) raise FileNotFoundError(message) timeDepRelDF = pd.read_csv(timeDepRelCsv, index_col=False) - timeDepRelDF = timeDepRelDF.fillna(0.0) # delete empty spaces and write column names in low case timeDepRelDF.columns = timeDepRelDF.columns.str.strip().str.lower() @@ -1211,17 +1216,90 @@ def getTimeDepRelCsv(timeDepRelCsv): timeDepRelValues = { "timeStep": timeDepRelDF["timestep"].to_numpy(dtype=np.float64), "thickness": timeDepRelDF["thickness"].to_numpy(dtype=np.float64), - "velocity": timeDepRelDF["velocity"].to_numpy(dtype=np.float64), } + # TODO: rethink this actual status: we only allow velocity magnitude together with shape file location and only velocity components with x and y (from csv) locations + for component in ["velocityx", "velocityy", "velocityz"]: + if "x" not in timeDepRelDF.columns and component in timeDepRelDF.columns: + message = "If release location is defined by shape file, only the velocity magnitude can be provided in the csv file." + log.error(message) + raise ValueError(message) + if "x" in timeDepRelDF.columns and "velocity" in timeDepRelDF.columns: + message = "If release location is defined by coordinates in csv file, only the velocity components can be provided." + log.error(message) + raise ValueError(message) + + if "velocity" in timeDepRelDF.columns: + timeDepRelDF["velocity"] = timeDepRelDF["velocity"].fillna(0.0) + timeDepRelValues["velocity"] = timeDepRelDF["velocity"].to_numpy(dtype=np.float64) + else: + for component in ["velocityx", "velocityy", "velocityz"]: + print() + if component not in timeDepRelDF.columns: + message = "Please provide x, y, and z velocity component in the time dependent release values (csv file)." + log.error(message) + raise ValueError(message) + timeDepRelValues["velocityX"] = timeDepRelDF["velocityx"].to_numpy(dtype=np.float64) + timeDepRelValues["velocityY"] = timeDepRelDF["velocityy"].to_numpy(dtype=np.float64) + timeDepRelValues["velocityZ"] = timeDepRelDF["velocityz"].to_numpy(dtype=np.float64) + + if "x" in timeDepRelDF.columns: + timeDepRelValues["x"] = timeDepRelDF["x"].to_numpy(dtype=np.float64) + if "y" in timeDepRelDF.columns: + timeDepRelValues["y"] = timeDepRelDF["y"].to_numpy(dtype=np.float64) # check if some criterias are satisfied in the csv file checkTimeDepRelease(timeDepRelValues, timeDepRelCsv) return timeDepRelValues, timeDepRelDF +def timeDepRelCoordsToRaster(timeDepRelValues, index, demHeader, parameter="thickness"): + """Create a one-cell release raster at the csv file location for one timestep. + + The x and y coordintates are map coordinates of the DEM cell centres. Coordinates + between centres are assigned to their nearest raster cell. + + Parameters + ---------- + timeDepRelValues: dict + contains time dependent release values: timestep, thickness, velocity, x and y coordinates + index: int + index of timestep (row in csv file) + demHeader: dict + header of DEM file + parameter: str + parameter of timeDepRelValues that is written into the raster + + Returns + -------- + raster: np.array + thickness raster read from time dependent csv file (coordinates and thickness) + """ + if "x" not in timeDepRelValues: + return None + + cellsize = demHeader["cellsize"] + cols = np.rint((timeDepRelValues["x"][index] - demHeader["xllcenter"]) / cellsize).astype(int) + rows = np.rint((timeDepRelValues["y"][index] - demHeader["yllcenter"]) / cellsize).astype(int) + if not np.all( + (cols >= 0) & (cols < demHeader["ncols"]) & + (rows >= 0) & (rows < demHeader["nrows"]) + ): + message = ( + "A time dependent release location at timestep %.3f s is outside the DEM extent" + % (timeDepRelValues["timeStep"][index]) + ) + log.error(message) + raise ValueError(message) + + raster = np.zeros((demHeader["nrows"], demHeader["ncols"]), dtype=float) + np.add.at(raster, (rows, cols), timeDepRelValues[parameter][index]) + return raster + + def checkTimeDepRelease(timeDepRelValues, timeDepRelCsv): """ check if time dependent release values satisfy the following requirements: - - release - timesteps are unique + - if coordinates are provided: x and y coordinates (not just one) are provided in every row + - release - timesteps are unique (when no coordinates are provided) - the release - timesteps are not too close (that the particle density becomes too high) - provided release - thickness is larger than zero - provided velocity is zero or larger. @@ -1233,17 +1311,50 @@ def checkTimeDepRelease(timeDepRelValues, timeDepRelCsv): timeDepRelValues: dict contains time dependent release values: timestep, thickness, velocity """ + + xInCsv = "x" in timeDepRelValues + yInCsv = "y" in timeDepRelValues + if xInCsv != yInCsv: + message = "Time dependent release csv file %s must provide both x and y columns" % timeDepRelCsv + log.error(message) + raise ValueError(message) + if xInCsv and (np.any(np.isnan(timeDepRelValues["x"])) or np.any(np.isnan(timeDepRelValues["y"]))): + message = "Time dependent release csv file %s must provide x and y values for every timestep" % timeDepRelCsv + log.error(message) + raise ValueError(message) + + # if x,y coordinates are provided, check that for each coordinate pair + # there is only one row per timestep (no duplicate timesteps at same location) + if xInCsv: + xVals = np.asarray(timeDepRelValues["x"]) + yVals = np.asarray(timeDepRelValues["y"]) + tVals = np.asarray(timeDepRelValues["timeStep"]) + seenDict = {} + for xVal, yVal, tVal in zip(xVals, yVals, tVals): + key = (xVal, yVal) + if key not in seenDict: + seenDict[key] = set() + if tVal in seenDict[key]: + message = ( + "For coordinate (x=%s, y=%s) in %s, only one timestep is allowed" + % (xVal, yVal, timeDepRelCsv) + ) + log.error(message) + raise ValueError(message) + seenDict[key].add(tVal) + # check if timesteps are unique timeStepUnique = np.unique(timeDepRelValues["timeStep"]) - if timeStepUnique.ndim == 0: - if timeStepUnique != timeDepRelValues["timeStep"]: - message = "The provided time dependent release time steps in %s are not unique" % (timeDepRelCsv) + if not xInCsv: + if timeStepUnique.ndim == 0: + if timeStepUnique != timeDepRelValues["timeStep"]: + message = "The provided time dependent release time steps in %s are not unique" % (timeDepRelCsv) + log.error(message) + raise ValueError(message) + elif len(timeStepUnique) != len(timeDepRelValues["timeStep"]): + message = "The provided time dependent release timesteps in %s are not unique" % (timeDepRelCsv) log.error(message) raise ValueError(message) - elif len(timeStepUnique) != len(timeDepRelValues["timeStep"]): - message = "The provided time dependent release timesteps in %s are not unique" % (timeDepRelCsv) - log.error(message) - raise ValueError(message) # check if a timestep = 0 is provided if 0 not in timeStepUnique: @@ -1263,11 +1374,12 @@ def checkTimeDepRelease(timeDepRelValues, timeDepRelCsv): log.error(message) raise ValueError(message) - for vel in timeDepRelValues["velocity"]: - if vel < 0: - message = "The initial velocity provided in %s can not be negative." % (timeDepRelCsv) - log.error(message) - raise ValueError(message) + if "velocity" in timeDepRelValues.keys(): + for vel in timeDepRelValues["velocity"]: + if vel < 0: + message = "The initial velocity provided in %s can not be negative." % (timeDepRelCsv) + log.error(message) + raise ValueError(message) def preprocessAssets(avalancheDir, dem, cfg): diff --git a/avaframe/out3Plot/outCom1DFA.py b/avaframe/out3Plot/outCom1DFA.py index a142b3e9e..c10ca2efa 100644 --- a/avaframe/out3Plot/outCom1DFA.py +++ b/avaframe/out3Plot/outCom1DFA.py @@ -571,6 +571,16 @@ def plotReleaseScenarioView( relArea.plot(ax=ax, edgecolor="darkblue", linewidth=2, facecolor="none") relPatch = Patch(color="darkblue", label="release") handles.append(relPatch) + else: + relAreaPlot = np.where(releaseLine["rasterData"] > 0, releaseLine["rasterData"], np.nan) + ax.imshow( + relAreaPlot, + extent=extentCells, + cmap="Blues", + vmin=0, + vmax=1, + zorder=1000, + ) count = 1 if reportAreaInfo["resistance"] == "Yes": if inputSimLines["resLine"]["initializedFrom"] == "shapefile": diff --git a/avaframe/tests/test_com1DFA.py b/avaframe/tests/test_com1DFA.py index 83b8aadcc..77a041dff 100644 --- a/avaframe/tests/test_com1DFA.py +++ b/avaframe/tests/test_com1DFA.py @@ -403,6 +403,111 @@ def test_prepareInputData(tmp_path): assert inputSimLines["entLine"] is None +def test_prepareInputDataCoordinateCsvRelease(tmp_path): + """A coordinate time-dependent CSV supplies both release geometry and thickness.""" + sourceInputs = pathlib.Path(__file__).parents[0] / ".." / "data" / "avaParabola" / "Inputs" + avaDir = tmp_path / "avaCoordinateRelease" + shutil.copytree(sourceInputs, avaDir / "Inputs") + demHeader = getInput.initializeDEM(avaDir, demPath="DEM_PF_Topo.asc")["header"] + coordinateCsv = avaDir / "Inputs" / "REL" / "coordinateRelease.csv" + coordinateCsv.write_text( + "timestep,thickness,velocityX,velocityY,velocityZ,x,y\n" + "%s,1.5,3.0,-4.0,5.0,%s,%s\n" + % (0.0, demHeader["xllcenter"], demHeader["yllcenter"]) + ) + inputSimFiles = { + "releaseScenario": coordinateCsv, + "entResInfo": {"flagEnt": "No", "flagRes": "No", "flagSecondaryRelease": "No"}, + "muFile": None, + "xiFile": None, + "kFile": None, + "tauCFile": None, + } + cfg = configparser.ConfigParser() + cfg["GENERAL"] = { + "avalancheDir": str(avaDir), + "timeDependentRelease": "True", + "secRelArea": "False", + "simTypeActual": "null", + "dam": "False", + "relThFromFile": "False", + } + cfg["INPUT"] = { + "DEM": "DEM_PF_Topo.asc", + "releaseScenario": coordinateCsv.stem, + "relThFile": "", + "timeDepRelCsv": str(coordinateCsv), + } + + demOri, inputSimLines = com1DFA.prepareInputData(inputSimFiles, cfg) + + releaseLine = inputSimLines["releaseLine"] + assert releaseLine["initializedFrom"] == "csvfile" + assert releaseLine["Name"] == [coordinateCsv.stem] + assert releaseLine["thickness"] == 1.5 + assert releaseLine["rasterData"].shape == (demOri["header"]["nrows"], demOri["header"]["ncols"]) + assert releaseLine["rasterData"][0, 0] == 1.5 + assert inputSimLines["relThField"][0, 0] == 1.5 + assert np.array_equal(releaseLine["timeDepRelValues"]["x"], np.array([demHeader["xllcenter"]])) + assert np.array_equal(releaseLine["timeDepRelValues"]["y"], np.array([demHeader["yllcenter"]])) + assert np.array_equal(releaseLine["timeDepRelValues"]["velocityX"], np.array([3.0])) + assert np.array_equal(releaseLine["timeDepRelValues"]["velocityY"], np.array([-4.0])) + assert np.array_equal(releaseLine["timeDepRelValues"]["velocityZ"], np.array([5.0])) + + # All coordinate CSV rows are retained and timestep-zero rows form the initial raster. + sourceInputs = pathlib.Path(__file__).parents[0] / ".." / "data" / "avaParabola" / "Inputs" + avaDir = tmp_path / "avaCoordinateReleaseMultipleRows" + shutil.copytree(sourceInputs, avaDir / "Inputs") + demHeader = getInput.initializeDEM(avaDir, demPath="DEM_PF_Topo.asc")["header"] + coordinateCsv = avaDir / "Inputs" / "REL" / "coordinateRelease.csv" + x0 = demHeader["xllcenter"] + y0 = demHeader["yllcenter"] + cellsize = demHeader["cellsize"] + coordinateCsv.write_text( + "timestep,thickness,velocityX,velocityY,velocityZ,x,y\n" + "%s,1.5,3.0,-4.0,5.0,%s,%s\n" + "%s,2.5,-1.0,2.0,0.0,%s,%s\n" + "%s,3.0,0.0,0.0,-2.0,%s,%s\n" + % (0.0, x0, y0, 0.0, x0 + cellsize, y0 + cellsize, 10.0, x0 + 2 * cellsize, y0 + 2 * cellsize) + ) + inputSimFiles = { + "releaseScenario": coordinateCsv, + "entResInfo": {"flagEnt": "No", "flagRes": "No", "flagSecondaryRelease": "No"}, + "muFile": None, + "xiFile": None, + "kFile": None, + "tauCFile": None, + } + cfg = configparser.ConfigParser() + cfg["GENERAL"] = { + "avalancheDir": str(avaDir), + "timeDependentRelease": "True", + "secRelArea": "False", + "simTypeActual": "null", + "dam": "False", + "relThFromFile": "False", + } + cfg["INPUT"] = { + "DEM": "DEM_PF_Topo.asc", + "releaseScenario": coordinateCsv.stem, + "relThFile": "", + "timeDepRelCsv": str(coordinateCsv), + } + + _, inputSimLines = com1DFA.prepareInputData(inputSimFiles, cfg) + + releaseLine = inputSimLines["releaseLine"] + assert releaseLine["initializedFrom"] == "csvfile" + assert releaseLine["thickness"] == 2.0 + assert releaseLine["rasterData"][0, 0] == 1.5 + assert releaseLine["rasterData"][1, 1] == 2.5 + assert np.sum(releaseLine["rasterData"]) == 4.0 + assert np.array_equal(releaseLine["timeDepRelValues"]["timeStep"], np.array([0.0, 0.0, 10.0])) + assert np.array_equal(releaseLine["timeDepRelValues"]["thickness"], np.array([1.5, 2.5, 3.0])) + assert np.array_equal(releaseLine["timeDepRelValues"]["x"], np.array([x0, x0 + cellsize, x0 + 2 * cellsize])) + assert np.array_equal(releaseLine["timeDepRelValues"]["y"], np.array([y0, y0 + cellsize, y0 + 2 * cellsize])) + + def test_prepareReleaseEntrainment(tmp_path): """test preparing release areas""" @@ -1523,6 +1628,7 @@ def test_initializeParticles(): "thickness": [1.0], "rasterData": relRaster, "type": "Release", + "initializedFrom": "shapefile", } releaseLine["header"] = demHeader diff --git a/avaframe/tests/test_debrisFunctions.py b/avaframe/tests/test_debrisFunctions.py index 92c01f8a0..696ddc1c6 100644 --- a/avaframe/tests/test_debrisFunctions.py +++ b/avaframe/tests/test_debrisFunctions.py @@ -14,6 +14,7 @@ def test_addReleaseParticles(): "Start": np.asarray([0.0]), "Length": np.asarray([5]), "type": "time dependent Release", + "initializedFrom": "shapefile", "x": np.asarray( [ 0, @@ -29,8 +30,7 @@ def test_addReleaseParticles(): "thickness": 1, } } - thickness = inputSimLines["releaseLine"]["thickness"] - velocityMag = 0 + timeDepRelValues = {"thickness": np.array([1.0]), "velocity": np.array([0.0])} demHeader = {} demHeader["xllcenter"] = 0 @@ -119,7 +119,7 @@ def test_addReleaseParticles(): zPartArray0Test = np.ones(particlesTest["nPart"]) particlesNewRel, zPartArray0NewRel = debF.addReleaseParticles( - cfg, particles, inputSimLines, thickness, velocityMag, dem, zPartArray0 + cfg, particles, inputSimLines, timeDepRelValues, dem, zPartArray0, timeDepRelIndex=0 ) assert np.all(np.equal(zPartArray0NewRel, zPartArray0Test)) @@ -136,7 +136,7 @@ def test_addReleaseParticles(): cfg["GENERAL"]["thresholdMassSplit"] = "1.5" particlesNewRel, zPartArray0NewRel = debF.addReleaseParticles( - cfg, particles, inputSimLines, thickness, velocityMag, dem, zPartArray0 + cfg, particles, inputSimLines, timeDepRelValues, dem, zPartArray0, timeDepRelIndex=0 ) assert particlesNewRel["nPart"] == 16 + 3 for key in ["ux", "uy", "uz", "velocityMag", "x", "y", "z"]: @@ -147,11 +147,111 @@ def test_addReleaseParticles(): particles["y"] = np.array([5, 3, 30]) with pytest.raises(ValueError): - debF.addReleaseParticles(cfg, particles, inputSimLines, thickness, velocityMag, dem, zPartArray0) + debF.addReleaseParticles(cfg, particles, inputSimLines, timeDepRelValues, dem, zPartArray0, timeDepRelIndex=0) + + +def test_addReleaseParticlesFromCoordinateCsv(): + """Coordinate CSV releases initialize particles and project vector velocity.""" + header = {"xllcenter": 0.0, "yllcenter": 0.0, "cellsize": 5.0, "nrows": 7, "ncols": 7} + dem = { + "header": header, + "originalHeader": header, + "rasterData": np.ones((7, 7)), + "areaRaster": np.ones((7, 7)), + "Nx": np.zeros((7, 7)), + "Ny": np.zeros((7, 7)), + "Nz": np.zeros((7, 7)), + } + cfg = configparser.ConfigParser() + cfg["GENERAL"] = { + "resType": "ppr|pft|pfv", + "rho": "1000.", + "gravAcc": "9.81", + "cpIce": "2050", + "TIni": "-10", + "avalancheDir": "data/avaParabola", + "massPerParticleDeterminationMethod": "MPPDIR", + "interpOption": "2", + "initialiseParticlesFromFile": "False", + "iniStep": "False", + "seed": "12345", + "sphKernelRadius": "1", + "deltaTh": "1", + "initPartDistType": "uniform", + "thresholdPointInPoly": "0.001", + "massPerPart": "100000", + "thresholdPointInRel": "0", + } + timeDepRelValues = { + "timeStep": np.array([5.0, 5.0]), + "thickness": np.array([1.5, 2.5]), + "velocityX": np.array([3.0, -1.0]), + "velocityY": np.array([-4.0, 2.0]), + "velocityZ": np.array([5.0, 0.0]), + "x": np.array([0.0, 5.0]), + "y": np.array([0.0, 5.0]), + } + inputSimLines = { + "releaseLine": { + "initializedFrom": "csvfile", + "timeDepRelValues": timeDepRelValues, + "rasterData": np.zeros((3, 3)), + } + } + particles = { + "nPart": 3, + "x": np.array([12.0, 20.0, 30.0]), + "y": np.array([5.0, 10.0, 30.0]), + "z": np.array([1.0, 1.0, 1.0]), + "m": np.array([1000.0, 1000.0, 1000.0]), + "idFixed": np.zeros(3), + "t": 4.0, + "dt": 0.5, + "massPerPart": 1000.0, + "mTot": 3000.0, + "totalEnthalpy": np.full(3, -20490.19), + "tPlot": 0, + "h": np.ones(3), + "ux": np.zeros(3), + "uy": np.zeros(3), + "uz": np.zeros(3), + "uAcc": np.zeros(3), + "velocityMag": np.zeros(3), + "trajectoryLengthXY": np.zeros(3), + "trajectoryLengthXYCor": np.zeros(3), + "trajectoryLengthXYZ": np.zeros(3), + "trajectoryAngle": np.zeros(3), + "stoppCriteria": False, + "peakForceSPH": 0.0, + "forceSPHIni": 0.0, + "peakMassFlowing": 0, + "xllcenter": 0.0, + "yllcenter": 0.0, + "nExitedParticles": 0.0, + "dmDet": np.zeros(3), + "dmEnt": np.zeros(3), + } + particles, zPartArray0 = debF.addReleaseParticles( + cfg, + particles, + inputSimLines, + timeDepRelValues, + dem, + np.ones(3), + timeDepRelIndex=np.array([0, 1]), + ) + + assert particles["nPart"] == 5 + assert np.array_equal(particles["x"], np.array([12.0, 20.0, 30.0, 0.0, 5.0])) + assert np.array_equal(particles["y"], np.array([5.0, 10.0, 30.0, 0.0, 5.0])) + assert np.array_equal(particles["ux"], np.array([0.0, 0.0, 0.0, 3.0, -1.0])) + assert np.array_equal(particles["uy"], np.array([0.0, 0.0, 0.0, -4.0, 2.0])) + assert np.array_equal(particles["uz"], np.array([0.0, 0.0, 0.0, 5.0, 0.0])) + assert np.array_equal(zPartArray0, np.ones(5)) """ -Test does not word because: When calling pytest, executing DFAfunctionsCython.upfateFieldsC() raises an error ("Fatal Python error: Aborted") +Test does not work because: When calling pytest, executing DFAfunctionsCython.upfateFieldsC() raises an error ("Fatal Python error: Aborted") (see issue #1002?) ------------------------------ diff --git a/avaframe/tests/test_getInput.py b/avaframe/tests/test_getInput.py index 2f504c3c8..7a3acb989 100644 --- a/avaframe/tests/test_getInput.py +++ b/avaframe/tests/test_getInput.py @@ -203,6 +203,28 @@ def test_getInputDataCom1DFA(tmp_path): assert inputSimFiles["relThFile"] == None +def test_getInputDataCom1DFAUsesCoordinateCsvWithoutReleaseGeometry(tmp_path): + """A coordinate time-dependent release CSV is a release scenario on its own.""" + dirPath = pathlib.Path(__file__).parents[0] + avaDir = tmp_path / "avaCoordinateRelease" + shutil.copytree(dirPath / ".." / "data" / "avaHockeyChannel" / "Inputs", avaDir / "Inputs") + releaseDir = avaDir / "Inputs" / "REL" + for releaseFile in releaseDir.glob("*.shp"): + releaseFile.unlink() + coordinateRelease = releaseDir / "coordinateRelease.csv" + coordinateRelease.write_text( + "timestep,thickness,velocityX,velocityY,velocityZ,x,y\n" + "0,1,0,0,0,1000,-5000\n" + ) + + inputSimFiles = getInput.getInputDataCom1DFA(avaDir) + + assert inputSimFiles["relFiles"] == [coordinateRelease] + assert inputSimFiles["timeDepRelCsv"] == [coordinateRelease] + assert inputSimFiles["entResInfo"]["relThFileType"] == ".csv" + assert inputSimFiles["entResInfo"]["timeDepRelCsvAvailable"] == "Yes" + + def test_getAndCheckInputFiles(tmp_path): """test fetching input files and checking if exist""" @@ -1583,10 +1605,6 @@ def test_getTimeDepRelCsv(): assert np.all(timeDepRelValues["velocity"] == np.array([5, 3, 0])) testDir = pathlib.Path(__file__).parents[0] - timeDepRelCsv = testDir / "data" / "testTimeDepRel" / "rel.csv" - - with pytest.raises(ValueError): - timeDepRelValues, timeDepRelValuesTxt = getInput.getTimeDepRelCsv(timeDepRelCsv) timeDepRelCsv = testDir / "data" / "testTimeDepRel" / "rel_notSorted.csv" timeDepRelValues, timeDepRelValuesTxt = getInput.getTimeDepRelCsv(timeDepRelCsv) @@ -1594,6 +1612,67 @@ def test_getTimeDepRelCsv(): assert np.all(timeDepRelValues["thickness"] == np.array([3, 1, 1])) +def test_getTimeDepRelCsvWithCoordinatesAndVectorVelocity(tmp_path): + """Coordinate releases may share a timestep and retain velocity components.""" + timeDepRelCsv = tmp_path / "coordinateRelease.csv" + timeDepRelCsv.write_text( + "timestep, thickness, velocityX, velocityY, velocityZ, x, y\n" + "0, 1.5, 3.0, -4.0, 5.0, 101.0, 202.0\n" + "0, 2.5, -1.0, 2.0, 0.0, 119.0, 219.0\n" + "10, 3.0, 0.0, 0.0, -2.0, 110.0, 210.0\n" + ) + + values, table = getInput.getTimeDepRelCsv(timeDepRelCsv) + + assert table.shape == (3, 7) + assert np.array_equal(values["timeStep"], np.array([0.0, 0.0, 10.0])) + assert np.array_equal(values["velocityX"], np.array([3.0, -1.0, 0.0])) + assert np.array_equal(values["velocityY"], np.array([-4.0, 2.0, 0.0])) + assert np.array_equal(values["velocityZ"], np.array([5.0, 0.0, -2.0])) + + +def test_timeDepRelCoordsToRaster(): + """Coordinate release values are assigned to their nearest DEM cells.""" + values = { + "timeStep": np.array([0.0, 0.0, 10.0]), + "thickness": np.array([1.5, 2.5, 3.0]), + "velocityX": np.array([3.0, -1.0, 0.0]), + "x": np.array([101.0, 119.0, 110.0]), + "y": np.array([202.0, 219.0, 210.0]), + } + header = {"xllcenter": 100.0, "yllcenter": 200.0, "cellsize": 10.0, "nrows": 3, "ncols": 3} + thicknessRaster = getInput.timeDepRelCoordsToRaster(values, np.array([0, 1]), header) + velocityXRaster = getInput.timeDepRelCoordsToRaster(values, 1, header, parameter="velocityX") + + assert np.array_equal(thicknessRaster, np.array([[1.5, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 2.5]])) + assert np.array_equal(velocityXRaster, np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, -1.0]])) + assert getInput.timeDepRelCoordsToRaster({"thickness": np.array([1.0])}, 0, header) is None + + values["x"][2] = 1000.0 + with pytest.raises(ValueError, match="outside the DEM extent"): + getInput.timeDepRelCoordsToRaster(values, 2, header) + + +def test_checkTimeDepReleaseCoordinateValidation(): + """Coordinate releases require complete coordinates and unique events per location.""" + csvPath = "coordinateRelease.csv" + values = { + "timeStep": np.array([0.0, 10.0]), + "thickness": np.array([1.0, 1.0]), + "velocityX": np.array([0.0, 0.0]), + "velocityY": np.array([0.0, 0.0]), + "velocityZ": np.array([0.0, 0.0]), + "x": np.array([100.0, 100.0]), + } + with pytest.raises(ValueError, match="both x and y columns"): + getInput.checkTimeDepRelease(values, csvPath) + + values["y"] = np.array([200.0, 200.0]) + values["timeStep"] = np.array([0.0, 0.0]) + with pytest.raises(ValueError, match="only one timestep is allowed"): + getInput.checkTimeDepRelease(values, csvPath) + + def test_checkTimeDepRelease(): timeDepRelValues = { "timeStep": np.array([0, 10, 20, 30, 35]), diff --git a/docs/moduleCom1DFA.rst b/docs/moduleCom1DFA.rst index f7b7414be..27487245a 100644 --- a/docs/moduleCom1DFA.rst +++ b/docs/moduleCom1DFA.rst @@ -179,17 +179,21 @@ input file (shape file or raster file) or 2) through the :py:mod:`com1DFA` confi - if the flag `timeDependentRelease` is True, in various provided time steps flowing mass is initialized (`relThFromFile` is also set to True, currently the only option to read time dependent thickness is from csv file) - - additional to a .shp file (raster file does not work yet), at least one csv file is provided in the `REL` folder, that contains: + - provide at least one csv file in the `REL` folder. A release polygon can still be used; alternatively, + include location columns in the CSV and no polygon is required. - a header (first line) - the following columns with the respective column names: - timestep values (column name: "timestep") - thickness values (column name: "thickness") - - initial velocity values (column name: "velocity") + - initial velocity values (column name: "velocity") or velocity components + (column names: "velocityX", "velocityY", "velocityZ") + - optional release location values (column names: "x" and "y"). Both columns must be provided + together; each row initializes particles in the DEM grid cell nearest to that coordinate. - the delimiter is: "," - In each provided timestep, particles are initialized with the provided corresponding thickness - in the release area. If a velocity is provided, the particles have this initial velocity + in the release area, or in the csv file defined grid cell when ``x``/``y`` are present. If a velocity components are provided, the particles have this initial velocity in the direction of the steepest descent.