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 

26 

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

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

29 R = 6371.0 # km. 

30 lat_rad = np.radians(original_lat) 

31 lon_rad = np.radians(original_lon) 

32 bearing_rad = np.radians(bearing) 

33 

34 new_lat = np.arcsin( 

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

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

37 ) 

38 new_lon = lon_rad + np.arctan2( 

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

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

41 ) 

42 

43 new_lat = np.degrees(new_lat) 

44 new_lon = np.degrees(new_lon) 

45 new_lon = np.where(new_lon > 360.0, new_lon - 360.0, new_lon) 

46 

47 return new_lat, new_lon 

48 

49 

50def global_curv(central, radius, num_radial_points=16, tol=0): 

51 """Calculate the CURV diagnostic for a global model.""" 

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

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

54 else: 

55 raise ValueError( 

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

57 ) 

58 

59 logging.info( 

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

61 ) 

62 

63 surroundings = iris.cube.CubeList([]) 

64 for b in bearing: 

65 surround = central.copy() 

66 # Need to try and avoid this loop for whole domain if possible, also need to raise time and realization bits as well. 

67 for lat_number, lat in enumerate(central.slices_over("latitude")): 

68 for lon_number, lon in enumerate(lat.slices_over("longitude")): 

69 new_lat, new_lon = _lat_lon_identifier( 

70 lon.coord("latitude").points, 

71 lon.coord("longitude").points, 

72 radius, 

73 b, 

74 ) 

75 # Need to get interpolation better. 

76 surround.data[lat_number, lon_number] = central.data[ 

77 np.abs(central.coord("latitude").points - new_lat[0]).argmin(), 

78 np.abs(central.coord("longitude").points - new_lon[0]).argmin(), 

79 ] 

80 surround.add_aux_coord( 

81 iris.coords.DimCoord(b, long_name="bearing", units="degrees") 

82 ) 

83 surroundings.append(surround) 

84 surroundings.merge() 

85 print(surroundings) 

86 

87 curv = surroundings[0].copy() 

88 curv -= central 

89 

90 curv.data[curv.data > tol] = -1.0 

91 curv.data[curv.data < tol] = 1.0 

92 

93 curv.collapsed("bearing", iris.analysis.SUM) 

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

95 curv.units = "1" 

96 return curv