diff --git a/avaframe/com1DFA/deriveParameterSet.py b/avaframe/com1DFA/deriveParameterSet.py index 8499d6d1a..c80bcc8ca 100644 --- a/avaframe/com1DFA/deriveParameterSet.py +++ b/avaframe/com1DFA/deriveParameterSet.py @@ -426,11 +426,16 @@ def checkThicknessSettings(cfg, thName, inputSimFiles): ) log.error(message) raise AssertionError(message) - if (cfg["GENERAL"].getboolean("timeDependentRelease") and inputSimFiles["entResInfo"]["relTh" + "FileType"] in [".asc", ".tif"]): + if cfg["GENERAL"].getboolean("timeDependentRelease") and inputSimFiles["entResInfo"][ + "relTh" + "FileType" + ] in [".asc", ".tif"]: message = "When release is time dependent, release file needs to be provided as shape file" log.error(message) raise FileNotFoundError(message) - if (cfg["GENERAL"].getboolean("timeDependentRelease") and cfg["GENERAL"].getboolean("relThFromFile") is False): + if ( + cfg["GENERAL"].getboolean("timeDependentRelease") + and cfg["GENERAL"].getboolean("relThFromFile") is False + ): message = "When release is time dependent, relThFromFile needs to be set to True" log.error(message) raise FileNotFoundError(message) @@ -890,11 +895,8 @@ def appendThicknessToCfg(cfg): for count, id in enumerate(idList): thNameId = thType + id if thNameId in cfg["GENERAL"].keys(): - log.info( - "Thickness value for %s already set in initial config file, \ - read from there not from shp file" - % thNameId - ) + log.info("Thickness value for %s already set in initial config file, \ + read from there not from shp file" % thNameId) else: cfgGen[thNameId] = str(float(thicknessList[count])) @@ -968,6 +970,17 @@ def checkExtentAndCellSize(cfg, inputFile, dem, fileType): """ inputField = IOf.readRaster(inputFile) + + # Check for NaN values before remeshing - not allowed in non-DEM input rasters + if fileType.upper() != "DEM": + if np.isnan(inputField["rasterData"]).any(): + message = "In %s file (%s) nan values found - this is not allowed" % ( + fileType, + inputFile.name, + ) + log.error(message) + raise AssertionError(message) + cellSizeOld = inputField["header"]["cellsize"] demHeader = dem["header"] diff --git a/avaframe/in3Utils/geoTrans.py b/avaframe/in3Utils/geoTrans.py index eb415640e..fb3ca157a 100644 --- a/avaframe/in3Utils/geoTrans.py +++ b/avaframe/in3Utils/geoTrans.py @@ -201,7 +201,15 @@ def resizeData(raster, rasterRef): X, Y = np.meshgrid(xgrid, ygrid) Points = {"x": X, "y": Y} Points, _ = projectOnRaster(raster, Points, interp="bilinear") - raster["rasterData"] = Points["z"] + bilinearData = Points["z"] + + if np.isnan(bilinearData).any(): + PointsNearest = {"x": X, "y": Y} + PointsNearest, _ = projectOnRaster(raster, PointsNearest, interp="nearest") + nanMask = np.isnan(bilinearData) + bilinearData[nanMask] = PointsNearest["z"][nanMask] + + raster["rasterData"] = bilinearData return raster["rasterData"], rasterRef["rasterData"] diff --git a/avaframe/tests/test_deriveParameterSet.py b/avaframe/tests/test_deriveParameterSet.py index a6d6ab749..3c1315c88 100644 --- a/avaframe/tests/test_deriveParameterSet.py +++ b/avaframe/tests/test_deriveParameterSet.py @@ -740,15 +740,30 @@ def test_checkExtentAndCellSize(tmp_path): headerInput["transform"] = IOf.transformFromASCHeader(headerInput) headerInput["crs"] = rasterio.crs.CRS() + # Test that NaN values in non-DEM input files are not allowed + 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" + with pytest.raises(AssertionError) as e: + dP.checkExtentAndCellSize(cfg, inputFile, dem, "rel") + assert "nan values found - this is not allowed" in str(e.value) + # Test with interior NaN (not just on border) inputFile = inDirR / "inputFile4.asc" inField = np.ones((4, 5)) inField[2, 2] = 10.0 @@ -757,15 +772,16 @@ def test_checkExtentAndCellSize(tmp_path): 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) + dP.checkExtentAndCellSize(cfg, inputFile, dem, "rel") + assert "nan values found - this is not allowed" in str(e.value) + # Test that DEM files ARE allowed to have NaN values inputFile = inDirR / "inputFile5.asc" headerInput = { "nrows": 4, "ncols": 5, - "xllcenter": 1.3, - "yllcenter": 4.2, + "xllcenter": 1, + "yllcenter": 5, "cellsize": 1, "nodata_value": -9999, "driver": "AAIGrid", @@ -778,33 +794,10 @@ def test_checkExtentAndCellSize(tmp_path): 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) + # fileType="DEM" should allow NaNs + testFile5, outFile5, remeshedFlag5 = dP.checkExtentAndCellSize(cfg, inputFile, dem, "DEM") + assert remeshedFlag5 == "No" - 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/avaframe/tests/test_geoTrans.py b/avaframe/tests/test_geoTrans.py index 6f20facad..95e9b578b 100644 --- a/avaframe/tests/test_geoTrans.py +++ b/avaframe/tests/test_geoTrans.py @@ -110,6 +110,20 @@ def test_resizeData(): testRes = np.allclose(np.nan_to_num(data1 - data2), zSol, atol=tol) assert testRes + raster3 = { + "header": {"ncols": 3, "nrows": 3, "xllcenter": 0.0, "yllcenter": 0.0, "cellsize": 10.0}, + "rasterData": np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]), + } + raster4 = { + "header": {"ncols": 5, "nrows": 5, "xllcenter": 0.0, "yllcenter": 0.0, "cellsize": 5.0}, + "rasterData": np.zeros((5, 5)), + } + + data3, _ = geoTrans.resizeData(raster3, raster4) + + assert np.isnan(data3).sum() == 0 + assert np.isclose(data3[2, 2], 5.0) + # def test_findAngleProfile(capfd): def test_findAngleProfile(): diff --git a/docs/moduleCom1DFA.rst b/docs/moduleCom1DFA.rst index 509f16d47..81b12941c 100644 --- a/docs/moduleCom1DFA.rst +++ b/docs/moduleCom1DFA.rst @@ -63,7 +63,7 @@ or GeoTIFF format, or shape files and specified for the respective input type be - 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) + - negative values and no-data values are not allowed - 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. @@ -80,7 +80,7 @@ at least two results are generated: the *null* variant and the variant with entr - 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) + - negative values and no-data values are not allowed * **one resistance area as (multi-) polygon shapefile OR raster file (in Inputs/RES)** - marks the (multiple) areas where resistance is considered @@ -89,7 +89,7 @@ at least two results are generated: the *null* variant and the variant with entr - 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) + - negative values and no-data values are not allowed - if the cellsize does not match the requested meshCellSize, the file is remeshed if within `resizeThreshold`