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/cset_workflow/meta/diagnostics/rose-meta.conf b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf index e3fbaace3..1c53be724 100644 --- a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf @@ -1327,6 +1327,34 @@ 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 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/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 diff --git a/src/CSET/loaders/spatial_field.py b/src/CSET/loaders/spatial_field.py index ba35ca1e2..aa91b79ab 100644 --- a/src/CSET/loaders/spatial_field.py +++ b/src/CSET/loaders/spatial_field.py @@ -572,6 +572,25 @@ def load(conf: Config): aggregation=False, ) + if conf.GLOBAL_CURV: + 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 + 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/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", diff --git a/src/CSET/operators/curvature.py b/src/CSET/operators/curvature.py new file mode 100644 index 000000000..834fd6014 --- /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 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 diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index c6a44e09b..dd55f62b8 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 "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 + + 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]) 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..846fc575d --- /dev/null +++ b/src/CSET/recipes/derived_diagnostics/daily_weather/curv_spatial_plot.yaml @@ -0,0 +1,35 @@ +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. + +$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_mean_sea_level + cell_methods_constraint: + operator: constraints.generate_cell_methods_constraint + cell_methods: [] + 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 + + - operator: write.write_cube_to_nc + overwrite: True