Fix average filtering rasters with nodata pixels

ocsmesh/raster.py

@@ -965,8 +965,13 @@ def average_filter(
                     mask_below = ma.getdata(outband) < drop_below
                     outband[mask_below] = np.nan
 
+                # Since the nbmean checks for isnan, other fill values
+                # will be treated as normal numbers
                 outband_new = generic_filter(
-                    outband, LowLevelCallable(nbmean.ctypes), size=size)
+                    outband.filled(np.nan),
+                    LowLevelCallable(nbmean.ctypes),
+                    size=size
+                )
 
                 # Mask raster based on result of filter?
                 if drop_above is not None:
@@ -978,6 +983,7 @@ def average_filter(
 
                 outband_new[mask] = dst.nodatavals[bnd_idx]
                 outband_ma = ma.masked_array(outband_new, mask=mask)
+                ma.set_fill_value(outband_ma, dst.nodatavals[bnd_idx])
                 dst.write_band(i, outband_ma)
 
 

tests/api/raster.py

@@ -0,0 +1,58 @@
+import shutil
+import tempfile
+import unittest
+from pathlib import Path
+
+import numpy as np
+
+import ocsmesh
+from ocsmesh.utils import raster_from_numpy
+
+
+class Raster(unittest.TestCase):
+    def setUp(self):
+        self.tdir = Path(tempfile.mkdtemp())
+        self.rast1 = self.tdir / 'rast_1.tif'
+        self.rast2 = self.tdir / 'rast_2.tif'
+        self.rast3 = self.tdir / 'rast_3.tif'
+
+        rast_xy = np.mgrid[-1:1:0.01, -0.7:0.7:0.01]
+        rast_z_1 = np.ones_like(rast_xy[0]) * 10
+        rast_z_2 = np.ma.MaskedArray(
+            np.ones_like(rast_xy[0]) * 10,
+            mask=np.random.random(size=rast_xy[0].shape) < 0.2,
+            fill_value=np.nan
+        )
+        rast_z_3 = np.ma.MaskedArray(
+            np.ones_like(rast_xy[0]) * 10,
+            mask=np.random.random(size=rast_xy[0].shape) < 0.2,
+            fill_value=-1e+15
+        )
+
+        raster_from_numpy(self.rast1, rast_z_1, rast_xy, 4326)
+        raster_from_numpy(self.rast2, rast_z_2, rast_xy, 4326)
+        raster_from_numpy(self.rast3, rast_z_3, rast_xy, 4326)
+
+
+    def tearDown(self):
+        shutil.rmtree(self.tdir)
+
+
+    def test_avg_filter_nomask(self):
+        rast = ocsmesh.Raster(self.rast1)
+        rast.average_filter(size=17)
+        self.assertTrue(np.all(rast.get_values() == 10))
+
+
+    def test_avg_filter_masked_nanfill(self):
+        rast = ocsmesh.Raster(self.rast2)
+        rast.average_filter(size=17)
+        self.assertTrue(
+            np.all(rast.values[~np.isnan(rast.values)] == 10))
+
+
+    def test_avg_filter_masked_nonnanfill(self):
+        rast = ocsmesh.Raster(self.rast3)
+        rast.average_filter(size=17)
+        self.assertTrue(
+            np.all(rast.values[rast.values != rast.nodata] == 10))