diff --git a/avaframe/ana3AIMEC/aimecTools.py b/avaframe/ana3AIMEC/aimecTools.py index 9aac43451..93520b9be 100644 --- a/avaframe/ana3AIMEC/aimecTools.py +++ b/avaframe/ana3AIMEC/aimecTools.py @@ -1,9 +1,10 @@ """ - Main logic for AIMEC post processing +Main logic for AIMEC post processing """ import logging import pathlib +import re import numpy as np import pandas as pd import copy @@ -30,8 +31,8 @@ # ----------------------------------------------------------- -def readAIMECinputs(avalancheDir, pathDict, defineRunoutArea, dirName='com1DFA'): - """ Read inputs for AIMEC postprocessing +def readAIMECinputs(avalancheDir, pathDict, defineRunoutArea, dirName="com1DFA"): + """Read inputs for AIMEC postprocessing Reads the required geometry files location for AIMEC postpocessing given an avalanche directory; avalanche path, DEM @@ -54,53 +55,60 @@ def readAIMECinputs(avalancheDir, pathDict, defineRunoutArea, dirName='com1DFA') updated dictionary with path to geometry input data """ - refDir = pathlib.Path(avalancheDir, 'Inputs', 'LINES') - profileLayer = list(refDir.glob('*aimec*.shp')) + refDir = pathlib.Path(avalancheDir, "Inputs", "LINES") + profileLayer = list(refDir.glob("*aimec*.shp")) try: - message = ('There should be exactly one path_aimec.shp file containing the avalanche path in %s/Inputs/LINES/' % - avalancheDir) + message = ( + "There should be exactly one path_aimec.shp file containing the avalanche path in %s/Inputs/LINES/" + % avalancheDir + ) assert len(profileLayer) == 1, message except AssertionError: raise - pathDict['profileLayer'] = profileLayer[0] + pathDict["profileLayer"] = profileLayer[0] if defineRunoutArea: - refDir = pathlib.Path(avalancheDir, 'Inputs', 'POINTS') - splitPointLayer = list(refDir.glob('*.shp')) + refDir = pathlib.Path(avalancheDir, "Inputs", "POINTS") + splitPointLayer = list(refDir.glob("*.shp")) try: - message = 'There should be exactly one .shp file containing the split point in %s/Inputs/POINTS/' % avalancheDir + message = ( + "There should be exactly one .shp file containing the split point in %s/Inputs/POINTS/" + % avalancheDir + ) assert len(splitPointLayer) == 1, message except AssertionError: raise - pathDict['splitPointSource'] = splitPointLayer[0] + pathDict["splitPointSource"] = splitPointLayer[0] else: - pathDict['splitPointSource'] = None + pathDict["splitPointSource"] = None - refDir = pathlib.Path(avalancheDir, 'Inputs') + refDir = pathlib.Path(avalancheDir, "Inputs") # check for DEM - if 'demFileName' not in pathDict.keys(): - demSource = list(refDir.glob('*.asc')) + list(refDir.glob('*.tif')) - elif pathDict['demFileName'] == '': - demSource = list(refDir.glob('*.asc')) + list(refDir.glob('*.tif')) + if "demFileName" not in pathDict.keys(): + demSource = list(refDir.glob("*.asc")) + list(refDir.glob("*.tif")) + elif pathDict["demFileName"] == "": + demSource = list(refDir.glob("*.asc")) + list(refDir.glob("*.tif")) else: - demSource = list(refDir.glob(pathDict['demFileName'])) + demSource = list(refDir.glob(pathDict["demFileName"])) try: - assert len(demSource) == 1, 'There should be exactly one topography .asc file in %s/Inputs/' % avalancheDir + assert len(demSource) == 1, ( + "There should be exactly one topography .asc file in %s/Inputs/" % avalancheDir + ) except AssertionError: raise - pathDict['demSource'] = demSource[0] + pathDict["demSource"] = demSource[0] - if dirName != '': - pathResult = pathlib.Path(avalancheDir, 'Outputs', 'ana3AIMEC', dirName) + if dirName != "": + pathResult = pathlib.Path(avalancheDir, "Outputs", "ana3AIMEC", dirName) else: - pathResult = pathlib.Path(avalancheDir, 'Outputs', 'ana3AIMEC') - pathDict['pathResult'] = pathResult + pathResult = pathlib.Path(avalancheDir, "Outputs", "ana3AIMEC") + pathDict["pathResult"] = pathResult # check for reference data - referenceDir= pathlib.Path(avalancheDir, 'Inputs') + referenceDir = pathlib.Path(avalancheDir, "Inputs") # file names need to contain: _LINE, _POLY, _POINT - referenceTypes = {"referenceLine": "LINE", "referencePolygon": "POLY",'referencePoint': 'POINT'} + referenceTypes = {"referenceLine": "LINE", "referencePolygon": "POLY", "referencePoint": "POINT"} for refType in referenceTypes: referenceFile, availableFile, _ = gI.getAndCheckInputFiles( referenceDir, @@ -109,114 +117,124 @@ def readAIMECinputs(avalancheDir, pathDict, defineRunoutArea, dirName='com1DFA') fileExt="shp", fileSuffix=("_" + referenceTypes[refType]), ) - if availableFile == 'Yes': + if availableFile == "Yes": # add file paths to pathDict pathDict[refType] = [referenceFile] else: pathDict[refType] = [] projectName = pathlib.Path(avalancheDir).stem - pathDict['projectName'] = projectName + pathDict["projectName"] = projectName pathName = profileLayer[0].stem - pathDict['pathName'] = pathName - pathDict['avalancheDir'] = avalancheDir + pathDict["pathName"] = pathName + pathDict["avalancheDir"] = avalancheDir return pathDict -def fetchReferenceSimNo(avaDir, inputsDF, comModule, cfg, inputDir=''): - """ Define reference simulation used for aimec analysis. - - if the configuration files are available and a varParList is provided, the simulations - are ordered and the referenceSimValue is used to define the reference. - - otherwise, if a referenceSimName is provided, the reference is based on this name - - otherwise, the first file found is used as reference - - Parameters - ----------- - avaDir : str - path to avalanche directory - inputsDF: dataFrame - dataFrame with simulations to analyze and path to result files - comModule: str - computational module used to produce the results to analyze - cfg: configParser object - configuration for aimec - referenceSimValue, varParList used here - inputDir: str or pathlib path - optional- directory where peak files are located, if '', - avaDir/Outputs/comMod/peakFiles is set - - Returns - -------- - refSimRowHash: str - dataframe hash (not simHash) of the simulation used as reference - refSimName: str - name of the simulation used as reference - inputsDF:dataFrame - dataFrame with simulations to analyze and path to result files - If com1DFA = comModule and a variation parameter was specified, the - com1DFA configuration is merged to the inputsDF - colorVariation: boolean - True if a color variation should be applied in the plots - valRef: str - value of vapParList[0] used to define reference sim +def fetchReferenceSimNo(avaDir, inputsDF, comModule, cfg, inputDir=""): + """Define reference simulation used for aimec analysis. + + if the configuration files are available and a varParList is provided, the simulations + are ordered and the referenceSimValue is used to define the reference. + + otherwise, if a referenceSimName is provided, the reference is based on this name + + otherwise, the first file found is used as reference + + Parameters + ----------- + avaDir : str + path to avalanche directory + inputsDF: dataFrame + dataFrame with simulations to analyze and path to result files + comModule: str + computational module used to produce the results to analyze + cfg: configParser object + configuration for aimec - referenceSimValue, varParList used here + inputDir: str or pathlib path + optional- directory where peak files are located, if '', + avaDir/Outputs/comMod/peakFiles is set + + Returns + -------- + refSimRowHash: str + dataframe hash (not simHash) of the simulation used as reference + refSimName: str + name of the simulation used as reference + inputsDF:dataFrame + dataFrame with simulations to analyze and path to result files + If com1DFA = comModule and a variation parameter was specified, the + com1DFA configuration is merged to the inputsDF + colorVariation: boolean + True if a color variation should be applied in the plots + valRef: str + value of vapParList[0] used to define reference sim """ - cfgSetup = cfg['AIMECSETUP'] - if inputDir != '': + cfgSetup = cfg["AIMECSETUP"] + if inputDir != "": inputDir = pathlib.Path(inputDir) else: - inputDir = pathlib.Path(avaDir, 'Outputs', comModule, 'peakFiles') + inputDir = pathlib.Path(avaDir, "Outputs", comModule, "peakFiles") if not inputDir.is_dir(): - message = 'Input directory %s does not exist - check anaMod' % inputDir + message = "Input directory %s does not exist - check anaMod" % inputDir log.error(message) raise NotADirectoryError(message) - referenceSimValue = cfgSetup['referenceSimValue'] - referenceSimName = cfgSetup['referenceSimName'] + referenceSimValue = cfgSetup["referenceSimValue"] + referenceSimName = cfgSetup["referenceSimName"] colorVariation = False # look for a configuration try: # load dataFrame for all configurations configurationDF = cfgUtils.createConfigurationInfo(avaDir, comModule=comModule) # Merge inputsDF with the configurationDF. Make sure to keep the indexing from inputs and to merge on 'simName' - inputsDF = inputsDF.reset_index().merge(configurationDF, on=['simName', 'modelType']).set_index('index') + inputsDF = ( + inputsDF.reset_index().merge(configurationDF, on=["simName", "modelType"]).set_index("index") + ) configFound = True except (NotADirectoryError, FileNotFoundError) as e: - if cfgSetup['varParList'] != '' and (any(item in inputsDF.columns.tolist() for item in cfgSetup['varParList'].split('|')) == False): - message = ('No configuration directory found. This is needed for sorting simulation according to ' - '%s' % cfgSetup['varParList']) + if cfgSetup["varParList"] != "" and ( + any(item in inputsDF.columns.tolist() for item in cfgSetup["varParList"].split("|")) == False + ): + message = ( + "No configuration directory found. This is needed for sorting simulation according to " + "%s" % cfgSetup["varParList"] + ) raise (message) - elif cfgSetup['varParList'] != '' and any(item in inputsDF.columns.tolist() for item in cfgSetup['varParList'].split('|')): + elif cfgSetup["varParList"] != "" and any( + item in inputsDF.columns.tolist() for item in cfgSetup["varParList"].split("|") + ): configFound = True else: configFound = False # filter simulations - if cfg.has_section('FILTER'): - parametersDict = fU.getFilterDict(cfg, 'FILTER') - simNameList = cfgHandling.filterSims(avaDir, parametersDict, specDir='') - inputsDF = inputsDF[inputsDF['simName'].isin(simNameList)] - log.info('%d simulations found matching filter criteria' % inputsDF.shape[0]) + if cfg.has_section("FILTER"): + parametersDict = fU.getFilterDict(cfg, "FILTER") + simNameList = cfgHandling.filterSims(avaDir, parametersDict, specDir="") + inputsDF = inputsDF[inputsDF["simName"].isin(simNameList)] + log.info("%d simulations found matching filter criteria" % inputsDF.shape[0]) # set reference value to empty string - valRef = '' + valRef = "" # if the simulations come from com1DFA, it is possible to order the files and define a reference - if configFound and cfgSetup['varParList'] != '': + if configFound and cfgSetup["varParList"] != "": # fetch parameters that shall be used for ordering - varParList = cfgSetup['varParList'].split('|') + varParList = cfgSetup["varParList"].split("|") # order simulations - varParList, inputsDF = cfgHandling.orderSimulations(varParList, cfgSetup.getboolean('ascendingOrder'), inputsDF) + varParList, inputsDF = cfgHandling.orderSimulations( + varParList, cfgSetup.getboolean("ascendingOrder"), inputsDF + ) colorVariation = True # now look for the reference - if referenceSimValue != '': + if referenceSimValue != "": # if a referenceSimValue is provided, find the corresponding simulation # get the value of the first parameter used for ordering (this will be usefull for colorcoding in plots) refSimRowHash, refSimName, valRef = defineRefOnSimValue(referenceSimValue, varParList, inputsDF) - elif cfgSetup['referenceSimName'] != '': + elif cfgSetup["referenceSimName"] != "": # else if a referenceSimName is provided, find the corresponding simulation - set # simulation with referenceSimName in name as referene simulation refSimRowHash, refSimName = defineRefOnSimName(referenceSimName, inputsDF) @@ -224,12 +242,16 @@ def fetchReferenceSimNo(avaDir, inputsDF, comModule, cfg, inputDir=''): # if no referenceSimValue is provided, we assume the first simulation after reordering is the reference # reference simulation refSimRowHash = inputsDF.index[0] - refSimName = inputsDF.loc[refSimRowHash, 'simName'] - log.info(('Reference simulation is the first simulation after reordering according to %s in ascending order' - '= %s and corresponds to simulation %s') - % (varParList, cfgSetup.getboolean('ascendingOrder'), refSimName)) + refSimName = inputsDF.loc[refSimRowHash, "simName"] + log.info( + ( + "Reference simulation is the first simulation after reordering according to %s in ascending order" + "= %s and corresponds to simulation %s" + ) + % (varParList, cfgSetup.getboolean("ascendingOrder"), refSimName) + ) - elif referenceSimName != '': + elif referenceSimName != "": # else if a referenceSimName is provided, find the corresponding simulation - set # simulation with referenceSimName in name as referene simulation refSimRowHash, refSimName = defineRefOnSimName(referenceSimName, inputsDF) @@ -238,16 +260,18 @@ def fetchReferenceSimNo(avaDir, inputsDF, comModule, cfg, inputDir=''): # if no ordering is done, no referenceSimValue nor referenceSimName are given, take the first simulation # as reference refSimRowHash = inputsDF.index[0] - refSimName = inputsDF.loc[refSimRowHash, 'simName'] - message = ('No information on how to define the reference was provided, taking simulation %s as reference.' - % refSimName) + refSimName = inputsDF.loc[refSimRowHash, "simName"] + message = ( + "No information on how to define the reference was provided, taking simulation %s as reference." + % refSimName + ) log.warning(message) return refSimRowHash, refSimName, inputsDF, colorVariation, valRef def defineRefOnSimValue(referenceSimValue, varParList, inputsDF): - """ Search for reference based on configuration parameter and value + """Search for reference based on configuration parameter and value Raise an error if no simulation is found @@ -282,8 +306,8 @@ def defineRefOnSimValue(referenceSimValue, varParList, inputsDF): valRef = sortingParameter[indexRef] elif typeCP in [float, int]: # check if thicknessValue read from shp - if np.isnan(sortingParameter).any() and varParList[0] in ['relTh', 'entTh', 'secondaryRelTh']: - sortingParameter = inputsDF[(varParList[0] + '0')].to_list() + sortingParameter + if np.isnan(sortingParameter).any() and varParList[0] in ["relTh", "entTh", "secondaryRelTh"]: + sortingParameter = inputsDF[(varParList[0] + "0")].to_list() + sortingParameter colorValues = np.asarray(sortingParameter) # look for closest value indexRef = np.nanargmin(np.abs(colorValues - typeCP(referenceSimValue))) @@ -294,18 +318,22 @@ def defineRefOnSimValue(referenceSimValue, varParList, inputsDF): # there might be multiple simulations matching the referenceSimValue, we take the first one refSimRowHash = inputsDF[inputsDF[varParList[0]] == valRef].index refSimRowHash, refSimName = checkMultipleSimFound(refSimRowHash, inputsDF) - log.info(('Reference simulation is based on %s = %s - closest value ' - 'found is: %s and corresponds to simulation %s') - % (varParList[0], referenceSimValue, str(valRef), refSimName)) + log.info( + ( + "Reference simulation is based on %s = %s - closest value " + "found is: %s and corresponds to simulation %s" + ) + % (varParList[0], referenceSimValue, str(valRef), refSimName) + ) except ValueError: - message = 'Did not find any simulation matching %s = %s.' % (varParList[0], referenceSimValue) + message = "Did not find any simulation matching %s = %s." % (varParList[0], referenceSimValue) log.error(message) raise ValueError(message) return refSimRowHash, refSimName, str(valRef) def defineRefOnSimName(referenceSimName, inputsDF): - """ Search for reference based on simulation name + """Search for reference based on simulation name Raise an error if no simulation is found @@ -323,20 +351,21 @@ def defineRefOnSimName(referenceSimName, inputsDF): refSimName : str name of THE reference simulation """ - refSimRowHash = inputsDF.index[inputsDF['simName'].str.contains('%s' % referenceSimName).to_list()] + refSimRowHash = inputsDF.index[inputsDF["simName"].str.contains("%s" % referenceSimName).to_list()] if len(refSimRowHash) == 0: - message = ('Found no simulation matching the provided referenceSimName = %s.' - % referenceSimName) + message = "Found no simulation matching the provided referenceSimName = %s." % referenceSimName log.error(message) raise IndexError(message) refSimRowHash, refSimName = checkMultipleSimFound(refSimRowHash, inputsDF) - log.info('Reference simulation is based on provided simName: %s and corresponds to simulation %s' - % (referenceSimName, refSimName)) + log.info( + "Reference simulation is based on provided simName: %s and corresponds to simulation %s" + % (referenceSimName, refSimName) + ) return refSimRowHash, refSimName def checkMultipleSimFound(refSimRowHash, inputsDF, error=False): - """ Check if more than one file can be chosen as reference + """Check if more than one file can be chosen as reference Log warning or error if it is the case @@ -358,90 +387,214 @@ def checkMultipleSimFound(refSimRowHash, inputsDF, error=False): """ if len(refSimRowHash) > 1: if not error: - message = ('Found multiple simulations (%s) matching the reference criterion, ' - 'taking the first one as reference' % len(refSimRowHash)) + message = ( + "Found multiple simulations (%s) matching the reference criterion, " + "taking the first one as reference" % len(refSimRowHash) + ) log.warning(message) else: - message = ('Found multiple simulations (%s) matching the reference criterion, ' - 'there should be only one reference' % len(refSimRowHash)) + message = ( + "Found multiple simulations (%s) matching the reference criterion, " + "there should be only one reference" % len(refSimRowHash) + ) log.error(message) raise NameError(message) refSimRowHash = refSimRowHash[0] - refSimName = inputsDF.loc[refSimRowHash, 'simName'] + refSimName = inputsDF.loc[refSimRowHash, "simName"] return refSimRowHash, refSimName -def checkAIMECinputs(cfgSetup, pathDict): - """ Check inputs before running AIMEC postprocessing +def resolveResTypeColumn(row, baseResType, layer=""): + """Resolve base resType to actual DataFrame column for a simulation row + + For single-layer sims: returns baseResType (e.g. 'ppr') + For multi-layer sims with layer set: returns layer-suffixed column (e.g. 'ppr_l2') + Returns None if no matching column has data. + + Parameters + ---------- + row : pandas Series + a row from inputsDF or resAnalysisDF (one row per simulation with info on + the simulation, available result files and optional configuration) + baseResType : str + base result type name (e.g. 'ppr', 'pfv') + layer : str + layer identifier (e.g. 'L2'). Empty string for single-layer. + + Returns + ------- + str or None + actual column name, or None if not resolvable + """ + if layer: + layeredCol = "%s_%s" % (baseResType, layer.lower()) + if layeredCol in row.index and pd.notna(row[layeredCol]): + return layeredCol + # single-layer sim or no layer configured + if baseResType in row.index and pd.notna(row[baseResType]): + return baseResType + return None + + +def computeBaseResTypeList(inputsDF, runoutLayer): + """Compute base resType names available for all simulations + + Discovers all base resType names from DataFrame columns (stripping _l\\d+ suffixes), + then checks that every row can resolve each base name via resolveResTypeColumn. + Returns only base names resolvable for ALL sims. + + Parameters + ---------- + inputsDF : pandas DataFrame + simulation dataframe with resType columns (base and/or layer-suffixed) + runoutLayer : str + layer identifier (e.g. 'L2'). Empty string for single-layer. + + Returns + ------- + list + sorted list of base resType names available for all simulations + """ + # metadata columns to exclude from resType discovery + metadataCols = { + "simName", + "releaseArea", + "simHash", + "simModified", + "simType", + "modelType", + "cellSize", + "layers", + } + + # discover base resType names by stripping _l\d+ suffixes + baseNames = set() + for col in inputsDF.columns: + if col in metadataCols: + continue + baseName = re.sub(r"_l\d+$", "", col) + baseNames.add(baseName) + + # check per-row resolvability for each base name + resolvableForAll = [] + for baseName in sorted(baseNames): + allResolve = True + for idx in inputsDF.index: + if resolveResTypeColumn(inputsDF.loc[idx], baseName, layer=runoutLayer) is None: + allResolve = False + break + if allResolve: + resolvableForAll.append(baseName) + + return resolvableForAll + + +def checkAIMECinputs(cfgSetup, pathDict, inputsDF): + """Check inputs before running AIMEC postprocessing - Make sure that the available data satisfies what is required in the ini file + Make sure that the available data satisfies what is required in the ini file. + Uses per-row resolution to validate that base resType names are available for + all simulations (handling both single-layer and multi-layer results). + Errors if multi-layer data is detected but no runoutLayer is specified. Parameters ---------- cfgSetup : configParser - aimec configuration + aimec configuration, reads runoutResType, resTypes and runoutLayer pathDict: dict aimec input dictionary (path to inputs and refSimName and hash) + inputsDF: pandas DataFrame + simulation dataframe with resType columns Returns ------- pathDict : dict - aimec input dictionary updated with the resTypeList (list of result file types to take into account) + aimec input dictionary updated with: + - resTypeList: set of base result type names to analyze + - runoutResType: base runout result type name + - runoutLayer: layer identifier for downstream resolution + - displayRunoutResType: display name for plot titles """ - # check that the resTypes and runoutResType asked for in the ini file are available for all simulations - # if not, raise an error - # what we have for all simulations - resTypeList = pathDict['resTypeList'] - # what we need for all simulations - if cfgSetup['resTypes'] != '': - # required res types - resTypesWanted = cfgSetup['resTypes'].split('|') + # runoutLayer is optional: not needed for single-layer sims, required for multi-layer sims + runoutLayer = cfgSetup.get("runoutLayer", "") + runoutResType = cfgSetup["runoutResType"] + + # detect multi-layer data from inputsDF columns + if not runoutLayer: + layerColumns = [col for col in inputsDF.columns if re.search(r"_l\d+$", col)] + if layerColumns: + message = ( + "Multi-layer result files detected (%s) but no runoutLayer is set. " + "Set runoutLayer (e.g. L1) in ana3AIMECCfg.ini to select a layer for analysis." + % ", ".join(sorted(layerColumns)) + ) + log.error(message) + raise FileNotFoundError(message) + + # compute base resType names available for all sims + baseResTypeList = computeBaseResTypeList(inputsDF, runoutLayer) + + # determine base resTypes used for analysis, if not provided in aimec cfg take base res types available for all sims + if cfgSetup["resTypes"] != "": + resTypesWanted = cfgSetup["resTypes"].split("|") else: - resTypesWanted = copy.deepcopy(pathDict['resTypeList']) + resTypesWanted = list(baseResTypeList) + # add the runout res type to what we need - resTypesWanted.append(cfgSetup['runoutResType']) - resTypesWanted = set(resTypesWanted) - # check for differences - match = resTypesWanted & set(resTypeList) # what is in both - diff = list(resTypesWanted.difference(match)) # what is in resTypesWanted but not in both - if diff != []: - message = '%s result type should be available for all simulations to analyze.' % diff + if runoutResType not in resTypesWanted: + resTypesWanted.append(runoutResType) + + # check that all wanted base types are available for all sims + diff = [rt for rt in resTypesWanted if rt not in baseResTypeList] + if diff: + message = "%s result type should be available for all simulations to analyze." % diff log.error(message) raise FileNotFoundError(message) + + pathDict["resTypeList"] = sorted(set(resTypesWanted)) + pathDict["runoutResType"] = runoutResType + pathDict["runoutLayer"] = runoutLayer + # display name for plot titles: "ppr (if multilayer: L2)" for multi-layer, "ppr" for single-layer + if runoutLayer: + pathDict["displayRunoutResType"] = "%s (if multilayer: %s)" % (runoutResType, runoutLayer) else: - pathDict['resTypeList'] = resTypesWanted - log.info('Analyzing %s results types. Runout based on %s result type' % - (pathDict['resTypeList'], cfgSetup['runoutResType'])) + pathDict["displayRunoutResType"] = runoutResType + log.info( + "Analyzing %s results types. Runout based on %s result type" + % (pathDict["resTypeList"], runoutResType) + ) return pathDict def computeCellSizeSL(cfgSetup, refCellSize): - """ Get the new (s, l) coordinate cell size - read by default the reference result file cell size. - If a 'cellSizeSL' is specified in cfgSetup then use this one - - Parameters - ----------- - refCellSize: float - cell size of reference simulation - cfgSetup: configParser object - configuration for aimec - with field cellSizeSL - - Returns - -------- - cellSizeSL: float - cell size to be used for the (s, l) coordinates + """Get the new (s, l) coordinate cell size + read by default the reference result file cell size. + If a 'cellSizeSL' is specified in cfgSetup then use this one + + Parameters + ----------- + refCellSize: float + cell size of reference simulation + cfgSetup: configParser object + configuration for aimec - with field cellSizeSL + + Returns + -------- + cellSizeSL: float + cell size to be used for the (s, l) coordinates """ - if cfgSetup['cellSizeSL'] == '': + if cfgSetup["cellSizeSL"] == "": cellSizeSL = float(refCellSize) - log.info('cellSizeSL is read from the reference header and is : %.2f m' % cellSizeSL) + log.info("cellSizeSL is read from the reference header and is : %.2f m" % cellSizeSL) else: try: - cellSizeSL = cfgSetup.getfloat('cellSizeSL') - log.info('cellSizeSL is read from the configuration file and is : %.2f m' % cellSizeSL) + cellSizeSL = cfgSetup.getfloat("cellSizeSL") + log.info("cellSizeSL is read from the configuration file and is : %.2f m" % cellSizeSL) except ValueError: - message = ('cellSizeSL is read from the configuration file but should be a number, you provided: %s' - % cfgSetup['cellSizeSL']) + message = ( + "cellSizeSL is read from the configuration file but should be a number, you provided: %s" + % cfgSetup["cellSizeSL"] + ) log.error(message) raise ValueError(message) @@ -449,7 +602,7 @@ def computeCellSizeSL(cfgSetup, refCellSize): def makeDomainTransfo(pathDict, dem, refCellSize, cfgSetup): - """ Make domain transformation + """Make domain transformation This function returns the information about the domain transformation Data given on a regular grid is projected on a nonuniform grid following @@ -491,17 +644,17 @@ def makeDomainTransfo(pathDict, dem, refCellSize, cfgSetup): index for start of the runout area (in s) if defineRunoutArea is False - indStartOfRunout=0 (start of thalweg) """ - w = cfgSetup.getfloat('domainWidth') + w = cfgSetup.getfloat("domainWidth") # get the cell size for the (s, l) raster cellSizeSL = computeCellSizeSL(cfgSetup, refCellSize) # Initialize transformation dictionary rasterTransfo = {} - rasterTransfo['domainWidth'] = w - rasterTransfo['cellSizeSL'] = cellSizeSL + rasterTransfo["domainWidth"] = w + rasterTransfo["cellSizeSL"] = cellSizeSL # read avaPath avaPath, splitPoint = setAvaPath(pathDict, dem) - rasterTransfo['avaPath'] = avaPath - rasterTransfo['splitPoint'] = splitPoint + rasterTransfo["avaPath"] = avaPath + rasterTransfo["splitPoint"] = splitPoint # Get new Domain Boundaries DB # input: ava path @@ -534,13 +687,13 @@ def makeDomainTransfo(pathDict, dem, refCellSize, cfgSetup): rasterTransfo = findStartOfRunoutArea(dem, rasterTransfo, cfgSetup, splitPoint) # add dem info to rasterTransfo - rasterTransfo['dem'] = dem + rasterTransfo["dem"] = dem return rasterTransfo def splitSection(DB, i): - """ Splits the ith segment of domain boundary DB in the s direction + """Splits the ith segment of domain boundary DB in the s direction (direction of the path) Parameters @@ -569,27 +722,27 @@ def splitSection(DB, i): number of elements on the new segments """ # left edge - xl0 = DB['DBXl'][i] - xl1 = DB['DBXl'][i+1] - yl0 = DB['DBYl'][i] - yl1 = DB['DBYl'][i+1] + xl0 = DB["DBXl"][i] + xl1 = DB["DBXl"][i + 1] + yl0 = DB["DBYl"][i] + yl1 = DB["DBYl"][i + 1] dxl = xl1 - xl0 dyl = yl1 - yl0 Vl = np.array((dxl, dyl)) zl = np.linalg.norm(Vl) # right edge - xr0 = DB['DBXr'][i] - xr1 = DB['DBXr'][i+1] - yr0 = DB['DBYr'][i] - yr1 = DB['DBYr'][i+1] + xr0 = DB["DBXr"][i] + xr1 = DB["DBXr"][i + 1] + yr0 = DB["DBYr"][i] + yr1 = DB["DBYr"][i + 1] dxr = xr1 - xr0 dyr = yr1 - yr0 Vr = np.array((dxr, dyr)) zr = np.linalg.norm(Vr) # number of segments - m = int(max(np.ceil(zl), np.ceil(zr))+1) + m = int(max(np.ceil(zl), np.ceil(zr)) + 1) # make left segment bxl = np.linspace(xl0, xl1, m) byl = np.linspace(yl0, yl1, m) @@ -601,7 +754,7 @@ def splitSection(DB, i): def makeTransfoMat(rasterTransfo): - """ Make transformation matrix. + """Make transformation matrix. Takes a Domain Boundary and finds the (x,y) coordinates of the new raster point (the one following the path) @@ -632,30 +785,30 @@ def makeTransfoMat(rasterTransfo): - l: 1D numpy array, new coord system in the cross direction """ - w = rasterTransfo['domainWidth'] - cellSize = rasterTransfo['cellSizeSL'] + w = rasterTransfo["domainWidth"] + cellSize = rasterTransfo["cellSizeSL"] # number of points describing the avaPath - n_pnt = np.shape(rasterTransfo['DBXr'])[0] + n_pnt = np.shape(rasterTransfo["DBXr"])[0] # Working with no dimentions # (the cellsize scaling will be readded at the end) # lcoord is the distance from the polyline (cross section) # the maximum step size should be smaller then the cellsize - nTot = np.ceil(w/cellSize) + nTot = np.ceil(w / cellSize) # take the next odd integer. This ensures that the lcoord = 0 exists - nTot = int(nTot+1) if ((nTot % 2) == 0) else int(nTot) - n2Tot = int(np.floor(nTot/2)) + nTot = int(nTot + 1) if ((nTot % 2) == 0) else int(nTot) + n2Tot = int(np.floor(nTot / 2)) lcoord = np.linspace(-n2Tot, n2Tot, nTot) # this way, 0 is in lcoord # initialize new_rasters newGridRasterX = np.array([]) # xcoord of the points of the new raster newGridRasterY = np.array([]) # ycoord of the points of the new raster # loop on each section of the path - for i in range(n_pnt-1): + for i in range(n_pnt - 1): # split edges in segments bxl, byl, bxr, byr, m = splitSection(rasterTransfo, i) # bxl, byl, bxr, byr reprensent the s direction (olong path) # loop on segments of section - for j in range(m-1): + for j in range(m - 1): # this is the cross section segment (l direction) x = np.linspace(bxl[j], bxr[j], nTot) # line coordinates x y = np.linspace(byl[j], byr[j], nTot) # line coordinates y @@ -664,26 +817,24 @@ def makeTransfoMat(rasterTransfo): newGridRasterX = x.reshape(1, nTot) newGridRasterY = y.reshape(1, nTot) else: - newGridRasterX = np.append(newGridRasterX, x.reshape(1, nTot), - axis=0) - newGridRasterY = np.append(newGridRasterY, y.reshape(1, nTot), - axis=0) + newGridRasterX = np.append(newGridRasterX, x.reshape(1, nTot), axis=0) + newGridRasterY = np.append(newGridRasterY, y.reshape(1, nTot), axis=0) # add last column - x = np.linspace(bxl[m-1], bxr[m-1], nTot) # line coordinates x - y = np.linspace(byl[m-1], byr[m-1], nTot) # line coordinates y + x = np.linspace(bxl[m - 1], bxr[m - 1], nTot) # line coordinates x + y = np.linspace(byl[m - 1], byr[m - 1], nTot) # line coordinates y newGridRasterX = np.append(newGridRasterX, x.reshape(1, nTot), axis=0) newGridRasterY = np.append(newGridRasterY, y.reshape(1, nTot), axis=0) - rasterTransfo['l'] = lcoord - rasterTransfo['gridx'] = newGridRasterX - rasterTransfo['gridy'] = newGridRasterY + rasterTransfo["l"] = lcoord + rasterTransfo["gridx"] = newGridRasterX + rasterTransfo["gridy"] = newGridRasterY return rasterTransfo def getSArea(rasterTransfo, dem): - """ Get the s curvilinear coordinate and area on the new raster + """Get the s curvilinear coordinate and area on the new raster Find the scoord corresponding to the transformation and the Area of the cells of the new raster @@ -710,8 +861,8 @@ def getSArea(rasterTransfo, dem): - rasterArea: 2D numpy array, real area of the cells of the new raster """ - xcoord = rasterTransfo['gridx'] - ycoord = rasterTransfo['gridy'] + xcoord = rasterTransfo["gridx"] + ycoord = rasterTransfo["gridy"] # add ghost lines and columns to the coord matrix # in order to perform dx and dy calculation n, m = np.shape(xcoord) @@ -722,45 +873,51 @@ def getSArea(rasterTransfo, dem): ycoord = np.append(ycoord, ycoord[-2, :].reshape(1, m), axis=0) n, m = np.shape(xcoord) # calculate dx and dy for each point in the s direction - dxs = xcoord[1:n, 0:m-1]-xcoord[0:n-1, 0:m-1] - dys = ycoord[1:n, 0:m-1]-ycoord[0:n-1, 0:m-1] + dxs = xcoord[1:n, 0 : m - 1] - xcoord[0 : n - 1, 0 : m - 1] + dys = ycoord[1:n, 0 : m - 1] - ycoord[0 : n - 1, 0 : m - 1] # deduce the distance in s direction - Vs2 = (dxs*dxs + dys*dys) + Vs2 = dxs * dxs + dys * dys Vs = np.sqrt(Vs2) # Method 1 # calculate area of each cell using Gauss's area formula # sum_i{x_(i)*y_(i+1) + y_(i)*x_(i+1)} # i = 1 : top left in the matrix : [0:n-1, 0:m-1] - x1 = xcoord[0:n-1, 0:m-1] - y1 = ycoord[0:n-1, 0:m-1] + x1 = xcoord[0 : n - 1, 0 : m - 1] + y1 = ycoord[0 : n - 1, 0 : m - 1] # i = 2 : top right in the matrix : [1:n, 0:m-1] - x2 = xcoord[1:n, 0:m-1] - y2 = ycoord[1:n, 0:m-1] + x2 = xcoord[1:n, 0 : m - 1] + y2 = ycoord[1:n, 0 : m - 1] # i = 3 : bottom right in the matrix : [1:n, 1:m] x3 = xcoord[1:n, 1:m] y3 = ycoord[1:n, 1:m] # i = 4 : bottom left in the matrix : [0:n-1, 1:m] - x4 = xcoord[0:n-1, 1:m] - y4 = ycoord[0:n-1, 1:m] + x4 = xcoord[0 : n - 1, 1:m] + y4 = ycoord[0 : n - 1, 1:m] - area = np.abs(x1*y2-y1*x2 + x2*y3-y2*x3 + x3*y4-y3*x4 + x4*y1-y4*x1)/2 + area = np.abs(x1 * y2 - y1 * x2 + x2 * y3 - y2 * x3 + x3 * y4 - y3 * x4 + x4 * y1 - y4 * x1) / 2 # save Area matrix - demCellSize = dem['header']['cellsize'] - xllcenter = dem['header']['xllcenter'] - yllcenter = dem['header']['yllcenter'] + demCellSize = dem["header"]["cellsize"] + xllcenter = dem["header"]["xllcenter"] + yllcenter = dem["header"]["yllcenter"] # area corection coef due to slope (1 if cell is horizontal, >1 if sloped) - demAreaCoef = dem['areaRaster']/(demCellSize*demCellSize) + demAreaCoef = dem["areaRaster"] / (demCellSize * demCellSize) # project on ney grid - areaCoef, _ = geoTrans.projectOnGrid(rasterTransfo['gridx']*demCellSize, rasterTransfo['gridy']*demCellSize, - demAreaCoef, csz=demCellSize, xllc=xllcenter, yllc=yllcenter) + areaCoef, _ = geoTrans.projectOnGrid( + rasterTransfo["gridx"] * demCellSize, + rasterTransfo["gridy"] * demCellSize, + demAreaCoef, + csz=demCellSize, + xllc=xllcenter, + yllc=yllcenter, + ) # real area - rasterTransfo['rasterArea'] = area * areaCoef + rasterTransfo["rasterArea"] = area * areaCoef # get scoord - ds = Vs[:, int(np.floor(m/2))-1] - scoord = np.cumsum(ds)-ds[0] - rasterTransfo['s'] = scoord + ds = Vs[:, int(np.floor(m / 2)) - 1] + scoord = np.cumsum(ds) - ds[0] + rasterTransfo["s"] = scoord return rasterTransfo @@ -790,20 +947,22 @@ def transform(data, name, rasterTransfo, interpMethod, dem=False): the new raster """ # read tranformation info - newGridRasterX = rasterTransfo['gridx'] - newGridRasterY = rasterTransfo['gridy'] + newGridRasterX = rasterTransfo["gridx"] + newGridRasterY = rasterTransfo["gridy"] n, m = np.shape(newGridRasterX) xx = newGridRasterX yy = newGridRasterY Points = {} - Points['x'] = xx.flatten() - Points['y'] = yy.flatten() - iib = len(Points['x']) + Points["x"] = xx.flatten() + Points["y"] = yy.flatten() + iib = len(Points["x"]) Points, ioob = geoTrans.projectOnRaster(data, Points, interp=interpMethod) - newData = Points['z'].reshape(n, m) - log.debug('Data-file: %s - %d raster values transferred - %d out of original raster' - 'bounds!' % (name, iib-ioob, ioob)) + newData = Points["z"].reshape(n, m) + log.debug( + "Data-file: %s - %d raster values transferred - %d out of original raster" + "bounds!" % (name, iib - ioob, ioob) + ) # # TODO: only required if check of folding on coordinate system where data is to be analysed # # check if data overlap in s, l coordinate system @@ -824,7 +983,7 @@ def transform(data, name, rasterTransfo, interpMethod, dem=False): def assignData(fnames, rasterTransfo, interpMethod): - """ Transfer data from old raster to new raster + """Transfer data from old raster to new raster Loops through paths in fnames and calls transfom Transform affects values to the points of the new raster @@ -849,7 +1008,7 @@ def assignData(fnames, rasterTransfo, interpMethod): maxtopo = len(fnames) avalData = np.array(([None] * maxtopo)) - log.debug('Transfer data of %d file(s) from old to new raster' % maxtopo) + log.debug("Transfer data of %d file(s) from old to new raster" % maxtopo) for i in range(maxtopo): fname = fnames[i] name = fname.name @@ -860,7 +1019,7 @@ def assignData(fnames, rasterTransfo, interpMethod): def analyzeMass(fnameMass, simRowHash, refSimRowHash, resAnalysisDF, time=None): - """ analyze Mass data + """analyze Mass data Parameters ---------- @@ -899,47 +1058,51 @@ def analyzeMass(fnameMass, simRowHash, refSimRowHash, resAnalysisDF, time=None): time: 1D numpy array time array for mass analysis """ - log.debug('Analyzing mass') - log.debug('{: <10} {: <10} {: <10}'.format('Sim number ', 'GI ', 'GR ')) + log.debug("Analyzing mass") + log.debug("{: <10} {: <10} {: <10}".format("Sim number ", "GI ", "GR ")) # analyze mass - releasedMass, entrainedMass, finalMass, grIndex, grGrad, entMassFlow, totalMass, time = readWrite(fnameMass, time) - resAnalysisDF.loc[simRowHash, 'relMass'] = releasedMass - resAnalysisDF.loc[simRowHash, 'finalMass'] = finalMass - resAnalysisDF.loc[simRowHash, 'entMass'] = entrainedMass - resAnalysisDF.loc[simRowHash, 'growthIndex'] = grIndex - resAnalysisDF.loc[simRowHash, 'growthGrad'] = grGrad - - releasedMassRef = resAnalysisDF.loc[refSimRowHash, 'relMass'] - finalMassRef = resAnalysisDF.loc[refSimRowHash, 'finalMass'] - relativMassDiff = (finalMass-finalMassRef)/finalMassRef*100 + releasedMass, entrainedMass, finalMass, grIndex, grGrad, entMassFlow, totalMass, time = readWrite( + fnameMass, time + ) + resAnalysisDF.loc[simRowHash, "relMass"] = releasedMass + resAnalysisDF.loc[simRowHash, "finalMass"] = finalMass + resAnalysisDF.loc[simRowHash, "entMass"] = entrainedMass + resAnalysisDF.loc[simRowHash, "growthIndex"] = grIndex + resAnalysisDF.loc[simRowHash, "growthGrad"] = grGrad + + releasedMassRef = resAnalysisDF.loc[refSimRowHash, "relMass"] + finalMassRef = resAnalysisDF.loc[refSimRowHash, "finalMass"] + relativMassDiff = (finalMass - finalMassRef) / finalMassRef * 100 if not (releasedMass == releasedMassRef): - log.warning('Release masses differ between simulations!') - log.info('{: <10} {:<10.4f} {:<10.4f}'.format(*[resAnalysisDF.loc[simRowHash, 'simName'], grIndex, grGrad])) - resAnalysisDF.loc[simRowHash, 'relativMassDiff'] = relativMassDiff + log.warning("Release masses differ between simulations!") + log.info( + "{: <10} {:<10.4f} {:<10.4f}".format(*[resAnalysisDF.loc[simRowHash, "simName"], grIndex, grGrad]) + ) + resAnalysisDF.loc[simRowHash, "relativMassDiff"] = relativMassDiff if simRowHash == refSimRowHash: - resAnalysisDF = pd.concat([resAnalysisDF, - pd.DataFrame({'entMassFlowArray': np.nan}, - index=resAnalysisDF.index)], - axis=1).copy() - resAnalysisDF['entMassFlowArray'] = resAnalysisDF['entMassFlowArray'].astype(object) - resAnalysisDF = pd.concat([resAnalysisDF, - pd.DataFrame({'totalMassArray': np.nan}, - index=resAnalysisDF.index)], - axis=1).copy() - resAnalysisDF['totalMassArray'] = resAnalysisDF['totalMassArray'].astype(object) - - resAnalysisDF.at[simRowHash, 'entMassFlowArray'] = entMassFlow - resAnalysisDF.at[simRowHash, 'totalMassArray'] = totalMass + resAnalysisDF = pd.concat( + [resAnalysisDF, pd.DataFrame({"entMassFlowArray": np.nan}, index=resAnalysisDF.index)], axis=1 + ).copy() + resAnalysisDF["entMassFlowArray"] = resAnalysisDF["entMassFlowArray"].astype(object) + resAnalysisDF = pd.concat( + [resAnalysisDF, pd.DataFrame({"totalMassArray": np.nan}, index=resAnalysisDF.index)], axis=1 + ).copy() + resAnalysisDF["totalMassArray"] = resAnalysisDF["totalMassArray"].astype(object) + + resAnalysisDF.at[simRowHash, "entMassFlowArray"] = entMassFlow + resAnalysisDF.at[simRowHash, "totalMassArray"] = totalMass return resAnalysisDF, time -def computeRunOut(cfgSetup, rasterTransfo, resAnalysisDF, transformedRasters, simRowHash): - """ Compute runout based on peak field results +def computeRunOut( + cfgSetup, rasterTransfo, resAnalysisDF, transformedRasters, simRowHash, runoutResType +): + """Compute runout based on peak field results Parameters ---------- - cfgSetup: confiParser + cfgSetup: configParser aimec analysis configuration rasterTransfo: dict transformation information @@ -953,6 +1116,8 @@ def computeRunOut(cfgSetup, rasterTransfo, resAnalysisDF, transformedRasters, si dict with transformed dem and peak results simRowHash: str simulation dataframe hash + runoutResType: str + resolved runout result type column name (e.g. "ppr" or "ppr_l1" for multi-layer). Returns ------- @@ -983,31 +1148,30 @@ def computeRunOut(cfgSetup, rasterTransfo, resAnalysisDF, transformedRasters, si """ # read inputs - scoord = rasterTransfo['s'] - lcoord = rasterTransfo['l'] - zThalweg = rasterTransfo['z'] + scoord = rasterTransfo["s"] + lcoord = rasterTransfo["l"] + zThalweg = rasterTransfo["z"] n = np.shape(lcoord)[0] - n = int(np.floor(n/2)+1) - x = rasterTransfo['x'] - y = rasterTransfo['y'] - gridx = rasterTransfo['gridx'] - gridy = rasterTransfo['gridy'] - - runoutResType = cfgSetup['runoutResType'] - thresholdValue = cfgSetup.getfloat('thresholdValue') - transformedDEMRasters = transformedRasters['newRasterDEM'] - PResRasters = transformedRasters['newRaster' + runoutResType.upper()] - PResCrossMax = resAnalysisDF.loc[simRowHash, runoutResType + 'CrossMax'] - PResCrossMean = resAnalysisDF.loc[simRowHash, runoutResType + 'CrossMean'] - - log.debug('Computing runout') + n = int(np.floor(n / 2) + 1) + x = rasterTransfo["x"] + y = rasterTransfo["y"] + gridx = rasterTransfo["gridx"] + gridy = rasterTransfo["gridy"] + + thresholdValue = cfgSetup.getfloat("thresholdValue") + transformedDEMRasters = transformedRasters["newRasterDEM"] + PResRasters = transformedRasters["newRaster" + runoutResType.upper()] + PResCrossMax = resAnalysisDF.loc[simRowHash, runoutResType + "CrossMax"] + PResCrossMean = resAnalysisDF.loc[simRowHash, runoutResType + "CrossMean"] + + log.debug("Computing runout") lindex = np.nonzero(PResCrossMax > thresholdValue)[0] if lindex.any(): cUpper = min(lindex) cLower = max(lindex) runoutFound = True else: - log.warning('No max values > threshold found. threshold = %4.2f, too high?' % thresholdValue) + log.warning("No max values > threshold found. threshold = %4.2f, too high?" % thresholdValue) cUpper = 0 cLower = 0 runoutFound = False @@ -1018,7 +1182,7 @@ def computeRunOut(cfgSetup, rasterTransfo, resAnalysisDF, transformedRasters, si cLowerm = max(lindex) flagMeanFound = True else: - log.warning('No average values > threshold found. threshold = %4.2f, too high?' % thresholdValue) + log.warning("No average values > threshold found. threshold = %4.2f, too high?" % thresholdValue) cUpperm = 0 cLowerm = 0 flagMeanFound = False @@ -1039,7 +1203,7 @@ def computeRunOut(cfgSetup, rasterTransfo, resAnalysisDF, transformedRasters, si resAnalysisDF.loc[simRowHash, "yRunout"] = np.nan resAnalysisDF.loc[simRowHash, "runoutAngle"] = np.nan - resAnalysisDF.loc[simRowHash, 'deltaSXY'] = scoord[cLower] - scoord[cUpper] + resAnalysisDF.loc[simRowHash, "deltaSXY"] = scoord[cLower] - scoord[cUpper] if flagMeanFound: resAnalysisDF.loc[simRowHash, "zRunout"] = zThalweg[cLower] @@ -1053,50 +1217,63 @@ def computeRunOut(cfgSetup, rasterTransfo, resAnalysisDF, transformedRasters, si resAnalysisDF.loc[simRowHash, "yMeanRunout"] = np.nan resAnalysisDF.loc[simRowHash, "zRelease"] = zThalweg[cUpper] - resAnalysisDF.loc[simRowHash, 'elevRel'] = zThalweg[cUpper] - resAnalysisDF.loc[simRowHash, 'deltaZ'] = zThalweg[cUpper] - zThalweg[cLower] + resAnalysisDF.loc[simRowHash, "elevRel"] = zThalweg[cUpper] + resAnalysisDF.loc[simRowHash, "deltaZ"] = zThalweg[cUpper] - zThalweg[cLower] resAnalysisDF.loc[simRowHash, "runoutFound"] = runoutFound return resAnalysisDF -def computeRunoutLine(cfgSetup, rasterTransfo, transformedRasters, simRowHash, actualType, name='', basedOnMax=False): - """ compute the runout line as for each l coordinate the furthest affected s coordinate using the desired - resType and thresholdvalue, or if basedOnMax searching for the max value (used for reference line, point) - also add runout point based on max s value of runout line - - Parameters - ----------- - cfgSetup: configparser object - configuration settings, here: runoutResType and tresholdValue - rasterTransfo: dict - information on raster transformation - transformedRasters: dict - transformed rasters from simulation - simRowHash: hash - index of current simulation - actualType: str - options: simulation, line, point, poly - name: str - optional name to find raster in transfromedRasters dict - basedOnMax: bool - if not threshold of runoutResType is used but max value along s for each l - - Returns - --------- - runoutLine : dict - coordinates of runout line in s, l, and x, y, index of rasterTransfo['s']mname and type - coordinates of runout point (max s extent) in s, l - """ +def computeRunoutLine( + cfgSetup, + rasterTransfo, + transformedRasters, + simRowHash, + actualType, + name="", + basedOnMax=False, + runoutResType=None, +): + """compute the runout line as for each l coordinate the furthest affected s coordinate using the desired + resType and thresholdvalue, or if basedOnMax searching for the max value (used for reference line, point) + also add runout point based on max s value of runout line + + Parameters + ----------- + cfgSetup: configparser object + configuration settings, here: runoutResType and tresholdValue + rasterTransfo: dict + information on raster transformation + transformedRasters: dict + transformed rasters from simulation + simRowHash: hash + index of current simulation + actualType: str + options: simulation, line, point, poly + name: str + optional name to find raster in transfromedRasters dict + basedOnMax: bool + if not threshold of runoutResType is used but max value along s for each l + runoutResType: str or None + resolved runout result type column name (e.g. "ppr" or "ppr_l1" for multi-layer). + Required when name is empty (simulation runout). Not used when name is set (reference data). + + Returns + --------- + runoutLine : dict + coordinates of runout line in s, l, and x, y, index of rasterTransfo['s']mname and type + coordinates of runout point (max s extent) in s, l + """ # get info on coordinate grids - gridx = rasterTransfo['gridx'] - gridy = rasterTransfo['gridy'] + gridx = rasterTransfo["gridx"] + gridy = rasterTransfo["gridy"] # fetch raster data - if name == '': - runoutResType = cfgSetup['runoutResType'] - anaRaster = transformedRasters['newRaster' + runoutResType.upper()] + if name == "": + if runoutResType is None: + raise ValueError("runoutResType is required when computing simulation runout (name is empty)") + anaRaster = transformedRasters["newRaster" + runoutResType.upper()] name = simRowHash else: anaRaster = transformedRasters[name] @@ -1106,59 +1283,63 @@ def computeRunoutLine(cfgSetup, rasterTransfo, transformedRasters, simRowHash, a anaRasterFlip = np.flip(anaRaster, axis=0) # setup runout line - thresholdValue = cfgSetup.getfloat('thresholdValue') - runoutLine = {'s': np.zeros(rasterTransfo['l'].shape) * np.nan, 'l': np.zeros(rasterTransfo['l'].shape) * np.nan, - 'x': np.zeros(rasterTransfo['l'].shape) * np.nan, 'y': np.zeros(rasterTransfo['l'].shape) * np.nan, - 'runoutLineFound': np.zeros(rasterTransfo['l'].shape) * np.nan} - runoutLineIndex = np.zeros(rasterTransfo['l'].shape) * np.nan + thresholdValue = cfgSetup.getfloat("thresholdValue") + runoutLine = { + "s": np.zeros(rasterTransfo["l"].shape) * np.nan, + "l": np.zeros(rasterTransfo["l"].shape) * np.nan, + "x": np.zeros(rasterTransfo["l"].shape) * np.nan, + "y": np.zeros(rasterTransfo["l"].shape) * np.nan, + "runoutLineFound": np.zeros(rasterTransfo["l"].shape) * np.nan, + } + runoutLineIndex = np.zeros(rasterTransfo["l"].shape) * np.nan # loop over each l coordinate to find max s coordinate sMaxInd = np.nan lMaxInd = np.nan sMax = 0 oneIndexFound = False - for ind, lCoor in enumerate(rasterTransfo['l']): + for ind, lCoor in enumerate(rasterTransfo["l"]): indexFound = False - if basedOnMax and (len(np.nonzero(anaRaster[:,ind])[0]) > 0): - index2 = np.argmax(anaRasterFlip[:,ind]) - index1 = anaRaster.shape[0]-index2 + if basedOnMax and (len(np.nonzero(anaRaster[:, ind])[0]) > 0): + index2 = np.argmax(anaRasterFlip[:, ind]) + index1 = anaRaster.shape[0] - index2 indexFound = True - oneIndexFound =True - elif (basedOnMax == False) and (len(np.nonzero(anaRaster[:,ind]> thresholdValue)[0]) > 0): - index1 = max(np.nonzero(anaRaster[:,ind] > thresholdValue)[0]) + oneIndexFound = True + elif (basedOnMax == False) and (len(np.nonzero(anaRaster[:, ind] > thresholdValue)[0]) > 0): + index1 = max(np.nonzero(anaRaster[:, ind] > thresholdValue)[0]) indexFound = True oneIndexFound = True if indexFound: runoutLineIndex[ind] = index1 - runoutLine['s'][ind] = (rasterTransfo['s'][index1]) - runoutLine['l'][ind] = lCoor - runoutLine['x'][ind] = gridx[index1, ind] - runoutLine['y'][ind] = gridy[index1, ind] - runoutLine['runoutLineFound'][ind] = True - if rasterTransfo['s'][index1] > sMax: + runoutLine["s"][ind] = rasterTransfo["s"][index1] + runoutLine["l"][ind] = lCoor + runoutLine["x"][ind] = gridx[index1, ind] + runoutLine["y"][ind] = gridy[index1, ind] + runoutLine["runoutLineFound"][ind] = True + if rasterTransfo["s"][index1] > sMax: sMaxInd = index1 lMaxInd = ind - sMax = rasterTransfo['s'][index1] + sMax = rasterTransfo["s"][index1] else: - runoutLine['runoutLineFound'][ind] = False + runoutLine["runoutLineFound"][ind] = False # add info on name, type - runoutLine['index'] = runoutLineIndex - runoutLine['name'] = name - runoutLine['type'] = actualType + runoutLine["index"] = runoutLineIndex + runoutLine["name"] = name + runoutLine["type"] = actualType if oneIndexFound: - runoutLine['sRunout'] = np.nanmax(runoutLine['s']) - runoutLine['lRunout'] = rasterTransfo['l'][lMaxInd] - runoutLine['xRunout'] = gridx[sMaxInd, lMaxInd] - runoutLine['yRunout'] = gridy[sMaxInd, lMaxInd] + runoutLine["sRunout"] = np.nanmax(runoutLine["s"]) + runoutLine["lRunout"] = rasterTransfo["l"][lMaxInd] + runoutLine["xRunout"] = gridx[sMaxInd, lMaxInd] + runoutLine["yRunout"] = gridy[sMaxInd, lMaxInd] else: - log.warning('For new Raster %s no runout line is found' % name) + log.warning("For new Raster %s no runout line is found" % name) return runoutLine def analyzeField(simRowHash, rasterTransfo, transformedRaster, dataType, resAnalysisDF): - """ analyze transformed field + """analyze transformed field analyze transformed raster: compute the Max and Mean values in each cross section as well as the overall maximum @@ -1187,24 +1368,33 @@ def analyzeField(simRowHash, rasterTransfo, transformedRaster, dataType, resAnal containing the mean of the field in each cross section """ # read inputs - rasterArea = rasterTransfo['rasterArea'] + rasterArea = rasterTransfo["rasterArea"] # initialize Arrays - log.debug('Analyzing %s' % (dataType)) - log.debug('{: <10} {: <10}'.format('Sim number ', 'maxCrossMax ')) + log.debug("Analyzing %s" % (dataType)) + log.debug("{: <10} {: <10}".format("Sim number ", "maxCrossMax ")) # Max Mean in each Cross-Section for each field maxaCrossMax, aCrossMax, aCrossMean = getMaxMeanValues(transformedRaster, rasterArea) - log.debug('{: <10} {:<10.4f}'.format(*[simRowHash, maxaCrossMax])) + log.debug("{: <10} {:<10.4f}".format(*[simRowHash, maxaCrossMax])) - resAnalysisDF.loc[simRowHash, 'max' + dataType + 'CrossMax'] = maxaCrossMax - resAnalysisDF.at[simRowHash, dataType + 'CrossMax'] = aCrossMax - resAnalysisDF.at[simRowHash, dataType + 'CrossMean'] = aCrossMean + resAnalysisDF.loc[simRowHash, "max" + dataType + "CrossMax"] = maxaCrossMax + resAnalysisDF.at[simRowHash, dataType + "CrossMax"] = aCrossMax + resAnalysisDF.at[simRowHash, dataType + "CrossMean"] = aCrossMean return resAnalysisDF -def analyzeArea(rasterTransfo, resAnalysisDF, simRowHash, newRasters, cfg, pathDict, contourDict, referenceType='newRefRaster'): +def analyzeArea( + rasterTransfo, + resAnalysisDF, + simRowHash, + newRasters, + cfg, + pathDict, + contourDict, + referenceType="newRefRaster", +): """Compare area results to reference. Compute True positive, False negative... areas. @@ -1247,39 +1437,42 @@ def analyzeArea(rasterTransfo, resAnalysisDF, simRowHash, newRasters, cfg, pathD for thresholdValue - updated """ - cfgSetup = cfg['AIMECSETUP'] - cfgPlots = cfg['PLOTS'] - runoutResType = cfgSetup['runoutResType'] - refSimRowHash = pathDict['refSimRowHash'] - cellarea = rasterTransfo['rasterArea'] - indStartOfRunout = rasterTransfo['indStartOfRunout'] - thresholdValue = cfgSetup.getfloat('thresholdValue') - contourLevels = fU.splitIniValueToArraySteps(cfgSetup['contourLevels']) - simName = resAnalysisDF.loc[simRowHash, 'simName'] + cfgSetup = cfg["AIMECSETUP"] + cfgPlots = cfg["PLOTS"] + runoutResType = pathDict["runoutResType"] + refSimRowHash = pathDict["refSimRowHash"] + cellarea = rasterTransfo["rasterArea"] + indStartOfRunout = rasterTransfo["indStartOfRunout"] + thresholdValue = cfgSetup.getfloat("thresholdValue") + contourLevels = fU.splitIniValueToArraySteps(cfgSetup["contourLevels"]) + simName = resAnalysisDF.loc[simRowHash, "simName"] # rasterinfo nStart = indStartOfRunout # inputs for plot inputs = {} - inputs['runoutLength'] = resAnalysisDF.loc[refSimRowHash, 'sRunout'] - inputs['refData'] = newRasters[referenceType + runoutResType.upper()] - inputs['nStart'] = nStart - inputs['runoutResType'] = runoutResType + inputs["runoutLength"] = resAnalysisDF.loc[refSimRowHash, "sRunout"] + inputs["refData"] = newRasters[referenceType + runoutResType.upper()] + inputs["nStart"] = nStart + inputs["runoutResType"] = runoutResType + # base name without layer suffix for cosmetic lookups (units, names, colorMaps) + inputs["baseRunoutResType"] = cfgSetup["runoutResType"] + inputs["displayRunoutResType"] = pathDict["displayRunoutResType"] thresholdArray = contourLevels thresholdArray = np.append(thresholdArray, thresholdValue) - inputs['thresholdArray'] = thresholdArray - inputs['diffLim'] = cfgSetup.getfloat('diffLim') + inputs["thresholdArray"] = thresholdArray + inputs["diffLim"] = cfgSetup.getfloat("diffLim") - rasterdata = newRasters['newRaster' + runoutResType.upper()] + rasterdata = newRasters["newRaster" + runoutResType.upper()] # take first simulation as reference - refMask = copy.deepcopy(inputs['refData']) + refMask = copy.deepcopy(inputs["refData"]) # prepare mask for area resAnalysis refMask = np.where(np.isnan(refMask), 0, refMask) refMask = np.where(refMask <= thresholdValue, 0, refMask) refMask = np.where(refMask > thresholdValue, 1, refMask) # comparison rasterdata with mask - log.debug('{: <15} {: <15} {: <15} {: <15} {: <15}'.format('Sim number ', 'TP ', 'FN ', 'FP ', 'TN')) + log.debug("{: <15} {: <15} {: <15} {: <15} {: <15}".format("Sim number ", "TP ", "FN ", "FP ", "TN")) compRasterMask = copy.deepcopy(rasterdata) # prepare mask for area resAnalysis compRasterMask = np.where(np.isnan(compRasterMask), 0, compRasterMask) @@ -1301,22 +1494,27 @@ def analyzeArea(rasterTransfo, resAnalysisDF, simRowHash, newRasters, cfg, pathD areaSum = tp + fn if areaSum > 0: - log.debug('{: <15} {:<15.4f} {:<15.4f} {:<15.4f} {:<15.4f}'.format( - *[simName, tp/areaSum, fn/areaSum, fp/areaSum, tn/areaSum])) + log.debug( + "{: <15} {:<15.4f} {:<15.4f} {:<15.4f} {:<15.4f}".format( + *[simName, tp / areaSum, fn / areaSum, fp / areaSum, tn / areaSum] + ) + ) if tp + fp == 0: - log.warning('Simulation %s did not reach the run-out area for threshold %.2f %s' % - (simName, thresholdValue, pU.cfgPlotUtils['unit' + cfgSetup['runoutResType']])) + log.warning( + "Simulation %s did not reach the run-out area for threshold %.2f %s" + % (simName, thresholdValue, pU.cfgPlotUtils["unit" + cfgSetup["runoutResType"]]) + ) - resAnalysisDF.loc[simRowHash, 'TP'] = tp - resAnalysisDF.loc[simRowHash, 'TN'] = tn - resAnalysisDF.loc[simRowHash, 'FP'] = fp - resAnalysisDF.loc[simRowHash, 'FN'] = fn + resAnalysisDF.loc[simRowHash, "TP"] = tp + resAnalysisDF.loc[simRowHash, "TN"] = tn + resAnalysisDF.loc[simRowHash, "FP"] = fp + resAnalysisDF.loc[simRowHash, "FN"] = fn # inputs for plot - inputs['compData'] = rasterdata + inputs["compData"] = rasterdata # masked data for the dataThreshold given in the ini file - inputs['refRasterMask'] = refMask - inputs['compRasterMask'] = compRasterMask - inputs['simName'] = simName + inputs["refRasterMask"] = refMask + inputs["compRasterMask"] = compRasterMask + inputs["simName"] = simName compPlotPath = None if ( @@ -1332,35 +1530,37 @@ def analyzeArea(rasterTransfo, resAnalysisDF, simRowHash, newRasters, cfg, pathD % resAnalysisDF.loc[simRowHash, "simName"] ) # add contourlines to contourDict - contourDict = outAimec.fetchContourLines(rasterTransfo, inputs, cfgSetup.getfloat('thresholdValue'), contourDict) + contourDict = outAimec.fetchContourLines( + rasterTransfo, inputs, cfgSetup.getfloat("thresholdValue"), contourDict + ) return resAnalysisDF, compPlotPath, contourDict def readWrite(fname_ent, time): - """ Get mass balance information - Read mass balance files to get mass properties of the simulation - (total mass, entrained mass...). Checks for mass conservation - - Parameters - ---------- - fname_ent: list of str - list of pah to mass balance files - Returns - ------- - relMass: float - release mass - entMassentMassFlow: float - entrained mass flow (kg/s) - finalMass: float - final mass - growthIndex: float - growthGrad: float + """Get mass balance information + Read mass balance files to get mass properties of the simulation + (total mass, entrained mass...). Checks for mass conservation + + Parameters + ---------- + fname_ent: list of str + list of pah to mass balance files + Returns + ------- + relMass: float + release mass + entMassentMassFlow: float + entrained mass flow (kg/s) + finalMass: float + final mass + growthIndex: float + growthGrad: float """ # load data # time, total mass, entrained mass - massTime = np.loadtxt(fname_ent, delimiter=',', skiprows=1) + massTime = np.loadtxt(fname_ent, delimiter=",", skiprows=1) timeSimulation = massTime[:, 0] dt = timeSimulation[1:] - timeSimulation[:-1] timeResults = [massTime[0, 0], massTime[-1, 0]] @@ -1368,33 +1568,34 @@ def readWrite(fname_ent, time): relMass = totMassResults[0] if time is None: time = np.arange(0, int(timeResults[1]), 0.1) - entMassFlow = np.interp(time, massTime[1:, 0], massTime[1:, 2]/dt) + entMassFlow = np.interp(time, massTime[1:, 0], massTime[1:, 2] / dt) totalMass = np.interp(time, massTime[:, 0], massTime[:, 1]) entrainedMass = np.sum(massTime[:, 2]) finalMass = totMassResults[1] # check mass balance - log.info('Total mass change between first and last time step in sim %s is: %.1f kg' % - (fname_ent.stem, totMassResults[1] - relMass)) - log.info('Total entrained mass in sim %s is: %.1f kg' % - (fname_ent.stem, entrainedMass)) + log.info( + "Total mass change between first and last time step in sim %s is: %.1f kg" + % (fname_ent.stem, totMassResults[1] - relMass) + ) + log.info("Total entrained mass in sim %s is: %.1f kg" % (fname_ent.stem, entrainedMass)) if (totMassResults[1] - relMass) == 0: diff = np.abs((totMassResults[1] - relMass) - entrainedMass) if diff > 0: - log.warning('Conservation of mass is not satisfied') - log.warning('Total mass change and total entrained mass differ from %.4f kg' % (diff)) + log.warning("Conservation of mass is not satisfied") + log.warning("Total mass change and total entrained mass differ from %.4f kg" % (diff)) else: - log.info('Total mass change and total entrained mass differ from %.4f kg' % (diff)) + log.info("Total mass change and total entrained mass differ from %.4f kg" % (diff)) else: - diff = np.abs((totMassResults[1] - relMass) - entrainedMass)/(totMassResults[1] - relMass) - if diff*100 > 0.05: - log.warning('Conservation of mass is not satisfied') - log.warning('Total mass change and total entrained mass differ from %.4f %%' % (diff*100)) + diff = np.abs((totMassResults[1] - relMass) - entrainedMass) / (totMassResults[1] - relMass) + if diff * 100 > 0.05: + log.warning("Conservation of mass is not satisfied") + log.warning("Total mass change and total entrained mass differ from %.4f %%" % (diff * 100)) else: - log.info('Total mass change and total entrained mass differ from %.4f %%' % (diff*100)) + log.info("Total mass change and total entrained mass differ from %.4f %%" % (diff * 100)) # growth results - growthIndex = totMassResults[1]/totMassResults[0] + growthIndex = totMassResults[1] / totMassResults[0] growthGrad = (totMassResults[1] - totMassResults[0]) / (timeResults[1] - timeResults[0]) return relMass, entrainedMass, finalMass, growthIndex, growthGrad, entMassFlow, totalMass, time @@ -1426,7 +1627,7 @@ def getMaxMeanValues(rasterdataA, rasterArea): rasterArea = np.where(np.where(np.isnan(rasterdataA), 0, rasterdataA) > 0, rasterArea, 0) areaSum = np.nansum(rasterArea, axis=1) areaSum = np.where(areaSum > 0, areaSum, 1) - aCrossMean = np.nansum(rasterdataA*rasterArea, axis=1)/areaSum + aCrossMean = np.nansum(rasterdataA * rasterArea, axis=1) / areaSum # aCrossMean = np.nanmean(rasterdataA, axis=1) aCrossMax = np.nanmax(rasterdataA, 1) # maximum of field a @@ -1436,28 +1637,28 @@ def getMaxMeanValues(rasterdataA, rasterArea): def setAvaPath(pathDict, dem): - """ fetch path shapefile, prepare for AIMEC, set z-coordinate - - Parameters - ----------- - pathDict: dict - info on path to aimec path, split Point, - dem: dict - dictionary with DEM header and data - - Returns - -------- - avaPath: dict - info on aimec path coordinates, ... - splitPoint: dict - info on split point coordinates + """fetch path shapefile, prepare for AIMEC, set z-coordinate + + Parameters + ----------- + pathDict: dict + info on path to aimec path, split Point, + dem: dict + dictionary with DEM header and data + + Returns + -------- + avaPath: dict + info on aimec path coordinates, ... + splitPoint: dict + info on split point coordinates """ # fetch input parameters - profileLayer = pathDict['profileLayer'] - splitPointSource = pathDict['splitPointSource'] - defaultName = pathDict['projectName'] + profileLayer = pathDict["profileLayer"] + splitPointSource = pathDict["splitPointSource"] + defaultName = pathDict["projectName"] # read avaPath avaPath = shpConv.readLine(profileLayer, defaultName, dem) @@ -1471,45 +1672,45 @@ def setAvaPath(pathDict, dem): # reverse avaPath if necessary _, avaPath = geoTrans.checkProfile(avaPath, projSplitPoint=None) - log.debug('Creating new raster along polyline: %s' % profileLayer) + log.debug("Creating new raster along polyline: %s" % profileLayer) return avaPath, splitPoint def scalePathWithCellSize(rasterTransfo, cellSizeSL): - """ use desired cellsize to scale ava path and create s, l coordinates - and coordinates of the resampled polyline in the old coord systems x, y - - Parameters - ---------- - rasterTransfo: dict - domain transformation info - cellSizeSL: float - desired cell size of sl raster - - Returns - -------- - rasterTransfo: dict - updated dictionary with new coordinate system + """use desired cellsize to scale ava path and create s, l coordinates + and coordinates of the resampled polyline in the old coord systems x, y + + Parameters + ---------- + rasterTransfo: dict + domain transformation info + cellSizeSL: float + desired cell size of sl raster + + Returns + -------- + rasterTransfo: dict + updated dictionary with new coordinate system """ # put back the scale due to the desired cellsize - rasterTransfo['s'] = rasterTransfo['s']*cellSizeSL - rasterTransfo['l'] = rasterTransfo['l']*cellSizeSL - rasterTransfo['gridx'] = rasterTransfo['gridx']*cellSizeSL - rasterTransfo['gridy'] = rasterTransfo['gridy']*cellSizeSL - rasterTransfo['rasterArea'] = rasterTransfo['rasterArea']*cellSizeSL*cellSizeSL + rasterTransfo["s"] = rasterTransfo["s"] * cellSizeSL + rasterTransfo["l"] = rasterTransfo["l"] * cellSizeSL + rasterTransfo["gridx"] = rasterTransfo["gridx"] * cellSizeSL + rasterTransfo["gridy"] = rasterTransfo["gridy"] * cellSizeSL + rasterTransfo["rasterArea"] = rasterTransfo["rasterArea"] * cellSizeSL * cellSizeSL # (x,y) coordinates of the resampled avapth (centerline where l = 0) - n = np.shape(rasterTransfo['l'])[0] - indCenter = int(np.floor(n/2)) - rasterTransfo['x'] = rasterTransfo['gridx'][:, indCenter] - rasterTransfo['y'] = rasterTransfo['gridy'][:, indCenter] + n = np.shape(rasterTransfo["l"])[0] + indCenter = int(np.floor(n / 2)) + rasterTransfo["x"] = rasterTransfo["gridx"][:, indCenter] + rasterTransfo["y"] = rasterTransfo["gridy"][:, indCenter] return rasterTransfo def addSurfaceParalleCoord(rasterTransfo): - """ Add the surface parallel coordinate (along flow path taking elevation into account) + """Add the surface parallel coordinate (along flow path taking elevation into account) Parameters ---------- @@ -1521,97 +1722,117 @@ def addSurfaceParalleCoord(rasterTransfo): rasterTransfo: dict updated dictionary with the surface parallel coordinate 'sParallel' """ - dz = np.diff(rasterTransfo['z'], prepend=rasterTransfo['z'][0]) - ds = np.diff(rasterTransfo['s'], prepend=rasterTransfo['s'][0]) - dsParallel2 = (ds*ds + dz*dz) + dz = np.diff(rasterTransfo["z"], prepend=rasterTransfo["z"][0]) + ds = np.diff(rasterTransfo["s"], prepend=rasterTransfo["s"][0]) + dsParallel2 = ds * ds + dz * dz dsParallel = np.sqrt(dsParallel2) - sParallel = np.cumsum(dsParallel)-dsParallel[0] - rasterTransfo['sParallel'] = sParallel + sParallel = np.cumsum(dsParallel) - dsParallel[0] + rasterTransfo["sParallel"] = sParallel return rasterTransfo def findStartOfRunoutArea(dem, rasterTransfo, cfgSetup, splitPoint): - """ find start of runout area point using splitPoint - add info on x, y coordinates of point, angle, index - - if defineRunoutArea=False - start of runout area is equal to the start of the thalweg - -> in this case the entire SL domain represents the runout area - - Parameters - ----------- - dem: dict - dictionary with DEM header and data - rasterTransfo: dict - domain transformation info - cfgSetup: configparser object - configuration for aimec - splitPoint: dict - dictionary with split Point coordinates - - Returns - -------- - rasterTransfo: dict - updated dictionary with info on start of runout location,... + """find start of runout area point using splitPoint + add info on x, y coordinates of point, angle, index + - if defineRunoutArea=False - start of runout area is equal to the start of the thalweg + -> in this case the entire SL domain represents the runout area + + Parameters + ----------- + dem: dict + dictionary with DEM header and data + rasterTransfo: dict + domain transformation info + cfgSetup: configparser object + configuration for aimec + splitPoint: dict + dictionary with split Point coordinates + + Returns + -------- + rasterTransfo: dict + updated dictionary with info on start of runout location,... """ if splitPoint != None: # fetch input parameters from config - startOfRunoutAreaAngle = cfgSetup.getfloat('startOfRunoutAreaAngle') + startOfRunoutAreaAngle = cfgSetup.getfloat("startOfRunoutAreaAngle") # add 'z' coordinate to the centerline rasterTransfo, _ = geoTrans.projectOnRaster(dem, rasterTransfo) # find projection of split point on the centerline projPoint = geoTrans.findSplitPoint(rasterTransfo, splitPoint) - rasterTransfo['indSplit'] = projPoint['indSplit'] - rasterTransfo['projSplitPoint'] = projPoint + rasterTransfo["indSplit"] = projPoint["indSplit"] + rasterTransfo["projSplitPoint"] = projPoint # prepare find start of runout area points angle, tmp, ds = geoTrans.prepareAngleProfile(startOfRunoutAreaAngle, rasterTransfo) # find the runout point: first point under startOfRunoutAreaAngle - indStartOfRunout = geoTrans.findAngleProfile(tmp, ds, cfgSetup.getfloat('dsMin')) - rasterTransfo['startOfRunoutAreaAngle'] = angle[indStartOfRunout] - rasterTransfo['labelRunout'] = ('start of runout area: ' + (r'$\beta_{%.1f °}$' % - rasterTransfo['startOfRunoutAreaAngle'])) + indStartOfRunout = geoTrans.findAngleProfile(tmp, ds, cfgSetup.getfloat("dsMin")) + rasterTransfo["startOfRunoutAreaAngle"] = angle[indStartOfRunout] + rasterTransfo["labelRunout"] = "start of runout area: " + ( + r"$\beta_{%.1f °}$" % rasterTransfo["startOfRunoutAreaAngle"] + ) else: - log.info('DefineRunoutArea is set to False - start of runout area set to start of thalweg') + log.info("DefineRunoutArea is set to False - start of runout area set to start of thalweg") indStartOfRunout = 0 - rasterTransfo['startOfRunoutAreaAngle'] = np.nan - rasterTransfo['labelRunout'] = 'start of runout area' + rasterTransfo["startOfRunoutAreaAngle"] = np.nan + rasterTransfo["labelRunout"] = "start of runout area" # add info to rasterTransfo dict - rasterTransfo['indStartOfRunout'] = indStartOfRunout - rasterTransfo['xBetaPoint'] = rasterTransfo['x'][indStartOfRunout] - rasterTransfo['yBetaPoint'] = rasterTransfo['y'][indStartOfRunout] + rasterTransfo["indStartOfRunout"] = indStartOfRunout + rasterTransfo["xBetaPoint"] = rasterTransfo["x"][indStartOfRunout] + rasterTransfo["yBetaPoint"] = rasterTransfo["y"][indStartOfRunout] - log.info('Start of run-out area at the %.2f ° point of coordinates (%.2f, %.2f)' % - (rasterTransfo['startOfRunoutAreaAngle'], rasterTransfo['xBetaPoint'], rasterTransfo['yBetaPoint'])) + log.info( + "Start of run-out area at the %.2f ° point of coordinates (%.2f, %.2f)" + % (rasterTransfo["startOfRunoutAreaAngle"], rasterTransfo["xBetaPoint"], rasterTransfo["yBetaPoint"]) + ) return rasterTransfo def addFieldsToDF(inputsDF): - """ add fields that will be added in aimec analysis to dataframe - - Parameters - ----------- - inputsDF: pandas DataFrame - DataFrame where fields should be added as empty columns - fields are: - cfgSetup: configparser object - configuration settings, here used: includeReference if inputs from Inputs/REFDATA shall be included and - used to compare sim results to - - Returns - --------- - inputsDF: pandas DataFrame - updated DataFrame + """add fields that will be added in aimec analysis to dataframe + + Parameters + ----------- + inputsDF: pandas DataFrame + DataFrame where fields should be added as empty columns + fields are: + cfgSetup: configparser object + configuration settings, here used: includeReference if inputs from Inputs/REFDATA shall be included and + used to compare sim results to + + Returns + --------- + inputsDF: pandas DataFrame + updated DataFrame """ - nanFields = ['sRunout', 'lRunout', 'xRunout', 'yRunout', 'deltaSXY', 'runoutAngle', 'zRelease', 'zRunout', - 'sMeanRunout', 'xMeanRunout', 'yMeanRunout', 'elevRel', 'deltaZ', 'refSim_Diff_sRunout', - 'refSim_Diff_lRunout', 'dataType', 'runoutLineDiff_line_RMSE', 'runoutLineDiff_poly_RMSE'] + nanFields = [ + "sRunout", + "lRunout", + "xRunout", + "yRunout", + "deltaSXY", + "runoutAngle", + "zRelease", + "zRunout", + "sMeanRunout", + "xMeanRunout", + "yMeanRunout", + "elevRel", + "deltaZ", + "refSim_Diff_sRunout", + "refSim_Diff_lRunout", + "dataType", + "runoutLineDiff_line_RMSE", + "runoutLineDiff_poly_RMSE", + ] emptyStrFields = [ "runoutLineDiff_line_pointsNotFoundInSim", @@ -1625,47 +1846,59 @@ def addFieldsToDF(inputsDF): inputsDF = pd.concat([inputsDF, pd.DataFrame({item: np.nan}, index=inputsDF.index)], axis=1).copy() for item in emptyStrFields: - inputsDF = pd.concat([inputsDF, pd.DataFrame({item: ''}, index=inputsDF.index)], axis=1).copy() + inputsDF = pd.concat([inputsDF, pd.DataFrame({item: ""}, index=inputsDF.index)], axis=1).copy() # add that datatype is simulation - inputsDF['dataType'] = ['simulation']*len(inputsDF) + inputsDF["dataType"] = ["simulation"] * len(inputsDF) return inputsDF def createReferenceDF(pathDict): - """ create data frame with one row per reference data set found in Inputs/REFDATA - Parameters - ----------- - pathDict: dict - dictionary with info paths to reference datasets + """create data frame with one row per reference data set found in Inputs/REFDATA + Parameters + ----------- + pathDict: dict + dictionary with info paths to reference datasets - Returns - --------- - referenceDF: pandas DataFrame - DataFrame with one row per reference dataset + Returns + --------- + referenceDF: pandas DataFrame + DataFrame with one row per reference dataset """ # add parameters that provide info on reference data and analysis that will be performed on this data - nanFields = ['reference_Type', 'reference_resType', 'reference_sRunout', 'reference_lRunout', - 'reference_xRunout', 'reference_yRunout', 'reference_name', 'reference_filePath', - 'dataType'] + nanFields = [ + "reference_Type", + "reference_resType", + "reference_sRunout", + "reference_lRunout", + "reference_xRunout", + "reference_yRunout", + "reference_name", + "reference_filePath", + "dataType", + ] # create empty dataframe referenceDF = pd.DataFrame(columns=nanFields) # add one row for each reference dataset - refData = [pathDict['referenceLine'], pathDict['referencePolygon'], pathDict['referencePoint']] + refData = [pathDict["referenceLine"], pathDict["referencePolygon"], pathDict["referencePoint"]] for ref in refData: if ref != []: for refFile in ref: referenceName = refFile.stem - newLine = pd.DataFrame([[referenceName, refFile]], columns=['reference_name', 'reference_filePath'], index=[referenceName]) + newLine = pd.DataFrame( + [[referenceName, refFile]], + columns=["reference_name", "reference_filePath"], + index=[referenceName], + ) hashRef = pd.util.hash_pandas_object(newLine) newLine = newLine.set_index(hashRef) referenceDF = pd.concat([referenceDF, newLine], ignore_index=False) - referenceDF.loc[hashRef, 'dataType'] = 'reference' + referenceDF.loc[hashRef, "dataType"] = "reference" referenceDF.loc[hashRef, "reference_Type"] = refFile.stem.split("_")[1] # TODO add here if additional info read from shp or a textfile? @@ -1674,164 +1907,176 @@ def createReferenceDF(pathDict): def computeRunoutPointDiff(resAnalysisDF, refData, simRowHash): - """ compute difference between runout point of simulation and runout point of refData - - Parameters - ------------- - resAnalysisDF: data frame - one row per simulation (config info, sim results and analysis) - refData: dict - dictionary of reference data set - simRowHash: - index of row in data frame for current simulation - - Returns - -------- - resAnalysisDF: data frame - updated DF with runout point difference + """compute difference between runout point of simulation and runout point of refData + + Parameters + ------------- + resAnalysisDF: data frame + one row per simulation (config info, sim results and analysis) + refData: dict + dictionary of reference data set + simRowHash: + index of row in data frame for current simulation + + Returns + -------- + resAnalysisDF: data frame + updated DF with runout point difference """ - for item in ['sRunout', 'lRunout']: + for item in ["sRunout", "lRunout"]: simPoint = resAnalysisDF.loc[simRowHash, item] refPoint = refData[item] diffRunout = refPoint - simPoint - resAnalysisDF.loc[simRowHash, 'refSim_Diff_%s' % item] = diffRunout - log.info('%s Runout diff for %s is %.2f' % (item, resAnalysisDF.loc[simRowHash, 'simName'], diffRunout)) + resAnalysisDF.loc[simRowHash, "refSim_Diff_%s" % item] = diffRunout + log.info( + "%s Runout diff for %s is %.2f" % (item, resAnalysisDF.loc[simRowHash, "simName"], diffRunout) + ) return resAnalysisDF def findSLCoors(rasterTransfo, refData, referenceType): - """ find closest s,l coordinates to given x,y coordinates using transformation matrix + """find closest s,l coordinates to given x,y coordinates using transformation matrix - Parameters - ----------- - rasterTransfo - refData + Parameters + ----------- + rasterTransfo + refData """ - xx = rasterTransfo['gridx'] - yy = rasterTransfo['gridy'] + xx = rasterTransfo["gridx"] + yy = rasterTransfo["gridy"] nrows = xx.shape[0] ncols = xx.shape[1] sizeXX = xx.size indexArray = np.arange(sizeXX).reshape(nrows, ncols) - coors = np.column_stack((np.reshape(xx, sizeXX), np.reshape(yy, sizeXX))) + coors = np.column_stack((np.reshape(xx, sizeXX), np.reshape(yy, sizeXX))) # initialize index1Array index1Array = None - if referenceType.lower() == 'point': + if referenceType.lower() == "point": # arrange x,y coordinates of ref Data into a 2D array for finding closest coordinates of coors test = np.zeros((len(coors), 2)) - test[:, 0] = refData['x'] - test[:, 1] = refData['y'] + test[:, 0] = refData["x"] + test[:, 1] = refData["y"] indexBool = np.isclose(coors, test) - index1Array = [np.where(np.all(indexBool==True, axis=1) == True)[0][0]] + index1Array = [np.where(np.all(indexBool == True, axis=1) == True)[0][0]] - elif referenceType.lower() == 'line': + elif referenceType.lower() == "line": indRowArray = [] indColArray = [] index1Array = [] - for x1, y1 in zip(refData['x'], refData['y']): + for x1, y1 in zip(refData["x"], refData["y"]): # arrange x,y coordinates of ref Data into a 2D array for finding closest coordinates of coors test = np.zeros((len(coors), 2)) test[:, 0] = x1 test[:, 1] = y1 indexBool = np.isclose(coors, test) - if len(np.where(np.all(indexBool==True, axis=1))[0]) > 1: - index1 = np.where(np.all(indexBool==True, axis=1) == True)[0][0] - indRowArray.append(np.where(indexArray==index1)[0][0]) - indColArray.append(np.where(indexArray==index1)[1][0]) + if len(np.where(np.all(indexBool == True, axis=1))[0]) > 1: + index1 = np.where(np.all(indexBool == True, axis=1) == True)[0][0] + indRowArray.append(np.where(indexArray == index1)[0][0]) + indColArray.append(np.where(indexArray == index1)[1][0]) index1Array.append(index1) if index1Array is None: Lcoors = None Scoors = None else: - gridL, gridS = np.meshgrid(rasterTransfo['l'], rasterTransfo['s']) - slCoors = np.column_stack((np.reshape(gridL, gridL.size), np.reshape(gridS, gridS.size))) - Lcoors = slCoors[index1Array,0] - Scoors = slCoors[index1Array,1] + gridL, gridS = np.meshgrid(rasterTransfo["l"], rasterTransfo["s"]) + slCoors = np.column_stack((np.reshape(gridL, gridL.size), np.reshape(gridS, gridS.size))) + Lcoors = slCoors[index1Array, 0] + Scoors = slCoors[index1Array, 1] return Lcoors, Scoors def addReferenceAnalysisTODF(referenceDF, refFile, refDataDict): - """ add reference analysis info to referenceDF - - Parameters - ------------ - referenceDF: data frame - one row per reference dataset - refFile: pathlib Path - file path of reference dataset file - refDataDict: dict - dict with info on reference data set and computed runout point - - Returns - ---------- - referenceDF: data frame - updated data DF with info on runout point coordinate + """add reference analysis info to referenceDF + + Parameters + ------------ + referenceDF: data frame + one row per reference dataset + refFile: pathlib Path + file path of reference dataset file + refDataDict: dict + dict with info on reference data set and computed runout point + + Returns + ---------- + referenceDF: data frame + updated data DF with info on runout point coordinate """ - for item in ['sRunout', 'lRunout', 'xRunout', 'yRunout']: - referenceDF.loc[referenceDF['reference_name'] == refFile.stem, 'reference_%s' % item] = refDataDict[item] + for item in ["sRunout", "lRunout", "xRunout", "yRunout"]: + referenceDF.loc[referenceDF["reference_name"] == refFile.stem, "reference_%s" % item] = refDataDict[ + item + ] return referenceDF def analyzeDiffsRunoutLines(cfgSetup, runoutLine, refDataTransformed, resAnalysisDF, simRowHash, pathDict): - """ analyze difference between runout line derived from simRaster and reference datasets (lines, polygons) - - Parameters - ------------- - runoutLine: dict - info on runout line derived from simRaster - refDataTransformed: dict - info on runout line derived from reference datasets - resAnalysisDF: data frame - one row per simulation - simRowHash: ID - index of row for current sim in resAnalysisDF - pathDict: dict - dict with info on where to save plots - - Returns - ---------- - resAnalysisDF: data frame - updated DF + """analyze difference between runout line derived from simRaster and reference datasets (lines, polygons) + + Parameters + ------------- + runoutLine: dict + info on runout line derived from simRaster + refDataTransformed: dict + info on runout line derived from reference datasets + resAnalysisDF: data frame + one row per simulation + simRowHash: ID + index of row for current sim in resAnalysisDF + pathDict: dict + dict with info on where to save plots + + Returns + ---------- + resAnalysisDF: data frame + updated DF """ for item in refDataTransformed: refLine = refDataTransformed[item] - if refLine['type'] in ['line', 'poly']: + if refLine["type"] in ["line", "poly"]: # simName is - simName = resAnalysisDF.loc[simRowHash, 'simName'] + simName = resAnalysisDF.loc[simRowHash, "simName"] # add runoutLine differences as vector - resAnalysisDF.at[simRowHash, 'runoutLineDiff_%s' % refLine['type']] = np.asarray(runoutLine['s'] - refLine['s']) + resAnalysisDF.at[simRowHash, "runoutLineDiff_%s" % refLine["type"]] = np.asarray( + runoutLine["s"] - refLine["s"] + ) # compute difference between runoutLine from simulation and reference data - diffAll = runoutLine['s'] - refLine['s'] + diffAll = runoutLine["s"] - refLine["s"] # check where both lines have points to compute RMSE - diffNoNans = diffAll[np.where((np.isnan(refLine['s']) == False) & (np.isnan(runoutLine['s']) == False))] + diffNoNans = diffAll[ + np.where((np.isnan(refLine["s"]) == False) & (np.isnan(runoutLine["s"]) == False)) + ] # check number of points where runoutLine sim has no points but reference Line has points - runoutLineNoPoints = len(np.where(np.isnan(runoutLine['s'][np.where(np.isnan(refLine['s']) == False)]))[0]) - runoutLineAllPoints = len(np.where(np.isnan(runoutLine['s']) == False)[0]) + runoutLineNoPoints = len( + np.where(np.isnan(runoutLine["s"][np.where(np.isnan(refLine["s"]) == False)]))[0] + ) + runoutLineAllPoints = len(np.where(np.isnan(runoutLine["s"]) == False)[0]) # check number of points where refLine has no points but runout line sim has - refLineNoPoints = len(np.where(np.isnan(refLine['s'][np.where(np.isnan(runoutLine['s']) == False)]))[0]) - refLineAllPoints = len(np.where(np.isnan(refLine['s']) == False)[0]) - runoutStr = '%d/%d' % (runoutLineNoPoints, refLineAllPoints) - refLineStr = '%d/%d' % (refLineNoPoints, runoutLineAllPoints) + refLineNoPoints = len( + np.where(np.isnan(refLine["s"][np.where(np.isnan(runoutLine["s"]) == False)]))[0] + ) + refLineAllPoints = len(np.where(np.isnan(refLine["s"]) == False)[0]) + runoutStr = "%d/%d" % (runoutLineNoPoints, refLineAllPoints) + refLineStr = "%d/%d" % (refLineNoPoints, runoutLineAllPoints) # compute RMSE between runout line and refLine # set RMSE to nan if no runout line found for sim to compare to reference @@ -1851,12 +2096,16 @@ def analyzeDiffsRunoutLines(cfgSetup, runoutLine, refDataTransformed, resAnalysi % simName ) - resAnalysisDF.at[simRowHash, 'runoutLineDiff_%s_pointsNotFoundInSim' % refLine['type']] = runoutStr - resAnalysisDF.at[simRowHash, 'runoutLineDiff_%s_pointsNotFoundInRef' % refLine['type']] = refLineStr - resAnalysisDF.at[simRowHash, 'runoutLineDiff_%s_RMSE' % refLine['type']] = RMSE + resAnalysisDF.at[simRowHash, "runoutLineDiff_%s_pointsNotFoundInSim" % refLine["type"]] = ( + runoutStr + ) + resAnalysisDF.at[simRowHash, "runoutLineDiff_%s_pointsNotFoundInRef" % refLine["type"]] = ( + refLineStr + ) + resAnalysisDF.at[simRowHash, "runoutLineDiff_%s_RMSE" % refLine["type"]] = RMSE else: - log.info('For reference data type %s, runout line comparison is not available' % refLine['type']) + log.info("For reference data type %s, runout line comparison is not available" % refLine["type"]) return resAnalysisDF diff --git a/avaframe/ana3AIMEC/ana3AIMEC.py b/avaframe/ana3AIMEC/ana3AIMEC.py index 8af04b01c..e3954de4c 100644 --- a/avaframe/ana3AIMEC/ana3AIMEC.py +++ b/avaframe/ana3AIMEC/ana3AIMEC.py @@ -66,7 +66,7 @@ def fullAimecAnalysis(avalancheDir, cfg, inputDir='', demFileName=''): 'demFileName': demFileName} pathDict = aimecTools.readAIMECinputs(avalancheDir, pathDict, cfgSetup.getboolean('defineRunoutArea'), dirName=anaMod) - pathDict = aimecTools.checkAIMECinputs(cfgSetup, pathDict) + pathDict = aimecTools.checkAIMECinputs(cfgSetup, pathDict, inputsDF) log.info("Running ana3AIMEC model on test case DEM \n %s \n with profile \n %s ", pathDict['demSource'], pathDict['profileLayer']) # Run AIMEC postprocessing @@ -117,7 +117,10 @@ def mainAIMEC(pathDict, inputsDF, cfg): # get hash of the reference refSimRowHash = pathDict['refSimRowHash'] # read reference file and raster and config - refResultSource = inputsDF.loc[refSimRowHash, cfgSetup['runoutResType']] + runoutResType = pathDict['runoutResType'] + runoutLayer = pathDict.get('runoutLayer', '') + refResTypeCol = aimecTools.resolveResTypeColumn(inputsDF.loc[refSimRowHash], runoutResType, runoutLayer) + refResultSource = inputsDF.loc[refSimRowHash, refResTypeCol] refRaster = IOf.readRaster(refResultSource) refRasterData = refRaster['rasterData'] refHeader = refRaster['header'] @@ -209,7 +212,14 @@ def mainAIMEC(pathDict, inputsDF, cfg): compResTypes1 = cfgPlots['compResType1'].split('|') compResTypes2 = cfgPlots['compResType2'].split('|') for indRes, compResType in enumerate(compResTypes1): - outAimec.plotMaxValuesComp(pathDict, resAnalysisDF, compResType, compResTypes2[indRes], + compResType2 = compResTypes2[indRes] + missingCols = [c for c in [compResType, compResType2] if c not in resAnalysisDF.columns] + if missingCols: + log.warning("Skipping comparison plot %s vs %s: columns %s not available in results " + "(result type not in resTypeList for all simulations)", + compResType, compResType2, missingCols) + continue + outAimec.plotMaxValuesComp(pathDict, resAnalysisDF, compResType, compResType2, hue=cfgPlots['scenarioName']) return rasterTransfo, resAnalysisDF, plotDict, newRasters @@ -331,13 +341,15 @@ def postProcessAIMEC(cfg, rasterTransfo, pathDict, resAnalysisDF, newRasters, ti flagMass = cfgFlags.getboolean('flagMass') refSimRowHash = pathDict['refSimRowHash'] resTypeList = pathDict['resTypeList'] + runoutLayer = pathDict.get('runoutLayer', '') # apply domain transformation log.info('Analyzing data in path coordinate system') for resType in resTypeList: log.debug("Assigning %s data to deskewed raster" % resType) - inputFiles = resAnalysisDF.loc[simRowHash, resType] + resTypeCol = aimecTools.resolveResTypeColumn(resAnalysisDF.loc[simRowHash], resType, runoutLayer) + inputFiles = resAnalysisDF.loc[simRowHash, resTypeCol] if isinstance(inputFiles, pathlib.PurePath): rasterData = IOf.readRaster(inputFiles) newRaster = aimecTools.transform(rasterData, inputFiles, rasterTransfo, interpMethod) @@ -380,16 +392,19 @@ def postProcessAIMEC(cfg, rasterTransfo, pathDict, resAnalysisDF, newRasters, ti simRowHash, rasterTransfo, newRaster, resType, resAnalysisDF ) - # compute runout based on runoutResType - resAnalysisDF = aimecTools.computeRunOut(cfgSetup, rasterTransfo, resAnalysisDF, newRasters, simRowHash) - runoutLine = aimecTools.computeRunoutLine(cfgSetup, rasterTransfo, newRasters, simRowHash, 'simulation', name='') + # compute runout based on base runoutResType (resolver handles per-row column access) + runoutResType = pathDict['runoutResType'] + resAnalysisDF = aimecTools.computeRunOut(cfgSetup, rasterTransfo, resAnalysisDF, newRasters, simRowHash, + runoutResType=runoutResType) + runoutLine = aimecTools.computeRunoutLine(cfgSetup, rasterTransfo, newRasters, simRowHash, 'simulation', name='', + runoutResType=runoutResType) if 'refPoint' in refDataTransformed: # compute differences between runout points resAnalysisDF = aimecTools.computeRunoutPointDiff(resAnalysisDF, refDataTransformed['refPoint'], simRowHash) # plot comparison between runout lines - outAimec.compareRunoutLines(cfgSetup, refDataTransformed, newRasters['newRaster'+cfgSetup['runoutResType'].upper()], + outAimec.compareRunoutLines(cfgSetup, refDataTransformed, newRasters['newRaster'+runoutResType.upper()], runoutLine, rasterTransfo, resAnalysisDF.loc[simRowHash], pathDict) # analyze distribution of diffs between runout lines diff --git a/avaframe/ana3AIMEC/ana3AIMECCfg.ini b/avaframe/ana3AIMEC/ana3AIMECCfg.ini index 597a41bcc..e8e08b678 100644 --- a/avaframe/ana3AIMEC/ana3AIMECCfg.ini +++ b/avaframe/ana3AIMEC/ana3AIMECCfg.ini @@ -24,6 +24,10 @@ resTypes = ppr|pft|pfv # data result type for runout analysis (ppr, pft, pfv) runoutResType = ppr +# layer to use for analysis (e.g. L1, L2). Leave empty for single-layer modules. +# Required when analyzing multi-layer results (e.g. com8MoTPSA). +runoutLayer = + # show cross profile max or mean values on slComparisonstat plot (options: Max, Mean) runoutCrossType = Max diff --git a/avaframe/ana3AIMEC/dfa2Aimec.py b/avaframe/ana3AIMEC/dfa2Aimec.py index 115ac43f3..eaebf63d1 100644 --- a/avaframe/ana3AIMEC/dfa2Aimec.py +++ b/avaframe/ana3AIMEC/dfa2Aimec.py @@ -10,6 +10,7 @@ # local modules from avaframe.in3Utils import fileHandlerUtils as fU from avaframe.in3Utils import cfgUtils +from avaframe.ana3AIMEC import aimecTools # create local logger # change log level in calling module to DEBUG to see log messages @@ -147,8 +148,6 @@ def dfaBench2Aimec(avaDir, cfg, simNameRef='', simNameComp=''): # Load all infos on comparison module simulations compData, resTypeCompList = fU.makeSimFromResDF(avaDir, None, inputDir=inputDirComp, simName=simNameComp) - resTypeList = list(set(resTypeRefList).intersection(resTypeCompList)) - pathDict['resTypeList'] = resTypeList pathDict['colorParameter'] = colorParameter pathDict['refSimRowHash'] = refSimRowHash pathDict['refSimName'] = refSimName @@ -182,6 +181,12 @@ def dfaBench2Aimec(avaDir, cfg, simNameRef='', simNameComp=''): message = ('Multiple rows of the dataFrame have the same simulation name.') log.warning(message) + # compute base-name resTypeList and validate via checkAIMECinputs + # use raw column-name lists as initial resTypeList (will be overwritten by checkAIMECinputs) + resTypeList = list(set(resTypeRefList).intersection(resTypeCompList)) + pathDict['resTypeList'] = resTypeList + pathDict = aimecTools.checkAIMECinputs(cfgSetup, pathDict, inputsDF) + return inputsDF, pathDict diff --git a/avaframe/com8MoTPSA/com8MoTPSA.py b/avaframe/com8MoTPSA/com8MoTPSA.py index e7e49b7c4..5cdad5418 100644 --- a/avaframe/com8MoTPSA/com8MoTPSA.py +++ b/avaframe/com8MoTPSA/com8MoTPSA.py @@ -24,6 +24,15 @@ def com8MoTPSAMain(cfgMain, cfgInfo=None): + """Run the full MoT-PSA workflow: generate configs, run simulations in parallel, postprocess. + + Parameters + ---------- + cfgMain : configparser.ConfigParser + main AvaFrame configuration (avalancheDir, nCPU, plot flags) + cfgInfo : dict or None, optional + override configuration info passed to MoTGenerateConfigs + """ # Get all necessary information from the configuration files currentModule = sys.modules[__name__] simDict, inputSimFiles = mT.MoTGenerateConfigs(cfgMain, cfgInfo, currentModule) @@ -62,7 +71,56 @@ def com8MoTPSAMain(cfgMain, cfgInfo=None): com8MoTPSAPostprocess(simDict, cfgMain, inputSimFiles) +def copyRawToLayerPeakFiles(workDir, outputDirPeakFile): + """Rename and copy raw MoT-PSA output files to peakFiles with L1/L2 layer naming. + + MoT-PSA produces raw suffixes (p1/p2_max, h1/h2_max, s1/s2_max) that map + to AvaFrame result types (ppr, pfd, pfv). The digit in the raw suffix encodes + the layer: 1 -> L1 (dense flow), 2 -> L2 (powder snow). + + Example: simKey_null_psa_p1_max.asc -> simKey_null_psa_L1_ppr.asc + + Parameters + ---------- + workDir : pathlib.Path + simulation work directory containing raw MoT-PSA output files + outputDirPeakFile : pathlib.Path + target directory for renamed peak files + """ + # Each entry: (glob pattern, raw L1 suffix, raw L2 suffix, AvaFrame resType) + layerRenameMap = [ + ("*p?_max*", "p1_max", "p2_max", "ppr"), + ("*h?_max*", "h1_max", "h2_max", "pfd"), + ("*s?_max*", "s1_max", "s2_max", "pfv"), + ] + for globPattern, rawL1, rawL2, resType in layerRenameMap: + rawFiles = list(workDir.glob(globPattern)) + # Replace raw suffixes with layer naming (e.g. p1_max -> L1_ppr, p2_max -> L2_ppr) + targetFiles = [ + pathlib.Path(str(f.name).replace(rawL1, "L1_%s" % resType).replace(rawL2, "L2_%s" % resType)) + for f in rawFiles + ] + # Prepend output directory and copy + targetFiles = [outputDirPeakFile / f for f in targetFiles] + for source, target in zip(rawFiles, targetFiles): + shutil.copy2(source, target) + + def com8MoTPSAPostprocess(simDict, cfgMain, inputSimFiles): + """Postprocess MoT-PSA results: rename outputs to L1/L2 peak files, generate plots and reports. + + For each simulation, copies DataTime.txt and renames raw MoT-PSA output files + (p1/p2_max, h1/h2_max, s1/s2_max) to AvaFrame layer naming (L1/L2 + ppr/pfd/pfv). + + Parameters + ---------- + simDict : dict + simulation dictionary + cfgMain : configparser.ConfigParser + main AvaFrame configuration (avalancheDir, plot flags) + inputSimFiles : dict + input file paths, must contain "demFile" + """ avalancheDir = cfgMain["MAIN"]["avalancheDir"] # Copy max files to output directory @@ -73,55 +131,11 @@ def com8MoTPSAPostprocess(simDict, cfgMain, inputSimFiles): for key in simDict: workDir = pathlib.Path(avalancheDir) / "Work" / "com8MoTPSA" / str(key) - # identify simType - simType = simDict[key]["simType"] - # Copy DataTime.txt dataTimeFile = workDir / "DataTime.txt" shutil.copy2(dataTimeFile, outputDir / (str(key) + "_DataTime.txt")) - # TODO: functionize it - # Copy ppr files - pprFiles = list(workDir.glob("*p?_max*")) - targetFiles = [ - pathlib.Path(str(f.name).replace("%s_psa_p1_max" % simType, "%s_dfa_ppr" % simType)) - for f in pprFiles - ] - targetFiles = [ - pathlib.Path(str(f).replace("%s_psa_p2_max" % simType, "%s_psa_ppr" % simType)) - for f in targetFiles - ] - targetFiles = [outputDirPeakFile / f for f in targetFiles] - for source, target in zip(pprFiles, targetFiles): - shutil.copy2(source, target) - - # Copy pfd files - pfdFiles = list(workDir.glob("*h?_max*")) - targetFiles = [ - pathlib.Path(str(f.name).replace("%s_psa_h1_max" % simType, "%s_dfa_pfd" % simType)) - for f in pfdFiles - ] - targetFiles = [ - pathlib.Path(str(f).replace("%s_psa_h2_max" % simType, "%s_psa_pfd" % simType)) - for f in targetFiles - ] - targetFiles = [outputDirPeakFile / f for f in targetFiles] - for source, target in zip(pfdFiles, targetFiles): - shutil.copy2(source, target) - - # Copy pfv files - pfvFiles = list(workDir.glob("*s?_max*")) - targetFiles = [ - pathlib.Path(str(f.name).replace("%s_psa_s1_max" % simType, "%s_dfa_pfv" % simType)) - for f in pfvFiles - ] - targetFiles = [ - pathlib.Path(str(f).replace("%s_psa_s2_max" % simType, "%s_psa_pfv" % simType)) - for f in targetFiles - ] - targetFiles = [outputDirPeakFile / f for f in targetFiles] - for source, target in zip(pfvFiles, targetFiles): - shutil.copy2(source, target) + copyRawToLayerPeakFiles(workDir, outputDirPeakFile) # create plots and report modName = __name__.split(".")[-1] @@ -134,6 +148,18 @@ def com8MoTPSAPostprocess(simDict, cfgMain, inputSimFiles): def com8MoTPSATask(rcfFile): + """Run a single MoT-PSA simulation by invoking the MoT-PSA executable with an rcf file. + + Parameters + ---------- + rcfFile : str or pathlib.Path + path to the .rcf configuration file for this simulation + + Returns + ------- + list + the command that was executed (["./MoT-PSA", rcfFile]) + """ # TODO: Obvious... os.chdir(os.path.dirname(os.path.abspath(__file__))) command = ["./MoT-PSA", rcfFile] @@ -144,6 +170,26 @@ def com8MoTPSATask(rcfFile): def com8MoTPSAPreprocess(simDict, inputSimFiles, cfgMain): + """Prepare all MoT-PSA simulations: derive input rasters, set config paths, write rcf files. + + For each simulation in simDict, processes release/entrainment/bed shear/deposition + input data, configures MoT-PSA file paths and parameters, and writes the .rcf + configuration file needed by the MoT-PSA executable. + + Parameters + ---------- + simDict : dict + simulation dictionary keyed by simKey, each value contains "cfgSim" and "simType" + inputSimFiles : dict + input file paths (DEM, release scenarios, etc.) + cfgMain : configparser.ConfigParser + main AvaFrame configuration (avalancheDir) + + Returns + ------- + list of pathlib.Path + paths to the generated .rcf files, one per simulation + """ # Load avalanche directory from general configuration file avalancheDir = cfgMain["MAIN"]["avalancheDir"] # set inputsDir where original input data and remeshed rasters are stored diff --git a/avaframe/com8MoTPSA/com8MoTPSACfg.ini b/avaframe/com8MoTPSA/com8MoTPSACfg.ini index 4fb42939a..6c5a76d08 100644 --- a/avaframe/com8MoTPSA/com8MoTPSACfg.ini +++ b/avaframe/com8MoTPSA/com8MoTPSACfg.ini @@ -5,6 +5,10 @@ # model type - only for file naming (psa - powder snow avalanche) modelType = psa +# available layers: L1 (dense flow), L2 (powder snow) +# if empty, all available layers are output +layers = L1|L2 + # list of simulations that shall be performed (null, ent, res, entres, available (use all available input data)) simTypeList = null @@ -73,6 +77,11 @@ secondaryRelThDistVariation = # secondary area release thickness (only considered if secondaryRelThFromFile=False) [m] secondaryRelTh = +#+++++++++++++general start conditions: time dependent release +# if timeDependentRelease is True, provide the the timesteps, thickness and velocity +# for a releases in a csv-file in the REL folder +timeDependentRelease = False + #+++++++++++++Volume classes [m³] volClassSmall = 25000. volClassMedium = 60000. diff --git a/avaframe/in3Utils/cfgUtils.py b/avaframe/in3Utils/cfgUtils.py index 5bba71ad9..dafd133ba 100644 --- a/avaframe/in3Utils/cfgUtils.py +++ b/avaframe/in3Utils/cfgUtils.py @@ -530,9 +530,10 @@ def parseSimName(name): """Parse simulation name handling both old and new formats. Auto-detects: - - Old format: relName_simHash_defID_[frictIndi]_simType_modelType[_resType][_timeStep] - - New format: relName_simHash_modName_defID_[frictIndi]_simType_modelType[_resType][_timeStep] + - Old format: relName_simHash_defID_[frictIndi]_simType_modelType[_layer][_resType][_timeStep] + - New format: relName_simHash_modName_defID_[frictIndi]_simType_modelType[_layer][_resType][_timeStep] [ ] denotes optional items + Layer component matches pattern L followed by digits (e.g., L1, L2, L12) Parameters ---------- @@ -550,6 +551,7 @@ def parseSimName(name): - frictIndi: str | None (optional, values: "S", "M", "L") - simType: str (required) - modelType: str (required) + - layer: str | None (optional, e.g., "L1", "L2") - resType: str | None (optional, only in filenames) - timeStep: str | None (optional, only in filenames) @@ -619,15 +621,25 @@ def parseSimName(name): simType = remainingParts[offset] modelType = remainingParts[offset + 1] - # Step 6: Extract optional file components (resType, timeStep) + # Step 6: Extract optional file components (layer, resType, timeStep) + layer = None resType = None timeStep = None if len(remainingParts) > offset + 2: - resType = remainingParts[offset + 2] - - if len(remainingParts) > offset + 3: - timeStep = remainingParts[offset + 3] + candidate = remainingParts[offset + 2] + if re.match(r"^L\d+$", candidate): + # Layer component detected (e.g., L1, L2, L12) + layer = candidate + if len(remainingParts) > offset + 3: + resType = remainingParts[offset + 3] + if len(remainingParts) > offset + 4: + timeStep = remainingParts[offset + 4] + else: + # No layer + resType = candidate + if len(remainingParts) > offset + 3: + timeStep = remainingParts[offset + 3] # Step 7: Return structured dictionary return { @@ -638,6 +650,7 @@ def parseSimName(name): "frictIndi": frictIndi, "simType": simType, "modelType": modelType, + "layer": layer, "resType": resType, "timeStep": timeStep, } diff --git a/avaframe/in3Utils/fileHandlerUtils.py b/avaframe/in3Utils/fileHandlerUtils.py index 91f42ea75..f2ea2fcb7 100644 --- a/avaframe/in3Utils/fileHandlerUtils.py +++ b/avaframe/in3Utils/fileHandlerUtils.py @@ -669,7 +669,7 @@ def makeSimFromResDF(avaDir, comModule, inputDir="", simName=""): dataframe with for each simulation, the full file path, file name, release area scenario, simulation type (null, entres, etc.), model type (dfa, ref, etc.), simID, path to result files (ppr, pft, etc.), simulation name, - cell size and optional name of avalanche, optional time step + cell size and optional name of avalanche, optional time step, layer name resTypeListAll: list list of res types available for all simulations """ @@ -695,8 +695,21 @@ def makeSimFromResDF(avaDir, comModule, inputDir="", simName=""): datafiles = list(inputDir.glob(ascName)) + list(inputDir.glob(tifName)) - # build the result data frame - resTypeListFromFiles = list(set([file.stem.split("_")[-1] for file in datafiles])) + # Parse all filenames once to determine column names and cache results + parsedFiles = [] + colNamesFromFiles = set() + for file in datafiles: + parts = cfgUtils.parseSimName(file.stem) + resType = parts["resType"] if parts["resType"] else file.stem.split("_")[-1] + layer = parts["layer"] + if layer: + colName = "%s_%s" % (resType, layer.lower()) + else: + colName = resType + colNamesFromFiles.add(colName) + parsedFiles.append((file, parts, resType, layer, colName)) + + # Build the result data frame columnsList = [ "simName", "releaseArea", @@ -705,17 +718,13 @@ def makeSimFromResDF(avaDir, comModule, inputDir="", simName=""): "simType", "modelType", "cellSize", - ] + resTypeListFromFiles + "layers", + ] + sorted(colNamesFromFiles) dataDF = pd.DataFrame(columns=columnsList) resTypeListOne = [] - for file in datafiles: + for file, simNameParts, resType, layer, colName in parsedFiles: name = file.stem - # Parse the filename to extract components - simNameParts = cfgUtils.parseSimName(name) - - # Extract simName (without resType/timeStep) and resType - resType = simNameParts["resType"] if simNameParts["resType"] else name.split("_")[-1] # Reconstruct simName without resType and timeStep # Preserve _AF_ separator if present in original name @@ -753,11 +762,18 @@ def makeSimFromResDF(avaDir, comModule, inputDir="", simName=""): # add info about the cell size header = IOf.readRasterHeader(file) dataDF.loc[simName, "cellSize"] = header["cellsize"] - # add full path to resType - dataDF.loc[simName, resType] = pathlib.Path(file) - # list all res types found - if resType not in resTypeListOne: - resTypeListOne.append(resType) + # add full path to resType column (layer-aware) + dataDF.loc[simName, colName] = pathlib.Path(file) + # update layers metadata for multi-layer files + if layer: + existingLayers = dataDF.loc[simName, "layers"] + if pd.isna(existingLayers): + dataDF.loc[simName, "layers"] = layer + elif layer not in existingLayers: + dataDF.loc[simName, "layers"] = existingLayers + "|" + layer + # list all column names found + if colName not in resTypeListOne: + resTypeListOne.append(colName) # add a hash for each line of the DF and use as index - required for identifcation hash = pd.util.hash_pandas_object(dataDF) diff --git a/avaframe/out3Plot/outAIMEC.py b/avaframe/out3Plot/outAIMEC.py index d26b33f47..8c39f5570 100644 --- a/avaframe/out3Plot/outAIMEC.py +++ b/avaframe/out3Plot/outAIMEC.py @@ -57,6 +57,7 @@ def visuTransfo(rasterTransfo, inputData, cfgSetup, pathDict): #################################### # Get input data runoutResType = cfgSetup['runoutResType'] + displayRunoutResType = pathDict.get('displayRunoutResType', runoutResType) unit = pU.cfgPlotUtils['unit' + runoutResType] # read paths projectName = pathDict['projectName'] @@ -137,13 +138,13 @@ def visuTransfo(rasterTransfo, inputData, cfgSetup, pathDict): ax2.text( 0.5, 0.5, - "reference only 0 values for %s" % runoutResType, + "reference only 0 values for %s" % displayRunoutResType, horizontalalignment="center", verticalalignment="center", transform=ax2.transAxes, ) else: - pU.addColorBar(im, ax2, ticks, unit, title=runoutResType) + pU.addColorBar(im, ax2, ticks, unit, title=displayRunoutResType) outFileName = '_'.join([projectName, 'DomainTransformation']) pU.saveAndOrPlot(pathDict, outFileName, fig) @@ -249,6 +250,9 @@ def visuRunoutStat(rasterTransfo, inputsDF, resAnalysisDF, newRasters, cfgSetup, percentile = cfgSetup.getfloat('percentile') runoutResType = cfgSetup['runoutResType'] + # layer-resolved name for data access (e.g. ppr_l2); falls back to base name for single-layer + resolvedRunoutResType = pathDict.get('runoutResType', runoutResType) + displayRunoutResType = pathDict.get('displayRunoutResType', runoutResType) thresholdValue = cfgSetup['thresholdValue'] unit = pU.cfgPlotUtils['unit' + runoutResType] name = pU.cfgPlotUtils['name' + runoutResType] @@ -258,10 +262,10 @@ def visuRunoutStat(rasterTransfo, inputsDF, resAnalysisDF, newRasters, cfgSetup, s = rasterTransfo['s'] l = rasterTransfo['l'] indStartOfRunout = rasterTransfo['indStartOfRunout'] - rasterdataPres = newRasters['newRefRaster' + runoutResType.upper()] + rasterdataPres = newRasters['newRefRaster' + resolvedRunoutResType.upper()] runout = resAnalysisDF['sRunout'].to_numpy() crossValue = 'Cross' + cfgSetup['runoutCrossType'].upper()[0] + cfgSetup['runoutCrossType'].lower()[1:] - pprCrossMax = np.stack(resAnalysisDF[runoutResType.lower() + crossValue].to_numpy()) + pprCrossMax = np.stack(resAnalysisDF[resolvedRunoutResType.lower() + crossValue].to_numpy()) ############################################ # compute mean, median and percenti. of peak field cross max values and mask array with threshold pMean = np.mean(pprCrossMax, axis=0) @@ -275,7 +279,7 @@ def visuRunoutStat(rasterTransfo, inputsDF, resAnalysisDF, newRasters, cfgSetup, if np.all(np.isnan(runout)): log.warning( "No %s values exceeding threshold of %.2f found to analyze" - % (runoutResType, float(thresholdValue)) + % (displayRunoutResType, float(thresholdValue)) ) outFilePath = "" else: @@ -403,7 +407,7 @@ def visuRunoutStat(rasterTransfo, inputsDF, resAnalysisDF, newRasters, cfgSetup, norm=normSC, marker=pU.markers[0], label=( - "runout points (%s<%.1f%s)" % (runoutResType, cfgSetup.getfloat("thresholdValue"), unit) + "runout points (%s<%.1f%s)" % (displayRunoutResType, cfgSetup.getfloat("thresholdValue"), unit) ), ) if displayColorBar: @@ -415,7 +419,7 @@ def visuRunoutStat(rasterTransfo, inputsDF, resAnalysisDF, newRasters, cfgSetup, ax1.set_xlim([xMin, xMax]) ax1.set_title( "%s field (reference) for (%s>%.1f%s)" - % (name, runoutResType, cfgSetup.getfloat("thresholdValue"), unit) + % (name, displayRunoutResType, cfgSetup.getfloat("thresholdValue"), unit) ) ax1.set_aspect("equal") pU.putAvaNameOnPlot(ax1, projectName) @@ -425,7 +429,7 @@ def visuRunoutStat(rasterTransfo, inputsDF, resAnalysisDF, newRasters, cfgSetup, ax1.text( 0.5, 0.5, - "reference only 0 values for %s" % runoutResType, + "reference only 0 values for %s" % displayRunoutResType, horizontalalignment="center", verticalalignment="center", transform=ax1.transAxes, @@ -451,7 +455,7 @@ def visuRunoutStat(rasterTransfo, inputsDF, resAnalysisDF, newRasters, cfgSetup, ax2.set_xlabel("$S_{xy}$ (thalweg) [m]") ax2.set_xlim([s.min(), s.max()]) ax2.set_ylim(auto=True) - ax2.set_ylabel("$%s_{%s}$ [%s]" % (runoutResType, crossValue, unit)) + ax2.set_ylabel("$%s_{%s}$ [%s]" % (displayRunoutResType, crossValue, unit)) # add middle panel with cross max values along s # loop over all sims and compute colorbar value and add line plot @@ -466,7 +470,7 @@ def visuRunoutStat(rasterTransfo, inputsDF, resAnalysisDF, newRasters, cfgSetup, if resAnalysisRow["simName"] == pathDict["refSimName"]: ax3.plot( s, - resAnalysisRow[runoutResType.lower() + crossValue], + resAnalysisRow[resolvedRunoutResType.lower() + crossValue], c="k", label="reference", zorder=nSamples + 1, @@ -476,12 +480,12 @@ def visuRunoutStat(rasterTransfo, inputsDF, resAnalysisDF, newRasters, cfgSetup, cmapVal = cmapValues3[varParStrValues.index(resAnalysisRow[varParList[0]])] ax3.plot( s, - resAnalysisRow[runoutResType.lower() + crossValue], + resAnalysisRow[resolvedRunoutResType.lower() + crossValue], c=cmap3.to_rgba(cmapVal), label=resAnalysisRow[varParList[0]], ) else: - ax3.plot(s, resAnalysisRow[runoutResType.lower() + crossValue], c=cmapSC(cmapVal)) + ax3.plot(s, resAnalysisRow[resolvedRunoutResType.lower() + crossValue], c=cmapSC(cmapVal)) countSim = countSim + 1 # add colorbar @@ -500,7 +504,7 @@ def visuRunoutStat(rasterTransfo, inputsDF, resAnalysisDF, newRasters, cfgSetup, ax3.set_xlabel("$S_{xy}$ (thalweg) [m]") ax3.set_xlim([s.min(), s.max()]) ax3.set_ylim(auto=True) - ax3.set_ylabel("$%s_{%s}$ [%s]" % (runoutResType, crossValue, unit)) + ax3.set_ylabel("$%s_{%s}$ [%s]" % (displayRunoutResType, crossValue, unit)) outFileName = "_".join( [projectName, runoutResType, str(thresholdValue).replace(".", "p"), "slComparisonStat"] @@ -689,8 +693,11 @@ def visuComparison(rasterTransfo, inputs, pathDict): simName = inputs['simName'] refSimName = pathDict['refSimName'] runoutResType = inputs['runoutResType'] - unit = pU.cfgPlotUtils['unit' + runoutResType] - name = pU.cfgPlotUtils['name' + runoutResType] + # base name without layer suffix for cosmetic lookups (units, names, colorMaps) + baseRunoutResType = inputs.get('baseRunoutResType', runoutResType) + displayRunoutResType = inputs.get('displayRunoutResType', runoutResType) + unit = pU.cfgPlotUtils['unit' + baseRunoutResType] + name = pU.cfgPlotUtils['name' + baseRunoutResType] thresholdArray = inputs['thresholdArray'] thresholdValue = thresholdArray[-1] @@ -702,7 +709,7 @@ def visuComparison(rasterTransfo, inputs, pathDict): outFilePath = None else: cmapTF, _, ticks, normTF = pU.makeColorMap( - pU.colorMaps[runoutResType], np.nanmin((refData)), np.nanmax((refData)), continuous=pU.contCmap + pU.colorMaps[baseRunoutResType], np.nanmin((refData)), np.nanmax((refData)), continuous=pU.contCmap ) cmapTF.set_bad(color="w") @@ -788,7 +795,7 @@ def visuComparison(rasterTransfo, inputs, pathDict): # add legend associated to the contour plot handles, _ = contourRef.legend_elements() ax2.legend( - title=runoutResType + " contour lines [" + unit + "]", + title=displayRunoutResType + " contour lines [" + unit + "]", handles=handles, labels=labels, ncol=ncolLegend, @@ -862,7 +869,7 @@ def visuComparison(rasterTransfo, inputs, pathDict): pU.putAvaNameOnPlot(ax1, namePrint) ax1.set_title( - "%s difference (sim - reference) in runout area" % runoutResType + "%s difference (sim - reference) in runout area" % displayRunoutResType + "\n" + "blue = FN, red = FP" ) @@ -878,8 +885,8 @@ def visuComparison(rasterTransfo, inputs, pathDict): ax3, insert="False", title=[ - "%s diff histogram" % runoutResType, - "%s diff CDF (95%% and 99%% centiles)" % runoutResType, + "%s diff histogram" % displayRunoutResType, + "%s diff CDF (95%% and 99%% centiles)" % displayRunoutResType, ], ) @@ -887,7 +894,7 @@ def visuComparison(rasterTransfo, inputs, pathDict): divider = make_axes_locatable(ax2) cax = divider.append_axes("right", size="5%", pad=0.1) pU.addColorBar( - im3, ax2, None, None, title=runoutResType + (" [%s]" % unit), extend="both", cax=cax + im3, ax2, None, None, title=displayRunoutResType + (" [%s]" % unit), extend="both", cax=cax ) ax2.set_aspect("equal") ax1.set_aspect("equal") @@ -908,13 +915,13 @@ def visuComparison(rasterTransfo, inputs, pathDict): if pathDict["compType"][0] == "comModules": ax2.set_title( - "%s difference and contour lines" % runoutResType + "%s difference and contour lines" % displayRunoutResType + "\n" + "refMod = full, compMod = dashed line" ) else: ax2.set_title( - "%s difference and contour lines" % runoutResType + "\n" + "ref = full, sim = dashed line" + "%s difference and contour lines" % displayRunoutResType + "\n" + "ref = full, sim = dashed line" ) # fig.subplots_adjust(hspace=0.13, wspace=0.3) @@ -1008,6 +1015,7 @@ def resultVisu(cfgSetup, inputsDF, pathDict, cfgFlags, rasterTransfo, resAnalysi #################################### # Get input data runoutResType = cfgSetup['runoutResType'] + displayRunoutResType = pathDict.get('displayRunoutResType', runoutResType) varParList = cfgSetup['varParList'].split('|') unit = cfgSetup['unit'] paraVar = cfgSetup['varParList'].split('|')[0] @@ -1066,7 +1074,7 @@ def resultVisu(cfgSetup, inputsDF, pathDict, cfgFlags, rasterTransfo, resAnalysi cbar.ax.set_ylabel('hit rate density') else: sc = ax1.scatter(rFP, rTP, c=colorSC, cmap=cmapSC, norm=normSC, marker=pU.markers[0], - label=('runout points (%s<%.1f%s)' % (runoutResType, cfgSetup.getfloat('thresholdValue'), + label=('runout points (%s<%.1f%s)' % (displayRunoutResType, cfgSetup.getfloat('thresholdValue'), pU.cfgPlotUtils['unit' + runoutResType]))) if displayColorBar: @@ -1147,6 +1155,7 @@ def plotContoursTransformed(contourDict, pathDict, rasterTransfo, cfgSetup, inpu l = rasterTransfo['l'] indStartOfRunout = rasterTransfo['indStartOfRunout'] unit = pU.cfgPlotUtils['unit' + cfgSetup['runoutResType']] + displayRunoutResType = pathDict.get('displayRunoutResType', cfgSetup['runoutResType']) colorOrdering = False nSamples = inputsDF.shape[0] @@ -1193,7 +1202,7 @@ def plotContoursTransformed(contourDict, pathDict, rasterTransfo, cfgSetup, inpu # PANEL 1 show contour lines of all sims in contourDict for thresholdValue runoutResType ax1 = fig.add_subplot(gs[0, 0:2]) - ax1.set_title('%s %s %s contour lines' % (cfgSetup['runoutResType'], cfgSetup['thresholdValue'], + ax1.set_title('%s %s %s contour lines' % (displayRunoutResType, cfgSetup['thresholdValue'], unit)) ax1.set_xlabel('$S_{XY}$ (thalweg) [m]') ax1.set_ylabel('$L_{XY}$ (thalweg) [m]') @@ -1837,6 +1846,7 @@ def plotRunoutLineComparisonToReference(cfgSetup, refLine, runoutLine, pathDict, colors = {'line': 'lightcoral', 'point': 'gold','poly': 'purple'} unitResType = pU.cfgPlotUtils['unit' + cfgSetup['runoutResType']] + displayRunoutResType = pathDict.get('displayRunoutResType', cfgSetup['runoutResType']) # create figure fig = plt.figure(figsize=(pU.figW*1.5, pU.figH)) @@ -1845,7 +1855,7 @@ def plotRunoutLineComparisonToReference(cfgSetup, refLine, runoutLine, pathDict, # add panel one with lines and differences across thalweg Lxy ax1 = fig.add_subplot(gs[0:3, 0:2]) ax1.set_title('Runout line (sim %s > %s %s) difference (from reference %s)' % - (cfgSetup['runoutResType'], cfgSetup['thresholdValue'], unitResType, refLine['type'])) + (displayRunoutResType, cfgSetup['thresholdValue'], unitResType, refLine['type'])) ax2 = ax1.twinx() ax1.bar(refLine['l'], runoutLine['s']-refLine['s'], width=5, zorder=1) ax2.plot(runoutLine['l'], runoutLine['s'], label='sim', c='k', zorder=3) diff --git a/avaframe/tests/test_aimecTools.py b/avaframe/tests/test_aimecTools.py index 58be31f68..c3b8e006c 100644 --- a/avaframe/tests/test_aimecTools.py +++ b/avaframe/tests/test_aimecTools.py @@ -276,7 +276,8 @@ def test_computeRunoutLine(tmp_path): # call function to be tested runoutLine = aT.computeRunoutLine( - cfgSetup, rasterTransfo, transformedRasters, "", "line", name="", basedOnMax=False + cfgSetup, rasterTransfo, transformedRasters, "", "line", name="", basedOnMax=False, + runoutResType=cfgSetup["runoutResType"] ) # print('runoutLine', runoutLine) @@ -301,7 +302,8 @@ def test_computeRunoutLine(tmp_path): # call function to be tested runoutLine = aT.computeRunoutLine( - cfgSetup, rasterTransfo, transformedRasters, "", "line", name="", basedOnMax=False + cfgSetup, rasterTransfo, transformedRasters, "", "line", name="", basedOnMax=False, + runoutResType=cfgSetup["runoutResType"] ) # print('runoutLine', runoutLine) @@ -414,3 +416,224 @@ def test_findIntSectCoors(): # Verify no points are marked assert np.sum(intPointsArray) == 0 + + +def _makeAimecCfgSetup(runoutResType="ppr", resTypes="ppr|pft|pfv", runoutLayer=""): + """Create a minimal configparser section mimicking AIMECSETUP for checkAIMECinputs tests""" + cfg = configparser.ConfigParser() + cfg["AIMECSETUP"] = { + "runoutResType": runoutResType, + "resTypes": resTypes, + "runoutLayer": runoutLayer, + } + return cfg["AIMECSETUP"] + + +def test_checkAIMECinputs_singleLayer_unchanged(): + """Single-layer data with no runoutLayer: existing behavior preserved""" + cfgSetup = _makeAimecCfgSetup(runoutResType="ppr", resTypes="ppr|pft|pfv", runoutLayer="") + pathDict = {"resTypeList": ["ppr", "pft", "pfv"]} + inputsDF = pd.DataFrame({ + "simName": ["sim1", "sim2"], + "ppr": ["/path/ppr1.asc", "/path/ppr2.asc"], + "pft": ["/path/pft1.asc", "/path/pft2.asc"], + "pfv": ["/path/pfv1.asc", "/path/pfv2.asc"], + }) + + result = aT.checkAIMECinputs(cfgSetup, pathDict, inputsDF) + + assert sorted(result["resTypeList"]) == ["pft", "pfv", "ppr"] + assert result["runoutResType"] == "ppr" + assert result.get("runoutLayer", "") == "" + + +def test_checkAIMECinputs_multiLayer_withRunoutLayer(): + """Multi-layer data with runoutLayer=L1: resTypeList contains base names""" + cfgSetup = _makeAimecCfgSetup(runoutResType="ppr", resTypes="ppr|pft|pfv", runoutLayer="L1") + pathDict = {"resTypeList": ["ppr_l1", "ppr_l2", "pft_l1", "pft_l2", "pfv_l1", "pfv_l2"]} + inputsDF = pd.DataFrame({ + "simName": ["sim1"], + "ppr_l1": ["/path/L1_ppr.asc"], + "ppr_l2": ["/path/L2_ppr.asc"], + "pft_l1": ["/path/L1_pft.asc"], + "pft_l2": ["/path/L2_pft.asc"], + "pfv_l1": ["/path/L1_pfv.asc"], + "pfv_l2": ["/path/L2_pfv.asc"], + }) + + result = aT.checkAIMECinputs(cfgSetup, pathDict, inputsDF) + + # resTypeList should contain base names, not layer-suffixed + assert sorted(result["resTypeList"]) == ["pft", "pfv", "ppr"] + # runoutResType should be base name + assert result["runoutResType"] == "ppr" + assert result["runoutLayer"] == "L1" + + +def test_checkAIMECinputs_multiLayer_noRunoutLayer_errors(): + """Multi-layer data without runoutLayer: must error out""" + cfgSetup = _makeAimecCfgSetup(runoutResType="ppr", resTypes="ppr|pft|pfv", runoutLayer="") + pathDict = {"resTypeList": ["ppr_l1", "ppr_l2", "pft_l1", "pft_l2", "pfv_l1", "pfv_l2"]} + inputsDF = pd.DataFrame({ + "simName": ["sim1"], + "ppr_l1": ["/path/L1_ppr.asc"], + "ppr_l2": ["/path/L2_ppr.asc"], + "pft_l1": ["/path/L1_pft.asc"], + "pft_l2": ["/path/L2_pft.asc"], + "pfv_l1": ["/path/L1_pfv.asc"], + "pfv_l2": ["/path/L2_pfv.asc"], + }) + + with pytest.raises(FileNotFoundError, match="Multi-layer result files detected.*runoutLayer"): + aT.checkAIMECinputs(cfgSetup, pathDict, inputsDF) + + +def test_checkAIMECinputs_multiLayer_runoutLayer_L2(): + """Multi-layer data with runoutLayer=L2: base names in resTypeList""" + cfgSetup = _makeAimecCfgSetup(runoutResType="ppr", resTypes="ppr|pfv", runoutLayer="L2") + pathDict = {"resTypeList": ["ppr_l1", "ppr_l2", "pfv_l1", "pfv_l2"]} + inputsDF = pd.DataFrame({ + "simName": ["sim1"], + "ppr_l1": ["/path/L1_ppr.asc"], + "ppr_l2": ["/path/L2_ppr.asc"], + "pfv_l1": ["/path/L1_pfv.asc"], + "pfv_l2": ["/path/L2_pfv.asc"], + }) + + result = aT.checkAIMECinputs(cfgSetup, pathDict, inputsDF) + + assert sorted(result["resTypeList"]) == ["pfv", "ppr"] + assert result["runoutResType"] == "ppr" + assert result["runoutLayer"] == "L2" + assert result["displayRunoutResType"] == "ppr (if multilayer: L2)" + + +def test_checkAIMECinputs_mixed_modules(): + """Mixed single-layer + multi-layer: base names resolvable for all sims""" + cfgSetup = _makeAimecCfgSetup(runoutResType="ppr", resTypes="ppr|pfv", runoutLayer="L2") + pathDict = {"resTypeList": []} # empty from makeSimFromResDF in mixed case + inputsDF = pd.DataFrame({ + "simName": ["sim_com1", "sim_com8"], + "ppr": ["/path/com1_ppr.asc", np.nan], + "pfv": ["/path/com1_pfv.asc", np.nan], + "ppr_l1": [np.nan, "/path/L1_ppr.asc"], + "ppr_l2": [np.nan, "/path/L2_ppr.asc"], + "pfv_l1": [np.nan, "/path/L1_pfv.asc"], + "pfv_l2": [np.nan, "/path/L2_pfv.asc"], + }) + + result = aT.checkAIMECinputs(cfgSetup, pathDict, inputsDF) + + assert sorted(result["resTypeList"]) == ["pfv", "ppr"] + assert result["runoutResType"] == "ppr" + assert result["runoutLayer"] == "L2" + + +# --- resolveResTypeColumn tests --- + + +def test_resolveResTypeColumn_singleLayer(): + """Single-layer sim: returns base column name""" + row = pd.Series({"ppr": "/path/to/ppr.asc", "pfv": "/path/to/pfv.asc"}) + assert aT.resolveResTypeColumn(row, "ppr") == "ppr" + assert aT.resolveResTypeColumn(row, "pfv") == "pfv" + + +def test_resolveResTypeColumn_multiLayer_withLayer(): + """Multi-layer sim with layer set: returns layer-suffixed column""" + row = pd.Series({"ppr_l1": "/path/to/L1_ppr.asc", "ppr_l2": "/path/to/L2_ppr.asc"}) + assert aT.resolveResTypeColumn(row, "ppr", layer="L2") == "ppr_l2" + assert aT.resolveResTypeColumn(row, "ppr", layer="L1") == "ppr_l1" + + +def test_resolveResTypeColumn_singleLayer_withLayer_fallback(): + """Single-layer sim with layer set but no layer columns: falls back to base""" + row = pd.Series({"ppr": "/path/to/ppr.asc", "pfv": "/path/to/pfv.asc"}) + assert aT.resolveResTypeColumn(row, "ppr", layer="L2") == "ppr" + + +def test_resolveResTypeColumn_missing_returns_none(): + """No matching column at all: returns None""" + row = pd.Series({"pfv": "/path/to/pfv.asc"}) + assert aT.resolveResTypeColumn(row, "ppr", layer="L2") is None + + +def test_resolveResTypeColumn_nan_skipped(): + """Column exists but is NaN: treated as missing""" + row = pd.Series({"ppr": np.nan, "ppr_l2": "/path/to/L2_ppr.asc"}) + assert aT.resolveResTypeColumn(row, "ppr", layer="L2") == "ppr_l2" + # base column is NaN, no layer set — returns None + row2 = pd.Series({"ppr": np.nan}) + assert aT.resolveResTypeColumn(row2, "ppr") is None + + +# --- computeBaseResTypeList tests --- + + +def test_computeBaseResTypeList_singleLayerOnly(): + """All sims are single-layer: returns base resType names""" + inputsDF = pd.DataFrame({ + "simName": ["sim1", "sim2"], + "simHash": ["abc", "def"], + "layers": [np.nan, np.nan], + "cellSize": [5, 5], + "ppr": ["/path/sim1_ppr.asc", "/path/sim2_ppr.asc"], + "pfv": ["/path/sim1_pfv.asc", "/path/sim2_pfv.asc"], + "pft": ["/path/sim1_pft.asc", "/path/sim2_pft.asc"], + }) + result = aT.computeBaseResTypeList(inputsDF, runoutLayer="") + assert sorted(result) == ["pft", "pfv", "ppr"] + + +def test_computeBaseResTypeList_multiLayerOnly(): + """All sims are multi-layer: returns base names resolvable via runoutLayer""" + inputsDF = pd.DataFrame({ + "simName": ["sim1"], + "simHash": ["abc"], + "layers": ["L1|L2"], + "cellSize": [5], + "ppr_l1": ["/path/L1_ppr.asc"], + "ppr_l2": ["/path/L2_ppr.asc"], + "pfv_l1": ["/path/L1_pfv.asc"], + "pfv_l2": ["/path/L2_pfv.asc"], + }) + result = aT.computeBaseResTypeList(inputsDF, runoutLayer="L2") + assert sorted(result) == ["pfv", "ppr"] + + +def test_computeBaseResTypeList_mixed(): + """Mixed single-layer and multi-layer: returns only base names resolvable for ALL sims""" + inputsDF = pd.DataFrame({ + "simName": ["sim_com1", "sim_com8"], + "simHash": ["abc", "def"], + "layers": [np.nan, "L1|L2"], + "cellSize": [5, 5], + "ppr": ["/path/com1_ppr.asc", np.nan], + "pfv": ["/path/com1_pfv.asc", np.nan], + "pft": ["/path/com1_pft.asc", np.nan], + "ppr_l1": [np.nan, "/path/L1_ppr.asc"], + "ppr_l2": [np.nan, "/path/L2_ppr.asc"], + "pfv_l1": [np.nan, "/path/L1_pfv.asc"], + "pfv_l2": [np.nan, "/path/L2_pfv.asc"], + "pfd_l1": [np.nan, "/path/L1_pfd.asc"], + "pfd_l2": [np.nan, "/path/L2_pfd.asc"], + }) + result = aT.computeBaseResTypeList(inputsDF, runoutLayer="L2") + # ppr and pfv are resolvable for both sims; pft and pfd are not + assert sorted(result) == ["pfv", "ppr"] + + +def test_computeBaseResTypeList_mixed_noLayer_errors(): + """Mixed modules without runoutLayer: multi-layer sims can't resolve base names""" + inputsDF = pd.DataFrame({ + "simName": ["sim_com1", "sim_com8"], + "simHash": ["abc", "def"], + "layers": [np.nan, "L1|L2"], + "cellSize": [5, 5], + "ppr": ["/path/com1_ppr.asc", np.nan], + "ppr_l1": [np.nan, "/path/L1_ppr.asc"], + "ppr_l2": [np.nan, "/path/L2_ppr.asc"], + }) + # Without runoutLayer, multi-layer sim can't resolve 'ppr' (base column is NaN) + result = aT.computeBaseResTypeList(inputsDF, runoutLayer="") + assert result == [] diff --git a/avaframe/tests/test_ana3AIMEC.py b/avaframe/tests/test_ana3AIMEC.py index 2ac45ef43..a7851060e 100644 --- a/avaframe/tests/test_ana3AIMEC.py +++ b/avaframe/tests/test_ana3AIMEC.py @@ -62,6 +62,8 @@ def test_analyzeArea(tmp_path): pathDict["compType"] = ["singleModule", "com1DFA"] pathDict["contCmap"] = True pathDict["resTypeList"] = ["ppr", "pft", "pfv"] + pathDict["runoutResType"] = "ppr" + pathDict["displayRunoutResType"] = "ppr" cfg = cfgUtils.getModuleConfig(ana3AIMEC) cfgSetup = cfg["AIMECSETUP"] @@ -189,6 +191,8 @@ def test_makeDomainTransfo(tmp_path): pathDict["refSimName"] = "testAimec_0" pathDict["compType"] = ["singleModule", "com1DFA"] pathDict["resTypeList"] = ["ppr", "pft", "pfv"] + pathDict["runoutResType"] = "ppr" + pathDict["displayRunoutResType"] = "ppr" cfg = cfgUtils.getModuleConfig(ana3AIMEC) cfgSetup = cfg["AIMECSETUP"] @@ -302,6 +306,8 @@ def test_mainAIMEC(tmp_path): pathDict["refSimName"] = "testAimec_0" pathDict["compType"] = ["singleModule", "com1DFA"] pathDict["resTypeList"] = ["ppr", "pft", "pfv"] + pathDict["runoutResType"] = "ppr" + pathDict["displayRunoutResType"] = "ppr" pathDict["colorParameter"] = False cfg = cfgUtils.getModuleConfig(ana3AIMEC, onlyDefault=True) diff --git a/avaframe/tests/test_cfgUtils.py b/avaframe/tests/test_cfgUtils.py index 8c5e232a2..c63a45069 100644 --- a/avaframe/tests/test_cfgUtils.py +++ b/avaframe/tests/test_cfgUtils.py @@ -1227,4 +1227,112 @@ def test_parseSimName_realWorld_examples(): assert result3["defID"] == "D" assert result3["frictIndi"] == "M" assert result3["simType"] == "ent" - assert result3["modelType"] == "dfa" \ No newline at end of file + assert result3["modelType"] == "dfa" + + +# --- Multi-layer parseSimName tests --- + + +def test_parseSimName_layer_newFormat_resTypeOnly(): + """Test parsing new format with layer and resType, no timeStep""" + name = "release1_a1b2c3_com8_C_null_psa_L1_ppr" + result = cfgUtils.parseSimName(name) + + assert result["releaseName"] == "release1" + assert result["simHash"] == "a1b2c3" + assert result["modName"] == "com8" + assert result["defID"] == "C" + assert result["frictIndi"] is None + assert result["simType"] == "null" + assert result["modelType"] == "psa" + assert result["layer"] == "L1" + assert result["resType"] == "ppr" + assert result["timeStep"] is None + + +def test_parseSimName_layer_newFormat_withTimeStep(): + """Test parsing new format with layer, resType, and timeStep""" + name = "release1_a1b2c3_com8_C_null_psa_L2_pfv_t12.50" + result = cfgUtils.parseSimName(name) + + assert result["layer"] == "L2" + assert result["resType"] == "pfv" + assert result["timeStep"] == "t12.50" + assert result["modName"] == "com8" + assert result["modelType"] == "psa" + + +def test_parseSimName_layer_newFormat_allOptional(): + """Test parsing new format with defID, frictIndi, layer, resType, timeStep""" + name = "release1_a1b2c3_com8_D_S_ent_psa_L3_pft_t50.00" + result = cfgUtils.parseSimName(name) + + assert result["defID"] == "D" + assert result["frictIndi"] == "S" + assert result["simType"] == "ent" + assert result["modelType"] == "psa" + assert result["layer"] == "L3" + assert result["resType"] == "pft" + assert result["timeStep"] == "t50.00" + + +def test_parseSimName_layer_AF_separator(): + """Test parsing with _AF_ separator and layer""" + name = "myRelease_Test_AF_a1b2c3_com8_C_null_psa_L1_ppr" + result = cfgUtils.parseSimName(name) + + assert result["releaseName"] == "myRelease_Test" + assert result["simHash"] == "a1b2c3" + assert result["modName"] == "com8" + assert result["layer"] == "L1" + assert result["resType"] == "ppr" + + +def test_parseSimName_layer_multiDigit(): + """Test parsing with multi-digit layer number""" + name = "release1_a1b2c3_com8_C_null_psa_L12_ppr" + result = cfgUtils.parseSimName(name) + + assert result["layer"] == "L12" + assert result["resType"] == "ppr" + + +def test_parseSimName_noLayer_backwardCompat(): + """Test that existing names without layer still parse correctly with layer=None""" + # New format without layer + name1 = "release1_a1b2c3_com1_C_S_ent_dfa_ppr_t50.2" + result1 = cfgUtils.parseSimName(name1) + assert result1["layer"] is None + assert result1["resType"] == "ppr" + assert result1["timeStep"] == "t50.2" + + # New format, resType only + name2 = "release1_a1b2c3_com1_C_ent_dfa_pfv" + result2 = cfgUtils.parseSimName(name2) + assert result2["layer"] is None + assert result2["resType"] == "pfv" + assert result2["timeStep"] is None + + # Old format without layer + name3 = "release1_a1b2c3_C_S_ent_dfa_ppr" + result3 = cfgUtils.parseSimName(name3) + assert result3["layer"] is None + assert result3["resType"] == "ppr" + + # Bare simName (no resType, no layer) + name4 = "release1_a1b2c3_com1_C_null_dfa" + result4 = cfgUtils.parseSimName(name4) + assert result4["layer"] is None + assert result4["resType"] is None + + +def test_parseSimName_layer_onlyNoResType(): + """Test parsing with layer but no resType — layer appears at offset+2 position""" + # This is an edge case: simName_L1 with no resType after it + # The parser should recognize L1 as a layer, leaving resType=None + name = "release1_a1b2c3_com8_C_null_psa_L1" + result = cfgUtils.parseSimName(name) + + assert result["layer"] == "L1" + assert result["resType"] is None + assert result["timeStep"] is None \ No newline at end of file diff --git a/avaframe/tests/test_dfa2Aimec.py b/avaframe/tests/test_dfa2Aimec.py index b9de0259d..f1f5468c3 100644 --- a/avaframe/tests/test_dfa2Aimec.py +++ b/avaframe/tests/test_dfa2Aimec.py @@ -53,7 +53,8 @@ def test_dfaComp2Aimec(tmp_path): pathDataComp = testPath / 'Outputs' / 'com1DFAComp' / 'peakFiles' cfg = configparser.ConfigParser() cfg['AIMECSETUP'] = {'comModules': 'com1DFARef|com1DFAComp', 'testName': '', 'referenceSimValue': '', - 'referenceSimName': '', 'varParList': ''} + 'referenceSimName': '', 'varParList': '', + 'runoutResType': 'ppr', 'resTypes': 'ppr|pft|pfv', 'runoutLayer': ''} cfg['FLAGS'] = {'flagMass': 'True'} inputDF, pathDict = dfa2Aimec.dfaBench2Aimec(testPath, cfg, 'release1HX', 'release1HX') @@ -150,7 +151,8 @@ def test_dfaBench2Aimec(): # setup required input cfg = configparser.ConfigParser() cfg['AIMECSETUP'] = {'comModules': 'com1DFARef|com1DFAComp', 'testName': '', 'referenceSimValue': '', - 'referenceSimName': '', 'varParList': ''} + 'referenceSimName': '', 'varParList': '', + 'runoutResType': 'ppr', 'resTypes': 'ppr|pft|pfv', 'runoutLayer': ''} cfg['FLAGS'] = {'flagMass': 'True'} dirPath = pathlib.Path(__file__).parents[0] avaTestName = 'avaHelixChannelPytest' diff --git a/avaframe/tests/test_fileHandlerUtils.py b/avaframe/tests/test_fileHandlerUtils.py index 3f129a82e..6fbd75e27 100644 --- a/avaframe/tests/test_fileHandlerUtils.py +++ b/avaframe/tests/test_fileHandlerUtils.py @@ -386,3 +386,111 @@ def test_fetchFlowFields(): assert flowFields[0].stem == 'release1HS_0dcd58fc86_ent_dfa_pft' assert len(flowFields) == 6 + + +# --- Multi-layer makeSimFromResDF tests --- + +RASTER_HEADER = """ncols 10 +nrows 10 +xllcenter 0 +yllcenter 0 +cellsize 5 +NODATA_value -9999 +""" + + +def _createTestRaster(filepath): + """Create a minimal valid raster file for testing""" + with open(filepath, "w") as f: + f.write(RASTER_HEADER) + for _ in range(10): + f.write(" ".join(["0"] * 10) + "\n") + + +def test_makeSimFromResDF_multiLayer(tmp_path): + """Test that multi-layer files produce layer-suffixed columns""" + peakDir = tmp_path / "peakFiles" + peakDir.mkdir() + + # Create multi-layer result files: sim_L1_ppr, sim_L1_pfv, sim_L2_ppr, sim_L2_pfv + simBase = "release1_abc123_com8_C_null_psa" + for layer in ["L1", "L2"]: + for resType in ["ppr", "pfv"]: + _createTestRaster(peakDir / f"{simBase}_{layer}_{resType}.asc") + + dataDF, resTypeListAll = fU.makeSimFromResDF(str(tmp_path), "com8", inputDir=str(peakDir)) + + # Layer-suffixed columns should exist + assert "ppr_l1" in dataDF.columns + assert "ppr_l2" in dataDF.columns + assert "pfv_l1" in dataDF.columns + assert "pfv_l2" in dataDF.columns + + # Layer-suffixed columns should contain file paths (not NaN) + assert not dataDF["ppr_l1"].isnull().any() + assert not dataDF["ppr_l2"].isnull().any() + assert not dataDF["pfv_l1"].isnull().any() + assert not dataDF["pfv_l2"].isnull().any() + + # Layers metadata column should exist and be populated + assert "layers" in dataDF.columns + layersVal = dataDF["layers"].iloc[0] + assert "L1" in layersVal + assert "L2" in layersVal + + # resTypeListAll should contain layer-suffixed names + assert "ppr_l1" in resTypeListAll + assert "ppr_l2" in resTypeListAll + assert "pfv_l1" in resTypeListAll + assert "pfv_l2" in resTypeListAll + + +def test_makeSimFromResDF_singleLayer_unchanged(tmp_path): + """Test that single-layer files still produce standard columns (backward compat)""" + peakDir = tmp_path / "peakFiles" + peakDir.mkdir() + + # Create standard single-layer result files + simBase = "release1_abc123_com1_C_null_dfa" + for resType in ["ppr", "pfv", "pft"]: + _createTestRaster(peakDir / f"{simBase}_{resType}.asc") + + dataDF, resTypeListAll = fU.makeSimFromResDF(str(tmp_path), "com1", inputDir=str(peakDir)) + + # Standard columns should exist + assert "ppr" in dataDF.columns + assert "pfv" in dataDF.columns + assert "pft" in dataDF.columns + + # Standard columns should contain file paths + assert not dataDF["ppr"].isnull().any() + assert not dataDF["pfv"].isnull().any() + assert not dataDF["pft"].isnull().any() + + # Layers metadata column should be NaN for single-layer + if "layers" in dataDF.columns: + assert dataDF["layers"].isnull().all() + + # resTypeListAll should contain standard names + assert "ppr" in resTypeListAll + assert "pfv" in resTypeListAll + assert "pft" in resTypeListAll + + +def test_makeSimFromResDF_multiLayer_simName_reconstructed(tmp_path): + """Test that simName is correctly reconstructed without layer component""" + peakDir = tmp_path / "peakFiles" + peakDir.mkdir() + + simBase = "release1_abc123_com8_C_null_psa" + _createTestRaster(peakDir / f"{simBase}_L1_ppr.asc") + _createTestRaster(peakDir / f"{simBase}_L2_ppr.asc") + + dataDF, _ = fU.makeSimFromResDF(str(tmp_path), "com8", inputDir=str(peakDir)) + + # Should have exactly one row (one simulation, two layers) + assert len(dataDF) == 1 + # simName should NOT contain the layer + assert "L1" not in dataDF["simName"].iloc[0] + assert "L2" not in dataDF["simName"].iloc[0] + assert dataDF["simName"].iloc[0] == simBase