From b39300799c3038fcbb512af80cf9b351c655d463 Mon Sep 17 00:00:00 2001 From: daflack Date: Thu, 30 Apr 2026 12:25:51 +0100 Subject: [PATCH 01/17] Adds a global model application for CURV Fixes #2089 --- src/CSET/operators/curvature.py | 99 +++++++++++++++++++++++++++++++++ 1 file changed, 99 insertions(+) create mode 100644 src/CSET/operators/curvature.py diff --git a/src/CSET/operators/curvature.py b/src/CSET/operators/curvature.py new file mode 100644 index 000000000..24ae37de2 --- /dev/null +++ b/src/CSET/operators/curvature.py @@ -0,0 +1,99 @@ +# © 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. + +"""A module containing different implementations of CURV. + +The diagnostics here are separate implementations of CURV for both +global models and limited-area models. +""" + +import numpy as np + + +def lat_lon_identifier(original_lat, original_lon, distance, bearing): + """Find the latitude and longitude of new points based on great circle distance.""" + R = 6371.0 + lat_rad = np.radians(original_lat) + lon_rad = np.radians(original_lon) + bearing_rad = np.radians(bearing) + + new_lat = np.arcsin( + np.sin(lat_rad) * np.cos(distance / R) + + np.cos(lat_rad) * np.sin(distance / R) * np.cos(bearing_rad) + ) + new_lon = lon_rad + np.arctan2( + np.sin(bearing_rad) * np.sin(distance / R) * np.cos(lat_rad), + np.cos(distance / R) - np.sin(lat_rad) * np.sin(new_lat), + ) + + new_lon = np.where(new_lon > 180.0, new_lon - 360.0, new_lon) + + return new_lat, new_lon + + +# def curvature_function(central, surrounding, tolerance, radial_points): +# """Calculate the intermediate data for CURV diagnostic.""" +# shape_of_surrounding = surrounding.shape +# binary_difference = np.zeros(shape_of_surrounding) +# for rp in range(radial_points): +# binary_difference[:, rp] = np.where( +# surrounding[:, rp] < central[:] + tolerance, 1, 0 +# ) +# binary_difference[:, rp] = np.where( +# surroudnig[:, rp] > central[:] - tolerance, -1, binary_difference[:, rp] +# ) +# curvature = bindary_difference.sum(axis=1) +# return curvature + + +# def curv( +# radius, num_radial_points, north, south, east, west, tolerance, cube, outputfile +# ): +# """Calculate the CURV diagnostic.""" +# if (num_radial_points != 8) or (num_radial_points != 16): +# raise ValueError( +# f"Number of radial points should be 8 or 16. You have provided {num_radial_points}." +# ) +# else: +# bearing = np.linspace(0.0, 360.0, num_radial_points + 1) +# +# data_area = [north, west, south, east] +# +# logging.info( +# f"CURV calculated base on {num_radial_points} radial points, a {radius} km radius, with {tolerance} tolerance." +# ) +# logging.info(f"CURV calculated over area {data_area}.") +# +# extracted_data = cube.extract(lon) +# extract_lats = extract_data.extract(lat) +# +# points = np.size(extract_data.core_data()) +# +# curv = np.zeros(points) +# +# new_lats = np.zeros((points, num_radial_points)) +# new_lons = np.zeros((points, num_radial_points)) +# +# new_vals = np.zeros((points, num_radial_points)) +# +# for i in range(num_radial_points): +# new_lats[:, i], new_lons[:, i] = lat_lon_identifier( +# lats, lons, radius, bearing[i] +# ) +# +# new_vals[:, i] = cube.interpolate(new_lats[:, i], new_lons[:, i]) +# +# curv = curvature_function(cube, new_vals, tolerance, num_radial_points) +# +# return curv From 2e85ea50404bedfd4b13bfaaafff1fec8791ebcf Mon Sep 17 00:00:00 2001 From: daflack Date: Thu, 30 Apr 2026 12:30:16 +0100 Subject: [PATCH 02/17] Adds operators into CSET operators init and documentation --- docs/source/reference/operators.rst | 9 +++++++++ src/CSET/operators/__init__.py | 2 ++ 2 files changed, 11 insertions(+) diff --git a/docs/source/reference/operators.rst b/docs/source/reference/operators.rst index db6ab6aa1..a2da92d35 100644 --- a/docs/source/reference/operators.rst +++ b/docs/source/reference/operators.rst @@ -86,6 +86,15 @@ CSET.operators.convection .. automodule:: CSET.operators.convection :members: +Curvature Operators +~~~~~~~~~~~~~~~~~~~ + +CSET.operators.curvature +------------------------ + +.. automodule:: CSET.operators.curvature + :members: + Ensemble Operators ~~~~~~~~~~~~~~~~~~ diff --git a/src/CSET/operators/__init__.py b/src/CSET/operators/__init__.py index cf583c30d..b8ab7812e 100644 --- a/src/CSET/operators/__init__.py +++ b/src/CSET/operators/__init__.py @@ -32,6 +32,7 @@ collapse, constraints, convection, + curvature, ensembles, filters, humidity, @@ -56,6 +57,7 @@ "collapse", "constraints", "convection", + "curvature", "ensembles", "execute_recipe", "filters", From f86b5d9f2089b68a6ccee2d46412733d464ea683 Mon Sep 17 00:00:00 2001 From: daflack Date: Thu, 30 Apr 2026 16:02:22 +0100 Subject: [PATCH 03/17] Near full version of a slow CURV operator --- src/CSET/operators/curvature.py | 113 ++++++++++++++++---------------- 1 file changed, 55 insertions(+), 58 deletions(-) diff --git a/src/CSET/operators/curvature.py b/src/CSET/operators/curvature.py index 24ae37de2..4fd7511a1 100644 --- a/src/CSET/operators/curvature.py +++ b/src/CSET/operators/curvature.py @@ -18,12 +18,15 @@ global models and limited-area models. """ +import logging + +import iris import numpy as np -def lat_lon_identifier(original_lat, original_lon, distance, bearing): +def _lat_lon_identifier(original_lat, original_lon, distance, bearing): """Find the latitude and longitude of new points based on great circle distance.""" - R = 6371.0 + R = 6371.0 # km. lat_rad = np.radians(original_lat) lon_rad = np.radians(original_lon) bearing_rad = np.radians(bearing) @@ -37,63 +40,57 @@ def lat_lon_identifier(original_lat, original_lon, distance, bearing): np.cos(distance / R) - np.sin(lat_rad) * np.sin(new_lat), ) - new_lon = np.where(new_lon > 180.0, new_lon - 360.0, new_lon) + new_lat = np.degrees(new_lat) + new_lon = np.degrees(new_lon) + new_lon = np.where(new_lon > 360.0, new_lon - 360.0, new_lon) return new_lat, new_lon -# def curvature_function(central, surrounding, tolerance, radial_points): -# """Calculate the intermediate data for CURV diagnostic.""" -# shape_of_surrounding = surrounding.shape -# binary_difference = np.zeros(shape_of_surrounding) -# for rp in range(radial_points): -# binary_difference[:, rp] = np.where( -# surrounding[:, rp] < central[:] + tolerance, 1, 0 -# ) -# binary_difference[:, rp] = np.where( -# surroudnig[:, rp] > central[:] - tolerance, -1, binary_difference[:, rp] -# ) -# curvature = bindary_difference.sum(axis=1) -# return curvature - - -# def curv( -# radius, num_radial_points, north, south, east, west, tolerance, cube, outputfile -# ): -# """Calculate the CURV diagnostic.""" -# if (num_radial_points != 8) or (num_radial_points != 16): -# raise ValueError( -# f"Number of radial points should be 8 or 16. You have provided {num_radial_points}." -# ) -# else: -# bearing = np.linspace(0.0, 360.0, num_radial_points + 1) -# -# data_area = [north, west, south, east] -# -# logging.info( -# f"CURV calculated base on {num_radial_points} radial points, a {radius} km radius, with {tolerance} tolerance." -# ) -# logging.info(f"CURV calculated over area {data_area}.") -# -# extracted_data = cube.extract(lon) -# extract_lats = extract_data.extract(lat) -# -# points = np.size(extract_data.core_data()) -# -# curv = np.zeros(points) -# -# new_lats = np.zeros((points, num_radial_points)) -# new_lons = np.zeros((points, num_radial_points)) -# -# new_vals = np.zeros((points, num_radial_points)) -# -# for i in range(num_radial_points): -# new_lats[:, i], new_lons[:, i] = lat_lon_identifier( -# lats, lons, radius, bearing[i] -# ) -# -# new_vals[:, i] = cube.interpolate(new_lats[:, i], new_lons[:, i]) -# -# curv = curvature_function(cube, new_vals, tolerance, num_radial_points) -# -# return curv +def curv(central, radius, num_radial_points=16, tol=0): + """Calculate the CURV diagnostic.""" + if (num_radial_points == 8) or (num_radial_points == 16): + bearing = np.linspace(0.0, 360.0, num_radial_points + 1)[:-1] + else: + raise ValueError( + f"Number of radial points should be 8 or 16. You have provided {num_radial_points}." + ) + + logging.info( + f"CURV calculated base on {num_radial_points} radial points, a {radius} km radius and {tol} tolerance." + ) + + surroundings = iris.cube.CubeList([]) + for b in bearing: + surround = central.copy() + # Need to try and avoid this loop for whole domain if possible, also need to raise time and realization bits as well. + for lat_number, lat in enumerate(central.slices_over("latitude")): + for lon_number, lon in enumerate(lat.slices_over("longitude")): + new_lat, new_lon = _lat_lon_identifier( + lon.coord("latitude").points, + lon.coord("longitude").points, + radius, + b, + ) + # Need to get interpolation better. + surround.data[lat_number, lon_number] = central.data[ + np.abs(central.coord("latitude").points - new_lat[0]).argmin(), + np.abs(central.coord("longitude").points - new_lon[0]).argmin(), + ] + surround.add_aux_coord( + iris.coords.DimCoord(b, long_name="bearing", units="degrees") + ) + surroundings.append(surround) + surroundings.merge() + print(surroundings) + + # curv = surroundings.copy() + # curv -= central + + # curv.data[curv.data > tol]=-1.0 + # curv.data[curv.data < tol]=1.0 + + # curv.collapsed('bearing',iris.analysis.SUM) + # curv.rename("CURV") + # curv.units="1" + # return curv From 822af6e514cf998f72a071b3a360f145dde21baf Mon Sep 17 00:00:00 2001 From: daflack Date: Thu, 30 Apr 2026 16:05:39 +0100 Subject: [PATCH 04/17] Update to a comment for code if current testing works --- src/CSET/operators/curvature.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/operators/curvature.py b/src/CSET/operators/curvature.py index 4fd7511a1..2a7483f72 100644 --- a/src/CSET/operators/curvature.py +++ b/src/CSET/operators/curvature.py @@ -84,7 +84,7 @@ def curv(central, radius, num_radial_points=16, tol=0): surroundings.merge() print(surroundings) - # curv = surroundings.copy() + # curv = surroundings[0].copy() # curv -= central # curv.data[curv.data > tol]=-1.0 From c47033fe071d4a88cf408c8542876cba232f55bf Mon Sep 17 00:00:00 2001 From: daflack Date: Thu, 30 Apr 2026 16:09:17 +0100 Subject: [PATCH 05/17] Update to basic CURV operator --- src/CSET/operators/curvature.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/CSET/operators/curvature.py b/src/CSET/operators/curvature.py index 2a7483f72..2dca2fb06 100644 --- a/src/CSET/operators/curvature.py +++ b/src/CSET/operators/curvature.py @@ -84,13 +84,13 @@ def curv(central, radius, num_radial_points=16, tol=0): surroundings.merge() print(surroundings) - # curv = surroundings[0].copy() - # curv -= central + curv = surroundings[0].copy() + curv -= central - # curv.data[curv.data > tol]=-1.0 - # curv.data[curv.data < tol]=1.0 + curv.data[curv.data > tol] = -1.0 + curv.data[curv.data < tol] = 1.0 - # curv.collapsed('bearing',iris.analysis.SUM) - # curv.rename("CURV") - # curv.units="1" - # return curv + curv.collapsed("bearing", iris.analysis.SUM) + curv.rename(f"CURV_calculated_from_{num_radial_points}_radial_points") + curv.units = "1" + return curv From 4c8c0952a327ffa153e249a198defb63760d5b9a Mon Sep 17 00:00:00 2001 From: daflack Date: Thu, 30 Apr 2026 16:19:52 +0100 Subject: [PATCH 06/17] Update CURV operator to show it is for global CURV --- src/CSET/operators/curvature.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/CSET/operators/curvature.py b/src/CSET/operators/curvature.py index 2dca2fb06..893060bad 100644 --- a/src/CSET/operators/curvature.py +++ b/src/CSET/operators/curvature.py @@ -47,8 +47,8 @@ def _lat_lon_identifier(original_lat, original_lon, distance, bearing): return new_lat, new_lon -def curv(central, radius, num_radial_points=16, tol=0): - """Calculate the CURV diagnostic.""" +def global_curv(central, radius, num_radial_points=16, tol=0): + """Calculate the CURV diagnostic for a global model.""" if (num_radial_points == 8) or (num_radial_points == 16): bearing = np.linspace(0.0, 360.0, num_radial_points + 1)[:-1] else: From 3cd964b99b0ebb6fdf5d2d4df38a655eed859b69 Mon Sep 17 00:00:00 2001 From: daflack Date: Wed, 13 May 2026 12:13:16 +0100 Subject: [PATCH 07/17] Adds basic recipe and plumbing into workflow This is mainly for testing purposes to see if it works. --- .../meta/diagnostics/rose-meta.conf | 12 +++++++++ src/CSET/loaders/spatial_field.py | 15 +++++++++++ .../daily_weather/curv_spatial_plot.yaml | 26 +++++++++++++++++++ 3 files changed, 53 insertions(+) create mode 100644 src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml diff --git a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf index e3fbaace3..84834adc6 100644 --- a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf @@ -1327,6 +1327,18 @@ type=python_boolean compulsory=true sort-key=dw-fog6 +[template variables = GLOBAL_CURV] +ns=Diagnostics/Derived/Daily +title=Global CURV: spatial plot +description=PROCESS-BASED DIAGNOSTIC. + Determines the curvature of the MSLP field for the global model. + Requires air_pressure_at_mean_sea_level. + TEMPORARY LOCATION FOR TESTING PURPOSES. +help=This diagnostic can be particularly useful for grouping of the data if + used within a recipe. This provides the plotted data. +type=python_boolean +compulsory=true +sort-key=dw-curv ################################### # Aviation [Diagnostics/Aviation] diff --git a/src/CSET/loaders/spatial_field.py b/src/CSET/loaders/spatial_field.py index ba35ca1e2..8e24d2c2a 100644 --- a/src/CSET/loaders/spatial_field.py +++ b/src/CSET/loaders/spatial_field.py @@ -572,6 +572,21 @@ def load(conf: Config): aggregation=False, ) + if conf.GLOBAL_CURV: + for model in models: + yield RawRecipe( + recipe="curv_spatial_plot.yaml", + variables={ + "MODEL_NAME": model["name"], + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT + if conf.SELECT_SUBAREA + else None, + }, + model_ids=model["id"], + aggregation=False, + ) + # Screen-level temperature probabilities for model, condition, threshold in itertools.product( models, diff --git a/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml b/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml new file mode 100644 index 000000000..57b28a0da --- /dev/null +++ b/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml @@ -0,0 +1,26 @@ +category: Daily Weather +title: $MODEL_NAME CURV spatial plot +description: | + Generates the curvature (CURV) from MSLP and creates a spatila plot. +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS + model_names: $MODEL_NAME + constraint: + operator: constraints.combine_constraints + varname_constraint: + operator: constraints.generate_var_constraint + varname: air_pressure_at_mean_sea_level + cell_methods_constraint: + operator: constraints.generate_cell_methods_constraint + cell_methods: [] + subarea_type: $SUBAREA_TYPE + subarea_extent: $SUBAREA_EXTENT + + - operator: curvature.curv + radius: 1000 + + - operator: plot.spatial_pcolormesh_plot + + - operator: write.write_cube_to_nc + overwrite: True From bcfce1f87f5640bb4e913e8a768960f7479e5291 Mon Sep 17 00:00:00 2001 From: daflack Date: Wed, 13 May 2026 12:18:41 +0100 Subject: [PATCH 08/17] Adds CURV to example conf --- src/CSET/cset_workflow/rose-suite.conf.example | 1 + 1 file changed, 1 insertion(+) diff --git a/src/CSET/cset_workflow/rose-suite.conf.example b/src/CSET/cset_workflow/rose-suite.conf.example index f8a730649..1997057e5 100644 --- a/src/CSET/cset_workflow/rose-suite.conf.example +++ b/src/CSET/cset_workflow/rose-suite.conf.example @@ -42,6 +42,7 @@ EXTRACT_PLEVEL_TRANSECT=False FOG_PRESENCE_DOMAIN_MEAN_TIMESERIES=False FOG_PRESENCE_SPATIAL_DIFFERENCE=False FOG_PRESENCE_SPATIAL_PLOT=False +GLOBAL_CURV=False GROUND_FROST_PRESENCE_DOMAIN_MEAN_TIMESERIES=False GROUND_FROST_PRESENCE_SPATIAL_DIFFERENCE=False GROUND_FROST_PRESENCE_SPATIAL_PLOT=False From 46b8843197d977b1a4a25aab677257418a33fad8 Mon Sep 17 00:00:00 2001 From: daflack Date: Wed, 13 May 2026 12:20:15 +0100 Subject: [PATCH 09/17] Update rose-meta.conf --- src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf index 84834adc6..27b470a9a 100644 --- a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf @@ -1327,7 +1327,7 @@ type=python_boolean compulsory=true sort-key=dw-fog6 -[template variables = GLOBAL_CURV] +[template variables=GLOBAL_CURV] ns=Diagnostics/Derived/Daily title=Global CURV: spatial plot description=PROCESS-BASED DIAGNOSTIC. From abeb3aceecd79304bbba5b0e222995d191af1b0e Mon Sep 17 00:00:00 2001 From: daflack Date: Wed, 13 May 2026 15:44:25 +0100 Subject: [PATCH 10/17] Updates recipe to use the faster operator, and creates a tentative fast operator This commit contains the first working version of global_curv in global_curv_version2, global_curv is slow and has not yet produced output. --- src/CSET/operators/curvature.py | 122 +++++++++++++----- .../daily_weather/curv_spatial_plot.yaml | 3 +- 2 files changed, 92 insertions(+), 33 deletions(-) diff --git a/src/CSET/operators/curvature.py b/src/CSET/operators/curvature.py index 893060bad..ce009ca7e 100644 --- a/src/CSET/operators/curvature.py +++ b/src/CSET/operators/curvature.py @@ -23,6 +23,8 @@ import iris import numpy as np +from CSET._common import iter_maybe + def _lat_lon_identifier(original_lat, original_lon, distance, bearing): """Find the latitude and longitude of new points based on great circle distance.""" @@ -59,38 +61,94 @@ def global_curv(central, radius, num_radial_points=16, tol=0): logging.info( f"CURV calculated base on {num_radial_points} radial points, a {radius} km radius and {tol} tolerance." ) + curv_cubes = iris.cube.CubeList([]) + for cube in iter_maybe(central): + surroundings = iris.cube.CubeList([]) + for b in bearing: + surround = cube.copy() + # Need to try and avoid this loop for whole domain if possible, also need to raise time and realization bits as well. + for lat_number, lat in enumerate(cube.slices_over("latitude")): + for lon_number, lon in enumerate(lat.slices_over("longitude")): + new_lat, new_lon = _lat_lon_identifier( + lon.coord("latitude").points, + lon.coord("longitude").points, + radius, + b, + ) + # Need to get interpolation better. + surround.data[:, lat_number, lon_number] = cube.data[ + :, + np.abs(cube.coord("latitude").points - new_lat[0]).argmin(), + np.abs(cube.coord("longitude").points - new_lon[0]).argmin(), + ] + surround.add_aux_coord( + iris.coords.DimCoord(b, long_name="bearing", units="degrees") + ) + surroundings.append(surround) + surroundings.merge() + print(surroundings) + + curv = surroundings[0].copy() + curv -= cube + + curv.data[curv.data > tol] = -1.0 + curv.data[curv.data < tol] = 1.0 + + curv.collapsed("bearing", iris.analysis.SUM) + curv.rename(f"CURV_calculated_from_{num_radial_points}_radial_points") + curv.units = "1" + curv_cubes.append(curv) + if len(curv_cubes) == 1: + return curv_cubes[0] + else: + return curv_cubes - surroundings = iris.cube.CubeList([]) - for b in bearing: - surround = central.copy() - # Need to try and avoid this loop for whole domain if possible, also need to raise time and realization bits as well. - for lat_number, lat in enumerate(central.slices_over("latitude")): - for lon_number, lon in enumerate(lat.slices_over("longitude")): - new_lat, new_lon = _lat_lon_identifier( - lon.coord("latitude").points, - lon.coord("longitude").points, - radius, - b, - ) - # Need to get interpolation better. - surround.data[lat_number, lon_number] = central.data[ - np.abs(central.coord("latitude").points - new_lat[0]).argmin(), - np.abs(central.coord("longitude").points - new_lon[0]).argmin(), - ] - surround.add_aux_coord( - iris.coords.DimCoord(b, long_name="bearing", units="degrees") - ) - surroundings.append(surround) - surroundings.merge() - print(surroundings) - - curv = surroundings[0].copy() - curv -= central - curv.data[curv.data > tol] = -1.0 - curv.data[curv.data < tol] = 1.0 +def global_curv_versiontwo(central, radius, num_radial_points=16, tol=0.0): + """Calculate the CURV diagnostic for a global model VERSION 2 for faster calculation speed.""" + if (num_radial_points == 8) or (num_radial_points == 16): + bearing = np.linspace(0.0, 360.0, num_radial_points + 1)[:-1] + else: + raise ValueError( + f"Number of radial points should be 8 or 16. You have provided {num_radial_points}." + ) - curv.collapsed("bearing", iris.analysis.SUM) - curv.rename(f"CURV_calculated_from_{num_radial_points}_radial_points") - curv.units = "1" - return curv + logging.info( + f"CURV calculated base on {num_radial_points} radial points, a {radius} km radius and {tol} tolerance." + ) + curv_cubes = iris.cube.CubeList([]) + for cube in iter_maybe(central): + curv = cube.copy() + surroundings = [] + lats, lons = np.meshgrid( + cube.coord("latitude").points, cube.coord("longitude").points + ) + for b in bearing: + target_data = cube.copy() + new_lat, new_lon = _lat_lon_identifier( + lats, + lons, + radius, + b, + ) + new_lat_points = new_lat[0, :] + new_lon_points = new_lon[:, 0] + surround = target_data.interpolate( + [("latitude", new_lat_points), ("longitude", new_lon_points)], + iris.analysis.Nearest(), + ) + surroundings.append(surround.data) + surroundings = np.array(surroundings) + curv_data = surroundings - cube.data + curv_data_binary = np.zeros(curv_data.shape) + curv_data_binary[curv_data > tol] = -1.0 + curv_data_binary[curv_data < tol] = 1.0 + curv_data_collapsed = np.sum(curv_data_binary, 0) + curv.data = curv_data_collapsed + curv.rename(f"CURV_calculated_from_{num_radial_points}_radial_points") + curv.units = "1" + curv_cubes.append(curv) + if len(curv_cubes) == 1: + return curv_cubes[0] + else: + return curv_cubes diff --git a/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml b/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml index 57b28a0da..db97a3f16 100644 --- a/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml +++ b/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml @@ -17,8 +17,9 @@ steps: subarea_type: $SUBAREA_TYPE subarea_extent: $SUBAREA_EXTENT - - operator: curvature.curv + - operator: curvature.global_curv_versiontwo radius: 1000 + num_radial_points: 8 - operator: plot.spatial_pcolormesh_plot From db3fd66b1924a3640ebbc5392c3c39d7bee2e438 Mon Sep 17 00:00:00 2001 From: daflack Date: Thu, 14 May 2026 13:27:26 +0100 Subject: [PATCH 11/17] Switches to linear interpolation and adds options to change radius and number of points in gui --- .../meta/diagnostics/rose-meta.conf | 18 +++++++++++++++++- src/CSET/loaders/spatial_field.py | 6 +++++- src/CSET/operators/curvature.py | 2 +- .../daily_weather/curv_spatial_plot.yaml | 10 ++++++---- 4 files changed, 29 insertions(+), 7 deletions(-) diff --git a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf index 27b470a9a..1c53be724 100644 --- a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf @@ -1338,7 +1338,23 @@ help=This diagnostic can be particularly useful for grouping of the data if used within a recipe. This provides the plotted data. type=python_boolean compulsory=true -sort-key=dw-curv +sort-key=dw-curv3 + +[template variables=CURV_POINTS] +ns=Diagnostics/Derived/Daily +title=Number of points to use in CURV +description=The number of points to use in CURV calculation. This should be 8 or 16. +type=python_list +compulsory=true +sort-key=dw-curv1 + +[template variables=CURV_RADIUS] +ns=Diagnostics/Derived/Daily +title=Radii to use in CURV +description=The radii in km to use in CURV calculation. +type=python_list +compulsory=true +sort-key=dw-curv2 ################################### # Aviation [Diagnostics/Aviation] diff --git a/src/CSET/loaders/spatial_field.py b/src/CSET/loaders/spatial_field.py index 8e24d2c2a..aa91b79ab 100644 --- a/src/CSET/loaders/spatial_field.py +++ b/src/CSET/loaders/spatial_field.py @@ -573,11 +573,15 @@ def load(conf: Config): ) if conf.GLOBAL_CURV: - for model in models: + for model, points, radius in itertools.product( + models, conf.CURV_POINTS, conf.CURV_RADIUS + ): yield RawRecipe( recipe="curv_spatial_plot.yaml", variables={ "MODEL_NAME": model["name"], + "CURV_POINTS": points, + "CURV_RADIUS": radius, "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA diff --git a/src/CSET/operators/curvature.py b/src/CSET/operators/curvature.py index ce009ca7e..3b30c3f44 100644 --- a/src/CSET/operators/curvature.py +++ b/src/CSET/operators/curvature.py @@ -135,7 +135,7 @@ def global_curv_versiontwo(central, radius, num_radial_points=16, tol=0.0): new_lon_points = new_lon[:, 0] surround = target_data.interpolate( [("latitude", new_lat_points), ("longitude", new_lon_points)], - iris.analysis.Nearest(), + iris.analysis.Linear(), ) surroundings.append(surround.data) surroundings = np.array(surroundings) diff --git a/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml b/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml index db97a3f16..87e55b2d3 100644 --- a/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml +++ b/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml @@ -1,7 +1,9 @@ category: Daily Weather -title: $MODEL_NAME CURV spatial plot +title: $MODEL_NAME CURV spatial plot from $CURV_POINTS points and radius of $CURV_RADIUS km description: | - Generates the curvature (CURV) from MSLP and creates a spatila plot. + Generates the curvature (CURV) from MSLP and creates a spatial plot. + It uses $CURV_POINTS points and a radius of $CURV_RADIUS km. + +$CURV_POINTS is a full anticyclone and -$CURV_POINTS is a full cyclone. steps: - operator: read.read_cubes file_paths: $INPUT_PATHS @@ -18,8 +20,8 @@ steps: subarea_extent: $SUBAREA_EXTENT - operator: curvature.global_curv_versiontwo - radius: 1000 - num_radial_points: 8 + radius: $CURV_RADIUS + num_radial_points: $CURV_POINTS - operator: plot.spatial_pcolormesh_plot From 08a5049fef2d62be5601f2b67e7ed2a87d697913 Mon Sep 17 00:00:00 2001 From: daflack Date: Thu, 14 May 2026 13:29:53 +0100 Subject: [PATCH 12/17] Removes slow version of code, renames fast version, updates recipe --- src/CSET/operators/curvature.py | 57 +------------------ .../daily_weather/curv_spatial_plot.yaml | 2 +- 2 files changed, 2 insertions(+), 57 deletions(-) diff --git a/src/CSET/operators/curvature.py b/src/CSET/operators/curvature.py index 3b30c3f44..c29b75c27 100644 --- a/src/CSET/operators/curvature.py +++ b/src/CSET/operators/curvature.py @@ -49,62 +49,7 @@ def _lat_lon_identifier(original_lat, original_lon, distance, bearing): return new_lat, new_lon -def global_curv(central, radius, num_radial_points=16, tol=0): - """Calculate the CURV diagnostic for a global model.""" - if (num_radial_points == 8) or (num_radial_points == 16): - bearing = np.linspace(0.0, 360.0, num_radial_points + 1)[:-1] - else: - raise ValueError( - f"Number of radial points should be 8 or 16. You have provided {num_radial_points}." - ) - - logging.info( - f"CURV calculated base on {num_radial_points} radial points, a {radius} km radius and {tol} tolerance." - ) - curv_cubes = iris.cube.CubeList([]) - for cube in iter_maybe(central): - surroundings = iris.cube.CubeList([]) - for b in bearing: - surround = cube.copy() - # Need to try and avoid this loop for whole domain if possible, also need to raise time and realization bits as well. - for lat_number, lat in enumerate(cube.slices_over("latitude")): - for lon_number, lon in enumerate(lat.slices_over("longitude")): - new_lat, new_lon = _lat_lon_identifier( - lon.coord("latitude").points, - lon.coord("longitude").points, - radius, - b, - ) - # Need to get interpolation better. - surround.data[:, lat_number, lon_number] = cube.data[ - :, - np.abs(cube.coord("latitude").points - new_lat[0]).argmin(), - np.abs(cube.coord("longitude").points - new_lon[0]).argmin(), - ] - surround.add_aux_coord( - iris.coords.DimCoord(b, long_name="bearing", units="degrees") - ) - surroundings.append(surround) - surroundings.merge() - print(surroundings) - - curv = surroundings[0].copy() - curv -= cube - - curv.data[curv.data > tol] = -1.0 - curv.data[curv.data < tol] = 1.0 - - curv.collapsed("bearing", iris.analysis.SUM) - curv.rename(f"CURV_calculated_from_{num_radial_points}_radial_points") - curv.units = "1" - curv_cubes.append(curv) - if len(curv_cubes) == 1: - return curv_cubes[0] - else: - return curv_cubes - - -def global_curv_versiontwo(central, radius, num_radial_points=16, tol=0.0): +def global_curv(central, radius, num_radial_points=16, tol=0.0): """Calculate the CURV diagnostic for a global model VERSION 2 for faster calculation speed.""" if (num_radial_points == 8) or (num_radial_points == 16): bearing = np.linspace(0.0, 360.0, num_radial_points + 1)[:-1] diff --git a/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml b/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml index 87e55b2d3..14e09add8 100644 --- a/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml +++ b/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml @@ -19,7 +19,7 @@ steps: subarea_type: $SUBAREA_TYPE subarea_extent: $SUBAREA_EXTENT - - operator: curvature.global_curv_versiontwo + - operator: curvature.global_curv radius: $CURV_RADIUS num_radial_points: $CURV_POINTS From 74c13daf1f9e57afc6ab0796deb0ac67e1ba8412 Mon Sep 17 00:00:00 2001 From: daflack Date: Thu, 14 May 2026 17:01:52 +0100 Subject: [PATCH 13/17] Adds new colorbar for curv --- src/CSET/operators/plot.py | 63 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index c6a44e09b..e026a0aed 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -285,6 +285,9 @@ def _colorbar_map_levels(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = if any("aviation_colour_state" in name for name in varnames): cmap, levels, norm = _custom_colormap_aviation_colour_state(cube) return cmap, levels, norm + if any("CURV_" in name for name in varnames): + cmap, levels, norm = _custom_colormap_curv(cube) + return cmap, levels, norm # If no valid colormap has been defined, use defaults and return. if not cmap: @@ -2293,6 +2296,66 @@ def _custom_colormap_aviation_colour_state(cube: iris.cube.Cube): return cmap, levels, norm +def _custom_colormap_curv(cube: iris.cube.Cube): + """Return custom colourmap for curv. + + If "CURV_" appears anywhere in the name of a cube + this function will be called. + + 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. + """ + if "8" in cube.long_name(): + levels = [-17, -15, -13, -11, -9, -7, -5, -3, -1, 1, 3, 5, 7, 9, 11, 13, 15, 17] + colors = [ + "#ff0000", + "#a52a2a", + "#dc143c", + "#ff4500", + "#f97306", + "#ffa500", + "#fac205", + "#ffd700", + "#ffffff", + "#7fffd4", + "#00ffff", + "#069af3", + "#0323df", + "#0000ff", + "#00008b", + "#030764", + "#01153e", + ] + else: + levels = [-9, -7, -5, -3, -1, 1, 3, 5, 7, 9] + colors = [ + "#ff0000", + "#dc143c", + "#f97306", + "#fac205", + "#ffffff", + "#00ffff", + "#0323df", + "#00008b", + "#01153e", + ] + # Create a custom colormap + cmap = mcolors.ListedColormap(colors) + # Normalise the levels + norm = mcolors.BoundaryNorm(levels, cmap.N) + return cmap, levels, norm + + def _custom_colourmap_visibility_in_air(cube: iris.cube.Cube, cmap, levels, norm): """Return a custom colourmap for the current recipe.""" varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name]) From d3482a32946e2b95a3b6827b1d80b4bc195ff1cf Mon Sep 17 00:00:00 2001 From: daflack Date: Thu, 14 May 2026 17:07:26 +0100 Subject: [PATCH 14/17] Updates colorbar --- 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 e026a0aed..d7dd5aa91 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -2315,7 +2315,7 @@ def _custom_colormap_curv(cube: iris.cube.Cube): should be taken as the range. norm: BoundaryNorm. """ - if "8" in cube.long_name(): + if "8" in cube.long_name: levels = [-17, -15, -13, -11, -9, -7, -5, -3, -1, 1, 3, 5, 7, 9, 11, 13, 15, 17] colors = [ "#ff0000", From 0cc4925a38fdbbfae02f25b4893c18471c8521c0 Mon Sep 17 00:00:00 2001 From: daflack Date: Thu, 14 May 2026 17:14:06 +0100 Subject: [PATCH 15/17] Updates colorbar --- src/CSET/operators/plot.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index d7dd5aa91..3b7ec7402 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -2315,7 +2315,7 @@ def _custom_colormap_curv(cube: iris.cube.Cube): should be taken as the range. norm: BoundaryNorm. """ - if "8" in cube.long_name: + if "16" in cube.long_name: levels = [-17, -15, -13, -11, -9, -7, -5, -3, -1, 1, 3, 5, 7, 9, 11, 13, 15, 17] colors = [ "#ff0000", @@ -2349,10 +2349,10 @@ def _custom_colormap_curv(cube: iris.cube.Cube): "#00008b", "#01153e", ] - # Create a custom colormap - cmap = mcolors.ListedColormap(colors) - # Normalise the levels - norm = mcolors.BoundaryNorm(levels, cmap.N) + # Create a custom colormap + cmap = mcolors.ListedColormap(colors) + # Normalise the levels + norm = mcolors.BoundaryNorm(levels, cmap.N) return cmap, levels, norm From 22aa95e5dda40359ed3763bc3c53762150fbf265 Mon Sep 17 00:00:00 2001 From: daflack Date: Fri, 15 May 2026 15:39:56 +0100 Subject: [PATCH 16/17] Update colorbar and default tolerance --- src/CSET/operators/curvature.py | 6 ++--- src/CSET/operators/plot.py | 46 ++++++++++++++++----------------- 2 files changed, 26 insertions(+), 26 deletions(-) diff --git a/src/CSET/operators/curvature.py b/src/CSET/operators/curvature.py index c29b75c27..834fd6014 100644 --- a/src/CSET/operators/curvature.py +++ b/src/CSET/operators/curvature.py @@ -44,12 +44,12 @@ def _lat_lon_identifier(original_lat, original_lon, distance, bearing): new_lat = np.degrees(new_lat) new_lon = np.degrees(new_lon) - new_lon = np.where(new_lon > 360.0, new_lon - 360.0, new_lon) + new_lon = np.where(new_lon > 180.0, new_lon - 360.0, new_lon) return new_lat, new_lon -def global_curv(central, radius, num_radial_points=16, tol=0.0): +def global_curv(central, radius, num_radial_points=16, tol=400.0): """Calculate the CURV diagnostic for a global model VERSION 2 for faster calculation speed.""" if (num_radial_points == 8) or (num_radial_points == 16): bearing = np.linspace(0.0, 360.0, num_radial_points + 1)[:-1] @@ -87,7 +87,7 @@ def global_curv(central, radius, num_radial_points=16, tol=0.0): curv_data = surroundings - cube.data curv_data_binary = np.zeros(curv_data.shape) curv_data_binary[curv_data > tol] = -1.0 - curv_data_binary[curv_data < tol] = 1.0 + curv_data_binary[curv_data < -tol] = 1.0 curv_data_collapsed = np.sum(curv_data_binary, 0) curv.data = curv_data_collapsed curv.rename(f"CURV_calculated_from_{num_radial_points}_radial_points") diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index 3b7ec7402..dd55f62b8 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -2318,36 +2318,36 @@ def _custom_colormap_curv(cube: iris.cube.Cube): if "16" in cube.long_name: levels = [-17, -15, -13, -11, -9, -7, -5, -3, -1, 1, 3, 5, 7, 9, 11, 13, 15, 17] colors = [ + "#01153e", + "#030764", + "#00008b", + "#0000ff", + "#0323df", + "#069af3", + "#00ffff", + "#7fffd4", + "#ffffff", + "#ffd700", + "#fac205", + "#ffa500", + "#f97306", + "#ff4500", "#ff0000", - "#a52a2a", "#dc143c", - "#ff4500", - "#f97306", - "#ffa500", - "#fac205", - "#ffd700", - "#ffffff", - "#7fffd4", - "#00ffff", - "#069af3", - "#0323df", - "#0000ff", - "#00008b", - "#030764", - "#01153e", + "#a52a2a", ] else: levels = [-9, -7, -5, -3, -1, 1, 3, 5, 7, 9] colors = [ - "#ff0000", - "#dc143c", - "#f97306", - "#fac205", - "#ffffff", - "#00ffff", - "#0323df", - "#00008b", "#01153e", + "#00008b", + "#0323df", + "#00ffff", + "#ffffff", + "#fac205", + "#f97306", + "#ff0000", + "#a52a2a", ] # Create a custom colormap cmap = mcolors.ListedColormap(colors) From a4f633d794a7df690c40753e2543c1b3e4fbf9cb Mon Sep 17 00:00:00 2001 From: daflack Date: Fri, 15 May 2026 16:30:37 +0100 Subject: [PATCH 17/17] Adds regridding to recipe --- .../daily_weather/curv_spatial_plot.yaml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml b/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml index 14e09add8..846fc575d 100644 --- a/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml +++ b/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml @@ -19,9 +19,15 @@ steps: subarea_type: $SUBAREA_TYPE subarea_extent: $SUBAREA_EXTENT + - operator: regrid.regrid_onto_xyspacing + xspacing: 1.0 + yspacing: 1.0 + method: Linear + - operator: curvature.global_curv radius: $CURV_RADIUS num_radial_points: $CURV_POINTS + tol: 400.0 - operator: plot.spatial_pcolormesh_plot