From 72151b2fba4a1fdd2e2a05e2bf775cf04530f8b7 Mon Sep 17 00:00:00 2001 From: paulasp Date: Fri, 21 Feb 2025 11:07:41 +0100 Subject: [PATCH 1/2] use raster functions from avaframe --- avaframe/com4FlowPy/com4FlowPy.py | 124 ++++++++----------- avaframe/com4FlowPy/splitAndMerge.py | 5 +- avaframe/runCom4FlowPy.py | 171 +++++++++------------------ docs/moduleCom4FlowPy.rst | 25 +--- 4 files changed, 112 insertions(+), 213 deletions(-) diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py index 83e488597..a9f909400 100755 --- a/avaframe/com4FlowPy/com4FlowPy.py +++ b/avaframe/com4FlowPy/com4FlowPy.py @@ -24,7 +24,6 @@ import avaframe.in3Utils.geoTrans as gT # com4FlowPy Libraries -import avaframe.com4FlowPy.rasterIo as io import avaframe.com4FlowPy.flowCore as fc import avaframe.com4FlowPy.splitAndMerge as SPAM @@ -173,19 +172,9 @@ def com4FlowPyMain(cfgPath, cfgSetup): # get information on cellsize and nodata value from demHeader rasterAttributes = {} - - demExt = os.path.splitext(modelPaths["demPath"])[1] - if demExt in ['.asc', '.ASC']: - demHeader = IOf.readRasterHeader(modelPaths["demPath"]) - rasterAttributes["nodata"] = demHeader["nodata_value"] - elif demExt in ['.tif', '.tiff', '.TIF', '.TIFF']: - demHeader = io.read_header(modelPaths["demPath"]) - rasterAttributes["nodata"] = demHeader["noDataValue"] - else: - _errMsg = f"file format '{demExt}' for DEM is not supported, please provide '.tif' or '.asc'" - log.error(_errMsg) - raise ValueError(_errMsg) - + + demHeader = IOf.readRasterHeader(modelPaths["demPath"]) + rasterAttributes["nodata"] = demHeader["nodata_value"] rasterAttributes["cellsize"] = demHeader["cellsize"] # tile input layers and write tiles (pickled np.arrays) to temp Folder @@ -291,15 +280,8 @@ def checkInputLayerDimensions(modelParameters, modelPaths): ext = os.path.splitext(modelPaths["demPath"])[1] - if ext in ['.asc', '.ASC']: - _demHeader = IOf.readRasterHeader(modelPaths["demPath"]) - _relHeader = io.read_header(modelPaths["releasePathWork"]) - elif ext in ['.tif', '.tiff', '.TIF', '.TIFF']: - _demHeader = io.read_header(modelPaths["demPath"]) - _relHeader = io.read_header(modelPaths["releasePathWork"]) - else: - _errMsg = f"file format '{ext}' for DEM is not supported, please provide '.tif' or '.asc'" - raise ValueError(_errMsg) + _demHeader = IOf.readRasterHeader(modelPaths["demPath"]) + _relHeader = IOf.readRasterHeader(modelPaths["releasePathWork"]) if _demHeader["ncols"] == _relHeader["ncols"] and _demHeader["nrows"] == _relHeader["nrows"]: log.info("DEM and Release Layer ok!") @@ -308,7 +290,7 @@ def checkInputLayerDimensions(modelParameters, modelPaths): sys.exit(1) if modelParameters["infraBool"]: - _infraHeader = io.read_header(modelPaths["infraPath"]) + _infraHeader = IOf.readRasterHeader(modelPaths["infraPath"]) if _demHeader["ncols"] == _infraHeader["ncols"] and _demHeader["nrows"] == _infraHeader["nrows"]: log.info("Infra Layer ok!") else: @@ -316,7 +298,7 @@ def checkInputLayerDimensions(modelParameters, modelPaths): sys.exit(1) if modelParameters["forestBool"]: - _forestHeader = io.read_header(modelPaths["forestPath"]) + _forestHeader = IOf.readRasterHeader(modelPaths["forestPath"]) if _demHeader["ncols"] == _forestHeader["ncols"] and _demHeader["nrows"] == _forestHeader["nrows"]: log.info("Forest Layer ok!") else: @@ -324,7 +306,7 @@ def checkInputLayerDimensions(modelParameters, modelPaths): sys.exit(1) if modelParameters["varUmaxBool"]: - _varUmaxHeader = io.read_header(modelPaths["varUmaxPath"]) + _varUmaxHeader = IOf.readRasterHeader(modelPaths["varUmaxPath"]) if _demHeader["ncols"] == _varUmaxHeader["ncols"] and _demHeader["nrows"] == _varUmaxHeader["nrows"]: log.info("uMax Limit Layer ok!") else: @@ -332,7 +314,7 @@ def checkInputLayerDimensions(modelParameters, modelPaths): sys.exit(1) if modelParameters["varAlphaBool"]: - _varAlphaHeader = io.read_header(modelPaths["varAlphaPath"]) + _varAlphaHeader = IOf.readRasterHeader(modelPaths["varAlphaPath"]) if _demHeader["ncols"] == _varAlphaHeader["ncols"] and _demHeader["nrows"] == _varAlphaHeader["nrows"]: log.info("variable Alpha Layer ok!") else: @@ -340,7 +322,7 @@ def checkInputLayerDimensions(modelParameters, modelPaths): sys.exit(1) if modelParameters["varExponentBool"]: - _varExponentHeader = io.read_header(modelPaths["varExponentPath"]) + _varExponentHeader = IOf.readRasterHeader(modelPaths["varExponentPath"]) if _demHeader["ncols"] == _varExponentHeader["ncols"] and _demHeader["nrows"] == _varExponentHeader["nrows"]: log.info("variable exponent Layer ok!") else: @@ -385,7 +367,8 @@ def checkInputParameterValues(modelParameters, modelPaths): _checkVarParams = True if modelParameters["varAlphaBool"]: - rasterValues, header = io.read_raster(modelPaths["varAlphaPath"]) + data = IOf.readRaster(modelPaths["varAlphaPath"]) + rasterValues = data["rasterData"] rasterValues[rasterValues < 0] = np.nan # handle different noData values if np.any(rasterValues > 90, where=~np.isnan(rasterValues)): log.error("Error: Not all Alpha-raster values are within a physically sensible range ([0,90]),\ @@ -393,7 +376,8 @@ def checkInputParameterValues(modelParameters, modelPaths): _checkVarParams = False if modelParameters["varUmaxBool"]: - rasterValues, header = io.read_raster(modelPaths["varUmaxPath"]) + data = IOf.readRaster(modelPaths["varUmaxPath"]) + rasterValues = data["rasterData"] rasterValues[rasterValues < 0] = np.nan if modelParameters["varUmaxType"].lower() == 'umax': _maxVal = 1500 # ~sqrt(8848*2*9.81) @@ -549,44 +533,49 @@ def mergeAndWriteResults(modelPaths, modelOptions): _oF = modelPaths["outputFileFormat"] _ts = modelPaths["timeString"] + demHeader = IOf.readRasterHeader(modelPaths["demPath"]) + # TODO: should we keep the modelPaths["outputFileFormat"] or can we use the DEM - format for the outputs + # TODO: do we want -99909 as nodata values or the DEM - nodata + if 'flux' in _outputs: - io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / "com4_{}_{}_flux{}".format(_uid, _ts, _oF), - flux) + output = IOf.writeResultToRaster(demHeader, flux, + modelPaths["resDir"] / "com4_{}_{}_flux".format(_uid, _ts), flip=True) if 'zDelta' in _outputs: - io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / "com4_{}_{}_zdelta{}".format(_uid, _ts, _oF), - zDelta) + output = IOf.writeResultToRaster(demHeader, zDelta, + modelPaths["resDir"] / "com4_{}_{}_zdelta".format(_uid, _ts), flip=True) if 'cellCounts' in _outputs: - io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / "com4_{}_{}_cellCounts{}".format(_uid, _ts, _oF), - cellCounts) + output = IOf.writeResultToRaster(demHeader, cellCounts, + modelPaths["resDir"] / "com4_{}_{}_cellCounts".format(_uid, _ts), flip=True) if 'zDeltaSum' in _outputs: - io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / "com4_{}_{}_zDeltaSum{}".format(_uid, _ts, _oF), - zDeltaSum) + output = IOf.writeResultToRaster(demHeader, zDeltaSum, + modelPaths["resDir"] / "com4_{}_{}_zDeltaSum".format(_uid, _ts), flip=True) if 'routFluxSum' in _outputs: - io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / "com4_{}_{}_routFluxSum{}".format(_uid, _ts, _oF), - routFluxSum) + output = IOf.writeResultToRaster(demHeader, routFluxSum, + modelPaths["resDir"] / "com4_{}_{}_routFluxSum".format(_uid, _ts), flip=True) if 'depFluxSum' in _outputs: - io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / "com4_{}_{}_depFluxSum{}".format(_uid, _ts, _oF), - depFluxSum) + output = IOf.writeResultToRaster(demHeader, depFluxSum, + modelPaths["resDir"] / "com4_{}_{}_depFluxSum".format(_uid, _ts), flip=True) if 'fpTravelAngle' in _outputs: - io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / "com4_{}_{}_fpTravelAngle{}".format(_uid, - _ts, _oF), fpTa) + output = IOf.writeResultToRaster(demHeader, fpTa, + modelPaths["resDir"] / "com4_{}_{}_fpTravelAngle".format(_uid, _ts), flip=True) if 'slTravelAngle' in _outputs: - io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / "com4_{}_{}_slTravelAngle{}".format(_uid, - _ts, _oF), slTa) + output = IOf.writeResultToRaster(demHeader, slTa, + modelPaths["resDir"] / "com4_{}_{}_slTravelAngle".format(_uid, _ts), flip=True) if 'travelLength' in _outputs: - io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / "com4_{}_{}_travelLength{}".format(_uid, - _ts, _oF), travelLength) + output = IOf.writeResultToRaster(demHeader, travelLength, + modelPaths["resDir"] / "com4_{}_{}_travelLength".format(_uid, _ts), flip=True) # NOTE: # if not modelOptions["infraBool"]: # if no infra # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("cell_counts%s" %(output_format)),cell_counts) # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("z_delta_sum%s" %(output_format)),z_delta_sum) if modelOptions["infraBool"]: # if infra - io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / "com4_{}_{}_backcalculation{}".format(_uid, - _ts, _oF), backcalc) + output = IOf.writeResultToRaster(demHeader, backcalc, + modelPaths["resDir"] / "com4_{}_{}_backcalculation".format(_uid, _ts), flip=True) if modelOptions["forestInteraction"]: - io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / "com4_{}_{}_forestInteraction{}".format(_uid, - _ts, _oF), forestInteraction) + output = IOf.writeResultToRaster(demHeader, forestInteraction, + modelPaths["resDir"] / "com4_{}_{}_forestInteraction".format(_uid, _ts), flip=True) + del output def checkConvertReleaseShp2Tif(modelPaths): @@ -604,34 +593,25 @@ def checkConvertReleaseShp2Tif(modelPaths): """ # the release is a shp polygon, we need to convert it to a raster # releaseLine = shpConv.readLine(releasePath, 'releasePolygon', demDict) - if modelPaths["releasePath"].suffix == ".shp": - ext = os.path.splitext(modelPaths["demPath"])[1] - - if ext in ['.asc', '.ASC']: - dem = IOf.readRaster(modelPaths["demPath"]) - demHeader = IOf.readRasterHeader(modelPaths["demPath"]) - elif ext in ['.tif', '.tiff', '.TIF', '.TIFF']: - _errMsg = "using release area in '.shp' format currently only supported in combination with '.asc' DEMs" - log.error(_errMsg) - raise ValueError(_errMsg) - else: - _errMsg = f"file format '{ext}' for DEM is not supported, please provide '.tif' or '.asc'" - log.error(_errMsg) - raise ValueError(_errMsg) - - dem['originalHeader'] = demHeader + dem = IOf.readRaster(modelPaths["demPath"]) + demHeader = dem["header"] + dem['originalHeader'] = demHeader # TODO: why? releaseLine = shpConv.SHP2Array(modelPaths["releasePath"], "releasePolygon") thresholdPointInPoly = 0.01 releaseLine = gT.prepareArea(releaseLine, dem, thresholdPointInPoly, combine=True, checkOverlap=False) # give the same header as the dem - releaseAreaHeader = demHeader releaseArea = np.flipud(releaseLine["rasterData"]) - modelPaths["releasePathWork"] = modelPaths["workDir"] / "release.tif" - io.output_raster(modelPaths["demPath"], modelPaths["workDir"] / "release.asc", releaseArea) - io.output_raster(modelPaths["demPath"], modelPaths["workDir"] / "release.tif", releaseArea) + if demHeader["driver"] == "AAIGrid": + modelPaths["releasePathWork"] = modelPaths["workDir"] / "release.asc" + if demHeader["driver"] == "GTiff": + modelPaths["releasePathWork"] = modelPaths["workDir"] / "release.tif" + # now the release area is only saved in the same format of the DEM and nodata values + # are the same as in the DEM (instead of -9999) + # TODO: check if this is a problem explicitly + result = IOf.writeResultToRaster(demHeader, releaseArea, modelPaths["workDir"] / "release") del releaseLine else: modelPaths["releasePathWork"] = modelPaths["releasePath"] diff --git a/avaframe/com4FlowPy/splitAndMerge.py b/avaframe/com4FlowPy/splitAndMerge.py index e621d6241..46486d4b6 100644 --- a/avaframe/com4FlowPy/splitAndMerge.py +++ b/avaframe/com4FlowPy/splitAndMerge.py @@ -10,7 +10,7 @@ import pickle import gc import numpy as np -import avaframe.com4FlowPy.rasterIo as io +import avaframe.in2Trans.rasterUtils as IOf # create local logger log = logging.getLogger(__name__) @@ -42,7 +42,8 @@ def tileRaster(fNameIn, fNameOut, dirName, xDim, yDim, U, isInit=False): # os.makedirs(dirName) # largeRaster, largeHeader = iof.f_readASC(fNameIn, dType='float') - largeRaster, largeHeader = io.read_raster(fNameIn) + largeData = IOf.readRaster(fNameIn) + largeRaster = largeData["rasterData"] # einlesen des Rasters und der Header i, j, imax, jmax = 0, 0, 0, 0 diff --git a/avaframe/runCom4FlowPy.py b/avaframe/runCom4FlowPy.py index 67162e634..8fd93ded2 100644 --- a/avaframe/runCom4FlowPy.py +++ b/avaframe/runCom4FlowPy.py @@ -201,149 +201,86 @@ def readFlowPyinputs(avalancheDir, cfgFlowPy, log): cfgPath = {} avalancheDir = pathlib.Path(avalancheDir) - # read release area - releaseDir = avalancheDir / "Inputs" / "REL" - - # from shapefile - relFiles = sorted(list(releaseDir.glob("*.shp"))) - tryTif = False - if len(relFiles) == 0: - log.info("Found no *.shp file containing the release area in %s, trying with *.tif" % releaseDir) - tryTif = True - elif len(relFiles) > 1: - message = "There should be exactly one *.shp file containing the release area in %s" % releaseDir + + inputDir = avalancheDir / "Inputs" + relFile, available = gI.getAndCheckInputFiles(inputDir, "REL", "Release Area", fileExt="shp") + + if available == "No": + relFile, available = gI.getAndCheckInputFiles(inputDir, "REL", "Release Area", fileExt="raster") + if available == "No": + message = f"There is no release area file in supported format provided in {avalancheDir}/REL" log.error(message) raise AssertionError(message) - else: - log.info("Release area file is: %s" % relFiles[0]) - cfgPath["releasePath"] = relFiles[0] - - # from tif - if tryTif: - relFiles = sorted(list(releaseDir.glob("*.tif"))) - if len(relFiles) == 0: - message = ( - "You need to provide one *.shp file or one *.tif file containing the release area in %s" - % releaseDir - ) - log.error(message) - raise FileNotFoundError(message) - elif len(relFiles) > 1: - message = "There should be exactly one *.tif file containing the release area in %s" % releaseDir + log.info("Release area file is: %s" % relFile) + cfgPath["releasePath"] = relFile + + # TODO: also use the getAndCheckInputFiles to get the paths for the following files? + # read infra area + if cfgFlowPy.getboolean("GENERAL", "infra") is True: + infraPath, available = gI.getAndCheckInputFiles(inputDir, "INFRA", "Infra", fileExt="raster") + if available == "No": + message = f"There is no infra file in supported format provided in {avalancheDir}/INFRA" log.error(message) raise AssertionError(message) - else: - log.info("Release area file is: %s" % relFiles[0]) - cfgPath["releasePath"] = relFiles[0] - - # read infra area - infraDir = avalancheDir / "Inputs" / "INFRA" - infraPath = sorted(list(infraDir.glob("*.tif"))) - if len(infraPath) == 0 or cfgFlowPy.getboolean("GENERAL", "infra") is False: - infraPath = "" - elif len(infraPath) > 1: - message = "More than one Infrastructure file .%s file in %s not allowed" % (infraDir) - log.error(message) - raise AssertionError(message) - else: - infraPath = infraPath[0] log.info("Infrastructure area file is: %s" % infraPath) + else: + infraPath = "" cfgPath["infraPath"] = infraPath # read uMax Limit Raster - varUmaxDir = avalancheDir / "Inputs" / "UMAX" - varUmaxPath = sorted(list(varUmaxDir.glob("*.tif"))) - if len(varUmaxPath) == 0 or cfgFlowPy.getboolean("GENERAL", "variableUmaxLim") is False: - varUmaxPath = "" - elif len(varUmaxPath) > 1: - message = "More than one uMax Limit file .%s file in %s not allowed" % (varUmaxDir) - log.error(message) - raise AssertionError(message) - else: - varUmaxPath = varUmaxPath[0] + if cfgFlowPy.getboolean("GENERAL", "variableUmaxLim") is True: + varUmaxPath, available = gI.getAndCheckInputFiles(inputDir, "UMAX", "Umax", fileExt="raster") + if available == "No": + message = f"There is no variable UMAX file in supported format provided in {avalancheDir}/UMAX" + log.error(message) + raise AssertionError(message) log.info("uMax Limit file is: %s" % varUmaxPath) + else: + varUmaxPath = "" cfgPath["varUmaxPath"] = varUmaxPath # read variable Alpha Angle Raster - varAlphaDir = avalancheDir / "Inputs" / "ALPHA" - varAlphaPath = sorted(list(varAlphaDir.glob("*.tif"))) - if len(varAlphaPath) == 0 or cfgFlowPy.getboolean("GENERAL", "variableAlpha") is False: - varAlphaPath = "" - elif len(varAlphaPath) > 1: - message = "More than one variable alpha file .%s file in %s not allowed" % (varAlphaDir) - log.error(message) - raise AssertionError(message) - else: - varAlphaPath = varAlphaPath[0] + + if cfgFlowPy.getboolean("GENERAL", "variableAlpha") is True: + varAlphaPath, available = gI.getAndCheckInputFiles(inputDir, "ALPHA", "alpha", fileExt="raster") + if available == "No": + message = f"There is no variable ALPHA file in supported format provided in {avalancheDir}/ALPHA" + log.error(message) + raise AssertionError(message) log.info("variable Alpha file is: %s" % varAlphaPath) + else: + varAlphaPath = "" cfgPath["varAlphaPath"] = varAlphaPath # read variable Exponent Raster - varExponentDir = avalancheDir / "Inputs" / "EXP" - varExponentPath = sorted(list(varExponentDir.glob("*.tif"))) - if len(varExponentPath) == 0 or cfgFlowPy.getboolean("GENERAL", "variableExponent") is False: - varExponentPath = "" - elif len(varExponentPath) > 1: - message = "More than one variable exponent file .%s file in %s not allowed" % (varExponentDir) - log.error(message) - raise AssertionError(message) - else: - varExponentPath = varExponentPath[0] + + if cfgFlowPy.getboolean("GENERAL", "variableExponent") is True: + varExponentPath, available = gI.getAndCheckInputFiles(inputDir, "EXP", "exp", fileExt="raster") + if available == "No": + message = f"There is no variable EXPONENT file in supported format provided in {avalancheDir}/EXP" + log.error(message) + raise AssertionError(message) log.info("variable Exponent file is: %s" % varExponentPath) + else: + varExponentPath = "" cfgPath["varExponentPath"] = varExponentPath # check if forest should be used (assumed to be in the RES - 'RESISTANCE' directory) + # TODO: should we also allow forest as shp - file? - forestDir = avalancheDir / "Inputs" / "RES" # use directory for Resistance Layer for the forest layer - - if cfgFlowPy.getboolean("GENERAL", "forest") is False: - forestPath = "" - else: - patterns = ("*.tif", "*.asc", "*.TIF", "*.tiff", "*.TIFF", "*.ASC") - forestPath = [f for f in forestDir.iterdir() if any(f.match(p) for p in patterns)] - - if len(forestPath) == 0: - message = ( - "Please provide a Forest file in %s or set 'forest = False' in the .ini file" % forestDir - ) + if cfgFlowPy.getboolean("GENERAL", "forest") is True: + forestPath, available = gI.getAndCheckInputFiles(inputDir, "RES", "forest", fileExt="raster") + if available == "No": + message = f"There is no forest file in supported format provided in {avalancheDir}/RES" log.error(message) raise AssertionError(message) - elif len(forestPath) > 1: - message = ( - "Please provide exactly one Forest file in %s or set 'forest = False' in the .ini file" - % forestDir - ) - log.error(message) - raise AssertionError(message) - else: - forestPath = forestPath[0] - log.info("Forest File file is: %s" % forestPath) - + log.info("Forest file is: %s" % forestPath) + else: + forestPath = "" cfgPath["forestPath"] = forestPath # read DEM - demDir = avalancheDir / "Inputs" - patterns = ("*.tif", "*.asc", "*.TIF", "*.tiff", "*.TIFF", "*.ASC") - - _demPath = [f for f in demDir.iterdir() if any(f.match(p) for p in patterns)] - - if len(_demPath) == 0: - message = ( - "Please provide a DEM file in %s" % demDir - ) - log.error(message) - raise AssertionError(message) - elif len(_demPath) > 1: - message = ( - "Please provide exactly 1 (One!) DEM file in %s" % demDir - ) - log.error(message) - raise AssertionError(message) - else: - if os.path.splitext(_demPath[0])[1] in [".asc", ".ASC"]: - demPath = gI.getDEMPath(avalancheDir) - else: - demPath = _demPath[0] + demPath = gI.getDEMPath(avalancheDir) log.info("DEM file is: %s" % demPath) cfgPath["demPath"] = demPath diff --git a/docs/moduleCom4FlowPy.rst b/docs/moduleCom4FlowPy.rst index a5497d214..243d4db49 100644 --- a/docs/moduleCom4FlowPy.rst +++ b/docs/moduleCom4FlowPy.rst @@ -30,29 +30,10 @@ The motivational background and concepts behind the model, as well as a list of :ref:`theoryCom4FlowPy:com4FlowPy theory`. -Additional installation requirements and running the code +Running the code ---------------- -In order to run :py:mod:`com4FlowPy` you need to separately install the ``rasterio`` package, which is not included in the general AvaFrame requirements. -On Linux distributions inside the ``avaframe_env`` `conda environment -`_ it should usually work with: - - :: - - pip install rasterio - - -Running ``pip install rasterio`` inside the newly created ``avaframe_env`` `conda environment -`_ -is now (status: December 2024) also the recommended way to install ``rasterio`` on Windows. - -If this does not work, alternatively you can try installing ``rasterio`` via ``conda`` from the ``conda-forge`` channel. - -Additional (albeit not straight-forward) information can be found on: - - https://rasterio.readthedocs.io/en/stable/installation.html - - https://github.com/conda-forge/rasterio-feedstock - -Once the rasterio package is installed you can test :py:mod:`com4FlowPy` by: +You can test :py:mod:`com4FlowPy` by: 1. setting ``avalancheDir = data/avaFlowPy`` inside ``local_avaframeCfg.ini`` in the ``avaframe`` directory (if you don't already hava a ``local_avaFrameCfg.ini`` create it by copying and renaming the ``avaframeCfg.ini`` -- also see :ref:`moduleCom4FlowPy:Configuration`.) 2. running within the ``avaframe`` directory: @@ -270,4 +251,4 @@ If ``forestInteraction = True`` this layer will be written automatically (no nee .. - ``max_z``: :math:`z_{\delta}^{max}\;[\rm{m}]` might be defined based on observed max. velocities for different GMFs. .. b) forest module - .. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \ No newline at end of file + .. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ From ee4784d1cfbcdbc94bd774b4d82e2d1c0707fa67 Mon Sep 17 00:00:00 2001 From: paulasp Date: Mon, 24 Feb 2025 09:35:54 +0100 Subject: [PATCH 2/2] outputfile format set in ini - file --- avaframe/com4FlowPy/com4FlowPy.py | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py index a9f909400..2a9f08262 100755 --- a/avaframe/com4FlowPy/com4FlowPy.py +++ b/avaframe/com4FlowPy/com4FlowPy.py @@ -534,35 +534,40 @@ def mergeAndWriteResults(modelPaths, modelOptions): _ts = modelPaths["timeString"] demHeader = IOf.readRasterHeader(modelPaths["demPath"]) - # TODO: should we keep the modelPaths["outputFileFormat"] or can we use the DEM - format for the outputs + outputHeader = demHeader.copy() + if _oF == ".asc": + outputHeader["driver"] = "AAIGrid" + elif _oF == ".tif": + outputHeader["driver"] = "GTiff" + # TODO: do we want -99909 as nodata values or the DEM - nodata if 'flux' in _outputs: - output = IOf.writeResultToRaster(demHeader, flux, + output = IOf.writeResultToRaster(outputHeader, flux, modelPaths["resDir"] / "com4_{}_{}_flux".format(_uid, _ts), flip=True) if 'zDelta' in _outputs: - output = IOf.writeResultToRaster(demHeader, zDelta, + output = IOf.writeResultToRaster(outputHeader, zDelta, modelPaths["resDir"] / "com4_{}_{}_zdelta".format(_uid, _ts), flip=True) if 'cellCounts' in _outputs: - output = IOf.writeResultToRaster(demHeader, cellCounts, + output = IOf.writeResultToRaster(outputHeader, cellCounts, modelPaths["resDir"] / "com4_{}_{}_cellCounts".format(_uid, _ts), flip=True) if 'zDeltaSum' in _outputs: - output = IOf.writeResultToRaster(demHeader, zDeltaSum, + output = IOf.writeResultToRaster(outputHeader, zDeltaSum, modelPaths["resDir"] / "com4_{}_{}_zDeltaSum".format(_uid, _ts), flip=True) if 'routFluxSum' in _outputs: - output = IOf.writeResultToRaster(demHeader, routFluxSum, + output = IOf.writeResultToRaster(outputHeader, routFluxSum, modelPaths["resDir"] / "com4_{}_{}_routFluxSum".format(_uid, _ts), flip=True) if 'depFluxSum' in _outputs: - output = IOf.writeResultToRaster(demHeader, depFluxSum, + output = IOf.writeResultToRaster(outputHeader, depFluxSum, modelPaths["resDir"] / "com4_{}_{}_depFluxSum".format(_uid, _ts), flip=True) if 'fpTravelAngle' in _outputs: - output = IOf.writeResultToRaster(demHeader, fpTa, + output = IOf.writeResultToRaster(outputHeader, fpTa, modelPaths["resDir"] / "com4_{}_{}_fpTravelAngle".format(_uid, _ts), flip=True) if 'slTravelAngle' in _outputs: - output = IOf.writeResultToRaster(demHeader, slTa, + output = IOf.writeResultToRaster(outputHeader, slTa, modelPaths["resDir"] / "com4_{}_{}_slTravelAngle".format(_uid, _ts), flip=True) if 'travelLength' in _outputs: - output = IOf.writeResultToRaster(demHeader, travelLength, + output = IOf.writeResultToRaster(outputHeader, travelLength, modelPaths["resDir"] / "com4_{}_{}_travelLength".format(_uid, _ts), flip=True) # NOTE: @@ -570,10 +575,10 @@ def mergeAndWriteResults(modelPaths, modelOptions): # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("cell_counts%s" %(output_format)),cell_counts) # io.output_raster(modelPaths["demPath"], modelPaths["resDir"] / ("z_delta_sum%s" %(output_format)),z_delta_sum) if modelOptions["infraBool"]: # if infra - output = IOf.writeResultToRaster(demHeader, backcalc, + output = IOf.writeResultToRaster(outputHeader, backcalc, modelPaths["resDir"] / "com4_{}_{}_backcalculation".format(_uid, _ts), flip=True) if modelOptions["forestInteraction"]: - output = IOf.writeResultToRaster(demHeader, forestInteraction, + output = IOf.writeResultToRaster(outputHeader, forestInteraction, modelPaths["resDir"] / "com4_{}_{}_forestInteraction".format(_uid, _ts), flip=True) del output