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