Adds a global model application for CURV

docs/source/reference/operators.rst

@@ -92,6 +92,15 @@ CSET.operators.convection
 .. automodule:: CSET.operators.convection
    :members:
 
+Curvature Operators
+~~~~~~~~~~~~~~~~~~~
+
+CSET.operators.curvature
+------------------------
+
+.. automodule:: CSET.operators.curvature
+   :members:
+
 Ensemble Operators
 ~~~~~~~~~~~~~~~~~~
 

src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf

@@ -1327,6 +1327,43 @@ 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-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 the 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 the CURV calculation.
+type=python_list
+compulsory=true
+sort-key=dw-curv2
+
+[template variables=CURV_TOLERANCE]
+ns=Diagnostics/Derived/Daily
+title=Tolerance to use in CURV
+description=The tolerance in Pa to use in the CURV calculation.
+            A recommended tolerance to use is 200 Pa (2 hPa).
+type=python_list
+compulsory=true
+sort-key=dw-curv3
 ###################################
 # Aviation
 [Diagnostics/Aviation]

src/CSET/cset_workflow/rose-suite.conf.example

@@ -44,6 +44,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

src/CSET/loaders/spatial_field.py

@@ -704,6 +704,26 @@ def load(conf: Config):
                 aggregation=False,
             )
 
+    if conf.GLOBAL_CURV:
+        for model, points, radius, tol in itertools.product(
+            models, conf.CURV_POINTS, conf.CURV_RADIUS, conf.CURV_TOLERANCE
+        ):
+            yield RawRecipe(
+                recipe="curv_spatial_plot.yaml",
+                variables={
+                    "MODEL_NAME": model["name"],
+                    "CURV_POINTS": points,
+                    "CURV_RADIUS": radius,
+                    "CURV_TOLERANCE": tol,
+                    "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,
+            )
+
     # Directional shear across Moist Absolutely Unstable Layer (of deepest)
     if conf.DIRECTIONAL_SHEAR_ACROSS_MAUL:
         for model in models:

src/CSET/operators/__init__.py

@@ -32,6 +32,7 @@
     collapse,
     constraints,
     convection,
+    curvature,
     ensembles,
     feature,
     filters,
@@ -61,6 +62,7 @@
     "collapse",
     "constraints",
     "convection",
+    "curvature",
     "ensembles",
     "execute_recipe",
     "feature",

src/CSET/operators/_colormaps.py

@@ -207,6 +207,9 @@ def colorbar_map_levels(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = N
     if any("feature_" in name for name in varnames):
         cmap, levels, norm = custom_colormap_feature_tracking(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:
@@ -727,3 +730,63 @@ def custom_colormap_feature_tracking(cube: iris.cube.Cube):
     cmap = cmap.with_extremes(under="white")
 
     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 "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",
+            "#dc143c",
+            "#a52a2a",
+        ]
+    else:
+        levels = [-9, -7, -5, -3, -1, 1, 3, 5, 7, 9]
+        colors = [
+            "#01153e",
+            "#00008b",
+            "#0323df",
+            "#00ffff",
+            "#ffffff",
+            "#fac205",
+            "#f97306",
+            "#ff0000",
+            "#a52a2a",
+        ]
+    # Create a custom colormap
+    cmap = mcolors.ListedColormap(colors)
+    # Normalise the levels
+    norm = mcolors.BoundaryNorm(levels, cmap.N)
+    return cmap, levels, norm

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 logging
+
+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."""
+    R = 6371.0  # km.
+    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_lat = np.degrees(new_lat)
+    new_lon = np.degrees(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=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]
+    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):
+        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.Linear(),
+            )
+            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

src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml

@@ -0,0 +1,31 @@
+category: Daily Weather
+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 spatial plot.
+  It uses $CURV_POINTS points and a radius of $CURV_RADIUS km and a
+  tolerance of $CURV_TOLERANCE Pa.
+  +$CURV_POINTS is a full anticyclone and -$CURV_POINTS is a full cyclone.
+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_sea_level
+        cell_methods_constraint:
+            operator: constraints.generate_cell_methods_constraint
+            cell_methods: []
+    subarea_type: $SUBAREA_TYPE
+    subarea_extent: $SUBAREA_EXTENT
+
+  - operator: curvature.global_curv
+    radius: $CURV_RADIUS
+    num_radial_points: $CURV_POINTS
+    tol: $CURV_TOLERANCE
+
+  - operator: plot.spatial_pcolormesh_plot
+
+  - operator: write.write_cube_to_nc
+    overwrite: True