[Numpy2] Support for numpy==2.0.0 (and pin iris >=3.11) (#2395)

Co-authored-by: Manuel Schlund <[email protected]>
Co-authored-by: Manuel Schlund <[email protected]>

environment.yml

@@ -18,7 +18,7 @@ dependencies:
   - fire
   - geopy
   - humanfriendly
-  - iris >=3.10.0
+  - iris >=3.11  # 3.11 first to support Numpy 2 and Python 3.13
   - iris-esmf-regrid >=0.11.0
   - iris-grib >=0.20.0  # github.com/ESMValGroup/ESMValCore/issues/2535
   - isodate >=0.7.0  # incompatible with very old 0.6.1
@@ -27,7 +27,7 @@ dependencies:
   - nc-time-axis
   - nested-lookup
   - netcdf4
-  - numpy !=1.24.3,<2.0.0  # avoid pulling 2.0.0rcX
+  - numpy !=1.24.3
   - packaging
   - pandas
   - pillow

esmvalcore/preprocessor/_compare_with_refs.py

@@ -452,7 +452,11 @@ def _calculate_rmse(
     weights = get_weights(cube, coords) if weighted else None
     squared_error = (cube.core_data() - reference.core_data()) ** 2
     npx = get_array_module(squared_error)
-    rmse = npx.sqrt(npx.ma.average(squared_error, axis=axis, weights=weights))
+    mse = npx.ma.average(squared_error, axis=axis, weights=weights)
+    if isinstance(mse, da.Array):
+        rmse = da.reductions.safe_sqrt(mse)
+    else:
+        rmse = np.ma.sqrt(mse)
 
     # Metadata
     metadata = CubeMetadata(

esmvalcore/preprocessor/_derive/_shared.py

@@ -163,7 +163,7 @@ def _create_pressure_array(cube, ps_cube, top_limit):
     ps_4d_array = iris.util.broadcast_to_shape(ps_cube.data, shape, [0, 2, 3])
 
     # Set pressure levels below the surface pressure to NaN
-    pressure_4d = np.where((ps_4d_array - p_4d_array) < 0, np.NaN, p_4d_array)
+    pressure_4d = np.where((ps_4d_array - p_4d_array) < 0, np.nan, p_4d_array)
 
     # Make top_limit last pressure level
     top_limit_array = np.full(ps_cube.shape, top_limit, dtype=np.float32)

esmvalcore/preprocessor/_derive/co2s.py

@@ -17,6 +17,8 @@ def _get_first_unmasked_data(array, axis):
         *[da.arange(array.shape[i]) for i in range(array.ndim) if i != axis],
         indexing="ij",
     )
+
+    indices = list(indices)
     indices.insert(axis, indices_first_positive)
     first_unmasked_data = np.array(array)[tuple(indices)]
     return first_unmasked_data

esmvalcore/preprocessor/_regrid.py

@@ -209,7 +209,7 @@ def _generate_cube_from_dimcoords(latdata, londata, circular: bool = False):
 
     # Construct the resultant stock cube, with dummy data.
     shape = (latdata.size, londata.size)
-    dummy = np.empty(shape, dtype=np.dtype("int8"))
+    dummy = np.empty(shape, dtype=np.int32)
     coords_spec = [(lats, 0), (lons, 1)]
     cube = Cube(dummy, dim_coords_and_dims=coords_spec)
 

esmvalcore/preprocessor/_regrid_esmpy.py

@@ -283,8 +283,9 @@ def get_grid(
         num_peri_dims = 1
     else:
         num_peri_dims = 0
+
     grid = esmpy.Grid(
-        np.array(esmpy_lat.shape),
+        np.vstack(esmpy_lat.shape),
         num_peri_dims=num_peri_dims,
         staggerloc=[esmpy.StaggerLoc.CENTER],
     )

esmvalcore/preprocessor/_regrid_unstructured.py

@@ -170,8 +170,9 @@ def _get_weights_and_idx(
         src_points_with_convex_hull = self._add_convex_hull_twice(
             src_points, hull.vertices
         )
-        src_points_with_convex_hull[-2 * n_hull : -n_hull, 1] -= 360
-        src_points_with_convex_hull[-n_hull:, 1] += 360
+        lon_period = np.array(360, dtype=src_points_with_convex_hull.dtype)
+        src_points_with_convex_hull[-2 * n_hull : -n_hull, 1] -= lon_period
+        src_points_with_convex_hull[-n_hull:, 1] += lon_period
 
         # Actual weights calculation
         (weights, indices) = self._calculate_weights(

esmvalcore/preprocessor/_time.py

@@ -632,7 +632,9 @@ def seasonal_statistics(
             cube, "time", name="clim_season", seasons=seasons
         )
     else:
-        old_seasons = sorted(set(cube.coord("clim_season").points))
+        old_seasons = sorted(
+            {str(s) for s in cube.coord("clim_season").points}
+        )
         if not all(osea in seasons for osea in old_seasons):
             raise ValueError(
                 f"Seasons {seasons} do not match prior season extraction "
@@ -1597,7 +1599,7 @@ def _transform_to_lst_eager(
     """
     # Apart from the time index, all other dimensions will stay the same; this
     # is ensured with np.ogrid
-    idx = np.ogrid[tuple(slice(0, d) for d in data.shape)]
+    idx = list(np.ogrid[tuple(slice(0, d) for d in data.shape)])
     time_index = broadcast_to_shape(
         time_index, data.shape, (time_dim, lon_dim)
     )

pyproject.toml

@@ -50,7 +50,7 @@ dependencies = [
     "nc-time-axis",  # needed by iris.plot
     "nested-lookup",
     "netCDF4",
-    "numpy!=1.24.3,<2.0.0",  # avoid pulling 2.0.0rc1
+    "numpy!=1.24.3",
     "packaging",
     "pandas",
     "pillow",
@@ -61,7 +61,7 @@ dependencies = [
     "pyyaml",
     "requests",
     "scipy>=1.6",
-    "scitools-iris>=3.10.0",
+    "scitools-iris>=3.11",  # 3.11 first to support Numpy 2 and Python 3.13
     "shapely>=2.0.0",
     "stratify>=0.3",
     "yamale",

tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py

@@ -15,15 +15,15 @@
 from tests import PreprocessorFile
 
 
-def assert_allclose(array_1, array_2):
+def assert_allclose(array_1, array_2, rtol=1e-7):
     """Assert that (masked) array 1 is close to (masked) array 2."""
     if np.ma.is_masked(array_1) or np.ma.is_masked(array_2):
         mask_1 = np.ma.getmaskarray(array_1)
         mask_2 = np.ma.getmaskarray(array_2)
         np.testing.assert_equal(mask_1, mask_2)
-        np.testing.assert_allclose(array_1[~mask_1], array_2[~mask_2])
+        np.testing.assert_allclose(array_1[~mask_1], array_2[~mask_2], rtol)
     else:
-        np.testing.assert_allclose(array_1, array_2)
+        np.testing.assert_allclose(array_1, array_2, rtol)
 
 
 def products_set_to_dict(products):
@@ -473,7 +473,10 @@ def test_distance_metric(
     assert out_cube.shape == ()
     assert out_cube.dtype == np.float32
     assert not out_cube.has_lazy_data()
-    assert_allclose(out_cube.data, ref_data)
+    # an rtol=1e-6 is needed for numpy >=2.0
+    assert_allclose(
+        out_cube.data, np.array(ref_data, dtype=np.float32), rtol=1e-6
+    )
     assert out_cube.var_name == var_name
     assert out_cube.long_name == long_name
     assert out_cube.standard_name is None
@@ -684,7 +687,8 @@ def test_distance_metric_masked_data(
         expected_data = np.ma.masked_invalid(data)
     else:
         expected_data = np.array(data, dtype=np.float32)
-    assert_allclose(out_cube.data, expected_data)
+    # an rtol=1e-6 is needed for numpy >=2.0
+    assert_allclose(out_cube.data, expected_data, rtol=1e-6)
     assert out_cube.var_name == var_name
     assert out_cube.long_name == long_name
     assert out_cube.standard_name is None

tests/unit/preprocessor/_derive/test_shared.py

@@ -179,7 +179,7 @@ def test_low_lev_below_surf_press():
     """Test for lowest level below surface pressure."""
     plev = 970
     top_limit = 5
-    col = np.array([np.NaN, 900, 800])
+    col = np.array([np.nan, 900, 800])
     col = np.insert(col, 0, plev)
     col = np.append(col, top_limit)
     result = np.array([0, 120, 845])
@@ -197,7 +197,7 @@ def test_low_lev_below_surf_press():
         np.atleast_3d(result),
     )
 
-    col = np.array([np.NaN, np.NaN, 900, 800])
+    col = np.array([np.nan, np.nan, 900, 800])
     col = np.insert(col, 0, plev)
     col = np.append(col, top_limit)
     result = np.array([0, 0, 120, 845])

tests/unit/preprocessor/_regrid/__init__.py

@@ -11,7 +11,7 @@
 def _make_vcoord(data, dtype=None):
     """Create a synthetic test vertical coordinate."""
     if dtype is None:
-        dtype = np.dtype("int8")
+        dtype = np.int32
 
     if isinstance(data, int):
         data = np.arange(data, dtype=dtype)
@@ -45,7 +45,7 @@ def _make_cube(
 ):
     """Create a 3d synthetic test cube."""
     if dtype is None:
-        dtype = np.dtype("int8")
+        dtype = np.int32
 
     if not isinstance(data, np.ndarray):
         data = np.empty(data, dtype=dtype)

tests/unit/preprocessor/_regrid/test__create_cube.py

@@ -12,7 +12,7 @@
 class Test(tests.Test):
     def setUp(self):
         shape = (3, 2, 1)
-        self.dtype = np.dtype("int8")
+        self.dtype = np.int32
         self.cube = _make_cube(shape, dtype=self.dtype)
 
     def test_invalid_shape__data_mismatch_with_levels(self):

tests/unit/preprocessor/_regrid/test_extract_levels.py

@@ -24,7 +24,7 @@ class Test(tests.Test):
     def setUp(self):
         self.shape = (3, 2, 1)
         self.z = self.shape[0]
-        self.dtype = np.dtype("int8")
+        self.dtype = np.int32
         data = np.arange(np.prod(self.shape), dtype=self.dtype).reshape(
             self.shape
         )