Repository: astropy/astropy

Issue:
Coordinate method fails if `SkyCoord` contains `MaskedQuantity` fields
<!-- This comments are hidden when you submit the issue,
so you do not need to remove them! -->

<!-- Please be sure to check out our contributing guidelines,
https://github.com/astropy/astropy/blob/main/CONTRIBUTING.md .
Please be sure to check out our code of conduct,
https://github.com/astropy/astropy/blob/main/CODE_OF_CONDUCT.md . -->

<!-- Please have a search on our GitHub repository to see if a similar
issue has already been posted.
If a similar issue is closed, have a quick look to see if you are satisfied
by the resolution.
If not please go ahead and open an issue! -->

<!-- Please check that the development version still produces the same bug.
You can install development version with
pip install git+https://github.com/astropy/astropy
command. -->

### Description
<!-- Provide a general description of the bug. -->
`apply_space_motion` on a coordinate fails, if that Coordinate contains `astropy.utils.masked.core.MaskedQuantity` fields. If, instead, the coordinate is only made with normal quantities everything works.
I admit, that the setup in the minimal example below looks a little obscure. That's only because I'm trying to give an example here that does not rely on other packages.  In practice, I'm getting the coordinates and proper motions from an astroquery call to Simbad. That results in a masked table. Then, I convert that result to a `QTable` because I need to do some unit-aware arithmetic on columns that are not directly related to this problem and so I end up with `MaskedQuantity` after conversion to the `QTable`. 

This particular problem in new in 5.0 because `MAskedColumns` were not converted  automatically to `MaskedQuantities` before.

### Expected behavior
<!-- What did you expect to happen. -->
`apply_space_motion` should work for any valid coordinate object, independent on how exactly that object was initialized. Note that in the minimal example below, I have two coordinates. In the first row, all numbers are provided, while in the second row some are masked. Yet, even a coordinate transform on the first row only fails. I would not be surprised to see the second row fail (although I personally would prefer a to assume that a masked quantity is treated as a no-op because that way I can run `apply_space_motion` on my entire catalog even if not all objects have known proper motions), but at the very least the first one should work.

### Actual behavior
<!-- What actually happened. -->
<!-- Was the output confusing or poorly described? -->
```
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [115], in <cell line: 1>()
----> 1 example[0]['coord'].apply_space_motion(new_obstime=Time('2020-10-22T12:23:45'))

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/astropy/coordinates/sky_coordinate.py:791, in SkyCoord.apply_space_motion(self, new_obstime, dt)
    788 except u.UnitConversionError:  # No RV
    789     rv = 0.
--> 791 starpm = erfa.pmsafe(icrsrep.lon.radian, icrsrep.lat.radian,
    792                      icrsvel.d_lon.to_value(u.radian/u.yr),
    793                      icrsvel.d_lat.to_value(u.radian/u.yr),
    794                      plx, rv, t1.jd1, t1.jd2, t2.jd1, t2.jd2)
    796 if parallax_zero:
    797     new_distance = None

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/erfa/core.py:4991, in pmsafe(ra1, dec1, pmr1, pmd1, px1, rv1, ep1a, ep1b, ep2a, ep2b)
   4855 def pmsafe(ra1, dec1, pmr1, pmd1, px1, rv1, ep1a, ep1b, ep2a, ep2b):
   4856     """
   4857     Star proper motion:  update star catalog data for space motion, with
   4858     special handling to handle the zero parallax case.
   (...)
   4989 
   4990     """
-> 4991     ra2, dec2, pmr2, pmd2, px2, rv2, c_retval = ufunc.pmsafe(
   4992         ra1, dec1, pmr1, pmd1, px1, rv1, ep1a, ep1b, ep2a, ep2b)
   4993     check_errwarn(c_retval, 'pmsafe')
   4994     return ra2, dec2, pmr2, pmd2, px2, rv2

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/astropy/utils/masked/core.py:743, in MaskedNDArray.__array_ufunc__(self, ufunc, method, *inputs, **kwargs)
    739         mask = result_masks if len(result_masks) > 1 else result_masks[0]
    741 elif method == '__call__':
    742     # Regular ufunc call.
--> 743     mask = self._combine_masks(masks, out=out_mask)
    745 elif method == 'outer':
    746     # Must have two arguments; adjust masks as will be done for data.
    747     assert len(masks) == 2

File /nfs/melkor/d1/guenther/soft/mambaforge/envs/ciao-4.14/lib/python3.9/site-packages/astropy/utils/masked/core.py:661, in MaskedNDArray._combine_masks(self, masks, out)
    659 out = np.logical_or(masks[0], masks[1], out=out)
    660 for mask in masks[2:]:
--> 661     np.logical_or(out, mask, out=out)
    662 return out

TypeError: return arrays must be of ArrayType


```

