[out1] add compute area indicators based on real area not transformed into thalweg following

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)

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)