Merge pull request #979 from MetOffice/950_spatial_perturbation_fields

Adds spatial perturbation operator

docs/source/reference/operators.rst

@@ -82,3 +82,12 @@ CSET.operators.ageofair
 
 .. automodule:: CSET.operators.ageofair
    :members:
+
+Mesoscale Operators
+~~~~~~~~~~~~~~~~~~~
+
+CSET.operators.mesoscale
+------------------------
+
+.. automodule:: CSET.operators.mesoscale
+   :members:

src/CSET/operators/__init__.py

@@ -33,6 +33,7 @@
     constraints,
     convection,
     filters,
+    mesoscale,
     misc,
     plot,
     read,
@@ -51,6 +52,7 @@
     "execute_recipe",
     "filters",
     "get_operator",
+    "mesoscale",
     "misc",
     "plot",
     "read",

src/CSET/operators/mesoscale.py

@@ -0,0 +1,117 @@
+# © Crown copyright, Met Office (2022-2024) and CSET contributors.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""A module containing different diagnostics for mesoscales.
+
+The diagnostics here are applicable at mesoscales and apply generally
+rather than for specific aspects of mesoscale meteorology (e.g. convection).
+For specific aspects, the user is referred to other modules available in
+CSET.
+
+"""
+
+import logging
+
+import iris
+from scipy.ndimage import gaussian_filter, uniform_filter
+
+from CSET.operators._utils import get_cube_yxcoordname
+
+
+def spatial_perturbation_field(
+    original_field: iris.cube.Cube,
+    apply_gaussian_filter: bool = True,
+    filter_scale: int = 40,
+) -> iris.cube.Cube:
+    """Calculate a spatial perturbation field.
+
+    Parameters
+    ----------
+    original_field: iris.cube.Cube
+        Iris cube containing data to smooth, supporting multiple dimensions
+        (at least two spatial dimensions must be supplied, i.e. 2D).
+    apply_gaussian_filter: boolean, optional
+        If set to True a Gaussian filter is applied; if set to False
+        a Uniform filter is applied.
+        Default is True.
+    filter_scale: int, optional
+        Scale at which to define the filter in grid boxes. If the
+        filter is a Gaussian convolution this value represents the
+        standard deviation of the Gaussian kernel.
+        Default is 40 grid boxes.
+
+    Returns
+    -------
+    pert_field: iris.cube.Cube
+        An iris cube of the spatial perturbation field.
+
+    Notes
+    -----
+    In mesoscale meteorology the perturbation field is more important than the
+    balanced flows for process understanding. This function is designed
+    to create spatial perturbation fields based on smoothing with a Gaussian
+    kernel or a uniform kernel.
+
+    The kernels are defined by the filter_scale, which for mesoscale
+    perturbations should be between an approximate cloud separation distance,
+    of 30 km, and synoptic scale variations (1000 km). In practice any
+    value between these ranges should provided broadly consistent results (e.g.
+    [Flacketal2016]_). The Gaussian kernel will give greater importance to
+    areas closer to the event and will produce a smooth perturbation field.
+    The uniform kernel will produce a smooth perturbation field but will not
+    give local features as much prominence.
+
+    Caution should be applied to boundaries, particularly if the domain is of
+    variable resolution, as some numerical artifacts could be introduced.
+
+    References
+    ----------
+    .. [Flacketal2016] Flack, D.L.A., Plant, R.S., Gray, S.L., Lean, H.W.,
+       Keil, C. and Craig, G.C. (2016) "Characterisation of Convective
+       Regimes over the British Isles." Quarterly Journal of the Royal
+       Meteorological Society, vol. 142, 1541-1553. doi:10.1002/qj.2758
+
+    Examples
+    --------
+    >>> Temperature_perturbation = meso.spatial_perturbation_fields(Temp,
+                                   gaussian_filter=True,filter_scale=40)
+    >>> iplt.pcolormesh(Temperature_perturabtion[0,:,:],cmap=mpl.cm.bwr)
+    >>> plt.gca().coastlines('10m')
+    >>> plt.clim(-5,5)
+    >>> plt.colorbar()
+    >>> plt.show()
+
+    """
+    pert_field = original_field.copy()
+    # find axes of spatial coordinates in field
+    coords = [coord.name() for coord in original_field.coords()]
+    # axes tuple containing latitude, longitude coordinate name.
+    axes = (
+        coords.index(get_cube_yxcoordname(original_field)[0]),
+        coords.index(get_cube_yxcoordname(original_field)[1]),
+    )
+    # apply convolution depending on type used
+    if apply_gaussian_filter:
+        filter_type = "Gaussian"
+        logging.info("Gaussian filter applied.")
+        pert_field.data -= gaussian_filter(original_field.data, filter_scale, axes=axes)
+    else:
+        logging.info("Uniform filter applied.")
+        filter_type = "Uniform"
+        pert_field.data -= uniform_filter(original_field.data, filter_scale, axes=axes)
+    # provide attributes to cube to indicate spatial perturbation field
+    pert_field.attributes["perturbation_field"] = (
+        f"{filter_type}_with_{str(filter_scale)}_grid_point_filter_scale"
+    )
+    return pert_field

src/CSET/recipes/Example_Gaussian_Spatial_Perturbation.yaml

@@ -0,0 +1,28 @@
+category: Surface Spatial Plot
+title: Mesoscale Gaussian Filter Spatial Perturbation Plot
+description: |
+  Creates a spatial map of a perturbation field based on a Gaussian filter
+  over a specified halfwidth in number of grid points.
+
+steps:
+  - operator: read.read_cubes
+
+  - operator: mesoscale.spatial_perturbation_field
+    original_field:
+        operator: filters.filter_cubes
+        constraint:
+            operator: constraints.combine_constraints
+            variable_constraint:
+                operator: constraints.generate_var_constraint
+                varname: surface_air_temperature
+            cell_method_constraint:
+                operator: constraints.generate_cell_methods_constraint
+                cell_methods: []
+    Gaussian_filter: True
+    filter_scale: 40
+
+  - operator: plot.spatial_pcolormesh_plot
+    sequence_coordinate: time
+
+  - operator: write.write_cube_to_nc
+    overwrite: True

src/CSET/recipes/Example_Uniform_Spatial_Perturbation.yaml

@@ -0,0 +1,28 @@
+category: Surface Spatial Plot
+title: Mesoscale Uniform Filter Spatial Perturbation Plot
+description: |
+  Creates a spatial map of a perturbation field based on a uniform filter
+  over a specified filter scale in number of grid points.
+
+steps:
+  - operator: read.read_cubes
+
+  - operator: mesoscale.spatial_perturbation_field
+    original_field:
+        operator: filters.filter_cubes
+        constraint:
+            operator: constraints.combine_constraints
+            variable_constraint:
+                operator: constraints.generate_var_constraint
+                varname: surface_air_temperature
+            cell_method_constraint:
+                operator: constraints.generate_cell_methods_constraint
+                cell_methods: []
+    Gaussian_filter: False
+    filter_scale: 40
+
+  - operator: plot.spatial_pcolormesh_plot
+    sequence_coordinate: time
+
+  - operator: write.write_cube_to_nc
+    overwrite: True

tests/operators/test_mesoscale.py

@@ -0,0 +1,75 @@
+# © Crown copyright, Met Office (2022-2024) and CSET contributors.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Test mesoscale operators."""
+
+import numpy as np
+from scipy.ndimage import gaussian_filter, uniform_filter
+
+import CSET.operators.mesoscale as mesoscale
+from CSET.operators._utils import get_cube_yxcoordname
+
+
+def test_spatial_perturbation_field_gaussian(cube):
+    """Test smoothing a cube with a Gaussian filter."""
+    calculated = cube.copy()
+    coords = [coord.name() for coord in cube.coords()]
+    axes = (
+        coords.index(get_cube_yxcoordname(cube)[0]),
+        coords.index(get_cube_yxcoordname(cube)[1]),
+    )
+    calculated.data -= gaussian_filter(cube.data, 40, axes=axes)
+    assert np.allclose(
+        calculated.data,
+        mesoscale.spatial_perturbation_field(cube).data,
+        rtol=1e-06,
+        atol=1e-02,
+    )
+
+    calculated_2 = cube.copy()
+    calculated_2.data -= gaussian_filter(cube.data, 100, axes=axes)
+    assert np.allclose(
+        calculated_2.data,
+        mesoscale.spatial_perturbation_field(cube, filter_scale=100).data,
+        rtol=1e-06,
+        atol=1e-02,
+    )
+
+
+def test_spatial_perturbation_field_uniform(cube):
+    """Test smoothing a cube with a uniform filter."""
+    calculated = cube.copy()
+    coords = [coord.name() for coord in cube.coords()]
+    axes = (
+        coords.index(get_cube_yxcoordname(cube)[0]),
+        coords.index(get_cube_yxcoordname(cube)[1]),
+    )
+    calculated.data -= uniform_filter(cube.data, 40, axes=axes)
+    assert np.allclose(
+        calculated.data,
+        mesoscale.spatial_perturbation_field(cube, apply_gaussian_filter=False).data,
+        rtol=1e-06,
+        atol=1e-02,
+    )
+
+    calculated_2 = cube.copy()
+    calculated_2.data -= uniform_filter(cube.data, 100, axes=axes)
+    assert np.allclose(
+        calculated_2.data,
+        mesoscale.spatial_perturbation_field(
+            cube, apply_gaussian_filter=False, filter_scale=100
+        ).data,
+        rtol=1e-06,
+        atol=1e-02,
+    )