Repository: ESMValGroup/ESMValCore

Issue:
Use iris' regridder caching for faster regridding?
Some regridding schemes of iris allow reusing the regridding weights (see [here ](https://scitools-iris.readthedocs.io/en/latest/further_topics/which_regridder_to_use.html)for an overview of them), which might significantly speed up our `regrid` preprocessor if multiple variables of the same dataset are regridded. This is described [here](https://scitools-iris.readthedocs.io/en/latest/userguide/interpolation_and_regridding.html#caching-a-regridder) and basically boils down to using the following pattern

```py
regridder = iris.analysis.Nearest().regridder(source_cube, target_cube)
for cube in list_of_cubes_on_source_grid:
    result = regridder(cube)
```

instead of

```py
for cube in list_of_cubes_on_source_grid:
    result = cube.regrid(target_cube, iris.analysis.Nearest())
```

This should be fairly easy to implement with a `dict` cache that uses the source and target latitude and longitude dimensions as keys (we cannot use `lru_cache` since the `hash` of a coordinate is simply its `id`, therefore this would never use cached results for different cubes).

One drawback here is that we have an additional comparison of the coordinates for each `regrid` call, which makes it slightly slower if just a single variable for many data sets is analyzed.

@ESMValGroup/esmvaltool-coreteam is this something we want? I can open a draft PR so that the discussion is easier.

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2344?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 94.12%. Comparing base [(`bb7866e`)](https://app.codecov.io/gh/ESMValGroup/ESMValCore/commit/bb7866e6a136f6e011bc6419c4e3ed9eec25721e?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) to head [(`c257945`)](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2344?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup).


<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #2344      +/-   ##
==========================================
+ Coverage   94.11%   94.12%   +0.01%     
==========================================
  Files         241      241              
  Lines       13376    13399      +23     
==========================================
+ Hits        12589    12612      +23     
  Misses        787      787              
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2344?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup).

schlunma: Thanks for reviewing V! Memory usage should be minimal: as mentioned [here](https://github.com/ESMValGroup/ESMValCore/issues/2341#issuecomment-1961400007), the weights should only be around 10 MiB for very high resolution grids (1000x1000), and much smaller for "normal" resolutions.
valeriupredoi: @bouweandela maybe you could have a look too, given that I took you are not 100% convinced about the need of caching (reading from the issue) - and pls merge if all's good by ye too :beer: 
bouweandela: Yes, a 10 percent increase in runtime in the best case doesn't seem like a huge gain, but it's nice to have of course. I am concerned about the size of the cache though, in the issue you talked about 1GB https://github.com/ESMValGroup/ESMValCore/issues/2341#issuecomment-1961400007 and I saw that there is some discussion planned at the workshop about re-using weights:
> High resolution model data: often the weights to interpolate grids can be reused, not only for different variables but even across different experiments, time periods, and simulations of the same model. Could ESMValTool support the reading and use of precalculated weights because their calculation is very time consuming? (@katjaweigel?)

I image they would be rather large arrays too if they are so expensive to compute.

To avoid this turning into a memory leak, would it be an option to:
1) add a method to the `regrid` preprocessor function to clear the cache, similar to how it's done if you're using lru_cache, so Python users can clear the cache when they are done with it?
1) make weights caching optional by adding a `cache_weights` or similar argument to the preprocessor function?
schlunma: > Yes, a 10 percent increase in runtime in the best case doesn't seem like a huge gain, but it's nice to have of course. I am concerned about the size of the cache though, in the issue you talked about 1GB [#2341 (comment)](https://github.com/ESMValGroup/ESMValCore/issues/2341#issuecomment-1961400007)

The 10% reduction in run time were only for this specific example, it's definitely not an upper bound. Also, the 1GiB is a really extreme example (e.g., a 0.1x0.1° grid would lead to a weights array of ~25 MiB).

> To avoid this turning into a memory leak, would it be an option to:
> 
> 1. add a method to the `regrid` preprocessor function to clear the cache, similar to how it's done if you're using lru_cache, so Python users can clear the cache when they are done with it?
> 
> 2. make weights caching optional by adding a `cache_weights` or similar argument to the preprocessor function, so it can be enabled from the recipe if it is relevant?

That sounds very reasonable, I will do that! 👍 
valeriupredoi: from me own experience it's always better to read from file/memory if possible, than it is to compute (yes OK call me Dr Obvious :rofl: ) - and given Manu's reassuring info on upper limits for mem intake, I think it's a go - but I do like Bouwe's suggestions too. My only concern, that I just thought of, is maybe this is better it sat in iris?
valeriupredoi: Also (sorry, I sat down did MO crap today, so now my brain is free at last) - Manu, maybe it'd be worth ploppoing a worst case scenario test in the tests ie build a super-high res netCDF file on the fly, and do the weights caching dance on it - we'd be able to monitor the memory using the test performance tool we have in Github Actions, and of course, decorate it with eg `@highmem` so we don't run it usually :beer: 
schlunma: > My only concern, that I just thought of, is maybe this is better it sat in iris?

Actually this entire functionality [is already in iris](https://scitools-iris.readthedocs.io/en/stable/userguide/interpolation_and_regridding.html#caching-a-regridder) in the form of `Regridder` classes that compute the weights and store them. With this PR, we are using those `Regridder` classes instead of the regridding *schemes* which do not offer weights caching.

> Also (sorry, I sat down did MO crap today, so now my brain is free at last) - Manu, maybe it'd be worth ploppoing a worst case scenario test in the tests ie build a super-high res netCDF file on the fly, and do the weights caching dance on it - we'd be able to monitor the memory using the test performance tool we have in Github Actions, and of course, decorate it with eg @highmem so we don't run it usually 🍺

I am honestly not sure if that's worth it, especially if we implement the option to turn off weights caching 😬 What do we want to achieve with such a test? Kill the CI machines? 😄 🤖 
schlunma: Added option to enable/disable (default: disabled) caching weights: https://github.com/ESMValGroup/ESMValCore/pull/2344/commits/b6eaade2d39c237136d8aa2da68333b1494da913

Added function to clear cache similar to `lru.cache_clear()`: `regrid.cache_clear()`, see https://github.com/ESMValGroup/ESMValCore/pull/2344/commits/52a536ccd0e14f56820ba745816c607e22015420
bouweandela: Thanks for making the changes @schlunma! @valeriupredoi Could you please do a final review to review the new code, since you reviewed the original PR, and then we can merge.
valeriupredoi: On it, bud 🍺 
valeriupredoi: > What do we want to achieve with such a test? Kill the CI machines?

Not a smart thing to do, especially since today Skynet's closer than ever :laughing: OK, leave the poor machines alone then

Patch:
diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst
index cd3566498..4d7ef23c5 100644
--- a/doc/recipe/preprocessor.rst
+++ b/doc/recipe/preprocessor.rst
@@ -780,10 +780,6 @@ regridding is based on the horizontal grid of another cube (the reference
 grid). If the horizontal grids of a cube and its reference grid are sufficiently
 the same, regridding is automatically and silently skipped for performance reasons.
 
-The underlying regridding mechanism in ESMValCore uses
-:obj:`iris.cube.Cube.regrid`
-from Iris.
-
 The use of the horizontal regridding functionality is flexible depending on
 what type of reference grid and what interpolation scheme is preferred. Below
 we show a few examples.
@@ -821,7 +817,7 @@ cell specification is oftentimes used when operating on localized data.
           target_grid: 2.5x2.5
           scheme: nearest
 
-In this case the ``NearestNeighbour`` interpolation scheme is used (see below
+In this case the nearest-neighbor interpolation scheme is used (see below
 for scheme definitions).
 
 When using a ``MxN`` type of grid it is possible to offset the grid cell
@@ -917,9 +913,6 @@ Built-in regridding schemes
   :class:`~esmvalcore.preprocessor.regrid_schemes.ESMPyAreaWeighted`.
   Source data on an unstructured grid is not supported, yet.
 
-See also :func:`esmvalcore.preprocessor.regrid`
-
-
 .. _generic regridding schemes:
 
 Generic regridding schemes
@@ -1017,6 +1010,37 @@ scheme available in :doc:`iris-esmf-regrid:index`:
           reference: esmf_regrid.schemes:regrid_rectilinear_to_rectilinear
           mdtol: 0.7
 
+.. _caching_regridding_weights:
+
+Reusing regridding weights
+--------------------------
+
+If desired, regridding weights can be cached to reduce run times (see `here
+<https://scitools-iris.readthedocs.io/en/latest/userguide/interpolation_and_regridding.html#caching-a-regridder>`__
+for technical details on this).
+This can speed up the regridding of different datasets with similar source and
+target grids massively, but may take up a lot of memory for extremely
+high-resolution data.
+By default, this feature is disabled; to enable it, use the option
+``cache_weights: true`` in the preprocessor definition:
+
+.. code-block:: yaml
+
+    preprocessors:
+      regrid_preprocessor:
+        regrid:
+          target_grid: 0.1x0.1
+          scheme: linear
+          cache_weights: true
+
+Not all regridding schemes support weights caching. An overview of those that
+do is given `here
+<https://scitools-iris.readthedocs.io/en/latest/further_topics/which_regridder_to_use.html#which-regridder-to-use>`__
+and in the docstrings :ref:`here <regridding_schemes>`.
+
+See also :func:`esmvalcore.preprocessor.regrid`
+
+
 .. _ensemble statistics:
 
 Ensemble statistics
diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py
index 4c84cedfc..c9bc0915d 100644
--- a/esmvalcore/preprocessor/_regrid.py
+++ b/esmvalcore/preprocessor/_regrid.py
@@ -500,8 +500,7 @@ def _get_target_grid_cube(
     elif isinstance(target_grid, (str, Path)) and os.path.isfile(target_grid):
         target_grid_cube = iris.load_cube(target_grid)
     elif isinstance(target_grid, str):
-        # Generate a target grid from the provided cell-specification,
-        # and cache the resulting stock cube for later use.
+        # Generate a target grid from the provided cell-specification
         target_grid_cube = _global_stock_cube(
             target_grid, lat_offset, lon_offset
         )
@@ -639,12 +638,82 @@ def _load_generic_scheme(scheme: dict):
     return loaded_scheme
 
 
+_CACHED_REGRIDDERS: dict[tuple, dict] = {}
+
+
+def _get_regridder(
+    src_cube: Cube,
+    tgt_cube: Cube,
+    scheme: str | dict,
+    cache_weights: bool,
+):
+    """Get regridder to actually perform regridding.
+
+    Note
+    ----
+    If possible, this uses an existing regridder to reduce runtime (see also
+    https://scitools-iris.readthedocs.io/en/latest/userguide/
+    interpolation_and_regridding.html#caching-a-regridder.)
+
+    """
+    # (1) Weights caching enabled
+    if cache_weights:
+        # To search for a matching regridder in the cache, first check the
+        # regridding scheme name and shapes of source and target coordinates.
+        # Only if these match, check coordinates themselves (this is much more
+        # expensive).
+        coord_key = _get_coord_key(src_cube, tgt_cube)
+        name_shape_key = _get_name_and_shape_key(src_cube, tgt_cube, scheme)
+        if name_shape_key in _CACHED_REGRIDDERS:
+            # We cannot simply do a test for `coord_key in
+            # _CACHED_REGRIDDERS[shape_key]` below since the hash() of a
+            # coordinate is simply its id() (thus, coordinates loaded from two
+            # different files would never be considered equal)
+            for (key, regridder) in _CACHED_REGRIDDERS[name_shape_key].items():
+                if key == coord_key:
+                    return regridder
+
+        # Regridder is not in cached -> return a new one and cache it
+        loaded_scheme = _load_scheme(src_cube, scheme)
+        regridder = loaded_scheme.regridder(src_cube, tgt_cube)
+        _CACHED_REGRIDDERS.setdefault(name_shape_key, {})
+        _CACHED_REGRIDDERS[name_shape_key][coord_key] = regridder
+
+    # (2) Weights caching disabled
+    else:
+        loaded_scheme = _load_scheme(src_cube, scheme)
+        regridder = loaded_scheme.regridder(src_cube, tgt_cube)
+
+    return regridder
+
+
+def _get_coord_key(src_cube: Cube, tgt_cube: Cube) -> tuple:
+    """Get dict key from coordinates."""
+    src_lat = src_cube.coord('latitude')
+    src_lon = src_cube.coord('longitude')
+    tgt_lat = tgt_cube.coord('latitude')
+    tgt_lon = tgt_cube.coord('longitude')
+    return (src_lat, src_lon, tgt_lat, tgt_lon)
+
+
+def _get_name_and_shape_key(
+    src_cube: Cube,
+    tgt_cube: Cube,
+    scheme: str | dict,
+) -> tuple:
+    """Get dict key from scheme name and coordinate shapes."""
+    name = str(scheme)
+    shapes = [c.shape for c in _get_coord_key(src_cube, tgt_cube)]
+    return (name, *shapes)
+
+
 def regrid(
     cube: Cube,
     target_grid: Cube | Dataset | Path | str | dict,
     scheme: str | dict,
     lat_offset: bool = True,
     lon_offset: bool = True,
+    cache_weights: bool = False,
 ) -> Cube:
     """Perform horizontal regridding.
 
@@ -691,6 +760,14 @@ def regrid(
         Offset the grid centers of the longitude coordinate w.r.t. Greenwich
         meridian by half a grid step. This argument is ignored if
         `target_grid` is a cube or file.
+    cache_weights:
+        If ``True``, cache regridding weights for later usage. This can speed
+        up the regridding of different datasets with similar source and target
+        grids massively, but may take up a lot of memory for extremely
+        high-resolution data. This option is ignored for schemes that do not
+        support weights caching. More details on this are given in the section
+        on :ref:`caching_regridding_weights`. To clear the cache, use
+        :func:`esmvalcore.preprocessor.regrid.cache_clear`.
 
     Returns
     -------
@@ -757,16 +834,26 @@ def regrid(
             )
         return cube
 
-    # Load scheme, rechunk and regrid
+    # Load scheme and reuse existing regridder if possible
     if isinstance(scheme, str):
         scheme = scheme.lower()
-    loaded_scheme = _load_scheme(cube, scheme)
+    regridder = _get_regridder(cube, target_grid_cube, scheme, cache_weights)
+
+    # Rechunk and actually perform the regridding
     cube = _rechunk(cube, target_grid_cube)
-    cube = cube.regrid(target_grid_cube, loaded_scheme)
+    cube = regridder(cube)
 
     return cube
 
 
+def _cache_clear():
+    """Clear regridding weights cache."""
+    _CACHED_REGRIDDERS.clear()
+
+
+regrid.cache_clear = _cache_clear  # type: ignore
+
+
 def _rechunk(cube: Cube, target_grid: Cube) -> Cube:
     """Re-chunk cube with optimal chunk sizes for target grid."""
     if not cube.has_lazy_data() or cube.ndim < 3:
diff --git a/esmvalcore/preprocessor/_regrid_esmpy.py b/esmvalcore/preprocessor/_regrid_esmpy.py
index 7da6dd63c..dd47b08e9 100755
--- a/esmvalcore/preprocessor/_regrid_esmpy.py
+++ b/esmvalcore/preprocessor/_regrid_esmpy.py
@@ -43,6 +43,8 @@ MASK_REGRIDDING_MASK_VALUE = {
 class ESMPyRegridder:
     """General ESMPy regridder.
 
+    Does not support lazy regridding nor weights caching.
+
     Parameters
     ----------
     src_cube:
diff --git a/esmvalcore/preprocessor/regrid_schemes.py b/esmvalcore/preprocessor/regrid_schemes.py
index 13dee49a0..f84284fd4 100644
--- a/esmvalcore/preprocessor/regrid_schemes.py
+++ b/esmvalcore/preprocessor/regrid_schemes.py
@@ -31,6 +31,9 @@ __all__ = [
 class GenericRegridder:
     r"""Generic function regridder.
 
+    Does support lazy regridding if `func` does. Does not support weights
+    caching.
+
     Parameters
     ----------
     src_cube:


Test patch:
diff --git a/tests/integration/preprocessor/_regrid/test_regrid.py b/tests/integration/preprocessor/_regrid/test_regrid.py
index c774433a4..8c272a1e4 100644
--- a/tests/integration/preprocessor/_regrid/test_regrid.py
+++ b/tests/integration/preprocessor/_regrid/test_regrid.py
@@ -107,19 +107,29 @@ class Test:
             aux_coords_and_dims=[(lats_2d, (0, 1)), (lons_2d, (0, 1))],
         )
 
-    def test_regrid__linear(self):
-        result = regrid(self.cube, self.grid_for_linear, 'linear')
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_regrid__linear(self, cache_weights):
+        result = regrid(
+            self.cube,
+            self.grid_for_linear,
+            'linear',
+            cache_weights=cache_weights,
+        )
         expected = np.array([[[1.5]], [[5.5]], [[9.5]]])
         assert_array_equal(result.data, expected)
 
-    def test_regrid__linear_file(self, tmp_path):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_regrid__linear_file(self, tmp_path, cache_weights):
         file = tmp_path / "file.nc"
         iris.save(self.grid_for_linear, target=file)
-        result = regrid(self.cube, file, 'linear')
+        result = regrid(
+            self.cube, file, 'linear', cache_weights=cache_weights
+        )
         expected = np.array([[[1.5]], [[5.5]], [[9.5]]])
         assert_array_equal(result.data, expected)
 
-    def test_regrid__linear_dataset(self, monkeypatch):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_regrid__linear_dataset(self, monkeypatch, cache_weights):
         monkeypatch.setattr(Dataset, 'files', ["file.nc"])
 
         def load(_):
@@ -129,11 +139,14 @@ class Test:
         dataset = Dataset(
             short_name='tas',
         )
-        result = regrid(self.cube, dataset, 'linear')
+        result = regrid(
+            self.cube, dataset, 'linear', cache_weights=cache_weights
+        )
         expected = np.array([[[1.5]], [[5.5]], [[9.5]]])
         assert_array_equal(result.data, expected)
 
-    def test_regrid__esmf_rectilinear(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_regrid__esmf_rectilinear(self, cache_weights):
         scheme_name = 'esmf_regrid.schemes:regrid_rectilinear_to_rectilinear'
         scheme = {
             'reference': scheme_name
@@ -141,11 +154,14 @@ class Test:
         result = regrid(
             self.cube,
             self.grid_for_linear,
-            scheme)
+            scheme,
+            cache_weights=cache_weights,
+        )
         expected = np.array([[[1.5]], [[5.5]], [[9.5]]])
         np.testing.assert_array_almost_equal(result.data, expected, decimal=1)
 
-    def test_regrid__regular_coordinates(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_regrid__regular_coordinates(self, cache_weights):
         data = np.ones((1, 1))
         lons = iris.coords.DimCoord([1.50000000000001],
                                     standard_name='longitude',
@@ -159,18 +175,30 @@ class Test:
                                     coord_system=self.cs)
         coords_spec = [(lats, 0), (lons, 1)]
         regular_grid = iris.cube.Cube(data, dim_coords_and_dims=coords_spec)
-        result = regrid(regular_grid, self.grid_for_linear, 'linear')
+        result = regrid(
+            regular_grid,
+            self.grid_for_linear,
+            'linear',
+            cache_weights=cache_weights,
+        )
         iris.common.resolve.Resolve(result, self.grid_for_linear)
 
-    def test_regrid__linear_do_not_preserve_dtype(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_regrid__linear_do_not_preserve_dtype(self, cache_weights):
         self.cube.data = self.cube.data.astype(int)
-        result = regrid(self.cube, self.grid_for_linear, 'linear')
+        result = regrid(
+            self.cube,
+            self.grid_for_linear,
+            'linear',
+            cache_weights=cache_weights,
+        )
         expected = np.array([[[1.5]], [[5.5]], [[9.5]]])
         assert_array_equal(result.data, expected)
         assert np.issubdtype(self.cube.dtype, np.integer)
         assert np.issubdtype(result.dtype, np.floating)
 
-    def test_regrid__linear_with_extrapolation(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_regrid__linear_with_extrapolation(self, cache_weights):
         data = np.empty((3, 3))
         lons = iris.coords.DimCoord([0, 1.5, 3],
                                     standard_name='longitude',
@@ -188,13 +216,14 @@ class Test:
             'reference': 'iris.analysis:Linear',
             'extrapolation_mode': 'extrapolate',
         }
-        result = regrid(self.cube, grid, scheme)
+        result = regrid(self.cube, grid, scheme, cache_weights=cache_weights)
         expected = [[[-3., -1.5, 0.], [0., 1.5, 3.], [3., 4.5, 6.]],
                     [[1., 2.5, 4.], [4., 5.5, 7.], [7., 8.5, 10.]],
                     [[5., 6.5, 8.], [8., 9.5, 11.], [11., 12.5, 14.]]]
         assert_array_equal(result.data, expected)
 
-    def test_regrid__linear_with_mask(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_regrid__linear_with_mask(self, cache_weights):
         data = np.empty((3, 3))
         grid = iris.cube.Cube(data)
         lons = iris.coords.DimCoord([0, 1.5, 3],
@@ -209,13 +238,14 @@ class Test:
                                     coord_system=self.cs)
         coords_spec = [(lats, 0), (lons, 1)]
         grid = iris.cube.Cube(data, dim_coords_and_dims=coords_spec)
-        result = regrid(self.cube, grid, 'linear')
+        result = regrid(self.cube, grid, 'linear', cache_weights=cache_weights)
         expected = ma.empty((3, 3, 3))
         expected.mask = ma.masked
         expected[:, 1, 1] = np.array([1.5, 5.5, 9.5])
         assert_array_equal(result.data, expected)
 
-    def test_regrid__nearest(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_regrid__nearest(self, cache_weights):
         data = np.empty((1, 1))
         lons = iris.coords.DimCoord([1.6],
                                     standard_name='longitude',
@@ -229,11 +259,14 @@ class Test:
                                     coord_system=self.cs)
         coords_spec = [(lats, 0), (lons, 1)]
         grid = iris.cube.Cube(data, dim_coords_and_dims=coords_spec)
-        result = regrid(self.cube, grid, 'nearest')
+        result = regrid(
+            self.cube, grid, 'nearest', cache_weights=cache_weights
+        )
         expected = np.array([[[3]], [[7]], [[11]]])
         assert_array_equal(result.data, expected)
 
-    def test_regrid__nearest_extrapolate_with_mask(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_regrid__nearest_extrapolate_with_mask(self, cache_weights):
         data = np.empty((3, 3))
         lons = iris.coords.DimCoord([0, 1.6, 3],
                                     standard_name='longitude',
@@ -247,13 +280,16 @@ class Test:
                                     coord_system=self.cs)
         coords_spec = [(lats, 0), (lons, 1)]
         grid = iris.cube.Cube(data, dim_coords_and_dims=coords_spec)
-        result = regrid(self.cube, grid, 'nearest')
+        result = regrid(
+            self.cube, grid, 'nearest', cache_weights=cache_weights
+        )
         expected = ma.empty((3, 3, 3))
         expected.mask = ma.masked
         expected[:, 1, 1] = np.array([3, 7, 11])
         assert_array_equal(result.data, expected)
 
-    def test_regrid__area_weighted(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_regrid__area_weighted(self, cache_weights):
         data = np.empty((1, 1))
         lons = iris.coords.DimCoord([1.6],
                                     standard_name='longitude',
@@ -267,11 +303,14 @@ class Test:
                                     coord_system=self.cs)
         coords_spec = [(lats, 0), (lons, 1)]
         grid = iris.cube.Cube(data, dim_coords_and_dims=coords_spec)
-        result = regrid(self.cube, grid, 'area_weighted')
+        result = regrid(
+            self.cube, grid, 'area_weighted', cache_weights=cache_weights
+        )
         expected = np.array([1.499886, 5.499886, 9.499886])
         np.testing.assert_array_almost_equal(result.data, expected, decimal=6)
 
-    def test_regrid__esmf_area_weighted(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_regrid__esmf_area_weighted(self, cache_weights):
         data = np.empty((1, 1))
         lons = iris.coords.DimCoord([1.6],
                                     standard_name='longitude',
@@ -288,15 +327,19 @@ class Test:
         scheme = {
             'reference': 'esmf_regrid.schemes:ESMFAreaWeighted'
         }
-        result = regrid(self.cube, grid, scheme)
+        result = regrid(self.cube, grid, scheme, cache_weights=cache_weights)
         expected = np.array([[[1.499886]], [[5.499886]], [[9.499886]]])
         np.testing.assert_array_almost_equal(result.data, expected, decimal=6)
 
-    def test_regrid_nearest_unstructured_grid_float(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_regrid_nearest_unstructured_grid_float(self, cache_weights):
         """Test `nearest` regridding with unstructured cube of floats."""
-        result = regrid(self.unstructured_grid_cube,
-                        self.tgt_grid_for_unstructured,
-                        'nearest')
+        result = regrid(
+            self.unstructured_grid_cube,
+            self.tgt_grid_for_unstructured,
+            'nearest',
+            cache_weights=cache_weights,
+        )
         expected = np.ma.array([[[3.0]], [[7.0]], [[11.0]]],
                                mask=[[[True]], [[False]], [[False]]])
         np.testing.assert_array_equal(result.data.mask, expected.mask)
@@ -308,16 +351,21 @@ class Test:
         assert self.unstructured_grid_cube.dtype == np.float32
         assert result.dtype == np.float32
 
-    def test_regrid_nearest_unstructured_grid_int(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_regrid_nearest_unstructured_grid_int(self, cache_weights):
         """Test `nearest` regridding with unstructured cube of ints."""
         self.unstructured_grid_cube.data = np.ones((3, 4), dtype=int)
-        result = regrid(self.unstructured_grid_cube,
-                        self.tgt_grid_for_unstructured,
-                        'nearest')
+        result = regrid(
+            self.unstructured_grid_cube,
+            self.tgt_grid_for_unstructured,
+            'nearest',
+            cache_weights=cache_weights,
+        )
         expected = np.array([[[1]], [[1]], [[1]]])
         np.testing.assert_array_equal(result.data, expected)
 
-    def test_invalid_scheme_for_unstructured_grid(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_invalid_scheme_for_unstructured_grid(self, cache_weights):
         """Test invalid scheme for unstructured cube."""
         msg = (
             "Regridding scheme 'invalid' does not support unstructured data, "
@@ -327,9 +375,11 @@ class Test:
                 self.unstructured_grid_cube,
                 self.tgt_grid_for_unstructured,
                 'invalid',
+                cache_weights=cache_weights,
             )
 
-    def test_invalid_scheme_for_irregular_grid(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_invalid_scheme_for_irregular_grid(self, cache_weights):
         """Test invalid scheme for irregular cube."""
         msg = (
             "Regridding scheme 'invalid' does not support irregular data, "
@@ -339,15 +389,18 @@ class Test:
                 self.irregular_grid,
                 self.tgt_grid_for_unstructured,
                 'invalid',
+                cache_weights=cache_weights,
             )
 
-    def test_deprecate_unstrucured_nearest(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_deprecate_unstrucured_nearest(self, cache_weights):
         """Test deprecation of `unstructured_nearest` regridding scheme."""
         with pytest.warns(ESMValCoreDeprecationWarning):
             result = regrid(
                 self.unstructured_grid_cube,
                 self.tgt_grid_for_unstructured,
                 'unstructured_nearest',
+                cache_weights=cache_weights,
             )
         expected = np.ma.array(
             [[[3.0]], [[7.0]], [[11.0]]],
@@ -356,7 +409,8 @@ class Test:
         np.testing.assert_array_equal(result.data.mask, expected.mask)
         np.testing.assert_array_almost_equal(result.data, expected, decimal=6)
 
-    def test_deprecate_linear_extrapolate(self):
+    @pytest.mark.parametrize('cache_weights', [True, False])
+    def test_deprecate_linear_extrapolate(self, cache_weights):
         """Test deprecation of `linear_extrapolate` regridding scheme."""
         data = np.empty((3, 3))
         lons = iris.coords.DimCoord([0, 1.5, 3],
@@ -373,7 +427,12 @@ class Test:
         grid = iris.cube.Cube(data, dim_coords_and_dims=coords_spec)
 
         with pytest.warns(ESMValCoreDeprecationWarning):
-            result = regrid(self.cube, grid, 'linear_extrapolate')
+            result = regrid(
+                self.cube,
+                grid,
+                'linear_extrapolate',
+                cache_weights=cache_weights,
+            )
 
         expected = [[[-3., -1.5, 0.], [0., 1.5, 3.], [3., 4.5, 6.]],
                     [[1., 2.5, 4.], [4., 5.5, 7.], [7., 8.5, 10.]],
diff --git a/tests/unit/preprocessor/_regrid/test_regrid.py b/tests/unit/preprocessor/_regrid/test_regrid.py
index 046da17c8..f7ffff122 100644
--- a/tests/unit/preprocessor/_regrid/test_regrid.py
+++ b/tests/unit/preprocessor/_regrid/test_regrid.py
@@ -1,145 +1,26 @@
 """Unit tests for the :func:`esmvalcore.preprocessor.regrid.regrid`
 function."""
-
-import unittest
-from unittest import mock
-
 import dask
 import dask.array as da
 import iris
 import numpy as np
 import pytest
 
-import tests
+import esmvalcore.preprocessor._regrid
 from esmvalcore.preprocessor import regrid
 from esmvalcore.preprocessor._regrid import (
-    HORIZONTAL_SCHEMES_REGULAR,
+    _get_regridder,
     _horizontal_grid_is_close,
     _rechunk,
 )
 
 
-class Test(tests.Test):
-
-    def _check(self, tgt_grid, scheme, spec=False):
-        expected_scheme = HORIZONTAL_SCHEMES_REGULAR[scheme]
-
-        if spec:
-            spec = tgt_grid
-            self.coord_system.asset_called_once()
-            expected_calls = [
-                mock.call(axis='x', dim_coords=True),
-                mock.call(axis='y', dim_coords=True)
-            ]
-            self.assertEqual(self.tgt_grid_coord.mock_calls, expected_calls)
-            self.regrid.assert_called_once_with(self.tgt_grid, expected_scheme)
-        else:
-            self.regrid.assert_called_once_with(tgt_grid, expected_scheme)
-
-        # Reset the mocks to enable multiple calls per test-case.
-        for mocker in self.mocks:
-            mocker.reset_mock()
-
-    def setUp(self):
-        self.coord_system = mock.Mock(return_value=None)
-        self.coord = mock.sentinel.coord
-        self.coords = mock.Mock(return_value=[self.coord])
-        self.remove_coord = mock.Mock()
-        self.regridded_cube = mock.Mock()
-        self.regridded_cube.data = mock.sentinel.data
-        self.regridded_cube_data = mock.Mock()
-        self.regridded_cube.core_data.return_value = self.regridded_cube_data
-        self.regrid = mock.Mock(return_value=self.regridded_cube)
-        self.src_cube = mock.Mock(
-            spec=iris.cube.Cube,
-            coord_system=self.coord_system,
-            coords=self.coords,
-            remove_coord=self.remove_coord,
-            regrid=self.regrid,
-            dtype=float,
-        )
-        self.src_cube.ndim = 1
-        self.tgt_grid_coord = mock.Mock()
-        self.tgt_grid = mock.Mock(spec=iris.cube.Cube,
-                                  coord=self.tgt_grid_coord)
-        self.regrid_schemes = ['linear', 'nearest', 'area_weighted']
-
-        def _mock_horizontal_grid_is_close(src, tgt):
-            return False
-
-        self.patch('esmvalcore.preprocessor._regrid._horizontal_grid_is_close',
-                   side_effect=_mock_horizontal_grid_is_close)
-
-        def _return_mock_global_stock_cube(
-            spec,
-            lat_offset=True,
-            lon_offset=True,
-        ):
-            return self.tgt_grid
-
-        self.mock_stock = self.patch(
-            'esmvalcore.preprocessor._regrid._global_stock_cube',
-            side_effect=_return_mock_global_stock_cube)
-        self.mocks = [
-            self.coord_system, self.coords, self.regrid, self.src_cube,
-            self.tgt_grid_coord, self.tgt_grid, self.mock_stock
-        ]
-
-    def test_invalid_tgt_grid__unknown(self):
-        dummy = mock.sentinel.dummy
-        scheme = 'linear'
-        emsg = 'Expecting a cube'
-        with self.assertRaisesRegex(ValueError, emsg):
-            regrid(self.src_cube, dummy, scheme)
-
-    def test_invalid_scheme__unknown(self):
-        emsg = "Got invalid regridding scheme string 'wibble'"
-        with self.assertRaisesRegex(ValueError, emsg):
-            regrid(self.src_cube, self.src_cube, 'wibble')
-
-    def test_horizontal_schemes(self):
-        self.assertEqual(set(HORIZONTAL_SCHEMES_REGULAR.keys()),
-                         set(self.regrid_schemes))
-
-    def test_regrid__horizontal_schemes(self):
-        for scheme in self.regrid_schemes:
-            result = regrid(self.src_cube, self.tgt_grid, scheme)
-            self.assertEqual(result, self.regridded_cube)
-            self.assertEqual(result.data, mock.sentinel.data)
-            self._check(self.tgt_grid, scheme)
-
-    def test_regrid__cell_specification(self):
-        specs = ['1x1', '2x2', '3x3', '4x4', '5x5']
-        scheme = 'linear'
-        for spec in specs:
-            result = regrid(self.src_cube, spec, scheme)
-            self.assertEqual(result, self.regridded_cube)
-            self.assertEqual(result.data, mock.sentinel.data)
-            self._check(spec, scheme, spec=True)
-
-    def test_regrid_generic_missing_reference(self):
-        emsg = "No reference specified for generic regridding."
-        with self.assertRaisesRegex(ValueError, emsg):
-            regrid(self.src_cube, '1x1', {})
-
-    def test_regrid_generic_invalid_reference(self):
-        emsg = "Could not import specified generic regridding module."
-        with self.assertRaisesRegex(ValueError, emsg):
-            regrid(self.src_cube, '1x1', {"reference": "this.does:not.exist"})
-
-    def test_regrid_generic_regridding(self):
-        regrid(self.src_cube, '1x1', {"reference": "iris.analysis:Linear"})
-
-    third_party_regridder = mock.Mock()
-
-    def test_regrid_generic_third_party(self):
-        regrid(
-            self.src_cube, '1x1', {
-                "reference": "tests.unit.preprocessor._regrid.test_regrid:"
-                "Test.third_party_regridder",
-                "method": "good",
-            })
-        self.third_party_regridder.assert_called_once_with(method="good")
+@pytest.fixture(autouse=True)
+def clear_regridder_cache(monkeypatch):
+    """Clear regridder cache before test runs."""
+    monkeypatch.setattr(
+        esmvalcore.preprocessor._regrid, '_CACHED_REGRIDDERS', {}
+    )
 
 
 def _make_coord(start: float, stop: float, step: int, *, name: str):
@@ -159,8 +40,9 @@ def _make_cube(*, lat: tuple, lon: tuple):
     lon_coord = _make_coord(*lon, name='longitude')
 
     return iris.cube.Cube(
-        np.empty([len(lat_coord.points),
-                  len(lon_coord.points)]),
+        np.zeros(
+            [len(lat_coord.points), len(lon_coord.points)], dtype=np.float32
+        ),
         dim_coords_and_dims=[(lat_coord, 0), (lon_coord, 1)],
     )
 
@@ -177,6 +59,116 @@ LON_SPEC2 = (5, 355, 35)
 LAT_SPEC3 = (-90, 90, 18)
 LON_SPEC3 = (0, 360, 36)
 
+# 30x30
+LAT_SPEC4 = (-75, 75, 30)
+LON_SPEC4 = (15, 345, 30)
+
+
+@pytest.fixture
+def cube_10x10():
+    """Test cube."""
+    return _make_cube(lat=LAT_SPEC1, lon=LON_SPEC1)
+
+
+@pytest.fixture
+def cube_30x30():
+    """Test cube."""
+    return _make_cube(lat=LAT_SPEC4, lon=LON_SPEC4)
+
+
+SCHEMES = ['area_weighted', 'linear', 'nearest']
+
+
+@pytest.mark.parametrize('cache_weights', [True, False])
+@pytest.mark.parametrize('scheme', SCHEMES)
+def test_builtin_regridding(scheme, cache_weights, cube_10x10, cube_30x30):
+    """Test `regrid.`"""
+    _cached_regridders = esmvalcore.preprocessor._regrid._CACHED_REGRIDDERS
+    assert _cached_regridders == {}
+
+    res = regrid(cube_10x10, cube_30x30, scheme, cache_weights=cache_weights)
+
+    assert res.coord('latitude') == cube_30x30.coord('latitude')
+    assert res.coord('longitude') == cube_30x30.coord('longitude')
+    assert res.dtype == np.float32
+    assert np.allclose(res.data, 0.0)
+
+    if cache_weights:
+        assert len(_cached_regridders) == 1
+        key = (scheme, (18,), (36,), (30,), (30,))
+        assert key in _cached_regridders
+    else:
+        assert not _cached_regridders
+
+
+@pytest.mark.parametrize('scheme', SCHEMES)
+def test_invalid_target_grid(scheme, cube_10x10, mocker):
+    """Test `regrid.`"""
+    target_grid = mocker.sentinel.target_grid
+    msg = "Expecting a cube"
+    with pytest.raises(ValueError, match=msg):
+        regrid(cube_10x10, target_grid, scheme)
+
+
+def test_invalid_scheme(cube_10x10, cube_30x30):
+    """Test `regrid.`"""
+    msg = "Got invalid regridding scheme string 'wibble'"
+    with pytest.raises(ValueError, match=msg):
+        regrid(cube_10x10, cube_30x30, 'wibble')
+
+
+def test_regrid_generic_missing_reference(cube_10x10, cube_30x30):
+    """Test `regrid.`"""
+    msg = "No reference specified for generic regridding."
+    with pytest.raises(ValueError, match=msg):
+        regrid(cube_10x10, cube_30x30, {})
+
+
+def test_regrid_generic_invalid_reference(cube_10x10, cube_30x30):
+    """Test `regrid.`"""
+    msg = "Could not import specified generic regridding module."
+    with pytest.raises(ValueError, match=msg):
+        regrid(cube_10x10, cube_30x30, {'reference': 'this.does:not.exist'})
+
+
+@pytest.mark.parametrize('cache_weights', [True, False])
+def test_regrid_generic_regridding(cache_weights, cube_10x10, cube_30x30):
+    """Test `regrid.`"""
+    _cached_regridders = esmvalcore.preprocessor._regrid._CACHED_REGRIDDERS
+    assert _cached_regridders == {}
+
+    cube_gen = regrid(
+        cube_10x10,
+        cube_30x30,
+        {
+            'reference': 'iris.analysis:Linear',
+            'extrapolation_mode': 'mask',
+        },
+        cache_weights=cache_weights,
+    )
+    cube_lin = regrid(
+        cube_10x10, cube_30x30, 'linear', cache_weights=cache_weights
+    )
+    assert cube_gen.dtype == np.float32
+    assert cube_lin.dtype == np.float32
+    assert cube_gen == cube_lin
+
+    if cache_weights:
+        assert len(_cached_regridders) == 2
+        key_1 = (
+            "{'reference': 'iris.analysis:Linear', 'extrapolation_mode': "
+            "'mask'}",
+            (18,),
+            (36,),
+            (30,),
+            (30,),
+        )
+        key_2 = ('linear', (18,), (36,), (30,), (30,))
+        assert key_1 in _cached_regridders
+        assert key_2 in _cached_regridders
+    else:
+        assert not _cached_regridders
+
 
 @pytest.mark.parametrize(
     'cube2_spec, expected',
@@ -361,5 +353,40 @@ def test_no_rechunk_unsupported_grid():
     assert result.core_data().chunks == expected_chunks
 
 
-if __name__ == '__main__':
-    unittest.main()
+@pytest.mark.parametrize('scheme', SCHEMES)
+def test_regridding_weights_use_cache(scheme, cube_10x10, cube_30x30, mocker):
+    """Test `regrid.`"""
+    _cached_regridders = esmvalcore.preprocessor._regrid._CACHED_REGRIDDERS
+    assert _cached_regridders == {}
+
+    src_lat = cube_10x10.coord('latitude')
+    src_lon = cube_10x10.coord('longitude')
+    tgt_lat = cube_30x30.coord('latitude')
+    tgt_lon = cube_30x30.coord('longitude')
+    key = (scheme, (18,), (36,), (30,), (30,))
+    _cached_regridders[key] = {}
+    _cached_regridders[key][(src_lat, src_lon, tgt_lat, tgt_lon)] = (
+        mocker.sentinel.regridder
+    )
+    mock_load_scheme = mocker.patch.object(
+        esmvalcore.preprocessor._regrid, '_load_scheme', autospec=True
+    )
+
+    reg = _get_regridder(cube_10x10, cube_30x30, scheme, cache_weights=True)
+
+    assert reg == mocker.sentinel.regridder
+
+    assert len(_cached_regridders) == 1
+    assert key in _cached_regridders
+
+    mock_load_scheme.assert_not_called()
+
+
+def test_clear_regridding_weights_cache():
+    """Test `regrid.cache_clear().`"""
+    _cached_regridders = esmvalcore.preprocessor._regrid._CACHED_REGRIDDERS
+    _cached_regridders['test'] = 'test'
+
+    regrid.cache_clear()
+
+    assert _cached_regridders == {}