Coverage for src/CSET/operators/feature.py: 96%
41 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-06 11:34 +0000
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"""Operators for identifying and tracking features."""
16import logging
17import os
19import iris
20import iris.cube
21import iris.util
22import numpy as np
23from simpletrack.track import Tracker
26def track(
27 cube: iris.cube.Cube,
28 threshold: float,
29 under_threshold: bool = False,
30 min_size: int = 4,
31 retain_lifetime_on_split: bool = True,
32 tracking_nbhood: int = 5,
33 overlap_threshold: float = 0.3,
34 save_data: bool = False,
35):
36 """Track features between subsequent timesteps.
38 Parameters
39 ----------
40 cube: iris.cube.Cube
41 The cube to identify features in. The cube must be 3D and contain a time coordinate
42 and horizontal coordinates of xy type (not latitude/longitude).
43 threshold: float
44 The threshold value for feature detection.
45 under_threshold: bool, optional
46 If set to True, features are identified where the data is below the threshold.
47 If set to False, features are identified where the data is above the threshold.
48 Default is False.
49 min_size: int, optional
50 The minimum number of contiguous grid points required for a feature to be tracked.
51 Default is 4.
52 retain_lifetime_on_split: bool, optional
53 If set to True, the lifetime of a feature is retained when it splits into
54 multiple features. If set to False, the lifetime is reset when a feature splits.
55 Default is True.
56 tracking_nbhood: int, optional
57 The size of the neighbourhood used for tracking features between timesteps.
58 This dictates the maximum pixel radius from a feature centroid at which new features could
59 reasonably be spawned.
60 Default is 5.
61 overlap_threshold: float, optional
62 The minimum overlap required between features in consecutive timesteps for
63 them to be considered the same feature.
64 Default is 0.3.
65 save_data: bool, optional
66 If set to True, all tracking data is saved to disk for further analysis (including csv
67 and txt files containing feature properties that are not returned in output cubes).
68 Default is False.
70 Returns
71 -------
72 tracking_cubes: iris.cube.CubeList
73 A list of iris cubes containing tracking data, including feature ID, lifetime,
74 and locations of initiating features.
76 Notes
77 -----
78 This operator uses the Simple-Track package to track features between timesteps. Simple-Track is a
79 data-agnostic, threshold-based object tracking algorithm for 2D data. Features are tracked between
80 consecutive frames of data by projecting feature fields onto common timeframes and matching
81 between them based on the degree of overlap. Matched features retain the same identification
82 between all tracked fields, while new features are assigned a unique label.
83 Thus, Simple-Track compiles comprehensive information about feature merging, splitting, accretion,
84 initiation and dissipation.
86 Currently outputs three cubes containing the following data:
87 "feature_id":
88 A 2D field containing the unique label assigned to each feature, which is retained
89 if the feature is tracked across multiple timesteps. This cube can be used as a mask
90 to identify the location of the tracked feature throughout the evaluation period.
91 "feature_lifetime":
92 A 2D field containing the lifetime of each feature in terms of the number of
93 timesteps it has been tracked for. This cube can be used to distinguish between
94 mature and fresh features.
95 "feature_init":
96 A 2D binary field indicating the location of newly initiated features at each timestep.
97 These features are identified as having a lifetime of 1 AND have initiated sufficiently
98 far from other, existing features that they are not considered to have spawned from them.
100 Links
101 ----------
102 .. https://github.com/ParaChute-UK/simple-track
104 Examples
105 --------
106 >>> tracking_cubes = feature.track(threshold=2)
107 >>> lifetime_cube = tracking_cubes.extract_cube("feature_lifetime")
108 # Plot the final timestep of lifetime cube. This will show
109 # the lifetime of features that have been tracked for multiple previous
110 # timesteps, as well as new features that have just been initiated.
111 >>> iplt.pcolormesh(lifetime_cube[-1,:,:],cmap=mpl.cm.bwr)
112 >>> plt.gca().coastlines('10m')
113 >>> plt.clim(-5,5)
114 >>> plt.colorbar()
115 >>> plt.show()
117 """
118 # Check that the input cube has horizontal coordinates of xy type, not latitude/longitude
119 _check_xy_coords(cube)
121 # Setup config
122 tracker_config = {
123 "FEATURE": {
124 "threshold": threshold,
125 "under_threshold": under_threshold,
126 "min_size": min_size,
127 },
128 "TRACKING": {
129 "retain_lifetime_on_split": retain_lifetime_on_split,
130 "overlap_nbhood": tracking_nbhood,
131 "overlap_threshold": overlap_threshold,
132 },
133 "OUTPUT": {
134 "save_data": save_data,
135 "experiment_name": "feature_tracking",
136 "path": f"{os.getcwd()}/tracking_data",
137 },
138 }
139 logging.debug(f"Tracker config: {tracker_config}")
141 # Get cube data into a dict to pass to Tracker
142 times = cube.coord("time").points
143 time_units = cube.coord("time").units
144 times_dt = [time_units.num2pydate(t) for t in times]
145 cube_dict = {
146 time: cube_slice.data
147 for time, cube_slice in zip(times_dt, cube.slices_over("time"), strict=True)
148 }
150 # Run tracking, returning Timeline object
151 timeline = Tracker(tracker_config).run(cube_dict)
152 logging.debug("Tracking completed")
154 # Use input cube as template to make returned cube
155 # By iterating over all cube times, this will ensure all data is present
156 # If a Frame at the given time is not contained in the timeline, error is raised
157 output_type_and_methods = {
158 "lifetime": {
159 "getter": "lifetime_field",
160 "cube_name": "feature_lifetime",
161 },
162 "feature": {
163 "getter": "feature_field",
164 "cube_name": "feature_id",
165 },
166 "init": {
167 "getter": "get_init_field",
168 "cube_name": "feature_init",
169 },
170 }
172 tracking_cubelist = iris.cube.CubeList()
173 for output_type in output_type_and_methods:
174 tracking_data = []
175 for time in times_dt:
176 frame = timeline.get_frame(time)
177 getter = getattr(frame, output_type_and_methods[output_type]["getter"])
178 if callable(getter):
179 tracking_data.append(getter())
180 else:
181 tracking_data.append(getter)
183 # Convert to numpy arrays
184 tracking_data = np.stack(tracking_data, axis=0)
186 # Create cubes
187 tracking_cube = cube.copy(data=tracking_data)
188 tracking_cube.long_name = output_type_and_methods[output_type]["cube_name"]
189 tracking_cube.standard_name = None
190 tracking_cube.var_name = None
191 tracking_cube.units = "1"
192 tracking_cubelist.append(tracking_cube)
194 return tracking_cubelist
197def _check_xy_coords(cube: iris.cube.Cube) -> None:
198 """Check that the input cube has horizontal coordinates of xy type, not latitude/longitude.
200 Parameters
201 ----------
202 cube: iris.cube.Cube
203 An iris cube containing horizontal coordinates.
205 Raises
206 ------
207 ValueError
208 If the input cube has horizontal coordinates of latitude/longitude type.
209 """
210 hzntl_coords = [
211 coord
212 for coord in cube.coords()
213 if iris.util.guess_coord_axis(coord) in ["X", "Y"]
214 ]
215 invalid_coord_names = ["latitude", "longitude", "grid_latitude", "grid_longitude"]
216 for coord in hzntl_coords:
217 if coord.name() in invalid_coord_names: 217 ↛ 218line 217 didn't jump to line 218 because the condition on line 217 was never true
218 raise ValueError(
219 f"Input cube {cube} has horizontal coordinate {coord}, "
220 "which is not of xy type. Please provide a cube with horizontal "
221 "coordinates of xy type."
222 )