Fix geotiff export geotransform values for numpy >= 2.4

docs/src/whatsnew/latest.rst

@@ -36,7 +36,10 @@ This document explains the changes made to Iris for this release
 🐛 Bugs Fixed
 =============
 
-#. N/A
+#. `@PraneelBhatia`_ fixed :func:`iris.experimental.raster.export_geotiff`,
+   which failed with a ``TypeError`` under numpy >= 2.4: the geotransform
+   values are now passed to GDAL as Python scalars instead of 1-element
+   arrays. (:issue:`7133`, :pull:`7140`)
 
 
 💣 Incompatible Changes
@@ -142,6 +145,7 @@ This document explains the changes made to Iris for this release
     core dev names are automatically included by the common_links.inc:
 
 .. _@hdyson: https://github.com/hdyson
+.. _@PraneelBhatia: https://github.com/PraneelBhatia
 
 
 .. comment

lib/iris/experimental/raster.py

@@ -173,7 +173,10 @@ def export_geotiff(cube, fname):
             )
         if not coord.is_contiguous():
             raise ValueError("Coordinate {!r} bounds must be contiguous.".format(name))
-        xy_step.append(np.diff(coord.bounds[0]))
+        # Take a Python scalar step value : GDAL's SetGeoTransform requires
+        # plain numbers, and numpy >= 2.4 no longer implicitly converts the
+        # 1-element array which np.diff produces here.
+        xy_step.append(np.diff(coord.bounds[0]).item())
         if not np.allclose(np.diff(coord.bounds), xy_step[-1]):
             msg = "Coordinate {!r} bounds must be regularly spaced.".format(name)
             raise ValueError(msg)
@@ -205,6 +208,6 @@ def export_geotiff(cube, fname):
             x_bounds = x_bounds.copy()
             x_bounds[big_indices] -= 360
 
-    x_min = np.min(x_bounds)
-    y_max = np.max(coord_y.bounds)
+    x_min = np.min(x_bounds).item()
+    y_max = np.max(coord_y.bounds).item()
     _gdal_write_array(x_min, x_step, y_max, y_step, coord_system, data, fname, "GTiff")

lib/iris/tests/unit/experimental/raster/test_export_geotiff.py

@@ -162,7 +162,7 @@ def test_ellipsoid(self, tmp_path):
 
 @_shared_utils.skip_gdal
 class TestGeoTransform:
-    def test_(self, tmp_path):
+    def _cube(self):
         data = np.arange(12).reshape(3, 4).astype(np.uint8)
         cube = Cube(data, "air_pressure_anomaly")
         coord = DimCoord([30, 40, 50], "latitude", units="degrees")
@@ -171,7 +171,23 @@ def test_(self, tmp_path):
         coord = DimCoord([-10, -5, 0, 5], "longitude", units="degrees")
         coord.guess_bounds()
         cube.add_dim_coord(coord, 1)
+        return cube
+
+    def test_geotransform(self, tmp_path):
+        cube = self._cube()
         temp_filename = str(tmp_path / "tmp.tif")
         export_geotiff(cube, temp_filename)
         dataset = gdal.Open(temp_filename, gdal.GA_ReadOnly)
         assert dataset.GetGeoTransform() == (-12.5, 5, 0, 55, 0, -10)
+
+    def test_scalar_transform_values(self, mocker, tmp_path):
+        # The geotransform values passed to GDAL must be Python scalars :
+        # numpy >= 2.4 no longer implicitly converts the 1-element arrays
+        # which were previously passed. See issue #7133.
+        write_array = mocker.patch("iris.experimental.raster._gdal_write_array")
+        export_geotiff(self._cube(), str(tmp_path / "tmp.tif"))
+        x_min, x_step, y_max, y_step = write_array.call_args.args[:4]
+        assert (x_min, x_step, y_max, y_step) == (-12.5, 5.0, 55.0, -10.0)
+        for value in (x_min, x_step, y_max, y_step):
+            assert isinstance(value, (int, float))
+            assert not isinstance(value, np.generic)