diff --git a/avaframe/com7Regional/__init__.py b/avaframe/com7Regional/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/avaframe/com7Regional/com7Regional.py b/avaframe/com7Regional/com7Regional.py new file mode 100644 index 000000000..90f8f393e --- /dev/null +++ b/avaframe/com7Regional/com7Regional.py @@ -0,0 +1,457 @@ +"""Module for handling regional avalanche simulations.""" + +import pathlib +import shutil +import logging +import numpy as np +from concurrent.futures import ProcessPoolExecutor, as_completed + +import avaframe.in3Utils.initializeProject as initProj +from avaframe.com1DFA import com1DFA +from avaframe.in3Utils import cfgUtils, cfgHandling +from avaframe.in3Utils import logUtils +from avaframe.in2Trans import rasterUtils +from avaframe.in3Utils import fileHandlerUtils as fU + +from rasterio.merge import merge + +# create local logger +log = logging.getLogger(__name__) + + +def com7RegionalMain(cfgMain, cfg): + """Run com7Regional with given configuration. + + This function processes multiple avalanche directories in parallel, running simulations + for each directory. + + Parameters + ---------- + cfgMain : configparser.ConfigParser + Main avaframe configuration settings + cfg : configparser.ConfigParser + Regional configuration settings with potential overrides + + Returns + ------- + allPeakFilesDir : pathlib.Path or None + Path to the directory containing all peak files, if copyPeakFiles is True + mergedRastersDir : pathlib.Path or None + Path to the directory containing merged rasters, if mergeOutput is True + + Notes + ----- + The function expects the following directory structure: + avalancheDir/ + └── regionalDir/ + ├── avalanche1/ + ├── avalanche2/ + └── ... + + Where: + - avalancheDir: Main directory specified in cfgMain + - regionalDir: Subdirectory specified in cfg['GENERAL']['regionalDir'] + """ + # Define the regional directory in relation to the avalanche directory + regionalDirFromCfg = str(cfg["GENERAL"]["regionalDir"]) + regionalDir = pathlib.Path(cfgMain["MAIN"]["avalancheDir"]) / regionalDirFromCfg + + # List valid avalanche directories within the regional directory + avaDirs = findAvaDirs(regionalDir) + + # Get total number of simulations + log.info(f"Getting total number of simulations to perform...") + with logUtils.silentLogger(): + totalSims = getTotalNumberOfSims(avaDirs, cfgMain, cfg) + log.info(f"Found {totalSims} (new) simulations to perform across {len(avaDirs)} directories") + + # Get number of processes based on number of avaDirs + nProcesses = cfgUtils.getNumberOfProcesses(cfgMain, len(avaDirs)) + + # Set nCPU for com1 to 1 to avoid nested parallelization + cfgMain["MAIN"]["nCPU"] = "1" + + # Track progress and results + completed = 0 + nSuccesses = 0 + + # Process avalanche directories within the regional folder concurrently + with ProcessPoolExecutor(max_workers=nProcesses) as executor: + # Submit each avalanche directory to the executor + futures = { + executor.submit(processAvaDirCom1Regional, cfgMain, cfg, avaDir): avaDir for avaDir in avaDirs + } + # Log results as each future completes + for future in as_completed(futures): + avaDir = futures[future] + try: + resultDir, status = future.result() + completed += 1 + + if status == "Success": + nSuccesses += 1 + + log.info( + f"{status} in directory: {resultDir.relative_to(pathlib.Path(regionalDir))} " + f"- Overall progress: {completed}/{len(avaDirs)}" + ) + except Exception as e: + log.error(f"Error processing {avaDir}: {e}") + + log.info(f"Processing complete. Success in {nSuccesses} out of {len(avaDirs)} directories.") + + # Copy/move peak files if configured + allPeakFilesDir = None + if cfg["GENERAL"].getboolean("copyPeakFiles"): + allPeakFilesDir = moveOrCopyPeakFiles(cfg, regionalDir) + + # Merge output rasters if configured + mergedRastersDir = None + if cfg["GENERAL"].getboolean("mergeOutput"): + mergedRastersDir = mergeOutputRasters(cfg, regionalDir) + + return allPeakFilesDir, mergedRastersDir + + +def findAvaDirs(Dir): + """Find all valid avalanche directories within a given directory. + + A directory is considered a valid avalanche directory if it contains an "Inputs" folder. + + Parameters + ---------- + Dir : pathlib.Path or str + Path to the directory to search in + + Returns + ------- + avaDirs : list + List of pathlib.Path objects pointing to valid avalanche directories + """ + avaDirs = [pathlib.Path(p).parent for p in pathlib.Path(Dir).glob("*/Inputs")] + log.info(f"Found a total of '{len(avaDirs)}' avalanche directories in: {Dir}:") + for avaDir in avaDirs: + log.info(f"'{avaDir.name}'") + + return avaDirs + + +def getTotalNumberOfSims(avaDirs, cfgMain, cfgCom7): + """Get total number of simulations across all avalanche directories. + + Parameters + ---------- + avaDirs : list + List of avalanche directories + cfgMain : configparser.ConfigParser + Main configuration + cfgCom7 : configparser.ConfigParser + Regional configuration with potential overrides + + Returns + ------- + int + Total number of simulations + """ + totalSims = 0 + for avaDir in avaDirs: + # Create copies of configs to avoid modifying originals + cfgMainCopy = cfgUtils.convertDictToConfigParser(cfgUtils.convertConfigParserToDict(cfgMain)) + cfgMainCopy["MAIN"]["avalancheDir"] = str(avaDir) + + # Get com1DFA config with regional overrides (same as in processAvaDirCom1Regional) + cfgCom1DFA = cfgUtils.getModuleConfig( + com1DFA, + fileOverride="", + toPrint=False, + onlyDefault=cfgCom7["com1DFA_com1DFA_override"].getboolean("defaultConfig"), + ) + cfgCom1DFA, _ = cfgHandling.applyCfgOverride(cfgCom1DFA, cfgCom7, com1DFA, addModValues=False) + + # Get simulations for this directory + try: + simDict, _, _, _ = com1DFA.com1DFAPreprocess(cfgMainCopy, "cfgFromObject", cfgCom1DFA) + totalSims += len(simDict) + except Exception as e: + log.warning(f"Could not get simulations for {avaDir}: {e}") + continue + + return totalSims + + +def processAvaDirCom1Regional(cfgMain, cfgCom7, avalancheDir): + """Run com1DFA simulation in a specific avalanche directory with regional settings. + + Note: This function calls com1DFA within each avalanche directory within the input directory. + If wanted it may be used as a template to call another operation within each directory, such as com2AB, ana5Utils, etc. + + Parameters + ---------- + cfgMain : configparser.ConfigParser + Main configuration settings + cfgCom7 : configparser.ConfigParser + Regional configuration settings with potential overrides + avalancheDir : pathlib.Path or str + Path to the avalanche directory to process + + Returns + ------- + avalancheDir : pathlib.Path or str + Path to the avalanche directory that was processed + status : str + Status of the simulation, "Success" if completed + """ + # Initialize log for each process + log = logUtils.initiateLogger(avalancheDir, logName="runCom1DFA") + log.info("COM1DFA PROCESS CALLED BY COM7REGIONAL RUN") + log.info("Current avalanche: %s", avalancheDir) + + # Update cfgMain setting to reflect the current avalancheDir + cfgMain["MAIN"]["avalancheDir"] = str(avalancheDir) + + # Clean input directory of old work and output files from module + initProj.cleanModuleFiles(avalancheDir, com1DFA, deleteOutput=True) + + # Create com1DFA configuration for the current avalanche directory and override with regional settings + cfgCom1DFA = cfgUtils.getModuleConfig( + com1DFA, + fileOverride="", + toPrint=False, + onlyDefault=cfgCom7["com1DFA_com1DFA_override"].getboolean("defaultConfig"), + ) + cfgCom1DFA, cfgCom7 = cfgHandling.applyCfgOverride(cfgCom1DFA, cfgCom7, com1DFA, addModValues=False) + + # Run com1DFA in the current avalanche directory + com1DFA.com1DFAMain(cfgMain, cfgInfo=cfgCom1DFA) + + return avalancheDir, "Success" + + +def moveOrCopyPeakFiles(cfg, avalancheDir): + """Collects peak files from multiple sub-avalanche directories. + + Creates directory allPeakFiles: Contains peak files from all avalanche directories + + Parameters + ---------- + cfg : configparser.ConfigParser + Configuration containing GENERAL settings: + - copyPeakFiles: If True, copy/move files; if False, do nothing + - moveInsteadOfCopy: If True, move files instead of copying + avalancheDir : pathlib.Path or str + Base directory where allPeakFiles will be created + + Returns + ------- + allPeakFilesDir : pathlib.Path or None + Path to the created allPeakFiles directory or None if copyPeakFiles is False + """ + if not cfg["GENERAL"].getboolean("copyPeakFiles"): + log.info("copyPeakFiles is False - no files will be copied or moved") + return None, None + + # Get avalanche directories + # with logUtils.silentLogger(): + avaDirs = findAvaDirs(avalancheDir) + if not avaDirs: + log.warning("No avalanche directories found to copy/move files from") + return None, None + + # Set up outdirs + allPeakFilesDir = pathlib.Path(avalancheDir, "allPeakFiles") + for dirPath in [allPeakFilesDir]: + if dirPath.exists(): + shutil.rmtree(str(dirPath)) + fU.makeADir(dirPath) + + # Get operation type + operation = shutil.move if cfg["GENERAL"].getboolean("moveInsteadOfCopy") else shutil.copy2 + operationType = "Moving" if operation == shutil.move else "Copying" + + # Process each avalanche directory + for avaDir in avaDirs: + log.info(f"{operationType} files from: {avaDir.name}") + + # Process peak files + peakFiles = list(avaDir.glob("Outputs/**/peakFiles/*.*")) + for peakFile in peakFiles: + targetPath = allPeakFilesDir / peakFile.name + operation(str(peakFile), str(targetPath)) + + return allPeakFilesDir + + +def getRasterBounds(rasterFiles): + """Get the union bounds and validate cell sizes of multiple rasters. + + Parameters + ---------- + rasterFiles : list + List of paths to raster files + + Returns + ------- + bounds : dict + Dictionary containing xMin, yMin, xMax, yMax of the union bounds + cellSize : float + Cell size of the rasters + + Raises + ------ + ValueError + If cell sizes of rasters differ + """ + # Read first raster to get cellSize and initialize bounds + firstRaster = rasterUtils.readRaster(rasterFiles[0]) + cellSize = float(firstRaster["header"]["cellsize"]) + bounds = { + "xMin": float("inf"), + "yMin": float("inf"), + "xMax": float("-inf"), + "yMax": float("-inf"), + } + + # Find bounds and validate cell sizes + for rasterFile in rasterFiles: + raster = rasterUtils.readRaster(rasterFile) + header = raster["header"] + + if float(header["cellsize"]) != cellSize: + raise ValueError(f"Different cell sizes found: {cellSize} vs {header['cellsize']}") + + # Update bounds + bounds["xMin"] = min(bounds["xMin"], float(header["xllcenter"])) + bounds["yMin"] = min(bounds["yMin"], float(header["yllcenter"])) + bounds["xMax"] = max( + bounds["xMax"], + float(header["xllcenter"]) + float(header["ncols"]) * cellSize, + ) + bounds["yMax"] = max( + bounds["yMax"], + float(header["yllcenter"]) + float(header["nrows"]) * cellSize, + ) + + return bounds, cellSize + + +def mergeRasters(rasterFiles, bounds, mergeMethod="max"): + """Merge multiple rasters into a single raster. + + Parameters + ---------- + rasterFiles : list + List of paths to raster files + bounds : dict + Dictionary containing xMin, yMin, xMax, yMax of the union bounds + mergeMethod : str, optional + Method to use for merging overlapping cells. Options: + - 'max': maximum value (default) + - 'min': minimum value + - 'sum': sum of values + - 'count': number of overlapping valid results per cell + + Returns + ------- + mergedHeader : dict + Header dictionary containing ncols, nrows, xllcenter, yllcenter, cellsize, nodata_value + mergedData : numpy.ndarray + 2D array containing the merged raster data + """ + + # Merge data with rasterio + # If something other than min/max is wanted, it is possible to provide a custom function to merge + mergedData, outputTransform = merge(rasterFiles, method=mergeMethod, masked=True) + + mergedData = np.squeeze(mergedData) + + # Calculate dimensions for merged raster; helps checking if merged raster is correct + nCols = int((bounds["xMax"] - bounds["xMin"]) / outputTransform[0]) + nRows = int((bounds["yMax"] - bounds["yMin"]) / outputTransform[0]) + # + # # Create merged raster header + exampleRaster = rasterUtils.readRaster(rasterFiles[0]) + mergedHeader = exampleRaster["header"] + mergedHeader["ncols"] = nCols + mergedHeader["nrows"] = nRows + mergedHeader["xllcenter"] = float(bounds["xMin"]) + mergedHeader["yllcenter"] = float(bounds["yMin"]) + mergedHeader["transform"] = outputTransform + + return mergedHeader, mergedData + + +def mergeOutputRasters(cfg, avalancheDir): + """Merge output rasters (peakFiles) from all avalanche simulations. + + Parameters + ---------- + cfg : configparser.ConfigParser + Configuration containing settings: + - GENERAL.mergeOutput: If True, merge rasters + - GENERAL.mergeTypes: Types of rasters to merge (e.g., 'ppr|pfv|pft') + - GENERAL.mergeMethods: Methods to use for merging (e.g., 'max') + avalancheDir : pathlib.Path or str + Base directory where merged files will be saved + + Returns + ------- + mergedRastersDir : pathlib.Path or None + Path to the directory containing merged rasters or None if mergeOutput is False + """ + if not cfg["GENERAL"].getboolean("mergeOutput", False): + log.info("mergeOutput is False - no rasters will be merged") + return None + + # Get all avalanche directories + # with logUtils.silentLogger(): + avaDirs = findAvaDirs(avalancheDir) + if not avaDirs: + log.warning("No avalanche directories found to merge") + return None + + # Set up merged rasters directory + mergedRastersDir = pathlib.Path(avalancheDir, "mergedRasters") + if mergedRastersDir.exists(): + shutil.rmtree(str(mergedRastersDir)) + mergedRastersDir.mkdir(parents=True, exist_ok=True) + + # Get types to merge + mergeTypes = cfg["GENERAL"].get("mergeTypes").split("|") + log.info(f"Merging raster types: {mergeTypes}") + + # Get merge methods + mergeMethods = cfg["GENERAL"].get("mergeMethods", "max").lower().split("|") + log.info(f"Using merge methods: {mergeMethods}") + + # Validate merge methods + validMethods = {"max", "min", "mean", "sum", "count"} + invalidMethods = set(mergeMethods) - validMethods + if invalidMethods: + raise ValueError(f"Invalid merge methods: {invalidMethods}. Valid options are: {validMethods}") + + # Process each raster type + for rasterType in mergeTypes: + # Find all files of this type across all avalanche directories + rasterFiles = [] + for avaDir in avaDirs: + peakFilesDir = avaDir / "Outputs" / "com1DFA" / "peakFiles" + if peakFilesDir.is_dir(): + rasterFiles.extend(list(peakFilesDir.glob(f"*_{rasterType}.*"))) + + if not rasterFiles: + log.warning(f"No {rasterType} rasters found to merge") + continue + + log.info(f"Merging {len(rasterFiles)} {rasterType} rasters") + + # Get bounds and validate cell sizes + bounds, cellSize = getRasterBounds(rasterFiles) + + # Merge and save rasters + for mergeMethod in mergeMethods: + mergedHeader, mergedData = mergeRasters(rasterFiles, bounds, mergeMethod=mergeMethod) + outputPath = mergedRastersDir / f"merged_{rasterType}_{mergeMethod}" + rasterUtils.writeResultToRaster(mergedHeader, mergedData, outputPath, flip=False) + log.info(f"Saved merged {rasterType} raster (method: {mergeMethod}) to: {outputPath}") + + return mergedRastersDir diff --git a/avaframe/com7Regional/com7RegionalCfg.ini b/avaframe/com7Regional/com7RegionalCfg.ini new file mode 100644 index 000000000..0818c30ee --- /dev/null +++ b/avaframe/com7Regional/com7RegionalCfg.ini @@ -0,0 +1,55 @@ +# Configuration settings for regional run +[GENERAL] +#-------Input------- +# define the path to the regional directory, which should contain multiple avalanche directories +# (by default, set to work directly with splitInputs output location: com7Regional) +regionalDir = com7Regional + +#-------Output------- +# Whether to copy peak files to an allPeakFiles directory +copyPeakFiles = True +# Whether to move files instead of copying +moveInsteadOfCopy = False + +# Whether to merge com1DFA peakfiles to output rasters +mergeOutput = True +# Types of rasters to merge (pipe-separated list). Available options are: ppr, pft, pfv, pta, FT, FV, P, FM, Vx, Vy, Vz, TA +mergeTypes = pfv +# Methods to use for merging overlapping cells (pipe-separated list). Available options: +# - max: use maximum value +# - min: use minimum value +# - sum: sum all values +# - count: number of overlapping rasters per cell +mergeMethods = max + +[com1DFA_com1DFA_override] +# in this section, settings may be overridden for com1DFA that are used for the regional run. any parameter that is not +# set here will be taken from the default com1DFA config (if defaultConfig = True). + +# time step [s] +dt = 0.5 + +# use default com1DFA config as base configuration (True) and override following parameters +# if False and local_com1DFACfg is available use local +defaultConfig = True + +# threshold under which no remeshing is done +# set to high value to avoid raster remeshing +meshCellSizeThreshold = 100 + + +# list of simulations that shall be performed (null, ent, res, entres, available (use all available input data)) +simTypeList = null + +#+++++++++++++ Output++++++++++++ +# desired result Parameters (ppr, pft, pfv, pta, FT, FV, P, FM, Vx, Vy, Vz, TA, particles) - separated by | +resType = pfv +# saving time step, i.e. time in seconds (first and last time step are always saved) +# option 1: give an interval with start:interval in seconds +# (tStep = 0:5 - this will save desired results every 5 seconds for the full simulation) +# option 2: explicitly list all desired time steps +# (closest to actual computational time step) separated by | (example tSteps = 1|50.2|100) +# NOTE: initial and last time step are always saved! +tSteps = + +# Consider also adjusting settings in avaframeCfg.ini, such as generation of reports and plots to reduce output file size \ No newline at end of file diff --git a/avaframe/com7Regional/splitInputs.py b/avaframe/com7Regional/splitInputs.py new file mode 100644 index 000000000..750d75e0e --- /dev/null +++ b/avaframe/com7Regional/splitInputs.py @@ -0,0 +1,692 @@ +"""Module for splitting and organizing regional avalanche input data.""" + +import logging +import shapefile # pyshp +from shapely.geometry import box +import pathlib +import time + +from avaframe.in2Trans import rasterUtils +from avaframe.in3Utils import fileHandlerUtils as fU +from avaframe.in1Data import getInput +from avaframe.in3Utils.initializeProject import initializeFolderStruct +from avaframe.in2Trans import shpConversion as shpConv +from avaframe.out3Plot.outCom7Regional import createReportPlot + +# create local logger +log = logging.getLogger(__name__) + + +def splitInputsMain(avalancheDir, outputDir, cfg, cfgMain): + """Process and organize avalanche input data into individual avalanche directories based + on release area's "group" and "scenario" attributes provided in the release area file. If no + "group" attribute is provided, one avalanche directory per feature will be created (scenario is + ignored in this case). + + Parameters + ---------- + avalancheDir : pathlib.Path object + Path to input directory containing release areas (REL) and DEM files + outputDir : pathlib.Path object + Path to output directory where organized folders will be created + cfg : dict + Configuration settings containing: + - GENERAL.bufferSize : float, buffer size for DEM clipping + cfgMain : dict + Configuration settings containing: + - FLAGS.createReport : bool, whether to write report + - FLAGS.savePlot : bool, whether to save plots + + Returns + ------- + none + + Notes + ----- + Expected input directory structure: + avalancheDir/ + └── Inputs/ + ├── REL/ + │ └── *.shp # all release areas + ├── ENT/ # all entrainment areas (optional) + │ └── *.shp + ├── RES/ # all resistance areas (optional) + │ └── *.shp + └── *.asc or *.tif # digital elevation model (DEM) + """ + # Fetch the necessary input + inputSimFilesAll = getInput.getInputDataCom1DFA(avalancheDir) + + # extract release shapefile, make sure only one exists + if len(inputSimFilesAll["relFiles"]) == 1: + inputShp = inputSimFilesAll["relFiles"][0] + else: + log.error(f"Expected only one release area file, found {len(inputSimFilesAll['relFiles'])}.") + return + + # Get the input DEM + inputDEM = getInput.getDEMPath(avalancheDir) + + # Create the output directory + fU.makeADir(outputDir) + + # Step 1: Create the directory list + log.info("Initializing folder structure for each group...") + dirListGrouped = createDirList(inputShp) + log.info("Finished creating folder list") + + # Step 2: Set up avalanche directories + log.info("Initializing folder structure for each entry...") + for entry in dirListGrouped: + dirName = entry["dirName"] + initializeFolderStruct(str(outputDir / dirName), removeExisting=True) + log.debug(f"Created folder structure for '{dirName}'.") + log.info("Finished folder initialization") + + # Step 3: Split and move release areas to each directory + log.info("Splitting and moving release areas...") + splitAndMoveReleaseAreas(dirListGrouped, inputShp, outputDir) + log.info("Finished splitting and moving release areas") + + # Step 4: Clip and move DEM + log.info("Clipping and moving DEM...") + groupExtents = clipDEMByReleaseGroup(dirListGrouped, inputDEM, outputDir, cfg) + log.info("Finished clipping and moving of DEM") + + # Step 5: Clip and move optional input (currently only ENT and RES) + log.info("Clipping and moving optional input...") + groupFeatures = clipAndMoveOptionalInput(inputSimFilesAll, outputDir, groupExtents) + log.info("Finished clipping and moving optional input") + + # Step 6: Divide release areas into scenarios + log.info("Separating release areas by scenarios...") + splitByScenarios(dirListGrouped, outputDir) + log.info("Finished separating by scenarios") + + # Step 7: Write reports + if cfgMain["FLAGS"].getboolean("createReport"): + log.info("Writing reports...") + writeScenarioReport(dirListGrouped, outputDir) + if cfgMain["FLAGS"].getboolean("savePlot"): + createReportPlots(dirListGrouped, inputDEM, outputDir, groupExtents, groupFeatures) + log.info("Finished writing reports") + + +def createDirList(inputShp): + """Create a list of entries from each feature in the input shapefile, grouped by the 'group' attribute. + + Parameters + ---------- + inputShp: pathlib.Path object + path to input shapefile + + Returns + ------- + dirListGrouped: list + list of dictionaries containing dirName (group name), properties list, and geometries list, + where features are grouped by their 'group' attribute + """ + fields, fieldNames, properties, geometries, srs = shpConv.readShapefile(inputShp) + + # Create dictionary to store groups + groups = {} + unnamedCount = 1 + + for props, geom in zip(properties, geometries): + propsLower = {key.lower(): value for key, value in props.items()} # Handle case sensitivity + + # Get group name from 'group' attribute, fallback to unnamed if not present + groupName = propsLower.get("group", "").strip() or f"{str(unnamedCount).zfill(5)}" + if not propsLower.get("group", "").strip(): + unnamedCount += 1 + log.info(f"No 'group' field or empty group found in {inputShp}, using '{groupName}'") + + # Initialize group if not exists + if groupName not in groups: + groups[groupName] = { + "dirName": groupName, + "properties": [], + "geometries": [], + } + + # Add feature to group + groups[groupName]["properties"].append(props) + groups[groupName]["geometries"].append(geom) + + # Convert dictionary to list and sort by dirName + dirListGrouped = list(groups.values()) + dirListGrouped.sort(key=lambda x: x["dirName"].lower()) + + # Log total number of features + totalFeatures = sum(len(group["geometries"]) for group in dirListGrouped) + log.info(f"Found '{totalFeatures}' features that were organized into '{len(dirListGrouped)}' groups") + + return dirListGrouped + + +def splitAndMoveReleaseAreas(dirList, inputShp, outputDir): + """Split release areas into individual shapefiles and write them to their respective folders. + + Parameters + ---------- + dirList: list + list of dictionaries containing dirName, properties list, and geometries list + inputShp: pathlib.Path object + path to input shapefile + outputDir: pathlib.Path object + path to output directory where folders will be created + + Returns + ------- + none + """ + # Read the input shapefile + fields, fieldNames, properties, geometries, srs = shpConv.readShapefile(inputShp) + + featuresByName = {} + for entry in dirList: + name = entry["dirName"] # Get release area name + # Group entries with the same name + if name not in featuresByName: + featuresByName[name] = [] + # add corresponding properties and geometries + for i, properties in enumerate(entry["properties"]): + featuresByName[name].append((properties, entry["geometries"][i])) + + # Write shapefiles to their respective folders + for name, features in featuresByName.items(): + shpOutPath = outputDir / name / "Inputs" / "REL" / name + shpConv.writeShapefile(shpOutPath, fields, fieldNames, features, srs) + log.debug(f"Saved release area to '{shpOutPath}'.") + + +def checkFeatureIsolation(geometries, properties, bufferSize, groupName): + """Check if any feature in the group is isolated from all others. + + A feature is considered isolated if its buffered bounding box does not overlap + with any other feature's buffered bounding box in the group. + + Parameters + ---------- + geometries: list + List of geometry objects to check + properties: list + List of dictionaries containing properties for each geometry + bufferSize: float + Buffer size to use when creating bounding boxes + groupName: str + Name of the group, used for error messages + + Raises + ------ + ValueError + If any feature is isolated from all others in the group + """ + # Skip check if only one feature + if len(geometries) <= 1: + log.debug(f"Group '{groupName}' has only one feature, proceeding without isolation check.") + return + + # Create buffered bounding boxes for each feature + boundingBoxes = [] + for geom in geometries: + center = geom.centroid + + # Calculate bounding box for this feature + currXMin = center.x - bufferSize + currYMin = center.y - bufferSize + currXMax = center.x + bufferSize + currYMax = center.y + bufferSize + + # Update group extent + boundingBoxes.append(box(currXMin, currYMin, currXMax, currYMax)) + + # Check each feature's bounding box against all others + for i, bbox in enumerate(boundingBoxes): + hasOverlap = False + for j, otherBbox in enumerate(boundingBoxes): + if i != j and bbox.intersects(otherBbox): + hasOverlap = True + break + + if not hasOverlap: + # Find feature name regardless of case (NAME, name, Name etc.) + featureProps = {key.lower(): value for key, value in properties[i].items()} + featureName = featureProps.get("name", f"unnamed feature {i+1}").strip() + + message = f"Feature '{featureName}' in group '{groupName}' is isolated from all other features - consider assigning it to a different group" + log.error(message) + raise ValueError(message) + + +def clipDEMByReleaseGroup(dirList, inputDEM, outputDir, cfg): + """Clip the DEM to include all features in each release group. Returns an error if any feature in a group is isolated. + + Parameters + ---------- + dirList : list + List of dictionaries containing dirName, and geometries list + inputDEM : pathlib.Path + Path to input DEM file + outputDir : pathlib.Path + Path to output directory where clipped DEMs will be saved + cfg : configparser object + Configuration settings containing: + - GENERAL.bufferSize : float + Size of buffer to add around release areas Configuration settings + + Returns + ------- + groupExtents : dict + Dictionary with dirName as key and (xMin, xMax, yMin, yMax) as value. + The extents are reduced by one pixel on each side to ensure DEM extents + are larger than clip extents of other input. + """ + # Read input DEM + demData = rasterUtils.readRaster(inputDEM) + header = demData["header"] + raster = demData["rasterData"] + cellSize = header["cellsize"] + xOrigin = header["xllcenter"] + yOrigin = header["yllcenter"] + nRows = header["nrows"] + nCols = header["ncols"] + + # Process each group + groupExtents = {} + for entry in dirList: + dirName = entry["dirName"] + geometries = entry["geometries"] + properties = entry["properties"] + + if not geometries: + message = f"No geometries found for {dirName}" + log.error(message) + raise ValueError(message) + + # Check if any features in the group are isolated + bufferSize = cfg["GENERAL"].getfloat("bufferSize") + checkFeatureIsolation(geometries, properties, bufferSize, dirName) + + # Get extent of all geometries in group + bounds = [geom.bounds for geom in geometries] + xMins, yMins, xMaxs, yMaxs = zip(*bounds) + + # Calculate extent with buffer + bufferSize = cfg["GENERAL"].getfloat("bufferSize") + xMin = min(xMins) - bufferSize + xMax = max(xMaxs) + bufferSize + yMin = min(yMins) - bufferSize + yMax = max(yMaxs) + bufferSize + groupExtents[dirName] = (xMin, xMax, yMin, yMax) # Store extent for this group + + # Convert extent to grid indices (using top-left origin) + colStart = max(0, int((xMin - xOrigin) / cellSize)) + colEnd = min(nCols, int((xMax - xOrigin) / cellSize) + 1) + + # Convert y-coordinates to row indices (flipped for bottom-left origin) + rowStart = max(0, int((yOrigin + nRows * cellSize - yMax) / cellSize)) + rowEnd = min(nRows, int((yOrigin + nRows * cellSize - yMin) / cellSize) + 1) + + # Ensure valid row indices + if rowEnd <= rowStart: + log.warning(f"Invalid row indices calculated for {dirName}: start={rowStart}, end={rowEnd}") + continue + + # Flip row indices for bottom-left origin + rowStart, rowEnd = nRows - rowEnd, nRows - rowStart + + # Clip the DEM data + clippedData = raster[rowStart:rowEnd, colStart:colEnd] + + # Create header for clipped DEM + clippedHeader = header.copy() + clippedHeader["ncols"] = colEnd - colStart + clippedHeader["nrows"] = rowEnd - rowStart + clippedHeader["xllcenter"] = xOrigin + colStart * cellSize + clippedHeader["yllcenter"] = yOrigin + rowStart * cellSize + + # Update transformation matrix for clipped DEM + clippedHeader["transform"] = rasterUtils.transformFromASCHeader(clippedHeader) + + # Write clipped DEM + outputDEM = outputDir / dirName / "Inputs" / f"{dirName}_DEM" + rasterUtils.writeResultToRaster(clippedHeader, clippedData, outputDEM, flip=True) + log.debug(f"Clipped DEM saved to: {outputDEM}") + + # Store final DEM extents (reduced by one pixel on each side to ensure DEM > clip extents of other input) + xMinDEM = clippedHeader["xllcenter"] + (cellSize * 0.5) + yMinDEM = clippedHeader["yllcenter"] + (cellSize * 0.5) + xMaxDEM = clippedHeader["xllcenter"] + (clippedHeader["ncols"] * cellSize) - (cellSize * 0.5) + yMaxDEM = clippedHeader["yllcenter"] + (clippedHeader["nrows"] * cellSize) - (cellSize * 0.5) + groupExtents[dirName] = (xMinDEM, xMaxDEM, yMinDEM, yMaxDEM) + + return groupExtents + + +def clipAndMoveOptionalInput(allSimInputFiles, outputDir, groupExtents): + """Clip and move ENT and RES files based on group DEM extent. + + Parameters + ---------- + allSimInputFiles: dict + With all input information for a com1DFA sim + outputDir : pathlib.Path + Path to output directory where clipped files will be saved + groupExtents : dict + Dictionary with dirName as key and (xMin, xMax, yMin, yMax) as value, + containing the DEM clipping extents for each group + + Returns + ------- + groupFeatures : dict + Dictionary containing clipped features for each group and type + {dirName: {'ENT': [...], 'RES': [...]}} + """ + groupFeatures = {} + # Process ENT and RES directories + for dirType in ["ENT", "RES"]: + + if dirType == "ENT": + if allSimInputFiles["entFile"]: + shpFile = allSimInputFiles["entFile"] + else: + log.info("No entrainment file found") + continue + + if dirType == "RES": + if allSimInputFiles["resFile"]: + shpFile = allSimInputFiles["resFile"] + else: + log.info("No resistance file found") + continue + + # Read shapefile + fields, fieldNames, properties, geometries, srs = shpConv.readShapefile(shpFile) + + # Process each output directory that has extents + for entry in outputDir.iterdir(): + if not entry.is_dir() or entry.name not in groupExtents: + continue + + # Get extent + xMin, xMax, yMin, yMax = groupExtents[entry.name] + scenarioBbox = box(xMin, yMin, xMax, yMax) + + # Initialize group in dictionary if not exists + if entry.name not in groupFeatures: + groupFeatures[entry.name] = {"ENT": [], "RES": []} + + # Clip geometries with groups DEM extent + clippedFeatures = [] + for prop, geom in zip(properties, geometries): + if geom.intersects(scenarioBbox): + clippedGeom = geom.intersection(scenarioBbox) + if not clippedGeom.is_empty: + clippedFeatures.append((prop, clippedGeom)) + groupFeatures[entry.name][dirType].append(clippedGeom) + + if not clippedFeatures: + log.debug(f"No {dirType} features intersect with DEM extent for {entry.name}") + continue + + # Create output directory and save clipped shp + targetDir = entry / "Inputs" / dirType + fU.makeADir(targetDir) + outputPath = targetDir / f"{entry.name}_{dirType}.shp" + shpConv.writeShapefile(outputPath, fields, fieldNames, clippedFeatures, srs) + log.debug(f"Clipped {dirType} shapefile saved to: {outputPath}") + + return groupFeatures + + +def getScenarioGroups(inputShp, fieldNames): + """Group shapefile records by their scenario attribute. + + Parameters + ---------- + inputShp : pathlib.Path + Path to input shapefile + fieldNames : list + List of field names in the shapefile + + Returns + ------- + scenarios: dict + Dictionary mapping scenario names to lists of shape records + """ + scenarios = {} + for shapeRecord in shapefile.Reader(str(inputShp)).iterShapeRecords(): + properties = {k.lower(): v for k, v in zip(fieldNames, shapeRecord.record)} + scenarioValues = properties.get("scenario", "").split(",") + for scenario in scenarioValues: + # Check if scenario value is empty and set flag + if scenario.strip() == "": + scenario = "NULL" + # If scenario is not in scenarios dict, add it + if scenario not in scenarios: + scenarios[scenario] = [] + scenarios[scenario].append(shapeRecord) + return scenarios + + +def writeScenarioShapefile(outputShp, records, fields, fieldNames, srs): + """Write a shapefile for a specific scenario. + + Parameters + ---------- + outputShp : pathlib.Path + Path where to write the shapefile + records : list + List of shape records for this scenario + fields : list + List of field definitions + fieldNames : list + List of field names + srs : str + Spatial reference system string + """ + # Filter out the scenario attribute + shapeFeatures = [(dict(zip(fieldNames, record.record)), record.shape) for record in records] + filteredFields = [field for field in fields if field[0].lower() != "scenario"] + filteredFieldNames = [name for name in fieldNames if name.lower() != "scenario"] + + # Write the shapefile + shpConv.writeShapefile(outputShp, filteredFields, filteredFieldNames, shapeFeatures, srs) + + +def splitByScenarios(dirList, outputDir): + """Split release areas into separate shapefiles based on their scenario attribute. + + Parameters + ---------- + dirList: list + list of dictionaries containing dirName and list of geometries + outputDir: pathlib.Path object + path to output directory + + Returns + ------- + none + + Notes + ----- + - If a feature has no scenario attribute or it's empty, it will be marked as 'NULL' and grouped together with other 'NULL' features + - Intermediate shapefiles are deleted or renamed after scenario splitting + """ + totalInputFiles = 0 + totalScenarioFiles = 0 + + # Loop through each folder + for folder in dirList: + inputShp = pathlib.Path(outputDir) / folder["dirName"] / "Inputs" / "REL" / folder["dirName"] + fields, fieldNames, properties, geometries, srs = shpConv.readShapefile(inputShp) + totalInputFiles += 1 + + # Get the scenario attribute values + if "scenario" in map(str.lower, fieldNames): + # Group records by scenario + scenarios = getScenarioGroups(inputShp, fieldNames) + + # Write a shapefile for each scenario + for scenario, records in scenarios.items(): + if all(scenario == "NULL" for scenario in scenarios): + outputShp = ( + pathlib.Path(outputDir) + / folder["dirName"] + / "Inputs" + / "REL" + / f"{folder['dirName']}_REL" + ) + elif scenario == "NULL": + outputShp = ( + pathlib.Path(outputDir) + / folder["dirName"] + / "Inputs" + / "REL" + / f"{folder['dirName']}_NULL" + ) + else: + outputShp = ( + pathlib.Path(outputDir) + / folder["dirName"] + / "Inputs" + / "REL" + / f"{folder['dirName']}_{scenario}" + ) + + writeScenarioShapefile(outputShp, records, fields, fieldNames, srs) + totalScenarioFiles += 1 + + # Delete the intermediate shapefile + for ext in [".shp", ".shx", ".dbf", ".prj"]: + if (inputShp.with_suffix(ext)).exists(): + (inputShp.with_suffix(ext)).unlink() + else: + # If no scenario attribute exists, rename the file (necessary for further processing) + outputShp = ( + pathlib.Path(outputDir) / folder["dirName"] / "Inputs" / "REL" / f"{folder['dirName']}_REL" + ) + shapeFeatures = [ + (dict(zip(fieldNames, record.record)), record.shape) + for record in shapefile.Reader(str(inputShp)).iterShapeRecords() + ] + shpConv.writeShapefile(outputShp, fields, fieldNames, shapeFeatures, srs) + for ext in [".shp", ".shx", ".dbf", ".prj"]: + if (inputShp.with_suffix(ext)).exists(): + (inputShp.with_suffix(ext)).unlink() + + if totalScenarioFiles == 0: + log.info("No 'scenario' attribute or only 'NULL' found in release area shapefiles, continuing") + else: + log.info(f"Split '{totalInputFiles}' release area shapefiles into '{totalScenarioFiles}' scenarios") + + +def writeScenarioReport(dirListGrouped, outputDir): + """Create a report in txt format listing all scenarios and their associated features. + + Parameters + ---------- + dirListGrouped : list + list of dictionaries containing dirName and list of geometries + outputDir : pathlib.Path + Path to output directory where the report will be saved + + Returns + ------- + none + """ + reportPath = outputDir / "splitInputs_scenarioReport.txt" + + with open(reportPath, "w") as f: + f.write("SCENARIO REPORT\n") + f.write("==============\n") + f.write(f"Generated: {time.strftime('%Y-%m-%d %H:%M:%S')}\n\n") + + # Process each group and their scenarios + for group in sorted(dirListGrouped, key=lambda x: x["dirName"].lower()): + dirName = group["dirName"] + f.write(f"Group: {dirName}\n") + f.write("-" * (len(dirName) + 7) + "\n\n") + + relPath = pathlib.Path(outputDir) / dirName / "Inputs" / "REL" + scenarioFiles = sorted(relPath.glob(f"{dirName}_*.shp"), key=lambda x: x.stem.split("_")[-1]) + + if not scenarioFiles: + f.write("No scenarios found\n\n") + continue + + # Write release areas for each scenario + for scenFile in scenarioFiles: + fields, fieldNames, properties, geometries, _ = shpConv.readShapefile(scenFile) + scenName = scenFile.stem.split("_")[-1] + + f.write(f"Scenario: {scenName}\n") + f.write(f"No. of release areas: {len(geometries)}\n") + + if "name" in map(str.lower, fieldNames): # Handle case sensitivity + nameIdx = [i for i, name in enumerate(fieldNames) if name.lower() == "name"][0] + with shapefile.Reader(str(scenFile)) as shp: + records = sorted(shp.records(), key=lambda x: x[nameIdx].lower()) + for record in records: + f.write(f"- {record[nameIdx]}\n") + else: + for i in range(len(geometries)): + f.write(f"- Release Area {i+1}\n") + f.write("\n") + + # Write total entrainment and resistance areas for the group + entPath = pathlib.Path(outputDir) / dirName / "Inputs" / "ENT" + resPath = pathlib.Path(outputDir) / dirName / "Inputs" / "RES" + + entFiles = list(entPath.glob(f"{dirName}_*.shp")) + if entFiles: + totalEnt = sum(len(shpConv.readShapefile(ef)[3]) for ef in entFiles) + if totalEnt > 0: + f.write(f"No. of entrainment areas: {totalEnt}\n") + + resFiles = list(resPath.glob(f"{dirName}_*.shp")) + if resFiles: + totalRes = sum(len(shpConv.readShapefile(rf)[3]) for rf in resFiles) + if totalRes > 0: + f.write(f"No. of resistance areas: {totalRes}\n") + + f.write("\n") + + log.info(f"Scenario report written to: {reportPath}") + + +def createReportPlots(dirListGrouped, inputDEM, outputDir, groupExtents, groupFeatures): + """Write visual reports summarizing the split inputs operation. + + Creates two visual reports in PNG format: + 1. Basic report showing DEM extent and release areas + 2. Optional features report showing DEM extent with entrainment and resistance areas + + Parameters + ---------- + dirListGrouped : list + List of dictionaries containing dirName and list of geometries + inputDEM : pathlib.Path + Path to input DEM file + outputDir : pathlib.Path + Path to output directory where reports will be saved + groupExtents : dict + Dictionary with dirName as key and (xMin, xMax, yMin, yMax) as value, + containing the DEM clipping extents for each group + groupFeatures : dict + Dictionary containing clipped features for each group and type + + Returns + ------- + none + """ + # Create basic features report + basicPath = createReportPlot(dirListGrouped, inputDEM, outputDir, groupExtents, groupFeatures, "basic") + log.info(f"Visual report (basic) written to: {basicPath}") + + # Create optional features report + optionalPath = createReportPlot( + dirListGrouped, inputDEM, outputDir, groupExtents, groupFeatures, "optional" + ) + log.info(f"Visual report (optional) written to: {optionalPath}") diff --git a/avaframe/com7Regional/splitInputsCfg.ini b/avaframe/com7Regional/splitInputsCfg.ini new file mode 100644 index 000000000..7c89a7349 --- /dev/null +++ b/avaframe/com7Regional/splitInputsCfg.ini @@ -0,0 +1,10 @@ +### Config File - This file contains the main settings for splitInputs +## Copy to local_splitInputsCfg.ini and set your parameters +# This file is part of Avaframe. + +[GENERAL] +# Splitting length around each release area group [m]. +# The value for ``bufferSize`` is added to each direction (+x, -x, +y, -y) from the borders of a union bounding box with +# the maximum x-y extent of all release features in the group, to define the total group extent. +# Consider adjusting value according to the expected maximum runout length of your avalanches. +bufferSize = 2500 \ No newline at end of file diff --git a/avaframe/data/avaMalRegional/Inputs/ENT/entMal.cpg b/avaframe/data/avaMalRegional/Inputs/ENT/entMal.cpg new file mode 100644 index 000000000..3ad133c04 --- /dev/null +++ b/avaframe/data/avaMalRegional/Inputs/ENT/entMal.cpg @@ -0,0 +1 @@ +UTF-8 \ No newline at end of file diff --git a/avaframe/data/avaMalRegional/Inputs/ENT/entMal.dbf b/avaframe/data/avaMalRegional/Inputs/ENT/entMal.dbf new file mode 100644 index 000000000..64c87950b Binary files /dev/null and b/avaframe/data/avaMalRegional/Inputs/ENT/entMal.dbf differ diff --git a/avaframe/data/avaMalRegional/Inputs/ENT/entMal.prj b/avaframe/data/avaMalRegional/Inputs/ENT/entMal.prj new file mode 100644 index 000000000..3d616c54e --- /dev/null +++ b/avaframe/data/avaMalRegional/Inputs/ENT/entMal.prj @@ -0,0 +1 @@ +PROJCS["MGI_Austria_Lambert",GEOGCS["GCS_MGI",DATUM["D_MGI",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",400000.0],PARAMETER["Central_Meridian",13.3333333333333],PARAMETER["Standard_Parallel_1",49.0],PARAMETER["Standard_Parallel_2",46.0],PARAMETER["Latitude_Of_Origin",47.5],UNIT["Meter",1.0]] \ No newline at end of file diff --git a/avaframe/data/avaMalRegional/Inputs/ENT/entMal.shp b/avaframe/data/avaMalRegional/Inputs/ENT/entMal.shp new file mode 100644 index 000000000..f2084d69d Binary files /dev/null and b/avaframe/data/avaMalRegional/Inputs/ENT/entMal.shp differ diff --git a/avaframe/data/avaMalRegional/Inputs/ENT/entMal.shx b/avaframe/data/avaMalRegional/Inputs/ENT/entMal.shx new file mode 100644 index 000000000..6ccfdd7be Binary files /dev/null and b/avaframe/data/avaMalRegional/Inputs/ENT/entMal.shx differ diff --git a/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.cpg b/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.cpg new file mode 100644 index 000000000..3ad133c04 --- /dev/null +++ b/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.cpg @@ -0,0 +1 @@ +UTF-8 \ No newline at end of file diff --git a/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.dbf b/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.dbf new file mode 100644 index 000000000..da336fad0 Binary files /dev/null and b/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.dbf differ diff --git a/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.prj b/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.prj new file mode 100644 index 000000000..3d616c54e --- /dev/null +++ b/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.prj @@ -0,0 +1 @@ +PROJCS["MGI_Austria_Lambert",GEOGCS["GCS_MGI",DATUM["D_MGI",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",400000.0],PARAMETER["Central_Meridian",13.3333333333333],PARAMETER["Standard_Parallel_1",49.0],PARAMETER["Standard_Parallel_2",46.0],PARAMETER["Latitude_Of_Origin",47.5],UNIT["Meter",1.0]] \ No newline at end of file diff --git a/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.shp b/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.shp new file mode 100644 index 000000000..5b85787b2 Binary files /dev/null and b/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.shp differ diff --git a/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.shx b/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.shx new file mode 100644 index 000000000..3dcd86c00 Binary files /dev/null and b/avaframe/data/avaMalRegional/Inputs/REL/relMalRegional.shx differ diff --git a/avaframe/data/avaMalRegional/Inputs/avaMal.tif b/avaframe/data/avaMalRegional/Inputs/avaMal.tif new file mode 100644 index 000000000..8757a0de9 Binary files /dev/null and b/avaframe/data/avaMalRegional/Inputs/avaMal.tif differ diff --git a/avaframe/in2Trans/shpConversion.py b/avaframe/in2Trans/shpConversion.py index b6493bfc1..620d81ced 100644 --- a/avaframe/in2Trans/shpConversion.py +++ b/avaframe/in2Trans/shpConversion.py @@ -1,14 +1,15 @@ """ - Conversion functions to read/ write Shape files +Conversion functions to read/ write Shape files """ - import pathlib import shapefile import copy import numpy as np import logging +from shapely.geometry import shape + # create local logger log = logging.getLogger(__name__) @@ -91,7 +92,6 @@ def SHP2Array(infile, defname=None): tilt_value = None direc_value = None offset_value = None - # get coordinate system sks = getSHPProjection(infile) @@ -116,13 +116,13 @@ def SHP2Array(infile, defname=None): maxdepthList = [] tiltList = [] direcList = [] - offsetList = [] - - Length = np.empty((0)) - Start = np.empty((0)) - Coordx = np.empty((0)) - Coordy = np.empty((0)) - Coordz = np.empty((0)) + offsetList = [] + + Length = np.empty(0) + Start = np.empty(0) + Coordx = np.empty(0) + Coordy = np.empty(0) + Coordz = np.empty(0) start = 0 nParts = [] @@ -297,7 +297,6 @@ def readThickness(infile, defname=None): """ # Input shapefile sf = shapefile.Reader(str(infile)) - infile = pathlib.Path(infile) thickness = None id = None @@ -547,20 +546,100 @@ def writePoint2SHPfile(pointDict, pointName, fileName): """ fileName = str(fileName) w = shapefile.Writer(fileName) - w.field('name', 'C') - if isinstance(pointDict['x'], (float, np.float64)) and isinstance(pointDict['y'], (float, np.float64)): - w.point(pointDict['x'], pointDict['y']) - elif isinstance(pointDict['x'], (np.ndarray)) and isinstance(pointDict['y'], (np.ndarray)): - if len(pointDict['x']) > 1 or len(pointDict['y']) > 1: - message = 'Length of pointDict is not allowed to exceed one' + w.field("name", "C") + if isinstance(pointDict["x"], (float, np.float64)) and isinstance(pointDict["y"], (float, np.float64)): + w.point(pointDict["x"], pointDict["y"]) + elif isinstance(pointDict["x"], (np.ndarray)) and isinstance(pointDict["y"], (np.ndarray)): + if len(pointDict["x"]) > 1 or len(pointDict["y"]) > 1: + message = "Length of pointDict is not allowed to exceed one" log.error(message) raise ValueError(message) else: - w.point(pointDict['x'][0], pointDict['y'][0]) + w.point(pointDict["x"][0], pointDict["y"][0]) else: - message = 'Format of point is neither float nor array of length 1' + message = "Format of point is neither float nor array of length 1" log.error(message) raise TypeError(message) w.record(pointName) w.close() return fileName + + +def writeShapefile(outputPath, fields, fieldNames, features, srs=None): + """Write features to a shapefile with given fields and properties. + + Parameters + ---------- + outputPath: pathlib.Path + path where the shapefile will be written + fields: list + the fields of the shapefile + fieldNames: list + the names of the fields + features: list + list of tuples containing (properties, geometry) for each feature + srs: str, optional + the spatial reference system for the .prj file + + Returns + ------- + none + """ + with shapefile.Writer(str(outputPath)) as dst: + for field in fields: + dst.field(*field) + for properties, geometry in features: + dst.shape(geometry) + record = [properties.get(fieldName, "") for fieldName in fieldNames] + dst.record(*record) + + if srs is not None: + prjOutPath = outputPath.with_suffix(".prj") + with open(prjOutPath, "w") as prjFile: + prjFile.write(srs) + + +def readShapefile(inputShp): + """Read the fields, properties, geometries, and spatial reference of an input shapefile. + To be used in combination with shapefile.Reader. Could be expanded upon to get e.g. + shapeTypes, bounds, numFeatures and metadata if needed. This does pure reading, no special + checks for valid attributes are performed. + + Parameters + ---------- + inputShp : pathlib.Path + Path to input shapefile + + Returns + ------- + fields: list + the fields of the input shapefile + fieldNames: list + the names of the fields + properties: list + a list of dictionaries containing the properties of each feature + geometries: list + a list of geometry objects + srs: str + the spatial reference system fetched from eventual .prj file + """ + with shapefile.Reader(str(inputShp)) as src: + fields = src.fields[1:] # Skip deletion flag + fieldNames = [field[0] for field in fields] + properties = [] + geometries = [] + for shapeRecord in src.iterShapeRecords(): + if shapeRecord.shape.shapeType == shapefile.NULL: + log.warning(f"Skipping NULL shape in {inputShp}") + continue + properties.append(dict(zip(fieldNames, shapeRecord.record))) + geometries.append(shape(shapeRecord.shape.__geo_interface__)) + + srs = None + # Check if .prj file exists and read it + srsfile = inputShp.with_suffix(".prj") + if srsfile.is_file(): + with open(srsfile, "r") as f: + srs = f.read().strip() + + return fields, fieldNames, properties, geometries, srs diff --git a/avaframe/in3Utils/cfgUtils.py b/avaframe/in3Utils/cfgUtils.py index 8a62f6192..69cd0a5e5 100644 --- a/avaframe/in3Utils/cfgUtils.py +++ b/avaframe/in3Utils/cfgUtils.py @@ -927,7 +927,7 @@ def getNumberOfProcesses(cfgMain, nSims): # if number of sims is lower than nCPU nCPU = min(nCPU, nSims) - log.info("Number of simulations to perform: %s " % nSims) + log.info("Number of tasks to perform: %s " % nSims) log.info("Taking %s cpu cores out of maximum of %s cores." % (nCPU, maxCPU)) return nCPU @@ -985,4 +985,3 @@ def cfgToRcf(cfg, fileName): key = key.strip() f.write(f"{key:<40}{value}\n") f.write("#\n") - diff --git a/avaframe/in3Utils/logUtils.py b/avaframe/in3Utils/logUtils.py index e186a9776..08e7971bd 100644 --- a/avaframe/in3Utils/logUtils.py +++ b/avaframe/in3Utils/logUtils.py @@ -1,44 +1,44 @@ """ - Defining a writable object to write the config to the log file +Defining a writable object to write the config to the log file - This file is part of Avaframe. +This file is part of Avaframe. """ + import logging import logging.config -import os import pathlib -from pathlib import Path, PureWindowsPath from datetime import datetime import avaframe.version as gv +from contextlib import contextmanager -def writeCfg2Log(cfg, cfgName='Unspecified'): - """ write a configparser object to log file""" +def writeCfg2Log(cfg, cfgName="Unspecified"): + """write a configparser object to log file""" - log = logging.getLogger('avaframe') + log = logging.getLogger("avaframe") - log.debug('Writing cfg for: %s', cfgName) + log.debug("Writing cfg for: %s", cfgName) for section in cfg.sections(): - log.info('\t%s', section) + log.info("\t%s", section) for key in dict(cfg[section]): - log.info('\t\t%s : %s', key, cfg[section][key]) + log.info("\t\t%s : %s", key, cfg[section][key]) -def initiateLogger(targetDir, logName='runLog', modelInfo=''): - """ Initiates logger object based on setup in logging.conf +def initiateLogger(targetDir, logName="runLog", modelInfo=""): + """Initiates logger object based on setup in logging.conf - Parameters - ---------- - targetDir: str - folder to save log file to - logName : str - filename of log file; optional; defaults to runLog.log - modelInfo: str - if not emtpy add info on modelInfo to log + Parameters + ---------- + targetDir: str + folder to save log file to + logName : str + filename of log file; optional; defaults to runLog.log + modelInfo: str + if not emtpy add info on modelInfo to log - Returns - ------- - log : logging object + Returns + ------- + log : logging object """ @@ -46,29 +46,58 @@ def initiateLogger(targetDir, logName='runLog', modelInfo=''): now = datetime.now() dtString = now.strftime("%Y%m%d_%Hh%Mm%Ss") - logFileName = pathlib.Path(targetDir, logName+'_'+dtString+'.log') + logFileName = pathlib.Path(targetDir, logName + "_" + dtString + ".log") # get path of module and generate logging.conf file path logConfPath = pathlib.Path(__file__).parents[0] - logConfFile = logConfPath / 'logging.conf' + logConfFile = logConfPath / "logging.conf" # clean logFileName for special characters # very hacky way, but no better way found working for both linux and windows # Add further special characters here as needed... cleanName = str(logFileName) - cleanName = cleanName.replace('\'','\\\'') - cleanName = cleanName.replace('\\', '/') + cleanName = cleanName.replace("'", "\\'") + cleanName = cleanName.replace("\\", "/") - logging.config.fileConfig(fname=logConfFile, - defaults={'logfilename': cleanName}, - disable_existing_loggers=False) - log = logging.getLogger('avaframe') + logging.config.fileConfig( + fname=logConfFile, + defaults={"logfilename": cleanName}, + disable_existing_loggers=False, + ) + log = logging.getLogger("avaframe") - log.info('Started logging at: %s', dtString) - log.info('Also logging to: %s', logFileName) - log.info('AvaFrame Version: %s', gv.getVersion()) + log.info("Started logging at: %s", dtString) + log.info("Also logging to: %s", logFileName) + log.info("AvaFrame Version: %s", gv.getVersion()) - if modelInfo != '': - log.info('Used by: %s' % modelInfo) + if modelInfo != "": + log.info("Used by: %s" % modelInfo) return log + + +@contextmanager +def silentLogger(loggerName="avaframe"): + """Context manager to temporarily silence a logger. + + Parameters + ---------- + loggerName : str + Name of the logger to silence. Defaults to 'avaframe' + + Example + ------- + >>> with silentLogger(): + >>> # This code block will not produce any log output + >>> function_with_noisy_logs() + """ + logger = logging.getLogger(loggerName) + # Store the original level + originalLevel = logger.level + # Temporarily set the level to ERROR (will suppress INFO and DEBUG messages) + logger.setLevel(logging.ERROR) + try: + yield + finally: + # Restore the original level + logger.setLevel(originalLevel) diff --git a/avaframe/out3Plot/outCom7Regional.py b/avaframe/out3Plot/outCom7Regional.py new file mode 100644 index 000000000..92c84a9ec --- /dev/null +++ b/avaframe/out3Plot/outCom7Regional.py @@ -0,0 +1,206 @@ +import matplotlib as mpl +from matplotlib import pyplot as plt +from matplotlib.patches import Rectangle, Patch +from shapely import MultiPolygon + +from avaframe.in2Trans import rasterUtils +from avaframe.out3Plot import plotUtils as pU + + +def createReportPlot(dirListGrouped, inputDEM, outputDir, groupExtents, groupFeatures, reportType): + """Create a visual report showing the DEM extent with either basic or optional inputs. + + Parameters + ---------- + dirListGrouped : list + List of dictionaries containing dirName and list of geometries + inputDEM : pathlib.Path + Path to input DEM file + outputDir : pathlib.Path + Path to output directory where the report will be saved + groupExtents : dict + Dictionary with dirName as key and (xMin, xMax, yMin, yMax) as value, + containing the DEM clipping extents for each group + groupFeatures : dict + Dictionary containing clipped features for each group and type + reportType : str + Type of report to create, either 'basic' or 'optional' + - 'basic': Shows DEM extent and release areas only + - 'optional': Shows DEM extent with entrainment and resistance areas + + Returns + ------- + reportPath: pathlib.Path + Path to the generated report image + """ + # Set up figure + plt.figure(figsize=(10, 8)) + ax = plt.subplot(1, 1, 1) + + # Read and plot DEM + demData = rasterUtils.readRaster(inputDEM) + header = demData["header"] + cellSize = header["cellsize"] + xMin = header["xllcenter"] + yMin = header["yllcenter"] + xMax = xMin + cellSize * header["ncols"] + yMax = yMin + cellSize * header["nrows"] + im = ax.imshow( + demData["rasterData"], + extent=[xMin, xMax, yMin, yMax], + cmap=pU.cmapDEM.reversed(), + alpha=1, + origin="lower", + zorder=1, + ) + + # Custom color scheme for groups + colors = [ + "#0EF8EA", + "#FFA500", + "#C71585", + "#00FF00", + "#FF4500", + "#800080", + "#ADFF2F", + "#FF6347", + "#8A2BE2", + "#FFFF00", + "#FF0000", + ] + + # Plot groups + for idx, group in enumerate(dirListGrouped): + dirName = group["dirName"] + color = colors[idx % len(colors)] + # Plot DEM extent using groupExtents + if dirName in groupExtents: + xMin, xMax, yMin, yMax = groupExtents[dirName] + rect = Rectangle( + (xMin, yMin), + xMax - xMin, + yMax - yMin, + fill=False, + linestyle="--", + color=color, + ) + ax.add_patch(rect) + + if reportType == "basic": + # Plot release areas + for geom in group["geometries"]: + for x, y in getExteriorCoords(geom): + plt.fill(x, y, alpha=1.0, color=color) + else: + # Plot optional features + if dirName in groupFeatures: + for geom in groupFeatures[dirName].get("ENT", []): + for x, y in getExteriorCoords(geom): + plt.fill(x, y, alpha=0.3, color=color, edgecolor="none") + for geom in groupFeatures[dirName].get("RES", []): + for x, y in getExteriorCoords(geom): + plt.fill( + x, + y, + alpha=0.5, + color=color, + hatch="xxxx", + fill=False, + edgecolor=color, + linewidth=0.5, + ) + + # Place group label using groupExtents + plt.text( + xMin, + yMax, + dirName, + color=color, + fontsize=8, + transform=mpl.transforms.offset_copy(ax.transData, x=1, y=-7, units="points", fig=ax.figure), + ) + + # Create legend + mapElements = [ + Rectangle( + (0, 0), + 1, + 1, + fill=False, + linestyle="--", + color="black", + label="Clipped DEM Extent", + ) + ] + if reportType == "basic": + mapElements.append(Patch(facecolor="black", alpha=1.0, label="Release Areas")) + else: # optional + mapElements.extend( + [ + Patch(facecolor="black", alpha=0.3, label="Entrainment Areas"), + Patch( + facecolor="none", + alpha=0.3, + hatch="xxxx", + label="Resistance Areas", + edgecolor="black", + linewidth=0.5, + ), + ] + ) + + plt.legend(handles=mapElements, title="Legend", loc="center left", bbox_to_anchor=(1, 0.5)) + + # Add DEM colorbar + cax = ax.inset_axes([1.015, 0, 0.375, 0.02]) # [x, y, width, height] + plt.colorbar(im, cax=cax, orientation="horizontal", label="Elevation [m]") + + # Format plot; add title and labels + ax.set_aspect("equal") + ax.xaxis.set_major_locator(mpl.ticker.MaxNLocator(5)) + ax.yaxis.set_major_locator(mpl.ticker.MaxNLocator(5)) + ax.yaxis.set_major_formatter(mpl.ticker.ScalarFormatter(useOffset=False)) + regionalDir = inputDEM.parent.parent.name + reportName = "Basic" if reportType == "basic" else "Optional" + plt.title(f"Split Inputs Report - {reportName} - {regionalDir}") + plt.xlabel("X Coordinate") + plt.ylabel("Y Coordinate") + plt.grid(True, linestyle="--", alpha=0.3) + + # Set axis limits based on DEM extents with a small margin + xMins = [ext[0] for ext in groupExtents.values()] + xMaxs = [ext[1] for ext in groupExtents.values()] + yMins = [ext[2] for ext in groupExtents.values()] + yMaxs = [ext[3] for ext in groupExtents.values()] + xMin, xMax = min(xMins), max(xMaxs) + yMin, yMax = min(yMins), max(yMaxs) + margin = 0.01 + dx = (xMax - xMin) * margin + dy = (yMax - yMin) * margin + ax.set_xlim(xMin - dx, xMax + dx) + ax.set_ylim(yMin - dy, yMax + dy) + + reportPath = outputDir / f"splitInputs_visualReport_{reportType}.png" + plt.savefig(reportPath, dpi=200, bbox_inches="tight") + plt.close() + + return reportPath + + +def getExteriorCoords(geom): + """Get the exterior coordinates of a shapely geometry to handle both single and multi-polygon geometries. + + Parameters + ---------- + geom : shapely.geometry + The shapely geometry to get the exterior coordinates from. + + Returns + ------- + list + A list of tuples containing the x and y coordinates of the geometry exterior. + """ + if isinstance(geom, MultiPolygon): + return [poly.exterior.xy for poly in geom.geoms] + else: + return [geom.exterior.xy] diff --git a/avaframe/runCom7Regional.py b/avaframe/runCom7Regional.py new file mode 100644 index 000000000..56f13f942 --- /dev/null +++ b/avaframe/runCom7Regional.py @@ -0,0 +1,67 @@ +""" +Run script for running simulations in parallel based on input from a regional folder containing +multiple avaFolders +""" + +import time +import argparse + +from avaframe.com7Regional import com7Regional as com7 +from avaframe.in3Utils import cfgUtils +from avaframe.in3Utils import logUtils + + +def runCom7Regional(avalancheDir=""): + """Run regional avalanche simulations in parallel. + + Parameters + ---------- + avalancheDir : str, optional + Path to the main avalanche directory. If not provided, uses the path from general configuration. + + Returns + ------- + allPeakFilesDir : pathlib.Path or None + Path to the directory containing the collected peak files from all sub-avalanche directories, if enabled + allTimeStepsDir : pathlib.Path or None + Path to the directory containing consolidated time step files, if enabled + mergedRastersDir : pathlib.Path or None + Path to the directory containing merged output rasters, if enabled + """ + + # Time the whole routine + startTime = time.time() + + # Load the avalanche directory from command input or the general configuration file + cfgMain = cfgUtils.getGeneralConfig() + if avalancheDir != "": + cfgMain["MAIN"]["avalancheDir"] = avalancheDir + else: + avalancheDir = cfgMain["MAIN"]["avalancheDir"] + + # Start logging + log = logUtils.initiateLogger(str(avalancheDir), logName="runCom7Regional") + log.info("MAIN SCRIPT") + + # Load module configuration + cfg = cfgUtils.getModuleConfig(com7, fileOverride="", toPrint=False, onlyDefault=False) + + # Call main function + allPeakFilesDir, mergedRastersDir = com7.com7RegionalMain(cfgMain, cfg) + + # Print time needed + endTime = time.time() + log.info("Regional process finished after %.1f seconds." % (endTime - startTime)) + + return allPeakFilesDir, mergedRastersDir + + +if __name__ == "__main__": + + parser = argparse.ArgumentParser(description="Run regional workflow") + parser.add_argument( + "avadir", metavar="avalancheDir", type=str, nargs="?", default="", help="the avalanche directory" + ) + + args = parser.parse_args() + runCom7Regional(str(args.avadir)) diff --git a/avaframe/runScripts/runSplitInputs.py b/avaframe/runScripts/runSplitInputs.py new file mode 100644 index 000000000..076567077 --- /dev/null +++ b/avaframe/runScripts/runSplitInputs.py @@ -0,0 +1,57 @@ +"""runScript for splitting avalanche input data into multiple folders based on individual release areas""" + +import time +import pathlib +import argparse + +from avaframe.in3Utils import cfgUtils, logUtils +from avaframe.com7Regional import splitInputs + +def runSplitInputs(avalancheDir=''): + """Main function to split input data for avalanche scenarios. + + Parameters + ---------- + avalancheDir : str + Path to the avalanche directory. If not provided, the function will + retrieve it from the general configuration file. + + Returns + ------- + None + + """ + # Time the whole routine + startTime = time.time() + + # Load the avalanche directory from general configuration file + cfgMain = cfgUtils.getGeneralConfig() + if avalancheDir != '': + cfgMain['MAIN']['avalancheDir'] = avalancheDir + else: + avalancheDir = cfgMain['MAIN']['avalancheDir'] + + # Start logging + log = logUtils.initiateLogger(avalancheDir, logName='runSplitInputs') + log.info('MAIN SCRIPT') + + # Define the output directory + outputDir = pathlib.Path(avalancheDir) / 'com7Regional' + + # Load module configuration + cfg = cfgUtils.getModuleConfig(splitInputs) + + # Run splitting process + splitInputs.splitInputsMain(pathlib.Path(avalancheDir), outputDir, cfg, cfgMain) + + # Print time needed + endTime = time.time() + log.info(f"Completed splitting input data after "f"{endTime - startTime:.1f} seconds.") + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description="Run split inputs for avalanche directories") + parser.add_argument('avalancheDir', type=str, nargs='?', default='', + help="Directory containing the main avalanche data") + + args = parser.parse_args() + runSplitInputs(args.avalancheDir) diff --git a/avaframe/tests/data/testIn4Region/Inputs/ENT/entrainment.dbf b/avaframe/tests/data/testIn4Region/Inputs/ENT/entrainment.dbf new file mode 100644 index 000000000..0fe285f26 Binary files /dev/null and b/avaframe/tests/data/testIn4Region/Inputs/ENT/entrainment.dbf differ diff --git a/avaframe/tests/data/testIn4Region/Inputs/ENT/entrainment.prj b/avaframe/tests/data/testIn4Region/Inputs/ENT/entrainment.prj new file mode 100644 index 000000000..0e46868f1 --- /dev/null +++ b/avaframe/tests/data/testIn4Region/Inputs/ENT/entrainment.prj @@ -0,0 +1 @@ +PROJCS["WGS_1984_UTM_Zone_32N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]] \ No newline at end of file diff --git a/avaframe/tests/data/testIn4Region/Inputs/ENT/entrainment.shp b/avaframe/tests/data/testIn4Region/Inputs/ENT/entrainment.shp new file mode 100644 index 000000000..5c4acf024 Binary files /dev/null and b/avaframe/tests/data/testIn4Region/Inputs/ENT/entrainment.shp differ diff --git a/avaframe/tests/data/testIn4Region/Inputs/ENT/entrainment.shx b/avaframe/tests/data/testIn4Region/Inputs/ENT/entrainment.shx new file mode 100644 index 000000000..b58ddc8df Binary files /dev/null and b/avaframe/tests/data/testIn4Region/Inputs/ENT/entrainment.shx differ diff --git a/avaframe/tests/data/testIn4Region/Inputs/REL/release_areas.dbf b/avaframe/tests/data/testIn4Region/Inputs/REL/release_areas.dbf new file mode 100644 index 000000000..ced69ef0e Binary files /dev/null and b/avaframe/tests/data/testIn4Region/Inputs/REL/release_areas.dbf differ diff --git a/avaframe/tests/data/testIn4Region/Inputs/REL/release_areas.prj b/avaframe/tests/data/testIn4Region/Inputs/REL/release_areas.prj new file mode 100644 index 000000000..0e46868f1 --- /dev/null +++ b/avaframe/tests/data/testIn4Region/Inputs/REL/release_areas.prj @@ -0,0 +1 @@ +PROJCS["WGS_1984_UTM_Zone_32N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]] \ No newline at end of file diff --git a/avaframe/tests/data/testIn4Region/Inputs/REL/release_areas.shp b/avaframe/tests/data/testIn4Region/Inputs/REL/release_areas.shp new file mode 100644 index 000000000..f7e2c6764 Binary files /dev/null and b/avaframe/tests/data/testIn4Region/Inputs/REL/release_areas.shp differ diff --git a/avaframe/tests/data/testIn4Region/Inputs/REL/release_areas.shx b/avaframe/tests/data/testIn4Region/Inputs/REL/release_areas.shx new file mode 100644 index 000000000..deb7e0636 Binary files /dev/null and b/avaframe/tests/data/testIn4Region/Inputs/REL/release_areas.shx differ diff --git a/avaframe/tests/data/testIn4Region/Inputs/RES/resistance.dbf b/avaframe/tests/data/testIn4Region/Inputs/RES/resistance.dbf new file mode 100644 index 000000000..8dbbca9cb Binary files /dev/null and b/avaframe/tests/data/testIn4Region/Inputs/RES/resistance.dbf differ diff --git a/avaframe/tests/data/testIn4Region/Inputs/RES/resistance.prj b/avaframe/tests/data/testIn4Region/Inputs/RES/resistance.prj new file mode 100644 index 000000000..0e46868f1 --- /dev/null +++ b/avaframe/tests/data/testIn4Region/Inputs/RES/resistance.prj @@ -0,0 +1 @@ +PROJCS["WGS_1984_UTM_Zone_32N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]] \ No newline at end of file diff --git a/avaframe/tests/data/testIn4Region/Inputs/RES/resistance.shp b/avaframe/tests/data/testIn4Region/Inputs/RES/resistance.shp new file mode 100644 index 000000000..e89dab088 Binary files /dev/null and b/avaframe/tests/data/testIn4Region/Inputs/RES/resistance.shp differ diff --git a/avaframe/tests/data/testIn4Region/Inputs/RES/resistance.shx b/avaframe/tests/data/testIn4Region/Inputs/RES/resistance.shx new file mode 100644 index 000000000..b1a7d30cb Binary files /dev/null and b/avaframe/tests/data/testIn4Region/Inputs/RES/resistance.shx differ diff --git a/avaframe/tests/data/testIn4Region/Inputs/dem.tif b/avaframe/tests/data/testIn4Region/Inputs/dem.tif new file mode 100644 index 000000000..c6f8e2d75 Binary files /dev/null and b/avaframe/tests/data/testIn4Region/Inputs/dem.tif differ diff --git a/avaframe/tests/test_com7Regional.py b/avaframe/tests/test_com7Regional.py new file mode 100644 index 000000000..c9940b3c3 --- /dev/null +++ b/avaframe/tests/test_com7Regional.py @@ -0,0 +1,163 @@ +"""Test functions for com7Regional module.""" + +import pathlib +import configparser +import shutil +from avaframe.com7Regional import splitInputs +from avaframe.com7Regional import com7Regional +from avaframe.in3Utils import generateTopo, getReleaseArea, cfgUtils, initializeProject + + +def setup_test_directory(tmp_path, avaName): + """Set up test directory with required files.""" + avaDir = tmp_path / "RegionalAvalanches" / avaName + initializeProject.createFolderStruct(avaDir) + + # Create topography + cfgTopo = cfgUtils.getModuleConfig(generateTopo) + cfgTopo["TOPO"]["dx"] = "5" + cfgTopo["TOPO"]["xEnd"] = "500" + cfgTopo["TOPO"]["yEnd"] = "500" + cfgTopo["TOPO"]["meanAlpha"] = "30" + cfgTopo["TOPO"]["demType"] = "FP" + cfgTopo["TOPO"]["z0"] = "500" + cfgTopo["DEMDATA"] = { + "xl": "0", + "yl": "0", + "zEdit": "", + "nodata_value": "-9999", + "demName": "ava_topo", + } + + generateTopo.generateTopo(cfgTopo, avaDir) + + # Create release area + cfgRelease = cfgUtils.getModuleConfig(getReleaseArea) + cfgRelease["FILE"]["relNo"] = "1" + cfgRelease["FILE"]["relName"] = "release1" + cfgRelease["GENERAL"]["dh"] = "1.0" + cfgRelease["FILE"]["outputtxt"] = "False" + cfgRelease["GEOMETRY"] = {"x0": "25", "y0": "25", "width": "10", "length": "10"} + + getReleaseArea.getReleaseArea(cfgTopo, cfgRelease, avaDir) + + return avaDir + + +def create_test_configs(): + """Create test configurations.""" + # Create main config + cfgMainDict = {"MAIN": {"avalancheDir": "", "nCPU": "2", "flagDev": "True"}} + + # Create com7Regional config + cfgDict = { + "GENERAL": { + "regionalDir": "RegionalAvalanches", + "copyPeakFiles": "True", + "moveInsteadOfCopy": "False", + "mergeOutput": "True", + "mergeTypes": "pfv", + "mergeMethods": "max|count", + }, + "com1DFA_com1DFA_override": { + "defaultConfig": "False", + "dt": "0.5", + "simTypeList": "null", + "resType": "ppr|pft|pfv", + "tSteps": "1", + "tEnd": "2", + "meshCellSizeThreshold": "0.001", + }, + "INPUT": {"releaseScenario": "release1"}, + "com1DFA": {"simTypeList": "null", "releaseScenario": "release1"}, + "FLAGS": {"debugPlot": "False"}, + } + + cfgMain = cfgUtils.convertDictToConfigParser(cfgMainDict) + cfg = cfgUtils.convertDictToConfigParser(cfgDict) + + return cfgMain, cfg + + +def test_com7RegionalMain(tmp_path): + """Test the main function of com7Regional module.""" + # Create avadirs with valid inputs + avaNames = ["avaTest1", "avaTest2", "avaTest3"] + avaDirs = [] + for avaName in avaNames: + avaDir = setup_test_directory(tmp_path, avaName) + avaDirs.append(avaDir) + + cfgMain, cfg = create_test_configs() + cfgMain["MAIN"]["avalancheDir"] = str(tmp_path) + + allPeakFilesDir, mergedRastersDir = com7Regional.com7RegionalMain(cfgMain, cfg) + + # Check for expected output files and directories + assert allPeakFilesDir is not None + assert allPeakFilesDir.exists() + assert mergedRastersDir is not None + assert mergedRastersDir.exists() + + for avaDir in avaDirs: + outputDir = avaDir / "Outputs" / "com1DFA" + assert outputDir.exists() + assert (outputDir / "peakFiles").exists() + + +def test_getTotalNumberOfSims(tmp_path): + """Test the getTotalNumberOfSims function.""" + avaDir = setup_test_directory(tmp_path, "testAva") + + cfgMain, cfg = create_test_configs() + cfgMain["MAIN"]["avalancheDir"] = str(tmp_path) + + totalSims = com7Regional.getTotalNumberOfSims([avaDir], cfgMain, cfg) + + assert totalSims > 0 + + +def test_splitInputsMain(tmp_path): + """Test splitInputsMain function using pre-generated test data.""" + # Set up test data + test_data_dir = pathlib.Path(__file__).parent / "data" / "testIn4Region" + inputDir = tmp_path / "avalancheDir" + shutil.copytree(test_data_dir, inputDir) + outputDir = tmp_path / "avalancheDir" / "com7Regional" + + # Configure test parameters + cfg = configparser.ConfigParser() + cfg["GENERAL"] = {"bufferSize": "2500"} + cfgMain = configparser.ConfigParser() + cfgMain["FLAGS"] = {"createReport": "True", "savePlot": "True"} + + # Run function + splitInputs.splitInputsMain(inputDir, outputDir, cfg, cfgMain) + + # Verify outputs + assert outputDir.exists() + + # Check group directories + groupDirs = list(outputDir.glob("group*")) + assert len(groupDirs) == 2 + + for groupDir in groupDirs: + # Check directory structure + assert (groupDir / "Inputs").exists() + assert (groupDir / "Inputs" / "REL").exists() + assert len(list((groupDir / "Inputs").glob("*.tif"))) == 1 + + # Check release areas were split by scenarios + relDir = groupDir / "Inputs" / "REL" + assert len(list(relDir.glob("*.shp"))) == 2 + + # Check optional inputs + assert (groupDir / "Inputs" / "ENT").exists() + assert len(list((groupDir / "Inputs" / "ENT").glob("*.shp"))) == 1 + assert (groupDir / "Inputs" / "RES").exists() + assert len(list((groupDir / "Inputs" / "RES").glob("*.shp"))) == 1 + + # Check reports + assert (outputDir / "splitInputs_scenarioReport.txt").exists() + assert (outputDir / "splitInputs_visualReport_basic.png").exists() + assert (outputDir / "splitInputs_visualReport_optional.png").exists() diff --git a/docs/_static/splitInputs_visualReport_basic.png b/docs/_static/splitInputs_visualReport_basic.png new file mode 100644 index 000000000..84ad3e252 Binary files /dev/null and b/docs/_static/splitInputs_visualReport_basic.png differ diff --git a/docs/_static/splitInputs_visualReport_optional.png b/docs/_static/splitInputs_visualReport_optional.png new file mode 100644 index 000000000..22942def5 Binary files /dev/null and b/docs/_static/splitInputs_visualReport_optional.png differ diff --git a/docs/api.rst b/docs/api.rst index 33d60005c..b0af8b82c 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -16,6 +16,7 @@ Computational Modules com4FlowPy com5SnowSlide com6RockAvalanche + com7Regional Input/Transformation Modules ============================ diff --git a/docs/dataSources.rst b/docs/dataSources.rst index 4332f099b..aa0672b4b 100644 --- a/docs/dataSources.rst +++ b/docs/dataSources.rst @@ -13,5 +13,5 @@ Here is a list of external data we used, with links and sources if available: **Datenquelle: Land Tirol - data.tirol.gv.at** - https://www.data.gv.at/katalog/dataset/land-tirol_tirolgelnde + `Open data Austria Land Tirol `_ diff --git a/docs/index.rst b/docs/index.rst index 08b4d11ed..1a8df07a0 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -52,6 +52,7 @@ Computational modules * :doc:`moduleCom4FlowPy` * :doc:`moduleCom5SnowSlide` * :doc:`moduleCom6RockAvalanche` + * :doc:`moduleCom7Regional` * :doc:`moduleCom8MoTPSA` * :doc:`moduleCom8MoTVoellmy` @@ -66,6 +67,7 @@ Computational modules moduleCom4FlowPy.rst moduleCom5SnowSlide.rst moduleCom6RockAvalanche.rst + moduleCom7Regional.rst moduleCom8MoTPSA.rst moduleCom9MoTVoellmy.rst diff --git a/docs/moduleCom7Regional.rst b/docs/moduleCom7Regional.rst new file mode 100644 index 000000000..73defbaf2 --- /dev/null +++ b/docs/moduleCom7Regional.rst @@ -0,0 +1,320 @@ +###################################### +com7Regional: Regional Modeling +###################################### + +.. note:: + This module is still under development and may contain bugs or incomplete features. + +:py:mod:`com7Regional` is an experimental module targeted at the application of the :ref:`com1DFA ` +kernel on a regional scale. +In essence, it is a wrapper module that allows for the concurrent execution of com1DFA within multiple avalanche directories, +with the aim of reducing the overall computational load through parallelization. Com1DFA is designed with one single avalanche track as +the basis for computation. It includes parallel computation of different scenarios, but always on one, reasonably sized DEM. +Com7Regional allows for the splitting of a setup with multiple avalanche tracks, for example one whole valley, into more managable +setups and then the subsequent results aggregation of scenarios or regions. + +The module therefor provides functions for: + +* preparing and splitting input files from a master DEM and master release shapefile +* managing and aggregating outputs after a regional simulation, e.g. merge com1DFA peak files to rasters of maximum / minimum values. + +The log output gives information about the processing status of each individual avalanche directory, as well as a summary +of how many simulations were processed successfully. + +.. note:: + Experience with :ref:`com1DFA ` is recommended before using this module. + +Regional Input Management +===================================== +:py:mod:`com7Regional.splitInputs` is a module for organizing large amounts of avalanche input data into multiple avalanche directories, based +on the AvaFrame data structure, thus enabling efficient and controlled input management of larger datasets for further processing. +Importantly, the module also provides a simple, automatic method for clipping the input DEM around each release area group, +which then forms the basis for the organization of other input. + +Input +----- +The module is currently compatible with the following input file types: + +**Required:** + * Digital elevation model (as .asc or .tif file). *Please note that the maximum possible size of this DEM depends on your + compute hardware; You might run into problems if it is too big and maxes out your resources*. + * Release areas: ONE shapefile in ``Inputs/REL`` directory, with additional attributes ``group`` and ``scenario``. Please + also provide attributes required by com1DFA. + +**Optional:** + * ONE file with entrainment areas (shapefile in ``Inputs/ENT`` directory) + * ONE file with resistance areas (shapefile in ``Inputs/RES`` directory) + + +Where the expected input directory structure is:: + + avalancheDir/ + └── Inputs/ + ├── REL/ + │ └── *.shp # release areas (one file) + ├── ENT/ # entrainment areas (one file; optional) + │ └── *.shp + ├── RES/ # resistance areas (one file; optional) + │ └── *.shp + └── *.asc or *.tif # digital elevation model (DEM) + +Group and scenario creation +--------------------------- +Input data organization is based on two key concepts: + +* **Groups**: Collections of avalanche release areas (single polygon features) that are located in the same spatial domain and may be + wanted to be simulated together. We recommend setting one group per avalanche path/track. + +* **Scenarios**: Collections of release area features WITHIN each group, that are simulated together in :ref:`com1DFA `, + i.e. particles within one scenario can interact. + +These are defined through two new attributes in the input release area shapefile attribute table: + +* ``group``: expected format: string (e.g. *"group1"*) + +* ``scenario``: expected format: comma-separated list of strings, without spaces (e.g. *"small,large"* or *"10y,30y,100y"*) + +For example: + +.. list-table:: + :header-rows: 1 + :widths: 40 40 40 + + * - name + - group + - scenario + * - rel1 + - avaPath1 + - 1,all + * - rel2 + - avaPath1 + - 2,all + * - rel3 + - avaPath2 + - 1 + +In this example, four release scenarios (shapefiles) would be created, in two separate sub-avalanche directories:: + + avalancheDir/ + ├── ... + └── com7Regional/ + ├── avaPath1/ + │ └── REL/ + │ ├── avaPath1_1.shp - containing [1] feature: rel1 + │ ├── avaPath1_2.shp - containing [1] feature: rel2 + │ └── avaPath1_all.shp - containing [2] features: rel1, rel2 + └── avaPath2/ + └── REL/ + └── avaPath2_1.shp - containing [1] feature: rel3 + +In the case that scenarios are defined for only some release features within a group, the rest will be grouped together as a single 'NULL' scenario. + +In the case that no attributes or values for ``group`` or ``scenario`` are provided, the procedure will create groups based on each release feature +- so one avalanche directory, including a DEM and release scenario shp file, per release feature (polygon). This may be +intended to simulate each release feature separately in an automated manner. However, keep in mind, that this approach may result in +a large amount of duplicated DEM data, in the case that release areas are located in close proximity to each other. + +Output +------ +Running ``runScripts/runSplitInputs.py`` with valid input data will result in the following output in ``/com7Regional``: + +1. Individual avalanche directories for each group containing: + + - Clipped DEM file + - Scenario-specific release area shapefiles + - Optional: Clipped entrainment and resistance areas + +2. Two visual reports (see :numref:`fig-splitInputs-basic` and :numref:`fig-splitInputs-optional`) + +3. Scenario report in txt format (see example below) + +.. list-table:: + :widths: 50 50 + + * - .. _fig-splitInputs-basic: + .. figure:: /_static/splitInputs_visualReport_basic.png + :width: 100% + :alt: Basic visual report + + Example of basic inputs report displaying resulting groups and their extent + - .. _fig-splitInputs-optional: + .. figure:: /_static/splitInputs_visualReport_optional.png + :width: 100% + :alt: Optional inputs report + + Example of optional inputs report displaying RES and ENT areas for each group + +Example Scenario Report: + +.. code-block:: text + + SCENARIO REPORT + ============== + Generated: 2025-02-04 10:58:56 + + Group: group1 + ------------ + + Scenario: sce1 + No. of release areas: 2 + - rel1 + - rel2 + + No. of entrainment areas: 1 + No. of resistance areas: 3 + + Group: group2 + ------------ + ... + +Configuration +------------- +Settings are controlled through ``splitInputsCfg.ini``, in which the ``bufferSize`` for the group extent is defined (which is used for DEM, RES, and ENT clipping into +smaller chunks). By default, this value is set to 2500 m. For each group, a bounding box is created from the maximum x-y extent of all release features in the group. +The value for ``bufferSize`` is then added to each direction (+x, -x, +y, -y). This buffer may be adjusted according to the expected maximum runout length of your avalanches - +a larger value will ensure that no simulation will exit its domain, but with the drawback of producing larger output files and longer run times. + +Algorithm +--------- +The ``splitInputsMain`` function, which is called in ``runScripts/runSplitInputs.py``, performs the following steps: + +1. Create central avalanche directory list +2. Set up avalanche directories +3. Split and write release areas to each directory +4. Clip and write DEM to each directory +5. Clip and write optional input to each directory (currently includes RES and ENT) +6. Divide release areas into scenarios +7. Write reports + +To Run - Script based +--------------------- +1. Prepare inputs in your ``/Inputs`` +2. If you want to adjust settings, copy ``splitInputsCfg.ini`` to ``local_splitInputsCfg.ini`` and adjust values in there +3. Either set path to avalanche directory in ``avaframeCfg.ini`` (or local version ``local_avaframeCfg.ini``) or + call command below with the avalanche directory as argument +4. Execute from ``AvaFrame/avaframe`` directory: + +.. code-block:: bash + + python runScripts/runSplitInputs.py + +--------------- + +Running multiple avalanche directories +====================================== + +Input +----- +A directory structure containing pre-configured avalanche directories (containing an ``Inputs`` folder) is required. For input preparation use +:ref:`moduleCom7Regional:Regional Input Management`, which splits merged input data into standard :ref:`com1DFA ` +inputs across multiple avalanche directories. + +Example of a valid directory structure (as produced by the regional input management above):: + + avalancheDir + ├── Inputs/ #NOT being used for running; optional + └── com7Regional/ #This is the default name, can be changed via .ini setup + ├── sub_avalanche1/ + │ └── Inputs/ + │ ├── REL/*.shp + │ └── *.asc or *.tif + ├── sub_avalanche2/ + │ └── Inputs/ + │ ├── REL/*.shp + │ └── *.asc or *.tif + └── ... + + +Output +------ +Outputs are organized in two levels: + +**1. Merged rasters** and **2. Individual outputs (per sub_avalanche directory)** + +Merged rasters +^^^^^^^^^^^^^^ +Configure in ``com7RegionalCfg.ini`` (or local): + +.. code-block:: ini + + [GENERAL] + mergeOutput = True + mergeTypes = pfv # Available options: [ppr|pfv|pft|pta|FT|FV|P|FM|Vx|Vy|Vz|TA] + mergeMethod = max # Available options: [max|min|sum|count] + +Produces merged rasters of all peakFile results found within the sub-avalanche directories, for each +``mergeTypes`` and ``mergeMethod`` configured, in ``/com7Regional/mergedRasters/``. The merged +raster combines output from ALL sub-avalanches. + +Creates:: + + avalancheDir + ├── .... + └── com7Regional/ + ├── sub_avalanche1/ + ... + └── mergedRasters/ <- this one is created + +Individual outputs +^^^^^^^^^^^^^^^^^^ +After running com7 with the given module (currently only :ref:`com1DFA `), the standard output is located +within each of the sub-avalanche directories within e.g. ``/com7Regional//Outputs/com1DFA``. +Additionally, com7Regional provides the option of aggregating all output peakFiles and tSteps results into a single directory +for easier management, either through copying or moving the files after an executed run. + +Configure in ``com7RegionalCfg.ini`` (or local): + +.. code-block:: ini + + [GENERAL] + copyPeakFiles = True + moveInsteadOfCopy = False + +Creates:: + + avalancheDir + ├── .... + └── com7Regional/ + ├── sub_avalanche1/ + ... + └── allPeakFiles/ <- this one is created + +Configuration +------------- +Three configuration files are used: + +1. Main configuration (``avaframeCfg.ini``) + - To set nCPUs for handling the amount of avalanche directories processed in parallel + - Handle plot and other output generation + +2. com7Regional configuration (``com7RegionalCfg.ini``) + - Manages output aggregation and merged raster creation + - Overrides com1DFA parameters if specified + +3. com1DFA configuration (``com1DFACfg.ini``) + - Standard simulation parameters + +Processing +---------- +Parallelization is handled through the concurrent.futures library, specifically the +`ProcessPoolExecutor class `_. +In essence, tasks are executed concurrently within the input regional directory (by default com7Regional), based on the number of currently available CPUs. +So each sub-avalanche directory is one task. +The maximum number of CPUs is set by the ``nCPU`` parameter in ``avaframeCfg.ini``. By default, to avoid nested parallelization, each +avalanche directory is assigned a single CPU , essentially meaning that any variations (e.g. through different +scenarios, parameter variations, etc.) within each avalanche directory are handled sequentially. As a consequence, if the number of +variations is high, and the number of avaDirs to process is lower than ``nCPU``, it may be more efficient to run simulations with the +standard ``runCom1DFA.py`` instead, to utilize its parallel processing of variations. Alternatively, advanced users may want to adjust the nCPU for +variations in ``com7Regional.py``. + +To Run +------ +1. Prepare input directories, we recommend using the regional input management above +2. If you want to adjust settings, copy ``com7RegionalCfg.ini`` to ``local_com7RegionalCfg.ini`` and adjust values in there +3. Set path to avalanche directory in ``avaframeCfg.ini`` (or local version ``local_avaframeCfg.ini``) + or supply the directory as argument to the command below. +4. Execute from AvaFrame/avaframe directory: + +.. code-block:: bash + + python runCom7Regional.py diff --git a/tinyHelper/extractDEM.sh b/tinyHelper/extractDEM.sh index 893262d01..803bed33e 100755 --- a/tinyHelper/extractDEM.sh +++ b/tinyHelper/extractDEM.sh @@ -2,7 +2,7 @@ # Steps to take to generate a new topo manually (instead of script below): # 1. go to https://www.data.gv.at/katalog/dataset/land-tirol_tirolgelnde#resources and get DGM 5m Tirol -# 2. Load zip into QGis, make sure the correct projection is used +# 2. Extract zip and load tif into QGis, make sure the correct projection is used # 3. Use Raster - Warp to reproject to epsg31287 (use resampling method: cubic) -> save as tiff # 4. Draw the extend in a shapefile with epsg31287 # 5. Use Raster - Clip raster by mask layer to cut out the new DGM -> set .asc as output right away