### Steps to Reproduce
<!-- Ideally a code example could be provided so we can run it ourselves. -->
<!-- If you are pasting code, use triple backticks (```) around
your code snippet. -->
<!-- If necessary, sanitize your screen output to be pasted so you do not
reveal secrets like tokens and passwords. -->

Note that the error does not appear, if I set `masked=False` in the creation of the table. In that case, the proper motion columns are converted `astropy.units.quantity.Quantity` in the `QTable` and masked fields are simply filled with nan.

```python
import numpy as np
from astropy.table import Table, QTable
from astropy.coordinates import SkyCoord
from astropy.time import Time
import astropy.units as u
example = Table({'ra': [1, 2] *u.deg, 'dec':[3,4] * u.deg,
                 'pm_ra_cosdec': np.ma.array([5, np.ma.masked])*u.mas/u.yr, 
                 'pm_dec': np.ma.array([1, np.ma.masked])*u.mas/u.yr,
                 'obstime': Time(['2000-10-22T12:23:45', '2000-10-22T12:23:45'])}, masked=True)
example = QTable(example)
example['coord'] = SkyCoord(**example)
example[0]['coord'].apply_space_motion(new_obstime=Time('2020-10-22T12:23:45'))
```

### System Details
<!-- Even if you do not think this is necessary, it is useful information for the maintainers.
Please run the following snippet and paste the output below:
import platform; print(platform.platform())
import sys; print("Python", sys.version)
import numpy; print("Numpy", numpy.__version__)
import erfa; print("pyerfa", erfa.__version__)
import astropy; print("astropy", astropy.__version__)
import scipy; print("Scipy", scipy.__version__)
import matplotlib; print("Matplotlib", matplotlib.__version__)
-->
Linux-3.10.0-1062.12.1.el7.x86_64-x86_64-with-glibc2.17
Python 3.9.12 | packaged by conda-forge | (main, Mar 24 2022, 23:25:59) 
[GCC 10.3.0]
Numpy 1.22.3
pyerfa 2.0.0.1
astropy 5.0.3
Scipy 1.8.0
Matplotlib 3.5.1

Discussion:
github-actions[bot]: Thank you for your contribution to Astropy! 🌌 This checklist is meant to remind the package maintainers who will review this pull request of some common things to look for.

  - [ ] Do the proposed changes actually accomplish desired goals?
  - [ ] Do the proposed changes follow the [Astropy coding guidelines](https://docs.astropy.org/en/latest/development/codeguide.html)?
  - [ ] Are tests added/updated as required? If so, do they follow the [Astropy testing guidelines](https://docs.astropy.org/en/latest/development/testguide.html)?
  - [ ] Are docs added/updated as required? If so, do they follow the [Astropy documentation guidelines](https://docs.astropy.org/en/latest/development/docguide.html#astropy-documentation-rules-and-guidelines)?
  - [ ] Is rebase and/or squash necessary? If so, please provide the author with appropriate instructions. Also see instructions for [rebase](https://docs.astropy.org/en/latest/development/workflow/development_workflow.html#rebase-if-necessary) and [squash](https://docs.astropy.org/en/latest/development/workflow/development_workflow.html#squash-if-necessary).
  - [ ] Did the CI pass? If no, are the failures related? If you need to run daily and weekly cron jobs as part of the PR, please apply the "Extra CI" label. Codestyle issues can be fixed by the [bot](https://docs.astropy.org/en/latest/development/workflow/development_workflow.html#pre-commit).
  - [ ] Is a change log needed? If yes, did the change log check pass? If no, add the "no-changelog-entry-needed" label. If this is a manual backport, use the "skip-changelog-checks" label unless special changelog handling is necessary.
  - [ ] Is this a big PR that makes a "What's new?" entry worthwhile and if so, is (1) a "what's new" entry included in this PR and (2) the "whatsnew-needed" label applied?
  - [ ] At the time of adding the milestone, if the milestone set requires a backport to release branch(es), apply the appropriate "backport-X.Y.x" label(s) *before* merge.
neutrinoceros: Sorry I've been a bit out of the loop last week, so I'm just now seeing your review request. I'll do it now !
mhvk: @neutrinoceros - I think I have addressed all your comments - with the PR clearly the better for it, thanks!

Patch:
diff --git a/astropy/utils/masked/core.py b/astropy/utils/masked/core.py
index 89b0f528b1..22efe15fe1 100644
--- a/astropy/utils/masked/core.py
+++ b/astropy/utils/masked/core.py
@@ -1017,6 +1017,18 @@ class MaskedNDArray(Masked, np.ndarray, base_cls=np.ndarray, data_cls=np.ndarray
         # For inplace, the mask will have been set already.
         return out
 
+    def __array_wrap__(self, obj, context=None, return_scalar=False):
+        if context is None:
+            # Functions like np.ediff1d call __array_wrap__ to turn the array
+            # into self's subclass.
+            return self.from_unmasked(*self._get_data_and_mask(obj))
+
+        raise NotImplementedError(
+            "__array_wrap__ should not be used with a context any more since all use "
+            "should go through array_function. Please raise an issue on "
+            "https://github.com/astropy/astropy"
+        )
+
     # Below are ndarray methods that need to be overridden as masked elements
     # need to be skipped and/or an initial value needs to be set.
     def _reduce_defaults(self, kwargs, initial_func=None):
@@ -1291,6 +1303,14 @@ class MaskedNDArray(Masked, np.ndarray, base_cls=np.ndarray, data_cls=np.ndarray
         else:
             return string
 
+    def __hash__(self):
+        # Try to be somewhat like a numpy array scalar if possible.
+        if self.ndim == 0 and not self.mask:
+            return hash(self.unmasked[()])
+
+        # Will raise regular ndarray error.
+        return hash((self.unmasked, self.mask))
+
 
 class MaskedRecarray(np.recarray, MaskedNDArray, data_cls=np.recarray):
     # Explicit definition since we need to override some methods.
diff --git a/astropy/utils/masked/function_helpers.py b/astropy/utils/masked/function_helpers.py
index c56415cc67..de84a3394e 100644
--- a/astropy/utils/masked/function_helpers.py
+++ b/astropy/utils/masked/function_helpers.py
@@ -2,7 +2,7 @@
 """Helpers for letting numpy functions interact with Masked arrays.
 
 The module supplies helper routines for numpy functions that propagate
-masks appropriately., for use in the ``__array_function__``
+masks appropriately, for use in the ``__array_function__``
 implementation of `~astropy.utils.masked.MaskedNDArray`.  They are not
 very useful on their own, but the ones with docstrings are included in
 the documentation so that there is a place to find out how the mask is
@@ -129,6 +129,8 @@ MASKED_SAFE_FUNCTIONS |= {
     np.fix, np.isneginf, np.isposinf,
     # np.lib._function_base_impl
     np.angle, np.i0,
+    # np.lib._arraysetops_impl
+    np.intersect1d, np.setxor1d, np.union1d, np.unique,
 }  # fmt: skip
 
 
@@ -141,7 +143,10 @@ if NUMPY_LT_2_0:
     MASKED_SAFE_FUNCTIONS |= {np.trapz}
 else:
     # new in numpy 2.0
-    MASKED_SAFE_FUNCTIONS |= {np.astype, np.trapezoid}
+    MASKED_SAFE_FUNCTIONS |= {
+        np.astype, np.trapezoid,
+        np.unique_all, np.unique_counts, np.unique_inverse, np.unique_values,
+    }  # fmt: skip
 
 IGNORED_FUNCTIONS = {
     # I/O - useless for Masked, since no way to store the mask.
@@ -164,16 +169,6 @@ IGNORED_FUNCTIONS |= {
     np.einsum, np.einsum_path,
 }  # fmt: skip
 
-# Really should do these...
-if NUMPY_LT_2_0:
-    from numpy.lib import arraysetops
-else:
-    # Public set operations have been moved to the top-level namespace in numpy 2.0
-    # (numpy/numpy#24507), raising an AttributeError when accessed through np.lib.arraysetops.
-    from numpy.lib import _arraysetops_impl as arraysetops
-
-IGNORED_FUNCTIONS |= {getattr(np, setopsname) for setopsname in arraysetops.__all__}
-
 # Explicitly unsupported functions
 UNSUPPORTED_FUNCTIONS |= {
     np.unravel_index,
@@ -1334,6 +1329,93 @@ for nanfuncname in _nplibnanfunctions.__all__:
     )
 
 
+@dispatched_function
+def ediff1d(ary, to_end=None, to_begin=None):
+    from astropy.utils.masked import Masked
+
+    # ediff1d works fine if ary is Masked, but not if it is not (and
+    # we got here because to_end and/or to_begin are Masked).
+    if not isinstance(ary, Masked):
+        ary = Masked(ary)
+
+    return np.ediff1d.__wrapped__(ary, to_end, to_begin), None, None
+
+
+def _in1d(ar1, ar2, assume_unique=False, invert=False, *, kind=None):
+    # Copy sorting implementation from _arraysetops_impl._in1d;
+    # the others cannot work in the presence of a mask.
+    if kind == "table":
+        raise ValueError(
+            "The 'table' method is not supported for Masked arrays."
+            "Please select 'sort' or None for kind."
+        )
+
+    # Straight copy from here on.
+    if not assume_unique:
+        ar1, rev_idx = np.unique(ar1, return_inverse=True)
+        ar2 = np.unique(ar2)
+
+    ar = np.concatenate((ar1, ar2))
+    # We need this to be a stable sort, so always use 'mergesort'
+    # here. The values from the first array should always come before
+    # the values from the second array.
+    order = ar.argsort(kind="mergesort")
+    sar = ar[order]
+    if invert:
+        bool_ar = sar[1:] != sar[:-1]
+    else:
+        bool_ar = sar[1:] == sar[:-1]
+    flag = np.concatenate((bool_ar, [invert]))
+    ret = np.empty(ar.shape, dtype=bool)
+    ret[order] = flag
+    return ret[: len(ar1)] if assume_unique else ret[rev_idx]
+
+
+def _copy_of_mask(a):
+    mask = getattr(a, "mask", None)
+    return mask.copy() if mask is not None else False
+
+
+if NUMPY_LT_1_24:  # "kind" argument introduced in 1.24.
+
+    @dispatched_function
+    def in1d(ar1, ar2, assume_unique=False, invert=False):
+        mask = _copy_of_mask(ar1).ravel()
+        return _in1d(ar1, ar2, assume_unique, invert), mask, None
+
+    @dispatched_function
+    def isin(element, test_elements, assume_unique=False, invert=False):
+        element = np.asanyarray(element)
+        result = _in1d(element, test_elements, assume_unique, invert)
+        result.shape = element.shape
+        return result, _copy_of_mask(element), None
+
+else:
+
+    @dispatched_function
+    def in1d(ar1, ar2, assume_unique=False, invert=False, *, kind=None):
+        mask = _copy_of_mask(ar1).ravel()
+        return _in1d(ar1, ar2, assume_unique, invert, kind=kind), mask, None
+
+    @dispatched_function
+    def isin(element, test_elements, assume_unique=False, invert=False, *, kind=None):
+        element = np.asanyarray(element)
+        result = _in1d(element, test_elements, assume_unique, invert, kind=kind)
+        result.shape = element.shape
+        return result, _copy_of_mask(element), None
+
+
+@dispatched_function
+def setdiff1d(ar1, ar2, assume_unique=False):
+    # Again, mostly just to avoid an asarray call.
+    if assume_unique:
+        ar1 = np.asanyarray(ar1).ravel()
+    else:
+        ar1 = np.unique(ar1)
+        ar2 = np.unique(ar2)
+    return ar1[np.isin(ar1, ar2, assume_unique=True, invert=True)], None, None
+
+
 # Add any dispatched or helper function that has a docstring to
 # __all__, so they will be typeset by sphinx. The logic is that for
 # those presumably the use of the mask is not entirely obvious.
diff --git a/docs/changes/time/16224.bugfix.rst b/docs/changes/time/16224.bugfix.rst
new file mode 100644
index 0000000000..67cee83f64
--- /dev/null
+++ b/docs/changes/time/16224.bugfix.rst
@@ -0,0 +1,3 @@
+Scalar ``Time`` instances are now hashable if they are not masked, also if one
+uses ``Masked`` internally, matching the behaviour prior to astropy 6.0 (and
+the current behaviour when masking using ``np.ma.MaskedArray``).
diff --git a/docs/changes/utils/16224.api.rst b/docs/changes/utils/16224.api.rst
new file mode 100644
index 0000000000..8f2830838e
--- /dev/null
+++ b/docs/changes/utils/16224.api.rst
@@ -0,0 +1,3 @@
+Unmasked ``Masked`` scalar instances are now considered hashable, to match the
+implicit behaviour of regular arrays, where if an operation leads to a scalar,
+a hashable array scalar is returned.
diff --git a/docs/changes/utils/16224.feature.rst b/docs/changes/utils/16224.feature.rst
new file mode 100644
index 0000000000..33f0d71ac3
--- /dev/null
+++ b/docs/changes/utils/16224.feature.rst
@@ -0,0 +1,2 @@
+``Masked`` instances now support the various numpy array set operations, such
+as ``np.unique`` and ``np.isin``.


Test patch:
diff --git a/astropy/time/tests/test_basic.py b/astropy/time/tests/test_basic.py
index a953974e74..990d4827f4 100644
--- a/astropy/time/tests/test_basic.py
+++ b/astropy/time/tests/test_basic.py
@@ -2394,16 +2394,13 @@ def test_hash_masked_time(location, masked_array_type):
     t = Time([1, 1, 2, 3], format="cxcsec", location=location)
     t[3] = np.ma.masked
     with conf.set_temp("masked_array_type", masked_array_type):
-        if masked_array_type == "numpy":
-            h1 = hash(t[0])
-            h2 = hash(t[1])
-            h3 = hash(t[2])
-            assert h1 == h2
-            assert h1 != h3
-        else:
-            with pytest.raises(TypeError, match="value is masked"):
-                hash(t[0])
-
+        # Unmasked scalars should always be fine.
+        h1 = hash(t[0])
+        h2 = hash(t[1])
+        assert h2 == h1
+        h3 = hash(t[2])
+        assert h3 != h1
+        # But arrays and masked elements cannot be hashed
         with pytest.raises(
             TypeError, match=r"unhashable type: 'Time' \(must be scalar\)"
         ):
@@ -2420,16 +2417,13 @@ def test_hash_time_delta_masked(masked_array_type):
     t = TimeDelta([1, 1, 2, 3], format="sec")
     t[3] = np.ma.masked
     with conf.set_temp("masked_array_type", masked_array_type):
-        if masked_array_type == "numpy":
-            h1 = hash(t[0])
-            h2 = hash(t[1])
-            h3 = hash(t[2])
-            assert h1 == h2
-            assert h1 != h3
-        else:
-            with pytest.raises(TypeError, match=r"'TimeDelta' \(value is masked\)"):
-                hash(t[0])
-
+        # Unmasked scalars should always be fine.
+        h1 = hash(t[0])
+        h2 = hash(t[1])
+        h3 = hash(t[2])
+        assert h2 == h1
+        assert h3 != h1
+        # But arrays and masked elements cannot be hashed
         with pytest.raises(TypeError, match=r"'TimeDelta' \(must be scalar\)"):
             hash(t)
 
diff --git a/astropy/utils/masked/tests/test_containers.py b/astropy/utils/masked/tests/test_containers.py
index c5c2e48ef0..0ac5ea624f 100644
--- a/astropy/utils/masked/tests/test_containers.py
+++ b/astropy/utils/masked/tests/test_containers.py
@@ -81,6 +81,30 @@ class TestSkyCoord:
         assert_array_equal(mgcrs.data.lon.unmasked, gcrs.data.lon)
         assert_array_equal(mgcrs.data.lat.unmasked, gcrs.data.lat)
 
+    @pytest.mark.filterwarnings("ignore:.*ERFA.*distance overridden.*")
+    def test_apply_space_motion(self):
+        # Regression test for gh-13041. It is important that the
+        # distance is missing here, so that a warning is raised.
+        # Before gh-16224, the ERFA warning machinery let to an error.
+        kwargs = dict(
+            ra=[1, 2] * u.deg,
+            dec=[3, 4] * u.deg,
+            pm_ra_cosdec=[5, 6] * u.mas / u.yr,
+            pm_dec=[1, 2] * u.mas / u.yr,
+            obstime=Time(["2000-10-22T12:23:45", "2000-10-22T12:23:45"]),
+        )
+        sc = SkyCoord(**kwargs)
+        mask = np.array([False, True])
+        kwargs["pm_ra_cosdec"] = Masked(kwargs["pm_ra_cosdec"], mask=mask)
+        kwargs["pm_dec"] = Masked(kwargs["pm_dec"], mask=mask)
+        msc = SkyCoord(**kwargs)
+        new_time = Time("2020-10-22T12:23:45")
+        sc2 = sc.apply_space_motion(new_obstime=new_time)
+        msc2 = msc.apply_space_motion(new_obstime=new_time)
+        assert_array_equal(msc2.data.lon.mask, [False, True])
+        assert_array_equal(msc2.data.lon.unmasked, sc2.data.lon)
+        assert_array_equal(msc2.data.lat.unmasked, sc2.data.lat)
+
 
 class TestTime:
     def setup_class(self):
diff --git a/astropy/utils/masked/tests/test_function_helpers.py b/astropy/utils/masked/tests/test_function_helpers.py
index 73086bf3f7..3d797ffb39 100644
--- a/astropy/utils/masked/tests/test_function_helpers.py
+++ b/astropy/utils/masked/tests/test_function_helpers.py
@@ -1454,6 +1454,198 @@ class TestNaNFunctions:
         self.check(np.nanpercentile, q=50)
 
 
+class TestArraySetOps:
+    """Tests based on those from numpy.ma.tests.test_extras.
+
+    Adjusted to take into account that comparing masked values should
+    result in masked equality.
+
+    """
+
+    @classmethod
+    def setup_class(self):
+        # Setup for unique (names as in unique_all NamedTuple)
+        # input data, unique values, indices in data to those,
+        # inverse indices in values to reconstruct data, counts.
+        self.data = Masked([1, 1, 1, 2, 2, 3], mask=[0, 0, 1, 0, 1, 0])
+        self.values = Masked([1, 2, 3, 1, 2], mask=[0, 0, 0, 1, 1])
+        self.indices = np.array([0, 3, 5, 2, 4])
+        self.inverse_indices = np.array([0, 0, 3, 1, 4, 2])
+        self.counts = np.array([2, 1, 1, 1, 1])
+
+    @pytest.mark.parametrize("dtype", [int, float, object])
+    def test_unique(self, dtype):
+        values, indices, inverse_indices = np.unique(
+            self.data.astype(dtype), return_index=True, return_inverse=True
+        )
+        assert_masked_equal(values, self.values.astype(dtype))
+        assert_array_equal(indices, self.indices)
+        assert_array_equal(inverse_indices, self.inverse_indices)
+        # All masked
+        data2 = Masked([2, 1, 3], mask=True)
+        values2, indices2, inverse_indices2 = np.unique(
+            data2.astype(dtype), return_index=True, return_inverse=True
+        )
+        expected_values2 = Masked([1, 2, 3], mask=True)
+        assert_masked_equal(values2, expected_values2.astype(dtype))
+        assert_array_equal(indices2, [1, 0, 2])
+        assert_array_equal(inverse_indices2, [1, 0, 2])
+
+    @pytest.mark.skipif(NUMPY_LT_2_0, reason="new in numpy 2.0")
+    def check_unique(self, test):
+        for name in test._fields:
+            assert_array_equal(getattr(test, name), getattr(self, name))
+
+    @pytest.mark.skipif(NUMPY_LT_2_0, reason="new in numpy 2.0")
+    def test_unique_all(self):
+        test = np.unique_all(self.data)
+        assert len(test) == 4
+        self.check_unique(test)
+
+    @pytest.mark.skipif(NUMPY_LT_2_0, reason="new in numpy 2.0")
+    def test_unique_counts(self):
+        test = np.unique_counts(self.data)
+        assert len(test) == 2
+        self.check_unique(test)
+
+    @pytest.mark.skipif(NUMPY_LT_2_0, reason="new in numpy 2.0")
+    def test_unique_inverse(self):
+        test = np.unique_inverse(self.data)
+        assert len(test) == 2
+        self.check_unique(test)
+
+    @pytest.mark.skipif(NUMPY_LT_2_0, reason="new in numpy 2.0")
+    def test_unique_values(self):
+        test = np.unique_values(self.data)
+        assert isinstance(test, Masked)
+        assert_array_equal(test, self.values)
+
+    def test_ediff1d(self):
+        x = Masked(np.arange(5), mask=[1, 0, 0, 0, 1])
+        control = Masked([1, 1, 1, 1], mask=[1, 0, 0, 1])
+        test = np.ediff1d(x)
+        assert_masked_equal(test, control)
+        # Test ediff1d w/ to_begin
+        test2 = np.ediff1d(x, to_begin=Masked(10, mask=True))
+        control2 = Masked([10, 1, 1, 1, 1], mask=[1, 1, 0, 0, 1])
+        assert_masked_equal(test2, control2)
+        test3 = np.ediff1d(x, to_begin=[1, 2, 3])
+        control3 = Masked([1, 2, 3, 1, 1, 1, 1], mask=[0, 0, 0, 1, 0, 0, 1])
+        assert_masked_equal(test3, control3)
+        # Test ediff1d w/ to_end
+        test4 = np.ediff1d(x, to_end=Masked(10, mask=True))
+        control4 = Masked([1, 1, 1, 1, 10], mask=[1, 0, 0, 1, 1])
+        assert_masked_equal(test4, control4)
+        test5 = np.ediff1d(x, to_end=[1, 2, 3])
+        control5 = Masked([1, 1, 1, 1, 1, 2, 3], mask=[1, 0, 0, 1, 0, 0, 0])
+        assert_masked_equal(test5, control5)
+        # Test ediff1d w/ to_begin and to_end
+        test6 = np.ediff1d(
+            x, to_end=Masked(10, mask=True), to_begin=Masked(20, mask=True)
+        )
+        control6 = Masked([20, 1, 1, 1, 1, 10], mask=[1, 1, 0, 0, 1, 1])
+        assert_masked_equal(test6, control6)
+        test7 = np.ediff1d(x, to_end=[1, 2, 3], to_begin=Masked(10, mask=True))
+        control7 = Masked([10, 1, 1, 1, 1, 1, 2, 3], mask=[1, 1, 0, 0, 1, 0, 0, 0])
+        assert_masked_equal(test7, control7)
+        # Test ediff1d w/ a ndarray
+        test8 = np.ediff1d(
+            np.arange(5), to_end=Masked(10, mask=True), to_begin=Masked(20, mask=True)
+        )
+        control8 = Masked([20, 1, 1, 1, 1, 10], mask=[1, 0, 0, 0, 0, 1])
+        assert_masked_equal(test8, control8)
+
+    def test_intersect1d(self):
+        x = Masked([1, 3, 3, 3, 4], mask=[0, 0, 0, 1, 1])
+        y = Masked([3, 1, 1, 1, 4], mask=[0, 0, 0, 1, 1])
+        test = np.intersect1d(x, y)
+        control = Masked([1, 3, 4], mask=[0, 0, 1])
+        assert_masked_equal(test, control)
+
+    def test_setxor1d(self):
+        a = Masked([1, 2, 5, 7, -1], mask=[0, 0, 0, 0, 1])
+        b = Masked([1, 2, 3, 4, 5, -1], mask=[0, 0, 0, 0, 0, 1])
+        test = np.setxor1d(a, b)
+        assert_masked_equal(test, Masked([3, 4, 7]))
+        a = Masked([1, 2, 5, 7, -1], mask=[0, 0, 0, 0, 1])
+        b = [1, 2, 3, 4, 5]
+        test = np.setxor1d(a, b)
+        assert_masked_equal(test, Masked([3, 4, 7, -1], mask=[0, 0, 0, 1]))
+        a = Masked([1, 8, 2, 3], mask=[0, 1, 0, 0])
+        b = Masked([6, 5, 4, 8], mask=[0, 0, 0, 1])
+        test = np.setxor1d(a, b)
+        assert_masked_equal(test, Masked([1, 2, 3, 4, 5, 6]))
+        #
+        assert_masked_equal(np.setxor1d(Masked([]), []), Masked([]))
+
+    @pytest.mark.parametrize("dtype", [int, float, object])
+    def test_isin(self, dtype):
+        a = np.arange(24).reshape((2, 3, 4))
+        mask = np.zeros(a.shape, bool)
+        mask[1, 2, 0] = 1  # 20
+        mask[1, 2, 1] = 1  # 21
+        a = Masked(a, mask=mask)
+        b = Masked([0, 10, 20, 30, 1, 3, 11, 21, 33], mask=[0, 1, 0, 1, 0, 1, 0, 1, 0])
+        # unmasked 0, 20, 1, 11, 33, masked 10, 30, 3, 21
+        ec = np.zeros((2, 3, 4), dtype=bool)
+        ec[0, 0, 0] = True  # 0
+        ec[0, 0, 1] = True  # 1
+        ec[0, 2, 3] = True  # 11
+        ec[1, 2, 1] = True  # masked 21
+        ec = Masked(ec, mask)
+        c = np.isin(a.astype(dtype), b.astype(dtype))
+        assert_masked_equal(c, ec)
+
+    @pytest.mark.filterwarnings("ignore:in1d.*deprecated")  # not NUMPY_LT_2_0
+    def test_in1d(self):
+        # Once we require numpy>=2.0, these tests should be joined with np.isin.
+        a = Masked([1, 2, 5, -2, -1], mask=[0, 0, 0, 1, 1])
+        b = Masked([1, 2, 3, 4, 5, -2], mask=[0, 0, 0, 0, 0, 1])
+        test = np.in1d(a, b)
+        assert_masked_equal(test, Masked([True, True, True, True, False], mask=a.mask))
+        assert_array_equal(np.in1d(a, b, invert=True), ~test)
+
+        a = Masked([5, 5, 2, -2, -1], mask=[0, 0, 0, 1, 1])
+        b = Masked([1, 5, -1], mask=[0, 0, 1])
+        test = np.in1d(a, b)
+        assert_masked_equal(test, Masked([True, True, False, False, True], mask=a.mask))
+
+        assert_masked_equal(np.in1d(Masked([]), []), Masked([]))
+        assert_masked_equal(np.in1d(Masked([]), [], invert=True), Masked([]))
+
+    @pytest.mark.skipif(NUMPY_LT_1_24, reason="kind introduced in numpy 1.24")
+    def test_in1d_kind_table_error(self):
+        with pytest.raises(ValueError, match="'table' method is not supported"):
+            np.in1d(Masked([1, 2, 3]), [4, 5], kind="table")
+
+    @pytest.mark.parametrize("dtype", [int, float, object])
+    def test_union1d(self, dtype):
+        a = Masked([1, 2, 5, 7, 5, 5], mask=[0, 0, 0, 0, 0, 1])
+        b = Masked([1, 2, 3, 4, 5, 6], mask=[0, 0, 0, 0, 0, 1])
+        control = Masked([1, 2, 3, 4, 5, 7, 5, 6], mask=[0, 0, 0, 0, 0, 0, 1, 1])
+        test = np.union1d(a.astype(dtype), b.astype(dtype))
+        assert_masked_equal(test, control.astype(dtype))
+
+        assert_masked_equal(np.union1d(Masked([]), []), Masked([]))
+
+    def test_setdiff1d(self):
+        a = Masked([6, 5, 4, 7, 7, 1, 2, 1], mask=[0, 0, 0, 0, 0, 0, 0, 1])
+        b = np.array([2, 4, 3, 3, 2, 1, 5])
+        test = np.setdiff1d(a, b)
+        assert_masked_equal(test, Masked([6, 7, 1], mask=[0, 0, 1]))
+        b2 = Masked(b, mask=[1, 1, 1, 1, 0, 0, 0])
+        test2 = np.setdiff1d(a, b2)
+        assert_masked_equal(test2, Masked([4, 6, 7, 1], mask=[0, 0, 0, 1]))
+
+        a = Masked(np.array([], dtype=np.uint32), mask=[])
+        assert np.setdiff1d(a, []).dtype == np.uint32
+
+        a = Masked(["a", "b", "c"], mask=[0, 1, 1])
+        b = Masked(["a", "b", "s"], mask=[0, 1, 1])
+        test3 = np.setdiff1d(a, b, assume_unique=True)
+        assert_masked_equal(test3, Masked(["c"], True))
+
+
 # Get wrapped and covered functions.
 all_wrapped_functions = get_wrapped_functions(np)
 tested_functions = get_covered_functions(locals())
diff --git a/astropy/utils/masked/tests/test_functions.py b/astropy/utils/masked/tests/test_functions.py
index 9a22a681d6..1ed5b8d555 100644
--- a/astropy/utils/masked/tests/test_functions.py
+++ b/astropy/utils/masked/tests/test_functions.py
@@ -17,6 +17,7 @@ from numpy.testing import assert_allclose, assert_array_equal
 
 from astropy import units as u
 from astropy.units import Quantity
+from astropy.utils import minversion
 from astropy.utils.compat.numpycompat import NUMPY_LT_1_25
 from astropy.utils.masked.core import Masked
 
@@ -598,15 +599,24 @@ class TestStructuredUfuncs:
 
 def test_erfa_no_warnings_on_masked_entries():
     # Erfa warns for invalid inputs for some routines.
-    ihour1 = [25, 10]
-    with pytest.warns(erfa.ErfaWarning, match="ihour outside range"):
+    msg = 'ERFA function "tf2d" yielded {count} of "ihour outside range 0-23"'
+    ihour1 = [25, 26, 10]
+    with pytest.warns(erfa.ErfaWarning, match=msg.format(count=2)):
         res1 = erfa.tf2d("+", ihour1, 0, 0.0)
     # But will not if they are masked.
-    mask = [True, False]
+    mask = [True, True, False]
     ihour2 = Masked(ihour1, mask)
     res2 = erfa.tf2d("+", ihour2, 0, 0.0)
     assert_array_equal(res2.unmasked, res1)
     assert_array_equal(res2.mask, mask)
+    # And will count correctly.
+    mask = [True, False, False]
+    ihour3 = Masked(ihour1, mask)
+    count = 1 if minversion(erfa, "2.0.1.1", inclusive=False) else "—"
+    with pytest.warns(erfa.ErfaWarning, match=msg.format(count=count)):
+        res3 = erfa.tf2d("+", ihour3, 0, 0.0)
+    assert_array_equal(res3.unmasked, res1)
+    assert_array_equal(res3.mask, mask)
 
 
 def test_erfa_no_exceptions_on_masked_entries():
diff --git a/astropy/utils/masked/tests/test_masked.py b/astropy/utils/masked/tests/test_masked.py
index 00b27b7526..b7ae99dcd5 100644
--- a/astropy/utils/masked/tests/test_masked.py
+++ b/astropy/utils/masked/tests/test_masked.py
@@ -662,6 +662,15 @@ class MaskedItemTests(MaskedArraySetup):
         assert_array_equal(base.unmasked, self.a)
         assert_array_equal(base.mask, expected_mask)
 
+    @pytest.mark.parametrize("item", VARIOUS_ITEMS)
+    def test_hash(self, item):
+        ma_part = self.ma[item]
+        if ma_part.ndim > 0 or ma_part.mask.any():
+            with pytest.raises(TypeError, match="unhashable"):
+                hash(ma_part)
+        else:
+            assert hash(ma_part) == hash(ma_part.unmasked)
+
 
 class TestMaskedArrayItems(MaskedItemTests):
     @classmethod
@@ -883,6 +892,16 @@ class TestMaskedArrayMethods(MaskedArraySetup):
         assert_array_equal(ma_sum.unmasked, expected_data)
         assert_array_equal(ma_sum.mask, expected_mask)
 
+    def test_sum_hash(self):
+        ma_sum = self.ma.sum()
+        assert ma_sum.mask
+        # Masked scalars cannot be hashed.
+        with pytest.raises(TypeError, match="unhashable"):
+            hash(ma_sum)
+        ma_sum2 = Masked(self.a).sum()
+        # But an unmasked scalar can.
+        assert hash(ma_sum2) == hash(self.a.sum())
+
     @pytest.mark.parametrize("axis", (0, 1, None))
     def test_cumsum(self, axis):
         ma_sum = self.ma.cumsum(axis)
@@ -953,6 +972,10 @@ class TestMaskedArrayMethods(MaskedArraySetup):
         assert_array_equal(ma_mean.unmasked, expected_data)
         assert_array_equal(ma_mean.mask, expected_mask)
 
+    def test_mean_hash(self):
+        ma_mean = self.ma.mean()
+        assert hash(ma_mean) == hash(ma_mean.unmasked[()])
+
     @pytest.mark.filterwarnings("ignore:.*encountered in.*divide")
     @pytest.mark.parametrize("axis", (0, 1, None))
     def test_var(self, axis):