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

src/CSET/operators/misc.py

@@ -703,3 +703,133 @@ def differentiate(
         return new_cubelist[0]
     else:
         return new_cubelist
+
+
+def _mask_fill_cube(
+    cube: iris.cube.Cube,
+    ulp_factor: int = 10,
+) -> iris.cube.Cube:
+    """
+    Replace masked and fill-value data with NaNs.
+
+    Parameters
+    ----------
+    cube : iris.cube.Cube
+        Input cube to clean.
+
+    ulp_factor : int, optional
+        Number of floating-point ULPs (units in the last place) used when
+        comparing values against known fill values. Larger values are more
+        tolerant of floating-point rounding differences. Default is 10.
+
+    Returns
+    -------
+    iris.cube.Cube
+        Cube with masked and fill-value points replaced by ``NaN``.
+
+    Notes
+    -----
+     The returned cube preserves metadata and coordinates.
+     Data are processed lazily using Dask where possible.
+     If no masked or fill values are detected, the original cube is
+      returned unchanged.
+     Genuine NetCDF ``_FillValue`` handling is normally performed by Iris
+      during file loading, but this routine additionally converts masked
+      points to NaNs and handles known sentinel values.
+    """
+    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 fill values
+    # - 1e10 observed as NetCDF _FillValue in some archived variables.
+    # - 1e11 documented as the data value for missing/bad core data flags.
+    fill_values.extend([1e10, 1e11])
+
+    # Defensive fallback: other NumPy masked-array default fill values.
+    fill_values.extend([999999, -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 * abs(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: iris.cube.Cube | CubeList,
+    ulp_factor: int = 10,
+) -> CubeList:
+    """
+    Replace masked and fill-value data with NaNs in one or more cubes.
+
+    Applies :func:`_mask_fill_cube` to every cube in the supplied
+    CubeList. This is primarily intended for observational
+    datasets where a combination of NetCDF ``_FillValue`` masking and
+    dataset-specific sentinel values are used to represent missing or
+    unusable data.
+
+    Parameters
+    ----------
+    cubes : iris.cube.Cube or iris.cube.CubeList
+        Cube or CubeList to clean.
+
+    ulp_factor : int, optional
+        Number of floating-point ULPs used when comparing values against
+        known fill values. Passed directly to
+        :func:`_mask_fill_cube`. Default is 10.
+
+    Returns
+    -------
+    iris.cube.CubeList
+        CubeList containing cleaned cubes.
+
+    Notes
+    -----
+    This operator should generally be applied after reading observational
+    data and before plotting.
+
+    It may also be useful after operators that recreate or transform data,
+    where masked values could otherwise propagate into derived diagnostics
+    and appear as spurious spikes or extrema in plots.
+    """
+    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

@@ -567,39 +567,91 @@ def test_extract_common_points_nocommonpoints(vertical_profile_cube):
         )
 
 
-def test_remove_scalar_coord():
-    """Test that scalar coordinate be removed."""
-    # Create simple 1D cube
-    data = np.arange(5)
-    time = iris.coords.DimCoord(
-        np.arange(5), standard_name="time", units="hours since 1970-01-01"
+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"
     )
-    cube = iris.cube.Cube(data, dim_coords_and_dims=[(time, 0)])
-    # Add a scalar coord
-    realization = iris.coords.AuxCoord(1, long_name="realization")
-    cube.add_aux_coord(realization)
-    # Check it's present and scalar
-    assert cube.coords("realization")
-    assert cube.coord_dims("realization") == ()
-    # Run function
-    out = misc.remove_scalar_coords(cube, ["realization"])
-    # Check it’s removed
-    cube_out = out[0]
-    assert not cube_out.coords("realization")
-
-
-def test_not_remove_non_scalar_coord():
-    """Test that non-scalar coordinate is not removed."""
-    # Create 1D cube
-    data = np.arange(5)
-    time = iris.coords.DimCoord(
-        np.arange(5), standard_name="time", units="hours since 1970-01-01"
+    return iris.cube.Cube(
+        data,
+        dim_coords_and_dims=[(lat, 0), (lon, 1)],
+        units=units,
     )
-    cube = iris.cube.Cube(data, dim_coords_and_dims=[(time, 0)])
-    # Confirm it's non-scalar
-    assert cube.coord_dims("time") != ()
-    # Run function
-    out = misc.remove_scalar_coords(cube, ["time"])
-    # Check it is still present
-    cube_out = out[0]
-    assert cube_out.coords("time")
+
+
+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
+
+
+@pytest.mark.parametrize(
+    "sentinel",
+    [1e10, 1e11, 999999, -999999],
+)
+def test_mask_fill_value_sentinels(sentinel):
+    """Test known sentinel values are converted to NaN."""
+    cube = _make_cube([[1.0, sentinel, 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():
+    """Masked values should become NaNs."""
+    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
+    result = np.ma.filled(result, np.nan)
+
+    assert result[0, 0] == 1.0
+    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)