Add mask fill operator so that bad fill values are set to NaN

src/CSET/operators/misc.py

@@ -615,3 +615,73 @@ def differentiate(
         return new_cubelist[0]
     else:
         return new_cubelist
+
+
+def _mask_fill_cube(cube: iris.cube.Cube, ulp_factor=10):
+    """
+    Avoid plotting data flagged as bad/missing.
+
+    Force masked data and known fill values to np.nan
+    so they are not plotted.
+
+    """
+    import dask.array as da
+
+    x = cube.lazy_data()
+    fill_values = []
+    # NetCDF-style fill value (if present)
+    try:
+        fv = getattr(x._meta, "fill_value", None)
+        if fv is not None:
+            fill_values.append(fv)
+    except AttributeError:
+        pass  # x has no _meta (plain ndarray)
+
+    # Known Cardington fill values
+    fill_values.extend([1e10, 1e11, 999999])
+
+    if np.ma.isMaskedArray(x):
+        x_data = np.ma.getdata(x)
+        x_mask = np.ma.getmaskarray(x)
+    else:
+        x_data = x
+        x_mask = None
+
+    data = da.asarray(x_data, dtype=np.float32)
+
+    if x_mask is not None:
+        m0 = da.asarray(x_mask, dtype=bool)
+        # Convert masked elements into NaN immediately
+        data = da.where(m0, np.nan, data)
+    else:
+        m0 = da.zeros(data.shape, dtype=bool, chunks=data.chunks)
+
+    # Build mask
+    m_fill = da.zeros(data.shape, dtype=bool, chunks=data.chunks)
+    for fv in fill_values:
+        ulp = ulp_factor * np.spacing(np.float32(fv))
+        m_fill |= da.isclose(data, np.float32(fv), rtol=0, atol=ulp)
+
+    if not da.any(m0 | m_fill).compute():
+        return cube  # nothing to clean
+
+    masked = da.ma.masked_array(data, mask=(m0 | m_fill))
+    y = da.ma.filled(masked, np.nan)
+
+    return cube.copy(data=y)
+
+
+def mask_fill_values(cubes, ulp_factor=10):
+    """
+    Apply _mask_fill_value to every cube in the CubeList.
+
+    This must be run AFTER any operator that recreates data.
+    """
+    if not isinstance(cubes, CubeList):
+        cubes = CubeList([cubes])
+
+    cleaned = CubeList()
+    for cube in cubes:
+        cleaned.append(_mask_fill_cube(cube, ulp_factor=ulp_factor))
+
+    return cleaned

tests/operators/test_misc.py

@@ -541,3 +541,80 @@ def test_extract_common_points_nocommonpoints(vertical_profile_cube):
         misc.extract_common_points(
             cubes=iris.cube.CubeList([cube1, cube2]), coordinate="pressure"
         )
+
+
+def _make_cube(data, units="1"):
+    """Tiny cube generator."""
+    data = np.asanyarray(data)
+    if data.ndim != 2:
+        raise ValueError(f"Expected 2D data, got shape {data.shape}")
+
+    ny, nx = data.shape
+    lat = iris.coords.DimCoord(np.arange(ny), standard_name="latitude", units="degrees")
+    lon = iris.coords.DimCoord(
+        np.arange(nx), standard_name="longitude", units="degrees"
+    )
+    return iris.cube.Cube(
+        data,
+        dim_coords_and_dims=[(lat, 0), (lon, 1)],
+        units=units,
+    )
+
+
+def test_mask_fill_value_no_change():
+    """Test mask fill value with no change."""
+    cube = _make_cube([[1.0, 2.0]])
+    out = misc._mask_fill_cube(cube)
+    # return same object
+    assert out is cube
+
+
+def test_mask_fill_value_sentinel():
+    """Test mask with known fill value."""
+    cube = _make_cube([[1.0, 1e10, 3.0]])
+    out = misc._mask_fill_cube(cube)
+    data = out.data.compute() if hasattr(out.data, "compute") else out.data
+    assert np.isnan(data[0, 1])
+    assert np.allclose(data[0, [0, 2]], [1.0, 3.0])
+
+
+def test_mask_fill_value_masked_array():
+    """Test with masked array input."""
+    data = np.ma.array([[1.0, 2.0]], mask=[[False, True]])
+    cube = _make_cube(data)
+    out = misc._mask_fill_cube(cube)
+    result = out.data.compute() if hasattr(out.data, "compute") else out.data
+    assert not np.isnan(result[0, 0])
+    result = result.filled(np.nan) if np.ma.isMaskedArray(result) else result
+    assert np.isnan(result[0, 1])
+
+
+def test_mask_fill_value_ulp():
+    """Test with ulp_factor input."""
+    fv = np.float32(1e10)
+    near_fv = fv + np.spacing(fv) * 5  # within tolerance
+    cube = _make_cube([[near_fv]])
+    out = misc._mask_fill_cube(cube, ulp_factor=10)
+    data = out.data.compute() if hasattr(out.data, "compute") else out.data
+    assert np.isnan(data[0, 0])
+
+
+def test_mask_fill_value_combined():
+    """Test with combined mask and real input."""
+    data = np.ma.array([[1e10, 2.0]], mask=[[False, True]])
+    cube = _make_cube(data)
+    out = misc._mask_fill_cube(cube)
+    result = out.data.compute() if hasattr(out.data, "compute") else out.data
+    assert np.isnan(result[0, 0])  # sentinel
+    assert np.isnan(result[0, 1])  # masked
+
+
+def test_mask_fill_values_cubelist():
+    """Test with cubelist input."""
+    cubes = iris.cube.CubeList([_make_cube([[1e10]]), _make_cube([[2.0]])])
+    out = misc.mask_fill_values(cubes)
+    assert isinstance(out, iris.cube.CubeList)
+    data0 = out[0].data.compute() if hasattr(out[0].data, "compute") else out[0].data
+    data1 = out[1].data.compute() if hasattr(out[1].data, "compute") else out[1].data
+    assert np.isnan(data0[0, 0])
+    assert np.allclose(data1[0, 0], 2.0)