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 

31 

32def spatial_perturbation_field( 

33 original_field: iris.cube.Cube, 

34 apply_gaussian_filter: bool = True, 

35 filter_scale: int = 40, 

36) -> iris.cube.Cube: 

37 """Calculate a spatial perturbation field. 

38 

39 Parameters 

40 ---------- 

41 original_field: iris.cube.Cube 

42 Iris cube containing data to smooth, supporting multiple dimensions 

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

44 apply_gaussian_filter: boolean, optional 

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

46 a Uniform filter is applied. 

47 Default is True. 

48 filter_scale: int, optional 

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

50 filter is a Gaussian convolution this value represents the 

51 standard deviation of the Gaussian kernel. 

52 Default is 40 grid boxes. 

53 

54 Returns 

55 ------- 

56 pert_field: iris.cube.Cube 

57 An iris cube of the spatial perturbation field. 

58 

59 Notes 

60 ----- 

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

62 balanced flows for process understanding. This function is designed 

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

64 kernel or a uniform kernel. 

65 

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

67 perturbations should be between an approximate cloud separation distance, 

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

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

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

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

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

73 give local features as much prominence. 

74 

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

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

77 

78 References 

79 ---------- 

80 .. [Flacketal2016] Flack, D.L.A., Plant, R.S., Gray, S.L., Lean, H.W., 

81 Keil, C. and Craig, G.C. (2016) "Characterisation of Convective 

82 Regimes over the British Isles." Quarterly Journal of the Royal 

83 Meteorological Society, vol. 142, 1541-1553. doi:10.1002/qj.2758 

84 

85 Examples 

86 -------- 

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

88 gaussian_filter=True,filter_scale=40) 

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

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

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

92 >>> plt.colorbar() 

93 >>> plt.show() 

94 

95 """ 

96 pert_field = original_field.copy() 

97 # find axes of spatial coordinates in field 

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

99 # axes tuple containing latitude, longitude coordinate name. 

100 axes = ( 

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

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

103 ) 

104 # apply convolution depending on type used 

105 if apply_gaussian_filter: 

106 filter_type = "Gaussian" 

107 logging.info("Gaussian filter applied.") 

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

109 else: 

110 logging.info("Uniform filter applied.") 

111 filter_type = "Uniform" 

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

113 # provide attributes to cube to indicate spatial perturbation field 

114 pert_field.attributes["perturbation_field"] = ( 

115 f"{filter_type}_with_{str(filter_scale)}_grid_point_filter_scale" 

116 ) 

117 return pert_field