Coverage for src/CSET/operators/curvature.py: 11%

1# © Crown copyright, Met Office (2022-2026) 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 implementations of CURV. 

16 

17The diagnostics here are separate implementations of CURV for both 

18global models and limited-area models. 

19""" 

20 

21import logging 

22 

23import iris 

24import numpy as np 

25 

26from CSET._common import iter_maybe 

27 

28 

29def _lat_lon_identifier(original_lat, original_lon, distance, bearing): 

30 """Find the latitude and longitude of new points based on great circle distance.""" 

31 R = 6371.0 # km. 

32 lat_rad = np.radians(original_lat) 

33 lon_rad = np.radians(original_lon) 

34 bearing_rad = np.radians(bearing) 

35 

36 new_lat = np.arcsin( 

37 np.sin(lat_rad) * np.cos(distance / R) 

38 + np.cos(lat_rad) * np.sin(distance / R) * np.cos(bearing_rad) 

39 ) 

40 new_lon = lon_rad + np.arctan2( 

41 np.sin(bearing_rad) * np.sin(distance / R) * np.cos(lat_rad), 

42 np.cos(distance / R) - np.sin(lat_rad) * np.sin(new_lat), 

43 ) 

44 

45 new_lat = np.degrees(new_lat) 

46 new_lon = np.degrees(new_lon) 

47 new_lon = np.where(new_lon > 180.0, new_lon - 360.0, new_lon) 

48 

49 return new_lat, new_lon 

50 

51 

52def global_curv(central, radius, num_radial_points=16, tol=400.0): 

53 """Calculate the CURV diagnostic for a global model VERSION 2 for faster calculation speed.""" 

54 if (num_radial_points == 8) or (num_radial_points == 16): 

55 bearing = np.linspace(0.0, 360.0, num_radial_points + 1)[:-1] 

56 else: 

57 raise ValueError( 

58 f"Number of radial points should be 8 or 16. You have provided {num_radial_points}." 

59 ) 

60 

61 logging.info( 

62 f"CURV calculated base on {num_radial_points} radial points, a {radius} km radius and {tol} tolerance." 

63 ) 

64 curv_cubes = iris.cube.CubeList([]) 

65 for cube in iter_maybe(central): 

66 curv = cube.copy() 

67 surroundings = [] 

68 lats, lons = np.meshgrid( 

69 cube.coord("latitude").points, cube.coord("longitude").points 

70 ) 

71 for b in bearing: 

72 target_data = cube.copy() 

73 new_lat, new_lon = _lat_lon_identifier( 

74 lats, 

75 lons, 

76 radius, 

77 b, 

78 ) 

79 new_lat_points = new_lat[0, :] 

80 new_lon_points = new_lon[:, 0] 

81 surround = target_data.interpolate( 

82 [("latitude", new_lat_points), ("longitude", new_lon_points)], 

83 iris.analysis.Linear(), 

84 ) 

85 surroundings.append(surround.data) 

86 surroundings = np.array(surroundings) 

87 curv_data = surroundings - cube.data 

88 curv_data_binary = np.zeros(curv_data.shape) 

89 curv_data_binary[curv_data > tol] = -1.0 

90 curv_data_binary[curv_data < -tol] = 1.0 

91 curv_data_collapsed = np.sum(curv_data_binary, 0) 

92 curv.data = curv_data_collapsed 

93 curv.rename(f"CURV_calculated_from_{num_radial_points}_radial_points") 

94 curv.units = "1" 

95 curv_cubes.append(curv) 

96 if len(curv_cubes) == 1: 

97 return curv_cubes[0] 

98 else: 

99 return curv_cubes