From 9d4b574458a72eb7b4bb8c16e3efebc633515026 Mon Sep 17 00:00:00 2001 From: "John M. Edwards" Date: Thu, 31 Jul 2025 09:30:09 +0100 Subject: [PATCH 1/7] Replicating changes in a new branch. --- src/CSET/cset_workflow/.gitignore | 1 + .../app/fetch_fcst/bin/fetch_data.py | 91 ++++++++++++++ .../cset_workflow/app/fetch_obs/rose-app.conf | 3 + src/CSET/cset_workflow/flow.cylc | 28 +++++ .../spatial_surface_model_obs_diference.cylc | 26 ++++ .../includes/spatial_surface_obs_field.cylc | 18 +++ .../meta/observations/rose-meta.conf | 116 ++++++++++++++++++ src/CSET/cset_workflow/meta/rose-meta.conf | 1 + .../cset_workflow/meta/rose-meta.conf.jinja2 | 11 ++ .../cset_workflow/rose-suite.conf.example | 8 ++ src/CSET/operators/plot.py | 85 ++++++++++++- src/CSET/operators/read.py | 11 ++ src/CSET/operators/regrid.py | 68 ++++++++++ ...eric_model_obs_difference_scatterplot.yaml | 44 +++++++ src/CSET/recipes/generic_obs_scatterplot.yaml | 20 +++ 15 files changed, 529 insertions(+), 2 deletions(-) create mode 100644 src/CSET/cset_workflow/app/fetch_obs/rose-app.conf create mode 100644 src/CSET/cset_workflow/includes/spatial_surface_model_obs_diference.cylc create mode 100644 src/CSET/cset_workflow/includes/spatial_surface_obs_field.cylc create mode 100644 src/CSET/cset_workflow/meta/observations/rose-meta.conf create mode 100644 src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml create mode 100644 src/CSET/recipes/generic_obs_scatterplot.yaml diff --git a/src/CSET/cset_workflow/.gitignore b/src/CSET/cset_workflow/.gitignore index 6932c7da3..a49976707 100644 --- a/src/CSET/cset_workflow/.gitignore +++ b/src/CSET/cset_workflow/.gitignore @@ -2,6 +2,7 @@ site/ opt/ demo_pointstat/ +*metdb* restricted* # User workflow configuration. diff --git a/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py b/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py index 39b4c92f2..dbdd168c5 100644 --- a/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py +++ b/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py @@ -15,6 +15,7 @@ """Retrieve the files from the filesystem for the current cycle point.""" import abc +import ast import glob import itertools import logging @@ -176,6 +177,25 @@ def _get_needed_environment_variables() -> dict: logging.debug("Environment variables loaded: %s", variables) return variables +def _get_needed_environment_variables_obs() -> dict: + """Load the needed variables from the environment.""" + variables = { + "subtype": os.environ["OBS_SUBTYPE"], + "data_time": datetime.fromisoformat(os.environ["CYLC_TASK_CYCLE_POINT"]), + "forecast_length": isodate.parse_duration(os.environ["ANALYSIS_LENGTH"]), + "obs_fields": ast.literal_eval(os.environ["SURFACE_SYNOP_FIELDS"]), + "model_identifier": "OBS", + "wmo_nmbrs": ast.literal_eval(os.environ.get("WMO_BLOCK_STTN_NMBRS")) \ + if len(os.environ.get("WMO_BLOCK_STTN_NMBRS")) > 0 else None, + "subarea_extent":ast.literal_eval(os.environ.get("SUBAREA_EXTENT")) \ + if len(os.environ.get("SUBAREA_EXTENT")) > 0 else None, + "obs_interval": isodate.parse_duration(os.environ["SURFACE_SYNOP_INTERVAL"]), + "obs_offset": isodate.parse_duration(os.environ["SURFACE_SYNOP_OFFSET"]), + "rose_datac": os.environ["ROSE_DATAC"], + } + logging.debug("Environment variables loaded: %s", variables) + return variables + def _template_file_path( raw_path: str, @@ -270,3 +290,74 @@ def fetch_data(file_retriever: FileRetrieverABC): # exiting the with block. if not files_found: raise FileNotFoundError("No files found for model!") + + +def fetch_obs(obs_retriever: FileRetrieverABC = FilesystemFileRetriever): + """Fetch the observations corresponding to a model run. + + The following environment variables need to be set: + * ANALYSIS_OFFSET + * ANALYSIS_LENGTH + * CYLC_TASK_CYCLE_POINT + * DATA_PATH + * DATA_PERIOD + * DATE_TYPE + * MODEL_IDENTIFIER + * ROSE_DATAC + + Parameters + ---------- + obs_retriever: ObsRetriever + ObsRetriever implementation to use. Defaults to FilesystemFileRetriever. + + Raises + ------ + FileNotFound: + If no observations are available. + """ + v = _get_needed_environment_variables_obs() + + # Prepare output directory. + cycle_obs_dir = f"{v['rose_datac']}/data/OBS" + os.makedirs(cycle_obs_dir, exist_ok=True) + logging.debug("Output directory: %s", cycle_obs_dir) + + # We will get just one file for now, but follow the templating + # syntax for the model for consistency. + obs_base_path = ( + v["subtype"] + + "_" + + "%Y%m%dT%H%MZ_dt_" + + str(int(v["forecast_length"].total_seconds() // 3600)).zfill(3) + + ".nc" + ) + paths = _template_file_path( + obs_base_path, + "initiation", + v["data_time"], + v["forecast_length"], + timedelta(seconds=0), + v["obs_interval"], + ) + logging.info("Retrieving paths:\n%s", "\n".join(paths)) + + # Use obs retriever to transfer data with multiple threads. + # We shouldn't need to iterate as we do for the forecast data + # because these files will be smaller. Passing None is a temporary + # measure until we can work out why python doesn't think this is a + # class method. + try: + obs_retriever.get_file( + paths[0], + v["subtype"], + v["obs_fields"], + v["data_time"], + v["obs_offset"], + v["forecast_length"], + v["obs_interval"], + cycle_obs_dir, + wmo_nmbrs=v["wmo_nmbrs"], + subarea_extent=v["subarea_extent"], + ) + except: + raise FileNotFoundError("No observations available.") diff --git a/src/CSET/cset_workflow/app/fetch_obs/rose-app.conf b/src/CSET/cset_workflow/app/fetch_obs/rose-app.conf new file mode 100644 index 000000000..ec3d1e950 --- /dev/null +++ b/src/CSET/cset_workflow/app/fetch_obs/rose-app.conf @@ -0,0 +1,3 @@ +[command] +default=echo "Please set ROSE_APP_COMMAND_KEY to your database."; false +metdb=app_env_wrapper fetch-obs-metdb.py diff --git a/src/CSET/cset_workflow/flow.cylc b/src/CSET/cset_workflow/flow.cylc index 9fad81ad1..835ec33d2 100644 --- a/src/CSET/cset_workflow/flow.cylc +++ b/src/CSET/cset_workflow/flow.cylc @@ -45,6 +45,7 @@ URL = https://metoffice.github.io/CSET R1/{{date}} = """ setup_complete[^] => FETCH_DATA:succeed-all => + FETCH_OBS:succeed-all => fetch_complete => PROCESS:{% if LOGLEVEL == "DEBUG" %}succeed-all{% else %}finish-all{% endif %} => cycle_complete @@ -55,6 +56,7 @@ URL = https://metoffice.github.io/CSET {{CSET_TRIAL_CYCLE_PERIOD}} = """ setup_complete[^] => FETCH_DATA:succeed-all => + FETCH_OBS:succeed-all => fetch_complete => PROCESS:{% if LOGLEVEL == "DEBUG" %}succeed-all{% else %}finish-all{% endif %} => cycle_complete @@ -112,6 +114,8 @@ URL = https://metoffice.github.io/CSET execution retry delays = PT1M, PT15M, PT1H [[[directives]]] --mem=10000 + [[[environment]]] + PLOTTING_PROJECTION = {{PLOTTING_PROJECTION | default("")}} [[PROCESS_CASE_AGGREGATION]] script = rose task-run -v --app-key=run_cset_recipe @@ -128,6 +132,12 @@ URL = https://metoffice.github.io/CSET [[[environment]]] ANALYSIS_LENGTH = {{ANALYSIS_LENGTH}} + [[FETCH_OBS]] + script = rose task-run -v --app-key=fetch_obs + execution time limit = PT1H + [[[environment]]] + ANALYSIS_LENGTH = {{ANALYSIS_LENGTH}} + [[METPLUS]] [[[environment]]] {% if METPLUS_GRID_STAT|default(False) %} @@ -153,6 +163,9 @@ URL = https://metoffice.github.io/CSET [[cycle_complete]] inherit = DUMMY_TASK + [[dummy_fetch_obs]] + inherit = DUMMY_TASK, FETCH_OBS + [[dummy_process]] inherit = DUMMY_TASK, PROCESS @@ -188,6 +201,21 @@ URL = https://metoffice.github.io/CSET ANALYSIS_OFFSET = {{model["analysis_offset"]}} {% endfor %} + {% if SURFACE_SYNOP_OBS %} + [[fetch_obs]] + # Fetch observations from the MetDB. + inherit = FETCH_OBS + [[[environment]]] + MODEL_IDENTIFIER = "OBS" + ROSE_APP_COMMAND_KEY = "metdb" + SURFACE_SYNOP_FIELDS = {{SURFACE_SYNOP_FIELDS}} + SURFACE_SYNOP_INTERVAL = {{SURFACE_SYNOP_INTERVAL}} + SURFACE_SYNOP_OFFSET = {{SURFACE_SYNOP_OFFSET}} + WMO_BLOCK_STTN_NMBRS = {{WMO_BLOCK_STTN_NMBRS | default("")}} + SUBAREA_EXTENT = {{SUBAREA_EXTENT | default("")}} + OBS_SUBTYPE = "LNDSYB" + {% endif %} + [[housekeeping]] # Housekeep unprocessed data files and processed intermediate files. [[[environment]]] diff --git a/src/CSET/cset_workflow/includes/spatial_surface_model_obs_diference.cylc b/src/CSET/cset_workflow/includes/spatial_surface_model_obs_diference.cylc new file mode 100644 index 000000000..2b5b92be0 --- /dev/null +++ b/src/CSET/cset_workflow/includes/spatial_surface_model_obs_diference.cylc @@ -0,0 +1,26 @@ +{% if SURFACE_SYNOP_OBS and SURFACE_SYNOP_DIFFS %} +{% for model in models %} +{% for obs_field in (SURFACE_SYNOP_FIELDS | list) %} +{# {% set model_field = (SURFACE_SYNOP_MODEL_FIELDS | list)[SURFACE_SYNOP_FIELDS.index(obs_fields)] %} #} +[runtime] + [[obs_model_difference_scatterplot_m{{model["id"]}}_{{obs_field}}]] + inherit = PROCESS + execution time limit = PT6H + [[[environment]]] + CSET_RECIPE_NAME = "generic_model_obs_difference_scatterplot.yaml" + CSET_ADDOPTS = """ + --OBSVARNAME={{obs_field}} + --VARNAME={{(SURFACE_SYNOP_MODEL_FIELDS | list)[loop.index0]}} + --BASE_MODEL='{{model["id"]}}' + --OBSERVATIONS='OBS' + --MODEL_NAME='{{model["name"]}}' + --DIFFERENCE=True + --SUBAREA_TYPE={{SUBAREA_TYPE | default("")}} + --SUBAREA_EXTENT={{SUBAREA_TYPE | default("")}} + --PLOTTING_PROJECTION={{PLOTTING_PROJECTION | default("")}} + """ + MODEL_IDENTIFIERS = 'OBS {{model["id"]}}' + +{% endfor %} +{% endfor %} +{% endif %} diff --git a/src/CSET/cset_workflow/includes/spatial_surface_obs_field.cylc b/src/CSET/cset_workflow/includes/spatial_surface_obs_field.cylc new file mode 100644 index 000000000..f0ce429d9 --- /dev/null +++ b/src/CSET/cset_workflow/includes/spatial_surface_obs_field.cylc @@ -0,0 +1,18 @@ +{% if SURFACE_SYNOP_OBS %} +{% for obs_field in SURFACE_SYNOP_FIELDS %} +[runtime] + [[generic_obs_scatterplot_{{obs_field}}]] + inherit = PROCESS + execution time limit = PT1H + [[[environment]]] + CSET_RECIPE_NAME = "generic_obs_scatterplot.yaml" + CSET_ADDOPTS = """ + --OBSVARNAME={{obs_field}} + --SUBAREA_TYPE={{SUBAREA_TYPE | default("")}} + --SUBAREA_EXTENT={{SUBAREA_TYPE | default("")}} + --PLOTTING_PROJECTION={{PLOTTING_PROJECTION | default("")}} + """ + MODEL_IDENTIFIERS = "OBS" + +{% endfor %} +{% endif %} diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf new file mode 100644 index 000000000..b25be25d4 --- /dev/null +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -0,0 +1,116 @@ +################################################################################ +# Observations +################################################################################ + +[Observations] +ns=Observations +sort-key=sec-d +title=Observational Comparisons + +################################################################################ +# Spatial (2D map) fields +################################################################################ + +[Observations/Types] +ns=Observations/Types +sort-key=sec-d1 +title=Types of observations + +# Surface synoptic fields. +[template variables=SURFACE_SYNOP_OBS] +ns=Observations/Types +description=Compare with surface synoptic observations +type=python_boolean +compulsory=true +sort-key=0synop1a +trigger=template variables=SURFACE_SYNOP_FIELDS: True; + template variables=USE_WMO_STATION_NUMBERS: True; + template variables=SURFACE_SYNOP_INTERVAL: True; + template variables=SURFACE_SYNOP_OFFSET: True; + +[template variables=PLACEHOLDER_OBS] +ns=Observations/Types +description=Compare with observations from an undefined source. +help=This variable exists as a placeholder to show the intended + structure the GUI. Additional sources of observations should be + added using one window for each source. + +type=python_boolean +compulsory=true +sort-key=0yopp1a + +################################################################ + +[Observations/Types/Synoptic] +ns=Observations/Types/Synoptic +sort-key=sec-d2 +title=Synoptic observations + +[template variables=SURFACE_SYNOP_FIELDS] +ns=Observations/Types/Synoptic +title=Surface synoptic observations +description=Synoptic observations +help=Name of field in the style of the MetDB. +compulsory=true +type=python_list +sort-key=0synop1b + +[template variables=SURFACE_SYNOP_DIFFS] +ns=Observations/Types/Synoptic +description=Calculate the corresponding differences. +type=python_boolean +compulsory=true +sort-key=0synop1c +trigger=template variables=SURFACE_SYNOP_MODEL_FIELDS: True; + +[template variables=SURFACE_SYNOP_MODEL_FIELDS] +ns=Observations/Types/Synoptic +title=Corresponding model variables +description=Synoptic observations +help=Name of field in the model: the order must match that of the synoptic observations +compulsory=true +type=python_list +sort-key=0synop1d + +[template variables=USE_WMO_STATION_NUMBERS] +ns=Observations/Types/Synoptic +description=Use WMO numbers +help=Select to set stations by WMO station or block number. +compulsory=true +type=python_boolean +sort-key=0synop1e +trigger=template variables=WMO_BLOCK_STTN_NMBRS: True; + +[template variables=WMO_BLOCK_STTN_NMBRS] +ns=Observations/Types/Synoptic +title=WMO Block or station numbers. +description=WMO numbers +help=WMO block numbers (2 digits) or station numbers (5 digits). +compulsory=true +type=python_list +sort-key=0synop1f + +[template variables=SURFACE_SYNOP_INTERVAL] +ns=Observations/Types/Synoptic +title=Interval between synoptic observations +description=Synoptic observations +help=The time interval between rerievals of synoptic observations. +compulsory=true +type=quoted +sort-key=0synop1g + +[template variables=SURFACE_SYNOP_OFFSET] +ns=Observations/Types/Synoptic +title=Offset for synoptic observations +description=Offset for observations +help=The offset of the first retrieved observations relative to the initial time. +compulsory=true +type=quoted +sort-key=0synop1h + +################################################################ + +[Observations/Types/Placeholder] +ns=Observations/Types/Placeholder +sort-key=sec-e2 +title=Observations from an undefined source. diff --git a/src/CSET/cset_workflow/meta/rose-meta.conf b/src/CSET/cset_workflow/meta/rose-meta.conf index 99a67b25a..d73c5f324 100644 --- a/src/CSET/cset_workflow/meta/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/rose-meta.conf @@ -1,5 +1,6 @@ # Diagnostics settings are split into a separate file. import=meta/diagnostics + =meta/observations [ns=Setup] sort-key=sec-a diff --git a/src/CSET/cset_workflow/meta/rose-meta.conf.jinja2 b/src/CSET/cset_workflow/meta/rose-meta.conf.jinja2 index a2ec063ac..ec33e2e9f 100644 --- a/src/CSET/cset_workflow/meta/rose-meta.conf.jinja2 +++ b/src/CSET/cset_workflow/meta/rose-meta.conf.jinja2 @@ -158,6 +158,17 @@ type=quoted compulsory=true sort-key=setup-h-out2 +[template variables=PLOTTING_PROJECTION] +ns=Setup +title=Projection for the plots +description=Geographical projection of plots +help=This is is the projection used for the final plots, specified as a string. + If no value is given, native latitudes and longitudes are used. +compulsory=true +values="Unspecified", "NP_Stereo" +value-titles=Unspecified, North polar stereographic +sort-key=setup-h-out3a + [template variables=PLOT_RESOLUTION] ns=Setup title=Plot resolution diff --git a/src/CSET/cset_workflow/rose-suite.conf.example b/src/CSET/cset_workflow/rose-suite.conf.example index 2c6b93a72..a7a419cec 100644 --- a/src/CSET/cset_workflow/rose-suite.conf.example +++ b/src/CSET/cset_workflow/rose-suite.conf.example @@ -53,9 +53,11 @@ MODEL_LEVEL_FIELDS=[] !!MODULES_LIST= !!MODULES_PURGE=True !!ONE_TO_ONE=False +PLACEHOLDER_OBS=False !!PLEVEL_TRANSECT_FINISHCOORDS= !!PLEVEL_TRANSECT_STARTCOORDS= PLOT_RESOLUTION=100 +!!PLOTTING_PROJECTION="" PRESSURE_LEVELS=[] PRESSURE_LEVEL_FIELDS=[] PROFILE_MLEVEL=False @@ -85,6 +87,11 @@ SPATIAL_SURFACE_FIELD_METHOD=[""] !!SUBAREA_EXTENT=0,0,0,0 !!SUBAREA_TYPE="realworld" SURFACE_FIELDS=[] +SURFACE_SYNOP_OBS=False +!!SURFACE_SYNOP_FIELDS=[] +!!WMO_BLOCK_STTN_NMBRS=[] +!!SURFACE_SYNOP_INTERVAL="" +!!SURFACE_SYNOP_OFFSET="" SURFACE_SINGLE_POINT_TIME_SERIES=False TIMESERIES_MLEVEL_FIELD=False TIMESERIES_MLEVEL_FIELD_AGGREGATION=False,False,False,False @@ -92,6 +99,7 @@ TIMESERIES_PLEVEL_FIELD=False TIMESERIES_PLEVEL_FIELD_AGGREGATION=False,False,False,False TIMESERIES_SURFACE_FIELD=False TIMESERIES_SURFACE_FIELD_AGGREGATION=False,False,False,False +!!USE_WMO_STATION_NUMBERS=False !!VERTICAL_COORDINATE_A=[] !!VERTICAL_COORDINATE_B=[] WEB_ADDR="" diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index b3a6db8b3..eeb604cb9 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -25,6 +25,7 @@ import sys from typing import Literal +import cartopy.crs as ccrs import iris import iris.coords import iris.cube @@ -1183,12 +1184,79 @@ def _plot_and_save_postage_stamps_in_single_plot_histogram_series( plt.close(fig) +def _plot_and_save_scattermap_plot( + cube: iris.cube.Cube, filename: str, title: str, projection=None: str, **kwargs +): + """Plot and save a geographical scatter plot. + + Parameters + ---------- + cube: Cube + 1 dimensional Cube of the data points with auxiliary latitude and + longitude coordinates, + filename: str + Filename of the plot to write. + title: str + Plot title. + projection: str + Mapping projection to be used by cartopy. + """ + # Setup plot details, size, resolution, etc. + fig = plt.figure(figsize=(15, 15), facecolor="w", edgecolor="k") + if (projection is not None): + # Apart from the default, the only projection we currently support is + # a stereographic projection over the North Pole. + if (projection == 'NP_Stereo'): + axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0)) + else: + Raise ValueError(f"Unknown projection: {projection}") + else: + axes = plt.axes(projection=ccrs.PlateCarree()) + + # Filled contour plot of the field. + mrk_size = int(np.sqrt(5000000.0 / len(cube.data))) + klon = None + klat = None + for kc in range(len(cube.aux_coords)): + if cube.aux_coords[kc].standard_name == "latitude": + klat = kc + elif cube.aux_coords[kc].standard_name == "longitude": + klon = kc + scatter_map = iplt.scatter( + cube.aux_coords[klon], + cube.aux_coords[klat], + c=cube.data[:], + s=mrk_size, + cmap="jet", + edgecolors="k", + ) + + # Add coastlines. + try: + axes.coastlines() + except AttributeError: + pass + + # Add title. + axes.set_title(title, fontsize=16) + + # Add colour bar. + cbar = fig.colorbar(scatter_map) + cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20) + + # Save plot. + fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution()) + logging.info("Saved geographical scatter plot to %s", filename) + plt.close(fig) + + def _spatial_plot( method: Literal["contourf", "pcolormesh"], cube: iris.cube.Cube, filename: str | None, sequence_coordinate: str, stamp_coordinate: str, + **kwargs, ): """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube. @@ -1239,6 +1307,14 @@ def _spatial_plot( except iris.exceptions.CoordinateNotFoundError: pass + # Produce a geographical scatter plot if the data have a + # dimension called observation or model_obs_error + if any( + crd.var_name == "station" or crd.var_name == "model_obs_error" + for crd in cube.coords() + ): + plotting_func = _plot_and_save_scattermap_plot + # Must have a sequence coordinate. try: cube.coord(sequence_coordinate) @@ -1266,6 +1342,7 @@ def _spatial_plot( stamp_coordinate=stamp_coordinate, title=title, method=method, + **kwargs, ) plot_index.append(plot_filename) @@ -1439,7 +1516,9 @@ def spatial_contour_plot( TypeError If the cube isn't a single cube. """ - _spatial_plot("contourf", cube, filename, sequence_coordinate, stamp_coordinate) + _spatial_plot( + "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs + ) return cube @@ -1488,7 +1567,9 @@ def spatial_pcolormesh_plot( TypeError If the cube isn't a single cube. """ - _spatial_plot("pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate) + _spatial_plot( + "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs + ) return cube diff --git a/src/CSET/operators/read.py b/src/CSET/operators/read.py index 2f5bb7002..1a21dbc58 100644 --- a/src/CSET/operators/read.py +++ b/src/CSET/operators/read.py @@ -510,6 +510,17 @@ def _lfric_time_coord_fix_callback(cube: iris.cube.Cube, field, filename): not isinstance(time_coord, iris.coords.DimCoord) and len(cube.coord_dims(time_coord)) == 1 ): + # Fudge the bounds to foil checking for strict monotonicity. + if time_coord.has_bounds(): + if (time_coord.bounds[-1][0] - time_coord.bounds[0][0]) < 1.0e-8: + time_coord.bounds = [ + [ + time_coord.bounds[i][0] + 1.0e-8 * float(i), + time_coord.bounds[i][1], + ] + for i in range(len(time_coord.bounds)) + ] + iris.util.promote_aux_coord_to_dim_coord(cube, time_coord) # Force single-valued coordinates to be scalar coordinates. diff --git a/src/CSET/operators/regrid.py b/src/CSET/operators/regrid.py index c24cac73d..9722122e4 100644 --- a/src/CSET/operators/regrid.py +++ b/src/CSET/operators/regrid.py @@ -364,3 +364,71 @@ def transform_lat_long_points(lon, lat, cube): lat_rot = rot_coords[1] return lon_rot, lat_rot + +def interpolate_to_point_cube(fld: iris.cube.Cube | iris.cube.CubeList, + point_cube: iris.cube.Cube, + **kwargs) -> iris.cube.Cube | iris.cube.CubeList: + # + # As a basis, create a copy of the point cube. + fld_point_cube = point_cube.copy() + # Get indices of positional coordinates. We assume that the + # point cube is unrotated. + klon = None + klat = None + for kc in range(len(fld_point_cube.aux_coords)): + if fld_point_cube.aux_coords[kc].standard_name == "latitude": + klat = kc + elif fld_point_cube.aux_coords[kc].standard_name == "longitude": + klon = kc + # + # The input may have a rotated coordinate system. + if len(fld.coords("grid_latitude")) > 0: + # Interpolate in rotated coordinates. + rot_csyst = fld.coords("grid_latitude")[0].coord_system + rotpt = iris.analysis.cartography.rotate_pole( + fld_point_cube.aux_coords[klon].points, + fld_point_cube.aux_coords[klat].points, + rot_csyst.grid_north_pole_longitude, + rot_csyst.grid_north_pole_latitude, + ) + # Add other interpolation options later. + fld_interpolator = iris.analysis.Linear(extrapolation_mode="mask").interpolator( + fld, ["time", "grid_latitude", "grid_longitude"] + ) + for jt in range(len(fld_point_cube.coords("time")[0].points)): + fld_point_cube.data[jt, :] = np.ma.masked_invalid( + [ + fld_interpolator( + [ + fld_point_cube.coord("time").points[jt], + rotpt[1][k], + rotpt[0][k], + ] + ).data + if ~point_cube.data.mask[jt][k] + else np.nan + for k in range(len(rotpt[0])) + ] + ) + else: + # Add other interpolation options later. + fld_interpolator = iris.analysis.Linear(extrapolation_mode="mask").interpolator( + fld, ["time", "latitude", "longitude"] + ) + for jt in range(len(fld_point_cube.coords("time")[0].points)): + fld_point_cube.data[jt, :] = np.ma.masked_invalid( + [ + fld_interpolator( + [ + fld_point_cube.coords("time")[0].points[jt], + fld_point_cube.coord("latitude").points[k], + fld_point_cube.coord("longitude").points[k], + ] + ).data + if ~point_cube.data.mask[jt][k] + else np.nan + for k in range(fld_point_cube.coord("latitude").points) + ] + ) + # + return fld_point_cube diff --git a/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml b/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml new file mode 100644 index 000000000..580b0b31e --- /dev/null +++ b/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml @@ -0,0 +1,44 @@ +category: Quick Look +title: Difference between modelled field of $VARNAME in model $MODEL_NAME and observation. +description: | + Plot some observations. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS + + - operator: misc.subtraction + minuend: + operator: regrid.interpolate_to_point_cube + fld: + operator: filters.filter_cubes + constraint: + operator: constraints.combine_constraints + var_constraint: + operator: constraints.generate_var_constraint + varname: $VARNAME + cell_method_constraint: + operator: constraints.generate_cell_methods_constraint + cell_methods: [] + # This will require more work. Several fields are stored with muliple + # cell methods, so generally we want the instantaneous one, but in + # other cases, the field may be available only with another cell method. + # In effect, we shall need to match cell methods to diagnostics. + point_cube: + operator: filters.filter_cubes + constraint: + operator: constraints.generate_var_constraint + varname: $OBSVARNAME + subtrahend: + operator: filters.filter_cubes + constraint: + operator: constraints.generate_var_constraint + varname: $OBSVARNAME + + + - operator: plot.spatial_contour_plot + sequence_coordinate: time + projection: $PLOTTING_PROJECTION + + - operator: write.write_cube_to_nc + overwrite: True diff --git a/src/CSET/recipes/generic_obs_scatterplot.yaml b/src/CSET/recipes/generic_obs_scatterplot.yaml new file mode 100644 index 000000000..b2260a382 --- /dev/null +++ b/src/CSET/recipes/generic_obs_scatterplot.yaml @@ -0,0 +1,20 @@ +category: Quick Look +title: Observed values of $OBSVARNAME. +description: | + Plot some observations. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS + constraint: + operator: constraints.combine_constraints + variable_constraint: + operator: constraints.generate_var_constraint + varname: $OBSVARNAME + + - operator: plot.spatial_contour_plot + sequence_coordinate: time + projection: $PLOTTING_PROJECTION + + - operator: write.write_cube_to_nc + overwrite: True From b6cb7e465d77694955881aed3aaa40b6aa517cbb Mon Sep 17 00:00:00 2001 From: "John M. Edwards" Date: Mon, 4 Aug 2025 11:48:28 +0100 Subject: [PATCH 2/7] Committing tested version. --- .../app/fetch_fcst/bin/fetch_data.py | 21 +++++------ src/CSET/operators/plot.py | 10 +++--- src/CSET/operators/read.py | 35 ++++++++++--------- ...eric_model_obs_difference_scatterplot.yaml | 29 +++++++-------- src/CSET/recipes/generic_obs_scatterplot.yaml | 8 ++--- 5 files changed, 54 insertions(+), 49 deletions(-) diff --git a/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py b/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py index dbdd168c5..d46e2934c 100644 --- a/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py +++ b/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py @@ -177,6 +177,7 @@ def _get_needed_environment_variables() -> dict: logging.debug("Environment variables loaded: %s", variables) return variables + def _get_needed_environment_variables_obs() -> dict: """Load the needed variables from the environment.""" variables = { @@ -185,10 +186,12 @@ def _get_needed_environment_variables_obs() -> dict: "forecast_length": isodate.parse_duration(os.environ["ANALYSIS_LENGTH"]), "obs_fields": ast.literal_eval(os.environ["SURFACE_SYNOP_FIELDS"]), "model_identifier": "OBS", - "wmo_nmbrs": ast.literal_eval(os.environ.get("WMO_BLOCK_STTN_NMBRS")) \ - if len(os.environ.get("WMO_BLOCK_STTN_NMBRS")) > 0 else None, - "subarea_extent":ast.literal_eval(os.environ.get("SUBAREA_EXTENT")) \ - if len(os.environ.get("SUBAREA_EXTENT")) > 0 else None, + "wmo_nmbrs": ast.literal_eval(os.environ.get("WMO_BLOCK_STTN_NMBRS")) + if len(os.environ.get("WMO_BLOCK_STTN_NMBRS")) > 0 + else None, + "subarea_extent": ast.literal_eval(os.environ.get("SUBAREA_EXTENT")) + if len(os.environ.get("SUBAREA_EXTENT")) > 0 + else None, "obs_interval": isodate.parse_duration(os.environ["SURFACE_SYNOP_INTERVAL"]), "obs_offset": isodate.parse_duration(os.environ["SURFACE_SYNOP_OFFSET"]), "rose_datac": os.environ["ROSE_DATAC"], @@ -292,7 +295,7 @@ def fetch_data(file_retriever: FileRetrieverABC): raise FileNotFoundError("No files found for model!") -def fetch_obs(obs_retriever: FileRetrieverABC = FilesystemFileRetriever): +def fetch_obs(obs_retriever: FileRetrieverABC): """Fetch the observations corresponding to a model run. The following environment variables need to be set: @@ -343,9 +346,7 @@ def fetch_obs(obs_retriever: FileRetrieverABC = FilesystemFileRetriever): # Use obs retriever to transfer data with multiple threads. # We shouldn't need to iterate as we do for the forecast data - # because these files will be smaller. Passing None is a temporary - # measure until we can work out why python doesn't think this is a - # class method. + # because these files will be smaller. try: obs_retriever.get_file( paths[0], @@ -359,5 +360,5 @@ def fetch_obs(obs_retriever: FileRetrieverABC = FilesystemFileRetriever): wmo_nmbrs=v["wmo_nmbrs"], subarea_extent=v["subarea_extent"], ) - except: - raise FileNotFoundError("No observations available.") + except Exception as exc: + raise ValueError("No observations available.") from exc diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index eeb604cb9..fe6503214 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -1185,7 +1185,7 @@ def _plot_and_save_postage_stamps_in_single_plot_histogram_series( def _plot_and_save_scattermap_plot( - cube: iris.cube.Cube, filename: str, title: str, projection=None: str, **kwargs + cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs ): """Plot and save a geographical scatter plot. @@ -1203,13 +1203,13 @@ def _plot_and_save_scattermap_plot( """ # Setup plot details, size, resolution, etc. fig = plt.figure(figsize=(15, 15), facecolor="w", edgecolor="k") - if (projection is not None): + if projection is not None: # Apart from the default, the only projection we currently support is # a stereographic projection over the North Pole. - if (projection == 'NP_Stereo'): + if projection == "NP_Stereo": axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0)) else: - Raise ValueError(f"Unknown projection: {projection}") + raise ValueError(f"Unknown projection: {projection}") else: axes = plt.axes(projection=ccrs.PlateCarree()) @@ -1233,7 +1233,7 @@ def _plot_and_save_scattermap_plot( # Add coastlines. try: - axes.coastlines() + axes.coastlines(resolution="10m") except AttributeError: pass diff --git a/src/CSET/operators/read.py b/src/CSET/operators/read.py index 1a21dbc58..a878378b5 100644 --- a/src/CSET/operators/read.py +++ b/src/CSET/operators/read.py @@ -851,30 +851,33 @@ def _fix_um_winds(cubes: iris.cube.CubeList): # this will be biased low in general because the components will mostly # be time averages. For simplicity, we do this only if there is just one # cube of a component. A more complicated approach would be to consider - # the cell methods, but it may not be warranted. - u_constr = iris.AttributeConstraint(STASH="m01s03i225") - v_constr = iris.AttributeConstraint(STASH="m01s03i226") - speed_constr = iris.AttributeConstraint(STASH="m01s03i227") + # the cell methods, but it may not be warranted. As LFRic does not + # have STASH attributes, we need a try...except... construct. try: - if cubes.extract(u_constr) and cubes.extract(v_constr): - if len(cubes.extract(u_constr)) == 1 and not cubes.extract(speed_constr): - _add_wind_speed_um(cubes) - # Convert winds in the UM to be relative to true east and true north. + u10m = cubes.extract(iris.AttributeConstraint(STASH="m01s03i225")) + v10m = cubes.extract(iris.AttributeConstraint(STASH="m01s03i226")) + m10m = cubes.extract(iris.AttributeConstraint(STASH="m01s03i227")) + if len(u10m) == 1 and len(v10m) == 1 and len(m10m) == 0: + for i in range(len(u10m)): + cubes.append(_add_wind_speed_um(u10m[i], v10m[i])) + # Convert winds in the UM to be relative to true east and + # true north. _convert_wind_true_dirn_um(cubes) - except (KeyError, AttributeError): + except ValueError: pass -def _add_wind_speed_um(cubes: iris.cube.CubeList): - """Add windspeeds to cubes from the UM.""" - wspd10 = ( - cubes.extract_cube(iris.AttributeConstraint(STASH="m01s03i225"))[0] ** 2 - + cubes.extract_cube(iris.AttributeConstraint(STASH="m01s03i226"))[0] ** 2 - ) ** 0.5 +def _add_wind_speed_um(u10m: iris.cube.Cube, v10m: iris.cube.Cube) -> iris.cube.Cube: + """Calculate windspeed at 10m from components. + + This is intended only for use with the UM and the result is + given a STASH attribute. + """ + wspd10 = (u10m**2 + v10m**2) ** 0.5 wspd10.attributes["STASH"] = "m01s03i227" wspd10.standard_name = "wind_speed" wspd10.long_name = "wind_speed_at_10m" - cubes.append(wspd10) + return wspd10 def _convert_wind_true_dirn_um(cubes: iris.cube.CubeList): diff --git a/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml b/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml index 580b0b31e..4c8ad1448 100644 --- a/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml +++ b/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml @@ -11,24 +11,25 @@ steps: minuend: operator: regrid.interpolate_to_point_cube fld: - operator: filters.filter_cubes - constraint: - operator: constraints.combine_constraints - var_constraint: - operator: constraints.generate_var_constraint - varname: $VARNAME - cell_method_constraint: - operator: constraints.generate_cell_methods_constraint - cell_methods: [] - # This will require more work. Several fields are stored with muliple - # cell methods, so generally we want the instantaneous one, but in - # other cases, the field may be available only with another cell method. - # In effect, we shall need to match cell methods to diagnostics. + operator: filters.filter_cubes + constraint: + operator: constraints.combine_constraints + var_constraint: + operator: constraints.generate_var_constraint + varname: $VARNAME + cell_method_constraint: + operator: constraints.generate_cell_methods_constraint + cell_methods: [] + # This will require more work. Several fields are stored + # with multiple cell methods, so generally we want the + # instantaneous one, but in other cases, the field may + # be available only with another cell method. In effect, + # we shall need to match cell methods to diagnostics. point_cube: operator: filters.filter_cubes constraint: operator: constraints.generate_var_constraint - varname: $OBSVARNAME + varname: $OBSVARNAME subtrahend: operator: filters.filter_cubes constraint: diff --git a/src/CSET/recipes/generic_obs_scatterplot.yaml b/src/CSET/recipes/generic_obs_scatterplot.yaml index b2260a382..3bf32a73a 100644 --- a/src/CSET/recipes/generic_obs_scatterplot.yaml +++ b/src/CSET/recipes/generic_obs_scatterplot.yaml @@ -7,10 +7,10 @@ steps: - operator: read.read_cubes file_paths: $INPUT_PATHS constraint: - operator: constraints.combine_constraints - variable_constraint: - operator: constraints.generate_var_constraint - varname: $OBSVARNAME + operator: constraints.combine_constraints + variable_constraint: + operator: constraints.generate_var_constraint + varname: $OBSVARNAME - operator: plot.spatial_contour_plot sequence_coordinate: time From c82020937778fe0c4a760625fafa94f5b90223de Mon Sep 17 00:00:00 2001 From: "John M. Edwards" Date: Mon, 4 Aug 2025 16:53:56 +0100 Subject: [PATCH 3/7] Adding a doc-string. --- src/CSET/operators/regrid.py | 27 ++++++++++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/src/CSET/operators/regrid.py b/src/CSET/operators/regrid.py index 9722122e4..8ea27ee1d 100644 --- a/src/CSET/operators/regrid.py +++ b/src/CSET/operators/regrid.py @@ -365,9 +365,30 @@ def transform_lat_long_points(lon, lat, cube): return lon_rot, lat_rot -def interpolate_to_point_cube(fld: iris.cube.Cube | iris.cube.CubeList, - point_cube: iris.cube.Cube, - **kwargs) -> iris.cube.Cube | iris.cube.CubeList: + +def interpolate_to_point_cube( + fld: iris.cube.Cube | iris.cube.CubeList, point_cube: iris.cube.Cube, **kwargs +) -> iris.cube.Cube | iris.cube.CubeList: + """Interpolate from a 2D field to a set of points. + + Interpolate the 2D field in fld to the set of points + specified in point_cube. + + Parameters + ---------- + fld: Cube + An iris cube containing a two-dimensional field. + point_cube: Cube + An iris cube specifying the point to which the data + will be interpolated. + + Returns + ------- + fld_point_cube: Cube + An iris cube containing interpolated values at the points + specified by the point cube. + + """ # # As a basis, create a copy of the point cube. fld_point_cube = point_cube.copy() From 18f8a0fc4709f59e9c619ad3e8e14c315db7258a Mon Sep 17 00:00:00 2001 From: "John M. Edwards" Date: Thu, 31 Jul 2025 09:30:09 +0100 Subject: [PATCH 4/7] Replicating changes in a new branch. --- src/CSET/cset_workflow/.gitignore | 1 + .../app/fetch_fcst/bin/fetch_data.py | 91 ++++++++++++++ .../cset_workflow/app/fetch_obs/rose-app.conf | 3 + src/CSET/cset_workflow/flow.cylc | 28 +++++ .../spatial_surface_model_obs_diference.cylc | 26 ++++ .../includes/spatial_surface_obs_field.cylc | 18 +++ .../meta/observations/rose-meta.conf | 116 ++++++++++++++++++ src/CSET/cset_workflow/meta/rose-meta.conf | 1 + .../cset_workflow/meta/rose-meta.conf.jinja2 | 11 ++ .../cset_workflow/rose-suite.conf.example | 8 ++ src/CSET/operators/plot.py | 85 ++++++++++++- src/CSET/operators/read.py | 11 ++ src/CSET/operators/regrid.py | 68 ++++++++++ ...eric_model_obs_difference_scatterplot.yaml | 44 +++++++ src/CSET/recipes/generic_obs_scatterplot.yaml | 20 +++ 15 files changed, 529 insertions(+), 2 deletions(-) create mode 100644 src/CSET/cset_workflow/app/fetch_obs/rose-app.conf create mode 100644 src/CSET/cset_workflow/includes/spatial_surface_model_obs_diference.cylc create mode 100644 src/CSET/cset_workflow/includes/spatial_surface_obs_field.cylc create mode 100644 src/CSET/cset_workflow/meta/observations/rose-meta.conf create mode 100644 src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml create mode 100644 src/CSET/recipes/generic_obs_scatterplot.yaml diff --git a/src/CSET/cset_workflow/.gitignore b/src/CSET/cset_workflow/.gitignore index 6932c7da3..a49976707 100644 --- a/src/CSET/cset_workflow/.gitignore +++ b/src/CSET/cset_workflow/.gitignore @@ -2,6 +2,7 @@ site/ opt/ demo_pointstat/ +*metdb* restricted* # User workflow configuration. diff --git a/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py b/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py index 39b4c92f2..dbdd168c5 100644 --- a/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py +++ b/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py @@ -15,6 +15,7 @@ """Retrieve the files from the filesystem for the current cycle point.""" import abc +import ast import glob import itertools import logging @@ -176,6 +177,25 @@ def _get_needed_environment_variables() -> dict: logging.debug("Environment variables loaded: %s", variables) return variables +def _get_needed_environment_variables_obs() -> dict: + """Load the needed variables from the environment.""" + variables = { + "subtype": os.environ["OBS_SUBTYPE"], + "data_time": datetime.fromisoformat(os.environ["CYLC_TASK_CYCLE_POINT"]), + "forecast_length": isodate.parse_duration(os.environ["ANALYSIS_LENGTH"]), + "obs_fields": ast.literal_eval(os.environ["SURFACE_SYNOP_FIELDS"]), + "model_identifier": "OBS", + "wmo_nmbrs": ast.literal_eval(os.environ.get("WMO_BLOCK_STTN_NMBRS")) \ + if len(os.environ.get("WMO_BLOCK_STTN_NMBRS")) > 0 else None, + "subarea_extent":ast.literal_eval(os.environ.get("SUBAREA_EXTENT")) \ + if len(os.environ.get("SUBAREA_EXTENT")) > 0 else None, + "obs_interval": isodate.parse_duration(os.environ["SURFACE_SYNOP_INTERVAL"]), + "obs_offset": isodate.parse_duration(os.environ["SURFACE_SYNOP_OFFSET"]), + "rose_datac": os.environ["ROSE_DATAC"], + } + logging.debug("Environment variables loaded: %s", variables) + return variables + def _template_file_path( raw_path: str, @@ -270,3 +290,74 @@ def fetch_data(file_retriever: FileRetrieverABC): # exiting the with block. if not files_found: raise FileNotFoundError("No files found for model!") + + +def fetch_obs(obs_retriever: FileRetrieverABC = FilesystemFileRetriever): + """Fetch the observations corresponding to a model run. + + The following environment variables need to be set: + * ANALYSIS_OFFSET + * ANALYSIS_LENGTH + * CYLC_TASK_CYCLE_POINT + * DATA_PATH + * DATA_PERIOD + * DATE_TYPE + * MODEL_IDENTIFIER + * ROSE_DATAC + + Parameters + ---------- + obs_retriever: ObsRetriever + ObsRetriever implementation to use. Defaults to FilesystemFileRetriever. + + Raises + ------ + FileNotFound: + If no observations are available. + """ + v = _get_needed_environment_variables_obs() + + # Prepare output directory. + cycle_obs_dir = f"{v['rose_datac']}/data/OBS" + os.makedirs(cycle_obs_dir, exist_ok=True) + logging.debug("Output directory: %s", cycle_obs_dir) + + # We will get just one file for now, but follow the templating + # syntax for the model for consistency. + obs_base_path = ( + v["subtype"] + + "_" + + "%Y%m%dT%H%MZ_dt_" + + str(int(v["forecast_length"].total_seconds() // 3600)).zfill(3) + + ".nc" + ) + paths = _template_file_path( + obs_base_path, + "initiation", + v["data_time"], + v["forecast_length"], + timedelta(seconds=0), + v["obs_interval"], + ) + logging.info("Retrieving paths:\n%s", "\n".join(paths)) + + # Use obs retriever to transfer data with multiple threads. + # We shouldn't need to iterate as we do for the forecast data + # because these files will be smaller. Passing None is a temporary + # measure until we can work out why python doesn't think this is a + # class method. + try: + obs_retriever.get_file( + paths[0], + v["subtype"], + v["obs_fields"], + v["data_time"], + v["obs_offset"], + v["forecast_length"], + v["obs_interval"], + cycle_obs_dir, + wmo_nmbrs=v["wmo_nmbrs"], + subarea_extent=v["subarea_extent"], + ) + except: + raise FileNotFoundError("No observations available.") diff --git a/src/CSET/cset_workflow/app/fetch_obs/rose-app.conf b/src/CSET/cset_workflow/app/fetch_obs/rose-app.conf new file mode 100644 index 000000000..ec3d1e950 --- /dev/null +++ b/src/CSET/cset_workflow/app/fetch_obs/rose-app.conf @@ -0,0 +1,3 @@ +[command] +default=echo "Please set ROSE_APP_COMMAND_KEY to your database."; false +metdb=app_env_wrapper fetch-obs-metdb.py diff --git a/src/CSET/cset_workflow/flow.cylc b/src/CSET/cset_workflow/flow.cylc index ddf05ecfe..e725bd7dd 100644 --- a/src/CSET/cset_workflow/flow.cylc +++ b/src/CSET/cset_workflow/flow.cylc @@ -42,6 +42,7 @@ final cycle point = {{CSET_TRIAL_END_DATE}} R1/{{date}} = """ setup_complete[^] => FETCH_DATA:succeed-all => + FETCH_OBS:succeed-all => fetch_complete => PROCESS:{% if LOGLEVEL == "DEBUG" %}succeed-all{% else %}finish-all{% endif %} => cycle_complete @@ -52,6 +53,7 @@ final cycle point = {{CSET_TRIAL_END_DATE}} {{CSET_TRIAL_CYCLE_PERIOD}} = """ setup_complete[^] => FETCH_DATA:succeed-all => + FETCH_OBS:succeed-all => fetch_complete => PROCESS:{% if LOGLEVEL == "DEBUG" %}succeed-all{% else %}finish-all{% endif %} => cycle_complete @@ -109,6 +111,8 @@ final cycle point = {{CSET_TRIAL_END_DATE}} execution retry delays = PT1M, PT15M, PT1H [[[directives]]] --mem=10000 + [[[environment]]] + PLOTTING_PROJECTION = {{PLOTTING_PROJECTION | default("")}} [[PROCESS_CASE_AGGREGATION]] script = rose task-run -v --app-key=run_cset_recipe @@ -125,6 +129,12 @@ final cycle point = {{CSET_TRIAL_END_DATE}} [[[environment]]] ANALYSIS_LENGTH = {{ANALYSIS_LENGTH}} + [[FETCH_OBS]] + script = rose task-run -v --app-key=fetch_obs + execution time limit = PT1H + [[[environment]]] + ANALYSIS_LENGTH = {{ANALYSIS_LENGTH}} + [[METPLUS]] [[[environment]]] {% if METPLUS_GRID_STAT|default(False) %} @@ -150,6 +160,9 @@ final cycle point = {{CSET_TRIAL_END_DATE}} [[cycle_complete]] inherit = DUMMY_TASK + [[dummy_fetch_obs]] + inherit = DUMMY_TASK, FETCH_OBS + [[dummy_process]] inherit = DUMMY_TASK, PROCESS @@ -185,6 +198,21 @@ final cycle point = {{CSET_TRIAL_END_DATE}} ANALYSIS_OFFSET = {{model["analysis_offset"]}} {% endfor %} + {% if SURFACE_SYNOP_OBS %} + [[fetch_obs]] + # Fetch observations from the MetDB. + inherit = FETCH_OBS + [[[environment]]] + MODEL_IDENTIFIER = "OBS" + ROSE_APP_COMMAND_KEY = "metdb" + SURFACE_SYNOP_FIELDS = {{SURFACE_SYNOP_FIELDS}} + SURFACE_SYNOP_INTERVAL = {{SURFACE_SYNOP_INTERVAL}} + SURFACE_SYNOP_OFFSET = {{SURFACE_SYNOP_OFFSET}} + WMO_BLOCK_STTN_NMBRS = {{WMO_BLOCK_STTN_NMBRS | default("")}} + SUBAREA_EXTENT = {{SUBAREA_EXTENT | default("")}} + OBS_SUBTYPE = "LNDSYB" + {% endif %} + [[housekeeping]] # Housekeep input data files. [[[environment]]] diff --git a/src/CSET/cset_workflow/includes/spatial_surface_model_obs_diference.cylc b/src/CSET/cset_workflow/includes/spatial_surface_model_obs_diference.cylc new file mode 100644 index 000000000..2b5b92be0 --- /dev/null +++ b/src/CSET/cset_workflow/includes/spatial_surface_model_obs_diference.cylc @@ -0,0 +1,26 @@ +{% if SURFACE_SYNOP_OBS and SURFACE_SYNOP_DIFFS %} +{% for model in models %} +{% for obs_field in (SURFACE_SYNOP_FIELDS | list) %} +{# {% set model_field = (SURFACE_SYNOP_MODEL_FIELDS | list)[SURFACE_SYNOP_FIELDS.index(obs_fields)] %} #} +[runtime] + [[obs_model_difference_scatterplot_m{{model["id"]}}_{{obs_field}}]] + inherit = PROCESS + execution time limit = PT6H + [[[environment]]] + CSET_RECIPE_NAME = "generic_model_obs_difference_scatterplot.yaml" + CSET_ADDOPTS = """ + --OBSVARNAME={{obs_field}} + --VARNAME={{(SURFACE_SYNOP_MODEL_FIELDS | list)[loop.index0]}} + --BASE_MODEL='{{model["id"]}}' + --OBSERVATIONS='OBS' + --MODEL_NAME='{{model["name"]}}' + --DIFFERENCE=True + --SUBAREA_TYPE={{SUBAREA_TYPE | default("")}} + --SUBAREA_EXTENT={{SUBAREA_TYPE | default("")}} + --PLOTTING_PROJECTION={{PLOTTING_PROJECTION | default("")}} + """ + MODEL_IDENTIFIERS = 'OBS {{model["id"]}}' + +{% endfor %} +{% endfor %} +{% endif %} diff --git a/src/CSET/cset_workflow/includes/spatial_surface_obs_field.cylc b/src/CSET/cset_workflow/includes/spatial_surface_obs_field.cylc new file mode 100644 index 000000000..f0ce429d9 --- /dev/null +++ b/src/CSET/cset_workflow/includes/spatial_surface_obs_field.cylc @@ -0,0 +1,18 @@ +{% if SURFACE_SYNOP_OBS %} +{% for obs_field in SURFACE_SYNOP_FIELDS %} +[runtime] + [[generic_obs_scatterplot_{{obs_field}}]] + inherit = PROCESS + execution time limit = PT1H + [[[environment]]] + CSET_RECIPE_NAME = "generic_obs_scatterplot.yaml" + CSET_ADDOPTS = """ + --OBSVARNAME={{obs_field}} + --SUBAREA_TYPE={{SUBAREA_TYPE | default("")}} + --SUBAREA_EXTENT={{SUBAREA_TYPE | default("")}} + --PLOTTING_PROJECTION={{PLOTTING_PROJECTION | default("")}} + """ + MODEL_IDENTIFIERS = "OBS" + +{% endfor %} +{% endif %} diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf new file mode 100644 index 000000000..b25be25d4 --- /dev/null +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -0,0 +1,116 @@ +################################################################################ +# Observations +################################################################################ + +[Observations] +ns=Observations +sort-key=sec-d +title=Observational Comparisons + +################################################################################ +# Spatial (2D map) fields +################################################################################ + +[Observations/Types] +ns=Observations/Types +sort-key=sec-d1 +title=Types of observations + +# Surface synoptic fields. +[template variables=SURFACE_SYNOP_OBS] +ns=Observations/Types +description=Compare with surface synoptic observations +type=python_boolean +compulsory=true +sort-key=0synop1a +trigger=template variables=SURFACE_SYNOP_FIELDS: True; + template variables=USE_WMO_STATION_NUMBERS: True; + template variables=SURFACE_SYNOP_INTERVAL: True; + template variables=SURFACE_SYNOP_OFFSET: True; + +[template variables=PLACEHOLDER_OBS] +ns=Observations/Types +description=Compare with observations from an undefined source. +help=This variable exists as a placeholder to show the intended + structure the GUI. Additional sources of observations should be + added using one window for each source. + +type=python_boolean +compulsory=true +sort-key=0yopp1a + +################################################################ + +[Observations/Types/Synoptic] +ns=Observations/Types/Synoptic +sort-key=sec-d2 +title=Synoptic observations + +[template variables=SURFACE_SYNOP_FIELDS] +ns=Observations/Types/Synoptic +title=Surface synoptic observations +description=Synoptic observations +help=Name of field in the style of the MetDB. +compulsory=true +type=python_list +sort-key=0synop1b + +[template variables=SURFACE_SYNOP_DIFFS] +ns=Observations/Types/Synoptic +description=Calculate the corresponding differences. +type=python_boolean +compulsory=true +sort-key=0synop1c +trigger=template variables=SURFACE_SYNOP_MODEL_FIELDS: True; + +[template variables=SURFACE_SYNOP_MODEL_FIELDS] +ns=Observations/Types/Synoptic +title=Corresponding model variables +description=Synoptic observations +help=Name of field in the model: the order must match that of the synoptic observations +compulsory=true +type=python_list +sort-key=0synop1d + +[template variables=USE_WMO_STATION_NUMBERS] +ns=Observations/Types/Synoptic +description=Use WMO numbers +help=Select to set stations by WMO station or block number. +compulsory=true +type=python_boolean +sort-key=0synop1e +trigger=template variables=WMO_BLOCK_STTN_NMBRS: True; + +[template variables=WMO_BLOCK_STTN_NMBRS] +ns=Observations/Types/Synoptic +title=WMO Block or station numbers. +description=WMO numbers +help=WMO block numbers (2 digits) or station numbers (5 digits). +compulsory=true +type=python_list +sort-key=0synop1f + +[template variables=SURFACE_SYNOP_INTERVAL] +ns=Observations/Types/Synoptic +title=Interval between synoptic observations +description=Synoptic observations +help=The time interval between rerievals of synoptic observations. +compulsory=true +type=quoted +sort-key=0synop1g + +[template variables=SURFACE_SYNOP_OFFSET] +ns=Observations/Types/Synoptic +title=Offset for synoptic observations +description=Offset for observations +help=The offset of the first retrieved observations relative to the initial time. +compulsory=true +type=quoted +sort-key=0synop1h + +################################################################ + +[Observations/Types/Placeholder] +ns=Observations/Types/Placeholder +sort-key=sec-e2 +title=Observations from an undefined source. diff --git a/src/CSET/cset_workflow/meta/rose-meta.conf b/src/CSET/cset_workflow/meta/rose-meta.conf index 99a67b25a..d73c5f324 100644 --- a/src/CSET/cset_workflow/meta/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/rose-meta.conf @@ -1,5 +1,6 @@ # Diagnostics settings are split into a separate file. import=meta/diagnostics + =meta/observations [ns=Setup] sort-key=sec-a diff --git a/src/CSET/cset_workflow/meta/rose-meta.conf.jinja2 b/src/CSET/cset_workflow/meta/rose-meta.conf.jinja2 index a2ec063ac..ec33e2e9f 100644 --- a/src/CSET/cset_workflow/meta/rose-meta.conf.jinja2 +++ b/src/CSET/cset_workflow/meta/rose-meta.conf.jinja2 @@ -158,6 +158,17 @@ type=quoted compulsory=true sort-key=setup-h-out2 +[template variables=PLOTTING_PROJECTION] +ns=Setup +title=Projection for the plots +description=Geographical projection of plots +help=This is is the projection used for the final plots, specified as a string. + If no value is given, native latitudes and longitudes are used. +compulsory=true +values="Unspecified", "NP_Stereo" +value-titles=Unspecified, North polar stereographic +sort-key=setup-h-out3a + [template variables=PLOT_RESOLUTION] ns=Setup title=Plot resolution diff --git a/src/CSET/cset_workflow/rose-suite.conf.example b/src/CSET/cset_workflow/rose-suite.conf.example index 2c6b93a72..a7a419cec 100644 --- a/src/CSET/cset_workflow/rose-suite.conf.example +++ b/src/CSET/cset_workflow/rose-suite.conf.example @@ -53,9 +53,11 @@ MODEL_LEVEL_FIELDS=[] !!MODULES_LIST= !!MODULES_PURGE=True !!ONE_TO_ONE=False +PLACEHOLDER_OBS=False !!PLEVEL_TRANSECT_FINISHCOORDS= !!PLEVEL_TRANSECT_STARTCOORDS= PLOT_RESOLUTION=100 +!!PLOTTING_PROJECTION="" PRESSURE_LEVELS=[] PRESSURE_LEVEL_FIELDS=[] PROFILE_MLEVEL=False @@ -85,6 +87,11 @@ SPATIAL_SURFACE_FIELD_METHOD=[""] !!SUBAREA_EXTENT=0,0,0,0 !!SUBAREA_TYPE="realworld" SURFACE_FIELDS=[] +SURFACE_SYNOP_OBS=False +!!SURFACE_SYNOP_FIELDS=[] +!!WMO_BLOCK_STTN_NMBRS=[] +!!SURFACE_SYNOP_INTERVAL="" +!!SURFACE_SYNOP_OFFSET="" SURFACE_SINGLE_POINT_TIME_SERIES=False TIMESERIES_MLEVEL_FIELD=False TIMESERIES_MLEVEL_FIELD_AGGREGATION=False,False,False,False @@ -92,6 +99,7 @@ TIMESERIES_PLEVEL_FIELD=False TIMESERIES_PLEVEL_FIELD_AGGREGATION=False,False,False,False TIMESERIES_SURFACE_FIELD=False TIMESERIES_SURFACE_FIELD_AGGREGATION=False,False,False,False +!!USE_WMO_STATION_NUMBERS=False !!VERTICAL_COORDINATE_A=[] !!VERTICAL_COORDINATE_B=[] WEB_ADDR="" diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index 5f5d7b90a..c192c1080 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -25,6 +25,7 @@ import sys from typing import Literal +import cartopy.crs as ccrs import iris import iris.coords import iris.cube @@ -1218,12 +1219,79 @@ def _plot_and_save_postage_stamps_in_single_plot_histogram_series( plt.close(fig) +def _plot_and_save_scattermap_plot( + cube: iris.cube.Cube, filename: str, title: str, projection=None: str, **kwargs +): + """Plot and save a geographical scatter plot. + + Parameters + ---------- + cube: Cube + 1 dimensional Cube of the data points with auxiliary latitude and + longitude coordinates, + filename: str + Filename of the plot to write. + title: str + Plot title. + projection: str + Mapping projection to be used by cartopy. + """ + # Setup plot details, size, resolution, etc. + fig = plt.figure(figsize=(15, 15), facecolor="w", edgecolor="k") + if (projection is not None): + # Apart from the default, the only projection we currently support is + # a stereographic projection over the North Pole. + if (projection == 'NP_Stereo'): + axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0)) + else: + Raise ValueError(f"Unknown projection: {projection}") + else: + axes = plt.axes(projection=ccrs.PlateCarree()) + + # Filled contour plot of the field. + mrk_size = int(np.sqrt(5000000.0 / len(cube.data))) + klon = None + klat = None + for kc in range(len(cube.aux_coords)): + if cube.aux_coords[kc].standard_name == "latitude": + klat = kc + elif cube.aux_coords[kc].standard_name == "longitude": + klon = kc + scatter_map = iplt.scatter( + cube.aux_coords[klon], + cube.aux_coords[klat], + c=cube.data[:], + s=mrk_size, + cmap="jet", + edgecolors="k", + ) + + # Add coastlines. + try: + axes.coastlines() + except AttributeError: + pass + + # Add title. + axes.set_title(title, fontsize=16) + + # Add colour bar. + cbar = fig.colorbar(scatter_map) + cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20) + + # Save plot. + fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution()) + logging.info("Saved geographical scatter plot to %s", filename) + plt.close(fig) + + def _spatial_plot( method: Literal["contourf", "pcolormesh"], cube: iris.cube.Cube, filename: str | None, sequence_coordinate: str, stamp_coordinate: str, + **kwargs, ): """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube. @@ -1274,6 +1342,14 @@ def _spatial_plot( except iris.exceptions.CoordinateNotFoundError: pass + # Produce a geographical scatter plot if the data have a + # dimension called observation or model_obs_error + if any( + crd.var_name == "station" or crd.var_name == "model_obs_error" + for crd in cube.coords() + ): + plotting_func = _plot_and_save_scattermap_plot + # Must have a sequence coordinate. try: cube.coord(sequence_coordinate) @@ -1301,6 +1377,7 @@ def _spatial_plot( stamp_coordinate=stamp_coordinate, title=title, method=method, + **kwargs, ) plot_index.append(plot_filename) @@ -1474,7 +1551,9 @@ def spatial_contour_plot( TypeError If the cube isn't a single cube. """ - _spatial_plot("contourf", cube, filename, sequence_coordinate, stamp_coordinate) + _spatial_plot( + "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs + ) return cube @@ -1523,7 +1602,9 @@ def spatial_pcolormesh_plot( TypeError If the cube isn't a single cube. """ - _spatial_plot("pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate) + _spatial_plot( + "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs + ) return cube diff --git a/src/CSET/operators/read.py b/src/CSET/operators/read.py index a0fe8ad4c..64676f465 100644 --- a/src/CSET/operators/read.py +++ b/src/CSET/operators/read.py @@ -507,6 +507,17 @@ def _lfric_time_coord_fix_callback(cube: iris.cube.Cube, field, filename): not isinstance(time_coord, iris.coords.DimCoord) and len(cube.coord_dims(time_coord)) == 1 ): + # Fudge the bounds to foil checking for strict monotonicity. + if time_coord.has_bounds(): + if (time_coord.bounds[-1][0] - time_coord.bounds[0][0]) < 1.0e-8: + time_coord.bounds = [ + [ + time_coord.bounds[i][0] + 1.0e-8 * float(i), + time_coord.bounds[i][1], + ] + for i in range(len(time_coord.bounds)) + ] + iris.util.promote_aux_coord_to_dim_coord(cube, time_coord) # Force single-valued coordinates to be scalar coordinates. diff --git a/src/CSET/operators/regrid.py b/src/CSET/operators/regrid.py index c24cac73d..9722122e4 100644 --- a/src/CSET/operators/regrid.py +++ b/src/CSET/operators/regrid.py @@ -364,3 +364,71 @@ def transform_lat_long_points(lon, lat, cube): lat_rot = rot_coords[1] return lon_rot, lat_rot + +def interpolate_to_point_cube(fld: iris.cube.Cube | iris.cube.CubeList, + point_cube: iris.cube.Cube, + **kwargs) -> iris.cube.Cube | iris.cube.CubeList: + # + # As a basis, create a copy of the point cube. + fld_point_cube = point_cube.copy() + # Get indices of positional coordinates. We assume that the + # point cube is unrotated. + klon = None + klat = None + for kc in range(len(fld_point_cube.aux_coords)): + if fld_point_cube.aux_coords[kc].standard_name == "latitude": + klat = kc + elif fld_point_cube.aux_coords[kc].standard_name == "longitude": + klon = kc + # + # The input may have a rotated coordinate system. + if len(fld.coords("grid_latitude")) > 0: + # Interpolate in rotated coordinates. + rot_csyst = fld.coords("grid_latitude")[0].coord_system + rotpt = iris.analysis.cartography.rotate_pole( + fld_point_cube.aux_coords[klon].points, + fld_point_cube.aux_coords[klat].points, + rot_csyst.grid_north_pole_longitude, + rot_csyst.grid_north_pole_latitude, + ) + # Add other interpolation options later. + fld_interpolator = iris.analysis.Linear(extrapolation_mode="mask").interpolator( + fld, ["time", "grid_latitude", "grid_longitude"] + ) + for jt in range(len(fld_point_cube.coords("time")[0].points)): + fld_point_cube.data[jt, :] = np.ma.masked_invalid( + [ + fld_interpolator( + [ + fld_point_cube.coord("time").points[jt], + rotpt[1][k], + rotpt[0][k], + ] + ).data + if ~point_cube.data.mask[jt][k] + else np.nan + for k in range(len(rotpt[0])) + ] + ) + else: + # Add other interpolation options later. + fld_interpolator = iris.analysis.Linear(extrapolation_mode="mask").interpolator( + fld, ["time", "latitude", "longitude"] + ) + for jt in range(len(fld_point_cube.coords("time")[0].points)): + fld_point_cube.data[jt, :] = np.ma.masked_invalid( + [ + fld_interpolator( + [ + fld_point_cube.coords("time")[0].points[jt], + fld_point_cube.coord("latitude").points[k], + fld_point_cube.coord("longitude").points[k], + ] + ).data + if ~point_cube.data.mask[jt][k] + else np.nan + for k in range(fld_point_cube.coord("latitude").points) + ] + ) + # + return fld_point_cube diff --git a/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml b/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml new file mode 100644 index 000000000..580b0b31e --- /dev/null +++ b/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml @@ -0,0 +1,44 @@ +category: Quick Look +title: Difference between modelled field of $VARNAME in model $MODEL_NAME and observation. +description: | + Plot some observations. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS + + - operator: misc.subtraction + minuend: + operator: regrid.interpolate_to_point_cube + fld: + operator: filters.filter_cubes + constraint: + operator: constraints.combine_constraints + var_constraint: + operator: constraints.generate_var_constraint + varname: $VARNAME + cell_method_constraint: + operator: constraints.generate_cell_methods_constraint + cell_methods: [] + # This will require more work. Several fields are stored with muliple + # cell methods, so generally we want the instantaneous one, but in + # other cases, the field may be available only with another cell method. + # In effect, we shall need to match cell methods to diagnostics. + point_cube: + operator: filters.filter_cubes + constraint: + operator: constraints.generate_var_constraint + varname: $OBSVARNAME + subtrahend: + operator: filters.filter_cubes + constraint: + operator: constraints.generate_var_constraint + varname: $OBSVARNAME + + + - operator: plot.spatial_contour_plot + sequence_coordinate: time + projection: $PLOTTING_PROJECTION + + - operator: write.write_cube_to_nc + overwrite: True diff --git a/src/CSET/recipes/generic_obs_scatterplot.yaml b/src/CSET/recipes/generic_obs_scatterplot.yaml new file mode 100644 index 000000000..b2260a382 --- /dev/null +++ b/src/CSET/recipes/generic_obs_scatterplot.yaml @@ -0,0 +1,20 @@ +category: Quick Look +title: Observed values of $OBSVARNAME. +description: | + Plot some observations. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS + constraint: + operator: constraints.combine_constraints + variable_constraint: + operator: constraints.generate_var_constraint + varname: $OBSVARNAME + + - operator: plot.spatial_contour_plot + sequence_coordinate: time + projection: $PLOTTING_PROJECTION + + - operator: write.write_cube_to_nc + overwrite: True From 017e5e50cad100e078f5c39de79f9611e42e3af1 Mon Sep 17 00:00:00 2001 From: "John M. Edwards" Date: Mon, 4 Aug 2025 11:48:28 +0100 Subject: [PATCH 5/7] Committing tested version. --- .../app/fetch_fcst/bin/fetch_data.py | 21 +++++------ src/CSET/operators/plot.py | 10 +++--- src/CSET/operators/read.py | 35 ++++++++++--------- ...eric_model_obs_difference_scatterplot.yaml | 29 +++++++-------- src/CSET/recipes/generic_obs_scatterplot.yaml | 8 ++--- 5 files changed, 54 insertions(+), 49 deletions(-) diff --git a/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py b/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py index dbdd168c5..d46e2934c 100644 --- a/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py +++ b/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py @@ -177,6 +177,7 @@ def _get_needed_environment_variables() -> dict: logging.debug("Environment variables loaded: %s", variables) return variables + def _get_needed_environment_variables_obs() -> dict: """Load the needed variables from the environment.""" variables = { @@ -185,10 +186,12 @@ def _get_needed_environment_variables_obs() -> dict: "forecast_length": isodate.parse_duration(os.environ["ANALYSIS_LENGTH"]), "obs_fields": ast.literal_eval(os.environ["SURFACE_SYNOP_FIELDS"]), "model_identifier": "OBS", - "wmo_nmbrs": ast.literal_eval(os.environ.get("WMO_BLOCK_STTN_NMBRS")) \ - if len(os.environ.get("WMO_BLOCK_STTN_NMBRS")) > 0 else None, - "subarea_extent":ast.literal_eval(os.environ.get("SUBAREA_EXTENT")) \ - if len(os.environ.get("SUBAREA_EXTENT")) > 0 else None, + "wmo_nmbrs": ast.literal_eval(os.environ.get("WMO_BLOCK_STTN_NMBRS")) + if len(os.environ.get("WMO_BLOCK_STTN_NMBRS")) > 0 + else None, + "subarea_extent": ast.literal_eval(os.environ.get("SUBAREA_EXTENT")) + if len(os.environ.get("SUBAREA_EXTENT")) > 0 + else None, "obs_interval": isodate.parse_duration(os.environ["SURFACE_SYNOP_INTERVAL"]), "obs_offset": isodate.parse_duration(os.environ["SURFACE_SYNOP_OFFSET"]), "rose_datac": os.environ["ROSE_DATAC"], @@ -292,7 +295,7 @@ def fetch_data(file_retriever: FileRetrieverABC): raise FileNotFoundError("No files found for model!") -def fetch_obs(obs_retriever: FileRetrieverABC = FilesystemFileRetriever): +def fetch_obs(obs_retriever: FileRetrieverABC): """Fetch the observations corresponding to a model run. The following environment variables need to be set: @@ -343,9 +346,7 @@ def fetch_obs(obs_retriever: FileRetrieverABC = FilesystemFileRetriever): # Use obs retriever to transfer data with multiple threads. # We shouldn't need to iterate as we do for the forecast data - # because these files will be smaller. Passing None is a temporary - # measure until we can work out why python doesn't think this is a - # class method. + # because these files will be smaller. try: obs_retriever.get_file( paths[0], @@ -359,5 +360,5 @@ def fetch_obs(obs_retriever: FileRetrieverABC = FilesystemFileRetriever): wmo_nmbrs=v["wmo_nmbrs"], subarea_extent=v["subarea_extent"], ) - except: - raise FileNotFoundError("No observations available.") + except Exception as exc: + raise ValueError("No observations available.") from exc diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index c192c1080..f97fe0276 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -1220,7 +1220,7 @@ def _plot_and_save_postage_stamps_in_single_plot_histogram_series( def _plot_and_save_scattermap_plot( - cube: iris.cube.Cube, filename: str, title: str, projection=None: str, **kwargs + cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs ): """Plot and save a geographical scatter plot. @@ -1238,13 +1238,13 @@ def _plot_and_save_scattermap_plot( """ # Setup plot details, size, resolution, etc. fig = plt.figure(figsize=(15, 15), facecolor="w", edgecolor="k") - if (projection is not None): + if projection is not None: # Apart from the default, the only projection we currently support is # a stereographic projection over the North Pole. - if (projection == 'NP_Stereo'): + if projection == "NP_Stereo": axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0)) else: - Raise ValueError(f"Unknown projection: {projection}") + raise ValueError(f"Unknown projection: {projection}") else: axes = plt.axes(projection=ccrs.PlateCarree()) @@ -1268,7 +1268,7 @@ def _plot_and_save_scattermap_plot( # Add coastlines. try: - axes.coastlines() + axes.coastlines(resolution="10m") except AttributeError: pass diff --git a/src/CSET/operators/read.py b/src/CSET/operators/read.py index 64676f465..289f8bf50 100644 --- a/src/CSET/operators/read.py +++ b/src/CSET/operators/read.py @@ -848,30 +848,33 @@ def _fix_um_winds(cubes: iris.cube.CubeList): # this will be biased low in general because the components will mostly # be time averages. For simplicity, we do this only if there is just one # cube of a component. A more complicated approach would be to consider - # the cell methods, but it may not be warranted. - u_constr = iris.AttributeConstraint(STASH="m01s03i225") - v_constr = iris.AttributeConstraint(STASH="m01s03i226") - speed_constr = iris.AttributeConstraint(STASH="m01s03i227") + # the cell methods, but it may not be warranted. As LFRic does not + # have STASH attributes, we need a try...except... construct. try: - if cubes.extract(u_constr) and cubes.extract(v_constr): - if len(cubes.extract(u_constr)) == 1 and not cubes.extract(speed_constr): - _add_wind_speed_um(cubes) - # Convert winds in the UM to be relative to true east and true north. + u10m = cubes.extract(iris.AttributeConstraint(STASH="m01s03i225")) + v10m = cubes.extract(iris.AttributeConstraint(STASH="m01s03i226")) + m10m = cubes.extract(iris.AttributeConstraint(STASH="m01s03i227")) + if len(u10m) == 1 and len(v10m) == 1 and len(m10m) == 0: + for i in range(len(u10m)): + cubes.append(_add_wind_speed_um(u10m[i], v10m[i])) + # Convert winds in the UM to be relative to true east and + # true north. _convert_wind_true_dirn_um(cubes) - except (KeyError, AttributeError): + except ValueError: pass -def _add_wind_speed_um(cubes: iris.cube.CubeList): - """Add windspeeds to cubes from the UM.""" - wspd10 = ( - cubes.extract_cube(iris.AttributeConstraint(STASH="m01s03i225"))[0] ** 2 - + cubes.extract_cube(iris.AttributeConstraint(STASH="m01s03i226"))[0] ** 2 - ) ** 0.5 +def _add_wind_speed_um(u10m: iris.cube.Cube, v10m: iris.cube.Cube) -> iris.cube.Cube: + """Calculate windspeed at 10m from components. + + This is intended only for use with the UM and the result is + given a STASH attribute. + """ + wspd10 = (u10m**2 + v10m**2) ** 0.5 wspd10.attributes["STASH"] = "m01s03i227" wspd10.standard_name = "wind_speed" wspd10.long_name = "wind_speed_at_10m" - cubes.append(wspd10) + return wspd10 def _convert_wind_true_dirn_um(cubes: iris.cube.CubeList): diff --git a/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml b/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml index 580b0b31e..4c8ad1448 100644 --- a/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml +++ b/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml @@ -11,24 +11,25 @@ steps: minuend: operator: regrid.interpolate_to_point_cube fld: - operator: filters.filter_cubes - constraint: - operator: constraints.combine_constraints - var_constraint: - operator: constraints.generate_var_constraint - varname: $VARNAME - cell_method_constraint: - operator: constraints.generate_cell_methods_constraint - cell_methods: [] - # This will require more work. Several fields are stored with muliple - # cell methods, so generally we want the instantaneous one, but in - # other cases, the field may be available only with another cell method. - # In effect, we shall need to match cell methods to diagnostics. + operator: filters.filter_cubes + constraint: + operator: constraints.combine_constraints + var_constraint: + operator: constraints.generate_var_constraint + varname: $VARNAME + cell_method_constraint: + operator: constraints.generate_cell_methods_constraint + cell_methods: [] + # This will require more work. Several fields are stored + # with multiple cell methods, so generally we want the + # instantaneous one, but in other cases, the field may + # be available only with another cell method. In effect, + # we shall need to match cell methods to diagnostics. point_cube: operator: filters.filter_cubes constraint: operator: constraints.generate_var_constraint - varname: $OBSVARNAME + varname: $OBSVARNAME subtrahend: operator: filters.filter_cubes constraint: diff --git a/src/CSET/recipes/generic_obs_scatterplot.yaml b/src/CSET/recipes/generic_obs_scatterplot.yaml index b2260a382..3bf32a73a 100644 --- a/src/CSET/recipes/generic_obs_scatterplot.yaml +++ b/src/CSET/recipes/generic_obs_scatterplot.yaml @@ -7,10 +7,10 @@ steps: - operator: read.read_cubes file_paths: $INPUT_PATHS constraint: - operator: constraints.combine_constraints - variable_constraint: - operator: constraints.generate_var_constraint - varname: $OBSVARNAME + operator: constraints.combine_constraints + variable_constraint: + operator: constraints.generate_var_constraint + varname: $OBSVARNAME - operator: plot.spatial_contour_plot sequence_coordinate: time From f5cb2ec01d36557ccc9385f04d8149f9211b61ea Mon Sep 17 00:00:00 2001 From: "John M. Edwards" Date: Mon, 4 Aug 2025 16:53:56 +0100 Subject: [PATCH 6/7] Adding a doc-string. --- src/CSET/operators/regrid.py | 27 ++++++++++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/src/CSET/operators/regrid.py b/src/CSET/operators/regrid.py index 9722122e4..8ea27ee1d 100644 --- a/src/CSET/operators/regrid.py +++ b/src/CSET/operators/regrid.py @@ -365,9 +365,30 @@ def transform_lat_long_points(lon, lat, cube): return lon_rot, lat_rot -def interpolate_to_point_cube(fld: iris.cube.Cube | iris.cube.CubeList, - point_cube: iris.cube.Cube, - **kwargs) -> iris.cube.Cube | iris.cube.CubeList: + +def interpolate_to_point_cube( + fld: iris.cube.Cube | iris.cube.CubeList, point_cube: iris.cube.Cube, **kwargs +) -> iris.cube.Cube | iris.cube.CubeList: + """Interpolate from a 2D field to a set of points. + + Interpolate the 2D field in fld to the set of points + specified in point_cube. + + Parameters + ---------- + fld: Cube + An iris cube containing a two-dimensional field. + point_cube: Cube + An iris cube specifying the point to which the data + will be interpolated. + + Returns + ------- + fld_point_cube: Cube + An iris cube containing interpolated values at the points + specified by the point cube. + + """ # # As a basis, create a copy of the point cube. fld_point_cube = point_cube.copy() From 666a9146a107e57b51f11263f5469a660fd354ea Mon Sep 17 00:00:00 2001 From: "John M. Edwards" Date: Tue, 26 Aug 2025 15:25:03 +0100 Subject: [PATCH 7/7] Changes following the review. --- src/CSET/cset_workflow/flow.cylc | 1 + .../spatial_surface_model_obs_diference.cylc | 1 - .../cset_workflow/meta/observations/rose-meta.conf | 11 ----------- src/CSET/cset_workflow/meta/rose-meta.conf.jinja2 | 4 ++-- src/CSET/operators/plot.py | 12 +++++++++--- .../generic_model_obs_difference_scatterplot.yaml | 3 ++- src/CSET/recipes/generic_obs_scatterplot.yaml | 2 +- 7 files changed, 15 insertions(+), 19 deletions(-) diff --git a/src/CSET/cset_workflow/flow.cylc b/src/CSET/cset_workflow/flow.cylc index 835ec33d2..0dfa4f19f 100644 --- a/src/CSET/cset_workflow/flow.cylc +++ b/src/CSET/cset_workflow/flow.cylc @@ -123,6 +123,7 @@ URL = https://metoffice.github.io/CSET [[[directives]]] --mem=100000 [[[environment]]] + PLOTTING_PROJECTION = {{PLOTTING_PROJECTION | default("")}} # As this process is used for case aggregation we hard code to True. DO_CASE_AGGREGATION = True diff --git a/src/CSET/cset_workflow/includes/spatial_surface_model_obs_diference.cylc b/src/CSET/cset_workflow/includes/spatial_surface_model_obs_diference.cylc index 2b5b92be0..8c421e474 100644 --- a/src/CSET/cset_workflow/includes/spatial_surface_model_obs_diference.cylc +++ b/src/CSET/cset_workflow/includes/spatial_surface_model_obs_diference.cylc @@ -1,7 +1,6 @@ {% if SURFACE_SYNOP_OBS and SURFACE_SYNOP_DIFFS %} {% for model in models %} {% for obs_field in (SURFACE_SYNOP_FIELDS | list) %} -{# {% set model_field = (SURFACE_SYNOP_MODEL_FIELDS | list)[SURFACE_SYNOP_FIELDS.index(obs_fields)] %} #} [runtime] [[obs_model_difference_scatterplot_m{{model["id"]}}_{{obs_field}}]] inherit = PROCESS diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf index b25be25d4..22564a4f9 100644 --- a/src/CSET/cset_workflow/meta/observations/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -28,17 +28,6 @@ trigger=template variables=SURFACE_SYNOP_FIELDS: True; template variables=SURFACE_SYNOP_INTERVAL: True; template variables=SURFACE_SYNOP_OFFSET: True; -[template variables=PLACEHOLDER_OBS] -ns=Observations/Types -description=Compare with observations from an undefined source. -help=This variable exists as a placeholder to show the intended - structure the GUI. Additional sources of observations should be - added using one window for each source. - -type=python_boolean -compulsory=true -sort-key=0yopp1a - ################################################################ [Observations/Types/Synoptic] diff --git a/src/CSET/cset_workflow/meta/rose-meta.conf.jinja2 b/src/CSET/cset_workflow/meta/rose-meta.conf.jinja2 index ec33e2e9f..1e486ab29 100644 --- a/src/CSET/cset_workflow/meta/rose-meta.conf.jinja2 +++ b/src/CSET/cset_workflow/meta/rose-meta.conf.jinja2 @@ -165,8 +165,8 @@ description=Geographical projection of plots help=This is is the projection used for the final plots, specified as a string. If no value is given, native latitudes and longitudes are used. compulsory=true -values="Unspecified", "NP_Stereo" -value-titles=Unspecified, North polar stereographic +values="", "NP_Stereo" +value-titles="", "North polar stereographic" sort-key=setup-h-out3a [template variables=PLOT_RESOLUTION] diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index fe6503214..edb970b43 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -1202,7 +1202,7 @@ def _plot_and_save_scattermap_plot( Mapping projection to be used by cartopy. """ # Setup plot details, size, resolution, etc. - fig = plt.figure(figsize=(15, 15), facecolor="w", edgecolor="k") + fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k") if projection is not None: # Apart from the default, the only projection we currently support is # a stereographic projection over the North Pole. @@ -1213,8 +1213,14 @@ def _plot_and_save_scattermap_plot( else: axes = plt.axes(projection=ccrs.PlateCarree()) - # Filled contour plot of the field. - mrk_size = int(np.sqrt(5000000.0 / len(cube.data))) + # Scatter plot of the field. The marker size is chosen to give + # symbols that decrease in size as the number of observations + # increases, although the fraction of the figure covered by + # symbols increases roughly as N^(1/2), disregarding overlaps, + # and has been selected for the default figure size of (10, 10). + # Should this be changed, the marker size should be adjusted in + # proportion to the area of the figure. + mrk_size = int(np.sqrt(2500000.0 / len(cube.data))) klon = None klat = None for kc in range(len(cube.aux_coords)): diff --git a/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml b/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml index 4c8ad1448..62a69affc 100644 --- a/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml +++ b/src/CSET/recipes/generic_model_obs_difference_scatterplot.yaml @@ -1,7 +1,8 @@ category: Quick Look title: Difference between modelled field of $VARNAME in model $MODEL_NAME and observation. description: | - Plot some observations. + Plot the difference between the model's value of the field $VARNAME and the + observed value at the points of observation. steps: - operator: read.read_cubes diff --git a/src/CSET/recipes/generic_obs_scatterplot.yaml b/src/CSET/recipes/generic_obs_scatterplot.yaml index 3bf32a73a..3e24fa5c2 100644 --- a/src/CSET/recipes/generic_obs_scatterplot.yaml +++ b/src/CSET/recipes/generic_obs_scatterplot.yaml @@ -1,7 +1,7 @@ category: Quick Look title: Observed values of $OBSVARNAME. description: | - Plot some observations. + Plot observations of OBSVARNAME as points at the locations of observation. steps: - operator: read.read_cubes