diff --git a/docs/source/reference/operators.rst b/docs/source/reference/operators.rst index d57ad7d6e..5274c7917 100644 --- a/docs/source/reference/operators.rst +++ b/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: diff --git a/src/CSET/operators/__init__.py b/src/CSET/operators/__init__.py index 5a07151e9..0b70e4f32 100644 --- a/src/CSET/operators/__init__.py +++ b/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", diff --git a/src/CSET/operators/mesoscale.py b/src/CSET/operators/mesoscale.py new file mode 100644 index 000000000..900a12f36 --- /dev/null +++ b/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 diff --git a/src/CSET/recipes/Example_Gaussian_Spatial_Perturbation.yaml b/src/CSET/recipes/Example_Gaussian_Spatial_Perturbation.yaml new file mode 100644 index 000000000..78cbbbb25 --- /dev/null +++ b/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 diff --git a/src/CSET/recipes/Example_Uniform_Spatial_Perturbation.yaml b/src/CSET/recipes/Example_Uniform_Spatial_Perturbation.yaml new file mode 100644 index 000000000..c78033302 --- /dev/null +++ b/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 diff --git a/tests/operators/test_mesoscale.py b/tests/operators/test_mesoscale.py new file mode 100644 index 000000000..cbb1311d5 --- /dev/null +++ b/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, + )