diff --git a/.gitignore b/.gitignore index 5fe83c2a9..590739681 100644 --- a/.gitignore +++ b/.gitignore @@ -138,3 +138,6 @@ dmypy.json # NFS synchronisation files .nfs* + +# MacOS temp files +.DS_Store diff --git a/pyproject.toml b/pyproject.toml index dbfa7fa59..4e169090c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,6 +35,7 @@ dependencies = [ "scores", "dask", "xarray", + "simple-track", ] [project.urls] diff --git a/requirements/environment.yml b/requirements/environment.yml index afe766bf6..7ab75d024 100644 --- a/requirements/environment.yml +++ b/requirements/environment.yml @@ -20,6 +20,7 @@ dependencies: - dask-core # Dask with minimal dependencies. - proj = 9.7.1 # Newer versions break plotting, see issue #2052. - matplotlib-base = 3.10.9 # Newer versions break plotting, see issue #2205. + - simple-track # For feature tracking and cell stats operators # Build dependencies - setuptools>=64 diff --git a/requirements/locks/py312-lock-linux-64.txt b/requirements/locks/py312-lock-linux-64.txt index 7fd65e438..cff0fa24e 100644 --- a/requirements/locks/py312-lock-linux-64.txt +++ b/requirements/locks/py312-lock-linux-64.txt @@ -1,6 +1,6 @@ # Generated by conda-lock. # platform: linux-64 -# input_hash: b250a7e669442d882266a2b36dbaf02ac22a344cf68c51b4909773b01bc8ed4e +# input_hash: 18f693332aade4640da699cf921afdf4a4ddc6a3c923cf3a3b804f97de3cd919 @EXPLICIT https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-20_gnu.conda#a9f577daf3de00bca7c3c76c0ecbd1de https://conda.anaconda.org/conda-forge/linux-64/aom-3.14.1-pl5321h039972f_1.conda#5a78a69eb3b50f24b379e9d2a93163ae @@ -283,6 +283,7 @@ https://conda.anaconda.org/conda-forge/noarch/scores-2.5.0-pyhd8ed1ab_0.conda#b8 https://conda.anaconda.org/conda-forge/noarch/setuptools-82.0.1-pyh332efcf_0.conda#8e194e7b992f99a5015edbd4ebd38efd https://conda.anaconda.org/conda-forge/noarch/setuptools-scm-10.2.0-pyhcf101f3_0.conda#224418e442ea786882979fbd2b36061f https://conda.anaconda.org/conda-forge/noarch/setuptools_scm-10.2.0-he0c3440_0.conda#cc27f4be260e32227b41d00e87d16b32 +https://conda.anaconda.org/conda-forge/noarch/simple-track-2.2.1-pyhcf101f3_0.conda#9e6b5bb8d63c55901c51efcc534c7442 https://conda.anaconda.org/conda-forge/noarch/six-1.17.0-pyhe01879c_1.conda#3339e3b65d58accf4ca4fb8748ab16b3 https://conda.anaconda.org/conda-forge/noarch/snowballstemmer-3.1.1-pyhd8ed1ab_0.conda#46b6abe31482f6bca064b965696ae807 https://conda.anaconda.org/conda-forge/noarch/soupsieve-2.8.4-pyhd8ed1ab_0.conda#9e21f087f087f805debe877d88e00a14 diff --git a/requirements/locks/py313-lock-linux-64.txt b/requirements/locks/py313-lock-linux-64.txt index 9d20aca63..c411b375f 100644 --- a/requirements/locks/py313-lock-linux-64.txt +++ b/requirements/locks/py313-lock-linux-64.txt @@ -1,6 +1,6 @@ # Generated by conda-lock. # platform: linux-64 -# input_hash: 9f964ad66b3ccd0df5e1897de428803405a71203d6fa546bb16f1fae4455533b +# input_hash: a4926971a72ffdc7239622d764c9e59e4be6aa95391ad6dd37f42b5bacc22cae @EXPLICIT https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-20_gnu.conda#a9f577daf3de00bca7c3c76c0ecbd1de https://conda.anaconda.org/conda-forge/linux-64/aom-3.14.1-pl5321h039972f_1.conda#5a78a69eb3b50f24b379e9d2a93163ae @@ -281,6 +281,7 @@ https://conda.anaconda.org/conda-forge/noarch/scores-2.5.0-pyhd8ed1ab_0.conda#b8 https://conda.anaconda.org/conda-forge/noarch/setuptools-82.0.1-pyh332efcf_0.conda#8e194e7b992f99a5015edbd4ebd38efd https://conda.anaconda.org/conda-forge/noarch/setuptools-scm-10.2.0-pyhcf101f3_0.conda#224418e442ea786882979fbd2b36061f https://conda.anaconda.org/conda-forge/noarch/setuptools_scm-10.2.0-he0c3440_0.conda#cc27f4be260e32227b41d00e87d16b32 +https://conda.anaconda.org/conda-forge/noarch/simple-track-2.2.1-pyhcf101f3_0.conda#9e6b5bb8d63c55901c51efcc534c7442 https://conda.anaconda.org/conda-forge/noarch/six-1.17.0-pyhe01879c_1.conda#3339e3b65d58accf4ca4fb8748ab16b3 https://conda.anaconda.org/conda-forge/noarch/snowballstemmer-3.1.1-pyhd8ed1ab_0.conda#46b6abe31482f6bca064b965696ae807 https://conda.anaconda.org/conda-forge/noarch/soupsieve-2.8.4-pyhd8ed1ab_0.conda#9e21f087f087f805debe877d88e00a14 diff --git a/requirements/locks/py314-lock-linux-64.txt b/requirements/locks/py314-lock-linux-64.txt index b3cc1552b..5e0a57b1d 100644 --- a/requirements/locks/py314-lock-linux-64.txt +++ b/requirements/locks/py314-lock-linux-64.txt @@ -1,6 +1,6 @@ # Generated by conda-lock. # platform: linux-64 -# input_hash: ce6655aa6b01babac128d5359fbf8cc501c833139bc7a07ec3a3d72835aea244 +# input_hash: 78675034effda826d4616dd32d2b3cf2f7ed14e76b8b9f35c009a784a9e21f24 @EXPLICIT https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-20_gnu.conda#a9f577daf3de00bca7c3c76c0ecbd1de https://conda.anaconda.org/conda-forge/linux-64/aom-3.14.1-pl5321h039972f_1.conda#5a78a69eb3b50f24b379e9d2a93163ae @@ -282,6 +282,7 @@ https://conda.anaconda.org/conda-forge/noarch/scores-2.5.0-pyhd8ed1ab_0.conda#b8 https://conda.anaconda.org/conda-forge/noarch/setuptools-82.0.1-pyh332efcf_0.conda#8e194e7b992f99a5015edbd4ebd38efd https://conda.anaconda.org/conda-forge/noarch/setuptools-scm-10.2.0-pyhcf101f3_0.conda#224418e442ea786882979fbd2b36061f https://conda.anaconda.org/conda-forge/noarch/setuptools_scm-10.2.0-he0c3440_0.conda#cc27f4be260e32227b41d00e87d16b32 +https://conda.anaconda.org/conda-forge/noarch/simple-track-2.2.1-pyhcf101f3_0.conda#9e6b5bb8d63c55901c51efcc534c7442 https://conda.anaconda.org/conda-forge/noarch/six-1.17.0-pyhe01879c_1.conda#3339e3b65d58accf4ca4fb8748ab16b3 https://conda.anaconda.org/conda-forge/noarch/snowballstemmer-3.1.1-pyhd8ed1ab_0.conda#46b6abe31482f6bca064b965696ae807 https://conda.anaconda.org/conda-forge/noarch/soupsieve-2.8.4-pyhd8ed1ab_0.conda#9e21f087f087f805debe877d88e00a14 diff --git a/requirements/locks/sources b/requirements/locks/sources index 9824ff3c4..8fe059b97 100644 --- a/requirements/locks/sources +++ b/requirements/locks/sources @@ -1 +1 @@ -97e468837e6c5cfff450fcad0310548b6d211f1148e4370f1e11130c3eb0de45 requirements/environment.yml +ef70f8db8d41711660879766b126d46c6b03d4e416e6ebec3598b64622a526c4 requirements/environment.yml diff --git a/src/CSET/operators/__init__.py b/src/CSET/operators/__init__.py index d5ef4afdf..354a01a9e 100644 --- a/src/CSET/operators/__init__.py +++ b/src/CSET/operators/__init__.py @@ -33,6 +33,7 @@ constraints, convection, ensembles, + feature, filters, humidity, imageprocessing, @@ -61,6 +62,7 @@ "convection", "ensembles", "execute_recipe", + "feature", "filters", "humidity", "get_operator", diff --git a/src/CSET/operators/_colorbar_definition.json b/src/CSET/operators/_colorbar_definition.json index 757219d00..39528db41 100644 --- a/src/CSET/operators/_colorbar_definition.json +++ b/src/CSET/operators/_colorbar_definition.json @@ -558,6 +558,26 @@ "max": 1e-05, "min": -1e-05 }, + "precipitation_flux": { + "cmap": "cividis", + "levels": [ + 0, + 0.125, + 0.25, + 0.5, + 1, + 2, + 4, + 8, + 16, + 32, + 64, + 128, + 256 + ], + "ymax": 1.0, + "ymin": 0.0 + }, "radar_reflectivity_at_1km_above_the_surface": { "cmap": "cubehelix_r", "max": 70.0, diff --git a/src/CSET/operators/_colormaps.py b/src/CSET/operators/_colormaps.py index efe1d83bd..74e684aa9 100644 --- a/src/CSET/operators/_colormaps.py +++ b/src/CSET/operators/_colormaps.py @@ -240,6 +240,7 @@ def colorbar_map_levels(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = N cmap, levels, norm = custom_colormap_precipitation(cube, cmap, levels, norm) cmap, levels, norm = custom_colormap_visibility_in_air(cube, cmap, levels, norm) cmap, levels, norm = custom_colormap_celsius(cube, cmap, levels, norm) + cmap, levels, norm = custom_colormap_feature_tracking(cube, cmap, levels, norm) return cmap, levels, norm @@ -597,3 +598,67 @@ def custom_colormap_scores(cube: iris.cube.Cube): if any("RMSE_" in name for name in varnames): cmap = plt.get_cmap("PuRd", 51) return cmap, levels, norm + + +def custom_colormap_feature_tracking(cube: iris.cube.Cube, cmap, levels, norm): + """Return altered colormap for feature tracking. + + Parameters + ---------- + cube: Cube + Cube of variable for which the colorbar information is desired. + + Returns + ------- + cmap: Matplotlib colormap. + levels: List + List of levels to use for plotting. For continuous plots the min and max + should be taken as the range. + norm: BoundaryNorm. + """ + varnames = list(filter(None, [cube.long_name, cube.standard_name, cube.var_name])) + if ( + any("feature_id" in name for name in varnames) + and "difference" not in cube.long_name + and "mask" not in cube.long_name + ): + # Define the levels and colors + levels = np.linspace(1, np.ma.max(cube.data), 10) + cmap = plt.get_cmap("viridis") + # Normalize the levels + norm = mcolors.BoundaryNorm(levels, cmap.N) + logging.info("change colormap for feature id variable colorbar.") + elif ( + any("feature_lifetime" in name for name in varnames) + and "difference" not in cube.long_name + and "mask" not in cube.long_name + ): + # Define the levels and colors + levels = np.linspace(1, np.ma.max(cube.data), 10) + cmap = plt.get_cmap("YlGnBu") + # Normalize the levels + norm = mcolors.BoundaryNorm(levels, cmap.N) + logging.info("change colormap for feature lifetime variable colorbar.") + elif ( + any("feature_init" in name for name in varnames) + and "difference" not in cube.long_name + and "mask" not in cube.long_name + ): + # Define the levels and colors + levels = np.array([0.5, 1]) + cmap = plt.get_cmap("Blues") + # Normalize the levels + norm = mcolors.BoundaryNorm(levels, cmap.N) + logging.info("change colormap for feature init variable colorbar.") + + else: + # do nothing and keep existing colorbar attributes + cmap = cmap + levels = levels + norm = norm + + # Set all non-feature data to white + if any("feature" in name for name in varnames): + cmap.with_extremes(under="white") + + return cmap, levels, norm diff --git a/src/CSET/operators/feature.py b/src/CSET/operators/feature.py new file mode 100755 index 000000000..3eb594e4b --- /dev/null +++ b/src/CSET/operators/feature.py @@ -0,0 +1,222 @@ +# © Crown copyright, Met Office (2022-2026) 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 iris.cube +import iris.util +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 + ---------- + 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 + 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 feature 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 spawned 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() + + """ + # Check that the input cube has horizontal coordinates of xy type, not latitude/longitude + _check_xy_coords(cube) + + # 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 + + +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 horizontal coordinates. + + 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." + ) 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 000000000..240119047 --- /dev/null +++ b/src/CSET/recipes/example_recipes/example_feature_track.yaml @@ -0,0 +1,26 @@ +category: Quick Look +title: Precipitation flux feature lifetime spatial plot +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 diff --git a/tests/operators/test_colormaps.py b/tests/operators/test_colormaps.py index 6789c23f5..2d8689f78 100644 --- a/tests/operators/test_colormaps.py +++ b/tests/operators/test_colormaps.py @@ -403,3 +403,42 @@ def test_colorbar_map_scores_rmse(cube, tmp_working_dir): assert cmap == plt.get_cmap("PuRd", 51) assert levels is None assert norm is None + + +def test_colorbar_feature_tracking_id_cube(cube): + """Colorbar definition is found for a feature id cube.""" + cube.rename("feature_id") + cmap, levels, norm = _colormaps.custom_colormap_feature_tracking( + cube, None, None, None + ) + expected_levels = np.linspace(1, np.ma.max(cube.data), 10) + assert cmap == plt.get_cmap("viridis") + assert cmap.get_under() is not None + assert (levels == expected_levels).all() + assert (norm.boundaries == levels).all() + + +def test_colorbar_feature_tracking_lifetime_cube(cube): + """Colorbar definition is found for a feature lifetime cube.""" + cube.rename("feature_lifetime") + cmap, levels, norm = _colormaps.custom_colormap_feature_tracking( + cube, None, None, None + ) + expected_levels = np.linspace(1, np.ma.max(cube.data), 10) + assert cmap == plt.get_cmap("YlGnBu") + assert cmap.get_under() is not None + assert (levels == expected_levels).all() + assert (norm.boundaries == levels).all() + + +def test_colorbar_feature_tracking_init_cube(cube): + """Colorbar definition is found for a feature init cube.""" + cube.rename("feature_init") + cmap, levels, norm = _colormaps.custom_colormap_feature_tracking( + cube, None, None, None + ) + expected_levels = np.array([0.5, 1]) + assert cmap == plt.get_cmap("Blues") + assert cmap.get_under() is not None + assert (levels == expected_levels).all() + assert (norm.boundaries == levels).all() diff --git a/tests/operators/test_feature.py b/tests/operators/test_feature.py new file mode 100644 index 000000000..54e144f69 --- /dev/null +++ b/tests/operators/test_feature.py @@ -0,0 +1,165 @@ +# © Crown copyright, Met Office (2022-2026) 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 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_working_dir) -> None: + """Test that tracking data is saved when save_data is True.""" + 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 = tmp_working_dir / "tracking_data" + expected_file = output_directory / "lifetime_20100101_0000.field" + assert expected_file.is_file() + + # Check expected csv file is created in output directory + expected_file = output_directory / "frame_20100101_0000.csv" + 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_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 + data_arr = np.zeros((10, 10)) + + lat_coord = iris.coords.DimCoord( + points=np.linspace(-90, 90, 10), + standard_name="latitude", + var_name="latitude", + units="degrees", + ) + lon_coord = iris.coords.DimCoord( + points=np.linspace(-180, 180, 10), + 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", + ) + + with pytest.raises(ValueError): + feature._check_xy_coords(cube)