Allow masked data in coords (#6468)

* Allow masked data to be passed into a Cell, and by extensions a Coord.

* Handled case where integer array with masked values could not be created.

* Handle case where reading of scalar masked value from netCDF file
looses the original dtype.

* Removed some comments

* Removed fill_value keyword to masked_array creation - wasn't fixing the
issue with numpy warnings about NaNs.

* Single method for handling points arrays with masked data

* Refactored to add small inline helper function to add AuxCoord

* Refactored the masked array handling and now filling masked arrays with default fill type.

* Added test for Cell hashing

* test_Cell.py: Converted unittest to pytest

* Added tests for merging of scalar aux coords with mising data.

* Removed unnecessary __main__ entry point from test_Cell.py

* Missing hash test for bounded points

* More rigorous test on masked data of merged AuxCoord

* Migrated test_merge.py to pytest.

* Moved tests/test_merge.py to tests/unit/merge/test_merge.py

* What's new entry.

* Fixed some typos in comment

* Split out integration test from unit/merge/test_merge.py to
integration/merge/test_merge.py

Author: ukmo-ccbunney

docs/src/whatsnew/latest.rst

@@ -44,6 +44,9 @@ This document explains the changes made to Iris for this release
    so it also converts the values of the attributes ``"actual_range"``,
    ``"valid_max"``, ``"valid_min"``, and ``"valid_range"``. (:pull:`6416`)
 
+#. `@ukmo-ccbunney`_ fixed loading and merging of masked data in scalar ``AuxCoords``.
+   (:issue:`3584`, :pull:`6468`)
+
 💣 Incompatible Changes
 =======================
 

lib/iris/_merge.py

@@ -1406,10 +1406,37 @@ def axis_and_name(name):
                     # TODO: Consider appropriate sort order (ascending,
                     # descending) i.e. use CF positive attribute.
                     cells = sorted(indexes[name])
-                    points = np.array(
-                        [cell.point for cell in cells],
-                        dtype=metadata[name].points_dtype,
-                    )
+                    points = [cell.point for cell in cells]
+
+                    # If any points are masked then create a masked array type,
+                    # otherwise create a standard ndarray.
+                    if np.ma.masked in points:
+                        dtype = metadata[name].points_dtype
+
+                        # Create a pre-filled array with all elements set to `fill_value` for dtype
+                        # This avoids the following problems when trying to do `np.ma.masked_array(points, dtype...)`:
+                        #   - Underlying data of masked elements is arbitrary
+                        #   - Can't convert a np.ma.masked to an integer type
+                        #   - For floating point arrays, numpy raises a warning about "converting masked elements to NaN"
+                        fill_value = np.trunc(
+                            np.ma.default_fill_value(dtype), dtype=dtype
+                        )  # truncation needed to deal with silly default fill values in Numpy
+
+                        # create array of fill values; ensures we have consistent data under mask
+                        arr_points = np.ma.repeat(dtype.type(fill_value), len(points))
+
+                        # get mask index and filtered data then store in new array:
+                        mask = np.array([p is np.ma.masked for p in points])
+                        arr_points.mask = mask
+
+                        # Need another list comprehension to avoid numpy warning "converting masked elements to NaN":
+                        arr_points[~mask] = np.array(
+                            [p for p in points if p is not np.ma.masked]
+                        )
+                        points = arr_points
+                    else:
+                        points = np.array(points, dtype=metadata[name].points_dtype)
+
                     if cells[0].bound is not None:
                         bounds = np.array(
                             [cell.bound for cell in cells],
@@ -1594,13 +1621,19 @@ def _build_coordinates(self):
             # the bounds are not monotonic, so try building the coordinate,
             # and if it fails make the coordinate into an auxiliary coordinate.
             # This will ultimately make an anonymous dimension.
-            try:
-                coord = iris.coords.DimCoord(
-                    template.points, bounds=template.bounds, **template.kwargs
-                )
-                dim_coords_and_dims.append(_CoordAndDims(coord, template.dims))
-            except ValueError:
+
+            # If the points contain masked values, if definitely cannot be built
+            # as a dim coord, so add it to the _aux_templates immediately
+            if np.ma.is_masked(template.points):
                 self._aux_templates.append(template)
+            else:
+                try:
+                    coord = iris.coords.DimCoord(
+                        template.points, bounds=template.bounds, **template.kwargs
+                    )
+                    dim_coords_and_dims.append(_CoordAndDims(coord, template.dims))
+                except ValueError:
+                    self._aux_templates.append(template)
 
         # There is the potential that there are still anonymous dimensions.
         # Get a list of the dimensions which are not anonymous at this stage.
@@ -1609,26 +1642,34 @@ def _build_coordinates(self):
         ]
 
         # Build the auxiliary coordinates.
+        def _build_aux_coord_from_template(template):
+            # kwarg not applicable to AuxCoord.
+            template.kwargs.pop("circular", None)
+            coord = iris.coords.AuxCoord(
+                template.points, bounds=template.bounds, **template.kwargs
+            )
+            aux_coords_and_dims.append(_CoordAndDims(coord, template.dims))
+
         for template in self._aux_templates:
             # Attempt to build a DimCoord and add it to the cube. If this
             # fails e.g it's non-monontic or multi-dimensional or non-numeric,
             # then build an AuxCoord.
-            try:
-                coord = iris.coords.DimCoord(
-                    template.points, bounds=template.bounds, **template.kwargs
-                )
-                if len(template.dims) == 1 and template.dims[0] not in covered_dims:
-                    dim_coords_and_dims.append(_CoordAndDims(coord, template.dims))
-                    covered_dims.append(template.dims[0])
-                else:
-                    aux_coords_and_dims.append(_CoordAndDims(coord, template.dims))
-            except ValueError:
-                # kwarg not applicable to AuxCoord.
-                template.kwargs.pop("circular", None)
-                coord = iris.coords.AuxCoord(
-                    template.points, bounds=template.bounds, **template.kwargs
-                )
-                aux_coords_and_dims.append(_CoordAndDims(coord, template.dims))
+
+            # Check here whether points are masked? If so then it has to be an AuxCoord
+            if np.ma.is_masked(template.points):
+                _build_aux_coord_from_template(template)
+            else:
+                try:
+                    coord = iris.coords.DimCoord(
+                        template.points, bounds=template.bounds, **template.kwargs
+                    )
+                    if len(template.dims) == 1 and template.dims[0] not in covered_dims:
+                        dim_coords_and_dims.append(_CoordAndDims(coord, template.dims))
+                        covered_dims.append(template.dims[0])
+                    else:
+                        aux_coords_and_dims.append(_CoordAndDims(coord, template.dims))
+                except ValueError:
+                    _build_aux_coord_from_template(template)
 
         # Mix in the vector coordinates.
         for item, dims in zip(

lib/iris/coords.py

@@ -1237,6 +1237,9 @@ class Cell(namedtuple("Cell", ["point", "bound"])):
     # Make this class's comparison operators override those of numpy
     __array_priority__ = 100
 
+    # pre-computed hash for un-hashable `np.ma.masked` value
+    _MASKED_VALUE_HASH = hash("<<##MASKED_VALUE##>>")
+
     def __new__(cls, point=None, bound=None):
         """Construct a Cell from point or point-and-bound information."""
         if point is None:
@@ -1277,13 +1280,17 @@ def __add__(self, mod):
 
     def __hash__(self):
         # See __eq__ for the definition of when two cells are equal.
+        point = self.point
+        if np.ma.is_masked(point):
+            # `np.ma.masked` is unhashable
+            point = Cell._MASKED_VALUE_HASH
         if self.bound is None:
-            return hash(self.point)
+            return hash(point)
         bound = self.bound
         rbound = bound[::-1]
         if rbound < bound:
             bound = rbound
-        return hash((self.point, bound))
+        return hash((point, bound))
 
     def __eq__(self, other):
         """Compare Cell equality depending on the type of the object to be compared."""
@@ -2086,7 +2093,8 @@ def cell(self, index):
         """
         index = iris.util._build_full_slice_given_keys(index, self.ndim)
 
-        point = tuple(np.array(self.core_points()[index], ndmin=1).flatten())
+        # Use `np.asanyaray` to preserve any masked values:
+        point = tuple(np.asanyarray(self.core_points()[index]).flatten())
         if len(point) != 1:
             raise IndexError(
                 "The index %s did not uniquely identify a single "
@@ -2809,6 +2817,9 @@ def _values(self, points):
         # Check validity requirements for dimension-coordinate points.
         self._new_points_requirements(points)
         # Cast to a numpy array for masked arrays with no mask.
+
+        # NOTE: This is the point where any mask is lost on a coordinate if none of the
+        # values are actually masked. What if we wanted this to be an AuxCoord with a mask?
         points = np.array(points)
 
         super(DimCoord, self.__class__)._values.fset(self, points)

lib/iris/fileformats/netcdf/loader.py

@@ -253,6 +253,14 @@ def _get_cf_var_data(cf_var, filename):
         if total_bytes < _LAZYVAR_MIN_BYTES:
             # Don't make a lazy array, as it will cost more memory AND more time to access.
             result = cf_var[:]
+
+            # Special handling of masked scalar value; this will be returned as
+            # an `np.ma.masked` instance which will lose the original dtype.
+            # Workaround for this it return a 1-element masked array of the
+            # correct dtype. Note: this is not an issue for masked arrays,
+            # only masked scalar values.
+            if result is np.ma.masked:
+                result = np.ma.masked_all(1, dtype=cf_var.datatype)
         else:
             # Get lazy chunked data out of a cf variable.
             # Creates Dask wrappers around data arrays for any cube components which