Coverage for src/CSET/operators/mesoscale.py: 100%

1# © Crown copyright, Met Office (2022-2024) and CSET contributors. 

2# 

3# Licensed under the Apache License, Version 2.0 (the "License"); 

4# you may not use this file except in compliance with the License. 

5# You may obtain a copy of the License at 

6# 

7# http://www.apache.org/licenses/LICENSE-2.0 

8# 

9# Unless required by applicable law or agreed to in writing, software 

10# distributed under the License is distributed on an "AS IS" BASIS, 

11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 

12# See the License for the specific language governing permissions and 

13# limitations under the License. 

14 

15"""A module containing different diagnostics for mesoscales. 

16 

17The diagnostics here are applicable at mesoscales and apply generally 

18rather than for specific aspects of mesoscale meteorology (e.g. convection). 

19For specific aspects, the user is referred to other modules available in 

20CSET. 

21 

22""" 

23 

24import logging 

25 

26import iris 

27from scipy.ndimage import gaussian_filter, uniform_filter 

28 

29from CSET.operators._utils import get_cube_yxcoordname 

30 

31logger = logging.getLogger(__name__) 

32 

33 

34def spatial_perturbation_field( 

35 original_field: iris.cube.Cube, 

36 apply_gaussian_filter: bool = True, 

37 filter_scale: int = 40, 

38) -> iris.cube.Cube: 

39 """Calculate a spatial perturbation field. 

40 

41 Parameters 

42 ---------- 

43 original_field: iris.cube.Cube 

44 Iris cube containing data to smooth, supporting multiple dimensions 

45 (at least two spatial dimensions must be supplied, i.e. 2D). 

46 apply_gaussian_filter: boolean, optional 

47 If set to True a Gaussian filter is applied; if set to False 

48 a Uniform filter is applied. 

49 Default is True. 

50 filter_scale: int, optional 

51 Scale at which to define the filter in grid boxes. If the 

52 filter is a Gaussian convolution this value represents the 

53 standard deviation of the Gaussian kernel. 

54 Default is 40 grid boxes. 

55 

56 Returns 

57 ------- 

58 pert_field: iris.cube.Cube 

59 An iris cube of the spatial perturbation field. 

60 

61 Notes 

62 ----- 

63 In mesoscale meteorology the perturbation field is more important than the 

64 balanced flows for process understanding. This function is designed 

65 to create spatial perturbation fields based on smoothing with a Gaussian 

66 kernel or a uniform kernel. 

67 

68 The kernels are defined by the filter_scale, which for mesoscale 

69 perturbations should be between an approximate cloud separation distance, 

70 of 30 km, and synoptic scale variations (1000 km). In practice any 

71 value between these ranges should provided broadly consistent results (e.g. 

72 [Flacketal2016]_). The Gaussian kernel will give greater importance to 

73 areas closer to the event and will produce a smooth perturbation field. 

74 The uniform kernel will produce a smooth perturbation field but will not 

75 give local features as much prominence. 

76 

77 Caution should be applied to boundaries, particularly if the domain is of 

78 variable resolution, as some numerical artifacts could be introduced. 

79 

80 Examples 

81 -------- 

82 >>> Temperature_perturbation = meso.spatial_perturbation_fields(Temp, 

83 gaussian_filter=True,filter_scale=40) 

84 >>> iplt.pcolormesh(Temperature_perturabtion[0,:,:],cmap=mpl.cm.bwr) 

85 >>> plt.gca().coastlines('10m') 

86 >>> plt.clim(-5,5) 

87 >>> plt.colorbar() 

88 >>> plt.show() 

89 

90 """ 

91 pert_field = original_field.copy() 

92 # find axes of spatial coordinates in field 

93 coords = [coord.name() for coord in original_field.coords()] 

94 # axes tuple containing latitude, longitude coordinate name. 

95 axes = ( 

96 coords.index(get_cube_yxcoordname(original_field)[0]), 

97 coords.index(get_cube_yxcoordname(original_field)[1]), 

98 ) 

99 # apply convolution depending on type used 

100 if apply_gaussian_filter: 

101 filter_type = "Gaussian" 

102 logger.info("Gaussian filter applied.") 

103 pert_field.data -= gaussian_filter(original_field.data, filter_scale, axes=axes) 

104 else: 

105 logger.info("Uniform filter applied.") 

106 filter_type = "Uniform" 

107 pert_field.data -= uniform_filter(original_field.data, filter_scale, axes=axes) 

108 # provide attributes to cube to indicate spatial perturbation field 

109 pert_field.attributes["perturbation_field"] = ( 

110 f"{filter_type}_with_{filter_scale}_grid_point_filter_scale" 

111 ) 

112 return pert_field