From 6124541c8934def9ecd947b6e8cc442a65617aaa Mon Sep 17 00:00:00 2001 From: Anna Wirbel Date: Thu, 7 Aug 2025 10:10:45 +0200 Subject: [PATCH] add compute area refactor(out1Peak): improve code readability with consistent formatting and parameter structuring --- avaframe/out1Peak/outPlotAllPeakDiffs.py | 269 +++++++++++++++++++++ avaframe/runScripts/runPlotAreaRefDiffs.py | 113 +++++++++ 2 files changed, 382 insertions(+) create mode 100644 avaframe/out1Peak/outPlotAllPeakDiffs.py create mode 100644 avaframe/runScripts/runPlotAreaRefDiffs.py diff --git a/avaframe/out1Peak/outPlotAllPeakDiffs.py b/avaframe/out1Peak/outPlotAllPeakDiffs.py new file mode 100644 index 000000000..08a7cc85a --- /dev/null +++ b/avaframe/out1Peak/outPlotAllPeakDiffs.py @@ -0,0 +1,269 @@ +""" +Functions to plot 2D simulation results and area indicators + +""" + +import logging +import numpy as np +from matplotlib import pyplot as plt +import pathlib +from mpl_toolkits.axes_grid1 import make_axes_locatable +import avaframe.out3Plot.plotUtils as pU +import geopandas as gpd +import copy +from cmcrameri import cm as cmapCrameri + +# create local logger +log = logging.getLogger(__name__) + + +def computeAreaDiff( + refRaster, + simRaster, + thresholdValueReference, + thresholdValueSimulation, + dem, + cropToArea=None, +): + """compute difference between two rasters based on a thresholdValue + + Parameters + ------------ + refRaster: numpy ndarray + reference raster + simRaster: numpy ndarray + simulation result raster + thresholdValueReference: float + threshold value to mask reference raster + thresholdValueSimlation: float + threshold value to mask simulation raster + dem: dict + dictionary with info on DEM used for sims, here areaRaster is required to compute actual areas of indicators + cropToArea: None or numpy ndarray + None or raster that is used to define where comparison is computed (0 -no, 1-yes) + + Returns + --------- + refMask: numpy ndarray + array with 1 indicating reference data found exceeding threshold and 0 not found + compRaster: numpy ndarray + array with 1 indicating simulation data found exceeding threshold and 0 not found + indicatorDict: dict + dictionary with info on true positive, false positive, false negative + arrays with 1 indicating true positive, false positive, false negative areas and 0 all other entries + number of cells where tp, fp, np + total area for tp, fp, fn + """ + + refMask = copy.deepcopy(refRaster) + # set all nans to 0 values + refMask = np.where(np.isnan(refMask), 0, refMask) + # set to 0 when values <= threshold, and else if threshold exceeded + refMask = np.where(refMask <= thresholdValueReference, 0, refMask) + refMask = np.where(refMask > thresholdValueReference, 1, refMask) + + # set all nans to 0 values + compRasterMask = copy.deepcopy(simRaster) + # set to 0 when values <= threshold, and else if threshold exceeded + compRasterMask = np.where(np.isnan(compRasterMask), 0, compRasterMask) + compRasterMask = np.where(compRasterMask <= thresholdValueSimulation, 0, compRasterMask) + compRasterMask = np.where(compRasterMask > thresholdValueSimulation, 1, compRasterMask) + + # if array to crop data is provided, perform analysis only for this area + if isinstance(cropToArea, np.ndarray): + refMask = np.where(cropToArea == 0, 0, refMask) + compRasterMask = np.where(cropToArea == 0, 0, compRasterMask) + + # create true positive, false positive, false negative areas + tp = np.where((compRasterMask == 1) & (refMask == 1), 1, 0) + fp = np.where((compRasterMask == 1) & (refMask == 0), 1, 0) + fn = np.where((compRasterMask == 0) & (refMask == 1), 1, 0) + + indicatorDict = { + "truePositive": { + "array": tp, + "nCells": np.nansum(tp), + "areaSum": np.nansum(tp * dem["areaRaster"]), + }, + "falsePositive": { + "array": fp, + "nCells": np.nansum(fp), + "areaSum": np.nansum(fp * dem["areaRaster"]), + }, + "falseNegative": { + "array": fn, + "nCells": np.nansum(fn), + "areaSum": np.nansum(fn * dem["areaRaster"]), + }, + } + + return refMask, compRasterMask, indicatorDict + + +def plotAreaDiff( + refRaster, + refMask, + simRaster, + compRasterMask, + resType, + rasterHeader, + thresholdValueSimulation, + outDir, + indicatorDict, + simName, + cropFile=None, +): + """plot comparison of reference and simulation and area indicators + + Parameters + ----------- + refRaster, simRaster: numpy ndarray + reference raster, simulation raster + refMask, compRasterMask: numpy ndarray + array with 1 indicating reference (simulation) data found exceeding threshold and 0 not found + resType: str + result type of simulation + rasterHeader: dict + info on extent + thresholdValueSimulation: float + threshold value to mask simulation raster + outDir: pathlib path + path to save plot + indicatorDict: dict + dictionary with info on true positive, false positive, false negative + arrays with 1 indicating true positive, false positive, false negative areas and 0 all other entries + number of cells where tp, fp, np + total area for tp, fp, fn + simName: str + name of simulation + cropFile: None or pathlib path + None or path to cropfile with polygon used to perfrom analysis + + """ + + data = compRasterMask - refMask + # mask where both don't show values + data = np.ma.masked_where((refMask == 0) & (compRasterMask == 0), data) + + # create figure + fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(pU.figW * 3, pU.figH)) + # fetch extent of domain + extentCellCenters, extentCellCorners = pU.createExtentMinMax( + refRaster, rasterHeader, originLLCenter=True + ) + cmapTF, _, ticks, normTF = pU.makeColorMap(cmapCrameri.hawaii, -0.5, 0.5, continuous=pU.contCmap) + cmapTF.set_under(color="b") + cmapTF.set_over(color="r") + im1 = ax[0].imshow( + data, + cmap=cmapTF, + norm=normTF, + extent=extentCellCorners, + origin="lower", + aspect="equal", + zorder=4, + ) + # add area indicator info + ax[0].text( + 0.01, + 0.85, + ( + "false negative:\n %.0f cells, %.2f m2" + % ( + indicatorDict["falseNegative"]["nCells"], + indicatorDict["falseNegative"]["areaSum"], + ) + ), + horizontalalignment="left", + verticalalignment="center", + transform=ax[0].transAxes, + color="royalblue", + alpha=0.75, + zorder=5, + ) + ax[0].text( + 0.01, + 0.95, + ( + "false positive:\n %.0f cells, %.2f m2" + % ( + indicatorDict["falsePositive"]["nCells"], + indicatorDict["falsePositive"]["areaSum"], + ) + ), + horizontalalignment="left", + verticalalignment="center", + transform=ax[0].transAxes, + color="darkred", + alpha=0.75, + zorder=5, + ) + ax[0].text( + 0.01, + 0.75, + ( + "true positive:\n %.0f cells, %2.f m2" + % ( + indicatorDict["truePositive"]["nCells"], + indicatorDict["truePositive"]["areaSum"], + ) + ), + horizontalalignment="left", + verticalalignment="center", + transform=ax[0].transAxes, + color="darkgreen", + alpha=0.75, + zorder=5, + ) + title = str( + "Affected area based on %s %.2f %s" + % (resType, thresholdValueSimulation, pU.cfgPlotUtils["unit" + resType]) + ) + + if isinstance(cropFile, pathlib.Path): + focus = gpd.read_file(cropFile) + focus.plot( + ax=ax[0], + zorder=12, + edgecolor="darkgrey", + linewidth=2, + facecolor="none", + alpha=1.0, + ) + + # choose colormap for simulation result + cmap1, col1, ticks1, norm1 = pU.makeColorMap( + pU.colorMaps[resType], + np.nanmin(refRaster), + np.nanmax(refRaster), + continuous=True, + ) + cmap1.set_bad(alpha=0) + im2 = ax[1].imshow( + np.ma.masked_where(simRaster == 0.0, simRaster), + cmap=cmap1, + extent=extentCellCorners, + origin="lower", + aspect="equal", + ) + # plot reference raster + ax[2].imshow(refRaster, extent=extentCellCorners, origin="lower", aspect="equal") + + # add colorbar + divider = make_axes_locatable(ax[1]) + cax = divider.append_axes("right", size="5%", pad=0.05) + cbar = fig.colorbar(im2, cax=cax) + cbar.set_label(resType) + + # add labels + ax[0].set_xlabel("x [m]") + ax[0].set_ylabel("y [m]") + ax[0].set_title(title) + ax[1].set_xlabel("x [m]") + ax[1].set_ylabel("y [m]") + ax[1].set_title("Simulation %s" % resType) + ax[2].set_xlabel("x [m]") + ax[2].set_ylabel("y [m]") + ax[2].set_title("Reference mask") + + pU.saveAndOrPlot({"pathResult": outDir}, "areaIndicatorAnalysis_%s" % simName, fig) diff --git a/avaframe/runScripts/runPlotAreaRefDiffs.py b/avaframe/runScripts/runPlotAreaRefDiffs.py new file mode 100644 index 000000000..c5b941cd2 --- /dev/null +++ b/avaframe/runScripts/runPlotAreaRefDiffs.py @@ -0,0 +1,113 @@ +""" + Run script for plotting a comparison of simulation result to reference polygon +""" +# Load modules +# importing general python modules +import pathlib +import numpy as np + +# Local imports +from avaframe.in3Utils import cfgUtils +from avaframe.in3Utils import logUtils +from avaframe.in3Utils import fileHandlerUtils as fU +import avaframe.in2Trans.rasterUtils as IOf +import avaframe.out1Peak.outPlotAllPeakDiffs as oPD +import avaframe.in1Data.getInput as gI +import avaframe.in2Trans.shpConversion as shpConv +import avaframe.in3Utils.geoTrans as gT +import avaframe.com1DFA.DFAtools as DFAtls + + +################USER Input############# +resType = "ppr" +thresholdValueSimulation = 0.9 +modName = 'com1DFA' +############################################################ + +# Load avalanche directory from general configuration file +cfgMain = cfgUtils.getGeneralConfig() +avalancheDir = cfgMain["MAIN"]["avalancheDir"] +outDir = pathlib.Path(avalancheDir, 'Outputs', 'out1Peak') +fU.makeADir(outDir) + +# Start logging +logName = "plotAreaDiff_%s" % (resType) +# Start logging +log = logUtils.initiateLogger(avalancheDir, logName) +log.info("MAIN SCRIPT") +log.info("Current avalanche: %s", avalancheDir) + + +# initialize DEM from avalancheDir (used to perform simulations) +# TODO: if meshCellSize was changed - use actual simulation DEM +dem = gI.readDEM(avalancheDir) +# get normal vector of the grid mesh +dem = gT.getNormalMesh(dem, num=1) +# get real Area +dem = DFAtls.getAreaMesh(dem, 1) +dem['originalHeader'] = dem['header'] + +# read reference data set +inDir = pathlib.Path(avalancheDir, 'Inputs') +referenceFile, availableFile = gI.getAndCheckInputFiles(inDir, 'REFDATA', 'POLY', + fileExt="shp", fileSuffix='POLY') +# convert polygon to raster with value 1 inside polygon and 0 outside the polygon +referenceLine = shpConv.readLine(referenceFile, "reference", dem) +referenceLine= gT.prepareArea(referenceLine, dem, np.sqrt(2),combine=True, checkOverlap=False) + +# if available zoom into area provided by crop shp file in Inputs/CROPSHAPE +cropFile, cropInfo = gI.getAndCheckInputFiles( + inDir, "POLYGONS", "cropFile", fileExt="shp", fileSuffix="_cropshape" + ) +if cropInfo: + cropLine = shpConv.readLine(cropFile, "cropFile", dem) + cropLine = gT.prepareArea(cropLine, dem, np.sqrt(2), combine=True, checkOverlap=False) + +if modName == 'com1DFA': + # load dataFrame for all configurations of simulations in avalancheDir + simDF = cfgUtils.createConfigurationInfo(avalancheDir) + # create data frame that lists all available simulations and path to their result type result files + inputsDF, resTypeList = fU.makeSimFromResDF(avalancheDir, "com1DFA") + # merge parameters as columns to dataDF for matching simNames + dataDF = inputsDF.merge(simDF, left_on="simName", right_on="simName") + + ## loop over all simulations and load desired resType + for index, row in dataDF.iterrows(): + simFile = row[resType] + simData = IOf.readRaster(simFile) + + # compute referenceMask and simulationMask and true positive, false positive and false neg. arrays + # here thresholdValueReference is set to 0.9 as when converting the polygon to a raster, + # values inside polygon are set to 1 and outside to 0 + refMask, compMask, indicatorDict = oPD.computeAreaDiff(referenceLine['rasterData'], + simData['rasterData'], + 0.9, + thresholdValueSimulation, + dem, + cropToArea=cropLine['rasterData']) + + # plot differences + oPD.plotAreaDiff(referenceLine['rasterData'], refMask, simData['rasterData'], compMask, resType, simData['header'], + thresholdValueSimulation, outDir, + indicatorDict, row['simName'], cropFile=cropFile) +else: + # load all result files + resultDir = pathlib.Path(avalancheDir, 'Outputs', modName, 'peakFiles') + peakFilesList = list(resultDir.glob("*_%s.tif" % resType)) + list(resultDir.glob("*_%s.asc" % resType)) + for pF in peakFilesList: + simData = IOf.readRaster(pF) + simName = pF.stem + + # compute referenceMask and simulationMask and true positive, false positive and false neg. arrays + refMask, compMask, indicatorDict = oPD.computeAreaDiff(referenceLine['rasterData'], + simData['rasterData'], + 0.9, + thresholdValueSimulation, + dem, + cropToArea=cropLine['rasterData']) + + # plot differences + oPD.plotAreaDiff(referenceLine['rasterData'], refMask, simData['rasterData'], compMask, resType, + simData['header'], + thresholdValueSimulation, outDir, + indicatorDict, simName, cropFile=cropFile)