From f37be35825e10383897269dcf05244a339346ba0 Mon Sep 17 00:00:00 2001 From: paulasp Date: Thu, 20 Nov 2025 18:16:51 +0100 Subject: [PATCH 1/2] correct documetation todo comment update input data type com1DFA add comment about nodata values add either or tests/test_deriveParameterSet.py change layout add comments add more info --- avaframe/com1DFA/com1DFACfg.ini | 2 +- avaframe/com1DFA/deriveParameterSet.py | 41 ++++-- avaframe/tests/test_deriveParameterSet.py | 98 +++++++++++++- docs/moduleCom1DFA.rst | 153 ++++++++++++---------- 4 files changed, 211 insertions(+), 83 deletions(-) diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index 121e0ef4d..e2b3cd77d 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/avaframe/com1DFA/com1DFACfg.ini @@ -238,7 +238,7 @@ meshCellSize = 5 meshCellSizeThreshold = 0.001 # clean DEMremeshed directory to ensure remeshing if chosen meshCellsize is different from rasters in Inputs/ cleanRemeshedRasters = True -# resize files read from Inputs/RASTERS to be resized to extent of DEM as resizeThreshold x meshCellSize +# remesh if extent matches computational DEM within resizeThreshold x meshCellSize resizeThreshold = 3 # Normal computation on rectangular grid diff --git a/avaframe/com1DFA/deriveParameterSet.py b/avaframe/com1DFA/deriveParameterSet.py index 009591e65..c463fa5f3 100644 --- a/avaframe/com1DFA/deriveParameterSet.py +++ b/avaframe/com1DFA/deriveParameterSet.py @@ -948,15 +948,8 @@ def checkExtentAndCellSize(cfg, inputFile, dem, fileType): cellSizeOld = inputField["header"]["cellsize"] demHeader = dem["header"] - # check if negative values or nan values - if np.any(np.isnan(inputField["rasterData"])): - message = "In %s file (%s) nan values found - this is not allowed" % ( - fileType, - inputFile.name, - ) - log.error(message) - raise AssertionError(message) - elif np.any(inputField["rasterData"] < 0): + # check if negative values in raster file + if np.any(inputField["rasterData"] < 0): message = "In %s file (%s) negative values found - this is not allowed" % ( fileType, inputFile.name, @@ -1025,6 +1018,36 @@ def checkExtentAndCellSize(cfg, inputFile, dem, fileType): returnStr = str(pathlib.Path("remeshedRasters", outFile.name)) remeshedFlag = "Yes" + # check if no data values only where DEM also has no data values + # if remeshed - potentially nans at edges, if inputField has a different origin/extent + # for example if extent of DEM is larger -> nans in this region in remeshed input file + # first maks nans values of DEM as there, also inputField is allowed to have nans + nanDEMMasked = np.where(np.isnan(dem["rasterData"]), -9999, inputField["rasterData"]) + # search for indices where nans come from remeshing use difference in extent prior to remeshing + # if diff is negative on Left side DEM is smaller + # if diff is negative on right side DEM is larger + nNeglectRowsMin = int(np.ceil(abs(min(0, diffY0 / dem["header"]["cellsize"])))) + nNeglectRowsMax = int(np.ceil(abs(max(0, diffY1 / dem["header"]["cellsize"])))) + # if diff is negative on lower side DEM is smaller + # if diff is negative on right side DEM is larger + nNeglectColsMin = int(np.ceil(abs(min(0, diffX0 / dem["header"]["cellsize"])))) + nNeglectColsMax = int(np.ceil(abs(max(0, diffX1 / dem["header"]["cellsize"])))) + maxRows = dem["header"]["nrows"] - nNeglectRowsMin + maxCols = dem["header"]["ncols"] - nNeglectColsMin + + # add mask where nans come from remeshing + # change order because diff is with origin lower!! + nanDEMMaskedLimited = nanDEMMasked[nNeglectRowsMax:maxRows, nNeglectColsMax:maxCols] + + # check if nan values in raster file data (excluding nans from remeshing and where also DEM has nans) + if np.any(np.isnan(nanDEMMaskedLimited)): + message = "In %s file (%s) nan values found inside DEM extent - this is not allowed" % ( + fileType, + inputFile.name, + ) + log.error(message) + raise AssertionError(message) + return returnStr, outFile, remeshedFlag diff --git a/avaframe/tests/test_deriveParameterSet.py b/avaframe/tests/test_deriveParameterSet.py index 03f89f2db..289819d38 100644 --- a/avaframe/tests/test_deriveParameterSet.py +++ b/avaframe/tests/test_deriveParameterSet.py @@ -686,7 +686,7 @@ def test_checkExtentAndCellSize(tmp_path): assert newRaster2["header"]["xllcenter"] == 1.0 assert newRaster2["header"]["yllcenter"] == 5.0 - inputFile2 = inDirR / "inputFile1.asc" + inputFile2 = inDirR / "inputFile2.asc" headerInput2 = { "nrows": 5, "ncols": 5, @@ -706,6 +706,102 @@ def test_checkExtentAndCellSize(tmp_path): assert dP.checkExtentAndCellSize(cfg, inputFile2, dem, "mu") assert "Lower left center coordinates of DEM and " in str(e.value) + demField = np.ones((4, 5)) + dem = { + "header": { + "nrows": 4, + "ncols": 5, + "xllcenter": 1, + "yllcenter": 5, + "cellsize": 1, + "nodata_value": -9999, + "driver": "AAIGrid", + }, + "rasterData": demField, + } + + demField[0, :] = np.nan + dem["header"]["transform"] = IOf.transformFromASCHeader(dem["header"]) + dem["header"]["crs"] = rasterio.crs.CRS() + + inputFile = inDirR / "inputFile3.asc" + headerInput = { + "nrows": 4, + "ncols": 5, + "xllcenter": 1, + "yllcenter": 5, + "cellsize": 1, + "nodata_value": -9999, + "driver": "AAIGrid", + } + headerInput["transform"] = IOf.transformFromASCHeader(headerInput) + headerInput["crs"] = rasterio.crs.CRS() + + inField = np.ones((4, 5)) + inField[2, 2] = 10.0 + inField[0, :] = np.nan + IOf.writeResultToRaster(headerInput, inField, inputFile.parent / inputFile.stem, flip=True) + + testFile4, outFile4, remeshedFlag4 = dP.checkExtentAndCellSize(cfg, inputFile, dem, "rel") + + assert remeshedFlag4 == "No" + + inputFile = inDirR / "inputFile4.asc" + inField = np.ones((4, 5)) + inField[2, 2] = 10.0 + inField[0, :] = np.nan + inField[1, 1] = np.nan + IOf.writeResultToRaster(headerInput, inField, inputFile.parent / inputFile.stem, flip=True) + + with pytest.raises(AssertionError) as e: + assert dP.checkExtentAndCellSize(cfg, inputFile, dem, "rel") + assert "nan values found inside DEM extent - this is not allowed" in str(e.value) + + inputFile = inDirR / "inputFile5.asc" + headerInput = { + "nrows": 4, + "ncols": 5, + "xllcenter": 1.3, + "yllcenter": 4.2, + "cellsize": 1, + "nodata_value": -9999, + "driver": "AAIGrid", + } + headerInput["transform"] = IOf.transformFromASCHeader(headerInput) + headerInput["crs"] = rasterio.crs.CRS() + + inField = np.ones((4, 5)) + inField[2, 2] = 10.0 + inField[0, :] = np.nan + IOf.writeResultToRaster(headerInput, inField, inputFile.parent / inputFile.stem, flip=True) + + testFile5, outFile5, remeshedFlag5 = dP.checkExtentAndCellSize(cfg, inputFile, dem, "rel") + newRaster5 = IOf.readRaster((inDir / testFile5)) + assert remeshedFlag5 == "Yes" + assert np.isnan(newRaster5["rasterData"][0, :]).all() + assert np.isnan(newRaster5["rasterData"][-1, :]).all() + + inputFile = inDirR / "inputFile51.asc" + headerInput = { + "nrows": 4, + "ncols": 5, + "xllcenter": 1.3, + "yllcenter": 4.2, + "cellsize": 1, + "nodata_value": -9999, + "driver": "AAIGrid", + } + headerInput["transform"] = IOf.transformFromASCHeader(headerInput) + headerInput["crs"] = rasterio.crs.CRS() + + inField = np.ones((4, 5)) + inField[2, 2] = 10.0 + inField[0:2, :] = np.nan + IOf.writeResultToRaster(headerInput, inField, inputFile.parent / inputFile.stem, flip=True) + + with pytest.raises(AssertionError) as e: + assert dP.checkExtentAndCellSize(cfg, inputFile, dem, "rel") + assert "nan values found inside DEM extent - this is not allowed" in str(e.value) # Produced by AI (test): diff --git a/docs/moduleCom1DFA.rst b/docs/moduleCom1DFA.rst index 101675e58..c5f0222e6 100644 --- a/docs/moduleCom1DFA.rst +++ b/docs/moduleCom1DFA.rst @@ -48,47 +48,58 @@ folder structure described below. In the directory ``Inputs``, the following files are required. Be aware that ALL inputs have to be provided in the same -projection: +projection. Allowed file types are either raster files, i.e. `ESRI grid format `_ +or GeoTIFF format, or shape files and specified for the respective input type below: -* digital elevation model as raster file with either `ESRI grid format `_ - or GeoTIFF format. The format of the DEM determines the format of the output files. -* release area scenario as (multi-) polygon shapefile (in Inputs/REL; multiple features are possible) - - - the release area polygon must not contain any "holes" or inner rings - - the release area name should not contain an underscore, if so '_AF' is added. - - recommended attributes are *name*, *thickness* (see :ref:`moduleCom1DFA:Release-, entrainment thickness settings`) - and *ci95* (see :ref:`moduleAna4Stats:probAna - Probability maps`) - - ALL features within one shapefile are released at the same time (and interact), this is what we refer to as *scenario* - - if you want to simulate different scenarios with the same features, you have to copy them to separate shapefiles +* **digital elevation model as raster file. The format of the DEM determines the format of the output files.** +* release area scenario as (multi-) polygon shapefile OR raster file (in Inputs/REL; only shapefiles OR raster files) + - either polygon shapefile(s): + - the release area polygon must not contain any "holes" or inner rings + - multiple features are allowed + - recommended attributes are *name*, *thickness* (see :ref:`moduleCom1DFA:Release-, entrainment thickness settings`) and *ci95* (see :ref:`moduleAna4Stats:probAna - Probability maps`) + - or raster file(s): + - cells with non-zero values define the release area + - cell value can be read as thickness (measured normal to the slope) (see :ref:`moduleCom1DFA:Release-, entrainment thickness settings`) + - negative values are not allowed (specified no-data values are not considered) + - ALL features within one shapefile or all cells with non-zero values in a raster file are released at the same time (and interact), this is what we refer to as *scenario* + - if you want to simulate different scenarios with the same features, you have to copy them to separate shapefiles/raster files + - the release area name should not contain an underscore, if so '_AF' is added. and the following files are optional. Please note: in the standard configuration (i.e. ``simTypeList = available``) , the *null* variant is always run! I.e. if a resistance and/or an entrainment file is given (as described below), at least two results are generated: the *null* variant and the variant with entrainment and/or resistance. -* one entrainment area (multi-) polygon shapefile (in Inputs/ENT) - - - marks the (multiple) areas where entrainment can occur. - - attribute *thickness* (see :ref:`moduleCom1DFA:Release-, entrainment thickness settings`) - - must not contain any "holes" or inner rings - - -* one resistance area (multi-) polygon shapefile (in Inputs/RES) - - - marks the (multiple) areas where resistance is considered - - resistance areas must not contain any "holes" or inner rings - - please consider the information about resistance below :ref:`moduleCom1DFA:Resistance setup` - - -* one secondary release area (multi-) polygon shapefile (in Inputs/SECREL) +* **one entrainment area as (multi-) polygon shapefile OR raster file (in Inputs/ENT)** + - marks the (multiple) areas where entrainment can occur. + - either polygon shapefile: + - attribute *thickness* (see :ref:`moduleCom1DFA:Release-, entrainment thickness settings`) + - must not contain any "holes" or inner rings + - or raster file: + - cells with non-zero values define where entrainment can occur + - cell value can be read as thickness (measured normal to the slope) (see :ref:`moduleCom1DFA:Release-, entrainment thickness settings`) + - negative values are not allowed (specified no-data values are not considered) + +* **one resistance area as (multi-) polygon shapefile OR raster file (in Inputs/RES)** + - marks the (multiple) areas where resistance is considered + - please consider the information about resistance below :ref:`moduleCom1DFA:Resistance setup` + - either polygon shapefile: + - resistance areas must not contain any "holes" or inner rings + - or raster file: + - cells with non-zero values define where entrainment can occur + - negative values are not allowed (specified no-data values are not considered) + - if the cellsize does not match the requested meshCellSize, the file is + remeshed if within `resizeThreshold` + +* **one secondary release area (multi-) as polygon shapefile OR raster file (in Inputs/SECREL)** - can have multiple release areas, each as one feature - same setup as the release area scenario (see above) - features will release as soon as at least one particle enters its area - release area polygons must not contain any "holes" or inner rings -* raster files for the Voellmy friction parameters :math:`\mu` and :math:`\xi` (in Inputs/RASTERS) +* **raster files for the Voellmy friction parameters :math:`\mu` and :math:`\xi` (in Inputs/RASTERS)** - spatial field of :math:`\mu` and :math:`\xi` values with same extent as DEM - file names need to end with ``_mu.*`` and ``_xi.*`` @@ -96,7 +107,7 @@ at least two results are generated: the *null* variant and the variant with entr - if ``meshCellSize`` is different from simulation ``meshCellSize`` fields will be remeshed - only used if ``frictionModel`` is set to ``spatialVoellmy`` -* one ``_cropshape.shp`` shape file (in Inputs/POLYGONS) +* **one ``_cropshape.shp`` shape file (in Inputs/POLYGONS)** - provides a polygon located inside the DEM to define area for report plots of peak fields (bounds of polygon) - if not provided peak fields are shown for the extent where peak field values are nonzero @@ -110,53 +121,51 @@ Release-, entrainment thickness settings .. Note:: Thickness is unambiguous: it is measured normal to the slope. -Release, entrainment and secondary release thickness can be specified in two different ways: - -1. Via **shape file**: - - - add an attribute called `thickness` for each feature - - important: ALL features have to have a single thickness value, which can differ between features - - for entrainment area only: if the thickness value is missing, the thickness value is - taken from `entThIfMissingInShp` (default 0.3 m) in the configuration file. If multiple features are - in the entrainment file the thickness attribute has to be set either for ALL or NONE of the features. - - for backwards compatibility, the attribute 'd0' also works, but we suggest to use `thickness` in new projects - - set the flag `THICKNESSFromShp` (i.e. relThFromFile, entThFromFile, - secondaryRelthFromShp) to True in the configuration file (default is True) - - a parameter variation can be added with the `THICKNESSPercentVariation` - parameter in the configuration file in the form of - ``+-percentage$numberOfSteps``. Provided a `+` a positive variation will be - performed, if `-` is given, only a negative variation is performed. If no - sign is given: both directions will be used. Additionally, a variation can be - added with the `THICKNESSRangeVariation` parameter in the configuration file - in the form of ``+-range$numberOfSteps``. Provided a `+` a positive variation - will be performed, if `-` is given, only a negative variation is performed. - If no sign is given: both directions will be used. Furthermore, there is the - option to vary the thickness in a range of +- the 95% confidence interval - value, which is also read from the shape file (requires an attribute called - ci95). In order to use this variation, set the 'THICKESSRangeFromCiVariation' - to ``ci95$numberOfSteps``. +Release, entrainment and secondary release thickness can be specified in two different ways, 1) directly from the provided +input file (shape file or raster file) or 2) through the :py:mod:`com1DFA` configuration file: + +1. **Read from input file**: + +- Via **shape file**: + + - add an attribute called `thickness` for each feature + - important: ALL features have to have a single thickness value, which can differ between features + - for entrainment area only: if the thickness value is missing, the thickness value is + taken from `entThIfMissingInShp` (default 0.3 m) in the configuration file. If multiple features are + in the entrainment file the thickness attribute has to be set either for ALL or NONE of the features. + - for backwards compatibility, the attribute 'd0' also works, but we suggest to use `thickness` in new projects + - set the flag `THICKNESSFromFile` (i.e. relThFromFile, entThFromFile, + secondaryRelThFromFile) to True in the configuration file (default is True) + - a parameter variation can be added with the `THICKNESSPercentVariation` + parameter in the configuration file in the form of + ``+-percentage$numberOfSteps``. Provided a `+` a positive variation will be + performed, if `-` is given, only a negative variation is performed. If no + sign is given: both directions will be used. Additionally, a variation can be + added with the `THICKNESSRangeVariation` parameter in the configuration file + in the form of ``+-range$numberOfSteps``. Provided a `+` a positive variation + will be performed, if `-` is given, only a negative variation is performed. + If no sign is given: both directions will be used. Furthermore, there is the + option to vary the thickness in a range of +- the 95% confidence interval + value, which is also read from the shape file (requires an attribute called + ci95). In order to use this variation, set the 'THICKESSRangeFromCiVariation' + to ``ci95$numberOfSteps``. + +- Via **raster file**: + + - non-zero cell values are read as thickness + - set the flag `THICKNESSFromFile` (i.e. relThFromFile, entThFromFile, + secondaryRelThFromFile) to True in the configuration file (default is True) + - no thickness variation is available if input type is raster file and `THICKNESSFromFile` is *True* 2. Via **configuration file (ini)**: - - set the flag 'THICKNESSFromShp' to False - - provide your desired thickness value in the respective THICKNESS parameter (i.e. relTh, entTh or secondaryRelth) - - in addition to the `THICKNESSPercentVariation` and `THICKNESSRangeVariation` - options (see option 1) and the standard variation options in - :ref:`configuration:Configuration`, you can also directly set e.g. `relTh = - 1.$50$2`, ``referenceValue$+-percentage$numberOfSteps``, resulting in a - variation of relTh from 0.5 to 1.5m in two steps. - -Only available for release thickness: - -3. Via **release thickness file**: - - - set the flag 'relThFromFile' to False - - set the flag 'relThFromFile' to True - - save a raster file with info on release thickness as raster file in - ``Inputs/RELTH`` the number of rows and columns must match the DEM raster - with desired meshCellSize (recommended) - - if the cellsize does not match the requested meshCellSize, the file is - remeshed. + - set the flag 'THICKNESSFromFile' to False + - provide your desired thickness value in the respective THICKNESS parameter (i.e. relTh, entTh or secondaryRelth) + - in addition to the `THICKNESSPercentVariation` and `THICKNESSRangeVariation` + options (see option 1) and the standard variation options in + :ref:`configuration:Configuration`, you can also directly set e.g. `relTh = + 1.$50$2`, ``referenceValue$+-percentage$numberOfSteps``, resulting in a + variation of relTh from 0.5 to 1.5m in two steps. Friction parameters From eb9f7b4affe279161772628355d5ebcb925be7d6 Mon Sep 17 00:00:00 2001 From: Felix Oesterle Date: Tue, 9 Dec 2025 15:20:38 +0100 Subject: [PATCH 2/2] Pin cython due to problems with numpy 2.0 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index e5780533a..0791b0b29 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,7 +32,7 @@ dependencies = [ "scipy", "cmcrameri", "seaborn", - "cython", + "cython<3.0", # cython 3.x compiles code that needs numpy > 2.0 "pandas", "shapely", "configUpdater",