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

1# © Crown copyright, Met Office (2022-2025) 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""" 

16A module containing different diagnostics for ensembles. 

17 

18The diagnostics here are applicable to ensembles and apply generally to 

19ensembles. They are not just limited to considering ensemble spread. 

20""" 

21 

22import iris 

23import iris.cube 

24import numpy as np 

25 

26 

27def DKE(u: iris.cube.Cube, v: iris.cube.Cube) -> iris.cube.Cube: 

28 r"""Calculate the Difference Kinetic Energy (DKE). 

29 

30 Parameters 

31 ---------- 

32 u: iris.cube.Cube 

33 Iris cube of the u component of the wind field for all members. Must 

34 include a realization coordinate; realization must be the first 

35 coordinate. 

36 v: iris.cube.Cube 

37 Iris cube of the v component of the wind field same format as u. 

38 

39 Returns 

40 ------- 

41 DKE: iris.cube.Cube 

42 An iris cube of the DKE for each of the control - member comparisons. 

43 

44 Notes 

45 ----- 

46 The Difference Kinetic Energy, or DKE was first used in an ensemble sense 

47 by [Zhangetal2002]_. It is calculated as 

48 

49 .. math:: \frac{1}{2}u'u' + \frac{1}{2}v'v' 

50 

51 where 

52 

53 .. math:: x' = x_{control} - x_{perturbed} 

54 

55 for each member of the ensemble. 

56 

57 The DKE is used to show links between the perturbation growth (growth of 

58 ensemble spread) or errors with dynamical features. It can be particularly 

59 useful in understanding differences in physical mechanisms between ensemble 

60 members. The larger the DKE the greater the difference between the two 

61 members being considered. 

62 

63 The DKE can be viewed as a domain average, horizontal integration (to 

64 produce profiles), or vertical integration/subsampling (to provide spatial maps). 

65 Furthermore, weightings based on volume or mass can be applied to these 

66 integrations. However, initially the DKE per unit mass is implemented here. 

67 

68 The DKE is often considered in the form of power spectra to identify the 

69 presence of upscale error growth or the dominant scales of error growth. It 

70 is often plotted on a logarithmic scale. 

71 

72 Examples 

73 -------- 

74 >>> DKE = ensembles.DKE(u, v) 

75 """ 

76 for cube in [u, v]: 

77 # Use dim_map to store the dimensions and check the realization 

78 # coordinate is the first coordinate. 

79 dim_map = {} 

80 for coord in cube.dim_coords: 

81 dim_index = cube.coord_dims(coord.name())[0] 

82 dim_map[coord.name()] = dim_index 

83 if "realization" not in dim_map: 

84 raise ValueError("Cube should have a realization coordinate.") 

85 if not dim_map["realization"] == 0: 

86 raise ValueError("Realization should be the first coordinate in cube.") 

87 

88 # Check cubes are the same shape. 

89 if u.shape != v.shape: 

90 raise ValueError("Cubes should be the same shape.") 

91 

92 # Check coordinates are identical. 

93 for coord_u, coord_v in zip(u.dim_coords, v.dim_coords, strict=True): 

94 if not u.coord(coord_u.name()) == v.coord(coord_v.name()): 

95 raise ValueError( 

96 f"u and v should have matching coordinates for {coord_u.name()}." 

97 ) 

98 # Define control member and perturbed members. 

99 u_ctrl = u[np.where(u.coord("realization").points[:] == 0)[0][0], :] 

100 u_mem = u[np.where(u.coord("realization").points[:] != 0)[0][:], :] 

101 v_ctrl = v[np.where(v.coord("realization").points[:] == 0)[0][0], :] 

102 v_mem = v[np.where(v.coord("realization").points[:] != 0)[0][:], :] 

103 

104 # Calculate the DKE 

105 DKE = u_mem.copy() 

106 DKE.data[:] = 0.0 

107 DKE.data = ( 

108 0.5 * (u_ctrl.data - u_mem.data) ** 2 + 0.5 * (v_ctrl.data - v_mem.data) ** 2 

109 ) 

110 DKE.rename("Difference Kinetic Energy") 

111 return DKE