From f479cb9c4745e17a05fe8ac8b8062a671c71311b Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Tue, 28 Apr 2026 17:14:33 +0100 Subject: [PATCH 01/33] added simple-track operator and tests Co-authored-by: Copilot --- .gitignore | 3 + src/CSET/operators/__init__.py | 5 + src/CSET/operators/feature.py | 186 ++++++++++++++++++++++++++++++++ tests/operators/test_feature.py | 129 ++++++++++++++++++++++ 4 files changed, 323 insertions(+) create mode 100755 src/CSET/operators/feature.py create mode 100644 tests/operators/test_feature.py diff --git a/.gitignore b/.gitignore index 5fe83c2a92..f0f7da2f9e 100644 --- a/.gitignore +++ b/.gitignore @@ -138,3 +138,6 @@ dmypy.json # NFS synchronisation files .nfs* + +#MacOS temp files +.DS_Store \ No newline at end of file diff --git a/src/CSET/operators/__init__.py b/src/CSET/operators/__init__.py index ef40fdf22a..5f8cb5c8a3 100644 --- a/src/CSET/operators/__init__.py +++ b/src/CSET/operators/__init__.py @@ -33,6 +33,7 @@ constraints, convection, ensembles, + feature, filters, imageprocessing, mesoscale, @@ -55,6 +56,7 @@ "convection", "ensembles", "execute_recipe", + "feature", "filters", "get_operator", "imageprocessing", @@ -104,7 +106,10 @@ def get_operator(name: str): name_sections = name.split(".") operator = CSET.operators for section in name_sections: + logging.debug(f"operator: {operator}") + logging.debug(f"section: {operator}") operator = getattr(operator, section) + if callable(operator): return operator else: diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py new file mode 100755 index 0000000000..a3cfe0c721 --- /dev/null +++ b/src/CSET/operators/feature.py @@ -0,0 +1,186 @@ +# © Crown copyright, Met Office (2022-2025) and CSET contributors. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +"""Operators for identifying and tracking features.""" + +import logging +import os + +import iris +import numpy as np +from simpletrack.track import Tracker + + +def track( + cube: iris.cube.Cube, + threshold: float, + under_threshold: bool = False, + min_size: int = 4, + retain_lifetime_on_split: bool = True, + tracking_nbhood: int = 5, + overlap_threshold: float = 0.3, + save_data: bool = False, +): + """Track features between subsequent timesteps. + + Parameters + ---------- + threshold: float + The threshold value for feature detection. + under_threshold: bool, optional + If set to True, features are identified where the data is below the threshold. + If set to False, features are identified where the data is above the threshold. + Default is False. + min_size: int, optional + The minimum number of contiguous grid points required for a feature to be tracked. + Default is 4. + retain_lifetime_on_split: bool, optional + If set to True, the lifetime of a feature is retained when it splits into + multiple features. If set to False, the lifetime is reset when a feature splits. + Default is True. + tracking_nbhood: int, optional + The size of the neighbourhood used for tracking features between timesteps. + This dictates the maximum pixel radius from a feature centroid at which new features could + reasonably be spawned. + Default is 5. + overlap_threshold: float, optional + The minimum overlap required between features in consecutive timesteps for + them to be considered the same feature. + Default is 0.3. + save_data: bool, optional + If set to True, all tracking data is saved to disk for further analysis (including csv + and txt files containing feature properties that are not returned in output cubes). + Default is False. + + Returns + ------- + tracking_cubes: iris.cube.CubeList + A list of iris cubes containing tracking data, including feauture ID, lifetime, + and locations of initiating features. + + Notes + ----- + This operator uses the Simple-Track package to track features between timesteps. Simple-Track is a + data-agnostic, threshold-based object tracking algorithm for 2D data. Features are tracked between + consecutive frames of data by projecting feature fields onto common timeframes and matching + between them based on the degree of overlap. Matched features retain the same identification + between all tracked fields, while new features are assigned a unique label. + Thus, Simple-Track compiles comprehensive information about feature merging, splitting, accretion, + initiation and dissipation. + + Currently outputs three cubes containing the following data: + "feature_id": + A 2D field containing the unique label assigned to each feature, which is retained + if the feature is tracked across multiple timesteps. This cube can be used as a mask + to identify the location of the tracked feature throughout the evaluation period. + "feature_lifetime": + A 2D field containing the lifetime of each feature in terms of the number of + timesteps it has been tracked for. This cube can be used to distinguish between + mature and fresh features. + "feature_init": + A 2D binary field indicating the location of newly initiated features at each timestep. + These features are identified as having a lifetime of 1 AND have initiated sufficiently + far from other, existing features that they are not considered to have spawed from them. + + Links + ---------- + .. https://github.com/ParaChute-UK/simple-track + + Examples + -------- + >>> tracking_cubes = feature.track(threshold=2) + >>> lifetime_cube = tracking_cubes.extract_cube("feature_lifetime") + # Plot the final timestep of lifetime cube. This will show + # the lifetime of features that have been tracked for multiple previous + # timesteps, as well as new features that have just been initiated. + >>> iplt.pcolormesh(lifetime_cube[-1,:,:],cmap=mpl.cm.bwr) + >>> plt.gca().coastlines('10m') + >>> plt.clim(-5,5) + >>> plt.colorbar() + >>> plt.show() + + """ + # Setup config + tracker_config = { + "FEATURE": { + "threshold": threshold, + "under_threshold": under_threshold, + "min_size": min_size, + }, + "TRACKING": { + "retain_lifetime_on_split": retain_lifetime_on_split, + "overlap_nbhood": tracking_nbhood, + "overlap_threshold": overlap_threshold, + }, + "OUTPUT": { + "save_data": save_data, + "experiment_name": "feature_tracking", + "path": f"{os.getcwd()}/tracking_data", + }, + } + logging.debug(f"Tracker config: {tracker_config}") + + # Get cube data into a dict to pass to Tracker + times = cube.coord("time").points + time_units = cube.coord("time").units + times_dt = [time_units.num2pydate(t) for t in times] + cube_dict = { + time: cube_slice.data + for time, cube_slice in zip(times_dt, cube.slices_over("time"), strict=True) + } + + # Run tracking, returning Timeline object + timeline = Tracker(tracker_config).run(cube_dict) + logging.debug("Tracking completed") + + # Use input cube as template to make returned cube + # By iterating over all cube times, this will ensure all data is present + # If a Frame at the given time is not contained in the timeline, error is raised + output_type_and_methods = { + "lifetime": { + "getter": "lifetime_field", + "cube_name": "feature_lifetime", + }, + "feature": { + "getter": "feature_field", + "cube_name": "feature_id", + }, + "init": { + "getter": "get_init_field", + "cube_name": "feature_init", + }, + } + + tracking_cubelist = iris.cube.CubeList() + for output_type in output_type_and_methods: + tracking_data = [] + for time in times_dt: + frame = timeline.get_frame(time) + getter = getattr(frame, output_type_and_methods[output_type]["getter"]) + if callable(getter): + tracking_data.append(getter()) + else: + tracking_data.append(getter) + + # Convert to numpy arrays + tracking_data = np.stack(tracking_data, axis=0) + + # Create cubes + tracking_cube = cube.copy(data=tracking_data) + tracking_cube.long_name = output_type_and_methods[output_type]["cube_name"] + tracking_cube.standard_name = None + tracking_cube.var_name = None + tracking_cube.units = "1" + tracking_cubelist.append(tracking_cube) + + return tracking_cubelist diff --git a/tests/operators/test_feature.py b/tests/operators/test_feature.py new file mode 100644 index 0000000000..4d2565a7ca --- /dev/null +++ b/tests/operators/test_feature.py @@ -0,0 +1,129 @@ +# © Crown copyright, Met Office (2022-2025) and CSET contributors. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +"""Tests for feature operators.""" + +import datetime as dt +import os + +import cf_units +import iris +import iris.coords +import iris.cube +import numpy as np +import pytest + +from CSET.operators import feature + + +@pytest.fixture +def feature_cube() -> iris.cube.Cube: + """Set up three timesteps of data and place into cube.""" + data_arr = np.zeros((3, 10, 10)) + data_arr[0, 2:6, 2:6] = 1 + data_arr[1, 3:7, 3:7] = 1 + data_arr[2, 4:8, 4:8] = 1 + + time_units = cf_units.Unit("days since 2000-01-01 00:00:00", calendar="gregorian") + time_start = dt.datetime(2010, 1, 1, 0, 0, 0) + time_dt_points = [time_start + dt.timedelta(minutes=5 * idx) for idx in range(3)] + time_points = [time_units.date2num(time_point) for time_point in time_dt_points] + time_coord = iris.coords.DimCoord( + points=time_points, standard_name="time", units=time_units + ) + + coord_system = iris.coord_systems.TransverseMercator( + latitude_of_projection_origin=55, longitude_of_central_meridian=0 + ) + coord_range = np.arange(0, 100, 10) + proj_y_coord = iris.coords.DimCoord( + points=coord_range, + standard_name="projection_y_coordinate", + var_name="projection_y_coordinate", + units="m", + coord_system=coord_system, + ) + proj_x_coord = iris.coords.DimCoord( + points=coord_range, + standard_name="projection_x_coordinate", + var_name="projection_x_coordinate", + units="m", + coord_system=coord_system, + ) + + proj_y_coord.guess_bounds() + proj_x_coord.guess_bounds() + + coords = (time_coord, proj_y_coord, proj_x_coord) + dim_coords_and_dims = [(coord, dim) for dim, coord in enumerate(coords)] + cube = iris.cube.Cube( + data=data_arr, + dim_coords_and_dims=dim_coords_and_dims, + long_name="Precipitation test", + ) + return cube + + +def test_tracking_valid(feature_cube) -> None: + """ + Test feature tracking returns same cube shape as input cube. + + Further tracking tests handled by Simple-Track dependency + """ + test_threshold = 0.5 + min_size = 1 + tracking_cubelist = feature.track( + feature_cube, threshold=test_threshold, min_size=min_size + ) + outputs = ["feature_lifetime", "feature_id", "feature_init"] + for output in outputs: + tracking_cube = tracking_cubelist.extract_cube(output) + assert tracking_cube.shape == feature_cube.shape + + +def test_tracking_lifetime_values(feature_cube) -> None: + """Test feature tracking returns expected lifetime values.""" + test_threshold = 0.5 + min_size = 1 + tracking_cubelist = feature.track( + feature_cube, threshold=test_threshold, min_size=min_size + ) + tracking_cube = tracking_cubelist.extract_cube("feature_lifetime") + # Check lifetime field values are expected, based on feature_cube data + for time_slice_idx in range(3): + expected_lifetime_field = np.where( + feature_cube.data[time_slice_idx] > test_threshold, time_slice_idx + 1, 0 + ) + actual_lifetime_field = tracking_cube.data[time_slice_idx] + np.testing.assert_array_equal(actual_lifetime_field, expected_lifetime_field) + + +def test_save_data(feature_cube, tmp_path) -> None: + """Test that tracking data is saved when save_data is True.""" + os.chdir(tmp_path) + test_threshold = 0.5 + min_size = 1 + feature.track( + feature_cube, + threshold=test_threshold, + min_size=min_size, + save_data=True, + ) + # Check expected lifetime field is created in output directory + output_directory = f"{tmp_path}/tracking_data" + expected_file = f"{output_directory}/lifetime_20100101_0000.field" + assert os.path.isfile(expected_file) + + # Check expected csv file is created in output directory + expected_file = f"{output_directory}/frame_20100101_0000.csv" + assert os.path.isfile(expected_file) From 4f56efb5da70fcb0bfe7fba3927eee8ba39b2a31 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Thu, 30 Apr 2026 12:36:13 +0100 Subject: [PATCH 02/33] added example tracking recipe --- .../example_feature_track.yaml | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100755 src/CSET/recipes/example_recipes/example_feature_track.yaml diff --git a/src/CSET/recipes/example_recipes/example_feature_track.yaml b/src/CSET/recipes/example_recipes/example_feature_track.yaml new file mode 100755 index 0000000000..67d67a0296 --- /dev/null +++ b/src/CSET/recipes/example_recipes/example_feature_track.yaml @@ -0,0 +1,26 @@ +category: Quick Look +title: Example running cell tracking and plotting spatial plots +description: | + Uses the feature.track operator to identify and track features in a cube with "time" series coordinate, + and then plots the lifetime of the identified features as a spatial plot. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS + + - operator: filters.filter_cubes + constraint: + operator: constraints.generate_var_constraint + varname: precipitation_flux + + - operator: feature.track + threshold: 3 + save_data: False # Whether to save raw tracking data for further analysis + + # Filter tracking cubelist to just one of "feature_lifetime", "feature_id" or "feature_init" + - operator: filters.filter_cubes + constraint: + operator: constraints.generate_var_constraint + varname: feature_lifetime + + - operator: plot.spatial_pcolormesh_plot From 4a1058ada87b636f2a2e32d3d906f5f02c4a7246 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Thu, 30 Apr 2026 12:37:18 +0100 Subject: [PATCH 03/33] added cell stats operator and changed nan behaviour in plot_histogram_series Co-authored-by: Copilot --- src/CSET/operators/feature.py | 144 ++++++++++++++++++++++++++++++++++ src/CSET/operators/plot.py | 4 +- 2 files changed, 146 insertions(+), 2 deletions(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index a3cfe0c721..6f35b9a021 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -184,3 +184,147 @@ def track( tracking_cubelist.append(tracking_cube) return tracking_cubelist + + +def cell_stats( + cube: iris.cube.Cube, + threshold: float, + under_threshold: bool = False, + min_size: int = 4, + save_data: bool = False, +): + """Identify features in each timestep and output statistics. + + Parameters + ---------- + threshold: float + The threshold value for feature detection. + under_threshold: bool, optional + If set to True, features are identified where the data is below the threshold. + If set to False, features are identified where the data is above the threshold. + Default is False. + min_size: int, optional + The minimum number of contiguous grid points required for a feature to be tracked. + Default is 4. + save_data: bool, optional + If set to True, all tracking data is saved to disk for further analysis (including csv + and txt files containing feature properties that are not returned in output cubes). + Default is False. + + Returns + ------- + cell_stats_cubes: iris.cube.CubeList + An iris CubeList containing "feature_size", "feature_mean", and "feature_max" cubes. + + Notes + ----- + This operator uses the Simple-Track package with tracking disabled to identify features + in each timestep and compile cell statistics. Outputs cubes containing feature size (number + of grid points), mean value within features, and maximum value within features. + + Links + ---------- + .. https://github.com/ParaChute-UK/simple-track + + Examples + -------- + >>> cell_stats_cubes = feature.cell_stats(threshold=2) + >>> feature_size_cube = cell_stats_cubes.extract_cube("feature_size") + >>> plt.hist(feature_size_cube[-1]) + >>> plt.show() + + """ + # Setup config + tracker_config = { + "FEATURE": { + "threshold": threshold, + "under_threshold": under_threshold, + "min_size": min_size, + }, + "OUTPUT": { + "save_data": save_data, + "experiment_name": "feature_tracking", + "path": f"{os.getcwd()}/cell-stats_data", + "skip_tracking": True, + }, + } + logging.debug(f"Tracker config: {tracker_config}") + + # Get cube data into a dict to pass to Tracker + times = cube.coord("time").points + time_units = cube.coord("time").units + times_dt = [time_units.num2pydate(t) for t in times] + cube_dict = { + time: cube_slice.data + for time, cube_slice in zip(times_dt, cube.slices_over("time"), strict=True) + } + + # Run tracking, returning Timeline object + timeline = Tracker(tracker_config).run(cube_dict) + logging.debug("Tracking completed") + + # Get feature data from each frame of data, append to list of lists + # before conversion to numpy array (which requires knowledge of max + # number of features before constructing array shape) + size_data, mean_data, max_data = [], [], [] + number_of_features = [] + for time in times_dt: + frame = timeline.get_frame(time) + features = frame.features + size_data.append([feature.get_size() for feature in features.values()]) + mean_data.append([feature.mean for feature in features.values()]) + max_data.append([feature.max for feature in features.values()]) + number_of_features.append(len(features)) + + # Pad data with NaNs to create arrays of consistent shape (max number of features across + # all timesteps) + arr_size = max(number_of_features) + + # Size data is integer, but we need to pad with NaNs, so fill with invalid value first + size_data = np.array( + [ + np.pad(sizes, (0, arr_size - len(sizes)), constant_values=-100) + for sizes in size_data + ], + dtype=float, + ) + size_data[size_data == -100] = np.nan + + # Mean and max data are already float, so can be padded with NaNs directly. + mean_data = np.array( + [ + np.pad(means, (0, arr_size - len(means)), constant_values=np.nan) + for means in mean_data + ] + ) + max_data = np.array( + [ + np.pad(maxs, (0, arr_size - len(maxs)), constant_values=np.nan) + for maxs in max_data + ] + ) + + # Create cubes + time_coord = cube.coord("time").copy() + feature_coord = iris.coords.DimCoord( + np.arange(arr_size), + long_name="feature_number", + var_name="feature_number", + units="1", + ) + coords = [time_coord, feature_coord] + coords_and_dims = [(coord, i) for i, coord in enumerate(coords)] + + cell_stats_cubelist = iris.cube.CubeList() + cube_names = ["feature_size", "feature_mean", "feature_max"] + data_arrays = [size_data, mean_data, max_data] + for cube_name, data_array in zip(cube_names, data_arrays, strict=True): + cell_stats_cube = iris.cube.Cube( + data_array, + long_name=cube_name, + units="1", + dim_coords_and_dims=coords_and_dims, + ) + cell_stats_cubelist.append(cell_stats_cube) + + return cell_stats_cubelist diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index 5ac8966a66..a88626181d 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -3074,8 +3074,8 @@ def plot_histogram_series( break if levels is None: - vmin = min(cb.data.min() for cb in cubes) - vmax = max(cb.data.max() for cb in cubes) + vmin = min(np.nanmin(cb.data) for cb in cubes) + vmax = max(np.nanmax(cb.data) for cb in cubes) # Make postage stamp plots if stamp_coordinate exists and has more than a # single point. If single_plot is True: From 40f87c21012bcba52d4bd6d7c499cf3c22b4f878 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Thu, 30 Apr 2026 12:39:26 +0100 Subject: [PATCH 04/33] added example cell_stats recipe --- .../example_feature_cell_stats.yaml | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100755 src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml diff --git a/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml b/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml new file mode 100755 index 0000000000..728c16cd11 --- /dev/null +++ b/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml @@ -0,0 +1,26 @@ +category: Quick Look +title: Example running cell_stats and plotting histograms +description: | + Uses the feature.cell_stats operator to calculate cell size, mean cell values and max cell values for + identified features. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS + + - operator: filters.filter_cubes + constraint: + operator: constraints.generate_var_constraint + varname: precipitation_flux + + - operator: feature.cell_stats + threshold: 3 + + # Filter tracking cubelist to one of "feature_size", "feature_mean" or "feature_max" + - operator: filters.filter_cubes + constraint: + operator: constraints.generate_var_constraint + varname: feature_max + + - operator: plot.plot_histogram_series + From 74cb58e71d8d0fc4e8ba6981fd5dde0f47b3a1ac Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Fri, 1 May 2026 14:15:05 +0100 Subject: [PATCH 05/33] log-log plotting for features, and frequency on y-axis rather than density --- src/CSET/operators/plot.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index a88626181d..562c2d5b72 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -1368,7 +1368,10 @@ def _plot_and_save_histogram_series( # Set default that histograms will produce probability density function # at each bin (integral over range sums to 1). - density = True + if "feature" in cubes[0].long_name: + density = False + else: + density = True for cube in iter_maybe(cubes): # Easier to check title (where var name originates) @@ -1400,6 +1403,10 @@ def _plot_and_save_histogram_series( np.max(bins), ) + if "feature" in cube.long_name: + ax.set_yscale("log") + ax.set_xscale("log") + # Reshape cube data into a single array to allow for a single histogram. # Otherwise we plot xdim histograms stacked. cube_data_1d = (cube.data).flatten() @@ -1432,6 +1439,8 @@ def _plot_and_save_histogram_series( ax.set_ylabel( f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14 ) + if "feature" in cubes[0].long_name: + ax.set_ylabel("Frequency", fontsize=14) ax.set_xlim(vmin, vmax) ax.tick_params(axis="both", labelsize=12) From e1070751fc8cce6920a0e604f90bbd89d329f423 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Fri, 1 May 2026 14:25:33 +0100 Subject: [PATCH 06/33] used same log bins for "surface_microphysical" as "feature" cubes - may need revisiting for mean and max --- src/CSET/operators/plot.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index 562c2d5b72..9e0cca1abd 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -1377,7 +1377,7 @@ def _plot_and_save_histogram_series( # Easier to check title (where var name originates) # than seeing if long names exist etc. # Exception case, where distribution better fits log scales/bins. - if "surface_microphysical" in title: + if "surface_microphysical" in title or "feature" in cube.long_name: if "amount" in title: # Compute histogram following Klingaman et al. (2017): ASoP bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100)) From f75e6d9c357aeceda20ecd6b519f55e142b5b56a Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Fri, 1 May 2026 15:51:51 +0100 Subject: [PATCH 07/33] added feature_effective_radius as alternative to feature_size Co-authored-by: Copilot --- src/CSET/operators/feature.py | 50 ++++++++++++++++--- .../example_feature_cell_stats.yaml | 6 ++- 2 files changed, 47 insertions(+), 9 deletions(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index 6f35b9a021..6dd32c510b 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -17,6 +17,8 @@ import os import iris +import iris.cube +import iris.util import numpy as np from simpletrack.track import Tracker @@ -191,6 +193,7 @@ def cell_stats( threshold: float, under_threshold: bool = False, min_size: int = 4, + feature_size_unit: str = "feature_effective_radius", save_data: bool = False, ): """Identify features in each timestep and output statistics. @@ -206,6 +209,11 @@ def cell_stats( min_size: int, optional The minimum number of contiguous grid points required for a feature to be tracked. Default is 4. + feature_size_unit: str, optional + The unit to define feature size output. Options are "feature_effective_radius" + (the radius of a circle with the same area as the feature, in km) or " + feature_grid_points" (the number of grid points contained in the feature). + Default is "feature_effective_radius". save_data: bool, optional If set to True, all tracking data is saved to disk for further analysis (including csv and txt files containing feature properties that are not returned in output cubes). @@ -314,15 +322,43 @@ def cell_stats( ) coords = [time_coord, feature_coord] coords_and_dims = [(coord, i) for i, coord in enumerate(coords)] - cell_stats_cubelist = iris.cube.CubeList() - cube_names = ["feature_size", "feature_mean", "feature_max"] - data_arrays = [size_data, mean_data, max_data] - for cube_name, data_array in zip(cube_names, data_arrays, strict=True): + + # Set cube properties + cube_properties = { + "feature_size": { + "data": size_data, + "long_name": feature_size_unit, + "units": "1", + }, + "feature_mean": {"data": mean_data, "long_name": "feature_mean", "units": 1}, + "feature_max": {"data": max_data, "long_name": "feature_max", "units": 1}, + } + if feature_size_unit == "feature_effective_radius": + # Conert feature size to effective radius, assuming uniform grid + # Guess coord representing horizontal grid (choose first available) + hzntl_coord = [ + coord + for coord in cube.coords() + if iris.util.guess_coord_axis(coord) in ["X", "Y"] + ][0] + logging.debug(f"Attempting to convert to effective radius using {hzntl_coord}") + # Convert to km if possible + hzntl_coord.convert_units("km") + # Naive grid spacing estimate, correct for regular grids, likely wildly + # inaccruate for LFRic/irregular grids + grid_spacing = np.abs(np.diff(hzntl_coord.points).mean()) + # Convert feature size in grid points to effective radius in km + size_data = np.sqrt(size_data * grid_spacing**2 / np.pi) + # Reset cube properties + cube_properties["feature_size"]["data"] = size_data + cube_properties["feature_size"]["units"] = "km" + + for cb_props in cube_properties.values(): cell_stats_cube = iris.cube.Cube( - data_array, - long_name=cube_name, - units="1", + data=cb_props["data"], + long_name=cb_props["long_name"], + units=cb_props["units"], dim_coords_and_dims=coords_and_dims, ) cell_stats_cubelist.append(cell_stats_cube) diff --git a/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml b/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml index 728c16cd11..300417442e 100755 --- a/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml +++ b/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml @@ -14,13 +14,15 @@ steps: varname: precipitation_flux - operator: feature.cell_stats + feature_size_unit: "feature_effective_radius" threshold: 3 - # Filter tracking cubelist to one of "feature_size", "feature_mean" or "feature_max" + # Filter tracking cubelist to one of "feature_mean", "feature_max", or + # "feature_size/feature_effective_radius" (depending on feature_size_unit set above), - operator: filters.filter_cubes constraint: operator: constraints.generate_var_constraint - varname: feature_max + varname: feature_effective_radius - operator: plot.plot_histogram_series From e928ba892cf757b3875b9e5d42c9531a50e1aaac Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Fri, 1 May 2026 15:56:26 +0100 Subject: [PATCH 08/33] additional comments --- src/CSET/operators/feature.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index 6dd32c510b..c581517412 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -288,7 +288,8 @@ def cell_stats( # all timesteps) arr_size = max(number_of_features) - # Size data is integer, but we need to pad with NaNs, so fill with invalid value first + # Size data is integer, but we need to pad with NaNs (which is a float), so fill + # with invalid value first size_data = np.array( [ np.pad(sizes, (0, arr_size - len(sizes)), constant_values=-100) @@ -344,6 +345,7 @@ def cell_stats( ][0] logging.debug(f"Attempting to convert to effective radius using {hzntl_coord}") # Convert to km if possible + # TODO: fall back to feature_size if this conversion fails hzntl_coord.convert_units("km") # Naive grid spacing estimate, correct for regular grids, likely wildly # inaccruate for LFRic/irregular grids @@ -354,6 +356,7 @@ def cell_stats( cube_properties["feature_size"]["data"] = size_data cube_properties["feature_size"]["units"] = "km" + # Populate cubelist for cb_props in cube_properties.values(): cell_stats_cube = iris.cube.Cube( data=cb_props["data"], From 7d89cbb7eb788554f52efd292ed827fceb163d1c Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Fri, 1 May 2026 17:30:33 +0100 Subject: [PATCH 09/33] added temporary feature cbar definitons (limits to be set dynamically) --- src/CSET/operators/_colorbar_definition.json | 44 ++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/src/CSET/operators/_colorbar_definition.json b/src/CSET/operators/_colorbar_definition.json index 0c075b7381..581554c894 100644 --- a/src/CSET/operators/_colorbar_definition.json +++ b/src/CSET/operators/_colorbar_definition.json @@ -163,6 +163,50 @@ "max": 1.1, "min": 0.5 }, + "feature_id": { + "cmap": "viridis", + "levels": [ + 1, + 50, + 100, + 150, + 200, + 250, + 300, + 350, + 400 + ], + "ymax": 1.0, + "ymin": 0.0 + }, + "feature_lifetime": { + "cmap": "YlGnBu", + "levels": [ + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12 + ], + "ymax": 1.0, + "ymin": 0.0 + }, + "feature_init": { + "cmap": "Blues", + "levels": [ + 0.5, + 1 + ], + "ymax": 1.0, + "ymin": 0.0 + }, "fog_fraction_at_screen_level": { "cmap": "viridis", "max": 1, From 98f39d96967f2ac97f06cf3af867df8fa8141cdd Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Fri, 1 May 2026 17:58:45 +0100 Subject: [PATCH 10/33] added simple-track dependency to pyproject and env.yml --- pyproject.toml | 1 + requirements/environment.yml | 1 + 2 files changed, 2 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index e236720400..b5716c1b5f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -33,6 +33,7 @@ dependencies = [ "scipy", "scikit-image", "dask", + "simple-track", ] [project.urls] diff --git a/requirements/environment.yml b/requirements/environment.yml index 4485c1e61b..2d6f9157ba 100644 --- a/requirements/environment.yml +++ b/requirements/environment.yml @@ -17,6 +17,7 @@ dependencies: - scikit-image # For image processing techniques. - dask-core # Dask with minimal dependencies. - proj = 9.7.1 # Newer versions break plotting, see issue #2052. + - simple-track # For feature tracking and cell stats operators # Build dependencies - setuptools>=64 From 730c7aa1dc1be06b0e2009e74dd097431e898769 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Fri, 1 May 2026 18:00:19 +0100 Subject: [PATCH 11/33] added simple-track to dependencies --- pyproject.toml | 1 + requirements/environment.yml | 1 + 2 files changed, 2 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index e236720400..b5716c1b5f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -33,6 +33,7 @@ dependencies = [ "scipy", "scikit-image", "dask", + "simple-track", ] [project.urls] diff --git a/requirements/environment.yml b/requirements/environment.yml index af5ba93ac6..2fbcaa3cbc 100644 --- a/requirements/environment.yml +++ b/requirements/environment.yml @@ -16,6 +16,7 @@ dependencies: - scipy - scikit-image # For image processing techniques. - dask + - simple-track # Build dependencies - setuptools>=64 From a3bce7723a9a190c9b88f98a0dfb153612c56d5c Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Fri, 8 May 2026 13:25:03 +0100 Subject: [PATCH 12/33] added cube metadata copying to feature cubes, and error catching for feature_effective_radius conversion --- src/CSET/operators/feature.py | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index c581517412..a9ebbcec3d 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -346,7 +346,12 @@ def cell_stats( logging.debug(f"Attempting to convert to effective radius using {hzntl_coord}") # Convert to km if possible # TODO: fall back to feature_size if this conversion fails - hzntl_coord.convert_units("km") + # TODO: this will still add the name "feature_effective_radius" to the feature_size + # cube, even though the data is actually pixel size. Figure out elegant solution. + try: + hzntl_coord.convert_units("km") + except: + pass # Naive grid spacing estimate, correct for regular grids, likely wildly # inaccruate for LFRic/irregular grids grid_spacing = np.abs(np.diff(hzntl_coord.points).mean()) @@ -356,6 +361,25 @@ def cell_stats( cube_properties["feature_size"]["data"] = size_data cube_properties["feature_size"]["units"] = "km" + # Get list of coords to copy from input cube to output cubes + copyable_coord_names = [ + "realization", + "hour", + "forecast_period", + "forecast_reference_time", + "model_name", + ] + input_cube_coord_names = [] + for coord in cube.coords(): + input_cube_coord_names.append(coord.standard_name) + input_cube_coord_names.append(coord.long_name) + + coords_to_copy = [ + coord_name + for coord_name in copyable_coord_names + if coord_name in input_cube_coord_names + ] + # Populate cubelist for cb_props in cube_properties.values(): cell_stats_cube = iris.cube.Cube( @@ -364,6 +388,11 @@ def cell_stats( units=cb_props["units"], dim_coords_and_dims=coords_and_dims, ) + # Add other metadata from input cube + for coord_name in coords_to_copy: + cell_stats_cube.add_aux_coord(cube.coord(coord_name).copy()) + + # Add to cubelist cell_stats_cubelist.append(cell_stats_cube) return cell_stats_cubelist From c52c4302851076a8a75d79fe62c44d957d0d7357 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Fri, 8 May 2026 13:31:04 +0100 Subject: [PATCH 13/33] added cset_comparison_base to list of coord metadata to copy to feature cube (may not be necessary?) --- src/CSET/operators/feature.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index a9ebbcec3d..a33a385f4d 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -368,6 +368,7 @@ def cell_stats( "forecast_period", "forecast_reference_time", "model_name", + "cset_comparison_base", ] input_cube_coord_names = [] for coord in cube.coords(): From 54c3db580e803ee70f3d2a7aa86b26bd97e63e1e Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Wed, 13 May 2026 11:06:41 +0100 Subject: [PATCH 14/33] copying coord to cell_stats_cube now copies dimension information (assuming same cube structure) --- src/CSET/operators/feature.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index a33a385f4d..adefa9bf32 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -391,7 +391,13 @@ def cell_stats( ) # Add other metadata from input cube for coord_name in coords_to_copy: - cell_stats_cube.add_aux_coord(cube.coord(coord_name).copy()) + coord = cube.coord(coord_name).copy() + # Check if this coord represents a dimension of data + dims = cube.coord_dims(coord) + if len(dims) > 1: + cell_stats_cube.add_aux_coord(coord, dims) + else: + cell_stats_cube.add_aux_coord(coord) # Add to cubelist cell_stats_cubelist.append(cell_stats_cube) From b8b1f6a3d820b79859a9e509d0d7a22a5f036871 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Tue, 19 May 2026 13:26:31 +0100 Subject: [PATCH 15/33] removed unnecessary logging, added set_under option to feature plot --- src/CSET/operators/__init__.py | 2 -- src/CSET/operators/plot.py | 3 +++ 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/CSET/operators/__init__.py b/src/CSET/operators/__init__.py index 5f8cb5c8a3..44e65cd2ba 100644 --- a/src/CSET/operators/__init__.py +++ b/src/CSET/operators/__init__.py @@ -106,8 +106,6 @@ def get_operator(name: str): name_sections = name.split(".") operator = CSET.operators for section in name_sections: - logging.debug(f"operator: {operator}") - logging.debug(f"section: {operator}") operator = getattr(operator, section) if callable(operator): diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index c6a44e09bb..48a62582e3 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -647,6 +647,9 @@ def _plot_and_save_spatial_plot( # Specify the color bar cmap, levels, norm = _colorbar_map_levels(cube) + if "feature" in cube.long_name: + cmap.set_under("white") + # If overplotting, set required colorbars if overlay_cube: over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube) From dda9d48c57f747bbe0e8d445a7186f7ad72b63e0 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Tue, 19 May 2026 15:04:53 +0100 Subject: [PATCH 16/33] removed feature_effective_radius support for now, added support for multiple models in cell_stats operator --- src/CSET/operators/feature.py | 331 +++++++++--------- .../example_feature_cell_stats.yaml | 19 +- 2 files changed, 182 insertions(+), 168 deletions(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index adefa9bf32..f9ce014d18 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -22,6 +22,8 @@ import numpy as np from simpletrack.track import Tracker +from CSET._common import iter_maybe + def track( cube: iris.cube.Cube, @@ -67,7 +69,7 @@ def track( Returns ------- tracking_cubes: iris.cube.CubeList - A list of iris cubes containing tracking data, including feauture ID, lifetime, + A list of iris cubes containing tracking data, including feature ID, lifetime, and locations of initiating features. Notes @@ -92,7 +94,7 @@ def track( "feature_init": A 2D binary field indicating the location of newly initiated features at each timestep. These features are identified as having a lifetime of 1 AND have initiated sufficiently - far from other, existing features that they are not considered to have spawed from them. + far from other, existing features that they are not considered to have spawned from them. Links ---------- @@ -189,7 +191,7 @@ def track( def cell_stats( - cube: iris.cube.Cube, + cubes: iris.cube.Cube | iris.cube.CubeList, threshold: float, under_threshold: bool = False, min_size: int = 4, @@ -242,164 +244,179 @@ def cell_stats( >>> plt.show() """ - # Setup config - tracker_config = { - "FEATURE": { - "threshold": threshold, - "under_threshold": under_threshold, - "min_size": min_size, - }, - "OUTPUT": { - "save_data": save_data, - "experiment_name": "feature_tracking", - "path": f"{os.getcwd()}/cell-stats_data", - "skip_tracking": True, - }, - } - logging.debug(f"Tracker config: {tracker_config}") - - # Get cube data into a dict to pass to Tracker - times = cube.coord("time").points - time_units = cube.coord("time").units - times_dt = [time_units.num2pydate(t) for t in times] - cube_dict = { - time: cube_slice.data - for time, cube_slice in zip(times_dt, cube.slices_over("time"), strict=True) - } + # Check inputs + cubes = iter_maybe(cubes) + cell_stats_cubelist = iris.cube.CubeList() - # Run tracking, returning Timeline object - timeline = Tracker(tracker_config).run(cube_dict) - logging.debug("Tracking completed") + # Run tracking on all input data + for cube in cubes: + model_name = cube.attributes.get("model_name", None) + # Setup config + tracker_config = { + "FEATURE": { + "threshold": threshold, + "under_threshold": under_threshold, + "min_size": min_size, + }, + "OUTPUT": { + "save_data": save_data, + "experiment_name": "feature_tracking", + "path": f"{os.getcwd()}/{model_name}/cell-stats_data", + "skip_tracking": True, + }, + } + logging.debug(f"Tracker config: {tracker_config}") + + # Get cube data into a dict to pass to Tracker + times = cube.coord("time").points + time_units = cube.coord("time").units + times_dt = [time_units.num2pydate(t) for t in times] + cube_dict = { + time: cube_slice.data + for time, cube_slice in zip(times_dt, cube.slices_over("time"), strict=True) + } + + # Run tracking, returning Timeline object + timeline = Tracker(tracker_config).run(cube_dict) + + logging.debug(f"Tracking completed for {model_name}") + + # Get feature data from each frame of data, append to list of lists + # before conversion to numpy array (which requires knowledge of max + # number of features before constructing array shape) + size_data, mean_data, max_data = [], [], [] + number_of_features = [] + for time in times_dt: + frame = timeline.get_frame(time) + features = frame.features + size_data.append([feature.get_size() for feature in features.values()]) + mean_data.append([feature.mean for feature in features.values()]) + max_data.append([feature.max for feature in features.values()]) + number_of_features.append(len(features)) + + # Pad data with NaNs to create arrays of consistent shape (max number of features across + # all timesteps) + arr_size = max(number_of_features) + + # Size data is integer, but we need to pad with NaNs (which is a float), so fill + # with invalid value first + size_data = np.array( + [ + np.pad(sizes, (0, arr_size - len(sizes)), constant_values=-100) + for sizes in size_data + ], + dtype=float, + ) + size_data[size_data == -100] = np.nan + + # Mean and max data are already float, so can be padded with NaNs directly. + mean_data = np.array( + [ + np.pad(means, (0, arr_size - len(means)), constant_values=np.nan) + for means in mean_data + ] + ) + max_data = np.array( + [ + np.pad(maxs, (0, arr_size - len(maxs)), constant_values=np.nan) + for maxs in max_data + ] + ) - # Get feature data from each frame of data, append to list of lists - # before conversion to numpy array (which requires knowledge of max - # number of features before constructing array shape) - size_data, mean_data, max_data = [], [], [] - number_of_features = [] - for time in times_dt: - frame = timeline.get_frame(time) - features = frame.features - size_data.append([feature.get_size() for feature in features.values()]) - mean_data.append([feature.mean for feature in features.values()]) - max_data.append([feature.max for feature in features.values()]) - number_of_features.append(len(features)) - - # Pad data with NaNs to create arrays of consistent shape (max number of features across - # all timesteps) - arr_size = max(number_of_features) - - # Size data is integer, but we need to pad with NaNs (which is a float), so fill - # with invalid value first - size_data = np.array( - [ - np.pad(sizes, (0, arr_size - len(sizes)), constant_values=-100) - for sizes in size_data - ], - dtype=float, - ) - size_data[size_data == -100] = np.nan - - # Mean and max data are already float, so can be padded with NaNs directly. - mean_data = np.array( - [ - np.pad(means, (0, arr_size - len(means)), constant_values=np.nan) - for means in mean_data + # Create cubes + time_coord = cube.coord("time").copy() + feature_coord = iris.coords.DimCoord( + np.arange(arr_size), + long_name="feature_number", + var_name="feature_number", + units="1", + ) + coords = [time_coord, feature_coord] + coords_and_dims = [(coord, i) for i, coord in enumerate(coords)] + + # Set cube properties + cube_properties = { + "feature_size": { + "data": size_data, + "long_name": feature_size_unit, + "units": "1", + }, + "feature_mean": { + "data": mean_data, + "long_name": "feature_mean", + "units": 1, + }, + "feature_max": {"data": max_data, "long_name": "feature_max", "units": 1}, + } + # TODO: need to consider better implementation which doesn't overwrite feature_size cube + # TODO: consider iris.analysis.cartography.area_weights() or hzntl_coord.is_regular() + # to check for regular spacing + # TODO: Check for hzntl_coords is empty list and handle elegantly + # if feature_size_unit == "feature_effective_radius": + # # Convert feature size to effective radius, assuming uniform grid + # # Guess coord representing horizontal grid (choose first available) + # hzntl_coord = [ + # coord + # for coord in cube.coords() + # if iris.util.guess_coord_axis(coord) in ["X", "Y"] + # ][0] + # logging.debug(f"Attempting to convert to effective radius using {hzntl_coord}") + # # Convert to km if possible + # # TODO: fall back to feature_size if this conversion fails + # # TODO: this will still add the name "feature_effective_radius" to the feature_size + # # cube, even though the data is actually pixel size. Figure out elegant solution. + # try: + # hzntl_coord.convert_units("km") + # except: + # pass + # # Naive grid spacing estimate, correct for regular grids, likely wildly + # # inaccruate for LFRic/irregular grids + # grid_spacing = np.abs(np.diff(hzntl_coord.points).mean()) + # # Convert feature size in grid points to effective radius in km + # size_data = np.sqrt(size_data * grid_spacing**2 / np.pi) + # # Reset cube properties + # cube_properties["feature_size"]["data"] = size_data + # cube_properties["feature_size"]["units"] = "km" + + # Get list of coords to copy from input cube to output cubes + copyable_coord_names = [ + "realization", + "hour", + "forecast_period", + "forecast_reference_time", + "model_name", + "cset_comparison_base", ] - ) - max_data = np.array( - [ - np.pad(maxs, (0, arr_size - len(maxs)), constant_values=np.nan) - for maxs in max_data + input_cube_coord_names = [] + for coord in cube.coords(): + input_cube_coord_names.append(coord.standard_name) + input_cube_coord_names.append(coord.long_name) + + coords_to_copy = [ + coord_name + for coord_name in copyable_coord_names + if coord_name in input_cube_coord_names ] - ) - - # Create cubes - time_coord = cube.coord("time").copy() - feature_coord = iris.coords.DimCoord( - np.arange(arr_size), - long_name="feature_number", - var_name="feature_number", - units="1", - ) - coords = [time_coord, feature_coord] - coords_and_dims = [(coord, i) for i, coord in enumerate(coords)] - cell_stats_cubelist = iris.cube.CubeList() - - # Set cube properties - cube_properties = { - "feature_size": { - "data": size_data, - "long_name": feature_size_unit, - "units": "1", - }, - "feature_mean": {"data": mean_data, "long_name": "feature_mean", "units": 1}, - "feature_max": {"data": max_data, "long_name": "feature_max", "units": 1}, - } - if feature_size_unit == "feature_effective_radius": - # Conert feature size to effective radius, assuming uniform grid - # Guess coord representing horizontal grid (choose first available) - hzntl_coord = [ - coord - for coord in cube.coords() - if iris.util.guess_coord_axis(coord) in ["X", "Y"] - ][0] - logging.debug(f"Attempting to convert to effective radius using {hzntl_coord}") - # Convert to km if possible - # TODO: fall back to feature_size if this conversion fails - # TODO: this will still add the name "feature_effective_radius" to the feature_size - # cube, even though the data is actually pixel size. Figure out elegant solution. - try: - hzntl_coord.convert_units("km") - except: - pass - # Naive grid spacing estimate, correct for regular grids, likely wildly - # inaccruate for LFRic/irregular grids - grid_spacing = np.abs(np.diff(hzntl_coord.points).mean()) - # Convert feature size in grid points to effective radius in km - size_data = np.sqrt(size_data * grid_spacing**2 / np.pi) - # Reset cube properties - cube_properties["feature_size"]["data"] = size_data - cube_properties["feature_size"]["units"] = "km" - - # Get list of coords to copy from input cube to output cubes - copyable_coord_names = [ - "realization", - "hour", - "forecast_period", - "forecast_reference_time", - "model_name", - "cset_comparison_base", - ] - input_cube_coord_names = [] - for coord in cube.coords(): - input_cube_coord_names.append(coord.standard_name) - input_cube_coord_names.append(coord.long_name) - - coords_to_copy = [ - coord_name - for coord_name in copyable_coord_names - if coord_name in input_cube_coord_names - ] - - # Populate cubelist - for cb_props in cube_properties.values(): - cell_stats_cube = iris.cube.Cube( - data=cb_props["data"], - long_name=cb_props["long_name"], - units=cb_props["units"], - dim_coords_and_dims=coords_and_dims, - ) - # Add other metadata from input cube - for coord_name in coords_to_copy: - coord = cube.coord(coord_name).copy() - # Check if this coord represents a dimension of data - dims = cube.coord_dims(coord) - if len(dims) > 1: - cell_stats_cube.add_aux_coord(coord, dims) - else: - cell_stats_cube.add_aux_coord(coord) - # Add to cubelist - cell_stats_cubelist.append(cell_stats_cube) + # Populate cubelist + for cb_props in cube_properties.values(): + cell_stats_cube = iris.cube.Cube( + data=cb_props["data"], + long_name=cb_props["long_name"], + units=cb_props["units"], + dim_coords_and_dims=coords_and_dims, + ) + # Add other metadata from input cube + for coord_name in coords_to_copy: + coord = cube.coord(coord_name).copy() + # Check if this coord represents a dimension of data + dims = cube.coord_dims(coord) + if len(dims) > 1: + cell_stats_cube.add_aux_coord(coord, dims) + else: + cell_stats_cube.add_aux_coord(coord) + + # Add to cubelist + cell_stats_cubelist.append(cell_stats_cube) return cell_stats_cubelist diff --git a/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml b/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml index 300417442e..2f21438eb3 100755 --- a/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml +++ b/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml @@ -1,28 +1,25 @@ category: Quick Look title: Example running cell_stats and plotting histograms description: | - Uses the feature.cell_stats operator to calculate cell size, mean cell values and max cell values for - identified features. + Uses the feature.cell_stats operator to calculate cell size, mean cell values and max cell values for + identified features. steps: - operator: read.read_cubes file_paths: $INPUT_PATHS - - operator: filters.filter_cubes + - operator: filters.filter_multiple_cubes constraint: operator: constraints.generate_var_constraint - varname: precipitation_flux + varname: $VARNAME - operator: feature.cell_stats - feature_size_unit: "feature_effective_radius" threshold: 3 - # Filter tracking cubelist to one of "feature_mean", "feature_max", or - # "feature_size/feature_effective_radius" (depending on feature_size_unit set above), - - operator: filters.filter_cubes + # Filter tracking cubelist to one of "feature_mean", "feature_max", or "feature_size" + - operator: filters.filter_multiple_cubes constraint: operator: constraints.generate_var_constraint - varname: feature_effective_radius - - - operator: plot.plot_histogram_series + varname: feature_size + - operator: plot.plot_histogram_series From 65764a747af5ada054c250c700aa7be2512f17de Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Tue, 19 May 2026 15:52:42 +0100 Subject: [PATCH 17/33] changed example cell_stats recipe to include MODEL_NAME input arg, fixed removal of feature_effective_radius, input cube attrs now copied to cell_stats_cube --- src/CSET/operators/feature.py | 16 +++++++++------- .../example_feature_cell_stats.yaml | 1 + 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index f9ce014d18..006cdc7d38 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -195,7 +195,6 @@ def cell_stats( threshold: float, under_threshold: bool = False, min_size: int = 4, - feature_size_unit: str = "feature_effective_radius", save_data: bool = False, ): """Identify features in each timestep and output statistics. @@ -211,11 +210,11 @@ def cell_stats( min_size: int, optional The minimum number of contiguous grid points required for a feature to be tracked. Default is 4. - feature_size_unit: str, optional - The unit to define feature size output. Options are "feature_effective_radius" - (the radius of a circle with the same area as the feature, in km) or " - feature_grid_points" (the number of grid points contained in the feature). - Default is "feature_effective_radius". + # feature_size_unit: str, optional + # The unit to define feature size output. Options are "feature_effective_radius" + # (the radius of a circle with the same area as the feature, in km) or " + # feature_grid_points" (the number of grid points contained in the feature). + # Default is "feature_effective_radius". save_data: bool, optional If set to True, all tracking data is saved to disk for further analysis (including csv and txt files containing feature properties that are not returned in output cubes). @@ -338,7 +337,7 @@ def cell_stats( cube_properties = { "feature_size": { "data": size_data, - "long_name": feature_size_unit, + "long_name": "feature_size", "units": "1", }, "feature_mean": { @@ -416,6 +415,9 @@ def cell_stats( else: cell_stats_cube.add_aux_coord(coord) + # Copy over attributes + cell_stats_cube.attributes = cube.attributes + # Add to cubelist cell_stats_cubelist.append(cell_stats_cube) diff --git a/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml b/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml index 2f21438eb3..9acf8e3811 100755 --- a/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml +++ b/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml @@ -7,6 +7,7 @@ description: | steps: - operator: read.read_cubes file_paths: $INPUT_PATHS + model_names: $MODEL_NAME - operator: filters.filter_multiple_cubes constraint: From 7ed0592c66622800f4767c9a3f3887764057a0ca Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Wed, 27 May 2026 13:36:14 +0100 Subject: [PATCH 18/33] fixed bug for assigning coord dims --- src/CSET/operators/feature.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index 006cdc7d38..6d9e80dea3 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -410,7 +410,7 @@ def cell_stats( coord = cube.coord(coord_name).copy() # Check if this coord represents a dimension of data dims = cube.coord_dims(coord) - if len(dims) > 1: + if len(dims) > 0: cell_stats_cube.add_aux_coord(coord, dims) else: cell_stats_cube.add_aux_coord(coord) From 8554d53ecc02418edc6e2034896cad5f5300bce0 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Tue, 16 Jun 2026 14:06:22 +0100 Subject: [PATCH 19/33] added feature effective radius calculation --- src/CSET/operators/feature.py | 119 ++++++++++++++++++++-------------- 1 file changed, 72 insertions(+), 47 deletions(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index 6d9e80dea3..0eba0d72fe 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -210,11 +210,6 @@ def cell_stats( min_size: int, optional The minimum number of contiguous grid points required for a feature to be tracked. Default is 4. - # feature_size_unit: str, optional - # The unit to define feature size output. Options are "feature_effective_radius" - # (the radius of a circle with the same area as the feature, in km) or " - # feature_grid_points" (the number of grid points contained in the feature). - # Default is "feature_effective_radius". save_data: bool, optional If set to True, all tracking data is saved to disk for further analysis (including csv and txt files containing feature properties that are not returned in output cubes). @@ -223,13 +218,15 @@ def cell_stats( Returns ------- cell_stats_cubes: iris.cube.CubeList - An iris CubeList containing "feature_size", "feature_mean", and "feature_max" cubes. + An iris CubeList containing "feature_size", "feature_effective_radius", "feature_mean", + and "feature_max" cubes. Notes ----- This operator uses the Simple-Track package with tracking disabled to identify features in each timestep and compile cell statistics. Outputs cubes containing feature size (number - of grid points), mean value within features, and maximum value within features. + of grid points), effective radius (in km), mean value within features, and maximum + value within features. Links ---------- @@ -308,6 +305,11 @@ def cell_stats( ) size_data[size_data == -100] = np.nan + # Get effective radius from feature size, using horizontal coordinate of input cube to estimate grid spacing + effective_radius_data = _get_effective_radius_from_feature_size( + size_data=size_data, cube_with_hzntl_coord=cube + ) + # Mean and max data are already float, so can be padded with NaNs directly. mean_data = np.array( [ @@ -322,17 +324,6 @@ def cell_stats( ] ) - # Create cubes - time_coord = cube.coord("time").copy() - feature_coord = iris.coords.DimCoord( - np.arange(arr_size), - long_name="feature_number", - var_name="feature_number", - units="1", - ) - coords = [time_coord, feature_coord] - coords_and_dims = [(coord, i) for i, coord in enumerate(coords)] - # Set cube properties cube_properties = { "feature_size": { @@ -346,36 +337,23 @@ def cell_stats( "units": 1, }, "feature_max": {"data": max_data, "long_name": "feature_max", "units": 1}, + "feature_effective_radius": { + "data": effective_radius_data, + "long_name": "feature_effective_radius", + "units": "km", + }, } - # TODO: need to consider better implementation which doesn't overwrite feature_size cube - # TODO: consider iris.analysis.cartography.area_weights() or hzntl_coord.is_regular() - # to check for regular spacing - # TODO: Check for hzntl_coords is empty list and handle elegantly - # if feature_size_unit == "feature_effective_radius": - # # Convert feature size to effective radius, assuming uniform grid - # # Guess coord representing horizontal grid (choose first available) - # hzntl_coord = [ - # coord - # for coord in cube.coords() - # if iris.util.guess_coord_axis(coord) in ["X", "Y"] - # ][0] - # logging.debug(f"Attempting to convert to effective radius using {hzntl_coord}") - # # Convert to km if possible - # # TODO: fall back to feature_size if this conversion fails - # # TODO: this will still add the name "feature_effective_radius" to the feature_size - # # cube, even though the data is actually pixel size. Figure out elegant solution. - # try: - # hzntl_coord.convert_units("km") - # except: - # pass - # # Naive grid spacing estimate, correct for regular grids, likely wildly - # # inaccruate for LFRic/irregular grids - # grid_spacing = np.abs(np.diff(hzntl_coord.points).mean()) - # # Convert feature size in grid points to effective radius in km - # size_data = np.sqrt(size_data * grid_spacing**2 / np.pi) - # # Reset cube properties - # cube_properties["feature_size"]["data"] = size_data - # cube_properties["feature_size"]["units"] = "km" + + # Create cubes + time_coord = cube.coord("time").copy() + feature_coord = iris.coords.DimCoord( + np.arange(arr_size), + long_name="feature_number", + var_name="feature_number", + units="1", + ) + coords = [time_coord, feature_coord] + coords_and_dims = [(coord, i) for i, coord in enumerate(coords)] # Get list of coords to copy from input cube to output cubes copyable_coord_names = [ @@ -422,3 +400,50 @@ def cell_stats( cell_stats_cubelist.append(cell_stats_cube) return cell_stats_cubelist + + +def _get_effective_radius_from_feature_size( + size_data: np.ndarray, cube_with_hzntl_coord: iris.cube.Cube +) -> np.ndarray: + """Convert feature size in grid points to effective radius in km. + + Parameters + ---------- + size_data: np.ndarray + An array containing "feature_size" data, in units of grid points. + cube_with_hzntl_coord: iris.cube.Cube + An iris cube containing a horizontal coordinate, which is used to + estimate the grid spacing for the effective radius calculation. + + Returns + ------- + effective_radii_data: np.ndarray + An array containing "feature_effective_radius" data, in units of km. + + Notes + ----- + This function assumes that the input cube has a horizontal coordinate system that is regular and + that the grid spacing can be estimated from the horizontal coordinates. The effective radius is + calculated as the radius of a circle with the same area as the feature size in grid points. + + """ + # Guess coord representing horizontal grid (choose first available) + hzntl_coord = [ + coord + for coord in cube_with_hzntl_coord.coords() + if iris.util.guess_coord_axis(coord) in ["X", "Y"] + ][0] + logging.debug(f"Attempting to convert to effective radius using {hzntl_coord}") + + # Check coordinate is regular, but only warn if not, this is a naive estimate + # and will be inaccurate for irregular grids + if not iris.util.is_regular(hzntl_coord): + logging.warning( + f"Horizontal coordinate {hzntl_coord} is not regular. " + "Effective radius calculation may be inaccurate." + ) + + grid_spacing = np.abs(np.mean(np.diff(hzntl_coord.points))) + effective_radii_data = np.sqrt(size_data * grid_spacing**2 / np.pi) + + return effective_radii_data From 963d808f87599d3ef74c0fad3132e223c1b18085 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Tue, 16 Jun 2026 15:01:33 +0100 Subject: [PATCH 20/33] modularised cell_stats operator, added test --- src/CSET/operators/feature.py | 252 ++++++++++++++++++++------------ tests/operators/test_feature.py | 41 +++++- 2 files changed, 195 insertions(+), 98 deletions(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index 0eba0d72fe..31a36e9b1c 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -20,6 +20,7 @@ import iris.cube import iris.util import numpy as np +from simpletrack.frame import Timeline from simpletrack.track import Tracker from CSET._common import iter_maybe @@ -274,62 +275,24 @@ def cell_stats( # Run tracking, returning Timeline object timeline = Tracker(tracker_config).run(cube_dict) - logging.debug(f"Tracking completed for {model_name}") - # Get feature data from each frame of data, append to list of lists - # before conversion to numpy array (which requires knowledge of max - # number of features before constructing array shape) - size_data, mean_data, max_data = [], [], [] - number_of_features = [] - for time in times_dt: - frame = timeline.get_frame(time) - features = frame.features - size_data.append([feature.get_size() for feature in features.values()]) - mean_data.append([feature.mean for feature in features.values()]) - max_data.append([feature.max for feature in features.values()]) - number_of_features.append(len(features)) - - # Pad data with NaNs to create arrays of consistent shape (max number of features across - # all timesteps) - arr_size = max(number_of_features) - - # Size data is integer, but we need to pad with NaNs (which is a float), so fill - # with invalid value first - size_data = np.array( - [ - np.pad(sizes, (0, arr_size - len(sizes)), constant_values=-100) - for sizes in size_data - ], - dtype=float, + # Get feature data from each frame of data + size_data, mean_data, max_data = _get_cell_stats_arrays_from_timeline( + timeline=timeline, expected_frame_times=times_dt ) - size_data[size_data == -100] = np.nan # Get effective radius from feature size, using horizontal coordinate of input cube to estimate grid spacing effective_radius_data = _get_effective_radius_from_feature_size( size_data=size_data, cube_with_hzntl_coord=cube ) - # Mean and max data are already float, so can be padded with NaNs directly. - mean_data = np.array( - [ - np.pad(means, (0, arr_size - len(means)), constant_values=np.nan) - for means in mean_data - ] - ) - max_data = np.array( - [ - np.pad(maxs, (0, arr_size - len(maxs)), constant_values=np.nan) - for maxs in max_data - ] - ) - - # Set cube properties + # Set output cube properties cube_properties = { "feature_size": { "data": size_data, "long_name": "feature_size", - "units": "1", + "units": 1, }, "feature_mean": { "data": mean_data, @@ -344,62 +307,75 @@ def cell_stats( }, } - # Create cubes - time_coord = cube.coord("time").copy() - feature_coord = iris.coords.DimCoord( - np.arange(arr_size), - long_name="feature_number", - var_name="feature_number", - units="1", + # Create cubes, add to existing cubelist + cell_stats_cubelist.extend( + _add_cell_stats_data_to_cubes( + data_and_metadata_dict=cube_properties, template_cube=cube + ) ) - coords = [time_coord, feature_coord] - coords_and_dims = [(coord, i) for i, coord in enumerate(coords)] - - # Get list of coords to copy from input cube to output cubes - copyable_coord_names = [ - "realization", - "hour", - "forecast_period", - "forecast_reference_time", - "model_name", - "cset_comparison_base", + + return cell_stats_cubelist + + +def _get_cell_stats_arrays_from_timeline( + timeline: Timeline, expected_frame_times: list +) -> list[np.ndarray]: + """Extract cell statistics data from a Simple-Track Timeline object. + + Parameters + ---------- + timeline: Timeline + A Simple-Track Timeline object containing tracked features. + + Returns + ------- + size_data: np.ndarray + A numpy array containing the size of each feature in grid points. + mean_data: np.ndarray + A numpy array containing the mean value of each feature. + max_data: np.ndarray + A numpy array containing the maximum value of each feature. + """ + size_data, mean_data, max_data = [], [], [] + number_of_features = [] + for time in expected_frame_times: + frame = timeline.get_frame(time) + features = frame.features + size_data.append([feature.get_size() for feature in features.values()]) + mean_data.append([feature.mean for feature in features.values()]) + max_data.append([feature.max for feature in features.values()]) + number_of_features.append(len(features)) + + # Pad data with NaNs to create arrays of consistent shape (max number of features across + # all timesteps) + arr_size = max(number_of_features) + + # Size data is integer, but we need to pad with NaNs (which is a float), so fill + # with invalid value first + size_data = np.array( + [ + np.pad(sizes, (0, arr_size - len(sizes)), constant_values=-100) + for sizes in size_data + ], + dtype=float, + ) + size_data[size_data == -100] = np.nan + + # Mean and max data are already float, so can be padded with NaNs directly. + mean_data = np.array( + [ + np.pad(means, (0, arr_size - len(means)), constant_values=np.nan) + for means in mean_data ] - input_cube_coord_names = [] - for coord in cube.coords(): - input_cube_coord_names.append(coord.standard_name) - input_cube_coord_names.append(coord.long_name) - - coords_to_copy = [ - coord_name - for coord_name in copyable_coord_names - if coord_name in input_cube_coord_names + ) + max_data = np.array( + [ + np.pad(maxs, (0, arr_size - len(maxs)), constant_values=np.nan) + for maxs in max_data ] + ) - # Populate cubelist - for cb_props in cube_properties.values(): - cell_stats_cube = iris.cube.Cube( - data=cb_props["data"], - long_name=cb_props["long_name"], - units=cb_props["units"], - dim_coords_and_dims=coords_and_dims, - ) - # Add other metadata from input cube - for coord_name in coords_to_copy: - coord = cube.coord(coord_name).copy() - # Check if this coord represents a dimension of data - dims = cube.coord_dims(coord) - if len(dims) > 0: - cell_stats_cube.add_aux_coord(coord, dims) - else: - cell_stats_cube.add_aux_coord(coord) - - # Copy over attributes - cell_stats_cube.attributes = cube.attributes - - # Add to cubelist - cell_stats_cubelist.append(cell_stats_cube) - - return cell_stats_cubelist + return size_data, mean_data, max_data def _get_effective_radius_from_feature_size( @@ -447,3 +423,89 @@ def _get_effective_radius_from_feature_size( effective_radii_data = np.sqrt(size_data * grid_spacing**2 / np.pi) return effective_radii_data + + +def _add_cell_stats_data_to_cubes( + data_and_metadata_dict: dict, template_cube: iris.cube.Cube +) -> iris.cube.CubeList: + """Add data to cubes, using template cube for metadata. + + Parameters + ---------- + data_and_metadata_dict: dict + A dictionary containing data and metadata for each cube to be created. + The keys are the long names of the cubes, and the values are dictionaries + containing the data and units for each cube. + + template_cube: iris.cube.Cube + An iris cube to use as a template for the new cubes. The new cubes will + have the same attributes as the template cube. + + Returns + ------- + cubelist: iris.cube.CubeList + A list of iris cubes containing the added data. + + """ + cubelist = iris.cube.CubeList() + + # Construct coordinates for new cubes + time_coord = template_cube.coord("time").copy() + # To construct feature coordinate, look at the size of dimension 1 for each data + arr_size = max( + [data_and_metadata_dict[cb]["data"].shape[1] for cb in data_and_metadata_dict] + ) + feature_coord = iris.coords.DimCoord( + np.arange(arr_size), + long_name="feature_number", + var_name="feature_number", + units="1", + ) + coords = [time_coord, feature_coord] + coords_and_dims = [(coord, i) for i, coord in enumerate(coords)] + + # Get list of coords to copy from input cube to output cubes + copyable_coord_names = [ + "realization", + "hour", + "forecast_period", + "forecast_reference_time", + "model_name", + "cset_comparison_base", + ] + input_cube_coord_names = [] + for coord in template_cube.coords(): + input_cube_coord_names.append(coord.standard_name) + input_cube_coord_names.append(coord.long_name) + + coords_to_copy = [ + coord_name + for coord_name in copyable_coord_names + if coord_name in input_cube_coord_names + ] + + # Populate cubelist + for cb_props in data_and_metadata_dict.values(): + cell_stats_cube = iris.cube.Cube( + data=cb_props["data"], + long_name=cb_props["long_name"], + units=cb_props["units"], + dim_coords_and_dims=coords_and_dims, + ) + # Add other metadata from input cube + for coord_name in coords_to_copy: + coord = template_cube.coord(coord_name).copy() + # Check if this coord represents a dimension of data + dims = template_cube.coord_dims(coord) + if len(dims) > 0: + cell_stats_cube.add_aux_coord(coord, dims) + else: + cell_stats_cube.add_aux_coord(coord) + + # Copy over attributes + cell_stats_cube.attributes = template_cube.attributes + + # Add to cubelist + cubelist.append(cell_stats_cube) + + return cubelist diff --git a/tests/operators/test_feature.py b/tests/operators/test_feature.py index 4d2565a7ca..0ea13ebdff 100644 --- a/tests/operators/test_feature.py +++ b/tests/operators/test_feature.py @@ -30,9 +30,9 @@ def feature_cube() -> iris.cube.Cube: """Set up three timesteps of data and place into cube.""" data_arr = np.zeros((3, 10, 10)) - data_arr[0, 2:6, 2:6] = 1 - data_arr[1, 3:7, 3:7] = 1 - data_arr[2, 4:8, 4:8] = 1 + data_arr[0, 2:6, 2:6] = 5 + data_arr[1, 3:7, 3:7] = 10 + data_arr[2, 4:8, 4:8] = 20 time_units = cf_units.Unit("days since 2000-01-01 00:00:00", calendar="gregorian") time_start = dt.datetime(2010, 1, 1, 0, 0, 0) @@ -127,3 +127,38 @@ def test_save_data(feature_cube, tmp_path) -> None: # Check expected csv file is created in output directory expected_file = f"{output_directory}/frame_20100101_0000.csv" assert os.path.isfile(expected_file) + + +def test_cell_stats_operator(feature_cube): + """ + Test the cell_stats operator returns expected size, mean, and max values. + + The expected values are based on the feature_cube data and the threshold of 0.5. + """ + threshold = 0.5 + min_size = 1 + cubelist = feature.cell_stats( + cubes=feature_cube, threshold=threshold, min_size=min_size + ) + + # Extract data from cubelist, squeeze since there is only one feature per timestep + # in this test case + size_data = np.squeeze(cubelist.extract_cube("feature_size").data) + mean_data = np.squeeze(cubelist.extract_cube("feature_mean").data) + max_data = np.squeeze(cubelist.extract_cube("feature_max").data) + effective_radius_data = np.squeeze( + cubelist.extract_cube("feature_effective_radius").data + ) + + # Expected values based on the feature_cube data + expected_size_data = np.array([16, 16, 16]) # Each feature is a 4x4 square + expected_mean_data = np.array([5.0, 10.0, 20.0]) # Mean values of each feature + expected_max_data = np.array([5.0, 10.0, 20.0]) # Max values of each feature + + grid_spacing = 10 # Assuming grid spacing is 10 meters from test setup + expected_radius_data = np.sqrt(expected_size_data * grid_spacing**2 / np.pi) + + np.testing.assert_array_equal(size_data, expected_size_data) + np.testing.assert_array_equal(mean_data, expected_mean_data) + np.testing.assert_array_equal(max_data, expected_max_data) + np.testing.assert_array_equal(effective_radius_data, expected_radius_data) From 244a6e25f0364a4276e10bfc50163d5fdd0a5dc9 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Tue, 30 Jun 2026 15:11:44 +0100 Subject: [PATCH 21/33] added requirement that input cube is xy coord type --- src/CSET/operators/feature.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index 31a36e9b1c..1053f6abb4 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -202,6 +202,11 @@ def cell_stats( Parameters ---------- + cubes: iris.cube.Cube | iris.cube.CubeList + An iris cube (single model) or cubelist (multiple models) containing 2D data to be + analysed. Cube must have horizontal coordinates of xy type, not latitude/longitude. + The cube must also have a time coordinate, which is used to identify features in + each timestep. threshold: float The threshold value for feature detection. under_threshold: bool, optional @@ -243,6 +248,24 @@ def cell_stats( """ # Check inputs cubes = iter_maybe(cubes) + + # Require inputs to have horizontal coordinates of xy type, not latitude/longitude + for cube in cubes: + hzntl_coords = [ + coord + for coord in cube.coords() + if iris.util.guess_coord_axis(coord) in ["X", "Y"] + ] + invalid_coord_names = ["latitude", "longitude"] + for coord in hzntl_coords: + if coord.name() in invalid_coord_names: + raise ValueError( + f"Input cube {cube} has horizontal coordinate {coord}, " + "which is not of xy type. Please provide a cube with horizontal " + "coordinates of xy type." + ) + + # Setup containing cube list cell_stats_cubelist = iris.cube.CubeList() # Run tracking on all input data From e7f5fd6bcea47258dfec76d32e450ec38d4a3f1a Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Tue, 30 Jun 2026 15:19:39 +0100 Subject: [PATCH 22/33] made xy checking into function --- src/CSET/operators/feature.py | 42 ++++++++++++++++++++++++----------- 1 file changed, 29 insertions(+), 13 deletions(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index 1053f6abb4..9cf8503b1a 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -251,19 +251,7 @@ def cell_stats( # Require inputs to have horizontal coordinates of xy type, not latitude/longitude for cube in cubes: - hzntl_coords = [ - coord - for coord in cube.coords() - if iris.util.guess_coord_axis(coord) in ["X", "Y"] - ] - invalid_coord_names = ["latitude", "longitude"] - for coord in hzntl_coords: - if coord.name() in invalid_coord_names: - raise ValueError( - f"Input cube {cube} has horizontal coordinate {coord}, " - "which is not of xy type. Please provide a cube with horizontal " - "coordinates of xy type." - ) + _check_xy_coords(cube) # Setup containing cube list cell_stats_cubelist = iris.cube.CubeList() @@ -340,6 +328,34 @@ def cell_stats( return cell_stats_cubelist +def _check_xy_coords(cube: iris.cube.Cube) -> None: + """Check that the input cube has horizontal coordinates of xy type, not latitude/longitude. + + Parameters + ---------- + cube: iris.cube.Cube + An iris cube containing 2D data to be analysed. + + Raises + ------ + ValueError + If the input cube has horizontal coordinates of latitude/longitude type. + """ + hzntl_coords = [ + coord + for coord in cube.coords() + if iris.util.guess_coord_axis(coord) in ["X", "Y"] + ] + invalid_coord_names = ["latitude", "longitude", "grid_latitude", "grid_longitude"] + for coord in hzntl_coords: + if coord.name() in invalid_coord_names: + raise ValueError( + f"Input cube {cube} has horizontal coordinate {coord}, " + "which is not of xy type. Please provide a cube with horizontal " + "coordinates of xy type." + ) + + def _get_cell_stats_arrays_from_timeline( timeline: Timeline, expected_frame_times: list ) -> list[np.ndarray]: From f7b8a6a2fcc37604c946f656d8d4055ee495c05e Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Wed, 12 Aug 2026 11:34:51 +0100 Subject: [PATCH 23/33] nanmax when finding vmin, vmax --- src/CSET/operators/plot.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index 26ad0fd037..843b3c6378 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -443,8 +443,8 @@ def _set_axis_range(cubes): break if levels is None: - vmin = min(cb.data.min() for cb in cubes) - vmax = max(cb.data.max() for cb in cubes) + vmin = min(np.nanmin(cb.data) for cb in cubes) + vmax = max(np.nanmax(cb.data) for cb in cubes) return vmin, vmax From 9446b273de361a086c368777506f3f8ede7f6db2 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Fri, 14 Aug 2026 12:10:19 +0100 Subject: [PATCH 24/33] Adds misc.flatten operator, which can remove nans in cell stats cubes. Includes tests. Changes bins for different cell stats outputs. Adds condition in collapse.collapse for "forecast" coords before attempting to extract common times. Updates example recipe --- src/CSET/operators/collapse.py | 26 +++++--- src/CSET/operators/misc.py | 62 +++++++++++++++++++ src/CSET/operators/plot.py | 14 +++++ .../example_feature_cell_stats.yaml | 7 ++- tests/operators/test_misc.py | 52 ++++++++++++++++ 5 files changed, 150 insertions(+), 11 deletions(-) diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index dcd2994844..cae6eacd57 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -75,18 +75,24 @@ def collapse( raise ValueError("Must specify additional_percent") # Retain only common time points between different models if multiple model inputs. + # Do this only if "forecast_reference_time" and "forecast_period" are present in the cubes. if isinstance(cubes, iris.cube.CubeList) and len(cubes) > 1: - logging.debug( - "Extracting common time points as multiple model inputs detected." - ) - for cube in cubes: - cube.coord("forecast_reference_time").bounds = None - cube.coord("forecast_period").bounds = None - cubes = cubes.extract_overlapping( - ["forecast_reference_time", "forecast_period"] + fcst_ref_time_check = all( + "forecast_reference_time" in cube.coords() for cube in cubes ) - if len(cubes) == 0: - raise ValueError("No overlapping times detected in input cubes.") + fcst_period_check = all("forecast_period" in cube.coords() for cube in cubes) + if fcst_ref_time_check and fcst_period_check: + logging.debug( + "Extracting common time points as multiple model inputs detected." + ) + for cube in cubes: + cube.coord("forecast_reference_time").bounds = None + cube.coord("forecast_period").bounds = None + cubes = cubes.extract_overlapping( + ["forecast_reference_time", "forecast_period"] + ) + if len(cubes) == 0: + raise ValueError("No overlapping times detected in input cubes.") collapsed_cubes = iris.cube.CubeList([]) with warnings.catch_warnings(): diff --git a/src/CSET/operators/misc.py b/src/CSET/operators/misc.py index f667273eb8..b8a5b0d5f8 100644 --- a/src/CSET/operators/misc.py +++ b/src/CSET/operators/misc.py @@ -702,3 +702,65 @@ def differentiate( return new_cubelist[0] else: return new_cubelist + + +def flatten( + cubes: iris.cube.Cube | iris.cube.CubeList, remove_nans: bool = False +) -> iris.cube.Cube | iris.cube.CubeList: + """Flatten a cube or cubelist along all dimensions. + + Flattened cube contains a single dimension coordinate named "flattened_index". + + Parameters + ---------- + cubes : iris.cube.Cube or iris.cube.CubeList + The input Cube or CubeList to flatten. + remove_nans : bool, optional + If True, remove NaN values from the flattened data. Default is True. + + Returns + ------- + iris.cube.Cube or iris.cube.CubeList + The flattened cube or cubelist. + """ + if isinstance(cubes, iris.cube.Cube): + cubes = iris.cube.CubeList([cubes]) + + if not isinstance(cubes, iris.cube.CubeList): + raise TypeError("Input must be an iris.cube.Cube or iris.cube.CubeList.") + + flattened_cubes = iris.cube.CubeList() + for cube in cubes: + # Remove NaN if required + cube_data = cube.data + if remove_nans: + cube_data = cube_data[~np.isnan(cube_data)] + + # Flatten the data + flattened_data = cube_data.flatten() + + # Create a new cube with the flattened data and the remaining coordinates + flat_coord = iris.coords.DimCoord( + np.arange(flattened_data.size), long_name="flattened_index", units="1" + ) + flattened_cube = iris.cube.Cube( + flattened_data, + standard_name=cube.standard_name, + long_name=cube.long_name, + var_name=cube.var_name, + units=cube.units, + attributes=cube.attributes, + dim_coords_and_dims=[(flat_coord, 0)], + ) + + # Add single time point as a scalar coord if it exists + if cube.coords("time"): + time_coord = cube.coord("time") + flattened_cube.add_aux_coord(time_coord[0]) + + flattened_cubes.append(flattened_cube) + + if len(flattened_cubes) == 1: + return flattened_cubes[0] + else: + return flattened_cubes diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index 843b3c6378..3af7368ddd 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -1540,6 +1540,20 @@ def _plot_and_save_histogram_series( ax.set_xscale("log") elif "lightning" in title: bins = [0, 1, 2, 3, 4, 5] + elif "feature_size" in cube.long_name: + bins = np.linspace(0, 500, 51) + elif "feature_effective_radius" in cube.long_name: + # From RMED toolbox + bins = 10 ** (np.arange(0.0, 3.12, 0.12)) + bins = np.insert(bins, 0, 0) + vmin = bins[1] + vmax = bins[-1] + elif "feature_mean" in cube.long_name or "feature_max" in cube.long_name: + # From RMED toolbox + bins = 10 ** (np.arange(-1, 2.7, 0.12)) + bins = np.insert(bins, 0, 0) + vmin = bins[1] + vmax = bins[-1] else: bins = np.linspace(vmin, vmax, 51) logging.debug( diff --git a/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml b/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml index 9acf8e3811..4f0f46153c 100755 --- a/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml +++ b/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml @@ -2,7 +2,7 @@ category: Quick Look title: Example running cell_stats and plotting histograms description: | Uses the feature.cell_stats operator to calculate cell size, mean cell values and max cell values for - identified features. + identified features. Saves cell stats data to output and plots histograms steps: - operator: read.read_cubes @@ -16,6 +16,7 @@ steps: - operator: feature.cell_stats threshold: 3 + save_data: True # save raw tracking data for further analysis # Filter tracking cubelist to one of "feature_mean", "feature_max", or "feature_size" - operator: filters.filter_multiple_cubes @@ -23,4 +24,8 @@ steps: operator: constraints.generate_var_constraint varname: feature_size + # Flatten data across case study period + - operator: misc.flatten + remove_nans: True + - operator: plot.plot_histogram_series diff --git a/tests/operators/test_misc.py b/tests/operators/test_misc.py index 72fc8d1d57..e5e230f4bf 100644 --- a/tests/operators/test_misc.py +++ b/tests/operators/test_misc.py @@ -603,3 +603,55 @@ def test_not_remove_non_scalar_coord(): # Check it is still present cube_out = out[0] assert cube_out.coords("time") + + +def test_flatten_cube_no_nans(): + """Test misc.flatten without nans in single cube.""" + data_shape = (3, 4) + data = np.arange(12).reshape(data_shape) + coords = [ + iris.coords.DimCoord(np.arange(shape), long_name=f"test{shape}", units="1") + for shape in data_shape + ] + dim_coords_and_dims = [(coord, dim) for dim, coord in enumerate(coords)] + cube = iris.cube.Cube(data, dim_coords_and_dims=dim_coords_and_dims) + flattened_cube = misc.flatten(cube) + assert flattened_cube.shape == (12,) + assert np.allclose(flattened_cube.data, np.arange(12)) + + +def test_flatten_cubelist_no_nans(): + """Test misc.flatten without nans in CubeList.""" + data_shape = (3, 4) + data = np.arange(12).reshape(data_shape) + coords = [ + iris.coords.DimCoord(np.arange(shape), long_name=f"test{shape}", units="1") + for shape in data_shape + ] + dim_coords_and_dims = [(coord, dim) for dim, coord in enumerate(coords)] + cubes = iris.cube.CubeList( + [ + iris.cube.Cube(data, dim_coords_and_dims=dim_coords_and_dims) + for __ in range(3) + ] + ) + flattened_cubes = misc.flatten(cubes) + assert len(flattened_cubes) == 3 + for cube in flattened_cubes: + assert cube.shape == (12,) + assert np.allclose(cube.data, np.arange(12)) + + +def test_flatten_cube_nans_removed(): + """Test misc.flatten with nans removed in Cube.""" + data_shape = (3, 4) + data = np.arange(12, dtype=float).reshape(data_shape) + data[:, 0] = np.nan # 3 nans are inserted + coords = [ + iris.coords.DimCoord(np.arange(shape), long_name=f"test{shape}", units="1") + for shape in data_shape + ] + dim_coords_and_dims = [(coord, dim) for dim, coord in enumerate(coords)] + cube = iris.cube.Cube(data, dim_coords_and_dims=dim_coords_and_dims) + flattened_cube = misc.flatten(cube, remove_nans=True) + assert flattened_cube.shape == (9,) From a8e7498352ca40473afe733c511da06fe0760fb9 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Fri, 14 Aug 2026 13:11:34 +0100 Subject: [PATCH 25/33] Added grid_spacing as attribute to output cubes in cell_stats operator. Larger range of effective_radius bins in histogram operator --- src/CSET/operators/feature.py | 11 +++++++++-- src/CSET/operators/plot.py | 6 ++++-- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index 3c6b310e2f..307e8fc3aa 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -300,10 +300,14 @@ def cell_stats( ) # Get effective radius from feature size, using horizontal coordinate of input cube to estimate grid spacing - effective_radius_data = _get_effective_radius_from_feature_size( + effective_radius_data, grid_spacing = _get_effective_radius_from_feature_size( size_data=size_data, cube_with_hzntl_coord=cube ) + # Add grid_spacing as an attribute to the template_cube, so it is copied to + # output cubes in following function + cube.attributes["grid_spacing"] = grid_spacing + # Set output cube properties cube_properties = { "feature_size": { @@ -441,6 +445,9 @@ def _get_effective_radius_from_feature_size( effective_radii_data: np.ndarray An array containing "feature_effective_radius" data, in units of km. + grid_spacing: float + The estimated grid spacing in m, calculated from the horizontal coordinate of the input cube. + Notes ----- This function assumes that the input cube has a horizontal coordinate system that is regular and @@ -467,7 +474,7 @@ def _get_effective_radius_from_feature_size( grid_spacing = np.abs(np.mean(np.diff(hzntl_coord.points))) effective_radii_data = np.sqrt(size_data * grid_spacing**2 / np.pi) - return effective_radii_data + return effective_radii_data, grid_spacing def _add_cell_stats_data_to_cubes( diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index 3af7368ddd..ec0278ff11 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -1543,8 +1543,10 @@ def _plot_and_save_histogram_series( elif "feature_size" in cube.long_name: bins = np.linspace(0, 500, 51) elif "feature_effective_radius" in cube.long_name: - # From RMED toolbox - bins = 10 ** (np.arange(0.0, 3.12, 0.12)) + # TODO: use grid_spacing attribute in cubes to find min bin size + # for effective radius, rather than being hard coded + # Modified from RMED toolbox + bins = 10 ** (np.arange(0, 5.28, 0.12)) bins = np.insert(bins, 0, 0) vmin = bins[1] vmax = bins[-1] From 7158074d6aeb883e0624b644b32df3a5796b41e7 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Thu, 20 Aug 2026 14:43:53 +0100 Subject: [PATCH 26/33] added support for multiple thresholds (one per model) --- src/CSET/operators/feature.py | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index 307e8fc3aa..b7a0830a39 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -199,7 +199,7 @@ def track( def cell_stats( cubes: iris.cube.Cube | iris.cube.CubeList, - threshold: float, + threshold: float | list[float], under_threshold: bool = False, min_size: int = 4, save_data: bool = False, @@ -213,8 +213,10 @@ def cell_stats( analysed. Cube must have horizontal coordinates of xy type, not latitude/longitude. The cube must also have a time coordinate, which is used to identify features in each timestep. - threshold: float - The threshold value for feature detection. + threshold: float | list[float] + The threshold value(s) for feature detection. If a list is provided, each value + is used to identify features in the corresponding cube in the cubelist. Therefore, + the list should match the number of models. under_threshold: bool, optional If set to True, features are identified where the data is below the threshold. If set to False, features are identified where the data is above the threshold. @@ -262,13 +264,24 @@ def cell_stats( # Setup containing cube list cell_stats_cubelist = iris.cube.CubeList() + # If threshold is a list, check that it matches the number of cubes + if isinstance(threshold, list): + if len(threshold) != len(cubes): + raise ValueError( + f"Length of threshold list ({len(threshold)}) does not match " + f"number of cubes ({len(cubes)})." + ) + # else, make it iterable by repeating the same value for each cube + else: + threshold = [threshold] * len(cubes) + # Run tracking on all input data - for cube in cubes: + for cube, thresh in zip(cubes, threshold, strict=True): model_name = cube.attributes.get("model_name", None) # Setup config tracker_config = { "FEATURE": { - "threshold": threshold, + "threshold": thresh, "under_threshold": under_threshold, "min_size": min_size, }, From e3a300425894974455b5a9427f51307302cf24e3 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Thu, 27 Aug 2026 12:30:31 +0100 Subject: [PATCH 27/33] Relaxed requirement for input to be on xy grid, instead checks for uniform grid and converts effective radius to km. Updated tests --- src/CSET/operators/feature.py | 69 +++++++++++++++++++++++++-------- tests/operators/test_feature.py | 37 ++++++++++-------- 2 files changed, 73 insertions(+), 33 deletions(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index b7a0830a39..a8c1b1b954 100755 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -118,9 +118,6 @@ def track( >>> plt.show() """ - # Check that the input cube has horizontal coordinates of xy type, not latitude/longitude - _check_xy_coords(cube) - # Setup config tracker_config = { "FEATURE": { @@ -257,9 +254,9 @@ def cell_stats( # Check inputs cubes = iter_maybe(cubes) - # Require inputs to have horizontal coordinates of xy type, not latitude/longitude + # Require inputs to have a uniform grid for cube in cubes: - _check_xy_coords(cube) + _check_uniform_grid(cube) # Setup containing cube list cell_stats_cubelist = iris.cube.CubeList() @@ -351,32 +348,42 @@ def cell_stats( return cell_stats_cubelist -def _check_xy_coords(cube: iris.cube.Cube) -> None: - """Check that the input cube has horizontal coordinates of xy type, not latitude/longitude. +def _check_uniform_grid(cube: iris.cube.Cube) -> bool: + """Check that the input cube has approximately uniform horizontal grid spacing. + + Prints warning if cube does not have uniform grid. Parameters ---------- cube: iris.cube.Cube An iris cube containing horizontal coordinates. + Returns + ------- + bool + True if the cube has a uniform grid, False otherwise. + Raises ------ ValueError - If the input cube has horizontal coordinates of latitude/longitude type. + If the input cube does not have a uniform grid. """ hzntl_coords = [ coord for coord in cube.coords() if iris.util.guess_coord_axis(coord) in ["X", "Y"] ] - invalid_coord_names = ["latitude", "longitude", "grid_latitude", "grid_longitude"] + for coord in hzntl_coords: - if coord.name() in invalid_coord_names: - raise ValueError( - f"Input cube {cube} has horizontal coordinate {coord}, " - "which is not of xy type. Please provide a cube with horizontal " - "coordinates of xy type." + if not iris.util.is_regular(coord): + warning_msg = ( + f"Horizontal coordinate {coord} is not regular. " + "Feature statistics calculation may be inaccurate." ) + logging.warning(warning_msg) + print(warning_msg) + return False + return True def _get_cell_stats_arrays_from_timeline( @@ -484,9 +491,39 @@ def _get_effective_radius_from_feature_size( "Effective radius calculation may be inaccurate." ) - grid_spacing = np.abs(np.mean(np.diff(hzntl_coord.points))) - effective_radii_data = np.sqrt(size_data * grid_spacing**2 / np.pi) + # Get grid spacing in native coord units (degrees, m, km etc) + grid_spacing = iris.util.regular_step(hzntl_coord) + + if hzntl_coord.units == "m": + grid_spacing = grid_spacing / 1000 # Convert to km + # If grid spacing is in degrees, convert to km using approximate conversion factor + if hzntl_coord.units == "degrees": + # Get latitude for better conversion to km + lat_coords = [ + coord + for coord in cube_with_hzntl_coord.coords() + if iris.util.guess_coord_axis(coord) in ["Y"] + ] + for coord in lat_coords: + if coord.units == "degrees": + lat_coord_for_conversion = coord + else: + lat_coord_for_conversion = None + + # Calculate conversion factor using latitude correction if available, + # otherwise use naive 111 km per degree conversion + if lat_coord_for_conversion is not None: + mean_latitude = np.mean(lat_coord_for_conversion.points) + else: + logging.warning( + "No latitude coordinate found for conversion to km. " + "Using naive conversion factor of 111 km per degree." + ) + mean_latitude = 0 + grid_spacing = grid_spacing * 111 * np.cos(np.radians(mean_latitude)) + + effective_radii_data = np.sqrt(size_data * grid_spacing**2 / np.pi) return effective_radii_data, grid_spacing diff --git a/tests/operators/test_feature.py b/tests/operators/test_feature.py index 1fe1c81797..60e5a7880b 100644 --- a/tests/operators/test_feature.py +++ b/tests/operators/test_feature.py @@ -14,6 +14,7 @@ """Tests for feature operators.""" import datetime as dt +import os import cf_units import iris @@ -118,7 +119,7 @@ def test_save_data(feature_cube, tmp_working_dir) -> None: save_data=True, ) # Check expected lifetime field is created in output directory - output_directory = f"{tmp_path}/tracking_data" + output_directory = f"{tmp_working_dir}/tracking_data" expected_file = f"{output_directory}/lifetime_20100101_0000.field" assert os.path.isfile(expected_file) @@ -127,7 +128,7 @@ def test_save_data(feature_cube, tmp_working_dir) -> None: assert os.path.isfile(expected_file) -def test_cell_stats_operator(feature_cube): +def test_cell_stats_operator(feature_cube, tmp_working_dir): """ Test the cell_stats operator returns expected size, mean, and max values. @@ -136,7 +137,7 @@ def test_cell_stats_operator(feature_cube): threshold = 0.5 min_size = 1 cubelist = feature.cell_stats( - cubes=feature_cube, threshold=threshold, min_size=min_size + cubes=feature_cube, threshold=threshold, min_size=min_size, save_data=True ) # Extract data from cubelist, squeeze since there is only one feature per timestep @@ -155,12 +156,14 @@ def test_cell_stats_operator(feature_cube): grid_spacing = 10 # Assuming grid spacing is 10 meters from test setup expected_radius_data = np.sqrt(expected_size_data * grid_spacing**2 / np.pi) + # Convert to km + expected_radius_data = expected_radius_data / 1000 np.testing.assert_array_equal(size_data, expected_size_data) np.testing.assert_array_equal(mean_data, expected_mean_data) np.testing.assert_array_equal(max_data, expected_max_data) - np.testing.assert_array_equal(effective_radius_data, expected_radius_data) - output_directory = tmp_working_dir / "tracking_data" + np.testing.assert_array_almost_equal(effective_radius_data, expected_radius_data) + output_directory = tmp_working_dir / "None/cell-stats_data" expected_file = output_directory / "lifetime_20100101_0000.field" assert expected_file.is_file() @@ -169,21 +172,21 @@ def test_cell_stats_operator(feature_cube): assert expected_file.is_file() -def test_check_xy_coords_valid(feature_cube) -> None: - """Test that _check_xy_coords does not raise an error for valid xy coordinates.""" - try: - feature._check_xy_coords(feature_cube) - except ValueError: - pytest.fail("Unexpected ValueError raised for valid xy coordinates.") +def test_check_uniform_grid(feature_cube) -> None: + """Test that _check_uniform_grid does not raise an error for valid uniform grid.""" + result = feature._check_uniform_grid(feature_cube) + assert result is True -def test_check_xy_coords_invalid() -> None: - """Test that _check_xy_coords raises a ValueError for invalid latitude/longitude coordinates.""" - # Create a cube with latitude and longitude coordinates +def test_check_uniform_grid_invalid() -> None: + """Test that _check_uniform_grid raises a ValueError for invalid uniform grid.""" + # Create a cube with non-uniform grid data_arr = np.zeros((10, 10)) + lat_points = np.linspace(-90, 90, 10) + lat_points[-1] = 95 lat_coord = iris.coords.DimCoord( - points=np.linspace(-90, 90, 10), + points=lat_points, standard_name="latitude", var_name="latitude", units="degrees", @@ -203,5 +206,5 @@ def test_check_xy_coords_invalid() -> None: long_name="Precipitation test", ) - with pytest.raises(ValueError): - feature._check_xy_coords(cube) + result = feature._check_uniform_grid(cube) + assert result is False From f4e0d3e9bae4369b9547a815a4f2acf8ef4a1803 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Thu, 27 Aug 2026 13:12:19 +0100 Subject: [PATCH 28/33] fixed failing tests --- src/CSET/operators/collapse.py | 23 ++++++++--------------- src/CSET/operators/feature.py | 13 ++++++++----- 2 files changed, 16 insertions(+), 20 deletions(-) diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index 817dc4cec0..bea4e82dff 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -79,22 +79,15 @@ def collapse( # Retain only common time points between different models if multiple model inputs. # Do this only if "forecast_reference_time" and "forecast_period" are present in the cubes. if isinstance(cubes, iris.cube.CubeList) and len(cubes) > 1: - fcst_ref_time_check = all( - "forecast_reference_time" in cube.coords() for cube in cubes + logger.debug("Extracting common time points as multiple model inputs detected.") + for cube in cubes: + cube.coord("forecast_reference_time").bounds = None + cube.coord("forecast_period").bounds = None + cubes = cubes.extract_overlapping( + ["forecast_reference_time", "forecast_period"] ) - fcst_period_check = all("forecast_period" in cube.coords() for cube in cubes) - if fcst_ref_time_check and fcst_period_check: - logger.debug( - "Extracting common time points as multiple model inputs detected." - ) - for cube in cubes: - cube.coord("forecast_reference_time").bounds = None - cube.coord("forecast_period").bounds = None - cubes = cubes.extract_overlapping( - ["forecast_reference_time", "forecast_period"] - ) - if len(cubes) == 0: - raise ValueError("No overlapping times detected in input cubes.") + if len(cubes) == 0: + raise ValueError("No overlapping times detected in input cubes.") collapsed_cubes = iris.cube.CubeList([]) with warnings.catch_warnings(): diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index 8c0018b9f7..5046eb7ec2 100644 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -482,12 +482,15 @@ def _get_effective_radius_from_feature_size( """ # Guess coord representing horizontal grid (choose first available) hzntl_coord = next( - [ - coord - for coord in cube_with_hzntl_coord.coords() - if iris.util.guess_coord_axis(coord) in ["X", "Y"] - ] + iter( + [ + coord + for coord in cube_with_hzntl_coord.coords() + if iris.util.guess_coord_axis(coord) in ["X", "Y"] + ] + ) ) + logger.debug(f"Attempting to convert to effective radius using {hzntl_coord}") # Check coordinate is regular, but only warn if not, this is a naive estimate From 16bec607a34d7eb3455a8941b2c428337618d287 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Thu, 27 Aug 2026 14:02:22 +0100 Subject: [PATCH 29/33] added more feature tests --- src/CSET/operators/feature.py | 3 + tests/operators/test_feature.py | 221 +++++++++++++++++++++++++++++++- 2 files changed, 223 insertions(+), 1 deletion(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index 5046eb7ec2..0756c1f9ab 100644 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -401,6 +401,9 @@ def _get_cell_stats_arrays_from_timeline( timeline: Timeline A Simple-Track Timeline object containing tracked features. + expected_frame_times: list + A list of expected frame times to extract data for. + Returns ------- size_data: np.ndarray diff --git a/tests/operators/test_feature.py b/tests/operators/test_feature.py index 60e5a7880b..fb188e2ec5 100644 --- a/tests/operators/test_feature.py +++ b/tests/operators/test_feature.py @@ -22,11 +22,13 @@ import iris.cube import numpy as np import pytest +from simpletrack import Tracker +from simpletrack.frame import Timeline from CSET.operators import feature -@pytest.fixture +@pytest.fixture(scope="session") def feature_cube() -> iris.cube.Cube: """Set up three timesteps of data and place into cube.""" data_arr = np.zeros((3, 10, 10)) @@ -41,6 +43,9 @@ def feature_cube() -> iris.cube.Cube: time_coord = iris.coords.DimCoord( points=time_points, standard_name="time", units=time_units ) + fcst_ref_coord = iris.coords.AuxCoord( + points=time_start, standard_name="forecast_reference_time", units=time_units + ) coord_system = iris.coord_systems.TransverseMercator( latitude_of_projection_origin=55, longitude_of_central_meridian=0 @@ -66,14 +71,51 @@ def feature_cube() -> iris.cube.Cube: coords = (time_coord, proj_y_coord, proj_x_coord) dim_coords_and_dims = [(coord, dim) for dim, coord in enumerate(coords)] + attributes = {"title": "Precipitation test"} cube = iris.cube.Cube( data=data_arr, dim_coords_and_dims=dim_coords_and_dims, long_name="Precipitation test", + attributes=attributes, + aux_coords_and_dims=[(fcst_ref_coord, None)], ) return cube +@pytest.fixture(scope="session") +def cell_stats_timeline(feature_cube) -> Timeline: + """Run feature tracking on feature_cube and return Timeline object.""" + # Setup config + tracker_config = { + "FEATURE": { + "threshold": 1, + "under_threshold": False, + "min_size": 4, + }, + "OUTPUT": { + "save_data": False, + "experiment_name": "feature_tracking", + "path": f"{os.getcwd()}/model/cell-stats_data", + "skip_tracking": True, + }, + } + + # Get cube data into a dict to pass to Tracker + times = feature_cube.coord("time").points + time_units = feature_cube.coord("time").units + times_dt = [time_units.num2pydate(t) for t in times] + cube_dict = { + time: cube_slice.data + for time, cube_slice in zip( + times_dt, feature_cube.slices_over("time"), strict=True + ) + } + + # Run tracking, returning Timeline object + timeline = Tracker(tracker_config).run(cube_dict) + return timeline + + def test_tracking_valid(feature_cube) -> None: """ Test feature tracking returns same cube shape as input cube. @@ -172,6 +214,13 @@ def test_cell_stats_operator(feature_cube, tmp_working_dir): assert expected_file.is_file() +def test_cell_stats_invalid_threshold_list(feature_cube): + """Test that cell_stats raises a ValueError if passed a list of wrong size.""" + invalid_thresholds = [0.5, 1.0] + with pytest.raises(ValueError): + feature.cell_stats(cubes=feature_cube, threshold=invalid_thresholds, min_size=1) + + def test_check_uniform_grid(feature_cube) -> None: """Test that _check_uniform_grid does not raise an error for valid uniform grid.""" result = feature._check_uniform_grid(feature_cube) @@ -208,3 +257,173 @@ def test_check_uniform_grid_invalid() -> None: result = feature._check_uniform_grid(cube) assert result is False + + +def test_get_cell_stats_arrays_from_timeline(cell_stats_timeline): + """Test that _get_cell_stats_arrays_from_timeline returns expected arrays.""" + expected_frame_times = [ + dt.datetime(2010, 1, 1, 0, 0, 0), + dt.datetime(2010, 1, 1, 0, 5, 0), + dt.datetime(2010, 1, 1, 0, 10, 0), + ] + + size_array, mean_array, max_array = feature._get_cell_stats_arrays_from_timeline( + cell_stats_timeline, expected_frame_times + ) + + # Expected values based on the feature_cube data and the threshold of 0.5 + expected_size_array = np.array([[16], [16], [16]]) # Each feature is a 4x4 square + expected_mean_array = np.array([[5], [10], [20]]) # Mean values of each feature + expected_max_array = np.array([[5], [10], [20]]) # Max values of each feature + + np.testing.assert_array_equal(size_array, expected_size_array) + np.testing.assert_array_equal(mean_array, expected_mean_array) + np.testing.assert_array_equal(max_array, expected_max_array) + + +def test_get_effective_radius_from_feature_size(feature_cube): + """Test that _get_effective_radius_from_feature_size returns expected values.""" + # Use the same size data from above test + size_data = np.array([[16], [16], [16]]) # Each feature is a 4x4 square + + effective_radius_data, grid_spacing = ( + feature._get_effective_radius_from_feature_size(size_data, feature_cube) + ) + + # Expected values based on the feature_cube data + grid_spacing = 10 # Assuming grid spacing is 10 meters from test setup + expected_radius_data = np.sqrt(size_data * grid_spacing**2 / np.pi) + # Convert to km + expected_radius_data = expected_radius_data / 1000 + + np.testing.assert_array_almost_equal(effective_radius_data, expected_radius_data) + assert grid_spacing == 10 + + +def test_get_effective_radius_from_feature_size_km_input_cube(feature_cube): + """Test that _get_effective_radius_from_feature_size returns expected values for km input cube.""" + # Use the same size data from above test + size_data = np.array([[16], [16], [16]]) # Each feature is a 4x4 square + + # Convert cube to km units + feature_cube_km = feature_cube.copy() + feature_cube_km.coord("projection_x_coordinate").convert_units("km") + feature_cube_km.coord("projection_y_coordinate").convert_units("km") + + effective_radius_data, grid_spacing = ( + feature._get_effective_radius_from_feature_size(size_data, feature_cube_km) + ) + + # Expected values based on the feature_cube data + grid_spacing = 0.01 # Assuming grid spacing is 10 meters (0.01 km) from test setup + expected_radius_data = np.sqrt(size_data * grid_spacing**2 / np.pi) + + np.testing.assert_array_almost_equal(effective_radius_data, expected_radius_data) + assert grid_spacing == 0.01 + + +def test_get_effective_radius_from_feature_size_latlon_cube(): + """Test that _get_effective_radius_from_feature_size returns expected values for lat/lon cube.""" + # Create a cube with lat/lon coordinates + data_arr = np.zeros((10, 10)) + lat_points = np.linspace(-1, 1, 10) + lon_points = np.linspace(-1, 1, 10) + + lat_coord = iris.coords.DimCoord( + points=lat_points, + standard_name="latitude", + var_name="latitude", + units="degrees", + ) + lon_coord = iris.coords.DimCoord( + points=lon_points, + standard_name="longitude", + var_name="longitude", + units="degrees", + ) + + coords = (lat_coord, lon_coord) + dim_coords_and_dims = [(coord, dim) for dim, coord in enumerate(coords)] + cube = iris.cube.Cube( + data=data_arr, + dim_coords_and_dims=dim_coords_and_dims, + long_name="Precipitation test", + ) + + # Use the same size data from above test + size_data = np.array([[2], [2], [2]]) # Each feature is a 4x4 square + + effective_radius_data, grid_spacing = ( + feature._get_effective_radius_from_feature_size(size_data, cube) + ) + + # Expected values based on the feature_cube data + # For lat/lon cube, grid spacing is calculated based on the distance between points + # Expected grid spacing is 2 degrees (222 km) / 9 intervals = 24.666 km + expected_grid_spacing = 24.666 + expected_radius_data = np.sqrt(size_data * expected_grid_spacing**2 / np.pi) + + np.testing.assert_array_almost_equal( + effective_radius_data, expected_radius_data, decimal=3 + ) + np.testing.assert_almost_equal(grid_spacing, expected_grid_spacing, decimal=3) + + +def test_add_cell_stats_data_to_cubes(cell_stats_timeline, feature_cube): + """Test that _add_cell_stats_data_to_cubes adds expected data to cubelist.""" + expected_frame_times = [ + dt.datetime(2010, 1, 1, 0, 0, 0), + dt.datetime(2010, 1, 1, 0, 5, 0), + dt.datetime(2010, 1, 1, 0, 10, 0), + ] + + # Get feature data from each frame of data + size_data, mean_data, max_data = feature._get_cell_stats_arrays_from_timeline( + timeline=cell_stats_timeline, expected_frame_times=expected_frame_times + ) + + # Get effective radius from feature size, using horizontal coordinate of input cube to estimate grid spacing + effective_radius_data, __ = feature._get_effective_radius_from_feature_size( + size_data=size_data, cube_with_hzntl_coord=feature_cube + ) + + cube_properties = { + "feature_size": { + "data": size_data, + "long_name": "feature_size", + "units": 1, + }, + "feature_mean": { + "data": mean_data, + "long_name": "feature_mean", + "units": 1, + }, + "feature_max": {"data": max_data, "long_name": "feature_max", "units": 1}, + "feature_effective_radius": { + "data": effective_radius_data, + "long_name": "feature_effective_radius", + "units": "km", + }, + } + + cubelist = feature._add_cell_stats_data_to_cubes(cube_properties, feature_cube) + + # Test a cube of each name is produced + expected_cube_names = [ + "feature_size", + "feature_mean", + "feature_max", + "feature_effective_radius", + ] + + for cube_name in expected_cube_names: + assert any(cube.long_name == cube_name for cube in cubelist) + + # Check each cube contains a forecast_reference_time coordinate, copied from + # feature_cube + for cube in cubelist: + assert cube.coords("forecast_reference_time") + + # test cube attributes copied to cell stats cubes + for cube in cubelist: + assert cube.attributes == feature_cube.attributes From 40a245261cfd05adddaf0ad0dd273a151216a2ef Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Thu, 3 Sep 2026 16:35:52 +0100 Subject: [PATCH 30/33] fixed cmap undefined error for cell_stats outputs --- src/CSET/operators/_colormaps.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/CSET/operators/_colormaps.py b/src/CSET/operators/_colormaps.py index 9d37fc5db1..871c32c9cf 100644 --- a/src/CSET/operators/_colormaps.py +++ b/src/CSET/operators/_colormaps.py @@ -723,7 +723,11 @@ def custom_colormap_feature_tracking(cube: iris.cube.Cube): norm = mcolors.BoundaryNorm(levels, cmap.N) logger.info("change colormap for feature init variable colorbar.") - # Set all non-feature data to white. - cmap = cmap.with_extremes(under="white") + else: + cmap, levels, norm = None, None, None + + # Set all non-feature data to white + if cmap: + cmap = cmap.with_extremes(under="white") return cmap, levels, norm From f3f3a2276401ab75502f78edff1339724ff95ec4 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Thu, 3 Sep 2026 16:41:13 +0100 Subject: [PATCH 31/33] fixed inaccurate docstrings --- src/CSET/operators/feature.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index 0756c1f9ab..23cae2c159 100644 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -45,7 +45,6 @@ def track( ---------- cube: iris.cube.Cube The cube to identify features in. The cube must be 3D and contain a time coordinate - and horizontal coordinates of xy type (not latitude/longitude). threshold: float The threshold value for feature detection. under_threshold: bool, optional @@ -212,7 +211,7 @@ def cell_stats( ---------- cubes: iris.cube.Cube | iris.cube.CubeList An iris cube (single model) or cubelist (multiple models) containing 2D data to be - analysed. Cube must have horizontal coordinates of xy type, not latitude/longitude. + analysed. Cube must have horizontal coordinates on a regular grid. The cube must also have a time coordinate, which is used to identify features in each timestep. threshold: float | list[float] @@ -233,7 +232,7 @@ def cell_stats( Returns ------- - cell_stats_cubes: iris.cube.CubeList + cell_stats_cubelist: iris.cube.CubeList An iris CubeList containing "feature_size", "feature_effective_radius", "feature_mean", and "feature_max" cubes. From b878727e377a5233caf87aa4174d4030af0f626b Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Wed, 16 Sep 2026 11:25:17 +0100 Subject: [PATCH 32/33] changed effective_radius to effective_diameter --- src/CSET/operators/feature.py | 38 ++++++----- .../example_feature_cell_stats.yaml | 5 +- tests/operators/test_feature.py | 66 ++++++++++--------- 3 files changed, 60 insertions(+), 49 deletions(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index 23cae2c159..394b173a64 100644 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -233,14 +233,14 @@ def cell_stats( Returns ------- cell_stats_cubelist: iris.cube.CubeList - An iris CubeList containing "feature_size", "feature_effective_radius", "feature_mean", + An iris CubeList containing "feature_size", "feature_effective_diameter", "feature_mean", and "feature_max" cubes. Notes ----- This operator uses the Simple-Track package with tracking disabled to identify features in each timestep and compile cell statistics. Outputs cubes containing feature size (number - of grid points), effective radius (in km), mean value within features, and maximum + of grid points), effective diameter (in km), mean value within features, and maximum value within features. Links @@ -313,9 +313,11 @@ def cell_stats( timeline=timeline, expected_frame_times=times_dt ) - # Get effective radius from feature size, using horizontal coordinate of input cube to estimate grid spacing - effective_radius_data, grid_spacing = _get_effective_radius_from_feature_size( - size_data=size_data, cube_with_hzntl_coord=cube + # Get effective diameter from feature size, using horizontal coordinate of input cube to estimate grid spacing + effective_diameter_data, grid_spacing = ( + _get_effective_diameter_from_feature_size( + size_data=size_data, cube_with_hzntl_coord=cube + ) ) # Add grid_spacing as an attribute to the template_cube, so it is copied to @@ -335,9 +337,9 @@ def cell_stats( "units": 1, }, "feature_max": {"data": max_data, "long_name": "feature_max", "units": 1}, - "feature_effective_radius": { - "data": effective_radius_data, - "long_name": "feature_effective_radius", + "feature_effective_diameter": { + "data": effective_diameter_data, + "long_name": "feature_effective_diameter", "units": "km", }, } @@ -454,10 +456,10 @@ def _get_cell_stats_arrays_from_timeline( return size_data, mean_data, max_data -def _get_effective_radius_from_feature_size( +def _get_effective_diameter_from_feature_size( size_data: np.ndarray, cube_with_hzntl_coord: iris.cube.Cube ) -> np.ndarray: - """Convert feature size in grid points to effective radius in km. + """Convert feature size in grid points to effective diameter in km. Parameters ---------- @@ -469,8 +471,8 @@ def _get_effective_radius_from_feature_size( Returns ------- - effective_radii_data: np.ndarray - An array containing "feature_effective_radius" data, in units of km. + effective_diameters_data: np.ndarray + An array containing "feature_effective_diameter" data, in units of km. grid_spacing: float The estimated grid spacing in m, calculated from the horizontal coordinate of the input cube. @@ -478,8 +480,8 @@ def _get_effective_radius_from_feature_size( Notes ----- This function assumes that the input cube has a horizontal coordinate system that is regular and - that the grid spacing can be estimated from the horizontal coordinates. The effective radius is - calculated as the radius of a circle with the same area as the feature size in grid points. + that the grid spacing can be estimated from the horizontal coordinates. The effective diameter is + calculated as the diameter of a circle with the same area as the feature size in grid points. """ # Guess coord representing horizontal grid (choose first available) @@ -493,14 +495,14 @@ def _get_effective_radius_from_feature_size( ) ) - logger.debug(f"Attempting to convert to effective radius using {hzntl_coord}") + logger.debug(f"Attempting to convert to effective diameter using {hzntl_coord}") # Check coordinate is regular, but only warn if not, this is a naive estimate # and will be inaccurate for irregular grids if not iris.util.is_regular(hzntl_coord): logger.warning( f"Horizontal coordinate {hzntl_coord} is not regular. " - "Effective radius calculation may be inaccurate." + "Effective diameter calculation may be inaccurate." ) # Get grid spacing in native coord units (degrees, m, km etc) @@ -535,8 +537,8 @@ def _get_effective_radius_from_feature_size( mean_latitude = 0 grid_spacing = grid_spacing * 111 * np.cos(np.radians(mean_latitude)) - effective_radii_data = np.sqrt(size_data * grid_spacing**2 / np.pi) - return effective_radii_data, grid_spacing + effective_diameters_data = np.sqrt(size_data * grid_spacing**2 / np.pi) * 2 + return effective_diameters_data, grid_spacing def _add_cell_stats_data_to_cubes( diff --git a/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml b/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml index 4f0f46153c..1d7908ba85 100755 --- a/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml +++ b/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml @@ -18,11 +18,14 @@ steps: threshold: 3 save_data: True # save raw tracking data for further analysis - # Filter tracking cubelist to one of "feature_mean", "feature_max", or "feature_size" + # Filter tracking cubelist to one of "feature_mean", "feature_max", "feature_size", or "feature_effective_diameter" - operator: filters.filter_multiple_cubes constraint: operator: constraints.generate_var_constraint varname: feature_size + # varname: feature_mean + # varname: feature_max + # varname: feature_effective_diameter # Flatten data across case study period - operator: misc.flatten diff --git a/tests/operators/test_feature.py b/tests/operators/test_feature.py index fb188e2ec5..a125038fdf 100644 --- a/tests/operators/test_feature.py +++ b/tests/operators/test_feature.py @@ -187,8 +187,8 @@ def test_cell_stats_operator(feature_cube, tmp_working_dir): size_data = np.squeeze(cubelist.extract_cube("feature_size").data) mean_data = np.squeeze(cubelist.extract_cube("feature_mean").data) max_data = np.squeeze(cubelist.extract_cube("feature_max").data) - effective_radius_data = np.squeeze( - cubelist.extract_cube("feature_effective_radius").data + effective_diameter_data = np.squeeze( + cubelist.extract_cube("feature_effective_diameter").data ) # Expected values based on the feature_cube data @@ -197,14 +197,16 @@ def test_cell_stats_operator(feature_cube, tmp_working_dir): expected_max_data = np.array([5.0, 10.0, 20.0]) # Max values of each feature grid_spacing = 10 # Assuming grid spacing is 10 meters from test setup - expected_radius_data = np.sqrt(expected_size_data * grid_spacing**2 / np.pi) + expected_diameter_data = np.sqrt(expected_size_data * grid_spacing**2 / np.pi) * 2 # Convert to km - expected_radius_data = expected_radius_data / 1000 + expected_diameter_data = expected_diameter_data / 1000 np.testing.assert_array_equal(size_data, expected_size_data) np.testing.assert_array_equal(mean_data, expected_mean_data) np.testing.assert_array_equal(max_data, expected_max_data) - np.testing.assert_array_almost_equal(effective_radius_data, expected_radius_data) + np.testing.assert_array_almost_equal( + effective_diameter_data, expected_diameter_data + ) output_directory = tmp_working_dir / "None/cell-stats_data" expected_file = output_directory / "lifetime_20100101_0000.field" assert expected_file.is_file() @@ -281,27 +283,29 @@ def test_get_cell_stats_arrays_from_timeline(cell_stats_timeline): np.testing.assert_array_equal(max_array, expected_max_array) -def test_get_effective_radius_from_feature_size(feature_cube): - """Test that _get_effective_radius_from_feature_size returns expected values.""" +def test_get_effective_diameter_from_feature_size(feature_cube): + """Test that _get_effective_diameter_from_feature_size returns expected values.""" # Use the same size data from above test size_data = np.array([[16], [16], [16]]) # Each feature is a 4x4 square - effective_radius_data, grid_spacing = ( - feature._get_effective_radius_from_feature_size(size_data, feature_cube) + effective_diameter_data, grid_spacing = ( + feature._get_effective_diameter_from_feature_size(size_data, feature_cube) ) # Expected values based on the feature_cube data grid_spacing = 10 # Assuming grid spacing is 10 meters from test setup - expected_radius_data = np.sqrt(size_data * grid_spacing**2 / np.pi) + expected_diameter_data = np.sqrt(size_data * grid_spacing**2 / np.pi) * 2 # Convert to km - expected_radius_data = expected_radius_data / 1000 + expected_diameter_data = expected_diameter_data / 1000 - np.testing.assert_array_almost_equal(effective_radius_data, expected_radius_data) + np.testing.assert_array_almost_equal( + effective_diameter_data, expected_diameter_data + ) assert grid_spacing == 10 -def test_get_effective_radius_from_feature_size_km_input_cube(feature_cube): - """Test that _get_effective_radius_from_feature_size returns expected values for km input cube.""" +def test_get_effective_diameter_from_feature_size_km_input_cube(feature_cube): + """Test that _get_effective_diameter_from_feature_size returns expected values for km input cube.""" # Use the same size data from above test size_data = np.array([[16], [16], [16]]) # Each feature is a 4x4 square @@ -310,20 +314,22 @@ def test_get_effective_radius_from_feature_size_km_input_cube(feature_cube): feature_cube_km.coord("projection_x_coordinate").convert_units("km") feature_cube_km.coord("projection_y_coordinate").convert_units("km") - effective_radius_data, grid_spacing = ( - feature._get_effective_radius_from_feature_size(size_data, feature_cube_km) + effective_diameter_data, grid_spacing = ( + feature._get_effective_diameter_from_feature_size(size_data, feature_cube_km) ) # Expected values based on the feature_cube data grid_spacing = 0.01 # Assuming grid spacing is 10 meters (0.01 km) from test setup - expected_radius_data = np.sqrt(size_data * grid_spacing**2 / np.pi) + expected_diameter_data = np.sqrt(size_data * grid_spacing**2 / np.pi) * 2 - np.testing.assert_array_almost_equal(effective_radius_data, expected_radius_data) + np.testing.assert_array_almost_equal( + effective_diameter_data, expected_diameter_data + ) assert grid_spacing == 0.01 -def test_get_effective_radius_from_feature_size_latlon_cube(): - """Test that _get_effective_radius_from_feature_size returns expected values for lat/lon cube.""" +def test_get_effective_diameter_from_feature_size_latlon_cube(): + """Test that _get_effective_diameter_from_feature_size returns expected values for lat/lon cube.""" # Create a cube with lat/lon coordinates data_arr = np.zeros((10, 10)) lat_points = np.linspace(-1, 1, 10) @@ -353,18 +359,18 @@ def test_get_effective_radius_from_feature_size_latlon_cube(): # Use the same size data from above test size_data = np.array([[2], [2], [2]]) # Each feature is a 4x4 square - effective_radius_data, grid_spacing = ( - feature._get_effective_radius_from_feature_size(size_data, cube) + effective_diameter_data, grid_spacing = ( + feature._get_effective_diameter_from_feature_size(size_data, cube) ) # Expected values based on the feature_cube data # For lat/lon cube, grid spacing is calculated based on the distance between points # Expected grid spacing is 2 degrees (222 km) / 9 intervals = 24.666 km expected_grid_spacing = 24.666 - expected_radius_data = np.sqrt(size_data * expected_grid_spacing**2 / np.pi) + expected_diameter_data = np.sqrt(size_data * expected_grid_spacing**2 / np.pi) * 2 np.testing.assert_array_almost_equal( - effective_radius_data, expected_radius_data, decimal=3 + effective_diameter_data, expected_diameter_data, decimal=3 ) np.testing.assert_almost_equal(grid_spacing, expected_grid_spacing, decimal=3) @@ -382,8 +388,8 @@ def test_add_cell_stats_data_to_cubes(cell_stats_timeline, feature_cube): timeline=cell_stats_timeline, expected_frame_times=expected_frame_times ) - # Get effective radius from feature size, using horizontal coordinate of input cube to estimate grid spacing - effective_radius_data, __ = feature._get_effective_radius_from_feature_size( + # Get effective diameter from feature size, using horizontal coordinate of input cube to estimate grid spacing + effective_diameter_data, __ = feature._get_effective_diameter_from_feature_size( size_data=size_data, cube_with_hzntl_coord=feature_cube ) @@ -399,9 +405,9 @@ def test_add_cell_stats_data_to_cubes(cell_stats_timeline, feature_cube): "units": 1, }, "feature_max": {"data": max_data, "long_name": "feature_max", "units": 1}, - "feature_effective_radius": { - "data": effective_radius_data, - "long_name": "feature_effective_radius", + "feature_effective_diameter": { + "data": effective_diameter_data, + "long_name": "feature_effective_diameter", "units": "km", }, } @@ -413,7 +419,7 @@ def test_add_cell_stats_data_to_cubes(cell_stats_timeline, feature_cube): "feature_size", "feature_mean", "feature_max", - "feature_effective_radius", + "feature_effective_diameter", ] for cube_name in expected_cube_names: From caa1b089c7724f4567cf890c943f0d559f55ba96 Mon Sep 17 00:00:00 2001 From: Adam Gainford Date: Thu, 17 Sep 2026 10:06:05 +0100 Subject: [PATCH 33/33] updates from code review (additional comments, catching scenario with no hzntl_coord in _get_effective_diameter function) --- src/CSET/operators/feature.py | 28 ++++++++++++------- src/CSET/operators/misc.py | 2 +- .../example_feature_cell_stats.yaml | 3 +- 3 files changed, 20 insertions(+), 13 deletions(-) diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py index 394b173a64..abba858399 100644 --- a/src/CSET/operators/feature.py +++ b/src/CSET/operators/feature.py @@ -485,20 +485,28 @@ def _get_effective_diameter_from_feature_size( """ # Guess coord representing horizontal grid (choose first available) - hzntl_coord = next( - iter( - [ - coord - for coord in cube_with_hzntl_coord.coords() - if iris.util.guess_coord_axis(coord) in ["X", "Y"] - ] + try: + hzntl_coord = next( + iter( + [ + coord + for coord in cube_with_hzntl_coord.coords() + if iris.util.guess_coord_axis(coord) in ["X", "Y"] + ] + ) ) - ) + except StopIteration: + raise ValueError( + "No horizontal coordinate found in input cube. " + "Cannot calculate effective diameter." + ) from None logger.debug(f"Attempting to convert to effective diameter using {hzntl_coord}") # Check coordinate is regular, but only warn if not, this is a naive estimate # and will be inaccurate for irregular grids + # Additionally, the current function only checks for regularity in one horizontal + # coordinate. Again, this is a bit of naive estimate. if not iris.util.is_regular(hzntl_coord): logger.warning( f"Horizontal coordinate {hzntl_coord} is not regular. " @@ -537,8 +545,8 @@ def _get_effective_diameter_from_feature_size( mean_latitude = 0 grid_spacing = grid_spacing * 111 * np.cos(np.radians(mean_latitude)) - effective_diameters_data = np.sqrt(size_data * grid_spacing**2 / np.pi) * 2 - return effective_diameters_data, grid_spacing + effective_diameter_data = np.sqrt(size_data * grid_spacing**2 / np.pi) * 2 + return effective_diameter_data, grid_spacing def _add_cell_stats_data_to_cubes( diff --git a/src/CSET/operators/misc.py b/src/CSET/operators/misc.py index 535f5311fe..46056db9c6 100644 --- a/src/CSET/operators/misc.py +++ b/src/CSET/operators/misc.py @@ -717,7 +717,7 @@ def flatten( cubes : iris.cube.Cube or iris.cube.CubeList The input Cube or CubeList to flatten. remove_nans : bool, optional - If True, remove NaN values from the flattened data. Default is True. + If True, remove NaN values from the flattened data. Default is False. Returns ------- diff --git a/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml b/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml index 1d7908ba85..aab22a63c4 100755 --- a/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml +++ b/src/CSET/recipes/example_recipes/example_feature_cell_stats.yaml @@ -1,8 +1,7 @@ category: Quick Look title: Example running cell_stats and plotting histograms description: | - Uses the feature.cell_stats operator to calculate cell size, mean cell values and max cell values for - identified features. Saves cell stats data to output and plots histograms + Uses the feature.cell_stats operator to calculate cell size, mean cell values and max cell values for identified features. Saves cell stats data to output and plots histograms steps: - operator: read.read_cubes