[in3Utils] (geotrans) fill bilinear remesh edge NaNs with nearest fallback

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"]
 

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"]
 
 

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):
 

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():

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`