diff --git a/lib/iris/tests/_shared_utils.py b/lib/iris/tests/_shared_utils.py
index a23c77a622..09fffcbf23 100644
--- a/lib/iris/tests/_shared_utils.py
+++ b/lib/iris/tests/_shared_utils.py
@@ -311,7 +311,7 @@ def assert_CML_approx_data(
             fname[-1] = fname[-1][:-4]
         fname[-1] += ".data.%d.json" % i
         assert_data_almost_equal(cube.data, fname, **kwargs)
-    assert_CML(cubes, reference_filename, checksum=False)
+    assert_CML(request, cubes, reference_filename, checksum=False)
 
 
 def assert_CDL(
@@ -459,8 +459,8 @@ def assert_data_almost_equal(data, reference_filename, **kwargs):
         kwargs.setdefault("err_msg", "Reference file %s" % reference_path)
         with open(reference_path, "r") as reference_file:
             stats = json.load(reference_file)
-            assert stats.get("shape", []), list(data.shape)
-            assert stats.get("masked", False), ma.is_masked(data)
+            assert stats.get("shape", []) == list(data.shape)
+            assert stats.get("masked", False) == ma.is_masked(data)
             nstats = np.array(
                 (
                     stats.get("mean", 0.0),
diff --git a/lib/iris/tests/test_analysis.py b/lib/iris/tests/test_analysis.py
index 6ed02b8ad4..abc68c06b5 100644
--- a/lib/iris/tests/test_analysis.py
+++ b/lib/iris/tests/test_analysis.py
@@ -3,10 +3,6 @@
 # This file is part of Iris and is released under the BSD license.
 # See LICENSE in the root of the repository for full licensing details.
 
-
-# import iris tests first so that some things can be initialised before importing anything else
-import iris.tests as tests  # isort:skip
-
 import cf_units
 import dask.array as da
 import numpy as np
@@ -20,12 +16,17 @@
 import iris.coord_systems
 import iris.coords
 import iris.cube
+from iris.tests import _shared_utils, stock
 import iris.tests.stock
 import iris.util
 
 
-class TestAnalysisCubeCoordComparison(tests.IrisTest):
-    def assertComparisonDict(self, comparison_dict, reference_filename):
+class TestAnalysisCubeCoordComparison:
+    @pytest.fixture(autouse=True)
+    def _setup(self, request):
+        self.request = request
+
+    def assert_comparison_dict(self, comparison_dict, reference_filename):
         string = ""
         for key in sorted(comparison_dict):
             coord_groups = comparison_dict[key]
@@ -36,7 +37,7 @@ def assertComparisonDict(self, comparison_dict, reference_filename):
             ]
             string += str(sorted(names))
             string += "\n"
-        self.assertString(string, reference_filename)
+        _shared_utils.assert_string(self.request, string, reference_filename)
 
     def test_coord_comparison(self):
         cube1 = iris.cube.Cube(np.zeros((41, 41)))
@@ -109,67 +110,71 @@ def test_coord_comparison(self):
 
         coord_comparison = iris.analysis._dimensional_metadata_comparison
 
-        self.assertComparisonDict(
+        self.assert_comparison_dict(
             coord_comparison(cube1, cube1),
             ("analysis", "coord_comparison", "cube1_cube1.txt"),
         )
-        self.assertComparisonDict(
+        self.assert_comparison_dict(
             coord_comparison(cube1, cube2),
             ("analysis", "coord_comparison", "cube1_cube2.txt"),
         )
-        self.assertComparisonDict(
+        self.assert_comparison_dict(
             coord_comparison(cube1, cube3),
             ("analysis", "coord_comparison", "cube1_cube3.txt"),
         )
-        self.assertComparisonDict(
+        self.assert_comparison_dict(
             coord_comparison(cube1, cube4),
             ("analysis", "coord_comparison", "cube1_cube4.txt"),
         )
-        self.assertComparisonDict(
+        self.assert_comparison_dict(
             coord_comparison(cube1, cube5),
             ("analysis", "coord_comparison", "cube1_cube5.txt"),
         )
-        self.assertComparisonDict(
+        self.assert_comparison_dict(
             coord_comparison(cube2, cube3),
             ("analysis", "coord_comparison", "cube2_cube3.txt"),
         )
-        self.assertComparisonDict(
+        self.assert_comparison_dict(
             coord_comparison(cube2, cube4),
             ("analysis", "coord_comparison", "cube2_cube4.txt"),
         )
-        self.assertComparisonDict(
+        self.assert_comparison_dict(
             coord_comparison(cube2, cube5),
             ("analysis", "coord_comparison", "cube2_cube5.txt"),
         )
-        self.assertComparisonDict(
+        self.assert_comparison_dict(
             coord_comparison(cube3, cube4),
             ("analysis", "coord_comparison", "cube3_cube4.txt"),
         )
-        self.assertComparisonDict(
+        self.assert_comparison_dict(
             coord_comparison(cube3, cube5),
             ("analysis", "coord_comparison", "cube3_cube5.txt"),
         )
-        self.assertComparisonDict(
+        self.assert_comparison_dict(
             coord_comparison(cube4, cube5),
             ("analysis", "coord_comparison", "cube4_cube5.txt"),
         )
 
-        self.assertComparisonDict(
+        self.assert_comparison_dict(
             coord_comparison(cube1, cube1, cube1),
             ("analysis", "coord_comparison", "cube1_cube1_cube1.txt"),
         )
-        self.assertComparisonDict(
+        self.assert_comparison_dict(
             coord_comparison(cube1, cube2, cube1),
             ("analysis", "coord_comparison", "cube1_cube2_cube1.txt"),
         )
 
         # get a coord comparison result and check that we are getting back what was expected
         coord_group = coord_comparison(cube1, cube2)["grouped_coords"][0]
-        self.assertIsInstance(coord_group, iris.analysis._CoordGroup)
-        self.assertIsInstance(list(coord_group)[0], iris.coords.Coord)
+        assert isinstance(coord_group, iris.analysis._CoordGroup)
+        assert isinstance(list(coord_group)[0], iris.coords.Coord)
 
 
-class TestAnalysisWeights(tests.IrisTest):
+class TestAnalysisWeights:
+    @pytest.fixture(autouse=True)
+    def _setup(self, request):
+        self.request = request
+
     def test_weighted_mean_little(self):
         data = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=np.float32)
         weights = np.array([[9, 8, 7], [6, 5, 4], [3, 2, 1]], dtype=np.float32)
@@ -196,7 +201,9 @@ def test_weighted_mean_little(self):
             ),
             1,
         )
-        self.assertCML(cube, ("analysis", "weighted_mean_source.cml"))
+        _shared_utils.assert_CML(
+            self.request, cube, ("analysis", "weighted_mean_source.cml")
+        )
 
         a = cube.collapsed("lat", iris.analysis.MEAN, weights=weights)
         # np.ma.average doesn't apply type promotion rules in some versions,
@@ -206,15 +213,19 @@ def test_weighted_mean_little(self):
         if a.dtype > np.float32:
             cast_data = a.data.astype(np.float32)
             a.data = cast_data
-        self.assertCMLApproxData(a, ("analysis", "weighted_mean_lat.cml"))
+        _shared_utils.assert_CML_approx_data(
+            self.request, a, ("analysis", "weighted_mean_lat.cml")
+        )
 
         b = cube.collapsed(lon_coord, iris.analysis.MEAN, weights=weights)
         if b.dtype > np.float32:
             cast_data = b.data.astype(np.float32)
             b.data = cast_data
         b.data = np.asarray(b.data)
-        self.assertCMLApproxData(b, ("analysis", "weighted_mean_lon.cml"))
-        self.assertEqual(b.coord("dummy").shape, (1,))
+        _shared_utils.assert_CML_approx_data(
+            self.request, b, ("analysis", "weighted_mean_lon.cml")
+        )
+        assert b.coord("dummy").shape == (1,)
 
         # test collapsing multiple coordinates (and the fact that one of the coordinates isn't the same coordinate instance as on the cube)
         c = cube.collapsed(
@@ -223,22 +234,26 @@ def test_weighted_mean_little(self):
         if c.dtype > np.float32:
             cast_data = c.data.astype(np.float32)
             c.data = cast_data
-        self.assertCMLApproxData(c, ("analysis", "weighted_mean_latlon.cml"))
-        self.assertEqual(c.coord("dummy").shape, (1,))
+        _shared_utils.assert_CML_approx_data(
+            self.request, c, ("analysis", "weighted_mean_latlon.cml")
+        )
+        assert c.coord("dummy").shape == (1,)
 
         # Check new coord bounds - made from points
-        self.assertArrayEqual(c.coord("lat").bounds, [[1, 3]])
+        _shared_utils.assert_array_equal(c.coord("lat").bounds, [[1, 3]])
 
         # Check new coord bounds - made from bounds
         cube.coord("lat").bounds = [[0.5, 1.5], [1.5, 2.5], [2.5, 3.5]]
         c = cube.collapsed(["lat", "lon"], iris.analysis.MEAN, weights=weights)
-        self.assertArrayEqual(c.coord("lat").bounds, [[0.5, 3.5]])
+        _shared_utils.assert_array_equal(c.coord("lat").bounds, [[0.5, 3.5]])
         cube.coord("lat").bounds = None
 
         # Check there was no residual change
-        self.assertCML(cube, ("analysis", "weighted_mean_source.cml"))
+        _shared_utils.assert_CML(
+            self.request, cube, ("analysis", "weighted_mean_source.cml")
+        )
 
-    @tests.skip_data
+    @_shared_utils.skip_data
     def test_weighted_mean(self):
         # compare with pp_area_avg - which collapses both lat and lon
         #
@@ -247,8 +262,10 @@ def test_weighted_mean(self):
         #     print, pp_area_avg(pp, /box)  #287.927
         #     ;gives an answer of 287.927
         #
-        e = iris.tests.stock.simple_pp()
-        self.assertCML(e, ("analysis", "weighted_mean_original.cml"))
+        e = stock.simple_pp()
+        _shared_utils.assert_CML(
+            self.request, e, ("analysis", "weighted_mean_original.cml")
+        )
         e.coord("latitude").guess_bounds()
         e.coord("longitude").guess_bounds()
         area_weights = iris.analysis.cartography.area_weights(e)
@@ -259,9 +276,9 @@ def test_weighted_mean(self):
         )
         g = f.collapsed("longitude", iris.analysis.MEAN, weights=collapsed_area_weights)
         # check it's a 0d, scalar cube
-        self.assertEqual(g.shape, ())
+        assert g.shape == ()
         # check the value - pp_area_avg's result of 287.927 differs by factor of 1.00002959
-        np.testing.assert_approx_equal(g.data, 287.935, significant=5)
+        _shared_utils.assert_array_almost_equal(g.data, 287.935, decimal=2)
 
         # check we get summed weights even if we don't give any
         h, summed_weights = e.collapsed("latitude", iris.analysis.MEAN, returned=True)
@@ -270,10 +287,12 @@ def test_weighted_mean(self):
         # Check there was no residual change
         e.coord("latitude").bounds = None
         e.coord("longitude").bounds = None
-        self.assertCML(e, ("analysis", "weighted_mean_original.cml"))
+        _shared_utils.assert_CML(
+            self.request, e, ("analysis", "weighted_mean_original.cml")
+        )
 
         # Test collapsing of missing coord
-        self.assertRaises(
+        pytest.raises(
             iris.exceptions.CoordinateNotFoundError,
             e.collapsed,
             "platitude",
@@ -281,7 +300,7 @@ def test_weighted_mean(self):
         )
 
         # Test collapsing of non data coord
-        self.assertRaises(
+        pytest.raises(
             iris.exceptions.CoordinateCollapseError,
             e.collapsed,
             "pressure",
@@ -289,13 +308,15 @@ def test_weighted_mean(self):
         )
 
 
-@tests.skip_data
-class TestAnalysisBasic(tests.IrisTest):
-    def setUp(self):
-        file = tests.get_data_path(("PP", "aPProt1", "rotatedMHtimecube.pp"))
+@_shared_utils.skip_data
+class TestAnalysisBasic:
+    @pytest.fixture(autouse=True)
+    def _setup(self, request):
+        self.request = request
+        file = _shared_utils.get_data_path(("PP", "aPProt1", "rotatedMHtimecube.pp"))
         cubes = iris.load(file)
         self.cube = cubes[0]
-        self.assertCML(self.cube, ("analysis", "original.cml"))
+        _shared_utils.assert_CML(self.request, self.cube, ("analysis", "original.cml"))
 
     def _common(
         self,
@@ -307,15 +328,16 @@ def _common(
     ):
         self.cube.data = self.cube.data.astype(np.float64)
 
-        self.assertCML(self.cube, ("analysis", original_name))
+        _shared_utils.assert_CML(self.request, self.cube, ("analysis", original_name))
 
         a = self.cube.collapsed("grid_latitude", aggregate)
-        self.assertCMLApproxData(
-            a, ("analysis", "%s_latitude.cml" % name), *args, **kwargs
+        _shared_utils.assert_CML_approx_data(
+            self.request, a, ("analysis", "%s_latitude.cml" % name), *args, **kwargs
         )
 
         b = a.collapsed("grid_longitude", aggregate)
-        self.assertCMLApproxData(
+        _shared_utils.assert_CML_approx_data(
+            self.request,
             b,
             ("analysis", "%s_latitude_longitude.cml" % name),
             *args,
@@ -323,7 +345,8 @@ def _common(
         )
 
         c = self.cube.collapsed(["grid_latitude", "grid_longitude"], aggregate)
-        self.assertCMLApproxData(
+        _shared_utils.assert_CML_approx_data(
+            self.request,
             c,
             ("analysis", "%s_latitude_longitude_1call.cml" % name),
             *args,
@@ -331,7 +354,7 @@ def _common(
         )
 
         # Check there was no residual change
-        self.assertCML(self.cube, ("analysis", original_name))
+        _shared_utils.assert_CML(self.request, self.cube, ("analysis", original_name))
 
     def test_mean(self):
         self._common("mean", iris.analysis.MEAN, rtol=1e-05)
@@ -369,12 +392,13 @@ def test_rms(self):
         self._common("rms", iris.analysis.RMS)
 
     def test_duplicate_coords(self):
-        self.assertRaises(ValueError, tests.stock.track_1d, duplicate_x=True)
+        pytest.raises(ValueError, stock.track_1d, duplicate_x=True)
 
 
-class TestMissingData(tests.IrisTest):
-    def setUp(self):
-        self.cube_with_nan = tests.stock.simple_2d()
+class TestMissingData:
+    @pytest.fixture(autouse=True)
+    def _setup(self):
+        self.cube_with_nan = stock.simple_2d()
 
         data = self.cube_with_nan.data.astype(np.float32)
         self.cube_with_nan.data = data.copy()
@@ -382,36 +406,37 @@ def setUp(self):
         self.cube_with_nan.data[2, 2] = np.nan
         self.cube_with_nan.data[2, 3] = np.nan
 
-        self.cube_with_mask = tests.stock.simple_2d()
+        self.cube_with_mask = stock.simple_2d()
         self.cube_with_mask.data = ma.array(
             self.cube_with_nan.data, mask=np.isnan(self.cube_with_nan.data)
         )
 
     def test_max(self):
         cube = self.cube_with_nan.collapsed("foo", iris.analysis.MAX)
-        np.testing.assert_array_equal(cube.data, np.array([3, np.nan, np.nan]))
+        _shared_utils.assert_array_equal(cube.data, np.array([3, np.nan, np.nan]))
 
         cube = self.cube_with_mask.collapsed("foo", iris.analysis.MAX)
-        np.testing.assert_array_equal(cube.data, np.array([3, 7, 9]))
+        _shared_utils.assert_array_equal(cube.data, np.array([3, 7, 9]))
 
     def test_min(self):
         cube = self.cube_with_nan.collapsed("foo", iris.analysis.MIN)
-        np.testing.assert_array_equal(cube.data, np.array([0, np.nan, np.nan]))
+        _shared_utils.assert_array_equal(cube.data, np.array([0, np.nan, np.nan]))
 
         cube = self.cube_with_mask.collapsed("foo", iris.analysis.MIN)
-        np.testing.assert_array_equal(cube.data, np.array([0, 5, 8]))
+        _shared_utils.assert_array_equal(cube.data, np.array([0, 5, 8]))
 
     def test_sum(self):
         cube = self.cube_with_nan.collapsed("foo", iris.analysis.SUM)
-        np.testing.assert_array_equal(cube.data, np.array([6, np.nan, np.nan]))
+        _shared_utils.assert_array_equal(cube.data, np.array([6, np.nan, np.nan]))
 
         cube = self.cube_with_mask.collapsed("foo", iris.analysis.SUM)
-        np.testing.assert_array_equal(cube.data, np.array([6, 18, 17]))
+        _shared_utils.assert_array_equal(cube.data, np.array([6, 18, 17]))
 
 
-class TestAuxCoordCollapse(tests.IrisTest):
-    def setUp(self):
-        self.cube_with_aux_coord = tests.stock.simple_4d_with_hybrid_height()
+class TestAuxCoordCollapse:
+    @pytest.fixture(autouse=True)
+    def _setup(self):
+        self.cube_with_aux_coord = stock.simple_4d_with_hybrid_height()
 
         # Guess bounds to get the weights
         self.cube_with_aux_coord.coord("grid_latitude").guess_bounds()
@@ -419,12 +444,12 @@ def setUp(self):
 
     def test_max(self):
         cube = self.cube_with_aux_coord.collapsed("grid_latitude", iris.analysis.MAX)
-        np.testing.assert_array_equal(
+        _shared_utils.assert_array_equal(
             cube.coord("surface_altitude").points,
             np.array([112, 113, 114, 115, 116, 117]),
         )
 
-        np.testing.assert_array_equal(
+        _shared_utils.assert_array_equal(
             cube.coord("surface_altitude").bounds,
             np.array(
                 [
@@ -441,29 +466,30 @@ def test_max(self):
         # Check collapsing over the whole coord still works
         cube = self.cube_with_aux_coord.collapsed("altitude", iris.analysis.MAX)
 
-        np.testing.assert_array_equal(
+        _shared_utils.assert_array_equal(
             cube.coord("surface_altitude").points, np.array([114])
         )
 
-        np.testing.assert_array_equal(
+        _shared_utils.assert_array_equal(
             cube.coord("surface_altitude").bounds, np.array([[100, 129]])
         )
 
         cube = self.cube_with_aux_coord.collapsed("grid_longitude", iris.analysis.MAX)
 
-        np.testing.assert_array_equal(
+        _shared_utils.assert_array_equal(
             cube.coord("surface_altitude").points,
             np.array([102, 108, 114, 120, 126]),
         )
 
-        np.testing.assert_array_equal(
+        _shared_utils.assert_array_equal(
             cube.coord("surface_altitude").bounds,
             np.array([[100, 105], [106, 111], [112, 117], [118, 123], [124, 129]]),
         )
 
 
-class TestAggregator_mdtol_keyword(tests.IrisTest):
-    def setUp(self):
+class TestAggregator_mdtol_keyword:
+    @pytest.fixture(autouse=True)
+    def _setup(self):
         data = ma.array(
             [[1, 2], [4, 5]],
             dtype=np.float32,
@@ -483,7 +509,7 @@ def setUp(self):
     def test_single_coord_no_mdtol(self):
         collapsed = self.cube.collapsed(self.cube.coord("lat"), iris.analysis.MEAN)
         t = ma.array([2.5, 5.0], mask=[False, True])
-        self.assertMaskedArrayEqual(collapsed.data, t)
+        _shared_utils.assert_masked_array_almost_equal(collapsed.data, t)
 
     def test_single_coord_mdtol(self):
         self.cube.data.mask = np.array([[False, True], [False, False]])
@@ -491,7 +517,7 @@ def test_single_coord_mdtol(self):
             self.cube.coord("lat"), iris.analysis.MEAN, mdtol=0.5
         )
         t = ma.array([2.5, 5], mask=[False, False])
-        self.assertMaskedArrayEqual(collapsed.data, t)
+        _shared_utils.assert_masked_array_almost_equal(collapsed.data, t)
 
     def test_single_coord_mdtol_alt(self):
         self.cube.data.mask = np.array([[False, True], [False, False]])
@@ -499,7 +525,7 @@ def test_single_coord_mdtol_alt(self):
             self.cube.coord("lat"), iris.analysis.MEAN, mdtol=0.4
         )
         t = ma.array([2.5, 5], mask=[False, True])
-        self.assertMaskedArrayEqual(collapsed.data, t)
+        _shared_utils.assert_masked_array_almost_equal(collapsed.data, t)
 
     def test_multi_coord_no_mdtol(self):
         collapsed = self.cube.collapsed(
@@ -507,7 +533,7 @@ def test_multi_coord_no_mdtol(self):
             iris.analysis.MEAN,
         )
         t = np.array(2.5)
-        self.assertArrayEqual(collapsed.data, t)
+        _shared_utils.assert_array_equal(collapsed.data, t)
 
     def test_multi_coord_mdtol(self):
         collapsed = self.cube.collapsed(
@@ -516,10 +542,14 @@ def test_multi_coord_mdtol(self):
             mdtol=0.4,
         )
         t = ma.array(2.5, mask=True)
-        self.assertMaskedArrayEqual(collapsed.data, t)
+        _shared_utils.assert_masked_array_almost_equal(collapsed.data, t)
 
 
-class TestAggregators(tests.IrisTest):
+class TestAggregators:
+    @pytest.fixture(autouse=True)
+    def _setup(self, request):
+        self.request = request
+
     def _check_collapsed_percentile(
         self,
         cube,
@@ -537,30 +567,32 @@ def _check_collapsed_percentile(
             percent=percents,
             **kwargs,
         )
-        np.testing.assert_array_almost_equal(result.data, expected_result)
-        self.assertEqual(type(result.data), cube_data_type)
+        _shared_utils.assert_array_almost_equal(result.data, expected_result)
+        assert type(result.data) is cube_data_type
         if CML_filename is not None:
-            self.assertCML(result, ("analysis", CML_filename), checksum=False)
+            _shared_utils.assert_CML(
+                self.request, result, ("analysis", CML_filename), checksum=False
+            )
 
     def _check_percentile(self, data, axis, percents, expected_result, **kwargs):
         result = iris.analysis._percentile(data, axis, percents, **kwargs)
-        np.testing.assert_array_almost_equal(result, expected_result)
-        self.assertEqual(type(result), type(expected_result))
+        _shared_utils.assert_array_almost_equal(result, expected_result)
+        assert type(result) is type(expected_result)
 
     def test_percentile_1d_25_percent(self):
-        cube = tests.stock.simple_1d()
+        cube = stock.simple_1d()
         self._check_collapsed_percentile(
             cube, 25, "foo", 2.5, CML_filename="first_quartile_foo_1d.cml"
         )
 
     def test_percentile_1d_75_percent(self):
-        cube = tests.stock.simple_1d()
+        cube = stock.simple_1d()
         self._check_collapsed_percentile(
             cube, 75, "foo", 7.5, CML_filename="third_quartile_foo_1d.cml"
         )
 
     def test_fast_percentile_1d_25_percent(self):
-        cube = tests.stock.simple_1d()
+        cube = stock.simple_1d()
         self._check_collapsed_percentile(
             cube,
             25,
@@ -571,7 +603,7 @@ def test_fast_percentile_1d_25_percent(self):
         )
 
     def test_fast_percentile_1d_75_percent(self):
-        cube = tests.stock.simple_1d()
+        cube = stock.simple_1d()
         self._check_collapsed_percentile(
             cube,
             75,
@@ -582,7 +614,7 @@ def test_fast_percentile_1d_75_percent(self):
         )
 
     def test_fast_percentile_1d_75_percent_masked_type_no_mask(self):
-        cube = tests.stock.simple_1d()
+        cube = stock.simple_1d()
         cube.data = ma.MaskedArray(cube.data)
         self._check_collapsed_percentile(
             cube,
@@ -594,7 +626,7 @@ def test_fast_percentile_1d_75_percent_masked_type_no_mask(self):
         )
 
     def test_percentile_2d_single_coord(self):
-        cube = tests.stock.simple_2d()
+        cube = stock.simple_2d()
         self._check_collapsed_percentile(
             cube,
             25,
@@ -604,7 +636,7 @@ def test_percentile_2d_single_coord(self):
         )
 
     def test_percentile_2d_two_coords(self):
-        cube = tests.stock.simple_2d()
+        cube = stock.simple_2d()
         self._check_collapsed_percentile(
             cube,
             25,
@@ -614,7 +646,7 @@ def test_percentile_2d_two_coords(self):
         )
 
     def test_fast_percentile_2d_single_coord(self):
-        cube = tests.stock.simple_2d()
+        cube = stock.simple_2d()
         self._check_collapsed_percentile(
             cube,
             25,
@@ -625,7 +657,7 @@ def test_fast_percentile_2d_single_coord(self):
         )
 
     def test_fast_percentile_2d_two_coords(self):
-        cube = tests.stock.simple_2d()
+        cube = stock.simple_2d()
         self._check_collapsed_percentile(
             cube,
             25,
@@ -636,7 +668,7 @@ def test_fast_percentile_2d_two_coords(self):
         )
 
     def test_fast_percentile_2d_single_coord_masked_type_no_mask(self):
-        cube = tests.stock.simple_2d()
+        cube = stock.simple_2d()
         cube.data = ma.MaskedArray(cube.data)
         self._check_collapsed_percentile(
             cube,
@@ -648,7 +680,7 @@ def test_fast_percentile_2d_single_coord_masked_type_no_mask(self):
         )
 
     def test_fast_percentile_2d_two_coords_masked_type_no_mask(self):
-        cube = tests.stock.simple_2d()
+        cube = stock.simple_2d()
         cube.data = ma.MaskedArray(cube.data)
         self._check_collapsed_percentile(
             cube,
@@ -744,7 +776,7 @@ def test_fast_percentile_3d_axis_two_masked_type_no_mask(self):
         )
 
     def test_percentile_3d_masked(self):
-        cube = tests.stock.simple_3d_mask()
+        cube = stock.simple_3d_mask()
         expected_result = [
             [12.0, 13.0, 14.0, 15.0],
             [16.0, 17.0, 18.0, 19.0],
@@ -760,10 +792,10 @@ def test_percentile_3d_masked(self):
         )
 
     def test_fast_percentile_3d_masked_type_masked(self):
-        cube = tests.stock.simple_3d_mask()
+        cube = stock.simple_3d_mask()
         msg = "Cannot use fast np.percentile method with masked array."
 
-        with self.assertRaisesRegex(TypeError, msg):
+        with pytest.raises(TypeError, match=msg):
             cube.collapsed(
                 "wibble",
                 iris.analysis.PERCENTILE,
@@ -772,7 +804,7 @@ def test_fast_percentile_3d_masked_type_masked(self):
             )
 
     def test_percentile_3d_notmasked(self):
-        cube = tests.stock.simple_3d()
+        cube = stock.simple_3d()
         expected_result = [
             [9.0, 10.0, 11.0, 12.0],
             [13.0, 14.0, 15.0, 16.0],
@@ -788,7 +820,7 @@ def test_percentile_3d_notmasked(self):
         )
 
     def test_fast_percentile_3d_notmasked(self):
-        cube = tests.stock.simple_3d()
+        cube = stock.simple_3d()
         expected_result = [
             [9.0, 10.0, 11.0, 12.0],
             [13.0, 14.0, 15.0, 16.0],
@@ -805,42 +837,50 @@ def test_fast_percentile_3d_notmasked(self):
         )
 
     def test_proportion(self):
-        cube = tests.stock.simple_1d()
+        cube = stock.simple_1d()
         assert np.any(cube.data >= 5)
         gt5 = cube.collapsed(
             "foo", iris.analysis.PROPORTION, function=lambda val: val >= 5
         )
-        np.testing.assert_array_almost_equal(gt5.data, np.array([6 / 11.0]))
-        self.assertCML(gt5, ("analysis", "proportion_foo_1d.cml"), checksum=False)
+        _shared_utils.assert_array_almost_equal(gt5.data, np.array([6 / 11.0]))
+        _shared_utils.assert_CML(
+            self.request, gt5, ("analysis", "proportion_foo_1d.cml"), checksum=False
+        )
 
     def test_proportion_2d(self):
-        cube = tests.stock.simple_2d()
+        cube = stock.simple_2d()
 
         gt6 = cube.collapsed(
             "foo", iris.analysis.PROPORTION, function=lambda val: val >= 6
         )
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             gt6.data, np.array([0, 0.5, 1], dtype=np.float32)
         )
-        self.assertCML(gt6, ("analysis", "proportion_foo_2d.cml"), checksum=False)
+        _shared_utils.assert_CML(
+            self.request, gt6, ("analysis", "proportion_foo_2d.cml"), checksum=False
+        )
 
         gt6 = cube.collapsed(
             "bar", iris.analysis.PROPORTION, function=lambda val: val >= 6
         )
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             gt6.data, np.array([1 / 3, 1 / 3, 2 / 3, 2 / 3], dtype=np.float32)
         )
-        self.assertCML(gt6, ("analysis", "proportion_bar_2d.cml"), checksum=False)
+        _shared_utils.assert_CML(
+            self.request, gt6, ("analysis", "proportion_bar_2d.cml"), checksum=False
+        )
 
         gt6 = cube.collapsed(
             ("foo", "bar"),
             iris.analysis.PROPORTION,
             function=lambda val: val >= 6,
         )
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             gt6.data, np.array([0.5], dtype=np.float32)
         )
-        self.assertCML(gt6, ("analysis", "proportion_foo_bar_2d.cml"), checksum=False)
+        _shared_utils.assert_CML(
+            self.request, gt6, ("analysis", "proportion_foo_bar_2d.cml"), checksum=False
+        )
 
         # mask the data
         cube.data = ma.array(cube.data, mask=cube.data % 2)
@@ -848,7 +888,7 @@ def test_proportion_2d(self):
         gt6_masked = cube.collapsed(
             "bar", iris.analysis.PROPORTION, function=lambda val: val >= 6
         )
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             gt6_masked.data,
             ma.array(
                 [1 / 3, None, 1 / 2, None],
@@ -856,58 +896,71 @@ def test_proportion_2d(self):
                 dtype=np.float32,
             ),
         )
-        self.assertCML(
+        _shared_utils.assert_CML(
+            self.request,
             gt6_masked,
             ("analysis", "proportion_foo_2d_masked.cml"),
             checksum=False,
         )
 
     def test_count(self):
-        cube = tests.stock.simple_1d()
+        cube = stock.simple_1d()
         gt5 = cube.collapsed("foo", iris.analysis.COUNT, function=lambda val: val >= 5)
-        np.testing.assert_array_almost_equal(gt5.data, np.array([6]))
+        _shared_utils.assert_array_almost_equal(gt5.data, np.array([6]))
         gt5.data = gt5.data.astype("i8")
-        self.assertCML(gt5, ("analysis", "count_foo_1d.cml"), checksum=False)
+        _shared_utils.assert_CML(
+            self.request, gt5, ("analysis", "count_foo_1d.cml"), checksum=False
+        )
 
     def test_count_2d(self):
-        cube = tests.stock.simple_2d()
+        cube = stock.simple_2d()
 
         gt6 = cube.collapsed("foo", iris.analysis.COUNT, function=lambda val: val >= 6)
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             gt6.data, np.array([0, 2, 4], dtype=np.float32)
         )
         gt6.data = gt6.data.astype("i8")
-        self.assertCML(gt6, ("analysis", "count_foo_2d.cml"), checksum=False)
+        _shared_utils.assert_CML(
+            self.request, gt6, ("analysis", "count_foo_2d.cml"), checksum=False
+        )
 
         gt6 = cube.collapsed("bar", iris.analysis.COUNT, function=lambda val: val >= 6)
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             gt6.data, np.array([1, 1, 2, 2], dtype=np.float32)
         )
         gt6.data = gt6.data.astype("i8")
-        self.assertCML(gt6, ("analysis", "count_bar_2d.cml"), checksum=False)
+        _shared_utils.assert_CML(
+            self.request, gt6, ("analysis", "count_bar_2d.cml"), checksum=False
+        )
 
         gt6 = cube.collapsed(
             ("foo", "bar"), iris.analysis.COUNT, function=lambda val: val >= 6
         )
-        np.testing.assert_array_almost_equal(gt6.data, np.array([6], dtype=np.float32))
+        _shared_utils.assert_array_almost_equal(
+            gt6.data, np.array([6], dtype=np.float32)
+        )
         gt6.data = gt6.data.astype("i8")
-        self.assertCML(gt6, ("analysis", "count_foo_bar_2d.cml"), checksum=False)
+        _shared_utils.assert_CML(
+            self.request, gt6, ("analysis", "count_foo_bar_2d.cml"), checksum=False
+        )
 
     def test_max_run_1d(self):
-        cube = tests.stock.simple_1d()
+        cube = stock.simple_1d()
         # [ 0  1  2  3  4  5  6  7  8  9 10]
         result = cube.collapsed(
             "foo",
             iris.analysis.MAX_RUN,
             function=lambda val: np.isin(val, [0, 1, 4, 5, 6, 8, 9]),
         )
-        self.assertArrayEqual(result.data, np.array(3))
-        self.assertEqual(result.units, 1)
-        self.assertTupleEqual(result.cell_methods, ())
-        self.assertCML(result, ("analysis", "max_run_foo_1d.cml"), checksum=False)
+        _shared_utils.assert_array_equal(result.data, np.array(3))
+        assert result.units == 1
+        assert result.cell_methods == ()
+        _shared_utils.assert_CML(
+            self.request, result, ("analysis", "max_run_foo_1d.cml"), checksum=False
+        )
 
     def test_max_run_lazy(self):
-        cube = tests.stock.simple_1d()
+        cube = stock.simple_1d()
         # [ 0  1  2  3  4  5  6  7  8  9 10]
         # Make data lazy
         cube.data = da.from_array(cube.data)
@@ -916,16 +969,18 @@ def test_max_run_lazy(self):
             iris.analysis.MAX_RUN,
             function=lambda val: np.isin(val, [0, 1, 4, 5, 6, 8, 9]),
         )
-        self.assertTrue(result.has_lazy_data())
+        assert result.has_lazy_data()
         # Realise data
         _ = result.data
-        self.assertArrayEqual(result.data, np.array(3))
-        self.assertEqual(result.units, 1)
-        self.assertTupleEqual(result.cell_methods, ())
-        self.assertCML(result, ("analysis", "max_run_foo_1d.cml"), checksum=False)
+        _shared_utils.assert_array_equal(result.data, np.array(3))
+        assert result.units == 1
+        assert result.cell_methods == ()
+        _shared_utils.assert_CML(
+            self.request, result, ("analysis", "max_run_foo_1d.cml"), checksum=False
+        )
 
     def test_max_run_2d(self):
-        cube = tests.stock.simple_2d()
+        cube = stock.simple_2d()
         # [[ 0  1  2  3]
         #  [ 4  5  6  7]
         #  [ 8  9 10 11]]
@@ -934,22 +989,30 @@ def test_max_run_2d(self):
             iris.analysis.MAX_RUN,
             function=lambda val: np.isin(val, [0, 3, 4, 5, 7, 9, 11]),
         )
-        self.assertArrayEqual(foo_result.data, np.array([1, 2, 1], dtype=np.float32))
-        self.assertEqual(foo_result.units, 1)
-        self.assertTupleEqual(foo_result.cell_methods, ())
-        self.assertCML(foo_result, ("analysis", "max_run_foo_2d.cml"), checksum=False)
+        _shared_utils.assert_array_equal(
+            foo_result.data, np.array([1, 2, 1], dtype=np.float32)
+        )
+        assert foo_result.units == 1
+        assert foo_result.cell_methods == ()
+        _shared_utils.assert_CML(
+            self.request, foo_result, ("analysis", "max_run_foo_2d.cml"), checksum=False
+        )
 
         bar_result = cube.collapsed(
             "bar",
             iris.analysis.MAX_RUN,
             function=lambda val: np.isin(val, [0, 3, 4, 5, 7, 9, 11]),
         )
-        self.assertArrayEqual(bar_result.data, np.array([2, 2, 0, 3], dtype=np.float32))
-        self.assertEqual(bar_result.units, 1)
-        self.assertTupleEqual(bar_result.cell_methods, ())
-        self.assertCML(bar_result, ("analysis", "max_run_bar_2d.cml"), checksum=False)
+        _shared_utils.assert_array_equal(
+            bar_result.data, np.array([2, 2, 0, 3], dtype=np.float32)
+        )
+        assert bar_result.units == 1
+        assert bar_result.cell_methods == ()
+        _shared_utils.assert_CML(
+            self.request, bar_result, ("analysis", "max_run_bar_2d.cml"), checksum=False
+        )
 
-        with self.assertRaises(ValueError):
+        with pytest.raises(ValueError):
             _ = cube.collapsed(
                 ("foo", "bar"),
                 iris.analysis.MAX_RUN,
@@ -957,7 +1020,7 @@ def test_max_run_2d(self):
             )
 
     def test_max_run_masked(self):
-        cube = tests.stock.simple_2d()
+        cube = stock.simple_2d()
         # [[ 0  1  2  3]
         #  [ 4  5  6  7]
         #  [ 8  9 10 11]]
@@ -972,40 +1035,51 @@ def test_max_run_masked(self):
             iris.analysis.MAX_RUN,
             function=lambda val: np.isin(val, [0, 1, 4, 5, 6, 9, 10, 11]),
         )
-        self.assertArrayEqual(result.data, np.array([1, 1, 2, 0], dtype=np.float32))
-        self.assertEqual(result.units, 1)
-        self.assertTupleEqual(result.cell_methods, ())
-        self.assertCML(
-            result, ("analysis", "max_run_bar_2d_masked.cml"), checksum=False
+        _shared_utils.assert_array_equal(
+            result.data, np.array([1, 1, 2, 0], dtype=np.float32)
+        )
+        assert result.units == 1
+        assert result.cell_methods == ()
+        _shared_utils.assert_CML(
+            self.request,
+            result,
+            ("analysis", "max_run_bar_2d_masked.cml"),
+            checksum=False,
         )
 
     def test_weighted_sum_consistency(self):
         # weighted sum with unit weights should be the same as a sum
-        cube = tests.stock.simple_1d()
+        cube = stock.simple_1d()
         normal_sum = cube.collapsed("foo", iris.analysis.SUM)
         weights = np.ones_like(cube.data)
         weighted_sum = cube.collapsed("foo", iris.analysis.SUM, weights=weights)
-        self.assertArrayAlmostEqual(normal_sum.data, weighted_sum.data)
+        _shared_utils.assert_array_almost_equal(normal_sum.data, weighted_sum.data)
 
     def test_weighted_sum_1d(self):
         # verify 1d weighted sum is correct
-        cube = tests.stock.simple_1d()
+        cube = stock.simple_1d()
         weights = np.array([0.05, 0.05, 0.1, 0.1, 0.2, 0.3, 0.2, 0.1, 0.1, 0.05, 0.05])
         result = cube.collapsed("foo", iris.analysis.SUM, weights=weights)
-        self.assertAlmostEqual(result.data, 6.5)
-        self.assertCML(result, ("analysis", "sum_weighted_1d.cml"), checksum=False)
+        assert result.data == pytest.approx(6.5)
+        _shared_utils.assert_CML(
+            self.request, result, ("analysis", "sum_weighted_1d.cml"), checksum=False
+        )
 
     def test_weighted_sum_2d(self):
         # verify 2d weighted sum is correct
-        cube = tests.stock.simple_2d()
+        cube = stock.simple_2d()
         weights = np.array([0.3, 0.4, 0.3])
         weights = iris.util.broadcast_to_shape(weights, cube.shape, [0])
         result = cube.collapsed("bar", iris.analysis.SUM, weights=weights)
-        self.assertArrayAlmostEqual(result.data, np.array([4.0, 5.0, 6.0, 7.0]))
-        self.assertCML(result, ("analysis", "sum_weighted_2d.cml"), checksum=False)
+        _shared_utils.assert_array_almost_equal(
+            result.data, np.array([4.0, 5.0, 6.0, 7.0])
+        )
+        _shared_utils.assert_CML(
+            self.request, result, ("analysis", "sum_weighted_2d.cml"), checksum=False
+        )
 
     def test_weighted_rms(self):
-        cube = tests.stock.simple_2d()
+        cube = stock.simple_2d()
         # modify cube data so that the results are nice numbers
         cube.data = np.array(
             [[4, 7, 10, 8], [21, 30, 12, 24], [14, 16, 20, 8]],
@@ -1017,12 +1091,14 @@ def test_weighted_rms(self):
         )
         expected_result = np.array([8.0, 24.0, 16.0])
         result = cube.collapsed("foo", iris.analysis.RMS, weights=weights)
-        self.assertArrayAlmostEqual(result.data, expected_result)
-        self.assertCML(result, ("analysis", "rms_weighted_2d.cml"), checksum=False)
+        _shared_utils.assert_array_almost_equal(result.data, expected_result)
+        _shared_utils.assert_CML(
+            self.request, result, ("analysis", "rms_weighted_2d.cml"), checksum=False
+        )
 
 
-@tests.skip_data
-class TestRotatedPole(tests.IrisTest):
+@_shared_utils.skip_data
+class TestRotatedPole:
     def _check_both_conversions(self, cube, index):
         rlons, rlats = iris.analysis.cartography.get_xy_grids(cube)
         rcs = cube.coord_system("RotatedGeogCS")
@@ -1032,21 +1108,21 @@ def _check_both_conversions(self, cube, index):
             rcs.grid_north_pole_longitude,
             rcs.grid_north_pole_latitude,
         )
-        self.assertDataAlmostEqual(
+        _shared_utils.assert_data_almost_equal(
             x, ("analysis", "rotated_pole.{}.x.json".format(index))
         )
-        self.assertDataAlmostEqual(
+        _shared_utils.assert_data_almost_equal(
             y, ("analysis", "rotated_pole.{}.y.json".format(index))
         )
-        self.assertDataAlmostEqual(
+        _shared_utils.assert_data_almost_equal(
             rlons, ("analysis", "rotated_pole.{}.rlon.json".format(index))
         )
-        self.assertDataAlmostEqual(
+        _shared_utils.assert_data_almost_equal(
             rlats, ("analysis", "rotated_pole.{}.rlat.json".format(index))
         )
 
     def test_all(self):
-        path = tests.get_data_path(("PP", "ukVorog", "ukv_orog_refonly.pp"))
+        path = _shared_utils.get_data_path(("PP", "ukVorog", "ukv_orog_refonly.pp"))
         master_cube = iris.load_cube(path)
 
         # Check overall behaviour.
@@ -1067,8 +1143,8 @@ def test_unrotate_nd(self):
         solx = np.array([[-16.42176094, -14.85892262], [-16.71055023, -14.58434624]])
         soly = np.array([[46.00724251, 51.29188893], [46.98728486, 50.30706042]])
 
-        self.assertArrayAlmostEqual(resx, solx)
-        self.assertArrayAlmostEqual(resy, soly)
+        _shared_utils.assert_array_almost_equal(resx, solx)
+        _shared_utils.assert_array_almost_equal(resy, soly)
 
     def test_unrotate_1d(self):
         rlons = np.array([350.0, 352.0, 354.0, 356.0])
@@ -1082,8 +1158,8 @@ def test_unrotate_1d(self):
         solx = np.array([-16.42176094, -14.85892262, -12.88946157, -10.35078336])
         soly = np.array([46.00724251, 51.29188893, 56.55031485, 61.77015703])
 
-        self.assertArrayAlmostEqual(resx, solx)
-        self.assertArrayAlmostEqual(resy, soly)
+        _shared_utils.assert_array_almost_equal(resx, solx)
+        _shared_utils.assert_array_almost_equal(resy, soly)
 
     def test_rotate_nd(self):
         rlons = np.array([[350.0, 351.0], [352.0, 353.0]])
@@ -1095,8 +1171,8 @@ def test_rotate_nd(self):
         solx = np.array([[148.69672569, 149.24727087], [149.79067025, 150.31754368]])
         soly = np.array([[18.60905789, 23.67749384], [28.74419024, 33.8087963]])
 
-        self.assertArrayAlmostEqual(resx, solx)
-        self.assertArrayAlmostEqual(resy, soly)
+        _shared_utils.assert_array_almost_equal(resx, solx)
+        _shared_utils.assert_array_almost_equal(resy, soly)
 
     def test_rotate_1d(self):
         rlons = np.array([350.0, 351.0, 352.0, 353.0])
@@ -1110,20 +1186,25 @@ def test_rotate_1d(self):
         solx = np.array([148.69672569, 149.24727087, 149.79067025, 150.31754368])
         soly = np.array([18.60905789, 23.67749384, 28.74419024, 33.8087963])
 
-        self.assertArrayAlmostEqual(resx, solx)
-        self.assertArrayAlmostEqual(resy, soly)
+        _shared_utils.assert_array_almost_equal(resx, solx)
+        _shared_utils.assert_array_almost_equal(resy, soly)
 
 
-@tests.skip_data
-class TestAreaWeights(tests.IrisTest):
+@_shared_utils.skip_data
+class TestAreaWeights:
+    @pytest.fixture(autouse=True)
+    def _setup(self, request):
+        self.request = request
+
     def test_area_weights(self):
-        small_cube = iris.tests.stock.simple_pp()
+        small_cube = stock.simple_pp()
         # Get offset, subsampled region: small enough to test against literals
         small_cube = small_cube[10:, 35:]
         small_cube = small_cube[::8, ::8]
         small_cube = small_cube[:5, :4]
         # pre-check non-data properties
-        self.assertCML(
+        _shared_utils.assert_CML(
+            self.request,
             small_cube,
             ("analysis", "areaweights_original.cml"),
             checksum=False,
@@ -1143,83 +1224,85 @@ def test_area_weights(self):
             ],
             dtype=np.float64,
         )
-        self.assertArrayAllClose(area_weights, expected_results, rtol=1e-8)
+        _shared_utils.assert_array_all_close(area_weights, expected_results, rtol=1e-8)
 
         # Check there was no residual change
         small_cube.coord("latitude").bounds = None
         small_cube.coord("longitude").bounds = None
-        self.assertCML(
+        _shared_utils.assert_CML(
+            self.request,
             small_cube,
             ("analysis", "areaweights_original.cml"),
             checksum=False,
         )
 
 
-@tests.skip_data
-class TestAreaWeightGeneration(tests.IrisTest):
-    def setUp(self):
-        self.cube = iris.tests.stock.realistic_4d()
+@_shared_utils.skip_data
+class TestAreaWeightGeneration:
+    @pytest.fixture(autouse=True)
+    def _setup(self):
+        self.cube = stock.realistic_4d()
 
     def test_area_weights_std(self):
         # weights for stock 4d data
         weights = iris.analysis.cartography.area_weights(self.cube)
-        self.assertEqual(weights.shape, self.cube.shape)
+        assert weights.shape == self.cube.shape
 
     def test_area_weights_order(self):
         # weights for data with dimensions in a different order
         order = [3, 2, 1, 0]  # (lon, lat, level, time)
         self.cube.transpose(order)
         weights = iris.analysis.cartography.area_weights(self.cube)
-        self.assertEqual(weights.shape, self.cube.shape)
+        assert weights.shape == self.cube.shape
 
     def test_area_weights_non_adjacent(self):
         # weights for cube with non-adjacent latitude/longitude dimensions
         order = [0, 3, 1, 2]  # (time, lon, level, lat)
         self.cube.transpose(order)
         weights = iris.analysis.cartography.area_weights(self.cube)
-        self.assertEqual(weights.shape, self.cube.shape)
+        assert weights.shape == self.cube.shape
 
     def test_area_weights_scalar_latitude(self):
         # weights for cube with a scalar latitude dimension
         cube = self.cube[:, :, 0, :]
         weights = iris.analysis.cartography.area_weights(cube)
-        self.assertEqual(weights.shape, cube.shape)
+        assert weights.shape == cube.shape
 
     def test_area_weights_scalar_longitude(self):
         # weights for cube with a scalar longitude dimension
         cube = self.cube[:, :, :, 0]
         weights = iris.analysis.cartography.area_weights(cube)
-        self.assertEqual(weights.shape, cube.shape)
+        assert weights.shape == cube.shape
 
     def test_area_weights_scalar(self):
         # weights for cube with scalar latitude and longitude dimensions
         cube = self.cube[:, :, 0, 0]
         weights = iris.analysis.cartography.area_weights(cube)
-        self.assertEqual(weights.shape, cube.shape)
+        assert weights.shape == cube.shape
 
     def test_area_weights_singleton_latitude(self):
         # singleton (1-point) latitude dimension
         cube = self.cube[:, :, 0:1, :]
         weights = iris.analysis.cartography.area_weights(cube)
-        self.assertEqual(weights.shape, cube.shape)
+        assert weights.shape == cube.shape
 
     def test_area_weights_singleton_longitude(self):
         # singleton (1-point) longitude dimension
         cube = self.cube[:, :, :, 0:1]
         weights = iris.analysis.cartography.area_weights(cube)
-        self.assertEqual(weights.shape, cube.shape)
+        assert weights.shape == cube.shape
 
     def test_area_weights_singletons(self):
         # singleton (1-point) latitude and longitude dimensions
         cube = self.cube[:, :, 0:1, 0:1]
         weights = iris.analysis.cartography.area_weights(cube)
-        self.assertEqual(weights.shape, cube.shape)
+        assert weights.shape == cube.shape
 
     def test_area_weights_normalized(self):
         # normalized area weights must sum to one over lat/lon dimensions.
         weights = iris.analysis.cartography.area_weights(self.cube, normalize=True)
         sumweights = weights.sum(axis=3).sum(axis=2)  # sum over lon and lat
-        self.assertArrayAlmostEqual(sumweights, 1)
+        _shared_utils.assert_array_almost_equal(sumweights, 1)
 
     def test_area_weights_non_contiguous(self):
         # Slice the cube so that we have non-contiguous longitude
@@ -1228,23 +1311,24 @@ def test_area_weights_non_contiguous(self):
         cube = self.cube[..., ind]
         weights = iris.analysis.cartography.area_weights(cube)
         expected = iris.analysis.cartography.area_weights(self.cube)[..., ind]
-        self.assertArrayEqual(weights, expected)
+        _shared_utils.assert_array_equal(weights, expected)
 
     def test_area_weights_no_lon_bounds(self):
         self.cube.coord("grid_longitude").bounds = None
-        with self.assertRaises(ValueError):
+        with pytest.raises(ValueError):
             iris.analysis.cartography.area_weights(self.cube)
 
     def test_area_weights_no_lat_bounds(self):
         self.cube.coord("grid_latitude").bounds = None
-        with self.assertRaises(ValueError):
+        with pytest.raises(ValueError):
             iris.analysis.cartography.area_weights(self.cube)
 
 
-@tests.skip_data
-class TestLatitudeWeightGeneration(tests.IrisTest):
-    def setUp(self):
-        path = iris.tests.get_data_path(
+@_shared_utils.skip_data
+class TestLatitudeWeightGeneration:
+    @pytest.fixture(autouse=True)
+    def _setup(self):
+        path = _shared_utils.get_data_path(
             ["NetCDF", "rotated", "xyt", "small_rotPole_precipitation.nc"]
         )
         self.cube = iris.load_cube(path)
@@ -1274,52 +1358,58 @@ def test_cosine_latitude_weights_range(self):
         )
         cube.add_dim_coord(lat_coord, 0)
         weights = iris.analysis.cartography.cosine_latitude_weights(cube)
-        self.assertTrue(weights.max() <= 1)
-        self.assertTrue(weights.min() >= 0)
+        assert weights.max() <= 1
+        assert weights.min() >= 0
 
     def test_cosine_latitude_weights_0d(self):
         # 0d latitude dimension (scalar coordinate)
         weights = iris.analysis.cartography.cosine_latitude_weights(
             self.cube_dim_lat[:, 0, :]
         )
-        self.assertEqual(weights.shape, self.cube_dim_lat[:, 0, :].shape)
-        self.assertAlmostEqual(weights[0, 0], np.cos(np.deg2rad(self.lat1d[0])))
+        assert weights.shape == self.cube_dim_lat[:, 0, :].shape
+        assert weights[0, 0] == pytest.approx(np.cos(np.deg2rad(self.lat1d[0])))
 
     def test_cosine_latitude_weights_1d_singleton(self):
         # singleton (1-point) 1d latitude coordinate (time, lat, lon)
         cube = self.cube_dim_lat[:, 0:1, :]
         weights = iris.analysis.cartography.cosine_latitude_weights(cube)
-        self.assertEqual(weights.shape, cube.shape)
-        self.assertAlmostEqual(weights[0, 0, 0], np.cos(np.deg2rad(self.lat1d[0])))
+        assert weights.shape == cube.shape
+        assert weights[0, 0, 0] == pytest.approx(np.cos(np.deg2rad(self.lat1d[0])))
 
     def test_cosine_latitude_weights_1d(self):
         # 1d latitude coordinate (time, lat, lon)
         weights = iris.analysis.cartography.cosine_latitude_weights(self.cube_dim_lat)
-        self.assertEqual(weights.shape, self.cube.shape)
-        self.assertArrayAlmostEqual(weights[0, :, 0], np.cos(np.deg2rad(self.lat1d)))
+        assert weights.shape == self.cube.shape
+        _shared_utils.assert_array_almost_equal(
+            weights[0, :, 0], np.cos(np.deg2rad(self.lat1d))
+        )
 
     def test_cosine_latitude_weights_1d_latitude_first(self):
         # 1d latitude coordinate with latitude first (lat, time, lon)
         order = [1, 0, 2]  # (lat, time, lon)
         self.cube_dim_lat.transpose(order)
         weights = iris.analysis.cartography.cosine_latitude_weights(self.cube_dim_lat)
-        self.assertEqual(weights.shape, self.cube_dim_lat.shape)
-        self.assertArrayAlmostEqual(weights[:, 0, 0], np.cos(np.deg2rad(self.lat1d)))
+        assert weights.shape == self.cube_dim_lat.shape
+        _shared_utils.assert_array_almost_equal(
+            weights[:, 0, 0], np.cos(np.deg2rad(self.lat1d))
+        )
 
     def test_cosine_latitude_weights_1d_latitude_last(self):
         # 1d latitude coordinate with latitude last (time, lon, lat)
         order = [0, 2, 1]  # (time, lon, lat)
         self.cube_dim_lat.transpose(order)
         weights = iris.analysis.cartography.cosine_latitude_weights(self.cube_dim_lat)
-        self.assertEqual(weights.shape, self.cube_dim_lat.shape)
-        self.assertArrayAlmostEqual(weights[0, 0, :], np.cos(np.deg2rad(self.lat1d)))
+        assert weights.shape == self.cube_dim_lat.shape
+        _shared_utils.assert_array_almost_equal(
+            weights[0, 0, :], np.cos(np.deg2rad(self.lat1d))
+        )
 
     def test_cosine_latitude_weights_2d_singleton1(self):
         # 2d latitude coordinate with first dimension singleton
         cube = self.cube_aux_lat[:, 0:1, :]
         weights = iris.analysis.cartography.cosine_latitude_weights(cube)
-        self.assertEqual(weights.shape, cube.shape)
-        self.assertArrayAlmostEqual(
+        assert weights.shape == cube.shape
+        _shared_utils.assert_array_almost_equal(
             weights[0, :, :], np.cos(np.deg2rad(self.lat2d[0:1, :]))
         )
 
@@ -1327,8 +1417,8 @@ def test_cosine_latitude_weights_2d_singleton2(self):
         # 2d latitude coordinate with second dimension singleton
         cube = self.cube_aux_lat[:, :, 0:1]
         weights = iris.analysis.cartography.cosine_latitude_weights(cube)
-        self.assertEqual(weights.shape, cube.shape)
-        self.assertArrayAlmostEqual(
+        assert weights.shape == cube.shape
+        _shared_utils.assert_array_almost_equal(
             weights[0, :, :], np.cos(np.deg2rad(self.lat2d[:, 0:1]))
         )
 
@@ -1336,47 +1426,56 @@ def test_cosine_latitude_weights_2d_singleton3(self):
         # 2d latitude coordinate with both dimensions singleton
         cube = self.cube_aux_lat[:, 0:1, 0:1]
         weights = iris.analysis.cartography.cosine_latitude_weights(cube)
-        self.assertEqual(weights.shape, cube.shape)
-        self.assertArrayAlmostEqual(
+        assert weights.shape == cube.shape
+        _shared_utils.assert_array_almost_equal(
             weights[0, :, :], np.cos(np.deg2rad(self.lat2d[0:1, 0:1]))
         )
 
     def test_cosine_latitude_weights_2d(self):
         # 2d latitude coordinate (time, lat, lon)
         weights = iris.analysis.cartography.cosine_latitude_weights(self.cube_aux_lat)
-        self.assertEqual(weights.shape, self.cube_aux_lat.shape)
-        self.assertArrayAlmostEqual(weights[0, :, :], np.cos(np.deg2rad(self.lat2d)))
+        assert weights.shape == self.cube_aux_lat.shape
+        _shared_utils.assert_array_almost_equal(
+            weights[0, :, :], np.cos(np.deg2rad(self.lat2d))
+        )
 
     def test_cosine_latitude_weights_2d_latitude_first(self):
         # 2d latitude coordinate with latitude first (lat, time, lon)
         order = [1, 0, 2]  # (lat, time, lon)
         self.cube_aux_lat.transpose(order)
         weights = iris.analysis.cartography.cosine_latitude_weights(self.cube_aux_lat)
-        self.assertEqual(weights.shape, self.cube_aux_lat.shape)
-        self.assertArrayAlmostEqual(weights[:, 0, :], np.cos(np.deg2rad(self.lat2d)))
+        assert weights.shape == self.cube_aux_lat.shape
+        _shared_utils.assert_array_almost_equal(
+            weights[:, 0, :], np.cos(np.deg2rad(self.lat2d))
+        )
 
     def test_cosine_latitude_weights_2d_latitude_last(self):
         # 2d latitude coordinate with latitude last (time, lon, lat)
         order = [0, 2, 1]  # (time, lon, lat)
         self.cube_aux_lat.transpose(order)
         weights = iris.analysis.cartography.cosine_latitude_weights(self.cube_aux_lat)
-        self.assertEqual(weights.shape, self.cube_aux_lat.shape)
-        self.assertArrayAlmostEqual(weights[0, :, :], np.cos(np.deg2rad(self.lat2d.T)))
+        assert weights.shape == self.cube_aux_lat.shape
+        _shared_utils.assert_array_almost_equal(
+            weights[0, :, :], np.cos(np.deg2rad(self.lat2d.T))
+        )
 
     def test_cosine_latitude_weights_no_latitude(self):
         # no coordinate identified as latitude
         self.cube_dim_lat.remove_coord("grid_latitude")
-        with self.assertRaises(ValueError):
+        with pytest.raises(ValueError):
             _ = iris.analysis.cartography.cosine_latitude_weights(self.cube_dim_lat)
 
     def test_cosine_latitude_weights_multiple_latitude(self):
         # two coordinates identified as latitude
-        with self.assertRaises(ValueError):
+        with pytest.raises(ValueError):
             _ = iris.analysis.cartography.cosine_latitude_weights(self.cube)
 
 
-class TestRollingWindow(tests.IrisTest):
-    def setUp(self):
+class TestRollingWindow:
+    @pytest.fixture(autouse=True)
+    def _setup(self, request):
+        self.request = request
+
         # XXX Comes from test_aggregated_by
         cube = iris.cube.Cube(
             np.array([[6, 10, 12, 18], [8, 12, 14, 20], [18, 12, 10, 6]]),
@@ -1407,7 +1506,7 @@ def test_non_mean_operator(self):
         expected_result = np.array(
             [[10, 12, 18], [12, 14, 20], [18, 12, 10]], dtype=np.float64
         )
-        self.assertArrayEqual(expected_result, res_cube.data)
+        _shared_utils.assert_array_equal(expected_result, res_cube.data)
 
     def test_longitude_simple(self):
         res_cube = self.cube.rolling_window("longitude", iris.analysis.MEAN, window=2)
@@ -1417,11 +1516,15 @@ def test_longitude_simple(self):
             dtype=np.float64,
         )
 
-        self.assertArrayEqual(expected_result, res_cube.data)
+        _shared_utils.assert_array_equal(expected_result, res_cube.data)
 
-        self.assertCML(res_cube, ("analysis", "rolling_window", "simple_longitude.cml"))
+        _shared_utils.assert_CML(
+            self.request,
+            res_cube,
+            ("analysis", "rolling_window", "simple_longitude.cml"),
+        )
 
-        self.assertRaises(
+        pytest.raises(
             ValueError,
             self.cube.rolling_window,
             "longitude",
@@ -1450,12 +1553,12 @@ def test_longitude_masked(self):
             dtype=np.float64,
         )
 
-        self.assertMaskedArrayEqual(expected_result, res_cube.data)
+        _shared_utils.assert_masked_array_almost_equal(expected_result, res_cube.data)
 
     def test_longitude_circular(self):
         cube = self.cube
         cube.coord("longitude").circular = True
-        self.assertRaises(
+        pytest.raises(
             iris.exceptions.NotYetImplementedError,
             self.cube.rolling_window,
             "longitude",
@@ -1468,12 +1571,16 @@ def test_different_length_windows(self):
 
         expected_result = np.array([[11.5], [13.5], [11.5]], dtype=np.float64)
 
-        self.assertArrayEqual(expected_result, res_cube.data)
+        _shared_utils.assert_array_equal(expected_result, res_cube.data)
 
-        self.assertCML(res_cube, ("analysis", "rolling_window", "size_4_longitude.cml"))
+        _shared_utils.assert_CML(
+            self.request,
+            res_cube,
+            ("analysis", "rolling_window", "size_4_longitude.cml"),
+        )
 
         # Window too long:
-        self.assertRaises(
+        pytest.raises(
             ValueError,
             self.cube.rolling_window,
             "longitude",
@@ -1481,7 +1588,7 @@ def test_different_length_windows(self):
             window=6,
         )
         # Window too small:
-        self.assertRaises(
+        pytest.raises(
             ValueError,
             self.cube.rolling_window,
             "longitude",
@@ -1490,7 +1597,7 @@ def test_different_length_windows(self):
         )
 
     def test_bad_coordinate(self):
-        self.assertRaises(
+        pytest.raises(
             KeyError,
             self.cube.rolling_window,
             "wibble",
@@ -1506,9 +1613,13 @@ def test_latitude_simple(self):
             dtype=np.float64,
         )
 
-        self.assertArrayEqual(expected_result, res_cube.data)
+        _shared_utils.assert_array_equal(expected_result, res_cube.data)
 
-        self.assertCML(res_cube, ("analysis", "rolling_window", "simple_latitude.cml"))
+        _shared_utils.assert_CML(
+            self.request,
+            res_cube,
+            ("analysis", "rolling_window", "simple_latitude.cml"),
+        )
 
     def test_mean_with_weights_consistency(self):
         # equal weights should be the same as the mean with no weights
@@ -1519,7 +1630,7 @@ def test_mean_with_weights_consistency(self):
         expected_result = self.cube.rolling_window(
             "longitude", iris.analysis.MEAN, window=2
         )
-        self.assertArrayEqual(expected_result.data, res_cube.data)
+        _shared_utils.assert_array_equal(expected_result.data, res_cube.data)
 
     def test_mean_with_weights(self):
         # rolling window mean with weights
@@ -1531,10 +1642,10 @@ def test_mean_with_weights(self):
             [[10.2, 13.6], [12.2, 15.6], [12.0, 9.0]], dtype=np.float64
         )
         # use almost equal to compare floats
-        self.assertArrayAlmostEqual(expected_result, res_cube.data)
+        _shared_utils.assert_array_almost_equal(expected_result, res_cube.data)
 
 
-class TestCreateWeightedAggregatorFn(tests.IrisTest):
+class TestCreateWeightedAggregatorFn:
     @staticmethod
     def aggregator_fn(data, axis, **kwargs):
         return (data, axis, kwargs)
@@ -1544,20 +1655,20 @@ def test_no_weights_supplied(self):
             self.aggregator_fn, 42, test_kwarg="test"
         )
         output = aggregator_fn("dummy_array", None)
-        self.assertEqual(len(output), 3)
-        self.assertEqual(output[0], "dummy_array")
-        self.assertEqual(output[1], 42)
-        self.assertEqual(output[2], {"test_kwarg": "test"})
+        assert len(output) == 3
+        assert output[0] == "dummy_array"
+        assert output[1] == 42
+        assert output[2] == {"test_kwarg": "test"}
 
     def test_weights_supplied(self):
         aggregator_fn = iris.analysis.create_weighted_aggregator_fn(
             self.aggregator_fn, 42, test_kwarg="test"
         )
         output = aggregator_fn("dummy_array", "w")
-        self.assertEqual(len(output), 3)
-        self.assertEqual(output[0], "dummy_array")
-        self.assertEqual(output[1], 42)
-        self.assertEqual(output[2], {"test_kwarg": "test", "weights": "w"})
+        assert len(output) == 3
+        assert output[0] == "dummy_array"
+        assert output[1] == 42
+        assert output[2] == {"test_kwarg": "test", "weights": "w"}
 
     def test_weights_in_kwargs(self):
         kwargs = {"test_kwarg": "test", "weights": "ignored"}
@@ -1565,16 +1676,16 @@ def test_weights_in_kwargs(self):
             self.aggregator_fn, 42, **kwargs
         )
         output = aggregator_fn("dummy_array", "w")
-        self.assertEqual(len(output), 3)
-        self.assertEqual(output[0], "dummy_array")
-        self.assertEqual(output[1], 42)
-        self.assertEqual(output[2], {"test_kwarg": "test", "weights": "w"})
-        self.assertEqual(kwargs, {"test_kwarg": "test", "weights": "ignored"})
+        assert len(output) == 3
+        assert output[0] == "dummy_array"
+        assert output[1] == 42
+        assert output[2] == {"test_kwarg": "test", "weights": "w"}
+        assert kwargs == {"test_kwarg": "test", "weights": "ignored"}
 
 
 class TestWeights:
     @pytest.fixture(autouse=True)
-    def setup_test_data(self):
+    def _setup_test_data(self):
         self.array_lib = np
         self.target_type = np.ndarray
         self.create_test_data()
@@ -1629,28 +1740,28 @@ def test_init_with_str_dim_coord(self):
         # DimCoord always realizes points
         assert isinstance(weights.array, np.ndarray)
         assert isinstance(weights.units, cf_units.Unit)
-        np.testing.assert_array_equal(weights.array, [[0, 0, 0], [1, 1, 1]])
+        _shared_utils.assert_array_equal(weights.array, [[0, 0, 0], [1, 1, 1]])
         assert weights.units == "degrees"
 
     def test_init_with_str_aux_coord(self):
         weights = _Weights("auxcoord", self.cube)
         assert isinstance(weights.array, self.target_type)
         assert isinstance(weights.units, cf_units.Unit)
-        np.testing.assert_array_equal(weights.array, [[3, 3, 3], [4, 4, 4]])
+        _shared_utils.assert_array_equal(weights.array, [[3, 3, 3], [4, 4, 4]])
         assert weights.units == "s"
 
     def test_init_with_str_ancillary_variable(self):
         weights = _Weights("ancvar", self.cube)
         assert isinstance(weights.array, self.target_type)
         assert isinstance(weights.units, cf_units.Unit)
-        np.testing.assert_array_equal(weights.array, [[5, 6, 7], [5, 6, 7]])
+        _shared_utils.assert_array_equal(weights.array, [[5, 6, 7], [5, 6, 7]])
         assert weights.units == "kg"
 
     def test_init_with_str_cell_measure(self):
         weights = _Weights("cell_area", self.cube)
         assert isinstance(weights.array, self.target_type)
         assert isinstance(weights.units, cf_units.Unit)
-        np.testing.assert_array_equal(weights.array, self.data)
+        _shared_utils.assert_array_equal(weights.array, self.data)
         assert weights.units == "m2"
 
     def test_init_with_dim_coord(self):
@@ -1658,28 +1769,28 @@ def test_init_with_dim_coord(self):
         # DimCoord always realizes points
         assert isinstance(weights.array, np.ndarray)
         assert isinstance(weights.units, cf_units.Unit)
-        np.testing.assert_array_equal(weights.array, [[0, 0, 0], [1, 1, 1]])
+        _shared_utils.assert_array_equal(weights.array, [[0, 0, 0], [1, 1, 1]])
         assert weights.units == "degrees"
 
     def test_init_with_aux_coord(self):
         weights = _Weights(self.aux_coord, self.cube)
         assert isinstance(weights.array, self.target_type)
         assert isinstance(weights.units, cf_units.Unit)
-        np.testing.assert_array_equal(weights.array, [[3, 3, 3], [4, 4, 4]])
+        _shared_utils.assert_array_equal(weights.array, [[3, 3, 3], [4, 4, 4]])
         assert weights.units == "s"
 
     def test_init_with_ancillary_variable(self):
         weights = _Weights(self.ancillary_variable, self.cube)
         assert isinstance(weights.array, self.target_type)
         assert isinstance(weights.units, cf_units.Unit)
-        np.testing.assert_array_equal(weights.array, [[5, 6, 7], [5, 6, 7]])
+        _shared_utils.assert_array_equal(weights.array, [[5, 6, 7], [5, 6, 7]])
         assert weights.units == "kg"
 
     def test_init_with_cell_measure(self):
         weights = _Weights(self.cell_measure, self.cube)
         assert isinstance(weights.array, self.target_type)
         assert isinstance(weights.units, cf_units.Unit)
-        np.testing.assert_array_equal(weights.array, self.data)
+        _shared_utils.assert_array_equal(weights.array, self.data)
         assert weights.units == "m2"
 
     def test_init_with_list(self):
@@ -1695,7 +1806,7 @@ class TestWeightsLazy(TestWeights):
     """Repeat tests from ``TestWeights`` with lazy arrays."""
 
     @pytest.fixture(autouse=True)
-    def setup_test_data(self):
+    def _setup_test_data(self):
         self.array_lib = da
         self.target_type = da.core.Array
         self.create_test_data()
@@ -1713,7 +1824,7 @@ def test__Groupby_repr():
 
 
 @pytest.mark.parametrize(
-    "kwargs,expected",
+    ("kwargs", "expected"),
     [
         ({}, "kg m-2"),
         ({"test": "m"}, "kg m-2"),
@@ -1742,7 +1853,3 @@ def test_sum_units_func(kwargs, expected):
     # changed if the units have not changed (even when weights units are "1")
     if result == units:
         assert result.origin == expected
-
-
-if __name__ == "__main__":
-    tests.main()
diff --git a/lib/iris/tests/test_analysis_calculus.py b/lib/iris/tests/test_analysis_calculus.py
index 70c1077def..74e0f90d8e 100644
--- a/lib/iris/tests/test_analysis_calculus.py
+++ b/lib/iris/tests/test_analysis_calculus.py
@@ -3,12 +3,8 @@
 # This file is part of Iris and is released under the BSD license.
 # See LICENSE in the root of the repository for full licensing details.
 
-# import iris tests first so that some things can be initialised before importing anything else
-import iris.tests as tests  # isort:skip
-
-import unittest
-
 import numpy as np
+import pytest
 
 import iris
 import iris.analysis.calculus
@@ -16,18 +12,19 @@
 import iris.coords
 from iris.coords import DimCoord
 import iris.cube
+from iris.tests import _shared_utils
 import iris.tests.stock
 
 
-class TestCubeDelta(tests.IrisTest):
-    @tests.skip_data
+class TestCubeDelta:
+    @_shared_utils.skip_data
     def test_invalid(self):
         cube = iris.tests.stock.realistic_4d()
-        with self.assertRaises(iris.exceptions.CoordinateMultiDimError):
+        with pytest.raises(iris.exceptions.CoordinateMultiDimError):
             _ = iris.analysis.calculus.cube_delta(cube, "surface_altitude")
-        with self.assertRaises(iris.exceptions.CoordinateMultiDimError):
+        with pytest.raises(iris.exceptions.CoordinateMultiDimError):
             _ = iris.analysis.calculus.cube_delta(cube, "altitude")
-        with self.assertRaises(ValueError):
+        with pytest.raises(ValueError):
             _ = iris.analysis.calculus.cube_delta(cube, "forecast_period")
 
     def test_delta_coord_lookup(self):
@@ -44,13 +41,13 @@ def test_delta_coord_lookup(self):
         cube.add_dim_coord(coord, 0)
         delta = iris.analysis.calculus.cube_delta(cube, "projection_x_coordinate")
         delta_coord = delta.coord("projection_x_coordinate")
-        self.assertEqual(delta_coord, delta.coord(coord))
-        self.assertEqual(coord, cube.coord(delta_coord))
+        assert delta_coord == delta.coord(coord)
+        assert coord == cube.coord(delta_coord)
 
 
-class TestDeltaAndMidpoint(tests.IrisTest):
+class TestDeltaAndMidpoint:
     def _simple_filename(self, suffix):
-        return tests.get_result_path(
+        return _shared_utils.get_result_path(
             ("analysis", "delta_and_midpoint", "simple%s.cml" % suffix)
         )
 
@@ -61,13 +58,13 @@ def test_simple1_delta_midpoint(self):
             units="degrees",
             circular=True,
         )
-        self.assertXMLElement(a, self._simple_filename("1"))
+        _shared_utils.assert_XML_element(a, self._simple_filename("1"))
 
         delta = iris.analysis.calculus._construct_delta_coord(a)
-        self.assertXMLElement(delta, self._simple_filename("1_delta"))
+        _shared_utils.assert_XML_element(delta, self._simple_filename("1_delta"))
 
         midpoint = iris.analysis.calculus._construct_midpoint_coord(a)
-        self.assertXMLElement(midpoint, self._simple_filename("1_midpoint"))
+        _shared_utils.assert_XML_element(midpoint, self._simple_filename("1_midpoint"))
 
     def test_simple2_delta_midpoint(self):
         a = iris.coords.DimCoord(
@@ -76,13 +73,13 @@ def test_simple2_delta_midpoint(self):
             units="degrees",
             circular=True,
         )
-        self.assertXMLElement(a, self._simple_filename("2"))
+        _shared_utils.assert_XML_element(a, self._simple_filename("2"))
 
         delta = iris.analysis.calculus._construct_delta_coord(a)
-        self.assertXMLElement(delta, self._simple_filename("2_delta"))
+        _shared_utils.assert_XML_element(delta, self._simple_filename("2_delta"))
 
         midpoint = iris.analysis.calculus._construct_midpoint_coord(a)
-        self.assertXMLElement(midpoint, self._simple_filename("2_midpoint"))
+        _shared_utils.assert_XML_element(midpoint, self._simple_filename("2_midpoint"))
 
     def test_simple3_delta_midpoint(self):
         a = iris.coords.DimCoord(
@@ -92,13 +89,13 @@ def test_simple3_delta_midpoint(self):
             circular=True,
         )
         a.guess_bounds(0.5)
-        self.assertXMLElement(a, self._simple_filename("3"))
+        _shared_utils.assert_XML_element(a, self._simple_filename("3"))
 
         delta = iris.analysis.calculus._construct_delta_coord(a)
-        self.assertXMLElement(delta, self._simple_filename("3_delta"))
+        _shared_utils.assert_XML_element(delta, self._simple_filename("3_delta"))
 
         midpoint = iris.analysis.calculus._construct_midpoint_coord(a)
-        self.assertXMLElement(midpoint, self._simple_filename("3_midpoint"))
+        _shared_utils.assert_XML_element(midpoint, self._simple_filename("3_midpoint"))
 
     def test_simple4_delta_midpoint(self):
         a = iris.coords.AuxCoord(
@@ -108,13 +105,13 @@ def test_simple4_delta_midpoint(self):
         )
         a.guess_bounds()
         b = a.copy()
-        self.assertXMLElement(b, self._simple_filename("4"))
+        _shared_utils.assert_XML_element(b, self._simple_filename("4"))
 
         delta = iris.analysis.calculus._construct_delta_coord(b)
-        self.assertXMLElement(delta, self._simple_filename("4_delta"))
+        _shared_utils.assert_XML_element(delta, self._simple_filename("4_delta"))
 
         midpoint = iris.analysis.calculus._construct_midpoint_coord(b)
-        self.assertXMLElement(midpoint, self._simple_filename("4_midpoint"))
+        _shared_utils.assert_XML_element(midpoint, self._simple_filename("4_midpoint"))
 
     def test_simple5_not_degrees_delta_midpoint(self):
         # Not sure it makes sense to have a circular coordinate which does not have a modulus but test it anyway.
@@ -124,13 +121,13 @@ def test_simple5_not_degrees_delta_midpoint(self):
             units="meter",
             circular=True,
         )
-        self.assertXMLElement(a, self._simple_filename("5"))
+        _shared_utils.assert_XML_element(a, self._simple_filename("5"))
 
         delta = iris.analysis.calculus._construct_delta_coord(a)
-        self.assertXMLElement(delta, self._simple_filename("5_delta"))
+        _shared_utils.assert_XML_element(delta, self._simple_filename("5_delta"))
 
         midpoints = iris.analysis.calculus._construct_midpoint_coord(a)
-        self.assertXMLElement(midpoints, self._simple_filename("5_midpoint"))
+        _shared_utils.assert_XML_element(midpoints, self._simple_filename("5_midpoint"))
 
     def test_simple6_delta_midpoint(self):
         a = iris.coords.DimCoord(
@@ -140,7 +137,7 @@ def test_simple6_delta_midpoint(self):
             circular=True,
         )
         midpoints = iris.analysis.calculus._construct_midpoint_coord(a)
-        self.assertXMLElement(midpoints, self._simple_filename("6"))
+        _shared_utils.assert_XML_element(midpoints, self._simple_filename("6"))
 
     def test_singular_delta(self):
         # Test single valued coordinate mid-points when circular
@@ -149,7 +146,7 @@ def test_singular_delta(self):
         )
 
         r_expl = iris.analysis.calculus._construct_delta_coord(lon)
-        self.assertXMLElement(
+        _shared_utils.assert_XML_element(
             r_expl,
             (
                 "analysis",
@@ -160,7 +157,7 @@ def test_singular_delta(self):
 
         # Test single valued coordinate mid-points when not circular
         lon.circular = False
-        with self.assertRaises(ValueError):
+        with pytest.raises(ValueError):
             iris.analysis.calculus._construct_delta_coord(lon)
 
     def test_singular_midpoint(self):
@@ -170,7 +167,7 @@ def test_singular_midpoint(self):
         )
 
         r_expl = iris.analysis.calculus._construct_midpoint_coord(lon)
-        self.assertXMLElement(
+        _shared_utils.assert_XML_element(
             r_expl,
             (
                 "analysis",
@@ -181,12 +178,13 @@ def test_singular_midpoint(self):
 
         # Test single valued coordinate mid-points when not circular
         lon.circular = False
-        with self.assertRaises(ValueError):
+        with pytest.raises(ValueError):
             iris.analysis.calculus._construct_midpoint_coord(lon)
 
 
-class TestCoordTrig(tests.IrisTest):
-    def setUp(self):
+class TestCoordTrig:
+    @pytest.fixture(autouse=True)
+    def _setup(self):
         points = np.arange(20, dtype=np.float32) * 2.3
         bounds = np.concatenate([[points - 0.5 * 2.3], [points + 0.5 * 2.3]]).T
         self.lat = iris.coords.AuxCoord(
@@ -204,41 +202,41 @@ def test_sin(self):
         sin_of_coord_radians = iris.analysis.calculus._coord_sin(self.rlat)
 
         # Check the values are correct (within a tolerance)
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             np.sin(self.rlat.points), sin_of_coord.points
         )
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             np.sin(self.rlat.bounds), sin_of_coord.bounds
         )
 
         # Check that the results of the sin function are almost equal when operating on a coord with degrees and radians
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             sin_of_coord.points, sin_of_coord_radians.points
         )
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             sin_of_coord.bounds, sin_of_coord_radians.bounds
         )
 
-        self.assertEqual(sin_of_coord.name(), "sin(latitude)")
-        self.assertEqual(sin_of_coord.units, "1")
+        assert sin_of_coord.name() == "sin(latitude)"
+        assert sin_of_coord.units == "1"
 
     def test_cos(self):
         cos_of_coord = iris.analysis.calculus._coord_cos(self.lat)
         cos_of_coord_radians = iris.analysis.calculus._coord_cos(self.rlat)
 
         # Check the values are correct (within a tolerance)
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             np.cos(self.rlat.points), cos_of_coord.points
         )
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             np.cos(self.rlat.bounds), cos_of_coord.bounds
         )
 
         # Check that the results of the cos function are almost equal when operating on a coord with degrees and radians
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             cos_of_coord.points, cos_of_coord_radians.points
         )
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             cos_of_coord.bounds, cos_of_coord_radians.bounds
         )
 
@@ -248,15 +246,19 @@ def test_cos(self):
             points=np.array([1], dtype=np.float32)
         )
 
-        self.assertXMLElement(cos_of_coord, ("analysis", "calculus", "cos_simple.xml"))
-        self.assertXMLElement(
+        _shared_utils.assert_XML_element(
+            cos_of_coord, ("analysis", "calculus", "cos_simple.xml")
+        )
+        _shared_utils.assert_XML_element(
             cos_of_coord_radians,
             ("analysis", "calculus", "cos_simple_radians.xml"),
         )
 
 
-class TestCalculusSimple3(tests.IrisTest):
-    def setUp(self):
+class TestCalculusSimple3:
+    @pytest.fixture(autouse=True)
+    def _setup(self, request):
+        self.request = request
         data = np.arange(2500, dtype=np.float32).reshape(50, 50)
         cube = iris.cube.Cube(data, standard_name="x_wind", units="km/h")
 
@@ -285,15 +287,21 @@ def setUp(self):
     def test_diff_wrt_lon(self):
         t = iris.analysis.calculus.differentiate(self.cube, "longitude")
 
-        self.assertCMLApproxData(t, ("analysis", "calculus", "handmade2_wrt_lon.cml"))
+        _shared_utils.assert_CML_approx_data(
+            self.request, t, ("analysis", "calculus", "handmade2_wrt_lon.cml")
+        )
 
     def test_diff_wrt_lat(self):
         t = iris.analysis.calculus.differentiate(self.cube, "latitude")
-        self.assertCMLApproxData(t, ("analysis", "calculus", "handmade2_wrt_lat.cml"))
+        _shared_utils.assert_CML_approx_data(
+            self.request, t, ("analysis", "calculus", "handmade2_wrt_lat.cml")
+        )
 
 
-class TestCalculusSimple2(tests.IrisTest):
-    def setUp(self):
+class TestCalculusSimple2:
+    @pytest.fixture(autouse=True)
+    def _setup(self, request):
+        self.request = request
         data = np.array(
             [
                 [1, 2, 3, 4, 5],
@@ -344,47 +352,57 @@ def setUp(self):
 
     def test_diff_wrt_x(self):
         t = iris.analysis.calculus.differentiate(self.cube, "x")
-        self.assertCMLApproxData(t, ("analysis", "calculus", "handmade_wrt_x.cml"))
+        _shared_utils.assert_CML_approx_data(
+            self.request, t, ("analysis", "calculus", "handmade_wrt_x.cml")
+        )
 
     def test_diff_wrt_y(self):
         t = iris.analysis.calculus.differentiate(self.cube, "y")
-        self.assertCMLApproxData(t, ("analysis", "calculus", "handmade_wrt_y.cml"))
+        _shared_utils.assert_CML_approx_data(
+            self.request, t, ("analysis", "calculus", "handmade_wrt_y.cml")
+        )
 
     def test_diff_wrt_lon(self):
         t = iris.analysis.calculus.differentiate(self.cube, "longitude")
-        self.assertCMLApproxData(t, ("analysis", "calculus", "handmade_wrt_lon.cml"))
+        _shared_utils.assert_CML_approx_data(
+            self.request, t, ("analysis", "calculus", "handmade_wrt_lon.cml")
+        )
 
     def test_diff_wrt_lat(self):
         t = iris.analysis.calculus.differentiate(self.cube, "latitude")
-        self.assertCMLApproxData(t, ("analysis", "calculus", "handmade_wrt_lat.cml"))
+        _shared_utils.assert_CML_approx_data(
+            self.request, t, ("analysis", "calculus", "handmade_wrt_lat.cml")
+        )
 
     def test_delta_wrt_x(self):
         t = iris.analysis.calculus.cube_delta(self.cube, "x")
-        self.assertCMLApproxData(
-            t, ("analysis", "calculus", "delta_handmade_wrt_x.cml")
+        _shared_utils.assert_CML_approx_data(
+            self.request, t, ("analysis", "calculus", "delta_handmade_wrt_x.cml")
         )
 
     def test_delta_wrt_y(self):
         t = iris.analysis.calculus.cube_delta(self.cube, "y")
-        self.assertCMLApproxData(
-            t, ("analysis", "calculus", "delta_handmade_wrt_y.cml")
+        _shared_utils.assert_CML_approx_data(
+            self.request, t, ("analysis", "calculus", "delta_handmade_wrt_y.cml")
         )
 
     def test_delta_wrt_lon(self):
         t = iris.analysis.calculus.cube_delta(self.cube, "longitude")
-        self.assertCMLApproxData(
-            t, ("analysis", "calculus", "delta_handmade_wrt_lon.cml")
+        _shared_utils.assert_CML_approx_data(
+            self.request, t, ("analysis", "calculus", "delta_handmade_wrt_lon.cml")
         )
 
     def test_delta_wrt_lat(self):
         t = iris.analysis.calculus.cube_delta(self.cube, "latitude")
-        self.assertCMLApproxData(
-            t, ("analysis", "calculus", "delta_handmade_wrt_lat.cml")
+        _shared_utils.assert_CML_approx_data(
+            self.request, t, ("analysis", "calculus", "delta_handmade_wrt_lat.cml")
         )
 
 
-class TestCalculusSimple1(tests.IrisTest):
-    def setUp(self):
+class TestCalculusSimple1:
+    @pytest.fixture(autouse=True)
+    def _setup(self, request):
+        self.request = request
         data = np.array(
             [
                 [1, 2, 3, 4, 5],
@@ -410,14 +428,14 @@ def setUp(self):
 
     def test_diff_wrt_x(self):
         t = iris.analysis.calculus.differentiate(self.cube, "x")
-        self.assertCMLApproxData(
-            t, ("analysis", "calculus", "handmade_simple_wrt_x.cml")
+        _shared_utils.assert_CML_approx_data(
+            self.request, t, ("analysis", "calculus", "handmade_simple_wrt_x.cml")
         )
 
     def test_delta_wrt_x(self):
         t = iris.analysis.calculus.cube_delta(self.cube, "x")
-        self.assertCMLApproxData(
-            t, ("analysis", "calculus", "delta_handmade_simple_wrt_x.cml")
+        _shared_utils.assert_CML_approx_data(
+            self.request, t, ("analysis", "calculus", "delta_handmade_simple_wrt_x.cml")
         )
 
 
@@ -501,7 +519,11 @@ def build_cube(data, spherical=False):
     return cube
 
 
-class TestCalculusWKnownSolutions(tests.IrisTest):
+class TestCalculusWKnownSolutions:
+    @pytest.fixture(autouse=True)
+    def _setup(self, request):
+        self.request = request
+
     def get_coord_pts(self, cube):
         """Return (x_pts, x_ones, y_pts, y_ones, z_pts, z_ones) for the given cube."""
         x = cube.coord(axis="X")
@@ -568,7 +590,7 @@ def test_contrived_differential1(self):
         data = -sin_x_pts * y_ones
         result = df_dlon.copy(data=data)
 
-        np.testing.assert_array_almost_equal(result.data, df_dlon.data, decimal=3)
+        _shared_utils.assert_array_almost_equal(result.data, df_dlon.data, decimal=3)
 
     def test_contrived_differential2(self):
         # testing :
@@ -585,7 +607,7 @@ def test_contrived_differential2(self):
         x_pts, x_ones, y_pts, y_ones, z_pts, z_ones = self.get_coord_pts(r)
         result = r.copy(data=y_pts * 2.0 * x_ones * z_ones)
 
-        np.testing.assert_array_almost_equal(result.data, r.data, decimal=6)
+        _shared_utils.assert_array_almost_equal(result.data, r.data, decimal=6)
 
     def test_contrived_non_spherical_curl1(self):
         # testing :
@@ -604,15 +626,16 @@ def test_contrived_non_spherical_curl1(self):
         r = iris.analysis.calculus.curl(u, v)
 
         # Curl returns None when there is no components of Curl
-        self.assertEqual(r[0], None)
-        self.assertEqual(r[1], None)
+        assert r[0] is None
+        assert r[1] is None
         cube = r[2]
-        self.assertCML(
+        _shared_utils.assert_CML(
+            self.request,
             cube,
             ("analysis", "calculus", "grad_contrived_non_spherical1.cml"),
             checksum=False,
         )
-        self.assertTrue(np.all(np.abs(cube.data - (-1.0)) < 1.0e-7))
+        assert np.all(np.abs(cube.data - (-1.0)) < 1.0e-7)
 
     def test_contrived_non_spherical_curl2(self):
         # testing :
@@ -639,18 +662,19 @@ def test_contrived_non_spherical_curl2(self):
         #        result.data = y_pts * 2. * x_ones * z_ones
         #        print(repr(r[0].data[0:1, 0:5, 0:25:5]))
         #        print(repr(result.data[0:1, 0:5, 0:25:5]))
-        #        np.testing.assert_array_almost_equal(result.data, r[0].data, decimal=2)
+        #        _shared_utils.assert_array_almost_equal(result.data, r[0].data, decimal=2)
         #
         #        result = r[1].copy(data=True)
         #        x_pts, x_ones, y_pts, y_ones, z_pts, z_ones = self.get_coord_pts(result)
         #        result.data = pow(z_pts, 2) * x_ones * y_ones
-        #        np.testing.assert_array_almost_equal(result.data, r[1].data, decimal=6)
+        #        _shared_utils.assert_array_almost_equal(result.data, r[1].data, decimal=6)
 
         result = r[2].copy()
         result.data = result.data * 0 + 1
-        np.testing.assert_array_almost_equal(result.data, r[2].data, decimal=4)
+        _shared_utils.assert_array_almost_equal(result.data, r[2].data, decimal=4)
 
-        self.assertCML(
+        _shared_utils.assert_CML(
+            self.request,
             r,
             ("analysis", "calculus", "curl_contrived_cartesian2.cml"),
             checksum=False,
@@ -681,11 +705,14 @@ def test_contrived_spherical_curl1(self):
         result = r.copy(data=r.data * 0)
 
         # Note: This numerical comparison was created when the radius was 1000 times smaller
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             result.data[5:-5], r.data[5:-5] / 1000.0, decimal=1
         )
-        self.assertCML(
-            r, ("analysis", "calculus", "grad_contrived1.cml"), checksum=False
+        _shared_utils.assert_CML(
+            self.request,
+            r,
+            ("analysis", "calculus", "grad_contrived1.cml"),
+            checksum=False,
         )
 
     def test_contrived_spherical_curl2(self):
@@ -730,22 +757,25 @@ def test_contrived_spherical_curl2(self):
         result = r.copy(data=-2 * cos_x_pts * cos_y_pts)
 
         # Note: This numerical comparison was created when the radius was 1000 times smaller
-        np.testing.assert_array_almost_equal(
+        _shared_utils.assert_array_almost_equal(
             result.data[30:-30, :], r.data[30:-30, :] / 1000.0, decimal=1
         )
-        self.assertCML(
-            r, ("analysis", "calculus", "grad_contrived2.cml"), checksum=False
+        _shared_utils.assert_CML(
+            self.request,
+            r,
+            ("analysis", "calculus", "grad_contrived2.cml"),
+            checksum=False,
         )
 
 
-class TestCurlInterface(tests.IrisTest):
+class TestCurlInterface:
     def test_non_conformed(self):
         u = build_cube(np.empty((50, 20)), spherical=True)
 
         v = u.copy()
         y = v.coord("latitude")
         y.points = y.points + 5
-        self.assertRaises(ValueError, iris.analysis.calculus.curl, u, v)
+        pytest.raises(ValueError, iris.analysis.calculus.curl, u, v)
 
     def test_standard_name(self):
         nx = 20
@@ -758,26 +788,26 @@ def test_standard_name(self):
         w.rename("w_wind")
 
         r = iris.analysis.calculus.spatial_vectors_with_phenom_name(u, v)
-        self.assertEqual(r, (("u", "v", "w"), "wind"))
+        assert r == (("u", "v", "w"), "wind")
 
         r = iris.analysis.calculus.spatial_vectors_with_phenom_name(u, v, w)
-        self.assertEqual(r, (("u", "v", "w"), "wind"))
+        assert r == (("u", "v", "w"), "wind")
 
-        self.assertRaises(
+        pytest.raises(
             ValueError,
             iris.analysis.calculus.spatial_vectors_with_phenom_name,
             u,
             None,
             w,
         )
-        self.assertRaises(
+        pytest.raises(
             ValueError,
             iris.analysis.calculus.spatial_vectors_with_phenom_name,
             None,
             None,
             w,
         )
-        self.assertRaises(
+        pytest.raises(
             ValueError,
             iris.analysis.calculus.spatial_vectors_with_phenom_name,
             None,
@@ -789,22 +819,22 @@ def test_standard_name(self):
         v.rename("y foobar wibble")
         w.rename("z foobar wibble")
         r = iris.analysis.calculus.spatial_vectors_with_phenom_name(u, v)
-        self.assertEqual(r, (("x", "y", "z"), "foobar wibble"))
+        assert r == (("x", "y", "z"), "foobar wibble")
 
         r = iris.analysis.calculus.spatial_vectors_with_phenom_name(u, v, w)
-        self.assertEqual(r, (("x", "y", "z"), "foobar wibble"))
+        assert r == (("x", "y", "z"), "foobar wibble")
 
         u.rename("wibble foobar")
         v.rename("wobble foobar")
         w.rename("tipple foobar")
         #        r = iris.analysis.calculus.spatial_vectors_with_phenom_name(u, v, w) #should raise a Value Error...
-        self.assertRaises(
+        pytest.raises(
             ValueError,
             iris.analysis.calculus.spatial_vectors_with_phenom_name,
             u,
             v,
         )
-        self.assertRaises(
+        pytest.raises(
             ValueError,
             iris.analysis.calculus.spatial_vectors_with_phenom_name,
             u,
@@ -816,14 +846,14 @@ def test_standard_name(self):
         v.rename("northward_foobar")
         w.rename("upward_foobar")
         r = iris.analysis.calculus.spatial_vectors_with_phenom_name(u, v)
-        self.assertEqual(r, (("eastward", "northward", "upward"), "foobar"))
+        assert r == (("eastward", "northward", "upward"), "foobar")
 
         r = iris.analysis.calculus.spatial_vectors_with_phenom_name(u, v, w)
-        self.assertEqual(r, (("eastward", "northward", "upward"), "foobar"))
+        assert r == (("eastward", "northward", "upward"), "foobar")
 
         # Change it to have an inconsistent phenomenon
         v.rename("northward_foobar2")
-        self.assertRaises(
+        pytest.raises(
             ValueError,
             iris.analysis.calculus.spatial_vectors_with_phenom_name,
             u,
@@ -837,8 +867,4 @@ def test_rotated_pole(self):
         v.rename("v_wind")
 
         x, y, z = iris.analysis.calculus.curl(u, v)
-        self.assertEqual(z.coord_system(), u.coord_system())
-
-
-if __name__ == "__main__":
-    unittest.main()
+        assert z.coord_system() == u.coord_system()