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

environment.yml

@@ -18,7 +18,7 @@ dependencies:
   - fire
   - geopy
   - humanfriendly
-  - iris >=3.10.0
+  - iris >=3.11  # 3.11 first to support numpy2 and Python3.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_esmpy.py

@@ -283,8 +283,11 @@ def get_grid(
         num_peri_dims = 1
     else:
         num_peri_dims = 0
+
+    # previous to numpy 2.0, a np.array(esmpy_lat.shape) was used
+    # numpy>=2.0 throws ValueError: matrix transpose with ndim<2 is undefined
     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,12 @@ 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
+
+        # explicitly casting to int32 since without it, in Numpy 2.0
+        # one gets OverflowError: Python integer 360 out of bounds for int8
+        # see notes https://numpy.org/devdocs/numpy_2_0_migration_guide.html
+        src_points_with_convex_hull[-2 * n_hull : -n_hull, 1] -= np.int32(360)
+        src_points_with_convex_hull[-n_hull:, 1] += np.int32(360)
 
         # Actual weights calculation
         (weights, indices) = self._calculate_weights(

esmvalcore/preprocessor/_time.py

@@ -1597,7 +1597,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 numpy2 and Python3.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
@@ -740,6 +744,7 @@ def test_distance_metric_fully_masked_data(
     assert out_cube.dtype == np.float64
 
     expected_data = np.ma.masked_all(())
+    print("out/in/metric", out_cube.data, expected_data, metric)
     assert_allclose(out_cube.data, expected_data)
     assert out_cube.var_name == var_name
     assert out_cube.long_name == long_name

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/_time/test_time.py

@@ -1063,9 +1063,12 @@ def test_season_not_available(self):
             name="clim_season",
             seasons=["JFMAMJ", "JASOND"],
         )
+
+        # numpy>=2.0 these need to be explicitly cast to numpy strings
+        two_seasons = [np.str_("JASOND"), np.str_("JFMAMJ")]
         msg = (
             "Seasons ('DJF', 'MAM', 'JJA', 'SON') do not match prior season "
-            "extraction ['JASOND', 'JFMAMJ']."
+            f"extraction {two_seasons}."
         )
         with pytest.raises(ValueError, match=re.escape(msg)):
             seasonal_statistics(cube, "mean")