Repository: pydata/xarray

Issue:
AttributeError: `np.unicode_` was removed in the NumPy 2.0 release. Use `np.str_` instead
### What happened?

We've been getting this AttributeError on PyGMT's dev CI which tests nightly packages of numpy/xarray/etc. See https://github.com/GenericMappingTools/pygmt/actions/runs/6043991082/job/16401918043#step:16:40

### What did you expect to happen?

Running `import xarray as xr` should work with numpy 2.0.

### Minimal Complete Verifiable Example

```Python
mamba create --name xr_with_np2 python=3.11
mamba activate xr_with_np2
python -m pip install --pre --extra-index https://pypi.anaconda.org/scientific-python-nightly-wheels/simple numpy pandas xarray

python -c "import xarray as xr"
```


### MVCE confirmation

- [X] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
- [X] Complete example — the example is self-contained, including all data and the text of any traceback.
- [X] Verifiable example — the example copy & pastes into an IPython prompt or [Binder notebook](https://mybinder.org/v2/gh/pydata/xarray/main?urlpath=lab/tree/doc/examples/blank_template.ipynb), returning the result.
- [X] New issue — a search of GitHub Issues suggests this is not a duplicate.

### Relevant log output

```Python
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/home/user/mambaforge/envs/xr_with_np2/lib/python3.11/site-packages/xarray/__init__.py", line 1, in <module>
    from xarray import testing, tutorial
  File "/home/user/mambaforge/envs/xr_with_np2/lib/python3.11/site-packages/xarray/testing.py", line 10, in <module>
    from xarray.core import duck_array_ops, formatting, utils
  File "/home/user/mambaforge/envs/xr_with_np2/lib/python3.11/site-packages/xarray/core/duck_array_ops.py", line 36, in <module>
    from xarray.core import dask_array_ops, dtypes, nputils
  File "/home/user/mambaforge/envs/xr_with_np2/lib/python3.11/site-packages/xarray/core/dask_array_ops.py", line 3, in <module>
    from xarray.core import dtypes, nputils
  File "/home/user/mambaforge/envs/xr_with_np2/lib/python3.11/site-packages/xarray/core/dtypes.py", line 43, in <module>
    (np.bytes_, np.unicode_),  # numpy promotes to unicode
                ^^^^^^^^^^^
  File "/home/user/mambaforge/envs/xr_with_np2/lib/python3.11/site-packages/numpy/__init__.py", line 372, in __getattr__
    raise AttributeError(
AttributeError: `np.unicode_` was removed in the NumPy 2.0 release. Use `np.str_` instead.
```


### Anything else we need to know?

See https://numpy.org/devdocs/release/2.0.0-notes.html#numpy-2-0-python-api-removals. Relevant code at xarray is here https://github.com/pydata/xarray/blob/1043a9e13574e859ec08d19425341b2e359d2802/xarray/core/dtypes.py#L36-L44

### Environment

<details>

```
# packages in environment at /home/weiji/mambaforge/envs/xr_with_np2:
#
# Name                    Version                   Build  Channel
_libgcc_mutex             0.1                 conda_forge    conda-forge
_openmp_mutex             4.5                       2_gnu    conda-forge
bzip2                     1.0.8                h7f98852_4    conda-forge
ca-certificates           2023.7.22            hbcca054_0    conda-forge
ld_impl_linux-64          2.40                 h41732ed_0    conda-forge
libexpat                  2.5.0                hcb278e6_1    conda-forge
libffi                    3.4.2                h7f98852_5    conda-forge
libgcc-ng                 13.1.0               he5830b7_0    conda-forge
libgomp                   13.1.0               he5830b7_0    conda-forge
libnsl                    2.0.0                h7f98852_0    conda-forge
libsqlite                 3.43.0               h2797004_0    conda-forge
libuuid                   2.38.1               h0b41bf4_0    conda-forge
libzlib                   1.2.13               hd590300_5    conda-forge
ncurses                   6.4                  hcb278e6_0    conda-forge
numpy                     2.0.0.dev0               pypi_0    pypi
openssl                   3.1.2                hd590300_0    conda-forge
packaging                 23.1                     pypi_0    pypi
pandas                    2.2.0.dev0+141.g413d92d471          pypi_0    pypi
pip                       23.2.1             pyhd8ed1ab_0    conda-forge
python                    3.11.5          hab00c5b_0_cpython    conda-forge
python-dateutil           2.8.2                    pypi_0    pypi
pytz                      2023.3                   pypi_0    pypi
readline                  8.2                  h8228510_1    conda-forge
setuptools                68.1.2             pyhd8ed1ab_0    conda-forge
six                       1.16.0                   pypi_0    pypi
tk                        8.6.12               h27826a3_0    conda-forge
tzdata                    2023.3                   pypi_0    pypi
wheel                     0.41.2             pyhd8ed1ab_0    conda-forge
xarray                    2023.8.1.dev7+g1043a9e1          pypi_0    pypi
xz                        5.2.6                h166bdaf_0    conda-forge
```

</details>


Discussion:
keewis: not sure what to do with the `float_` removal, but otherwise this should be ready for merging. That way, at least the `import` check in the upstream-dev CI succeeds.
kmuehlbauer: > Isn't `np.float_` and `float` equivalent

Not for numpy 2.0, where is has been removed: https://numpy.org/devdocs/release/2.0.0-notes.html#numpy-2-0-python-api-removals.

Maybe just replacing `np.float_` with `np.float64` as suggested in the error message will fix it?
mathause: Right, thanks. I found no deprecation warning or error in our upstream-dev CI, though. Maybe that's still to come? Anyway, given the description that seems to be the right approach.
dcherian: Yeah we'll need to do something about the. mypy errors:
```
xarray/core/dataset.py:30: error: Cannot find implementation or library stub for module named "numpy.exceptions"  [import]
xarray/core/dataset.py:30: note: See https://mypy.readthedocs.io/en/stable/running_mypy.html#missing-imports
xarray/tests/test_dataset.py:20: error: Cannot find implementation or library stub for module named "numpy.exceptions"  [import]
xarray/tests/test_dataarray.py:18: error: Cannot find implementation or library stub for module named "numpy.exceptions"  [import]
```
keewis: > Maybe just replacing `np.float_` with `np.float64` as suggested in the error message will fix it?

We can't, as we cannot know the default floating-point size in advance: as far as I remember, 32bit OSs (windows in particular) still has `float32` as the default `float`. Instead, just using the python `float` as dtype might work?

> Yeah we'll need to do something about the mypy errors:

I'd be tempted to just add `numpy.exceptions.*` to the list of ignored modules in `pyproject.toml`, but `numpy` does have a `.pyi` file for that, so not sure how to deal with the transition.
keewis: this should be ready for merging now. Note that there are lots of failing upstream-dev mypy errors and `dask` seems to have some internal inconsistency in its `rolling` implementation, but these can be fixed / discussed in separate issues / PRs.
mathause: numpy suggests replacing `np.float_` with `np.float64`: https://numpy.org/devdocs/release/2.0.0-notes.html

Patch:
diff --git a/pyproject.toml b/pyproject.toml
index 4d63fd56..15c6bf19 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -64,6 +64,7 @@ module = [
   "sparse.*",
   "toolz.*",
   "zarr.*",
+  "numpy.exceptions.*",  # remove once support for `numpy<2.0` has been dropped
 ]
 
 [[tool.mypy.overrides]]
diff --git a/xarray/backends/netCDF4_.py b/xarray/backends/netCDF4_.py
index e3f97d7b..7cbfa5b5 100644
--- a/xarray/backends/netCDF4_.py
+++ b/xarray/backends/netCDF4_.py
@@ -207,7 +207,7 @@ def _ensure_fill_value_valid(data, attributes):
     # work around for netCDF4/scipy issue where _FillValue has the wrong type:
     # https://github.com/Unidata/netcdf4-python/issues/271
     if data.dtype.kind == "S" and "_FillValue" in attributes:
-        attributes["_FillValue"] = np.string_(attributes["_FillValue"])
+        attributes["_FillValue"] = np.bytes_(attributes["_FillValue"])
 
 
 def _force_native_endianness(var):
diff --git a/xarray/coding/strings.py b/xarray/coding/strings.py
index d10af65c..89ceaddd 100644
--- a/xarray/coding/strings.py
+++ b/xarray/coding/strings.py
@@ -100,7 +100,7 @@ def ensure_fixed_length_bytes(var):
     dims, data, attrs, encoding = unpack_for_encoding(var)
     if check_vlen_dtype(data.dtype) == bytes:
         # TODO: figure out how to handle this with dask
-        data = np.asarray(data, dtype=np.string_)
+        data = np.asarray(data, dtype=np.bytes_)
     return Variable(dims, data, attrs, encoding)
 
 
@@ -151,7 +151,7 @@ def bytes_to_char(arr):
 def _numpy_bytes_to_char(arr):
     """Like netCDF4.stringtochar, but faster and more flexible."""
     # ensure the array is contiguous
-    arr = np.array(arr, copy=False, order="C", dtype=np.string_)
+    arr = np.array(arr, copy=False, order="C", dtype=np.bytes_)
     return arr.reshape(arr.shape + (1,)).view("S1")
 
 
@@ -168,7 +168,7 @@ def char_to_bytes(arr):
 
     if not size:
         # can't make an S0 dtype
-        return np.zeros(arr.shape[:-1], dtype=np.string_)
+        return np.zeros(arr.shape[:-1], dtype=np.bytes_)
 
     if is_chunked_array(arr):
         chunkmanager = get_chunked_array_type(arr)
diff --git a/xarray/coding/times.py b/xarray/coding/times.py
index 4291d959..b531dc97 100644
--- a/xarray/coding/times.py
+++ b/xarray/coding/times.py
@@ -467,7 +467,7 @@ def convert_times(times, date_type, raise_on_invalid: bool = True) -> np.ndarray
     Useful to convert between calendars in numpy and cftime or between cftime calendars.
 
     If raise_on_valid is True (default), invalid dates trigger a ValueError.
-    Otherwise, the invalid element is replaced by np.NaN for cftime types and np.NaT for np.datetime64.
+    Otherwise, the invalid element is replaced by np.nan for cftime types and np.NaT for np.datetime64.
     """
     if date_type in (pd.Timestamp, np.datetime64) and not is_np_datetime_like(
         times.dtype
@@ -489,7 +489,7 @@ def convert_times(times, date_type, raise_on_invalid: bool = True) -> np.ndarray
                     f"{date_type(2000, 1, 1).calendar} calendar.  Reason: {e}."
                 )
             else:
-                dt = np.NaN
+                dt = np.nan
 
         new[i] = dt
     return new
diff --git a/xarray/core/accessor_str.py b/xarray/core/accessor_str.py
index 31028f10..aa6dc2c7 100644
--- a/xarray/core/accessor_str.py
+++ b/xarray/core/accessor_str.py
@@ -471,7 +471,7 @@ class StringAccessor(Generic[T_DataArray]):
         ... )
         >>> values_2 = np.array(3.4)
         >>> values_3 = ""
-        >>> values_4 = np.array("test", dtype=np.unicode_)
+        >>> values_4 = np.array("test", dtype=np.str_)
 
         Determine the separator to use
 
diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py
index 57757179..fd3ff60c 100644
--- a/xarray/core/dataarray.py
+++ b/xarray/core/dataarray.py
@@ -5772,8 +5772,8 @@ class DataArray(
         >>> array = xr.DataArray(
         ...     [
         ...         [2.0, 1.0, 2.0, 0.0, -2.0],
-        ...         [-4.0, np.NaN, 2.0, np.NaN, -2.0],
-        ...         [np.NaN, np.NaN, 1.0, np.NaN, np.NaN],
+        ...         [-4.0, np.nan, 2.0, np.nan, -2.0],
+        ...         [np.nan, np.nan, 1.0, np.nan, np.nan],
         ...     ],
         ...     dims=["y", "x"],
         ...     coords={"y": [-1, 0, 1], "x": np.arange(5.0) ** 2},
@@ -5868,8 +5868,8 @@ class DataArray(
         >>> array = xr.DataArray(
         ...     [
         ...         [2.0, 1.0, 2.0, 0.0, -2.0],
-        ...         [-4.0, np.NaN, 2.0, np.NaN, -2.0],
-        ...         [np.NaN, np.NaN, 1.0, np.NaN, np.NaN],
+        ...         [-4.0, np.nan, 2.0, np.nan, -2.0],
+        ...         [np.nan, np.nan, 1.0, np.nan, np.nan],
         ...     ],
         ...     dims=["y", "x"],
         ...     coords={"y": [-1, 0, 1], "x": np.arange(5.0) ** 2},
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index 4e5ca374..4e1fe621 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -24,6 +24,13 @@ from os import PathLike
 from typing import IO, TYPE_CHECKING, Any, Callable, Generic, Literal, cast, overload
 
 import numpy as np
+
+# remove once numpy 2.0 is the oldest supported version
+try:
+    from numpy.exceptions import RankWarning
+except ImportError:
+    from numpy import RankWarning
+
 import pandas as pd
 
 from xarray.coding.calendar_ops import convert_calendar, interp_calendar
@@ -79,10 +86,7 @@ from xarray.core.merge import (
 )
 from xarray.core.missing import get_clean_interp_index
 from xarray.core.options import OPTIONS, _get_keep_attrs
-from xarray.core.parallelcompat import (
-    get_chunked_array_type,
-    guess_chunkmanager,
-)
+from xarray.core.parallelcompat import get_chunked_array_type, guess_chunkmanager
 from xarray.core.pycompat import (
     array_type,
     is_chunked_array,
@@ -5275,34 +5279,31 @@ class Dataset(
 
         stacking_dims = tuple(dim for dim in self.dims if dim not in sample_dims)
 
-        for variable in self:
-            dims = self[variable].dims
-            dims_include_sample_dims = set(sample_dims) <= set(dims)
-            if not dims_include_sample_dims:
+        for key, da in self.data_vars.items():
+            missing_sample_dims = set(sample_dims) - set(da.dims)
+            if missing_sample_dims:
                 raise ValueError(
-                    "All variables in the dataset must contain the "
-                    f"dimensions {dims}."
+                    "Variables in the dataset must contain all ``sample_dims`` "
+                    f"({sample_dims!r}) but '{key}' misses {sorted(missing_sample_dims)}"
                 )
 
-        def ensure_stackable(val):
-            assign_coords = {variable_dim: val.name}
-            for dim in stacking_dims:
-                if dim not in val.dims:
-                    assign_coords[dim] = None
+        def stack_dataarray(da):
+            # add missing dims/ coords and the name of the variable
+
+            missing_stack_coords = {variable_dim: da.name}
+            for dim in set(stacking_dims) - set(da.dims):
+                missing_stack_coords[dim] = None
 
-            expand_dims = set(stacking_dims).difference(set(val.dims))
-            expand_dims.add(variable_dim)
-            # must be list for .expand_dims
-            expand_dims = list(expand_dims)
+            missing_stack_dims = list(missing_stack_coords)
 
             return (
-                val.assign_coords(**assign_coords)
-                .expand_dims(expand_dims)
+                da.assign_coords(**missing_stack_coords)
+                .expand_dims(missing_stack_dims)
                 .stack({new_dim: (variable_dim,) + stacking_dims})
             )
 
         # concatenate the arrays
-        stackable_vars = [ensure_stackable(self[key]) for key in self.data_vars]
+        stackable_vars = [stack_dataarray(da) for da in self.data_vars.values()]
         data_array = concat(stackable_vars, dim=new_dim)
 
         if name is not None:
@@ -8791,9 +8792,9 @@ class Dataset(
 
             with warnings.catch_warnings():
                 if full:  # Copy np.polyfit behavior
-                    warnings.simplefilter("ignore", np.RankWarning)
+                    warnings.simplefilter("ignore", RankWarning)
                 else:  # Raise only once per variable
-                    warnings.simplefilter("once", np.RankWarning)
+                    warnings.simplefilter("once", RankWarning)
 
                 coeffs, residuals = duck_array_ops.least_squares(
                     lhs, rhs.data, rcond=rcond, skipna=skipna_da
@@ -9083,8 +9084,8 @@ class Dataset(
         >>> array2 = xr.DataArray(
         ...     [
         ...         [2.0, 1.0, 2.0, 0.0, -2.0],
-        ...         [-4.0, np.NaN, 2.0, np.NaN, -2.0],
-        ...         [np.NaN, np.NaN, 1.0, np.NaN, np.NaN],
+        ...         [-4.0, np.nan, 2.0, np.nan, -2.0],
+        ...         [np.nan, np.nan, 1.0, np.nan, np.nan],
         ...     ],
         ...     dims=["y", "x"],
         ...     coords={"y": [-1, 0, 1], "x": ["a", "b", "c", "d", "e"]},
@@ -9180,8 +9181,8 @@ class Dataset(
         >>> array2 = xr.DataArray(
         ...     [
         ...         [2.0, 1.0, 2.0, 0.0, -2.0],
-        ...         [-4.0, np.NaN, 2.0, np.NaN, -2.0],
-        ...         [np.NaN, np.NaN, 1.0, np.NaN, np.NaN],
+        ...         [-4.0, np.nan, 2.0, np.nan, -2.0],
+        ...         [np.nan, np.nan, 1.0, np.nan, np.nan],
         ...     ],
         ...     dims=["y", "x"],
         ...     coords={"y": [-1, 0, 1], "x": ["a", "b", "c", "d", "e"]},
diff --git a/xarray/core/dtypes.py b/xarray/core/dtypes.py
index 7ac342e3..0762fa03 100644
--- a/xarray/core/dtypes.py
+++ b/xarray/core/dtypes.py
@@ -40,7 +40,7 @@ NINF = AlwaysLessThan()
 PROMOTE_TO_OBJECT: tuple[tuple[type[np.generic], type[np.generic]], ...] = (
     (np.number, np.character),  # numpy promotes to character
     (np.bool_, np.character),  # numpy promotes to character
-    (np.bytes_, np.unicode_),  # numpy promotes to unicode
+    (np.bytes_, np.str_),  # numpy promotes to unicode
 )
 
 
diff --git a/xarray/core/missing.py b/xarray/core/missing.py
index c6efaebc..137d689d 100644
--- a/xarray/core/missing.py
+++ b/xarray/core/missing.py
@@ -730,7 +730,7 @@ def interp_func(var, x, new_x, method: InterpOptions, kwargs):
         # scipy.interpolate.interp1d always forces to float.
         # Use the same check for blockwise as well:
         if not issubclass(var.dtype.type, np.inexact):
-            dtype = np.float_
+            dtype = float
         else:
             dtype = var.dtype
 
diff --git a/xarray/core/nputils.py b/xarray/core/nputils.py
index 1c5b0d3d..9efa5824 100644
--- a/xarray/core/nputils.py
+++ b/xarray/core/nputils.py
@@ -6,6 +6,12 @@ import numpy as np
 import pandas as pd
 from numpy.core.multiarray import normalize_axis_index  # type: ignore[attr-defined]
 
+# remove once numpy 2.0 is the oldest supported version
+try:
+    from numpy.exceptions import RankWarning
+except ImportError:
+    from numpy import RankWarning
+
 from xarray.core.options import OPTIONS
 from xarray.core.pycompat import is_duck_array
 
@@ -194,7 +200,7 @@ def _nanpolyfit_1d(arr, x, rcond=None):
 
 def warn_on_deficient_rank(rank, order):
     if rank != order:
-        warnings.warn("Polyfit may be poorly conditioned", np.RankWarning, stacklevel=2)
+        warnings.warn("Polyfit may be poorly conditioned", RankWarning, stacklevel=2)
 
 
 def least_squares(lhs, rhs, rcond=None, skipna=False):


Test patch:
diff --git a/xarray/tests/test_accessor_str.py b/xarray/tests/test_accessor_str.py
index 168d3232..dc325a84 100644
--- a/xarray/tests/test_accessor_str.py
+++ b/xarray/tests/test_accessor_str.py
@@ -279,22 +279,20 @@ def test_case_bytes() -> None:
 def test_case_str() -> None:
     # This string includes some unicode characters
     # that are common case management corner cases
-    value = xr.DataArray(["SOme wOrd Ǆ ß ᾛ ΣΣ ﬃ⁵Å Ç Ⅰ"]).astype(np.unicode_)
-
-    exp_capitalized = xr.DataArray(["Some word ǆ ß ᾓ σς ﬃ⁵å ç ⅰ"]).astype(np.unicode_)
-    exp_lowered = xr.DataArray(["some word ǆ ß ᾓ σς ﬃ⁵å ç ⅰ"]).astype(np.unicode_)
-    exp_swapped = xr.DataArray(["soME WoRD ǆ SS ᾛ σς FFI⁵å ç ⅰ"]).astype(np.unicode_)
-    exp_titled = xr.DataArray(["Some Word ǅ Ss ᾛ Σς Ffi⁵Å Ç Ⅰ"]).astype(np.unicode_)
-    exp_uppered = xr.DataArray(["SOME WORD Ǆ SS ἫΙ ΣΣ FFI⁵Å Ç Ⅰ"]).astype(np.unicode_)
-    exp_casefolded = xr.DataArray(["some word ǆ ss ἣι σσ ffi⁵å ç ⅰ"]).astype(
-        np.unicode_
-    )
-
-    exp_norm_nfc = xr.DataArray(["SOme wOrd Ǆ ß ᾛ ΣΣ ﬃ⁵Å Ç Ⅰ"]).astype(np.unicode_)
-    exp_norm_nfkc = xr.DataArray(["SOme wOrd DŽ ß ᾛ ΣΣ ffi5Å Ç I"]).astype(np.unicode_)
-    exp_norm_nfd = xr.DataArray(["SOme wOrd Ǆ ß ᾛ ΣΣ ﬃ⁵Å Ç Ⅰ"]).astype(np.unicode_)
+    value = xr.DataArray(["SOme wOrd Ǆ ß ᾛ ΣΣ ﬃ⁵Å Ç Ⅰ"]).astype(np.str_)
+
+    exp_capitalized = xr.DataArray(["Some word ǆ ß ᾓ σς ﬃ⁵å ç ⅰ"]).astype(np.str_)
+    exp_lowered = xr.DataArray(["some word ǆ ß ᾓ σς ﬃ⁵å ç ⅰ"]).astype(np.str_)
+    exp_swapped = xr.DataArray(["soME WoRD ǆ SS ᾛ σς FFI⁵å ç ⅰ"]).astype(np.str_)
+    exp_titled = xr.DataArray(["Some Word ǅ Ss ᾛ Σς Ffi⁵Å Ç Ⅰ"]).astype(np.str_)
+    exp_uppered = xr.DataArray(["SOME WORD Ǆ SS ἫΙ ΣΣ FFI⁵Å Ç Ⅰ"]).astype(np.str_)
+    exp_casefolded = xr.DataArray(["some word ǆ ss ἣι σσ ffi⁵å ç ⅰ"]).astype(np.str_)
+
+    exp_norm_nfc = xr.DataArray(["SOme wOrd Ǆ ß ᾛ ΣΣ ﬃ⁵Å Ç Ⅰ"]).astype(np.str_)
+    exp_norm_nfkc = xr.DataArray(["SOme wOrd DŽ ß ᾛ ΣΣ ffi5Å Ç I"]).astype(np.str_)
+    exp_norm_nfd = xr.DataArray(["SOme wOrd Ǆ ß ᾛ ΣΣ ﬃ⁵Å Ç Ⅰ"]).astype(np.str_)
     exp_norm_nfkd = xr.DataArray(["SOme wOrd DŽ ß ᾛ ΣΣ ffi5Å Ç I"]).astype(
-        np.unicode_
+        np.str_
     )
 
     res_capitalized = value.str.capitalize()
@@ -680,7 +678,7 @@ def test_extract_extractall_name_collision_raises(dtype) -> None:
 def test_extract_single_case(dtype) -> None:
     pat_str = r"(\w+)_Xy_\d*"
     pat_re: str | bytes = (
-        pat_str if dtype == np.unicode_ else bytes(pat_str, encoding="UTF-8")
+        pat_str if dtype == np.str_ else bytes(pat_str, encoding="UTF-8")
     )
     pat_compiled = re.compile(pat_re)
 
@@ -728,7 +726,7 @@ def test_extract_single_case(dtype) -> None:
 def test_extract_single_nocase(dtype) -> None:
     pat_str = r"(\w+)?_Xy_\d*"
     pat_re: str | bytes = (
-        pat_str if dtype == np.unicode_ else bytes(pat_str, encoding="UTF-8")
+        pat_str if dtype == np.str_ else bytes(pat_str, encoding="UTF-8")
     )
     pat_compiled = re.compile(pat_re, flags=re.IGNORECASE)
 
@@ -770,7 +768,7 @@ def test_extract_single_nocase(dtype) -> None:
 def test_extract_multi_case(dtype) -> None:
     pat_str = r"(\w+)_Xy_(\d*)"
     pat_re: str | bytes = (
-        pat_str if dtype == np.unicode_ else bytes(pat_str, encoding="UTF-8")
+        pat_str if dtype == np.str_ else bytes(pat_str, encoding="UTF-8")
     )
     pat_compiled = re.compile(pat_re)
 
@@ -810,7 +808,7 @@ def test_extract_multi_case(dtype) -> None:
 def test_extract_multi_nocase(dtype) -> None:
     pat_str = r"(\w+)_Xy_(\d*)"
     pat_re: str | bytes = (
-        pat_str if dtype == np.unicode_ else bytes(pat_str, encoding="UTF-8")
+        pat_str if dtype == np.str_ else bytes(pat_str, encoding="UTF-8")
     )
     pat_compiled = re.compile(pat_re, flags=re.IGNORECASE)
 
@@ -876,7 +874,7 @@ def test_extract_broadcast(dtype) -> None:
 def test_extractall_single_single_case(dtype) -> None:
     pat_str = r"(\w+)_Xy_\d*"
     pat_re: str | bytes = (
-        pat_str if dtype == np.unicode_ else bytes(pat_str, encoding="UTF-8")
+        pat_str if dtype == np.str_ else bytes(pat_str, encoding="UTF-8")
     )
     pat_compiled = re.compile(pat_re)
 
@@ -908,7 +906,7 @@ def test_extractall_single_single_case(dtype) -> None:
 def test_extractall_single_single_nocase(dtype) -> None:
     pat_str = r"(\w+)_Xy_\d*"
     pat_re: str | bytes = (
-        pat_str if dtype == np.unicode_ else bytes(pat_str, encoding="UTF-8")
+        pat_str if dtype == np.str_ else bytes(pat_str, encoding="UTF-8")
     )
     pat_compiled = re.compile(pat_re, flags=re.I)
 
@@ -937,7 +935,7 @@ def test_extractall_single_single_nocase(dtype) -> None:
 def test_extractall_single_multi_case(dtype) -> None:
     pat_str = r"(\w+)_Xy_\d*"
     pat_re: str | bytes = (
-        pat_str if dtype == np.unicode_ else bytes(pat_str, encoding="UTF-8")
+        pat_str if dtype == np.str_ else bytes(pat_str, encoding="UTF-8")
     )
     pat_compiled = re.compile(pat_re)
 
@@ -983,7 +981,7 @@ def test_extractall_single_multi_case(dtype) -> None:
 def test_extractall_single_multi_nocase(dtype) -> None:
     pat_str = r"(\w+)_Xy_\d*"
     pat_re: str | bytes = (
-        pat_str if dtype == np.unicode_ else bytes(pat_str, encoding="UTF-8")
+        pat_str if dtype == np.str_ else bytes(pat_str, encoding="UTF-8")
     )
     pat_compiled = re.compile(pat_re, flags=re.I)
 
@@ -1030,7 +1028,7 @@ def test_extractall_single_multi_nocase(dtype) -> None:
 def test_extractall_multi_single_case(dtype) -> None:
     pat_str = r"(\w+)_Xy_(\d*)"
     pat_re: str | bytes = (
-        pat_str if dtype == np.unicode_ else bytes(pat_str, encoding="UTF-8")
+        pat_str if dtype == np.str_ else bytes(pat_str, encoding="UTF-8")
     )
     pat_compiled = re.compile(pat_re)
 
@@ -1065,7 +1063,7 @@ def test_extractall_multi_single_case(dtype) -> None:
 def test_extractall_multi_single_nocase(dtype) -> None:
     pat_str = r"(\w+)_Xy_(\d*)"
     pat_re: str | bytes = (
-        pat_str if dtype == np.unicode_ else bytes(pat_str, encoding="UTF-8")
+        pat_str if dtype == np.str_ else bytes(pat_str, encoding="UTF-8")
     )
     pat_compiled = re.compile(pat_re, flags=re.I)
 
@@ -1097,7 +1095,7 @@ def test_extractall_multi_single_nocase(dtype) -> None:
 def test_extractall_multi_multi_case(dtype) -> None:
     pat_str = r"(\w+)_Xy_(\d*)"
     pat_re: str | bytes = (
-        pat_str if dtype == np.unicode_ else bytes(pat_str, encoding="UTF-8")
+        pat_str if dtype == np.str_ else bytes(pat_str, encoding="UTF-8")
     )
     pat_compiled = re.compile(pat_re)
 
@@ -1147,7 +1145,7 @@ def test_extractall_multi_multi_case(dtype) -> None:
 def test_extractall_multi_multi_nocase(dtype) -> None:
     pat_str = r"(\w+)_Xy_(\d*)"
     pat_re: str | bytes = (
-        pat_str if dtype == np.unicode_ else bytes(pat_str, encoding="UTF-8")
+        pat_str if dtype == np.str_ else bytes(pat_str, encoding="UTF-8")
     )
     pat_compiled = re.compile(pat_re, flags=re.I)
 
@@ -3419,12 +3417,12 @@ def test_cat_multi() -> None:
 
     values_4 = ""
 
-    values_5 = np.array("", dtype=np.unicode_)
+    values_5 = np.array("", dtype=np.str_)
 
     sep = xr.DataArray(
         [" ", ", "],
         dims=["ZZ"],
-    ).astype(np.unicode_)
+    ).astype(np.str_)
 
     expected = xr.DataArray(
         [
@@ -3440,7 +3438,7 @@ def test_cat_multi() -> None:
             ],
         ],
         dims=["X", "Y", "ZZ"],
-    ).astype(np.unicode_)
+    ).astype(np.str_)
 
     res = values_1.str.cat(values_2, values_3, values_4, values_5, sep=sep)
 
@@ -3561,7 +3559,7 @@ def test_format_scalar() -> None:
     values = xr.DataArray(
         ["{}.{Y}.{ZZ}", "{},{},{X},{X}", "{X}-{Y}-{ZZ}"],
         dims=["X"],
-    ).astype(np.unicode_)
+    ).astype(np.str_)
 
     pos0 = 1
     pos1 = 1.2
@@ -3574,7 +3572,7 @@ def test_format_scalar() -> None:
     expected = xr.DataArray(
         ["1.X.None", "1,1.2,'test','test'", "'test'-X-None"],
         dims=["X"],
-    ).astype(np.unicode_)
+    ).astype(np.str_)
 
     res = values.str.format(pos0, pos1, pos2, X=X, Y=Y, ZZ=ZZ, W=W)
 
@@ -3586,7 +3584,7 @@ def test_format_broadcast() -> None:
     values = xr.DataArray(
         ["{}.{Y}.{ZZ}", "{},{},{X},{X}", "{X}-{Y}-{ZZ}"],
         dims=["X"],
-    ).astype(np.unicode_)
+    ).astype(np.str_)
 
     pos0 = 1
     pos1 = 1.2
@@ -3608,7 +3606,7 @@ def test_format_broadcast() -> None:
             ["'test'-X-None", "'test'-X-None"],
         ],
         dims=["X", "YY"],
-    ).astype(np.unicode_)
+    ).astype(np.str_)
 
     res = values.str.format(pos0, pos1, pos2, X=X, Y=Y, ZZ=ZZ, W=W)
 
@@ -3620,7 +3618,7 @@ def test_mod_scalar() -> None:
     values = xr.DataArray(
         ["%s.%s.%s", "%s,%s,%s", "%s-%s-%s"],
         dims=["X"],
-    ).astype(np.unicode_)
+    ).astype(np.str_)
 
     pos0 = 1
     pos1 = 1.2
@@ -3629,7 +3627,7 @@ def test_mod_scalar() -> None:
     expected = xr.DataArray(
         ["1.1.2.2.3", "1,1.2,2.3", "1-1.2-2.3"],
         dims=["X"],
-    ).astype(np.unicode_)
+    ).astype(np.str_)
 
     res = values.str % (pos0, pos1, pos2)
 
@@ -3641,7 +3639,7 @@ def test_mod_dict() -> None:
     values = xr.DataArray(
         ["%(a)s.%(a)s.%(b)s", "%(b)s,%(c)s,%(b)s", "%(c)s-%(b)s-%(a)s"],
         dims=["X"],
-    ).astype(np.unicode_)
+    ).astype(np.str_)
 
     a = 1
     b = 1.2
@@ -3650,7 +3648,7 @@ def test_mod_dict() -> None:
     expected = xr.DataArray(
         ["1.1.1.2", "1.2,2.3,1.2", "2.3-1.2-1"],
         dims=["X"],
-    ).astype(np.unicode_)
+    ).astype(np.str_)
 
     res = values.str % {"a": a, "b": b, "c": c}
 
@@ -3662,7 +3660,7 @@ def test_mod_broadcast_single() -> None:
     values = xr.DataArray(
         ["%s_1", "%s_2", "%s_3"],
         dims=["X"],
-    ).astype(np.unicode_)
+    ).astype(np.str_)
 
     pos = xr.DataArray(
         ["2.3", "3.44444"],
@@ -3672,7 +3670,7 @@ def test_mod_broadcast_single() -> None:
     expected = xr.DataArray(
         [["2.3_1", "3.44444_1"], ["2.3_2", "3.44444_2"], ["2.3_3", "3.44444_3"]],
         dims=["X", "YY"],
-    ).astype(np.unicode_)
+    ).astype(np.str_)
 
     res = values.str % pos
 
@@ -3684,7 +3682,7 @@ def test_mod_broadcast_multi() -> None:
     values = xr.DataArray(
         ["%s.%s.%s", "%s,%s,%s", "%s-%s-%s"],
         dims=["X"],
-    ).astype(np.unicode_)
+    ).astype(np.str_)
 
     pos0 = 1
     pos1 = 1.2
@@ -3701,7 +3699,7 @@ def test_mod_broadcast_multi() -> None:
             ["1-1.2-2.3", "1-1.2-3.44444"],
         ],
         dims=["X", "YY"],
-    ).astype(np.unicode_)
+    ).astype(np.str_)
 
     res = values.str % (pos0, pos1, pos2)
 
diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py
index d54e1004..e2ae34f9 100644
--- a/xarray/tests/test_backends.py
+++ b/xarray/tests/test_backends.py
@@ -814,7 +814,7 @@ class DatasetIOBase:
     def test_dropna(self) -> None:
         # regression test for GH:issue:1694
         a = np.random.randn(4, 3)
-        a[1, 1] = np.NaN
+        a[1, 1] = np.nan
         in_memory = xr.Dataset(
             {"a": (("y", "x"), a)}, coords={"y": np.arange(4), "x": np.arange(3)}
         )
diff --git a/xarray/tests/test_dataarray.py b/xarray/tests/test_dataarray.py
index 2a28939d..1f4d259d 100644
--- a/xarray/tests/test_dataarray.py
+++ b/xarray/tests/test_dataarray.py
@@ -14,6 +14,12 @@ import pandas as pd
 import pytest
 from packaging.version import Version
 
+# remove once numpy 2.0 is the oldest supported version
+try:
+    from numpy.exceptions import RankWarning
+except ImportError:
+    from numpy import RankWarning
+
 import xarray as xr
 from xarray import (
     DataArray,
@@ -2846,7 +2852,7 @@ class TestDataArray:
     )
     def test_quantile(self, q, axis, dim, skipna) -> None:
         va = self.va.copy(deep=True)
-        va[0, 0] = np.NaN
+        va[0, 0] = np.nan
 
         actual = DataArray(va).quantile(q, dim=dim, keep_attrs=True, skipna=skipna)
         _percentile_func = np.nanpercentile if skipna in (True, None) else np.percentile
@@ -3124,10 +3130,10 @@ class TestDataArray:
         b = DataArray([1, 2], dims=["x"], coords={"x": ["b", "c"]})
 
         expected_a = DataArray(
-            [0, 1, np.NaN], dims=["x"], coords={"x": ["a", "b", "c"]}
+            [0, 1, np.nan], dims=["x"], coords={"x": ["a", "b", "c"]}
         )
         expected_b = DataArray(
-            [np.NaN, 1, 2], dims=["x"], coords={"x": ["a", "b", "c"]}
+            [np.nan, 1, 2], dims=["x"], coords={"x": ["a", "b", "c"]}
         )
 
         actual_a, actual_b = xr.align(a, b, join="outer")
@@ -4211,7 +4217,7 @@ class TestDataArray:
 
         # Full output and deficient rank
         with warnings.catch_warnings():
-            warnings.simplefilter("ignore", np.RankWarning)
+            warnings.simplefilter("ignore", RankWarning)
             out = da.polyfit("x", 12, full=True)
             assert out.polyfit_residuals.isnull().all()
 
@@ -4232,7 +4238,7 @@ class TestDataArray:
         np.testing.assert_almost_equal(out.polyfit_residuals, [0, 0])
 
         with warnings.catch_warnings():
-            warnings.simplefilter("ignore", np.RankWarning)
+            warnings.simplefilter("ignore", RankWarning)
             out = da.polyfit("x", 8, full=True)
             np.testing.assert_array_equal(out.polyfit_residuals.isnull(), [True, False])
 
@@ -4253,7 +4259,7 @@ class TestDataArray:
         ar = xr.DataArray([9], dims="x")
 
         actual = ar.pad(x=1)
-        expected = xr.DataArray([np.NaN, 9, np.NaN], dims="x")
+        expected = xr.DataArray([np.nan, 9, np.nan], dims="x")
         assert_identical(actual, expected)
 
         actual = ar.pad(x=1, constant_values=1.23456)
@@ -4261,7 +4267,7 @@ class TestDataArray:
         assert_identical(actual, expected)
 
         with pytest.raises(ValueError, match="cannot convert float NaN to integer"):
-            ar.pad(x=1, constant_values=np.NaN)
+            ar.pad(x=1, constant_values=np.nan)
 
     def test_pad_coords(self) -> None:
         ar = DataArray(
@@ -4699,10 +4705,10 @@ class TestReduce:
             np.array([0.0, 1.0, 2.0, 0.0, -2.0, -4.0, 2.0]), 5, 2, None, id="float"
         ),
         pytest.param(
-            np.array([1.0, np.NaN, 2.0, np.NaN, -2.0, -4.0, 2.0]), 5, 2, 1, id="nan"
+            np.array([1.0, np.nan, 2.0, np.nan, -2.0, -4.0, 2.0]), 5, 2, 1, id="nan"
         ),
         pytest.param(
-            np.array([1.0, np.NaN, 2.0, np.NaN, -2.0, -4.0, 2.0]).astype("object"),
+            np.array([1.0, np.nan, 2.0, np.nan, -2.0, -4.0, 2.0]).astype("object"),
             5,
             2,
             1,
@@ -4711,7 +4717,7 @@ class TestReduce:
             ),
             id="obj",
         ),
-        pytest.param(np.array([np.NaN, np.NaN]), np.NaN, np.NaN, 0, id="allnan"),
+        pytest.param(np.array([np.nan, np.nan]), np.nan, np.nan, 0, id="allnan"),
         pytest.param(
             np.array(
                 ["2015-12-31", "2020-01-02", "2020-01-01", "2016-01-01"],
@@ -4906,7 +4912,7 @@ class TestReduce1D(TestReduce):
 
         if hasna:
             coordarr1[...] = 1
-            fill_value_0 = np.NaN
+            fill_value_0 = np.nan
         else:
             fill_value_0 = 1
 
@@ -4920,7 +4926,7 @@ class TestReduce1D(TestReduce):
         assert_identical(result0, expected0)
 
         # Manually specify NaN fill_value
-        result1 = ar0.idxmin(fill_value=np.NaN)
+        result1 = ar0.idxmin(fill_value=np.nan)
         assert_identical(result1, expected0)
 
         # keep_attrs
@@ -5021,7 +5027,7 @@ class TestReduce1D(TestReduce):
 
         if hasna:
             coordarr1[...] = 1
-            fill_value_0 = np.NaN
+            fill_value_0 = np.nan
         else:
             fill_value_0 = 1
 
@@ -5035,7 +5041,7 @@ class TestReduce1D(TestReduce):
         assert_identical(result0, expected0)
 
         # Manually specify NaN fill_value
-        result1 = ar0.idxmax(fill_value=np.NaN)
+        result1 = ar0.idxmax(fill_value=np.nan)
         assert_identical(result1, expected0)
 
         # keep_attrs
@@ -5200,12 +5206,12 @@ class TestReduce1D(TestReduce):
             np.array(
                 [
                     [2.0, 1.0, 2.0, 0.0, -2.0, -4.0, 2.0],
-                    [-4.0, np.NaN, 2.0, np.NaN, -2.0, -4.0, 2.0],
-                    [np.NaN] * 7,
+                    [-4.0, np.nan, 2.0, np.nan, -2.0, -4.0, 2.0],
+                    [np.nan] * 7,
                 ]
             ),
-            [5, 0, np.NaN],
-            [0, 2, np.NaN],
+            [5, 0, np.nan],
+            [0, 2, np.nan],
             [None, 1, 0],
             id="nan",
         ),
@@ -5213,12 +5219,12 @@ class TestReduce1D(TestReduce):
             np.array(
                 [
                     [2.0, 1.0, 2.0, 0.0, -2.0, -4.0, 2.0],
-                    [-4.0, np.NaN, 2.0, np.NaN, -2.0, -4.0, 2.0],
-                    [np.NaN] * 7,
+                    [-4.0, np.nan, 2.0, np.nan, -2.0, -4.0, 2.0],
+                    [np.nan] * 7,
                 ]
             ).astype("object"),
-            [5, 0, np.NaN],
-            [0, 2, np.NaN],
+            [5, 0, np.nan],
+            [0, 2, np.nan],
             [None, 1, 0],
             marks=pytest.mark.filterwarnings(
                 "ignore:invalid value encountered in reduce:RuntimeWarning:"
@@ -5493,7 +5499,7 @@ class TestReduce2D(TestReduce):
         coordarr1[hasna, :] = 1
         minindex0 = [x if not np.isnan(x) else 0 for x in minindex]
 
-        nan_mult_0 = np.array([np.NaN if x else 1 for x in hasna])[:, None]
+        nan_mult_0 = np.array([np.nan if x else 1 for x in hasna])[:, None]
         expected0list = [
             (coordarr1 * nan_mult_0).isel(y=yi).isel(x=indi, drop=True)
             for yi, indi in enumerate(minindex0)
@@ -5508,7 +5514,7 @@ class TestReduce2D(TestReduce):
 
         # Manually specify NaN fill_value
         with raise_if_dask_computes(max_computes=max_computes):
-            result1 = ar0.idxmin(dim="x", fill_value=np.NaN)
+            result1 = ar0.idxmin(dim="x", fill_value=np.nan)
         assert_identical(result1, expected0)
 
         # keep_attrs
@@ -5635,7 +5641,7 @@ class TestReduce2D(TestReduce):
         coordarr1[hasna, :] = 1
         maxindex0 = [x if not np.isnan(x) else 0 for x in maxindex]
 
-        nan_mult_0 = np.array([np.NaN if x else 1 for x in hasna])[:, None]
+        nan_mult_0 = np.array([np.nan if x else 1 for x in hasna])[:, None]
         expected0list = [
             (coordarr1 * nan_mult_0).isel(y=yi).isel(x=indi, drop=True)
             for yi, indi in enumerate(maxindex0)
@@ -5650,7 +5656,7 @@ class TestReduce2D(TestReduce):
 
         # Manually specify NaN fill_value
         with raise_if_dask_computes(max_computes=max_computes):
-            result1 = ar0.idxmax(dim="x", fill_value=np.NaN)
+            result1 = ar0.idxmax(dim="x", fill_value=np.nan)
         assert_identical(result1, expected0)
 
         # keep_attrs
@@ -5909,31 +5915,31 @@ class TestReduce2D(TestReduce):
             np.array(
                 [
                     [[2.0, 1.0, 2.0, 0.0], [-2.0, -4.0, 2.0, 0.0]],
-                    [[-4.0, np.NaN, 2.0, np.NaN], [-2.0, -4.0, 2.0, 0.0]],
-                    [[np.NaN] * 4, [np.NaN] * 4],
+                    [[-4.0, np.nan, 2.0, np.nan], [-2.0, -4.0, 2.0, 0.0]],
+                    [[np.nan] * 4, [np.nan] * 4],
                 ]
             ),
             {"x": np.array([[1, 0, 0, 0], [0, 0, 0, 0]])},
             {
                 "y": np.array(
-                    [[1, 1, 0, 0], [0, 1, 0, 1], [np.NaN, np.NaN, np.NaN, np.NaN]]
+                    [[1, 1, 0, 0], [0, 1, 0, 1], [np.nan, np.nan, np.nan, np.nan]]
                 )
             },
-            {"z": np.array([[3, 1], [0, 1], [np.NaN, np.NaN]])},
+            {"z": np.array([[3, 1], [0, 1], [np.nan, np.nan]])},
             {"x": np.array([1, 0, 0, 0]), "y": np.array([0, 1, 0, 0])},
             {"x": np.array([1, 0]), "z": np.array([0, 1])},
-            {"y": np.array([1, 0, np.NaN]), "z": np.array([1, 0, np.NaN])},
+            {"y": np.array([1, 0, np.nan]), "z": np.array([1, 0, np.nan])},
             {"x": np.array(0), "y": np.array(1), "z": np.array(1)},
             {"x": np.array([[0, 0, 0, 0], [0, 0, 0, 0]])},
             {
                 "y": np.array(
-                    [[0, 0, 0, 0], [1, 1, 0, 1], [np.NaN, np.NaN, np.NaN, np.NaN]]
+                    [[0, 0, 0, 0], [1, 1, 0, 1], [np.nan, np.nan, np.nan, np.nan]]
                 )
             },
-            {"z": np.array([[0, 2], [2, 2], [np.NaN, np.NaN]])},
+            {"z": np.array([[0, 2], [2, 2], [np.nan, np.nan]])},
             {"x": np.array([0, 0, 0, 0]), "y": np.array([0, 0, 0, 0])},
             {"x": np.array([0, 0]), "z": np.array([2, 2])},
-            {"y": np.array([0, 0, np.NaN]), "z": np.array([0, 2, np.NaN])},
+            {"y": np.array([0, 0, np.nan]), "z": np.array([0, 2, np.nan])},
             {"x": np.array(0), "y": np.array(0), "z": np.array(0)},
             {"x": np.array([[2, 1, 2, 1], [2, 2, 2, 2]])},
             {
@@ -5952,31 +5958,31 @@ class TestReduce2D(TestReduce):
             np.array(
                 [
                     [[2.0, 1.0, 2.0, 0.0], [-2.0, -4.0, 2.0, 0.0]],
-                    [[-4.0, np.NaN, 2.0, np.NaN], [-2.0, -4.0, 2.0, 0.0]],
-                    [[np.NaN] * 4, [np.NaN] * 4],
+                    [[-4.0, np.nan, 2.0, np.nan], [-2.0, -4.0, 2.0, 0.0]],
+                    [[np.nan] * 4, [np.nan] * 4],
                 ]
             ).astype("object"),
             {"x": np.array([[1, 0, 0, 0], [0, 0, 0, 0]])},
             {
                 "y": np.array(
-                    [[1, 1, 0, 0], [0, 1, 0, 1], [np.NaN, np.NaN, np.NaN, np.NaN]]
+                    [[1, 1, 0, 0], [0, 1, 0, 1], [np.nan, np.nan, np.nan, np.nan]]
                 )
             },
-            {"z": np.array([[3, 1], [0, 1], [np.NaN, np.NaN]])},
+            {"z": np.array([[3, 1], [0, 1], [np.nan, np.nan]])},
             {"x": np.array([1, 0, 0, 0]), "y": np.array([0, 1, 0, 0])},
             {"x": np.array([1, 0]), "z": np.array([0, 1])},
-            {"y": np.array([1, 0, np.NaN]), "z": np.array([1, 0, np.NaN])},
+            {"y": np.array([1, 0, np.nan]), "z": np.array([1, 0, np.nan])},
             {"x": np.array(0), "y": np.array(1), "z": np.array(1)},
             {"x": np.array([[0, 0, 0, 0], [0, 0, 0, 0]])},
             {
                 "y": np.array(
-                    [[0, 0, 0, 0], [1, 1, 0, 1], [np.NaN, np.NaN, np.NaN, np.NaN]]
+                    [[0, 0, 0, 0], [1, 1, 0, 1], [np.nan, np.nan, np.nan, np.nan]]
                 )
             },
-            {"z": np.array([[0, 2], [2, 2], [np.NaN, np.NaN]])},
+            {"z": np.array([[0, 2], [2, 2], [np.nan, np.nan]])},
             {"x": np.array([0, 0, 0, 0]), "y": np.array([0, 0, 0, 0])},
             {"x": np.array([0, 0]), "z": np.array([2, 2])},
-            {"y": np.array([0, 0, np.NaN]), "z": np.array([0, 2, np.NaN])},
+            {"y": np.array([0, 0, np.nan]), "z": np.array([0, 2, np.nan])},
             {"x": np.array(0), "y": np.array(0), "z": np.array(0)},
             {"x": np.array([[2, 1, 2, 1], [2, 2, 2, 2]])},
             {
@@ -6522,12 +6528,12 @@ def test_isin(da) -> None:
 
 def test_raise_no_warning_for_nan_in_binary_ops() -> None:
     with assert_no_warnings():
-        xr.DataArray([1, 2, np.NaN]) > 0
+        xr.DataArray([1, 2, np.nan]) > 0
 
 
 @pytest.mark.filterwarnings("error")
 def test_no_warning_for_all_nan() -> None:
-    _ = xr.DataArray([np.NaN, np.NaN]).mean()
+    _ = xr.DataArray([np.nan, np.nan]).mean()
 
 
 def test_name_in_masking() -> None:
@@ -6567,7 +6573,7 @@ class TestIrisConversion:
         )
 
         # Set a bad value to test the masking logic
-        original.data[0, 2] = np.NaN
+        original.data[0, 2] = np.nan
 
         original.attrs["cell_methods"] = "height: mean (comment: A cell method)"
         actual = original.to_iris()
diff --git a/xarray/tests/test_dataset.py b/xarray/tests/test_dataset.py
index 226e2b6d..882285ac 100644
--- a/xarray/tests/test_dataset.py
+++ b/xarray/tests/test_dataset.py
@@ -15,6 +15,12 @@ import pandas as pd
 import pytest
 from pandas.core.indexes.datetimes import DatetimeIndex
 
+# remove once numpy 2.0 is the oldest supported version
+try:
+    from numpy.exceptions import RankWarning
+except ImportError:
+    from numpy import RankWarning
+
 import xarray as xr
 from xarray import (
     DataArray,
@@ -118,7 +124,7 @@ def create_append_test_data(seed=None) -> tuple[Dataset, Dataset, Dataset]:
                 ),
                 "unicode_var": xr.DataArray(
                     unicode_var, coords=[time1], dims=["time"]
-                ).astype(np.unicode_),
+                ).astype(np.str_),
                 "datetime_var": xr.DataArray(
                     datetime_var, coords=[time1], dims=["time"]
                 ),
@@ -141,7 +147,7 @@ def create_append_test_data(seed=None) -> tuple[Dataset, Dataset, Dataset]:
                 ),
                 "unicode_var": xr.DataArray(
                     unicode_var[:nt2], coords=[time2], dims=["time"]
-                ).astype(np.unicode_),
+                ).astype(np.str_),
                 "datetime_var": xr.DataArray(
                     datetime_var_to_append, coords=[time2], dims=["time"]
                 ),
@@ -2432,10 +2438,10 @@ class TestDataset:
         b = Dataset({"foo": ("x", [1, 2])}, coords={"x": ["b", "c"]})
 
         expected_a = Dataset(
-            {"foo": ("x", [0, 1, np.NaN])}, coords={"x": ["a", "b", "c"]}
+            {"foo": ("x", [0, 1, np.nan])}, coords={"x": ["a", "b", "c"]}
         )
         expected_b = Dataset(
-            {"foo": ("x", [np.NaN, 1, 2])}, coords={"x": ["a", "b", "c"]}
+            {"foo": ("x", [np.nan, 1, 2])}, coords={"x": ["a", "b", "c"]}
         )
 
         actual_a, actual_b = xr.align(a, b, join="outer")
@@ -3792,7 +3798,10 @@ class TestDataset:
             data_vars={"a": (("x", "y"), [[0, 1, 2], [3, 4, 5]]), "b": ("x", [6, 7])},
             coords={"y": ["u", "v", "w"]},
         )
-        with pytest.raises(ValueError):
+        with pytest.raises(
+            ValueError,
+            match=r"Variables in the dataset must contain all ``sample_dims`` \(\['y'\]\) but 'b' misses \['y'\]",
+        ):
             data.to_stacked_array("features", sample_dims=["y"])
 
     def test_to_stacked_array_name(self) -> None:
@@ -5502,7 +5511,7 @@ class TestDataset:
     @pytest.mark.parametrize("q", [0.25, [0.50], [0.25, 0.75]])
     def test_quantile(self, q, skipna) -> None:
         ds = create_test_data(seed=123)
-        ds.var1.data[0, 0] = np.NaN
+        ds.var1.data[0, 0] = np.nan
 
         for dim in [None, "dim1", ["dim1"]]:
             ds_quantile = ds.quantile(q, dim=dim, skipna=skipna)
@@ -6375,7 +6384,7 @@ class TestDataset:
         with warnings.catch_warnings(record=True) as ws:
             ds.var1.polyfit("dim2", 10, full=False)
             assert len(ws) == 1
-            assert ws[0].category == np.RankWarning
+            assert ws[0].category == RankWarning
             ds.var1.polyfit("dim2", 10, full=True)
             assert len(ws) == 1
 
@@ -6702,7 +6711,7 @@ def test_dir_unicode(ds) -> None:
 
 def test_raise_no_warning_for_nan_in_binary_ops() -> None:
     with assert_no_warnings():
-        Dataset(data_vars={"x": ("y", [1, 2, np.NaN])}) > 0
+        Dataset(data_vars={"x": ("y", [1, 2, np.nan])}) > 0
 
 
 @pytest.mark.filterwarnings("error")
diff --git a/xarray/tests/test_dtypes.py b/xarray/tests/test_dtypes.py
index 490520c8..3c2ee5e8 100644
--- a/xarray/tests/test_dtypes.py
+++ b/xarray/tests/test_dtypes.py
@@ -10,12 +10,12 @@ from xarray.core import dtypes
     "args, expected",
     [
         ([bool], bool),
-        ([bool, np.string_], np.object_),
+        ([bool, np.bytes_], np.object_),
         ([np.float32, np.float64], np.float64),
-        ([np.float32, np.string_], np.object_),
-        ([np.unicode_, np.int64], np.object_),
-        ([np.unicode_, np.unicode_], np.unicode_),
-        ([np.bytes_, np.unicode_], np.object_),
+        ([np.float32, np.bytes_], np.object_),
+        ([np.str_, np.int64], np.object_),
+        ([np.str_, np.str_], np.str_),
+        ([np.bytes_, np.str_], np.object_),
     ],
 )
 def test_result_type(args, expected) -> None:
diff --git a/xarray/tests/test_groupby.py b/xarray/tests/test_groupby.py
index 5d99eda1..b961abef 100644
--- a/xarray/tests/test_groupby.py
+++ b/xarray/tests/test_groupby.py
@@ -232,11 +232,11 @@ def test_da_groupby_quantile() -> None:
     assert_identical(expected, actual)
 
     array = xr.DataArray(
-        data=[np.NaN, 2, 3, 4, 5, 6], coords={"x": [1, 1, 1, 2, 2, 2]}, dims="x"
+        data=[np.nan, 2, 3, 4, 5, 6], coords={"x": [1, 1, 1, 2, 2, 2]}, dims="x"
     )
 
     for skipna in (True, False, None):
-        e = [np.NaN, 5] if skipna is False else [2.5, 5]
+        e = [np.nan, 5] if skipna is False else [2.5, 5]
 
         expected = xr.DataArray(data=e, coords={"x": [1, 2], "quantile": 0.5}, dims="x")
         actual = array.groupby("x").quantile(0.5, skipna=skipna)
@@ -346,12 +346,12 @@ def test_ds_groupby_quantile() -> None:
     assert_identical(expected, actual)
 
     ds = xr.Dataset(
-        data_vars={"a": ("x", [np.NaN, 2, 3, 4, 5, 6])},
+        data_vars={"a": ("x", [np.nan, 2, 3, 4, 5, 6])},
         coords={"x": [1, 1, 1, 2, 2, 2]},
     )
 
     for skipna in (True, False, None):
-        e = [np.NaN, 5] if skipna is False else [2.5, 5]
+        e = [np.nan, 5] if skipna is False else [2.5, 5]
 
         expected = xr.Dataset(
             data_vars={"a": ("x", e)}, coords={"quantile": 0.5, "x": [1, 2]}
diff --git a/xarray/tests/test_variable.py b/xarray/tests/test_variable.py
index 2ef34201..c4bac3b2 100644
--- a/xarray/tests/test_variable.py
+++ b/xarray/tests/test_variable.py
@@ -196,7 +196,7 @@ class VariableSubclassobjects(ABC):
             self._assertIndexedLikeNDArray(x, value, dtype)
 
     def test_index_0d_float(self):
-        for value, dtype in [(0.5, np.float_), (np.float32(0.5), np.float32)]:
+        for value, dtype in [(0.5, float), (np.float32(0.5), np.float32)]:
             x = self.cls(["x"], [value])
             self._assertIndexedLikeNDArray(x, value, dtype)
 
@@ -1127,9 +1127,9 @@ class TestVariable(VariableSubclassobjects):
         assert v.dtype == np.dtype("U3")
         assert v.values == "foo"
 
-        v = Variable([], np.string_("foo"))
+        v = Variable([], np.bytes_("foo"))
         assert v.dtype == np.dtype("S3")
-        assert v.values == bytes("foo", "ascii")
+        assert v.values == "foo".encode("ascii")
 
     def test_0d_datetime(self):
         v = Variable([], pd.Timestamp("2000-01-01"))
@@ -1466,10 +1466,10 @@ class TestVariable(VariableSubclassobjects):
 
     def test_index_0d_numpy_string(self):
         # regression test to verify our work around for indexing 0d strings
-        v = Variable([], np.string_("asdf"))
+        v = Variable([], np.bytes_("asdf"))
         assert_identical(v[()], v)
 
-        v = Variable([], np.unicode_("asdf"))
+        v = Variable([], np.str_("asdf"))
         assert_identical(v[()], v)
 
     def test_indexing_0d_unicode(self):
@@ -1810,7 +1810,7 @@ class TestVariable(VariableSubclassobjects):
     )
     def test_quantile(self, q, axis, dim, skipna):
         d = self.d.copy()
-        d[0, 0] = np.NaN
+        d[0, 0] = np.nan
 
         v = Variable(["x", "y"], d)
         actual = v.quantile(q, dim=dim, skipna=skipna)
@@ -2719,7 +2719,7 @@ class TestAsCompatibleData:
 
 def test_raise_no_warning_for_nan_in_binary_ops():
     with assert_no_warnings():
-        Variable("x", [1, 2, np.NaN]) > 0
+        Variable("x", [1, 2, np.nan]) > 0
 
 
 class TestBackendIndexing:
diff --git a/xarray/tests/test_weighted.py b/xarray/tests/test_weighted.py
index 628d6310..95fda3fa 100644
--- a/xarray/tests/test_weighted.py
+++ b/xarray/tests/test_weighted.py
@@ -608,7 +608,7 @@ def test_weighted_operations_3D(dim, add_nans, skipna):
     # add approximately 25 % NaNs (https://stackoverflow.com/a/32182680/3010700)
     if add_nans:
         c = int(data.size * 0.25)
-        data.ravel()[np.random.choice(data.size, c, replace=False)] = np.NaN
+        data.ravel()[np.random.choice(data.size, c, replace=False)] = np.nan
 
     data = DataArray(data, dims=dims, coords=coords)
 
@@ -631,7 +631,7 @@ def test_weighted_quantile_3D(dim, q, add_nans, skipna):
     # add approximately 25 % NaNs (https://stackoverflow.com/a/32182680/3010700)
     if add_nans:
         c = int(data.size * 0.25)
-        data.ravel()[np.random.choice(data.size, c, replace=False)] = np.NaN
+        data.ravel()[np.random.choice(data.size, c, replace=False)] = np.nan
 
     da = DataArray(data, dims=dims, coords=coords)
 
@@ -709,7 +709,7 @@ def test_weighted_operations_different_shapes(
     # add approximately 25 % NaNs
     if add_nans:
         c = int(data.size * 0.25)
-        data.ravel()[np.random.choice(data.size, c, replace=False)] = np.NaN
+        data.ravel()[np.random.choice(data.size, c, replace=False)] = np.nan
 
     data = DataArray(data)