Coverage for src/CSET/operators/transect.py: 95%

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"""Operators to extract a transect given a tuple of xy coords to start/finish.""" 

16 

17import logging 

18 

19import iris 

20import numpy as np 

21 

22from CSET.operators._utils import get_cube_yxcoordname 

23 

24 

25def _check_within_bounds(point: tuple[float, float], lat_coord, lon_coord): 

26 """Check if the point (lat, lon) is within the bounds of the data.""" 

27 lon_min = min(lon_coord.points) 

28 lon_max = max(lon_coord.points) 

29 lat_min = min(lat_coord.points) 

30 lat_max = max(lat_coord.points) 

31 if ( 

32 lat_min > point[0] 

33 or lat_max < point[0] 

34 or lon_min > point[1] 

35 or lon_max < point[1] 

36 ): 

37 raise IndexError( 

38 f"Point {point} not between {(lat_min, lon_min)} and {(lat_max, lon_max)}" 

39 ) 

40 

41 

42def calc_transect(cube: iris.cube.Cube, startcoords: tuple, endcoords: tuple): 

43 """Compute transect between startcoords and endcoords. 

44 

45 Computes a transect for a given cube containing at least latitude 

46 and longitude coordinates, using an appropriate sampling interval along the 

47 transect based on grid spacing. Also decides a suitable x coordinate to plot along 

48 the transect - see notes for details on this. 

49 

50 Arguments 

51 --------- 

52 

53 cube: Cube 

54 An iris cube containing latitude and longitude coordinate dimensions, 

55 to compute the transect on. 

56 startcoords: tuple 

57 A tuple containing the start coordinates for the transect using model coordinates, 

58 ordered (latitude,longitude). 

59 endcoords: tuple 

60 A tuple containing the end coordinates for the transect using model coordinates, 

61 ordered (latitude,longitude). 

62 

63 Returns 

64 ------- 

65 cube: Cube 

66 A cube containing at least pressure and the coordinate specified by coord, for 

67 the transect specified between startcoords and endcoords. 

68 

69 Notes 

70 ----- 

71 This operator uses the iris.Nearest method to interpolate the specific point along 

72 the transect. 

73 Identification of an appropriate coordinate to plot along the x axis is done by 

74 determining the largest distance out of latitude and longitude. For example, if 

75 the transect is 90 degrees (west to east), then delta latitude is zero, so it will 

76 always plot as a function of longitude. Vice versa if the transect is 180 degrees 

77 (south to north). If a transect is 45 degrees, and delta latitude and longitude 

78 are the same, it will default to plotting as a function of longitude. Note this 

79 doesn't affect the transect plot, its purely for interpretation with some appropriate 

80 x axis labelling/points. 

81 """ 

82 # Find out xy coord name 

83 lat_name, lon_name = get_cube_yxcoordname(cube) 

84 

85 lon_coord = cube.coord(lon_name) 

86 lat_coord = cube.coord(lat_name) 

87 

88 _check_within_bounds(startcoords, lat_coord, lon_coord) 

89 _check_within_bounds(endcoords, lat_coord, lon_coord) 

90 

91 # Compute vector distance between start and end points in degrees. 

92 dist_deg = np.sqrt( 

93 (startcoords[0] - endcoords[0]) ** 2 + (startcoords[1] - endcoords[1]) ** 2 

94 ) 

95 

96 # Compute minimum gap between x/y spatial coords. 

97 lon_min = np.abs(np.min(lon_coord.points[1:] - lon_coord.points[:-1])) 

98 lat_min = np.abs(np.min(lat_coord.points[1:] - lat_coord.points[:-1])) 

99 

100 # For scenarios where coord is at 90 degree to the grid (i.e. no 

101 # latitude/longitude change). Only xmin or ymin will be zero, not both 

102 # (otherwise startcoords and endcoords the same). 

103 if startcoords[1] == endcoords[1]: 

104 # Along latitude. 

105 lon_pnts = np.repeat(startcoords[1], int(dist_deg / lat_min)) 

106 lat_pnts = np.linspace(startcoords[0], endcoords[0], int(dist_deg / lat_min)) 

107 transect_coord = "latitude" 

108 elif startcoords[0] == endcoords[0]: 

109 # Along longitude. 

110 lon_pnts = np.linspace(startcoords[1], endcoords[1], int(dist_deg / lon_min)) 

111 lat_pnts = np.repeat(startcoords[0], int(dist_deg / lon_min)) 

112 transect_coord = "longitude" 

113 else: 

114 # Else use the smallest grid space in x or y direction. 

115 number_of_points = int(dist_deg / np.min([lon_min, lat_min])) 

116 lon_pnts = np.linspace(startcoords[1], endcoords[1], number_of_points) 

117 lat_pnts = np.linspace(startcoords[0], endcoords[0], number_of_points) 

118 

119 # If change in latitude larger than change in longitude: 

120 if abs(startcoords[0] - endcoords[0]) > abs(startcoords[1] - endcoords[1]): 

121 transect_coord = "latitude" 

122 elif abs(startcoords[0] - endcoords[0]) <= abs(startcoords[1] - endcoords[1]): 122 ↛ 126line 122 didn't jump to line 126 because the condition on line 122 was always true

123 transect_coord = "longitude" 

124 

125 # Create cubelist to store interpolated points along transect. 

126 interpolated_cubes = iris.cube.CubeList() 

127 

128 # Iterate over all points along transect, lon_pnts will be the same shape as 

129 # lat_pnts so we can use either to iterate over. 

130 for i in range(0, lon_pnts.shape[0]): 

131 logging.info("%s/%s", i + 1, lon_pnts.shape[0]) 

132 

133 # Get point along transect. 

134 cube_slice = cube.interpolate( 

135 [(lon_name, lon_pnts[i]), (lat_name, lat_pnts[i])], iris.analysis.Nearest() 

136 ) 

137 

138 # Remove existing coordinates ready to add one single map coordinate 

139 # Note latitude/longitude cubes may have additional AuxCoord to remove 

140 for coord_name in ["latitude", "longitude", "grid_latitude", "grid_longitude"]: 

141 try: 

142 cube_slice.remove_coord(coord_name) 

143 except iris.exceptions.CoordinateNotFoundError: 

144 pass 

145 

146 if transect_coord == "latitude": 

147 dist_coord = iris.coords.DimCoord( 

148 lat_pnts[i], long_name="latitude", units="degrees" 

149 ) 

150 cube_slice.add_aux_coord(dist_coord) 

151 cube_slice = iris.util.new_axis(cube_slice, scalar_coord="latitude") 

152 elif transect_coord == "longitude": 152 ↛ 159line 152 didn't jump to line 159 because the condition on line 152 was always true

153 dist_coord = iris.coords.DimCoord( 

154 lon_pnts[i], long_name="longitude", units="degrees" 

155 ) 

156 cube_slice.add_aux_coord(dist_coord) 

157 cube_slice = iris.util.new_axis(cube_slice, scalar_coord="longitude") 

158 

159 interpolated_cubes.append(cube_slice) 

160 

161 # Concatenate into single cube. 

162 interpolated_cubes = interpolated_cubes.concatenate() 

163 

164 # Add metadata to interpolated cubes showing coordinates. 

165 interpolated_cubes[0].attributes["transect_coords"] = ( 

166 f"{startcoords[0]}_{startcoords[1]}_{endcoords[0]}_{endcoords[1]}" 

167 ) 

168 # If concatenation successful, should be CubeList with one cube left. 

169 assert len(interpolated_cubes) == 1 

170 return interpolated_cubes[0]