From 9fb1eaae30acb855fbd0eecfec5bd4881cd214a8 Mon Sep 17 00:00:00 2001 From: James Warner Date: Mon, 13 Jul 2026 09:47:44 +0100 Subject: [PATCH 01/56] placeholder in rose-meta --- src/CSET/cset_workflow/meta/observations/rose-meta.conf | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf index 1bb00e1254..e145787abf 100644 --- a/src/CSET/cset_workflow/meta/observations/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -110,3 +110,10 @@ description=Create Cardington diagnostic variable timeseries type=python_boolean compulsory=true sort-key=0card2 + + +################################################################ +# Model analysis/reanalysis +[ns=Observations/Analysis] +sort-key=sec-d5 +title=Model Analysis/Reanalysis From abdde1dbf99b915f6f333e9446e9d5b3c7c6f9c5 Mon Sep 17 00:00:00 2001 From: James Warner Date: Mon, 13 Jul 2026 10:08:05 +0100 Subject: [PATCH 02/56] add meta to rose edit --- .../meta/observations/rose-meta.conf | 44 +++++++++++++++++-- 1 file changed, 41 insertions(+), 3 deletions(-) diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf index e145787abf..bb02945227 100644 --- a/src/CSET/cset_workflow/meta/observations/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -111,9 +111,47 @@ type=python_boolean compulsory=true sort-key=0card2 - ################################################################ # Model analysis/reanalysis -[ns=Observations/Analysis] +[Observations/Analysis_and_Reanalysis] +ns=Observations/Analysis sort-key=sec-d5 -title=Model Analysis/Reanalysis +title=Model Analysis and Reanalysis + +[template variables=USE_REANALYSIS] +ns=Observations/Analysis +title=Use analysis/reanalysis in CSET +description=If True, then plot analysis/reanalysis where possible. +type=python_boolean +compulsory=true +sort-key=opta + +[template variables=REANALYSIS_NAME] +ns=Observations/Analysis +title=A name to use for the analysis/reanalysis. +description=A string to help identify the analysis i.e. 'ERA5' or 'BARPA' +type=quoted +compulsory=true +sort-key=optb + +[template variables=REANALYSIS_PATH] +ns=Observations/Analysis +title=Data path to analysis/reanalysis +description=A path to the analysis, using wildcards to parse time (see help) +help=Full path (including file name) to the analysis data on your chosen storage + system. Can contain wildcards. No quotation marks required in rose edit, as + it is automatically quoted there. + + strftime format strings are supported, and will be replaced with the + desired case study date or trial time. E.g: `/data/%Y%m%d/model1/*.nc` + + Commonly useful placeholders: + %Y: Year, e.g: 2024 + %m: Month, e.g: 12 + %d: Day, e.g: 31 + %H: Hour, e.g: 23 + %M: Minute, e.g: 59 + +type=quoted +compulsory=true +sort-key=optc From 8f110e883e7b8295264f4254a8e42ae36c4ed9ae Mon Sep 17 00:00:00 2001 From: James Warner Date: Mon, 13 Jul 2026 10:08:23 +0100 Subject: [PATCH 03/56] add template variables --- src/CSET/cset_workflow/rose-suite.conf.example | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/CSET/cset_workflow/rose-suite.conf.example b/src/CSET/cset_workflow/rose-suite.conf.example index af47f9a382..2741df07cb 100644 --- a/src/CSET/cset_workflow/rose-suite.conf.example +++ b/src/CSET/cset_workflow/rose-suite.conf.example @@ -117,6 +117,8 @@ PROFILE_PLEVEL_AGGREGATION=False,False,False,False RAIN_PRESENCE_DOMAIN_MEAN_TIMESERIES=False RAIN_PRESENCE_SPATIAL_DIFFERENCE=False RAIN_PRESENCE_SPATIAL_PLOT=False +REANALYSIS_NAME = "" +REANALYSIS_PATH = "" RUN_METPLUS_GRID_STAT=False RUN_METPLUS_POINT_STAT=False SCORES_ALL=False @@ -200,6 +202,7 @@ TIMESERIES_SURFACE_FIELD_AGGREGATION=False,False,False,False TIMESERIES_SURFACE_FIELD_LAND_MASK=False TIMESERIES_SURFACE_FIELD_SEA_MASK=False !!USE_WMO_STATION_NUMBERS=False +USE_REANALYSIS=False VIOLENT_RAIN_PRESENCE_DOMAIN_MEAN_TIMESERIES=False VIOLENT_RAIN_PRESENCE_SPATIAL_DIFFERENCE=False VIOLENT_RAIN_PRESENCE_SPATIAL_PLOT=False From 48939295941c47794c1d7d7adf703066f5bc994d Mon Sep 17 00:00:00 2001 From: James Warner Date: Tue, 14 Jul 2026 17:04:24 +0100 Subject: [PATCH 04/56] initial code lodge for tidying --- src/CSET/operators/collapse.py | 147 +++++++++++++++++++++++++++++++++ 1 file changed, 147 insertions(+) diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index dcd2994844..449deb1ce8 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -31,6 +31,150 @@ from CSET.operators.aggregate import add_hour_coordinate +def _coord_is_effectively_scalar(coord): + """True if coordinate is scalar or all points have the same value.""" + points = np.asarray(coord.points) + + if points.size <= 1: + return True + + return np.all(points == points.flat[0]) + + +def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): + + + analysis_cubes = iris.cube.CubeList() + forecast_cubes = iris.cube.CubeList() + + for cube in cubes: + fp = cube.coord("forecast_period") + + if fp.points.max() == 0: + analysis_cubes.append(cube) + else: + forecast_cubes.append(cube) + + + longest_forecast_period = max( + cube.coord("forecast_period").points.max() + for cube in forecast_cubes + ) + + frts = { + cube.coord("forecast_reference_time").points[0] + for cube in forecast_cubes + } + + if len(frts) > 1: + raise ValueError( + "Forecast cubes have different forecast_reference_times" + ) + + frt_coord = forecast_cubes[0].coord("forecast_reference_time") + + forecast_reference_time = frt_coord.units.num2date( + frt_coord.points[0] +) + + longest_cube = max( + forecast_cubes, + key=lambda c: c.coord("forecast_period").points.max() + ) + + forecast_reference_time = ( + longest_cube.coord("forecast_reference_time") + .units.num2date( + longest_cube.coord("forecast_reference_time").points[0] + ) + ) + + forecast_end_time = ( + longest_cube.coord("time") + .units.num2date( + longest_cube.coord("time").points.max() + ) + ) + + time_constraint = iris.Constraint( + time=lambda cell: ( + forecast_reference_time <= cell.point <= forecast_end_time + ) + ) + + # manually slicing needed rather than just assigning cube = cube.extract... + for i, cube in enumerate(analysis_cubes): + analysis_cubes[i] = cube.extract(time_constraint) + + import cf_units + + fixed_cubes = iris.cube.CubeList() + + for cube in forecast_cubes: + fixed_cubes.append(cube) + + print(analysis_cubes) + print(analysis_cubes[0]) + print(analysis_cubes[0].coord('time')) + print(analysis_cubes[0].coord('forecast_reference_time')) + print(analysis_cubes[0].coord('forecast_period')) + print('pre') + + for cube in analysis_cubes: + # Get coords + time_coord = cube.coord("time") + frt_coord = cube.coord("forecast_reference_time") + fp_coord = cube.coord("forecast_period") + + + if (_coord_is_effectively_scalar(frt_coord) + and _coord_is_effectively_scalar(fp_coord)): + print("forecast_reference_time and forecast_period already static") + return cube + + + # + frt0_hours = frt_coord.points[0] + + fp_hours = time_coord.points - frt0_hours + + + # Remove old coords + cube.remove_coord("forecast_reference_time") + cube.remove_coord("forecast_period") + + # Add constant FRT back as scalar + new_frt = iris.coords.AuxCoord( + frt0_hours, + standard_name="forecast_reference_time", + units=time_coord.units, + ) + cube.add_aux_coord(new_frt) + + # Add forecast period varying along time dimension + time_dim = cube.coord_dims("time")[0] + + new_fp = iris.coords.AuxCoord( + fp_hours, + standard_name="forecast_period", + units="hours", + ) + + cube.add_aux_coord(new_fp, data_dims=(time_dim,)) + + fixed_cubes.append(cube) + + print(fixed_cubes) + print(fixed_cubes[0]) + print(fixed_cubes[0].coord('time')) + print(fixed_cubes[0].coord('forecast_reference_time')) + print(fixed_cubes[0].coord('forecast_period')) + + #something to check we are matching input, i.e. if cube, only return cube, not cubelist. + + return fixed_cubes + + def collapse( cubes: iris.cube.Cube | iris.cube.CubeList, coordinate: str | list[str], @@ -75,6 +219,9 @@ def collapse( raise ValueError("Must specify additional_percent") # Retain only common time points between different models if multiple model inputs. + + cubes = _fix_analysis_forecasttime(cubes) + if isinstance(cubes, iris.cube.CubeList) and len(cubes) > 1: logging.debug( "Extracting common time points as multiple model inputs detected." From b3aa6fe11cd5dbf9e88b7f91506dde966cb61e38 Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 15 Jul 2026 09:02:15 +0100 Subject: [PATCH 05/56] add notes to pickup --- src/CSET/operators/collapse.py | 77 +++++++++++++--------------------- 1 file changed, 29 insertions(+), 48 deletions(-) diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index 449deb1ce8..5b818e58a7 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -43,7 +43,6 @@ def _coord_is_effectively_scalar(coord): def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): - analysis_cubes = iris.cube.CubeList() forecast_cubes = iris.cube.CubeList() @@ -55,59 +54,43 @@ def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): else: forecast_cubes.append(cube) - longest_forecast_period = max( - cube.coord("forecast_period").points.max() - for cube in forecast_cubes + cube.coord("forecast_period").points.max() for cube in forecast_cubes ) - frts = { - cube.coord("forecast_reference_time").points[0] - for cube in forecast_cubes - } + frts = {cube.coord("forecast_reference_time").points[0] for cube in forecast_cubes} + # We might have multiple cubes, with different reference times, + # so we need to go through all possibilities and ensure analysis + # matches all of these. Maybe a list of forecast reference times and leadtimes + # to construct the appropriate cubes? if len(frts) > 1: - raise ValueError( - "Forecast cubes have different forecast_reference_times" - ) - + raise ValueError("Forecast cubes have different forecast_reference_times") + frt_coord = forecast_cubes[0].coord("forecast_reference_time") - forecast_reference_time = frt_coord.units.num2date( - frt_coord.points[0] -) + forecast_reference_time = frt_coord.units.num2date(frt_coord.points[0]) longest_cube = max( - forecast_cubes, - key=lambda c: c.coord("forecast_period").points.max() + forecast_cubes, key=lambda c: c.coord("forecast_period").points.max() ) - forecast_reference_time = ( - longest_cube.coord("forecast_reference_time") - .units.num2date( - longest_cube.coord("forecast_reference_time").points[0] - ) - ) + forecast_reference_time = longest_cube.coord( + "forecast_reference_time" + ).units.num2date(longest_cube.coord("forecast_reference_time").points[0]) - forecast_end_time = ( - longest_cube.coord("time") - .units.num2date( - longest_cube.coord("time").points.max() - ) + forecast_end_time = longest_cube.coord("time").units.num2date( + longest_cube.coord("time").points.max() ) time_constraint = iris.Constraint( - time=lambda cell: ( - forecast_reference_time <= cell.point <= forecast_end_time - ) + time=lambda cell: forecast_reference_time <= cell.point <= forecast_end_time ) # manually slicing needed rather than just assigning cube = cube.extract... for i, cube in enumerate(analysis_cubes): analysis_cubes[i] = cube.extract(time_constraint) - import cf_units - fixed_cubes = iris.cube.CubeList() for cube in forecast_cubes: @@ -115,10 +98,10 @@ def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): print(analysis_cubes) print(analysis_cubes[0]) - print(analysis_cubes[0].coord('time')) - print(analysis_cubes[0].coord('forecast_reference_time')) - print(analysis_cubes[0].coord('forecast_period')) - print('pre') + print(analysis_cubes[0].coord("time")) + print(analysis_cubes[0].coord("forecast_reference_time")) + print(analysis_cubes[0].coord("forecast_period")) + print("pre") for cube in analysis_cubes: # Get coords @@ -126,19 +109,17 @@ def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): frt_coord = cube.coord("forecast_reference_time") fp_coord = cube.coord("forecast_period") - - if (_coord_is_effectively_scalar(frt_coord) - and _coord_is_effectively_scalar(fp_coord)): + if _coord_is_effectively_scalar(frt_coord) and _coord_is_effectively_scalar( + fp_coord + ): print("forecast_reference_time and forecast_period already static") return cube - - # - frt0_hours = frt_coord.points[0] + # + frt0_hours = frt_coord.points[0] fp_hours = time_coord.points - frt0_hours - # Remove old coords cube.remove_coord("forecast_reference_time") cube.remove_coord("forecast_period") @@ -166,11 +147,11 @@ def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): print(fixed_cubes) print(fixed_cubes[0]) - print(fixed_cubes[0].coord('time')) - print(fixed_cubes[0].coord('forecast_reference_time')) - print(fixed_cubes[0].coord('forecast_period')) + print(fixed_cubes[0].coord("time")) + print(fixed_cubes[0].coord("forecast_reference_time")) + print(fixed_cubes[0].coord("forecast_period")) - #something to check we are matching input, i.e. if cube, only return cube, not cubelist. + # something to check we are matching input, i.e. if cube, only return cube, not cubelist. return fixed_cubes From b9f652ec1e6e616c7e5272e7e93b638d0371a143 Mon Sep 17 00:00:00 2001 From: James Warner Date: Thu, 16 Jul 2026 13:46:01 +0100 Subject: [PATCH 06/56] latest changes --- src/CSET/operators/collapse.py | 234 +++++++++++++++++++++++---------- 1 file changed, 163 insertions(+), 71 deletions(-) diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index 5b818e58a7..aa45e87c40 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -26,6 +26,7 @@ import iris.exceptions import iris.util import numpy as np +import datetime as dt from CSET._common import iter_maybe from CSET.operators.aggregate import add_hour_coordinate @@ -43,6 +44,7 @@ def _coord_is_effectively_scalar(coord): def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): + # Partition out forecast and analysis cubes based on forecast_period. analysis_cubes = iris.cube.CubeList() forecast_cubes = iris.cube.CubeList() @@ -54,106 +56,196 @@ def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): else: forecast_cubes.append(cube) - longest_forecast_period = max( - cube.coord("forecast_period").points.max() for cube in forecast_cubes - ) - - frts = {cube.coord("forecast_reference_time").points[0] for cube in forecast_cubes} + if len(analysis_cubes) > 1: + raise ValueError + if len(forecast_cubes) > 1: + raise ValueError + + analysis_cube = analysis_cubes[0] + forecast_cube = forecast_cubes[0] - # We might have multiple cubes, with different reference times, - # so we need to go through all possibilities and ensure analysis - # matches all of these. Maybe a list of forecast reference times and leadtimes - # to construct the appropriate cubes? - if len(frts) > 1: - raise ValueError("Forecast cubes have different forecast_reference_times") + # Analysis cube will think it has multiple realizations cause it duplicates + analysis_cube = analysis_cube.extract(iris.Constraint(realization = analysis_cube.coord('realization').points[0])) - frt_coord = forecast_cubes[0].coord("forecast_reference_time") + # Get forecast_reference_time + frt_coord = forecast_cube.coord("forecast_reference_time") - forecast_reference_time = frt_coord.units.num2date(frt_coord.points[0]) + # Get forecast length + fc_length = forecast_cube.coord("forecast_period").points[-1] - longest_cube = max( - forecast_cubes, key=lambda c: c.coord("forecast_period").points.max() - ) + # Work out minimum, maximum valid times in analysis + an_min = analysis_cube.coord("time").units.num2date(analysis_cube.coord("time").points[0]) + an_max = analysis_cube.coord("time").units.num2date(analysis_cube.coord("time").points[-1]) - forecast_reference_time = longest_cube.coord( - "forecast_reference_time" - ).units.num2date(longest_cube.coord("forecast_reference_time").points[0]) + fixed_analysis_cubes = iris.cube.CubeList() - forecast_end_time = longest_cube.coord("time").units.num2date( - longest_cube.coord("time").points.max() - ) + for t in frt_coord.points: + time = frt_coord.units.num2date(t) - time_constraint = iris.Constraint( - time=lambda cell: forecast_reference_time <= cell.point <= forecast_end_time + if time < an_min or (time+dt.timedelta(hours=fc_length)) > an_max: + print("analysis time range ",an_min, an_max) + print("forecast spans",time, time+dt.timedelta(hours=fc_length)) + print('analysis does not cover full forecast') + + time_constraint = iris.Constraint( + time=lambda cell: time <= cell.point <= time+dt.timedelta(hours=fc_length) ) - # manually slicing needed rather than just assigning cube = cube.extract... - for i, cube in enumerate(analysis_cubes): - analysis_cubes[i] = cube.extract(time_constraint) - - fixed_cubes = iris.cube.CubeList() - - for cube in forecast_cubes: - fixed_cubes.append(cube) - - print(analysis_cubes) - print(analysis_cubes[0]) - print(analysis_cubes[0].coord("time")) - print(analysis_cubes[0].coord("forecast_reference_time")) - print(analysis_cubes[0].coord("forecast_period")) - print("pre") - - for cube in analysis_cubes: - # Get coords - time_coord = cube.coord("time") - frt_coord = cube.coord("forecast_reference_time") - fp_coord = cube.coord("forecast_period") - - if _coord_is_effectively_scalar(frt_coord) and _coord_is_effectively_scalar( - fp_coord - ): - print("forecast_reference_time and forecast_period already static") - return cube - - # - frt0_hours = frt_coord.points[0] + analysis_slice = analysis_cube.extract(time_constraint) - fp_hours = time_coord.points - frt0_hours + print(analysis_slice) - # Remove old coords - cube.remove_coord("forecast_reference_time") - cube.remove_coord("forecast_period") + # Remove old coords + analysis_slice.remove_coord("forecast_reference_time") + analysis_slice.remove_coord("forecast_period") # Add constant FRT back as scalar new_frt = iris.coords.AuxCoord( - frt0_hours, + t, standard_name="forecast_reference_time", - units=time_coord.units, + units=frt_coord.units, ) - cube.add_aux_coord(new_frt) + analysis_slice.add_aux_coord(new_frt) # Add forecast period varying along time dimension time_dim = cube.coord_dims("time")[0] new_fp = iris.coords.AuxCoord( - fp_hours, + analysis_slice.coord('time').points - analysis_slice.coord('time').points[0], standard_name="forecast_period", units="hours", ) - cube.add_aux_coord(new_fp, data_dims=(time_dim,)) + print(time_dim) #producing value of 1, yet time dim should be zero, hence failing below + # when trying to link forecast period to the latitude dimension rather than time. + quit() + + analysis_slice.add_aux_coord(new_fp, data_dims=(time_dim,)) + + fixed_analysis_cubes.append(analysis_slice) + + fixed_analysis_cubes = fixed_analysis_cubes.merge() + + print(fixed_analysis_cubes) + for cube in fixed_analysis_cubes: + print(cube) + quit() + + # frt_coord = forecast_cubes[0].coord("forecast_reference_time") + + # forecast_reference_time = frt_coord.units.num2date(frt_coord.points[0]) + + # # Get all forecast reference times from forecast_cubes + # frt_coord = forecast_cubes.coord("forecast_reference_time").points + + + + + + # anaysis will be double loaded and givena realisation time if aggregation, so + # we'll need to extract the zero'th member. Need to check if that is the case if + # no aggregation. + + # longest_forecast_period = max( + # cube.coord("forecast_period").points.max() for cube in forecast_cubes + # ) + + # frts = {cube.coord("forecast_reference_time").points[0] for cube in forecast_cubes} + + # # We might have multiple cubes, with different reference times, + # # so we need to go through all possibilities and ensure analysis + # # matches all of these. Maybe a list of forecast reference times and leadtimes + # # to construct the appropriate cubes? + # if len(frts) > 1: + # raise ValueError("Forecast cubes have different forecast_reference_times") + + # frt_coord = forecast_cubes[0].coord("forecast_reference_time") + + # forecast_reference_time = frt_coord.units.num2date(frt_coord.points[0]) + + # longest_cube = max( + # forecast_cubes, key=lambda c: c.coord("forecast_period").points.max() + # ) + + # forecast_reference_time = longest_cube.coord( + # "forecast_reference_time" + # ).units.num2date(longest_cube.coord("forecast_reference_time").points[0]) + + # forecast_end_time = longest_cube.coord("time").units.num2date( + # longest_cube.coord("time").points.max() + # ) + + # time_constraint = iris.Constraint( + # time=lambda cell: forecast_reference_time <= cell.point <= forecast_end_time + # ) + + # # manually slicing needed rather than just assigning cube = cube.extract... + # for i, cube in enumerate(analysis_cubes): + # analysis_cubes[i] = cube.extract(time_constraint) + + # fixed_cubes = iris.cube.CubeList() + + # for cube in forecast_cubes: + # fixed_cubes.append(cube) + + # # print(analysis_cubes) + # # print(analysis_cubes[0]) + # # print(analysis_cubes[0].coord("time")) + # # print(analysis_cubes[0].coord("forecast_reference_time")) + # # print(analysis_cubes[0].coord("forecast_period")) + # # print("pre") + + # for cube in analysis_cubes: + # # Get coords + # time_coord = cube.coord("time") + # frt_coord = cube.coord("forecast_reference_time") + # fp_coord = cube.coord("forecast_period") + + # if _coord_is_effectively_scalar(frt_coord) and _coord_is_effectively_scalar( + # fp_coord + # ): + # print("forecast_reference_time and forecast_period already static") + # return cube + + # # + # frt0_hours = frt_coord.points[0] + + # fp_hours = time_coord.points - frt0_hours + + # # Remove old coords + # cube.remove_coord("forecast_reference_time") + # cube.remove_coord("forecast_period") + + # # Add constant FRT back as scalar + # new_frt = iris.coords.AuxCoord( + # frt0_hours, + # standard_name="forecast_reference_time", + # units=time_coord.units, + # ) + # cube.add_aux_coord(new_frt) + + # # Add forecast period varying along time dimension + # time_dim = cube.coord_dims("time")[0] + + # new_fp = iris.coords.AuxCoord( + # fp_hours, + # standard_name="forecast_period", + # units="hours", + # ) + + # cube.add_aux_coord(new_fp, data_dims=(time_dim,)) - fixed_cubes.append(cube) + # fixed_cubes.append(cube) - print(fixed_cubes) - print(fixed_cubes[0]) - print(fixed_cubes[0].coord("time")) - print(fixed_cubes[0].coord("forecast_reference_time")) - print(fixed_cubes[0].coord("forecast_period")) + # # print(fixed_cubes) + # # print(fixed_cubes[0]) + # # print(fixed_cubes[0].coord("time")) + # # print(fixed_cubes[0].coord("forecast_reference_time")) + # # print(fixed_cubes[0].coord("forecast_period")) - # something to check we are matching input, i.e. if cube, only return cube, not cubelist. + # # something to check we are matching input, i.e. if cube, only return cube, not cubelist. - return fixed_cubes + # return fixed_cubes def collapse( From ca8033c4cf168454d29c9e66560499dc9814d1f2 Mon Sep 17 00:00:00 2001 From: James Warner Date: Thu, 16 Jul 2026 15:16:30 +0100 Subject: [PATCH 07/56] lodge latest code --- src/CSET/operators/collapse.py | 190 +++++++++++++++++++++++++++------ 1 file changed, 157 insertions(+), 33 deletions(-) diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index aa45e87c40..c872da53ae 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -44,6 +44,10 @@ def _coord_is_effectively_scalar(coord): def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): + """ + TODO - in progress + """ + # Partition out forecast and analysis cubes based on forecast_period. analysis_cubes = iris.cube.CubeList() forecast_cubes = iris.cube.CubeList() @@ -59,6 +63,7 @@ def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): if len(analysis_cubes) > 1: raise ValueError if len(forecast_cubes) > 1: + print(forecast_cubes) raise ValueError analysis_cube = analysis_cubes[0] @@ -67,11 +72,12 @@ def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): # Analysis cube will think it has multiple realizations cause it duplicates analysis_cube = analysis_cube.extract(iris.Constraint(realization = analysis_cube.coord('realization').points[0])) - # Get forecast_reference_time + # Get forecast_reference_time and period information frt_coord = forecast_cube.coord("forecast_reference_time") - # Get forecast length - fc_length = forecast_cube.coord("forecast_period").points[-1] + fc_period_coord = forecast_cube.coord("forecast_period") + fc_periods = fc_period_coord.points + fc_length = fc_periods[-1] # Work out minimum, maximum valid times in analysis an_min = analysis_cube.coord("time").units.num2date(analysis_cube.coord("time").points[0]) @@ -79,59 +85,177 @@ def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): fixed_analysis_cubes = iris.cube.CubeList() + # Iterate over all forecast reference times for t in frt_coord.points: - time = frt_coord.units.num2date(t) + frt = frt_coord.units.num2date(t) + + if frt < an_min or ( + frt + dt.timedelta(hours=float(fc_length)) + ) > an_max: - if time < an_min or (time+dt.timedelta(hours=fc_length)) > an_max: - print("analysis time range ",an_min, an_max) - print("forecast spans",time, time+dt.timedelta(hours=fc_length)) - print('analysis does not cover full forecast') + raise ValueError( + f"Analysis does not cover forecast span " + f"{frt} -> " + f"{frt + dt.timedelta(hours=float(fc_length))}") + time_constraint = iris.Constraint( - time=lambda cell: time <= cell.point <= time+dt.timedelta(hours=fc_length) + time=lambda cell: frt <= cell.point <= frt+dt.timedelta(hours=fc_length) ) analysis_slice = analysis_cube.extract(time_constraint) - print(analysis_slice) - # Remove old coords analysis_slice.remove_coord("forecast_reference_time") analysis_slice.remove_coord("forecast_period") - # Add constant FRT back as scalar - new_frt = iris.coords.AuxCoord( - t, - standard_name="forecast_reference_time", + + time_coord = analysis_slice.coord("time") + time_dim = analysis_slice.coord_dims("time")[0] + + # Forecast periods relative to this FRT. + # + fp_points = ( + time_coord.points - time_coord.points[0] + ) + + fp_coord = iris.coords.DimCoord( + fp_points, + standard_name="forecast_period", + units='hours', + ) + + # Replace time dimension coord with forecast_period. + # + analysis_slice.remove_coord("time") + analysis_slice.add_dim_coord(fp_coord, time_dim) + + # Add FRT scalar coordinate. + # + analysis_slice.add_aux_coord( + iris.coords.AuxCoord( + t, + standard_name="forecast_reference_time", + units=frt_coord.units, + ) + ) + + fixed_analysis_cubes.append( + analysis_slice + ) + + if len(fixed_analysis_cubes) == 1: + + cube = fixed_analysis_cubes[0] + + fp_coord = cube.coord("forecast_period") + + time_points = ( + cube.coord( + "forecast_reference_time" + ).points[0] + + fp_coord.points + ) + + time_coord = iris.coords.AuxCoord( + time_points, + standard_name="time", units=frt_coord.units, ) - analysis_slice.add_aux_coord(new_frt) - # Add forecast period varying along time dimension - time_dim = cube.coord_dims("time")[0] + cube.add_aux_coord( + time_coord, + data_dims=( + cube.coord_dims( + "forecast_period" + )[0], + ), + ) - new_fp = iris.coords.AuxCoord( - analysis_slice.coord('time').points - analysis_slice.coord('time').points[0], - standard_name="forecast_period", - units="hours", + else: + # multiple FRTs + # + cube = fixed_analysis_cubes.merge_cube() + + fp_dim = cube.coord_dims( + "forecast_period" + )[0] + + frt_dim = cube.coord_dims( + "forecast_reference_time" + )[0] + + fp_points = cube.coord( + "forecast_period" + ).points + + frt_points = cube.coord( + "forecast_reference_time" + ).points + + time_points = ( + fp_points[:, None] + + frt_points[None, :] + ) + + time_coord = iris.coords.AuxCoord( + time_points, + standard_name="time", + units=frt_coord.units, + ) + + cube.add_aux_coord( + time_coord, + data_dims=(fp_dim, frt_dim), ) - print(time_dim) #producing value of 1, yet time dim should be zero, hence failing below - # when trying to link forecast period to the latitude dimension rather than time. - quit() + fixed_analysis_cubes = iris.cube.CubeList() + fixed_analysis_cubes.append(cube) - analysis_slice.add_aux_coord(new_fp, data_dims=(time_dim,)) + return iris.cube.CubeList([cube, forecast_cube]) - fixed_analysis_cubes.append(analysis_slice) - fixed_analysis_cubes = fixed_analysis_cubes.merge() - print(fixed_analysis_cubes) - for cube in fixed_analysis_cubes: - print(cube) - quit() - # frt_coord = forecast_cubes[0].coord("forecast_reference_time") +# # Keep valid time as auxiliary coordinate. +# # +# analysis_slice.add_aux_coord(time_coord, data_dims=(time_dim,)) + +# # + + +# fixed_analysis_cubes.append(analysis_slice) + +# if len(fixed_analysis_cubes) > 1: + + + + + + + + + + +# for cube in fixed_analysis_cubes: +# print(cube) + + +# print(fixed_analysis_cubes.concatenate()) +# print(fixed_analysis_cubes.merge()) + + +# try: +# merged = fixed_analysis_cubes.merge_cube() +# except Exception as e: +# print(type(e)) +# print(e) + + +# # print(fixed_analysis_cubes) +# quit() + +# # frt_coord = forecast_cubes[0].coord("forecast_reference_time") # forecast_reference_time = frt_coord.units.num2date(frt_coord.points[0]) From ae02f932ab614f8e592a6007627147c2eef57a2b Mon Sep 17 00:00:00 2001 From: James Warner Date: Thu, 16 Jul 2026 15:17:02 +0100 Subject: [PATCH 08/56] clear out old code --- src/CSET/operators/collapse.py | 157 --------------------------------- 1 file changed, 157 deletions(-) diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index c872da53ae..db933c7731 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -215,163 +215,6 @@ def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): return iris.cube.CubeList([cube, forecast_cube]) - - -# # Keep valid time as auxiliary coordinate. -# # -# analysis_slice.add_aux_coord(time_coord, data_dims=(time_dim,)) - -# # - - -# fixed_analysis_cubes.append(analysis_slice) - -# if len(fixed_analysis_cubes) > 1: - - - - - - - - - - -# for cube in fixed_analysis_cubes: -# print(cube) - - -# print(fixed_analysis_cubes.concatenate()) -# print(fixed_analysis_cubes.merge()) - - -# try: -# merged = fixed_analysis_cubes.merge_cube() -# except Exception as e: -# print(type(e)) -# print(e) - - -# # print(fixed_analysis_cubes) -# quit() - -# # frt_coord = forecast_cubes[0].coord("forecast_reference_time") - - # forecast_reference_time = frt_coord.units.num2date(frt_coord.points[0]) - - # # Get all forecast reference times from forecast_cubes - # frt_coord = forecast_cubes.coord("forecast_reference_time").points - - - - - - # anaysis will be double loaded and givena realisation time if aggregation, so - # we'll need to extract the zero'th member. Need to check if that is the case if - # no aggregation. - - # longest_forecast_period = max( - # cube.coord("forecast_period").points.max() for cube in forecast_cubes - # ) - - # frts = {cube.coord("forecast_reference_time").points[0] for cube in forecast_cubes} - - # # We might have multiple cubes, with different reference times, - # # so we need to go through all possibilities and ensure analysis - # # matches all of these. Maybe a list of forecast reference times and leadtimes - # # to construct the appropriate cubes? - # if len(frts) > 1: - # raise ValueError("Forecast cubes have different forecast_reference_times") - - # frt_coord = forecast_cubes[0].coord("forecast_reference_time") - - # forecast_reference_time = frt_coord.units.num2date(frt_coord.points[0]) - - # longest_cube = max( - # forecast_cubes, key=lambda c: c.coord("forecast_period").points.max() - # ) - - # forecast_reference_time = longest_cube.coord( - # "forecast_reference_time" - # ).units.num2date(longest_cube.coord("forecast_reference_time").points[0]) - - # forecast_end_time = longest_cube.coord("time").units.num2date( - # longest_cube.coord("time").points.max() - # ) - - # time_constraint = iris.Constraint( - # time=lambda cell: forecast_reference_time <= cell.point <= forecast_end_time - # ) - - # # manually slicing needed rather than just assigning cube = cube.extract... - # for i, cube in enumerate(analysis_cubes): - # analysis_cubes[i] = cube.extract(time_constraint) - - # fixed_cubes = iris.cube.CubeList() - - # for cube in forecast_cubes: - # fixed_cubes.append(cube) - - # # print(analysis_cubes) - # # print(analysis_cubes[0]) - # # print(analysis_cubes[0].coord("time")) - # # print(analysis_cubes[0].coord("forecast_reference_time")) - # # print(analysis_cubes[0].coord("forecast_period")) - # # print("pre") - - # for cube in analysis_cubes: - # # Get coords - # time_coord = cube.coord("time") - # frt_coord = cube.coord("forecast_reference_time") - # fp_coord = cube.coord("forecast_period") - - # if _coord_is_effectively_scalar(frt_coord) and _coord_is_effectively_scalar( - # fp_coord - # ): - # print("forecast_reference_time and forecast_period already static") - # return cube - - # # - # frt0_hours = frt_coord.points[0] - - # fp_hours = time_coord.points - frt0_hours - - # # Remove old coords - # cube.remove_coord("forecast_reference_time") - # cube.remove_coord("forecast_period") - - # # Add constant FRT back as scalar - # new_frt = iris.coords.AuxCoord( - # frt0_hours, - # standard_name="forecast_reference_time", - # units=time_coord.units, - # ) - # cube.add_aux_coord(new_frt) - - # # Add forecast period varying along time dimension - # time_dim = cube.coord_dims("time")[0] - - # new_fp = iris.coords.AuxCoord( - # fp_hours, - # standard_name="forecast_period", - # units="hours", - # ) - - # cube.add_aux_coord(new_fp, data_dims=(time_dim,)) - - # fixed_cubes.append(cube) - - # # print(fixed_cubes) - # # print(fixed_cubes[0]) - # # print(fixed_cubes[0].coord("time")) - # # print(fixed_cubes[0].coord("forecast_reference_time")) - # # print(fixed_cubes[0].coord("forecast_period")) - - # # something to check we are matching input, i.e. if cube, only return cube, not cubelist. - - # return fixed_cubes - - def collapse( cubes: iris.cube.Cube | iris.cube.CubeList, coordinate: str | list[str], From 6c3ed5ea2b40fb4087197450e05bd81156ad547e Mon Sep 17 00:00:00 2001 From: James Warner Date: Thu, 16 Jul 2026 15:51:23 +0100 Subject: [PATCH 09/56] return if no analysis --- src/CSET/operators/collapse.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index db933c7731..759ce6c7d1 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -65,6 +65,9 @@ def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): if len(forecast_cubes) > 1: print(forecast_cubes) raise ValueError + if len(analysis_cubes) == 0: + return cubes #i.e. no reanalysis so ignore this function + analysis_cube = analysis_cubes[0] forecast_cube = forecast_cubes[0] From c04e46ff15025d78431dc5a554e83967e5d38c6b Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 29 Jul 2026 09:33:39 +0100 Subject: [PATCH 10/56] back out rose meta changes --- .../meta/observations/rose-meta.conf | 70 ------------------- 1 file changed, 70 deletions(-) diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf index 6947167d74..8c7bc0d139 100644 --- a/src/CSET/cset_workflow/meta/observations/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -257,73 +257,3 @@ help=Switch to select gridded UK Nimrod 5 minute rainfall rate. compulsory=true type=python_boolean sort-key=nimrod-radar-7 - -################################################################ -# Cardington UK Research Observations -[Observations/Cardington] -ns=Observations/Cardington -sort-key=sec-d4 -title=Cardington UK observations - -# Air temperature timeseries -[template variables=CARDINGTON_AIR_TEMPERATURE_SINGLE_POINT_TIME_SERIES] -ns=Observations/Cardington -title=Cardington air temperature timeseries -description=Create Cardington diagnostic variable timeseries -type=python_boolean -compulsory=true -sort-key=0card1 - -# Relative humidity timeseries -[template variables=CARDINGTON_RELATIVE_HUMIDITY_SINGLE_POINT_TIME_SERIES] -ns=Observations/Cardington -title=Cardington relative humidity timeseries -description=Create Cardington diagnostic variable timeseries -type=python_boolean -compulsory=true -sort-key=0card2 - -################################################################ -# Model analysis/reanalysis -[Observations/Analysis_and_Reanalysis] -ns=Observations/Analysis -sort-key=sec-d5 -title=Model Analysis and Reanalysis - -[template variables=USE_REANALYSIS] -ns=Observations/Analysis -title=Use analysis/reanalysis in CSET -description=If True, then plot analysis/reanalysis where possible. -type=python_boolean -compulsory=true -sort-key=opta - -[template variables=REANALYSIS_NAME] -ns=Observations/Analysis -title=A name to use for the analysis/reanalysis. -description=A string to help identify the analysis i.e. 'ERA5' or 'BARPA' -type=quoted -compulsory=true -sort-key=optb - -[template variables=REANALYSIS_PATH] -ns=Observations/Analysis -title=Data path to analysis/reanalysis -description=A path to the analysis, using wildcards to parse time (see help) -help=Full path (including file name) to the analysis data on your chosen storage - system. Can contain wildcards. No quotation marks required in rose edit, as - it is automatically quoted there. - - strftime format strings are supported, and will be replaced with the - desired case study date or trial time. E.g: `/data/%Y%m%d/model1/*.nc` - - Commonly useful placeholders: - %Y: Year, e.g: 2024 - %m: Month, e.g: 12 - %d: Day, e.g: 31 - %H: Hour, e.g: 23 - %M: Minute, e.g: 59 - -type=quoted -compulsory=true -sort-key=optc From d16ae1233a163c02146c485bb5892ec633b80225 Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 29 Jul 2026 09:34:11 +0100 Subject: [PATCH 11/56] remove template variables for reanalysis --- src/CSET/cset_workflow/rose-suite.conf.example | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/CSET/cset_workflow/rose-suite.conf.example b/src/CSET/cset_workflow/rose-suite.conf.example index 71644aef7f..a43f7a5f63 100644 --- a/src/CSET/cset_workflow/rose-suite.conf.example +++ b/src/CSET/cset_workflow/rose-suite.conf.example @@ -142,8 +142,6 @@ PROFILE_PLEVEL_AGGREGATION=False,False,False,False RAIN_PRESENCE_DOMAIN_MEAN_TIMESERIES=False RAIN_PRESENCE_SPATIAL_DIFFERENCE=False RAIN_PRESENCE_SPATIAL_PLOT=False -REANALYSIS_NAME = "" -REANALYSIS_PATH = "" RUN_METPLUS_ENSEMBLE_STAT=False RUN_METPLUS_GRID_STAT=False RUN_METPLUS_POINT_STAT=False @@ -230,8 +228,6 @@ TIMESERIES_SURFACE_FIELD=False TIMESERIES_SURFACE_FIELD_AGGREGATION=False,False,False,False TIMESERIES_SURFACE_FIELD_LAND_MASK=False TIMESERIES_SURFACE_FIELD_SEA_MASK=False -!!USE_WMO_STATION_NUMBERS=False -USE_REANALYSIS=False VIOLENT_RAIN_PRESENCE_DOMAIN_MEAN_TIMESERIES=False VIOLENT_RAIN_PRESENCE_SPATIAL_DIFFERENCE=False VIOLENT_RAIN_PRESENCE_SPATIAL_PLOT=False From 12868ad25b9b5bb9be95066860afb5b5e998747b Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 29 Jul 2026 09:34:44 +0100 Subject: [PATCH 12/56] restore changes in collapse backing out code --- src/CSET/operators/collapse.py | 190 --------------------------------- 1 file changed, 190 deletions(-) diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index 759ce6c7d1..dcd2994844 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -26,198 +26,11 @@ import iris.exceptions import iris.util import numpy as np -import datetime as dt from CSET._common import iter_maybe from CSET.operators.aggregate import add_hour_coordinate -def _coord_is_effectively_scalar(coord): - """True if coordinate is scalar or all points have the same value.""" - points = np.asarray(coord.points) - - if points.size <= 1: - return True - - return np.all(points == points.flat[0]) - - -def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): - - """ - TODO - in progress - """ - - # Partition out forecast and analysis cubes based on forecast_period. - analysis_cubes = iris.cube.CubeList() - forecast_cubes = iris.cube.CubeList() - - for cube in cubes: - fp = cube.coord("forecast_period") - - if fp.points.max() == 0: - analysis_cubes.append(cube) - else: - forecast_cubes.append(cube) - - if len(analysis_cubes) > 1: - raise ValueError - if len(forecast_cubes) > 1: - print(forecast_cubes) - raise ValueError - if len(analysis_cubes) == 0: - return cubes #i.e. no reanalysis so ignore this function - - - analysis_cube = analysis_cubes[0] - forecast_cube = forecast_cubes[0] - - # Analysis cube will think it has multiple realizations cause it duplicates - analysis_cube = analysis_cube.extract(iris.Constraint(realization = analysis_cube.coord('realization').points[0])) - - # Get forecast_reference_time and period information - frt_coord = forecast_cube.coord("forecast_reference_time") - - fc_period_coord = forecast_cube.coord("forecast_period") - fc_periods = fc_period_coord.points - fc_length = fc_periods[-1] - - # Work out minimum, maximum valid times in analysis - an_min = analysis_cube.coord("time").units.num2date(analysis_cube.coord("time").points[0]) - an_max = analysis_cube.coord("time").units.num2date(analysis_cube.coord("time").points[-1]) - - fixed_analysis_cubes = iris.cube.CubeList() - - # Iterate over all forecast reference times - for t in frt_coord.points: - frt = frt_coord.units.num2date(t) - - if frt < an_min or ( - frt + dt.timedelta(hours=float(fc_length)) - ) > an_max: - - raise ValueError( - f"Analysis does not cover forecast span " - f"{frt} -> " - f"{frt + dt.timedelta(hours=float(fc_length))}") - - - time_constraint = iris.Constraint( - time=lambda cell: frt <= cell.point <= frt+dt.timedelta(hours=fc_length) - ) - - analysis_slice = analysis_cube.extract(time_constraint) - - # Remove old coords - analysis_slice.remove_coord("forecast_reference_time") - analysis_slice.remove_coord("forecast_period") - - - time_coord = analysis_slice.coord("time") - time_dim = analysis_slice.coord_dims("time")[0] - - # Forecast periods relative to this FRT. - # - fp_points = ( - time_coord.points - time_coord.points[0] - ) - - fp_coord = iris.coords.DimCoord( - fp_points, - standard_name="forecast_period", - units='hours', - ) - - # Replace time dimension coord with forecast_period. - # - analysis_slice.remove_coord("time") - analysis_slice.add_dim_coord(fp_coord, time_dim) - - # Add FRT scalar coordinate. - # - analysis_slice.add_aux_coord( - iris.coords.AuxCoord( - t, - standard_name="forecast_reference_time", - units=frt_coord.units, - ) - ) - - fixed_analysis_cubes.append( - analysis_slice - ) - - if len(fixed_analysis_cubes) == 1: - - cube = fixed_analysis_cubes[0] - - fp_coord = cube.coord("forecast_period") - - time_points = ( - cube.coord( - "forecast_reference_time" - ).points[0] - + fp_coord.points - ) - - time_coord = iris.coords.AuxCoord( - time_points, - standard_name="time", - units=frt_coord.units, - ) - - cube.add_aux_coord( - time_coord, - data_dims=( - cube.coord_dims( - "forecast_period" - )[0], - ), - ) - - else: - # multiple FRTs - # - cube = fixed_analysis_cubes.merge_cube() - - fp_dim = cube.coord_dims( - "forecast_period" - )[0] - - frt_dim = cube.coord_dims( - "forecast_reference_time" - )[0] - - fp_points = cube.coord( - "forecast_period" - ).points - - frt_points = cube.coord( - "forecast_reference_time" - ).points - - time_points = ( - fp_points[:, None] - + frt_points[None, :] - ) - - time_coord = iris.coords.AuxCoord( - time_points, - standard_name="time", - units=frt_coord.units, - ) - - cube.add_aux_coord( - time_coord, - data_dims=(fp_dim, frt_dim), - ) - - fixed_analysis_cubes = iris.cube.CubeList() - fixed_analysis_cubes.append(cube) - - return iris.cube.CubeList([cube, forecast_cube]) - - def collapse( cubes: iris.cube.Cube | iris.cube.CubeList, coordinate: str | list[str], @@ -262,9 +75,6 @@ def collapse( raise ValueError("Must specify additional_percent") # Retain only common time points between different models if multiple model inputs. - - cubes = _fix_analysis_forecasttime(cubes) - if isinstance(cubes, iris.cube.CubeList) and len(cubes) > 1: logging.debug( "Extracting common time points as multiple model inputs detected." From 89abe06f2784de9f7896b0ed29f0445fed649907 Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 29 Jul 2026 09:35:13 +0100 Subject: [PATCH 13/56] intial code move to utils --- utils/proc_reanalysis/README.md | 1 + utils/proc_reanalysis/process_reanalysis.py | 193 ++++++++++++++++++++ 2 files changed, 194 insertions(+) create mode 100644 utils/proc_reanalysis/README.md create mode 100644 utils/proc_reanalysis/process_reanalysis.py diff --git a/utils/proc_reanalysis/README.md b/utils/proc_reanalysis/README.md new file mode 100644 index 0000000000..30404ce4c5 --- /dev/null +++ b/utils/proc_reanalysis/README.md @@ -0,0 +1 @@ +TODO \ No newline at end of file diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py new file mode 100644 index 0000000000..e78a7b744b --- /dev/null +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -0,0 +1,193 @@ +""" +Some code + +MAINTAINER: james.warner@metoffice.gov.uk / jwarner8 +""" + + +import iris + +def _coord_is_effectively_scalar(coord): + """True if coordinate is scalar or all points have the same value.""" + points = np.asarray(coord.points) + + if points.size <= 1: + return True + + return np.all(points == points.flat[0]) + + +def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): + + """ + TODO - in progress + """ + + # Partition out forecast and analysis cubes based on forecast_period. + analysis_cubes = iris.cube.CubeList() + forecast_cubes = iris.cube.CubeList() + + for cube in cubes: + fp = cube.coord("forecast_period") + + if fp.points.max() == 0: + analysis_cubes.append(cube) + else: + forecast_cubes.append(cube) + + if len(analysis_cubes) > 1: + raise ValueError + if len(forecast_cubes) > 1: + print(forecast_cubes) + raise ValueError + if len(analysis_cubes) == 0: + return cubes #i.e. no reanalysis so ignore this function + + + analysis_cube = analysis_cubes[0] + forecast_cube = forecast_cubes[0] + + # Analysis cube will think it has multiple realizations cause it duplicates + analysis_cube = analysis_cube.extract(iris.Constraint(realization = analysis_cube.coord('realization').points[0])) + + # Get forecast_reference_time and period information + frt_coord = forecast_cube.coord("forecast_reference_time") + + fc_period_coord = forecast_cube.coord("forecast_period") + fc_periods = fc_period_coord.points + fc_length = fc_periods[-1] + + # Work out minimum, maximum valid times in analysis + an_min = analysis_cube.coord("time").units.num2date(analysis_cube.coord("time").points[0]) + an_max = analysis_cube.coord("time").units.num2date(analysis_cube.coord("time").points[-1]) + + fixed_analysis_cubes = iris.cube.CubeList() + + # Iterate over all forecast reference times + for t in frt_coord.points: + frt = frt_coord.units.num2date(t) + + if frt < an_min or ( + frt + dt.timedelta(hours=float(fc_length)) + ) > an_max: + + raise ValueError( + f"Analysis does not cover forecast span " + f"{frt} -> " + f"{frt + dt.timedelta(hours=float(fc_length))}") + + + time_constraint = iris.Constraint( + time=lambda cell: frt <= cell.point <= frt+dt.timedelta(hours=fc_length) + ) + + analysis_slice = analysis_cube.extract(time_constraint) + + # Remove old coords + analysis_slice.remove_coord("forecast_reference_time") + analysis_slice.remove_coord("forecast_period") + + + time_coord = analysis_slice.coord("time") + time_dim = analysis_slice.coord_dims("time")[0] + + # Forecast periods relative to this FRT. + # + fp_points = ( + time_coord.points - time_coord.points[0] + ) + + fp_coord = iris.coords.DimCoord( + fp_points, + standard_name="forecast_period", + units='hours', + ) + + # Replace time dimension coord with forecast_period. + # + analysis_slice.remove_coord("time") + analysis_slice.add_dim_coord(fp_coord, time_dim) + + # Add FRT scalar coordinate. + # + analysis_slice.add_aux_coord( + iris.coords.AuxCoord( + t, + standard_name="forecast_reference_time", + units=frt_coord.units, + ) + ) + + fixed_analysis_cubes.append( + analysis_slice + ) + + if len(fixed_analysis_cubes) == 1: + + cube = fixed_analysis_cubes[0] + + fp_coord = cube.coord("forecast_period") + + time_points = ( + cube.coord( + "forecast_reference_time" + ).points[0] + + fp_coord.points + ) + + time_coord = iris.coords.AuxCoord( + time_points, + standard_name="time", + units=frt_coord.units, + ) + + cube.add_aux_coord( + time_coord, + data_dims=( + cube.coord_dims( + "forecast_period" + )[0], + ), + ) + + else: + # multiple FRTs + # + cube = fixed_analysis_cubes.merge_cube() + + fp_dim = cube.coord_dims( + "forecast_period" + )[0] + + frt_dim = cube.coord_dims( + "forecast_reference_time" + )[0] + + fp_points = cube.coord( + "forecast_period" + ).points + + frt_points = cube.coord( + "forecast_reference_time" + ).points + + time_points = ( + fp_points[:, None] + + frt_points[None, :] + ) + + time_coord = iris.coords.AuxCoord( + time_points, + standard_name="time", + units=frt_coord.units, + ) + + cube.add_aux_coord( + time_coord, + data_dims=(fp_dim, frt_dim), + ) + + fixed_analysis_cubes = iris.cube.CubeList() + fixed_analysis_cubes.append(cube) + + return iris.cube.CubeList([cube, forecast_cube]) From f227f7eb6d015244b2ce71ab37ef1c54521a45cc Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 29 Jul 2026 09:39:06 +0100 Subject: [PATCH 14/56] add conf bits back in --- .../meta/observations/rose-meta.conf | 25 +++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf index 8c7bc0d139..bdfedfbe07 100644 --- a/src/CSET/cset_workflow/meta/observations/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -257,3 +257,28 @@ help=Switch to select gridded UK Nimrod 5 minute rainfall rate. compulsory=true type=python_boolean sort-key=nimrod-radar-7 + +################################################################ +# Cardington UK Research Observations +[Observations/Cardington] +ns=Observations/Cardington +sort-key=sec-d4 +title=Cardington UK observations + +# Air temperature timeseries +[template variables=CARDINGTON_AIR_TEMPERATURE_SINGLE_POINT_TIME_SERIES] +ns=Observations/Cardington +title=Cardington air temperature timeseries +description=Create Cardington diagnostic variable timeseries +type=python_boolean +compulsory=true +sort-key=0card1 + +# Relative humidity timeseries +[template variables=CARDINGTON_RELATIVE_HUMIDITY_SINGLE_POINT_TIME_SERIES] +ns=Observations/Cardington +title=Cardington relative humidity timeseries +description=Create Cardington diagnostic variable timeseries +type=python_boolean +compulsory=true +sort-key=0card2 From cca95ddb1eee150a8b34e6eb938bb501a3d39979 Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 29 Jul 2026 11:33:15 +0100 Subject: [PATCH 15/56] working code for UM and ERA5 --- utils/proc_reanalysis/process_reanalysis.py | 328 ++++++++++---------- 1 file changed, 167 insertions(+), 161 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index e78a7b744b..a7511fe3e6 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -4,190 +4,196 @@ MAINTAINER: james.warner@metoffice.gov.uk / jwarner8 """ - import iris +import iris.cube -def _coord_is_effectively_scalar(coord): - """True if coordinate is scalar or all points have the same value.""" - points = np.asarray(coord.points) - - if points.size <= 1: - return True - - return np.all(points == points.flat[0]) - +iris.FUTURE.date_microseconds = True +iris.FUTURE.save_split_attrs = True +import argparse +from datetime import datetime, timedelta -def _fix_analysis_forecasttime(cubes: iris.cube.CubeList): +def _create_forecasts(reanalysis, forecast_initialisations, forecastlength, outpath): """ TODO - in progress """ - - # Partition out forecast and analysis cubes based on forecast_period. - analysis_cubes = iris.cube.CubeList() - forecast_cubes = iris.cube.CubeList() - - for cube in cubes: - fp = cube.coord("forecast_period") - - if fp.points.max() == 0: - analysis_cubes.append(cube) - else: - forecast_cubes.append(cube) - - if len(analysis_cubes) > 1: - raise ValueError - if len(forecast_cubes) > 1: - print(forecast_cubes) - raise ValueError - if len(analysis_cubes) == 0: - return cubes #i.e. no reanalysis so ignore this function - - - analysis_cube = analysis_cubes[0] - forecast_cube = forecast_cubes[0] - - # Analysis cube will think it has multiple realizations cause it duplicates - analysis_cube = analysis_cube.extract(iris.Constraint(realization = analysis_cube.coord('realization').points[0])) - - # Get forecast_reference_time and period information - frt_coord = forecast_cube.coord("forecast_reference_time") - - fc_period_coord = forecast_cube.coord("forecast_period") - fc_periods = fc_period_coord.points - fc_length = fc_periods[-1] - - # Work out minimum, maximum valid times in analysis - an_min = analysis_cube.coord("time").units.num2date(analysis_cube.coord("time").points[0]) - an_max = analysis_cube.coord("time").units.num2date(analysis_cube.coord("time").points[-1]) - - fixed_analysis_cubes = iris.cube.CubeList() - - # Iterate over all forecast reference times - for t in frt_coord.points: - frt = frt_coord.units.num2date(t) - - if frt < an_min or ( - frt + dt.timedelta(hours=float(fc_length)) - ) > an_max: - - raise ValueError( - f"Analysis does not cover forecast span " - f"{frt} -> " - f"{frt + dt.timedelta(hours=float(fc_length))}") - - - time_constraint = iris.Constraint( - time=lambda cell: frt <= cell.point <= frt+dt.timedelta(hours=fc_length) - ) - - analysis_slice = analysis_cube.extract(time_constraint) - - # Remove old coords - analysis_slice.remove_coord("forecast_reference_time") - analysis_slice.remove_coord("forecast_period") - - - time_coord = analysis_slice.coord("time") - time_dim = analysis_slice.coord_dims("time")[0] - - # Forecast periods relative to this FRT. - # - fp_points = ( - time_coord.points - time_coord.points[0] + for forecast in forecast_initialisations: + print(f"Working on forecast initialisation {forecast} out to {forecastlength}H") + + start = forecast + end = forecast + timedelta(hours=forecastlength) + + cutouts = iris.cube.CubeList() + + for cube in reanalysis: + print(f"{cube.name()}...") + + # Work out minimum, maximum valid times in analysis + an_min = cube.coord("time").units.num2date(cube.coord("time").points[0]) + an_max = cube.coord("time").units.num2date(cube.coord("time").points[-1]) + + if start < an_min or end > an_max: + print("WARNING!! SOMEthING...") + else: + time_constraint = iris.Constraint( + time=lambda cell: start <= cell.point <= end + ) + + cube_slice = cube.extract(time_constraint) + + # Remove unnecessary coords and attributes + coords_attrs_to_remove = [ + "forecast_period", + "forecast_reference_time", + "originating_centre", + "source", + "um_version", + ] + for item in coords_attrs_to_remove: + if cube_slice.coords(item): + cube_slice.remove_coord(item) + if item in cube_slice.attributes: + del cube_slice.attributes[item] + + time_coord = cube_slice.coord("time") + time_dim = cube_slice.coord_dims("time")[0] + + # Forecast periods relative to this FRT. + fp_points = time_coord.points - time_coord.points[0] + + units_str = str(time_coord.units) + + if units_str.startswith("seconds since"): + fp_points = fp_points / 3600.0 + fp_units = "hours" + elif units_str.startswith("minutes since"): + fp_points = fp_points / 60.0 + fp_units = "hours" + elif units_str.startswith("hours since"): + fp_units = "hours" + else: + raise ValueError(f"Unhandled time units: {time_coord.units}") + + fp_coord = iris.coords.DimCoord( + fp_points, + standard_name="forecast_period", + units=fp_units, + ) + + cube_slice.remove_coord("time") + cube_slice.add_dim_coord(fp_coord, time_dim) + + cube_slice.add_aux_coord( + iris.coords.AuxCoord( + time_coord.units.date2num(start), + standard_name="forecast_reference_time", + units=time_coord.units, + ) + ) + + time_points = ( + cube_slice.coord("forecast_reference_time").points[0] + + fp_coord.points + ) + + new_time_coord = iris.coords.AuxCoord( + time_points, + standard_name="time", + units=time_coord.units, + ) + + cube_slice.add_aux_coord( + new_time_coord, + data_dims=(cube_slice.coord_dims("forecast_period")[0],), + ) + + print(cube_slice) + quit() + + cutouts.append(cube_slice) + print(f"{cube.name()}...done.") + + print(f"Saving {outpath + '/reanalysis_' + start.strftime('%Y%m%dT%H%MZ')}.nc") + iris.save( + cutouts, outpath + "/reanalysis_" + start.strftime("%Y%m%dT%H%MZ") + ".nc" ) - fp_coord = iris.coords.DimCoord( - fp_points, - standard_name="forecast_period", - units='hours', - ) - # Replace time dimension coord with forecast_period. - # - analysis_slice.remove_coord("time") - analysis_slice.add_dim_coord(fp_coord, time_dim) - - # Add FRT scalar coordinate. - # - analysis_slice.add_aux_coord( - iris.coords.AuxCoord( - t, - standard_name="forecast_reference_time", - units=frt_coord.units, - ) - ) - - fixed_analysis_cubes.append( - analysis_slice - ) +def _identify_number_of_cycles_required(cyclestart, cycleend, cyclefreq): + """ + TODO + """ + start_dt = datetime.strptime(cyclestart, "%Y%m%dT%H%MZ") + end_dt = datetime.strptime(cycleend, "%Y%m%dT%H%MZ") - if len(fixed_analysis_cubes) == 1: + forecast_initialisations = [] + current = start_dt - cube = fixed_analysis_cubes[0] + while current <= end_dt: + forecast_initialisations.append(current) + current += timedelta(hours=cyclefreq) - fp_coord = cube.coord("forecast_period") + return forecast_initialisations - time_points = ( - cube.coord( - "forecast_reference_time" - ).points[0] - + fp_coord.points - ) - time_coord = iris.coords.AuxCoord( - time_points, - standard_name="time", - units=frt_coord.units, - ) +def main(): - cube.add_aux_coord( - time_coord, - data_dims=( - cube.coord_dims( - "forecast_period" - )[0], - ), - ) + parser = argparse.ArgumentParser(description="Process forecast data.") - else: - # multiple FRTs - # - cube = fixed_analysis_cubes.merge_cube() + parser.add_argument("--filepath", required=True, help="Path to file(s)") + parser.add_argument( + "--cyclestart", type=str, required=True, help="First forecast initiation/cycle" + ) + parser.add_argument( + "--cycleend", type=str, required=True, help="Final forecast initiation/cycle" + ) + parser.add_argument( + "--cyclefreq", + type=int, + required=True, + help="Time between forecast initiationss/cycles", + ) + parser.add_argument( + "--forecastlength", + type=int, + required=True, + help="Forecast length in SI units i.e. PT48H", + ) + parser.add_argument( + "--outpath", type=str, required=True, help="Where to write output data" + ) - fp_dim = cube.coord_dims( - "forecast_period" - )[0] + args = parser.parse_args() - frt_dim = cube.coord_dims( - "forecast_reference_time" - )[0] + # Populate required variables + filepath = args.filepath + cyclestart = args.cyclestart + cycleend = args.cycleend + cyclefreq = args.cyclefreq + forecastlength = args.forecastlength + outpath = args.outpath - fp_points = cube.coord( - "forecast_period" - ).points + print("") + print("Starting process_reanalysis.py...") - frt_points = cube.coord( - "forecast_reference_time" - ).points + # Get all forecast initiations + forecast_initialisations = _identify_number_of_cycles_required( + cyclestart, cycleend, cyclefreq + ) - time_points = ( - fp_points[:, None] - + frt_points[None, :] - ) + # Load all reanalysis supplied + print(f"Loading reanalysis from {filepath}") + reanalysis = iris.load(filepath) + print("") + print("Found the following cubes...") + print(reanalysis) - time_coord = iris.coords.AuxCoord( - time_points, - standard_name="time", - units=frt_coord.units, - ) + print("") + print("Creating postprocessed files...") + _create_forecasts(reanalysis, forecast_initialisations, forecastlength, outpath) - cube.add_aux_coord( - time_coord, - data_dims=(fp_dim, frt_dim), - ) + print("Done") - fixed_analysis_cubes = iris.cube.CubeList() - fixed_analysis_cubes.append(cube) - return iris.cube.CubeList([cube, forecast_cube]) +if __name__ == "__main__": + main() From 8b3e66c777acfe21b35a4a5c9681c465428d666f Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 29 Jul 2026 11:47:42 +0100 Subject: [PATCH 16/56] add docstring to main --- utils/proc_reanalysis/process_reanalysis.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index a7511fe3e6..433d95fb3f 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -1,7 +1,10 @@ """ -Some code +Code that restructures reanalysis data to give it an effective forecast_period, so it can +be directly compared to model forecasts in CSET. +The code base currently supports UM and ERA5, and does not perform any additional metadata +correction beyond time axis and removing some surplus coords/attributes. -MAINTAINER: james.warner@metoffice.gov.uk / jwarner8 +Please see README for futher information on how to run the script. """ import iris @@ -136,7 +139,13 @@ def _identify_number_of_cycles_required(cyclestart, cycleend, cyclefreq): return forecast_initialisations -def main(): +def main() -> None: + """Generate forecast-like datasets from reanalysis data. + + Parse command-line arguments, create forecast initialisation times, + process the input reanalysis data, and write the resulting forecast + files to disk. + """ parser = argparse.ArgumentParser(description="Process forecast data.") From b9afb02c44d36b4c2e4c4a0b0f1af3d49cb7e56e Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 29 Jul 2026 12:02:17 +0100 Subject: [PATCH 17/56] tidy up function docstring, comments --- utils/proc_reanalysis/process_reanalysis.py | 32 ++++++++++++++++----- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index 433d95fb3f..4afb0302a6 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -122,16 +122,35 @@ def _create_forecasts(reanalysis, forecast_initialisations, forecastlength, outp ) -def _identify_number_of_cycles_required(cyclestart, cycleend, cyclefreq): - """ - TODO +def _identify_number_of_cycles_required( + cyclestart: str, cycleend: str, cyclefreq: str +) -> list: + """Generate forecast initialisation datetimes between two cycle bounds. + + Parameters + ---------- + cyclestart : str + First forecast cycle in YYYYMMDDTHHMMZ format. + cycleend : str + Last forecast cycle in YYYYMMDDTHHMMZ format. + cyclefreq : int + Frequency between forecast cycles in hours. + + Returns + ------- + forecast_initialisations: list + Forecast initialisation datetimes from cyclestart to cycleend, + inclusive, separated by cyclefreq hours. """ + # Identify the start and end times, and create datetime objects for these start_dt = datetime.strptime(cyclestart, "%Y%m%dT%H%MZ") end_dt = datetime.strptime(cycleend, "%Y%m%dT%H%MZ") + # To store initialisation times forecast_initialisations = [] current = start_dt + # Iterate over all initiations within the bounds, using the cyclefreq to determine interval. while current <= end_dt: forecast_initialisations.append(current) current += timedelta(hours=cyclefreq) @@ -146,7 +165,6 @@ def main() -> None: process the input reanalysis data, and write the resulting forecast files to disk. """ - parser = argparse.ArgumentParser(description="Process forecast data.") parser.add_argument("--filepath", required=True, help="Path to file(s)") @@ -182,7 +200,7 @@ def main() -> None: forecastlength = args.forecastlength outpath = args.outpath - print("") + print() print("Starting process_reanalysis.py...") # Get all forecast initiations @@ -193,11 +211,11 @@ def main() -> None: # Load all reanalysis supplied print(f"Loading reanalysis from {filepath}") reanalysis = iris.load(filepath) - print("") + print() print("Found the following cubes...") print(reanalysis) - print("") + print() print("Creating postprocessed files...") _create_forecasts(reanalysis, forecast_initialisations, forecastlength, outpath) From 9a978799185e4d0c39a92aed699246b7f2575f44 Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 29 Jul 2026 12:17:36 +0100 Subject: [PATCH 18/56] comments for _create_forecasts and docstring --- utils/proc_reanalysis/process_reanalysis.py | 70 ++++++++++++++++----- 1 file changed, 56 insertions(+), 14 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index 4afb0302a6..6ba373a02e 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -16,18 +16,54 @@ from datetime import datetime, timedelta -def _create_forecasts(reanalysis, forecast_initialisations, forecastlength, outpath): - """ - TODO - in progress +from datetime import datetime +from pathlib import Path + +import iris + + +def _create_forecasts( + reanalysis: iris.cube.CubeList, + forecast_initialisations: list, + forecastlength: int, + outpath: str, +) -> None: + """Create forecast files from reanalysis data. + + For each forecast initialisation time, extract the corresponding + analysis period from each input cube and convert it into a + forecast-style representation. This includes generating + forecast-period and forecast-reference-time coordinates and + writing the resulting cubes to disk. + + Parameters + ---------- + reanalysis: iris.cube.CubeList + Collection of reanalysis cubes from which forecast periods + will be extracted. + forecast_initialisations: list + Forecast initialisation times to process. + forecastlength: int + Forecast length in hours. + outpath: str + Directory to which the generated forecast files will be saved. + + Returns + ------- + None """ + + # Iterate over all forecast initialisations sequentially. for forecast in forecast_initialisations: print(f"Working on forecast initialisation {forecast} out to {forecastlength}H") + # Work out start and end time start = forecast end = forecast + timedelta(hours=forecastlength) cutouts = iris.cube.CubeList() + # For each cube (variable) loaded for cube in reanalysis: print(f"{cube.name()}...") @@ -35,13 +71,17 @@ def _create_forecasts(reanalysis, forecast_initialisations, forecastlength, outp an_min = cube.coord("time").units.num2date(cube.coord("time").points[0]) an_max = cube.coord("time").units.num2date(cube.coord("time").points[-1]) + # Check reanalysis spans what we are looking for time wise, otherwise ignore if start < an_min or end > an_max: - print("WARNING!! SOMEthING...") + print(f"Warning: Required time {start} {end} outside that found in analysis {an_min} {an_max}") else: + + # Generate time constraint object inclusive of time bounds. time_constraint = iris.Constraint( time=lambda cell: start <= cell.point <= end ) + # Extract required timeslice. cube_slice = cube.extract(time_constraint) # Remove unnecessary coords and attributes @@ -58,12 +98,15 @@ def _create_forecasts(reanalysis, forecast_initialisations, forecastlength, outp if item in cube_slice.attributes: del cube_slice.attributes[item] + # Get a copy of time coord, and dimension this corresponds to. time_coord = cube_slice.coord("time") + time_coord_points = time_coord.points.copy() time_dim = cube_slice.coord_dims("time")[0] - # Forecast periods relative to this FRT. + # Forecast periods relative to initialisation. fp_points = time_coord.points - time_coord.points[0] + # Work out units of forecast_period and adjust if necessary. units_str = str(time_coord.units) if units_str.startswith("seconds since"): @@ -77,15 +120,18 @@ def _create_forecasts(reanalysis, forecast_initialisations, forecastlength, outp else: raise ValueError(f"Unhandled time units: {time_coord.units}") + # Create forecast period dimension fp_coord = iris.coords.DimCoord( fp_points, standard_name="forecast_period", units=fp_units, ) + # Remove time dimension temporarily, as forecast_period will be lead dimension cube_slice.remove_coord("time") cube_slice.add_dim_coord(fp_coord, time_dim) + # Add auxiliary forecast initialisation dimension cube_slice.add_aux_coord( iris.coords.AuxCoord( time_coord.units.date2num(start), @@ -94,28 +140,24 @@ def _create_forecasts(reanalysis, forecast_initialisations, forecastlength, outp ) ) - time_points = ( - cube_slice.coord("forecast_reference_time").points[0] - + fp_coord.points - ) - + # Create auxiliary time dimension of valid time. new_time_coord = iris.coords.AuxCoord( - time_points, + time_coord_points, standard_name="time", units=time_coord.units, ) + # Add this dimsnion to the cube, tied to the forecast_period dimension. cube_slice.add_aux_coord( new_time_coord, data_dims=(cube_slice.coord_dims("forecast_period")[0],), ) - print(cube_slice) - quit() - + # Append slice to cutout list, ready for saving. cutouts.append(cube_slice) print(f"{cube.name()}...done.") + # Once all cubes processed, save to disk. print(f"Saving {outpath + '/reanalysis_' + start.strftime('%Y%m%dT%H%MZ')}.nc") iris.save( cutouts, outpath + "/reanalysis_" + start.strftime("%Y%m%dT%H%MZ") + ".nc" From b0d82d675e1b8061159f661593bf1e175a9ca3f2 Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 29 Jul 2026 12:23:14 +0100 Subject: [PATCH 19/56] do pre-commit fixes --- utils/proc_reanalysis/process_reanalysis.py | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index 6ba373a02e..62cd7aa33d 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -1,10 +1,11 @@ """ -Code that restructures reanalysis data to give it an effective forecast_period, so it can -be directly compared to model forecasts in CSET. +Code that restructures reanalysis data to give it an effective forecast_period. + +As a result, reanalysis can be directly compared to model forecasts in CSET and treated as another model. The code base currently supports UM and ERA5, and does not perform any additional metadata correction beyond time axis and removing some surplus coords/attributes. -Please see README for futher information on how to run the script. +Please see README for further information on how to run the script. """ import iris @@ -15,10 +16,6 @@ import argparse from datetime import datetime, timedelta - -from datetime import datetime -from pathlib import Path - import iris @@ -52,7 +49,6 @@ def _create_forecasts( ------- None """ - # Iterate over all forecast initialisations sequentially. for forecast in forecast_initialisations: print(f"Working on forecast initialisation {forecast} out to {forecastlength}H") @@ -73,12 +69,13 @@ def _create_forecasts( # Check reanalysis spans what we are looking for time wise, otherwise ignore if start < an_min or end > an_max: - print(f"Warning: Required time {start} {end} outside that found in analysis {an_min} {an_max}") + print( + f"Warning: Required time {start} {end} outside that found in analysis {an_min} {an_max}" + ) else: - # Generate time constraint object inclusive of time bounds. time_constraint = iris.Constraint( - time=lambda cell: start <= cell.point <= end + time=lambda cell, start=start, end=end: start <= cell.point <= end ) # Extract required timeslice. From b7b0752c8d57f69745ece8bba6cfb98b83933fec Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 29 Jul 2026 16:39:41 +0100 Subject: [PATCH 20/56] add tests placeholder --- .../test_process_reanalysis.py | 262 ++++++++++++++++++ 1 file changed, 262 insertions(+) create mode 100644 utils/proc_reanalysis/test_process_reanalysis.py diff --git a/utils/proc_reanalysis/test_process_reanalysis.py b/utils/proc_reanalysis/test_process_reanalysis.py new file mode 100644 index 0000000000..477f7786a6 --- /dev/null +++ b/utils/proc_reanalysis/test_process_reanalysis.py @@ -0,0 +1,262 @@ +""" +Unit tests - some formed from GitHub copilot for attribution purposes +""" +# conftest.py or in test module + +import numpy as np +import iris +from iris.coords import DimCoord +from iris.cube import Cube +import pytest +from datetime import datetime + + +def test_single_cycle(): + """Assert single datetime returned if one initialisation.""" + result = _identify_number_of_cycles_required( + "20240101T0000Z", + "20240101T0000Z", + 6, + ) + + assert result == [datetime(2024, 1, 1, 0, 0)] + + +def test_multiple_cycles(): + """Test handling of multiple cycles identified.""" + result = _identify_number_of_cycles_required( + "20240101T0000Z", + "20240101T1200Z", + 6, + ) + + assert result == [ + datetime(2024, 1, 1, 0), + datetime(2024, 1, 1, 6), + datetime(2024, 1, 1, 12), + ] + + +def test_non_divisible_interval(): + """Check end point is not exceeded.""" + result = _identify_number_of_cycles_required( + "20240101T0000Z", + "20240101T1000Z", + 6, + ) + + assert result == [ + datetime(2024, 1, 1, 0), + datetime(2024, 1, 1, 6), + ] + + +def test_invalid_date_format(): + """Raise problem with input.""" + with pytest.raises(ValueError): + _identify_number_of_cycles_required( + "2024-01-01", + "20240101T1200Z", + 6, + ) + + +def make_cube( + name="air_temperature", + units="hours since 2024-01-01 00:00:00", +): + """Create a minimal synthetic reanalysis cube.""" + + time = DimCoord( + np.arange(5), + standard_name="time", + units=units, + ) + + lat = DimCoord([50.0], standard_name="latitude", units="degrees") + lon = DimCoord([0.0], standard_name="longitude", units="degrees") + + data = np.arange(5).reshape(5, 1, 1) + + cube = Cube( + data, + standard_name=name, + dim_coords_and_dims=[ + (time, 0), + (lat, 1), + (lon, 2), + ], + ) + + return cube + + +def test_forecast_period_created(tmp_path): + """Check forecast period constructed correctly.""" + cube = make_cube() + + _create_forecasts( + iris.cube.CubeList([cube]), + [datetime(2024, 1, 1, 0)], + forecastlength=4, + outpath=str(tmp_path), + ) + + outfile = tmp_path / "reanalysis_20240101T0000Z.nc" + + cubes = iris.load(str(outfile)) + + result = cubes[0] + + fp = result.coord("forecast_period") + + assert fp.points.tolist() == [0, 1, 2, 3, 4] + assert str(fp.units) == "hours" + + +def test_seconds_converted_to_hours(tmp_path): + """Check conversion from seconds to hours in forecast period.""" + cube = make_cube( + time_units="seconds since 2024-01-01 00:00:00" + ) + + cube.coord("time").points = [0, 3600, 7200] + + _create_forecasts( + iris.cube.CubeList([cube]), + [datetime(2024, 1, 1)], + forecastlength=2, + outpath=str(tmp_path), + ) + + cubes = iris.load( + str(tmp_path / "reanalysis_20240101T0000Z.nc") + ) + + fp = cubes[0].coord("forecast_period") + + assert fp.points.tolist() == [0, 1, 2] + + +def test_minutes_converted_to_hours(tmp_path): + """Check conversion from minutes to hours in forecast period.""" + cube = make_cube( + time_units="minutes since 2024-01-01 00:00:00" + ) + + cube.coord("time").points = [0, 60, 120] + + _create_forecasts( + iris.cube.CubeList([cube]), + [datetime(2024, 1, 1)], + forecastlength=2, + outpath=str(tmp_path), + ) + + cubes = iris.load( + str(tmp_path / "reanalysis_20240101T0000Z.nc") + ) + + fp = cubes[0].coord("forecast_period") + + assert fp.points.tolist() == [0, 1, 2] + + +def test_unknown_time_units_raise(tmp_path): + """Check that error raised if time units unhandled.""" + cube = make_cube( + time_units="foobars since 2024-01-01 00:00:00" + ) + + with pytest.raises(ValueError, match="Unhandled time units"): + _create_forecasts( + iris.cube.CubeList([cube]), + [datetime(2024, 1, 1)], + forecastlength=1, + outpath=str(tmp_path), + ) + + +def test_forecast_reference_time_created(tmp_path): + """Check that forecast reference time has been set correctly.""" + init_time = datetime(2024, 1, 1, 0) + + cube = make_cube() + + _create_forecasts( + iris.cube.CubeList([cube]), + [init_time], + forecastlength=4, + outpath=str(tmp_path), + ) + + cubes = iris.load( + str(tmp_path / "reanalysis_20240101T0000Z.nc") + ) + + frt = cubes[0].coord("forecast_reference_time") + + recovered = frt.units.num2date(frt.points[0]) + + assert recovered.replace(tzinfo=None) == init_time + + +def test_forecast_attributes_removed(tmp_path): + """Check that common analysis attributes have been removed.""" + cube = make_cube() + + cube.attributes["source"] = "ERA5" + cube.attributes["um_version"] = "13.0" + + _create_forecasts( + iris.cube.CubeList([cube]), + [datetime(2024, 1, 1)], + forecastlength=4, + outpath=str(tmp_path), + ) + + cubes = iris.load( + str(tmp_path / "reanalysis_20240101T0000Z.nc") + ) + + attrs = cubes[0].attributes + + assert "source" not in attrs + assert "um_version" not in attrs + + +def test_cube_skipped_if_insufficient_data(tmp_path): + """Check that nothing returned if analysis doesn't overlap with target forecast.""" + cube = make_cube() + + _create_forecasts( + iris.cube.CubeList([cube]), + [datetime(2024, 1, 1)], + forecastlength=10, + outpath=str(tmp_path), + ) + + cubes = iris.load( + str(tmp_path / "reanalysis_20240101T0000Z.nc") + ) + + assert len(cubes) == 0 + + +def test_multiple_cubes_processed(tmp_path): + """Check working with multiple cubes.""" + cube1 = make_cube("air_temperature") + cube2 = make_cube("air_pressure") + + _create_forecasts( + iris.cube.CubeList([cube1, cube2]), + [datetime(2024, 1, 1)], + forecastlength=4, + outpath=str(tmp_path), + ) + + cubes = iris.load( + str(tmp_path / "reanalysis_20240101T0000Z.nc") + ) + + assert len(cubes) == 2 From f6df532599d9e66014f562cb5e37189f50369111 Mon Sep 17 00:00:00 2001 From: James Warner Date: Mon, 3 Aug 2026 10:53:23 +0100 Subject: [PATCH 21/56] add imports --- .../test_process_reanalysis.py | 73 +++++++------------ 1 file changed, 27 insertions(+), 46 deletions(-) diff --git a/utils/proc_reanalysis/test_process_reanalysis.py b/utils/proc_reanalysis/test_process_reanalysis.py index 477f7786a6..2522c69f66 100644 --- a/utils/proc_reanalysis/test_process_reanalysis.py +++ b/utils/proc_reanalysis/test_process_reanalysis.py @@ -1,19 +1,19 @@ -""" -Unit tests - some formed from GitHub copilot for attribution purposes -""" +"""Unit tests.""" # conftest.py or in test module -import numpy as np +from datetime import datetime + import iris +import numpy as np +import process_reanalysis as proc_reanalysis +import pytest from iris.coords import DimCoord from iris.cube import Cube -import pytest -from datetime import datetime def test_single_cycle(): """Assert single datetime returned if one initialisation.""" - result = _identify_number_of_cycles_required( + result = proc_reanalysis._identify_number_of_cycles_required( "20240101T0000Z", "20240101T0000Z", 6, @@ -24,7 +24,7 @@ def test_single_cycle(): def test_multiple_cycles(): """Test handling of multiple cycles identified.""" - result = _identify_number_of_cycles_required( + result = proc_reanalysis._identify_number_of_cycles_required( "20240101T0000Z", "20240101T1200Z", 6, @@ -39,7 +39,7 @@ def test_multiple_cycles(): def test_non_divisible_interval(): """Check end point is not exceeded.""" - result = _identify_number_of_cycles_required( + result = proc_reanalysis._identify_number_of_cycles_required( "20240101T0000Z", "20240101T1000Z", 6, @@ -54,7 +54,7 @@ def test_non_divisible_interval(): def test_invalid_date_format(): """Raise problem with input.""" with pytest.raises(ValueError): - _identify_number_of_cycles_required( + proc_reanalysis._identify_number_of_cycles_required( "2024-01-01", "20240101T1200Z", 6, @@ -66,7 +66,6 @@ def make_cube( units="hours since 2024-01-01 00:00:00", ): """Create a minimal synthetic reanalysis cube.""" - time = DimCoord( np.arange(5), standard_name="time", @@ -95,7 +94,7 @@ def test_forecast_period_created(tmp_path): """Check forecast period constructed correctly.""" cube = make_cube() - _create_forecasts( + proc_reanalysis._create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1, 0)], forecastlength=4, @@ -116,22 +115,18 @@ def test_forecast_period_created(tmp_path): def test_seconds_converted_to_hours(tmp_path): """Check conversion from seconds to hours in forecast period.""" - cube = make_cube( - time_units="seconds since 2024-01-01 00:00:00" - ) + cube = make_cube(time_units="seconds since 2024-01-01 00:00:00") cube.coord("time").points = [0, 3600, 7200] - _create_forecasts( + proc_reanalysis._create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], forecastlength=2, outpath=str(tmp_path), ) - cubes = iris.load( - str(tmp_path / "reanalysis_20240101T0000Z.nc") - ) + cubes = iris.load(str(tmp_path / "reanalysis_20240101T0000Z.nc")) fp = cubes[0].coord("forecast_period") @@ -140,22 +135,18 @@ def test_seconds_converted_to_hours(tmp_path): def test_minutes_converted_to_hours(tmp_path): """Check conversion from minutes to hours in forecast period.""" - cube = make_cube( - time_units="minutes since 2024-01-01 00:00:00" - ) + cube = make_cube(time_units="minutes since 2024-01-01 00:00:00") cube.coord("time").points = [0, 60, 120] - _create_forecasts( + proc_reanalysis._create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], forecastlength=2, outpath=str(tmp_path), ) - cubes = iris.load( - str(tmp_path / "reanalysis_20240101T0000Z.nc") - ) + cubes = iris.load(str(tmp_path / "reanalysis_20240101T0000Z.nc")) fp = cubes[0].coord("forecast_period") @@ -164,12 +155,10 @@ def test_minutes_converted_to_hours(tmp_path): def test_unknown_time_units_raise(tmp_path): """Check that error raised if time units unhandled.""" - cube = make_cube( - time_units="foobars since 2024-01-01 00:00:00" - ) + cube = make_cube(time_units="foobars since 2024-01-01 00:00:00") with pytest.raises(ValueError, match="Unhandled time units"): - _create_forecasts( + proc_reanalysis._create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], forecastlength=1, @@ -183,16 +172,14 @@ def test_forecast_reference_time_created(tmp_path): cube = make_cube() - _create_forecasts( + proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [init_time], forecastlength=4, outpath=str(tmp_path), ) - cubes = iris.load( - str(tmp_path / "reanalysis_20240101T0000Z.nc") - ) + cubes = iris.load(str(tmp_path / "reanalysis_20240101T0000Z.nc")) frt = cubes[0].coord("forecast_reference_time") @@ -208,16 +195,14 @@ def test_forecast_attributes_removed(tmp_path): cube.attributes["source"] = "ERA5" cube.attributes["um_version"] = "13.0" - _create_forecasts( + proc_reanalysis._create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], forecastlength=4, outpath=str(tmp_path), ) - cubes = iris.load( - str(tmp_path / "reanalysis_20240101T0000Z.nc") - ) + cubes = iris.load(str(tmp_path / "reanalysis_20240101T0000Z.nc")) attrs = cubes[0].attributes @@ -229,16 +214,14 @@ def test_cube_skipped_if_insufficient_data(tmp_path): """Check that nothing returned if analysis doesn't overlap with target forecast.""" cube = make_cube() - _create_forecasts( + proc_reanalysis._create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], forecastlength=10, outpath=str(tmp_path), ) - cubes = iris.load( - str(tmp_path / "reanalysis_20240101T0000Z.nc") - ) + cubes = iris.load(str(tmp_path / "reanalysis_20240101T0000Z.nc")) assert len(cubes) == 0 @@ -248,15 +231,13 @@ def test_multiple_cubes_processed(tmp_path): cube1 = make_cube("air_temperature") cube2 = make_cube("air_pressure") - _create_forecasts( + proc_reanalysis._create_forecasts( iris.cube.CubeList([cube1, cube2]), [datetime(2024, 1, 1)], forecastlength=4, outpath=str(tmp_path), ) - cubes = iris.load( - str(tmp_path / "reanalysis_20240101T0000Z.nc") - ) + cubes = iris.load(str(tmp_path / "reanalysis_20240101T0000Z.nc")) assert len(cubes) == 2 From c4863592c0d4594aabaa3af438ed24e505d4e013 Mon Sep 17 00:00:00 2001 From: James Warner Date: Mon, 3 Aug 2026 13:31:14 +0100 Subject: [PATCH 22/56] remove double import --- utils/proc_reanalysis/process_reanalysis.py | 28 +++++++++++++-------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index 62cd7aa33d..93c19cd3b2 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -16,8 +16,6 @@ import argparse from datetime import datetime, timedelta -import iris - def _create_forecasts( reanalysis: iris.cube.CubeList, @@ -97,15 +95,13 @@ def _create_forecasts( # Get a copy of time coord, and dimension this corresponds to. time_coord = cube_slice.coord("time") - time_coord_points = time_coord.points.copy() - time_dim = cube_slice.coord_dims("time")[0] - - # Forecast periods relative to initialisation. - fp_points = time_coord.points - time_coord.points[0] # Work out units of forecast_period and adjust if necessary. units_str = str(time_coord.units) + # Forecast periods relative to initialisation. + fp_points = time_coord.points - time_coord.points[0] + if units_str.startswith("seconds since"): fp_points = fp_points / 3600.0 fp_units = "hours" @@ -117,6 +113,10 @@ def _create_forecasts( else: raise ValueError(f"Unhandled time units: {time_coord.units}") + # Get copy of time coordinate + time_coord_points = time_coord.points.copy() + time_dim = cube_slice.coord_dims("time")[0] + # Create forecast period dimension fp_coord = iris.coords.DimCoord( fp_points, @@ -155,10 +155,16 @@ def _create_forecasts( print(f"{cube.name()}...done.") # Once all cubes processed, save to disk. - print(f"Saving {outpath + '/reanalysis_' + start.strftime('%Y%m%dT%H%MZ')}.nc") - iris.save( - cutouts, outpath + "/reanalysis_" + start.strftime("%Y%m%dT%H%MZ") + ".nc" - ) + if len(cutouts) > 0: + print( + f"Saving {outpath + '/reanalysis_' + start.strftime('%Y%m%dT%H%MZ')}.nc" + ) + iris.save( + cutouts, + outpath + "/reanalysis_" + start.strftime("%Y%m%dT%H%MZ") + ".nc", + ) + else: + raise ValueError("No suitable cubes found for saving!") def _identify_number_of_cycles_required( From c91ced8d5baf500ede42fae6deb8e092eb292248 Mon Sep 17 00:00:00 2001 From: James Warner Date: Mon, 3 Aug 2026 13:34:11 +0100 Subject: [PATCH 23/56] add tests --- .../test_process_reanalysis.py | 36 +++++++++---------- 1 file changed, 16 insertions(+), 20 deletions(-) diff --git a/utils/proc_reanalysis/test_process_reanalysis.py b/utils/proc_reanalysis/test_process_reanalysis.py index 2522c69f66..969748aae4 100644 --- a/utils/proc_reanalysis/test_process_reanalysis.py +++ b/utils/proc_reanalysis/test_process_reanalysis.py @@ -1,5 +1,4 @@ -"""Unit tests.""" -# conftest.py or in test module +"""Unit tests for process_reanalysis.py.""" from datetime import datetime @@ -115,9 +114,9 @@ def test_forecast_period_created(tmp_path): def test_seconds_converted_to_hours(tmp_path): """Check conversion from seconds to hours in forecast period.""" - cube = make_cube(time_units="seconds since 2024-01-01 00:00:00") + cube = make_cube(units="seconds since 2024-01-01 00:00:00") - cube.coord("time").points = [0, 3600, 7200] + cube.coord("time").points = [0, 3600, 7200, 10800, 14400] proc_reanalysis._create_forecasts( iris.cube.CubeList([cube]), @@ -135,9 +134,9 @@ def test_seconds_converted_to_hours(tmp_path): def test_minutes_converted_to_hours(tmp_path): """Check conversion from minutes to hours in forecast period.""" - cube = make_cube(time_units="minutes since 2024-01-01 00:00:00") + cube = make_cube(units="minutes since 2024-01-01 00:00:00") - cube.coord("time").points = [0, 60, 120] + cube.coord("time").points = [0, 60, 120, 180, 240] proc_reanalysis._create_forecasts( iris.cube.CubeList([cube]), @@ -155,7 +154,7 @@ def test_minutes_converted_to_hours(tmp_path): def test_unknown_time_units_raise(tmp_path): """Check that error raised if time units unhandled.""" - cube = make_cube(time_units="foobars since 2024-01-01 00:00:00") + cube = make_cube(units="days since 2024-01-01 00:00:00") with pytest.raises(ValueError, match="Unhandled time units"): proc_reanalysis._create_forecasts( @@ -172,7 +171,7 @@ def test_forecast_reference_time_created(tmp_path): cube = make_cube() - proc_reanalysis.create_forecasts( + proc_reanalysis._create_forecasts( iris.cube.CubeList([cube]), [init_time], forecastlength=4, @@ -185,7 +184,7 @@ def test_forecast_reference_time_created(tmp_path): recovered = frt.units.num2date(frt.points[0]) - assert recovered.replace(tzinfo=None) == init_time + assert str(recovered) == "2024-01-01 00:00:00" def test_forecast_attributes_removed(tmp_path): @@ -211,19 +210,16 @@ def test_forecast_attributes_removed(tmp_path): def test_cube_skipped_if_insufficient_data(tmp_path): - """Check that nothing returned if analysis doesn't overlap with target forecast.""" + """Check that nothing save if analysis doesn't overlap with target forecast.""" cube = make_cube() - proc_reanalysis._create_forecasts( - iris.cube.CubeList([cube]), - [datetime(2024, 1, 1)], - forecastlength=10, - outpath=str(tmp_path), - ) - - cubes = iris.load(str(tmp_path / "reanalysis_20240101T0000Z.nc")) - - assert len(cubes) == 0 + with pytest.raises(ValueError, match="No suitable cubes found for saving!"): + proc_reanalysis._create_forecasts( + iris.cube.CubeList([cube]), + [datetime(2024, 1, 1)], + forecastlength=10, + outpath=str(tmp_path), + ) def test_multiple_cubes_processed(tmp_path): From 710097c5ab1b997a4daf81b27f7fb5834b29f099 Mon Sep 17 00:00:00 2001 From: James Warner Date: Mon, 3 Aug 2026 13:48:53 +0100 Subject: [PATCH 24/56] lodge about section to readme --- utils/proc_reanalysis/README.md | 59 ++++++++++++++++++++++++++++++++- 1 file changed, 58 insertions(+), 1 deletion(-) diff --git a/utils/proc_reanalysis/README.md b/utils/proc_reanalysis/README.md index 30404ce4c5..a8848c4cf7 100644 --- a/utils/proc_reanalysis/README.md +++ b/utils/proc_reanalysis/README.md @@ -1 +1,58 @@ -TODO \ No newline at end of file +# process_reanalysis + +## About + +The script `process_reanalysis.py` is a utility for converting atmospheric reanalysis datasets into a forecast-style format that can be directly compared with numerical weather prediction (NWP) model forecasts. + +The primary motivation for this tool is to be able to analyse reanalysis alongside model forecasts in CSET, by choosing reanalysis to be the base model for verification and evaluation. CSET expects forecast data to contain forecast metadata such as: + +- `forecast_reference_time` +- `forecast_period` +- `time` + +Reanalysis datasets typically only contain valid time `time`, and are a series of files where the `forecast_reference_time` changes every 6 hours, and the `forecast_period` is always zero. + +This script resolves this issue by transforming reanalysis data into an effective forecast representation. Rather than treating reanalysis as a special data source in CSET, the transformed output can be ingested directly into the CSET workflow. This allows reanalysis to be treated as another "model" within verification systems such as CSET. + +The script currently supports datasets that can be loaded by Iris and has primarily been developed and tested using: + +- ERA5 reanalysis +- Unified Model (UM) analysis data + +Other model analyses may work, as only the time dimensions are manipulated, but this has not been tested. +No scientific changes are made to the meteorological fields themselves. + +For each requested forecast cycle, for each variable the script: +1. Extracts the required period of time from reanalysis data. +2. Treats the start of that extraction as a forecast initialisation. +3. Generates a forecast period coordinate. +4. Generates a forecast reference time coordinate. +5. Preserves the original valid time information. +6. Saves the result as a forecast-style NetCDF file. + +## Usage + +The python script requires the Iris package to be present within the python install. + +Run it with: + +``` +python process_reanalysis.py \ + --filepath "" \ + --cyclestart YYYYMMDDTHHMMZ \ + --cycleend YYYYMMDDTHHMMZ \ + --cyclefreq \ + --forecastlength \ + --outpath "" +``` + +[ncgen-docs]: https://docs.unidata.ucar.edu/nug/current/netcdf_utilities_guide.html + +## Owners + +> [!TIP] +> Utilities must have at least one named owner. + +The following people should be contacted for queries or issues with this utility: + +* [@jfrost-mo](https://github.com/jfrost-mo) From 13485332ebc8bddaa7a771ece04b1d4e2d466aaf Mon Sep 17 00:00:00 2001 From: James Warner Date: Mon, 3 Aug 2026 14:12:10 +0100 Subject: [PATCH 25/56] add req args --- utils/proc_reanalysis/README.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/utils/proc_reanalysis/README.md b/utils/proc_reanalysis/README.md index a8848c4cf7..0160163f9a 100644 --- a/utils/proc_reanalysis/README.md +++ b/utils/proc_reanalysis/README.md @@ -46,6 +46,10 @@ python process_reanalysis.py \ --outpath "" ``` +Required Arguments: + +--filepath: Path to the input reanalysis data. This can be a single file or wildcard expression understood by Iris. + [ncgen-docs]: https://docs.unidata.ucar.edu/nug/current/netcdf_utilities_guide.html ## Owners From 5d48bac5256a80ea03cb68e1a9889334d881dd31 Mon Sep 17 00:00:00 2001 From: James Warner Date: Mon, 3 Aug 2026 14:12:37 +0100 Subject: [PATCH 26/56] add req args --- utils/proc_reanalysis/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/proc_reanalysis/README.md b/utils/proc_reanalysis/README.md index 0160163f9a..47d61e2545 100644 --- a/utils/proc_reanalysis/README.md +++ b/utils/proc_reanalysis/README.md @@ -48,7 +48,7 @@ python process_reanalysis.py \ Required Arguments: ---filepath: Path to the input reanalysis data. This can be a single file or wildcard expression understood by Iris. +`--filepath`: Path to the input reanalysis data. This can be a single file or wildcard expression understood by Iris. [ncgen-docs]: https://docs.unidata.ucar.edu/nug/current/netcdf_utilities_guide.html From de3553f6b34db306bc85b5e8c2f080d665558e61 Mon Sep 17 00:00:00 2001 From: James Warner Date: Mon, 3 Aug 2026 14:20:03 +0100 Subject: [PATCH 27/56] additional readme args --- utils/proc_reanalysis/README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/utils/proc_reanalysis/README.md b/utils/proc_reanalysis/README.md index 47d61e2545..1993c1d843 100644 --- a/utils/proc_reanalysis/README.md +++ b/utils/proc_reanalysis/README.md @@ -49,6 +49,8 @@ python process_reanalysis.py \ Required Arguments: `--filepath`: Path to the input reanalysis data. This can be a single file or wildcard expression understood by Iris. +`--cyclestart`: First forecast initialisation time that you want the reanalysis to simulate, in format TZ. +--cycleend: Final forecast initialisation time, inclusive, that you want the reanalysis to simulate, in format TZ. [ncgen-docs]: https://docs.unidata.ucar.edu/nug/current/netcdf_utilities_guide.html From 03303804ae459b22b5fd22c2d5c0bdf0d1761ead Mon Sep 17 00:00:00 2001 From: James Warner Date: Mon, 3 Aug 2026 14:36:24 +0100 Subject: [PATCH 28/56] completed readme --- utils/proc_reanalysis/README.md | 48 ++++++++++++++++++++++++++++----- 1 file changed, 42 insertions(+), 6 deletions(-) diff --git a/utils/proc_reanalysis/README.md b/utils/proc_reanalysis/README.md index 1993c1d843..7c947f2eac 100644 --- a/utils/proc_reanalysis/README.md +++ b/utils/proc_reanalysis/README.md @@ -30,6 +30,9 @@ For each requested forecast cycle, for each variable the script: 5. Preserves the original valid time information. 6. Saves the result as a forecast-style NetCDF file. +> [!TIP] +> This script does not download reanalysis data - a user must fetch this first from archive. + ## Usage The python script requires the Iris package to be present within the python install. @@ -50,15 +53,48 @@ Required Arguments: `--filepath`: Path to the input reanalysis data. This can be a single file or wildcard expression understood by Iris. `--cyclestart`: First forecast initialisation time that you want the reanalysis to simulate, in format TZ. ---cycleend: Final forecast initialisation time, inclusive, that you want the reanalysis to simulate, in format TZ. +`--cycleend`: Final forecast initialisation time, inclusive, that you want the reanalysis to simulate, in format TZ. +`--cyclefreq`: Frequency between forecast cycles, in hours, as an integer. +`--forecastlength`: Length of the forecast you want the reanalysis to simulate, in hours, as an integer. +`--outpath`: Path of where to store the output data. The code will write a file per forecast initialisation, in the format of `reanalysis_%Y%m%dT%H%MZ_.nc`. -[ncgen-docs]: https://docs.unidata.ucar.edu/nug/current/netcdf_utilities_guide.html +## Examples -## Owners +1. A single forecast that goes out to 48h, initialised on the 1st January 2024 at 00Z. -> [!TIP] -> Utilities must have at least one named owner. +``` +python process_reanalysis.py \ + --filepath "/data/era5/*.nc" \ + --cyclestart 20240101T0000Z \ + --cycleend 20240101T0000Z \ + --cyclefreq 6 \ + --forecastlength 48 \ + --outpath /my/output/path/ +``` +Producing one file `my/output/path/reanalysis_20240101T0000Z.nc` + +2. Produce 6-hourly analysis across one day. + +``` +python process_reanalysis.py \ + --filepath "/data/era5/*.nc" \ + --cyclestart 20240101T0000Z \ + --cycleend 20240101T1800Z \ + --cyclefreq 6 \ + --forecastlength 48 \ + --outpath /my/output/path/ +``` +Produces + +``` +/my/output/path/reanalysis_20240101T0000Z.nc +/my/output/path/reanalysis_20240101T0600Z.nc +/my/output/path/reanalysis_20240101T1200Z.nc +/my/output/path/reanalysis_20240101T1800Z.nc +``` + +## Owners The following people should be contacted for queries or issues with this utility: -* [@jfrost-mo](https://github.com/jfrost-mo) +* [@jwarner8](https://github.com/jwarner8) From 11a9ffaa03fa68b170bf7c07ce2063f865d2caef Mon Sep 17 00:00:00 2001 From: James Warner Date: Mon, 3 Aug 2026 14:37:35 +0100 Subject: [PATCH 29/56] indent arg list --- utils/proc_reanalysis/README.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/utils/proc_reanalysis/README.md b/utils/proc_reanalysis/README.md index 7c947f2eac..722170ecf5 100644 --- a/utils/proc_reanalysis/README.md +++ b/utils/proc_reanalysis/README.md @@ -51,12 +51,12 @@ python process_reanalysis.py \ Required Arguments: -`--filepath`: Path to the input reanalysis data. This can be a single file or wildcard expression understood by Iris. -`--cyclestart`: First forecast initialisation time that you want the reanalysis to simulate, in format TZ. -`--cycleend`: Final forecast initialisation time, inclusive, that you want the reanalysis to simulate, in format TZ. -`--cyclefreq`: Frequency between forecast cycles, in hours, as an integer. -`--forecastlength`: Length of the forecast you want the reanalysis to simulate, in hours, as an integer. -`--outpath`: Path of where to store the output data. The code will write a file per forecast initialisation, in the format of `reanalysis_%Y%m%dT%H%MZ_.nc`. +- `--filepath`: Path to the input reanalysis data. This can be a single file or wildcard expression understood by Iris. +- `--cyclestart`: First forecast initialisation time that you want the reanalysis to simulate, in format TZ. +- `--cycleend`: Final forecast initialisation time, inclusive, that you want the reanalysis to simulate, in format TZ. +- `--cyclefreq`: Frequency between forecast cycles, in hours, as an integer. +- `--forecastlength`: Length of the forecast you want the reanalysis to simulate, in hours, as an integer. +- `--outpath`: Path of where to store the output data. The code will write a file per forecast initialisation, in the format of `reanalysis_%Y%m%dT%H%MZ_.nc`. ## Examples From bf3ab586d392441c8f828f9a22553f1b2a2b97f1 Mon Sep 17 00:00:00 2001 From: James Warner Date: Mon, 3 Aug 2026 15:56:01 +0100 Subject: [PATCH 30/56] minor text clarifications and tidy --- utils/proc_reanalysis/README.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/utils/proc_reanalysis/README.md b/utils/proc_reanalysis/README.md index 722170ecf5..8c4477fb45 100644 --- a/utils/proc_reanalysis/README.md +++ b/utils/proc_reanalysis/README.md @@ -10,9 +10,9 @@ The primary motivation for this tool is to be able to analyse reanalysis alongsi - `forecast_period` - `time` -Reanalysis datasets typically only contain valid time `time`, and are a series of files where the `forecast_reference_time` changes every 6 hours, and the `forecast_period` is always zero. +Reanalysis datasets typically only contain valid time `time`, and is a series of files where the `forecast_reference_time` changes every 6 hours, and the `forecast_period` is zero. -This script resolves this issue by transforming reanalysis data into an effective forecast representation. Rather than treating reanalysis as a special data source in CSET, the transformed output can be ingested directly into the CSET workflow. This allows reanalysis to be treated as another "model" within verification systems such as CSET. +This script resolves this issue by transforming reanalysis data into an effective forecast representation. Rather than treating reanalysis as a special data source in CSET, the transformed output can be ingested directly into the CSET workflow. This allows reanalysis to be treated as another "model" within CSET. The script currently supports datasets that can be loaded by Iris and has primarily been developed and tested using: @@ -22,7 +22,7 @@ The script currently supports datasets that can be loaded by Iris and has primar Other model analyses may work, as only the time dimensions are manipulated, but this has not been tested. No scientific changes are made to the meteorological fields themselves. -For each requested forecast cycle, for each variable the script: +For each requested forecast cycle, for each variable found in the reanalysis the script: 1. Extracts the required period of time from reanalysis data. 2. Treats the start of that extraction as a forecast initialisation. 3. Generates a forecast period coordinate. @@ -31,11 +31,11 @@ For each requested forecast cycle, for each variable the script: 6. Saves the result as a forecast-style NetCDF file. > [!TIP] -> This script does not download reanalysis data - a user must fetch this first from archive. +> This script does not download reanalysis data - a user must fetch this first from archive/api. ## Usage -The python script requires the Iris package to be present within the python install. +The python script requires the Iris package to be installed and available to python. Run it with: @@ -51,7 +51,7 @@ python process_reanalysis.py \ Required Arguments: -- `--filepath`: Path to the input reanalysis data. This can be a single file or wildcard expression understood by Iris. +- `--filepath`: Path to the input reanalysis data. This can be a single file or wildcard expression understood by Iris. If a wildcard is used, then quote the input to prevent the shell expanding the filelist as arguments to python. - `--cyclestart`: First forecast initialisation time that you want the reanalysis to simulate, in format TZ. - `--cycleend`: Final forecast initialisation time, inclusive, that you want the reanalysis to simulate, in format TZ. - `--cyclefreq`: Frequency between forecast cycles, in hours, as an integer. From 04c6a5a7f0ab34537a87abaae4a20481f72a0140 Mon Sep 17 00:00:00 2001 From: James Warner Date: Mon, 3 Aug 2026 15:56:47 +0100 Subject: [PATCH 31/56] refactor func location --- utils/proc_reanalysis/process_reanalysis.py | 74 ++++++++++----------- 1 file changed, 37 insertions(+), 37 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index 93c19cd3b2..a48938d1dc 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -17,6 +17,42 @@ from datetime import datetime, timedelta +def _identify_number_of_cycles_required( + cyclestart: str, cycleend: str, cyclefreq: str +) -> list: + """Generate forecast initialisation datetimes between two cycle bounds. + + Parameters + ---------- + cyclestart : str + First forecast cycle in YYYYMMDDTHHMMZ format. + cycleend : str + Last forecast cycle in YYYYMMDDTHHMMZ format. + cyclefreq : int + Frequency between forecast cycles in hours. + + Returns + ------- + forecast_initialisations: list + Forecast initialisation datetimes from cyclestart to cycleend, + inclusive, separated by cyclefreq hours. + """ + # Identify the start and end times, and create datetime objects for these + start_dt = datetime.strptime(cyclestart, "%Y%m%dT%H%MZ") + end_dt = datetime.strptime(cycleend, "%Y%m%dT%H%MZ") + + # To store initialisation times + forecast_initialisations = [] + current = start_dt + + # Iterate over all initiations within the bounds, using the cyclefreq to determine interval. + while current <= end_dt: + forecast_initialisations.append(current) + current += timedelta(hours=cyclefreq) + + return forecast_initialisations + + def _create_forecasts( reanalysis: iris.cube.CubeList, forecast_initialisations: list, @@ -167,42 +203,6 @@ def _create_forecasts( raise ValueError("No suitable cubes found for saving!") -def _identify_number_of_cycles_required( - cyclestart: str, cycleend: str, cyclefreq: str -) -> list: - """Generate forecast initialisation datetimes between two cycle bounds. - - Parameters - ---------- - cyclestart : str - First forecast cycle in YYYYMMDDTHHMMZ format. - cycleend : str - Last forecast cycle in YYYYMMDDTHHMMZ format. - cyclefreq : int - Frequency between forecast cycles in hours. - - Returns - ------- - forecast_initialisations: list - Forecast initialisation datetimes from cyclestart to cycleend, - inclusive, separated by cyclefreq hours. - """ - # Identify the start and end times, and create datetime objects for these - start_dt = datetime.strptime(cyclestart, "%Y%m%dT%H%MZ") - end_dt = datetime.strptime(cycleend, "%Y%m%dT%H%MZ") - - # To store initialisation times - forecast_initialisations = [] - current = start_dt - - # Iterate over all initiations within the bounds, using the cyclefreq to determine interval. - while current <= end_dt: - forecast_initialisations.append(current) - current += timedelta(hours=cyclefreq) - - return forecast_initialisations - - def main() -> None: """Generate forecast-like datasets from reanalysis data. @@ -210,7 +210,7 @@ def main() -> None: process the input reanalysis data, and write the resulting forecast files to disk. """ - parser = argparse.ArgumentParser(description="Process forecast data.") + parser = argparse.ArgumentParser(description="Process arguments.") parser.add_argument("--filepath", required=True, help="Path to file(s)") parser.add_argument( From 725c1b843716b111b8b6fcd5da14af7af4371bbc Mon Sep 17 00:00:00 2001 From: James Frost Date: Tue, 4 Aug 2026 18:35:26 +0100 Subject: [PATCH 32/56] Add utils folder to ruff linting and formatting --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index ee24b4a0fa..06615a1a3b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -83,7 +83,7 @@ skip = "build,*.css,*.ipynb,*.js,*.html,*.svg,*.xml,.git" [tool.ruff] line-length = 88 -src = ["src", "test"] +src = ["src", "test", "utils"] [tool.ruff.lint] extend-select = ["B", "D"] From 20615a5dc02a5924378e27744c106f911a438d1a Mon Sep 17 00:00:00 2001 From: James Frost Date: Tue, 4 Aug 2026 18:41:06 +0100 Subject: [PATCH 33/56] Add utils to testpaths This will cause any tests in the community utilities folder to be run. --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 06615a1a3b..1c264b000e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -66,7 +66,7 @@ markers = [ ] minversion = "7" pythonpath = ["src"] -testpaths = ["tests"] +testpaths = ["tests", "utils"] [tool.coverage.run] branch = true From 8c5187002b17b698a0a2d3728c75b6551eba71da Mon Sep 17 00:00:00 2001 From: James Warner Date: Fri, 7 Aug 2026 08:34:08 +0100 Subject: [PATCH 34/56] switch to using isoformat --- utils/proc_reanalysis/process_reanalysis.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index a48938d1dc..4841fa5ab3 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -38,8 +38,8 @@ def _identify_number_of_cycles_required( inclusive, separated by cyclefreq hours. """ # Identify the start and end times, and create datetime objects for these - start_dt = datetime.strptime(cyclestart, "%Y%m%dT%H%MZ") - end_dt = datetime.strptime(cycleend, "%Y%m%dT%H%MZ") + start_dt = datetime.fromisoformat(cyclestart) + end_dt = datetime.fromisoformat(cycleend) # To store initialisation times forecast_initialisations = [] From 54207477d2cc73c97461a4ee1eacb54648c653e8 Mon Sep 17 00:00:00 2001 From: James Warner Date: Fri, 7 Aug 2026 08:38:35 +0100 Subject: [PATCH 35/56] switch to passing datetime rather than str --- utils/proc_reanalysis/process_reanalysis.py | 25 +++++++++++---------- 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index 4841fa5ab3..bbfbad2865 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -18,29 +18,25 @@ def _identify_number_of_cycles_required( - cyclestart: str, cycleend: str, cyclefreq: str + start_dt: datetime, end_dt: datetime, cyclefreq: timedelta ) -> list: """Generate forecast initialisation datetimes between two cycle bounds. Parameters ---------- - cyclestart : str - First forecast cycle in YYYYMMDDTHHMMZ format. - cycleend : str - Last forecast cycle in YYYYMMDDTHHMMZ format. - cyclefreq : int + cyclestart : datetime + First forecast cycle time. + cycleend : datetime + Last forecast cycle time. + cyclefreq : timedelta Frequency between forecast cycles in hours. Returns ------- forecast_initialisations: list - Forecast initialisation datetimes from cyclestart to cycleend, + Forecast initialisation datetimes from start_dt to end_dt, inclusive, separated by cyclefreq hours. """ - # Identify the start and end times, and create datetime objects for these - start_dt = datetime.fromisoformat(cyclestart) - end_dt = datetime.fromisoformat(cycleend) - # To store initialisation times forecast_initialisations = [] current = start_dt @@ -248,9 +244,14 @@ def main() -> None: print() print("Starting process_reanalysis.py...") + # Identify the start and end times, and create datetime objects for these + start_dt = datetime.fromisoformat(cyclestart) + end_dt = datetime.fromisoformat(cycleend) + cyclefreq = timedelta(cyclefreq) + # Get all forecast initiations forecast_initialisations = _identify_number_of_cycles_required( - cyclestart, cycleend, cyclefreq + start_dt, end_dt, cyclefreq ) # Load all reanalysis supplied From 42580a9df0195f2258b81716bdafc195d0549f71 Mon Sep 17 00:00:00 2001 From: James Warner Date: Fri, 7 Aug 2026 08:42:48 +0100 Subject: [PATCH 36/56] add datetime parsing as command line arg --- utils/proc_reanalysis/process_reanalysis.py | 22 ++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index bbfbad2865..35d6218519 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -18,7 +18,7 @@ def _identify_number_of_cycles_required( - start_dt: datetime, end_dt: datetime, cyclefreq: timedelta + cyclestart: datetime, cycleend: datetime, cyclefreq: timedelta ) -> list: """Generate forecast initialisation datetimes between two cycle bounds. @@ -39,10 +39,10 @@ def _identify_number_of_cycles_required( """ # To store initialisation times forecast_initialisations = [] - current = start_dt + current = cyclestart # Iterate over all initiations within the bounds, using the cyclefreq to determine interval. - while current <= end_dt: + while current <= cycleend: forecast_initialisations.append(current) current += timedelta(hours=cyclefreq) @@ -210,16 +210,22 @@ def main() -> None: parser.add_argument("--filepath", required=True, help="Path to file(s)") parser.add_argument( - "--cyclestart", type=str, required=True, help="First forecast initiation/cycle" + "--cyclestart", + type=datetime.fromisoformat, + required=True, + help="First forecast initiation/cycle, in format %Y%m%dT%H%MZ", ) parser.add_argument( - "--cycleend", type=str, required=True, help="Final forecast initiation/cycle" + "--cycleend", + type=datetime.fromisoformat, + required=True, + help="Final forecast initiation/cycle, in format %Y%m%dT%H%MZ", ) parser.add_argument( "--cyclefreq", type=int, required=True, - help="Time between forecast initiationss/cycles", + help="Time between forecast initiations/cycles", ) parser.add_argument( "--forecastlength", @@ -245,13 +251,11 @@ def main() -> None: print("Starting process_reanalysis.py...") # Identify the start and end times, and create datetime objects for these - start_dt = datetime.fromisoformat(cyclestart) - end_dt = datetime.fromisoformat(cycleend) cyclefreq = timedelta(cyclefreq) # Get all forecast initiations forecast_initialisations = _identify_number_of_cycles_required( - start_dt, end_dt, cyclefreq + cyclestart, cycleend, cyclefreq ) # Load all reanalysis supplied From d5c82de4f193df8ea7fee3513570ffb7854d3cca Mon Sep 17 00:00:00 2001 From: James Warner Date: Fri, 7 Aug 2026 08:44:29 +0100 Subject: [PATCH 37/56] specify cycle freq units --- utils/proc_reanalysis/process_reanalysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index 35d6218519..2030361b98 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -225,7 +225,7 @@ def main() -> None: "--cyclefreq", type=int, required=True, - help="Time between forecast initiations/cycles", + help="Hours between forecast initiations/cycles", ) parser.add_argument( "--forecastlength", From 604388d6b2c62a289331f4a9927679f55546ec49 Mon Sep 17 00:00:00 2001 From: James Warner Date: Fri, 7 Aug 2026 08:45:48 +0100 Subject: [PATCH 38/56] directly convert str of cyclefreq to timedelta --- utils/proc_reanalysis/process_reanalysis.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index 2030361b98..2dd8234fe2 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -243,16 +243,13 @@ def main() -> None: filepath = args.filepath cyclestart = args.cyclestart cycleend = args.cycleend - cyclefreq = args.cyclefreq + cyclefreq = timedelta(hours=args.cyclefreq) forecastlength = args.forecastlength outpath = args.outpath print() print("Starting process_reanalysis.py...") - # Identify the start and end times, and create datetime objects for these - cyclefreq = timedelta(cyclefreq) - # Get all forecast initiations forecast_initialisations = _identify_number_of_cycles_required( cyclestart, cycleend, cyclefreq From fadf4d5a0467be02f433e57bb5214526a0f08a3b Mon Sep 17 00:00:00 2001 From: James Warner Date: Fri, 7 Aug 2026 08:49:11 +0100 Subject: [PATCH 39/56] tidyup from code review --- utils/proc_reanalysis/process_reanalysis.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index 2dd8234fe2..a8a7cd96e5 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -176,7 +176,7 @@ def _create_forecasts( units=time_coord.units, ) - # Add this dimsnion to the cube, tied to the forecast_period dimension. + # Add this dimension to the cube, tied to the forecast_period dimension. cube_slice.add_aux_coord( new_time_coord, data_dims=(cube_slice.coord_dims("forecast_period")[0],), @@ -188,13 +188,9 @@ def _create_forecasts( # Once all cubes processed, save to disk. if len(cutouts) > 0: - print( - f"Saving {outpath + '/reanalysis_' + start.strftime('%Y%m%dT%H%MZ')}.nc" - ) - iris.save( - cutouts, - outpath + "/reanalysis_" + start.strftime("%Y%m%dT%H%MZ") + ".nc", - ) + filename = f"{outpath}/reanalysis_{start.strftime('%Y%m%dT%H%MZ').nc}" + print(f"Saving {filename}") + iris.save(cutouts, filename) else: raise ValueError("No suitable cubes found for saving!") From 361d410bd403e56840a12df992fdb564e5461f87 Mon Sep 17 00:00:00 2001 From: James Warner Date: Fri, 7 Aug 2026 08:54:15 +0100 Subject: [PATCH 40/56] apply changes from code review --- utils/proc_reanalysis/process_reanalysis.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index a8a7cd96e5..90b1da1594 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -17,7 +17,7 @@ from datetime import datetime, timedelta -def _identify_number_of_cycles_required( +def identify_number_of_cycles_required( cyclestart: datetime, cycleend: datetime, cyclefreq: timedelta ) -> list: """Generate forecast initialisation datetimes between two cycle bounds. @@ -49,7 +49,7 @@ def _identify_number_of_cycles_required( return forecast_initialisations -def _create_forecasts( +def create_forecasts( reanalysis: iris.cube.CubeList, forecast_initialisations: list, forecastlength: int, @@ -125,7 +125,7 @@ def _create_forecasts( if item in cube_slice.attributes: del cube_slice.attributes[item] - # Get a copy of time coord, and dimension this corresponds to. + # Get a copy of time coord time_coord = cube_slice.coord("time") # Work out units of forecast_period and adjust if necessary. @@ -145,8 +145,8 @@ def _create_forecasts( else: raise ValueError(f"Unhandled time units: {time_coord.units}") - # Get copy of time coordinate - time_coord_points = time_coord.points.copy() + # Get time points and dimension that time corresponds to. + time_coord_points = time_coord.points time_dim = cube_slice.coord_dims("time")[0] # Create forecast period dimension @@ -247,7 +247,7 @@ def main() -> None: print("Starting process_reanalysis.py...") # Get all forecast initiations - forecast_initialisations = _identify_number_of_cycles_required( + forecast_initialisations = identify_number_of_cycles_required( cyclestart, cycleend, cyclefreq ) @@ -260,7 +260,7 @@ def main() -> None: print() print("Creating postprocessed files...") - _create_forecasts(reanalysis, forecast_initialisations, forecastlength, outpath) + create_forecasts(reanalysis, forecast_initialisations, forecastlength, outpath) print("Done") From 450467f68ba0039b6f3b78bec255c0d8beb0990d Mon Sep 17 00:00:00 2001 From: James Warner Date: Fri, 7 Aug 2026 08:56:05 +0100 Subject: [PATCH 41/56] update tests given changes to function names --- .../test_process_reanalysis.py | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/utils/proc_reanalysis/test_process_reanalysis.py b/utils/proc_reanalysis/test_process_reanalysis.py index 969748aae4..d2025739ca 100644 --- a/utils/proc_reanalysis/test_process_reanalysis.py +++ b/utils/proc_reanalysis/test_process_reanalysis.py @@ -12,7 +12,7 @@ def test_single_cycle(): """Assert single datetime returned if one initialisation.""" - result = proc_reanalysis._identify_number_of_cycles_required( + result = proc_reanalysis.identify_number_of_cycles_required( "20240101T0000Z", "20240101T0000Z", 6, @@ -23,7 +23,7 @@ def test_single_cycle(): def test_multiple_cycles(): """Test handling of multiple cycles identified.""" - result = proc_reanalysis._identify_number_of_cycles_required( + result = proc_reanalysis.identify_number_of_cycles_required( "20240101T0000Z", "20240101T1200Z", 6, @@ -38,7 +38,7 @@ def test_multiple_cycles(): def test_non_divisible_interval(): """Check end point is not exceeded.""" - result = proc_reanalysis._identify_number_of_cycles_required( + result = proc_reanalysis.identify_number_of_cycles_required( "20240101T0000Z", "20240101T1000Z", 6, @@ -53,7 +53,7 @@ def test_non_divisible_interval(): def test_invalid_date_format(): """Raise problem with input.""" with pytest.raises(ValueError): - proc_reanalysis._identify_number_of_cycles_required( + proc_reanalysis.identify_number_of_cycles_required( "2024-01-01", "20240101T1200Z", 6, @@ -93,7 +93,7 @@ def test_forecast_period_created(tmp_path): """Check forecast period constructed correctly.""" cube = make_cube() - proc_reanalysis._create_forecasts( + proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1, 0)], forecastlength=4, @@ -138,7 +138,7 @@ def test_minutes_converted_to_hours(tmp_path): cube.coord("time").points = [0, 60, 120, 180, 240] - proc_reanalysis._create_forecasts( + proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], forecastlength=2, @@ -157,7 +157,7 @@ def test_unknown_time_units_raise(tmp_path): cube = make_cube(units="days since 2024-01-01 00:00:00") with pytest.raises(ValueError, match="Unhandled time units"): - proc_reanalysis._create_forecasts( + proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], forecastlength=1, @@ -171,7 +171,7 @@ def test_forecast_reference_time_created(tmp_path): cube = make_cube() - proc_reanalysis._create_forecasts( + proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [init_time], forecastlength=4, @@ -194,7 +194,7 @@ def test_forecast_attributes_removed(tmp_path): cube.attributes["source"] = "ERA5" cube.attributes["um_version"] = "13.0" - proc_reanalysis._create_forecasts( + proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], forecastlength=4, @@ -214,7 +214,7 @@ def test_cube_skipped_if_insufficient_data(tmp_path): cube = make_cube() with pytest.raises(ValueError, match="No suitable cubes found for saving!"): - proc_reanalysis._create_forecasts( + proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], forecastlength=10, @@ -227,7 +227,7 @@ def test_multiple_cubes_processed(tmp_path): cube1 = make_cube("air_temperature") cube2 = make_cube("air_pressure") - proc_reanalysis._create_forecasts( + proc_reanalysis.create_forecasts( iris.cube.CubeList([cube1, cube2]), [datetime(2024, 1, 1)], forecastlength=4, From 9cd6c6f428719ae07f68640f0f253defa057eff8 Mon Sep 17 00:00:00 2001 From: James Warner Date: Fri, 7 Aug 2026 09:46:00 +0100 Subject: [PATCH 42/56] update refs to using new ISO datetime standard --- utils/proc_reanalysis/README.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/utils/proc_reanalysis/README.md b/utils/proc_reanalysis/README.md index 8c4477fb45..b6fbfe9357 100644 --- a/utils/proc_reanalysis/README.md +++ b/utils/proc_reanalysis/README.md @@ -42,8 +42,8 @@ Run it with: ``` python process_reanalysis.py \ --filepath "" \ - --cyclestart YYYYMMDDTHHMMZ \ - --cycleend YYYYMMDDTHHMMZ \ + --cyclestart "YYYY-MM-DD HH:MM:SS" \ + --cycleend "YYYY-MM-DD HH:MM:SS" \ --cyclefreq \ --forecastlength \ --outpath "" @@ -65,8 +65,8 @@ Required Arguments: ``` python process_reanalysis.py \ --filepath "/data/era5/*.nc" \ - --cyclestart 20240101T0000Z \ - --cycleend 20240101T0000Z \ + --cyclestart "2024-01-01 00:00:00" \ + --cycleend "2024-01-01 00:00:00" \ --cyclefreq 6 \ --forecastlength 48 \ --outpath /my/output/path/ @@ -78,8 +78,8 @@ Producing one file `my/output/path/reanalysis_20240101T0000Z.nc` ``` python process_reanalysis.py \ --filepath "/data/era5/*.nc" \ - --cyclestart 20240101T0000Z \ - --cycleend 20240101T1800Z \ + --cyclestart "2024-01-01 00:00:00" \ + --cycleend "2024-01-01 18:00:00" \ --cyclefreq 6 \ --forecastlength 48 \ --outpath /my/output/path/ From 6f9b667f80adf10b72f73f7f76ea28ed832f106e Mon Sep 17 00:00:00 2001 From: James Warner Date: Fri, 7 Aug 2026 09:47:01 +0100 Subject: [PATCH 43/56] update identify_number_of_cycles_required to use timedelta --- utils/proc_reanalysis/process_reanalysis.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index 90b1da1594..bdc582eb07 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -44,7 +44,7 @@ def identify_number_of_cycles_required( # Iterate over all initiations within the bounds, using the cyclefreq to determine interval. while current <= cycleend: forecast_initialisations.append(current) - current += timedelta(hours=cyclefreq) + current += cyclefreq return forecast_initialisations @@ -188,7 +188,7 @@ def create_forecasts( # Once all cubes processed, save to disk. if len(cutouts) > 0: - filename = f"{outpath}/reanalysis_{start.strftime('%Y%m%dT%H%MZ').nc}" + filename = f"{outpath}/reanalysis_{start.strftime('%Y%m%dT%H%MZ')}.nc" print(f"Saving {filename}") iris.save(cutouts, filename) else: @@ -204,18 +204,18 @@ def main() -> None: """ parser = argparse.ArgumentParser(description="Process arguments.") - parser.add_argument("--filepath", required=True, help="Path to file(s)") + parser.add_argument("--filepath", required=True, help="Path to file(s) quoted") parser.add_argument( "--cyclestart", type=datetime.fromisoformat, required=True, - help="First forecast initiation/cycle, in format %Y%m%dT%H%MZ", + help="First forecast initiation/cycle, in format YYYY-MM-DD HH:MM:SS, quoted", ) parser.add_argument( "--cycleend", type=datetime.fromisoformat, required=True, - help="Final forecast initiation/cycle, in format %Y%m%dT%H%MZ", + help="Final forecast initiation/cycle, in format YYYY-MM-DD HH:MM:SS, quoted", ) parser.add_argument( "--cyclefreq", From 7d2fc047e786bfcefbd27e7aacb86460dc3ae9c2 Mon Sep 17 00:00:00 2001 From: James Warner Date: Fri, 7 Aug 2026 09:47:27 +0100 Subject: [PATCH 44/56] update tests to use new datetime format --- .../test_process_reanalysis.py | 32 +++++++------------ 1 file changed, 11 insertions(+), 21 deletions(-) diff --git a/utils/proc_reanalysis/test_process_reanalysis.py b/utils/proc_reanalysis/test_process_reanalysis.py index d2025739ca..47e782ef85 100644 --- a/utils/proc_reanalysis/test_process_reanalysis.py +++ b/utils/proc_reanalysis/test_process_reanalysis.py @@ -1,6 +1,6 @@ """Unit tests for process_reanalysis.py.""" -from datetime import datetime +from datetime import datetime, timedelta import iris import numpy as np @@ -13,9 +13,9 @@ def test_single_cycle(): """Assert single datetime returned if one initialisation.""" result = proc_reanalysis.identify_number_of_cycles_required( - "20240101T0000Z", - "20240101T0000Z", - 6, + datetime.fromisoformat("2024-01-01 00:00:00"), + datetime.fromisoformat("2024-01-01 00:00:00"), + timedelta(hours=6), ) assert result == [datetime(2024, 1, 1, 0, 0)] @@ -24,9 +24,9 @@ def test_single_cycle(): def test_multiple_cycles(): """Test handling of multiple cycles identified.""" result = proc_reanalysis.identify_number_of_cycles_required( - "20240101T0000Z", - "20240101T1200Z", - 6, + datetime.fromisoformat("2024-01-01 00:00:00"), + datetime.fromisoformat("2024-01-01 12:00:00"), + timedelta(hours=6), ) assert result == [ @@ -39,9 +39,9 @@ def test_multiple_cycles(): def test_non_divisible_interval(): """Check end point is not exceeded.""" result = proc_reanalysis.identify_number_of_cycles_required( - "20240101T0000Z", - "20240101T1000Z", - 6, + datetime.fromisoformat("2024-01-01 00:00:00"), + datetime.fromisoformat("2024-01-01 10:00:00"), + timedelta(hours=6), ) assert result == [ @@ -50,16 +50,6 @@ def test_non_divisible_interval(): ] -def test_invalid_date_format(): - """Raise problem with input.""" - with pytest.raises(ValueError): - proc_reanalysis.identify_number_of_cycles_required( - "2024-01-01", - "20240101T1200Z", - 6, - ) - - def make_cube( name="air_temperature", units="hours since 2024-01-01 00:00:00", @@ -118,7 +108,7 @@ def test_seconds_converted_to_hours(tmp_path): cube.coord("time").points = [0, 3600, 7200, 10800, 14400] - proc_reanalysis._create_forecasts( + proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], forecastlength=2, From 8aa62455133887ef653d9ad6e256593fee1dc0f9 Mon Sep 17 00:00:00 2001 From: James Warner Date: Fri, 7 Aug 2026 09:50:49 +0100 Subject: [PATCH 45/56] update import for github actions --- utils/proc_reanalysis/test_process_reanalysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/proc_reanalysis/test_process_reanalysis.py b/utils/proc_reanalysis/test_process_reanalysis.py index 47e782ef85..3e892eda24 100644 --- a/utils/proc_reanalysis/test_process_reanalysis.py +++ b/utils/proc_reanalysis/test_process_reanalysis.py @@ -4,10 +4,10 @@ import iris import numpy as np -import process_reanalysis as proc_reanalysis import pytest from iris.coords import DimCoord from iris.cube import Cube +from utils.proc_reanalysis import process_reanalysis as proc_reanalysis def test_single_cycle(): From d4b5f1b900ee64402e963ac516733405c175c4d0 Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 19 Aug 2026 13:18:14 +0100 Subject: [PATCH 46/56] add new entry and generalise function name, tidy up docstring --- src/CSET/operators/read.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/CSET/operators/read.py b/src/CSET/operators/read.py index 96791d383d..fe64eecbc4 100644 --- a/src/CSET/operators/read.py +++ b/src/CSET/operators/read.py @@ -413,7 +413,7 @@ def _loading_callback(cube: iris.cube.Cube, field, filename: str) -> iris.cube.C _lfric_time_callback(cube) _lfric_forecast_period_callback(cube) cube = _fix_no_time_coords_callback(cube) - _normalise_ML_varname(cube) + _normalise_longname(cube) return cube @@ -1098,8 +1098,8 @@ def _fix_no_time_coords_callback(cube: iris.cube.Cube): return cube -def _normalise_ML_varname(cube: iris.cube.Cube): - """Fix plev variable names to standard names.""" +def _normalise_longname(cube: iris.cube.Cube): + """Normalise long_name to the LFRic standard list.""" if cube.coords("pressure"): if cube.name() == "x_wind": cube.long_name = "zonal_wind_at_pressure_levels" @@ -1116,6 +1116,8 @@ def _normalise_ML_varname(cube: iris.cube.Cube): cube.long_name = "eastward_wind_at_10m" if cube.name() == "y_wind" and cube.var_name == "v_wind_at_10m": cube.long_name = "northward_wind_at_10m" + if cube.name() == "air_pressure_at_sea_level": + cube.long_name = "air_pressure_at_mean_sea_level" def _check_combine_point_observations(cubes: iris.cube.CubeList): From 10183835385ecd4357a51ed86b129818d4ff65a5 Mon Sep 17 00:00:00 2001 From: James Warner <62252918+jwarner8@users.noreply.github.com> Date: Wed, 19 Aug 2026 13:28:00 +0100 Subject: [PATCH 47/56] Add shebang utils/proc_reanalysis/process_reanalysis.py Co-authored-by: James Frost --- utils/proc_reanalysis/process_reanalysis.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index bdc582eb07..f43929c0b4 100644 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -1,3 +1,5 @@ +#!/usr/bin/python3 + """ Code that restructures reanalysis data to give it an effective forecast_period. From ef09c9b94fd869ff2e87b1f38c58005ca353abbb Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 19 Aug 2026 13:29:46 +0100 Subject: [PATCH 48/56] make executable --- utils/proc_reanalysis/process_reanalysis.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 utils/proc_reanalysis/process_reanalysis.py diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py old mode 100644 new mode 100755 From 80860bfe055bd25c07b108acb5190c5eac24df97 Mon Sep 17 00:00:00 2001 From: James Warner <62252918+jwarner8@users.noreply.github.com> Date: Wed, 19 Aug 2026 13:30:45 +0100 Subject: [PATCH 49/56] Update utils/proc_reanalysis/process_reanalysis.py Co-authored-by: James Frost --- utils/proc_reanalysis/process_reanalysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index f43929c0b4..b9a4b9b865 100755 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -229,7 +229,7 @@ def main() -> None: "--forecastlength", type=int, required=True, - help="Forecast length in SI units i.e. PT48H", + help="Forecast length in hours, i.e. 48", ) parser.add_argument( "--outpath", type=str, required=True, help="Where to write output data" From fd184413c101eff1e66010971b2b0dbdc1d8ed6d Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 19 Aug 2026 13:32:55 +0100 Subject: [PATCH 50/56] use timedelta for forecast length --- utils/proc_reanalysis/process_reanalysis.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index f43929c0b4..97fe8bf348 100755 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -54,7 +54,7 @@ def identify_number_of_cycles_required( def create_forecasts( reanalysis: iris.cube.CubeList, forecast_initialisations: list, - forecastlength: int, + forecastlength: datetime.timedelta, outpath: str, ) -> None: """Create forecast files from reanalysis data. @@ -72,7 +72,7 @@ def create_forecasts( will be extracted. forecast_initialisations: list Forecast initialisation times to process. - forecastlength: int + forecastlength: datetime.timedelta Forecast length in hours. outpath: str Directory to which the generated forecast files will be saved. @@ -83,11 +83,13 @@ def create_forecasts( """ # Iterate over all forecast initialisations sequentially. for forecast in forecast_initialisations: - print(f"Working on forecast initialisation {forecast} out to {forecastlength}H") + print( + f"Working on forecast initialisation {forecast} out to {forecastlength.total_seconds() / 3600}H" + ) # Work out start and end time start = forecast - end = forecast + timedelta(hours=forecastlength) + end = forecast + forecastlength cutouts = iris.cube.CubeList() @@ -242,7 +244,7 @@ def main() -> None: cyclestart = args.cyclestart cycleend = args.cycleend cyclefreq = timedelta(hours=args.cyclefreq) - forecastlength = args.forecastlength + forecastlength = timedelta(hours=args.forecastlength) outpath = args.outpath print() From b9bebe14dfe2c5cddb1f92f6c5cc2736028e5bb2 Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 19 Aug 2026 13:35:07 +0100 Subject: [PATCH 51/56] update argument names --- utils/proc_reanalysis/README.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/utils/proc_reanalysis/README.md b/utils/proc_reanalysis/README.md index b6fbfe9357..cc2603db9c 100644 --- a/utils/proc_reanalysis/README.md +++ b/utils/proc_reanalysis/README.md @@ -41,7 +41,7 @@ Run it with: ``` python process_reanalysis.py \ - --filepath "" \ + --files "" \ --cyclestart "YYYY-MM-DD HH:MM:SS" \ --cycleend "YYYY-MM-DD HH:MM:SS" \ --cyclefreq \ @@ -51,7 +51,7 @@ python process_reanalysis.py \ Required Arguments: -- `--filepath`: Path to the input reanalysis data. This can be a single file or wildcard expression understood by Iris. If a wildcard is used, then quote the input to prevent the shell expanding the filelist as arguments to python. +- `--files`: Path to the input reanalysis data. This can be a single file or wildcard expression understood by Iris. If a wildcard is used, then quote the input to prevent the shell expanding the filelist as arguments to python. - `--cyclestart`: First forecast initialisation time that you want the reanalysis to simulate, in format TZ. - `--cycleend`: Final forecast initialisation time, inclusive, that you want the reanalysis to simulate, in format TZ. - `--cyclefreq`: Frequency between forecast cycles, in hours, as an integer. @@ -64,7 +64,7 @@ Required Arguments: ``` python process_reanalysis.py \ - --filepath "/data/era5/*.nc" \ + --files "/data/era5/*.nc" \ --cyclestart "2024-01-01 00:00:00" \ --cycleend "2024-01-01 00:00:00" \ --cyclefreq 6 \ @@ -77,7 +77,7 @@ Producing one file `my/output/path/reanalysis_20240101T0000Z.nc` ``` python process_reanalysis.py \ - --filepath "/data/era5/*.nc" \ + --files "/data/era5/*.nc" \ --cyclestart "2024-01-01 00:00:00" \ --cycleend "2024-01-01 18:00:00" \ --cyclefreq 6 \ From 53270d5b17837869e4ac0b301c57f67617175288 Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 19 Aug 2026 13:35:25 +0100 Subject: [PATCH 52/56] update argument names to script --- utils/proc_reanalysis/process_reanalysis.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index 97fe8bf348..933809718e 100755 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -208,7 +208,7 @@ def main() -> None: """ parser = argparse.ArgumentParser(description="Process arguments.") - parser.add_argument("--filepath", required=True, help="Path to file(s) quoted") + parser.add_argument("--files", required=True, help="Path to file(s) quoted") parser.add_argument( "--cyclestart", type=datetime.fromisoformat, @@ -240,7 +240,7 @@ def main() -> None: args = parser.parse_args() # Populate required variables - filepath = args.filepath + filepath = args.files cyclestart = args.cyclestart cycleend = args.cycleend cyclefreq = timedelta(hours=args.cyclefreq) From 218289a9d6ba3ab801b3f7f689139e55476f0062 Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 19 Aug 2026 13:45:33 +0100 Subject: [PATCH 53/56] use standard cycle time iso time format --- utils/proc_reanalysis/README.md | 12 ++++++------ utils/proc_reanalysis/process_reanalysis.py | 4 ++-- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/utils/proc_reanalysis/README.md b/utils/proc_reanalysis/README.md index cc2603db9c..7ea2571a6e 100644 --- a/utils/proc_reanalysis/README.md +++ b/utils/proc_reanalysis/README.md @@ -42,8 +42,8 @@ Run it with: ``` python process_reanalysis.py \ --files "" \ - --cyclestart "YYYY-MM-DD HH:MM:SS" \ - --cycleend "YYYY-MM-DD HH:MM:SS" \ + --cyclestart YYYYMMDDTHHMMZ" \ + --cycleend YYYYMMDDTHHMMZ" \ --cyclefreq \ --forecastlength \ --outpath "" @@ -65,8 +65,8 @@ Required Arguments: ``` python process_reanalysis.py \ --files "/data/era5/*.nc" \ - --cyclestart "2024-01-01 00:00:00" \ - --cycleend "2024-01-01 00:00:00" \ + --cyclestart "20240101T0000Z" \ + --cycleend "20240101T0000Z" \ --cyclefreq 6 \ --forecastlength 48 \ --outpath /my/output/path/ @@ -78,8 +78,8 @@ Producing one file `my/output/path/reanalysis_20240101T0000Z.nc` ``` python process_reanalysis.py \ --files "/data/era5/*.nc" \ - --cyclestart "2024-01-01 00:00:00" \ - --cycleend "2024-01-01 18:00:00" \ + --cyclestart "20240101T0000Z" \ + --cycleend "20240101T1800Z" \ --cyclefreq 6 \ --forecastlength 48 \ --outpath /my/output/path/ diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index d284c810e4..40e9cf1c76 100755 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -213,13 +213,13 @@ def main() -> None: "--cyclestart", type=datetime.fromisoformat, required=True, - help="First forecast initiation/cycle, in format YYYY-MM-DD HH:MM:SS, quoted", + help="First forecast initiation/cycle, in format YYYYMMDDTHHMMZ", ) parser.add_argument( "--cycleend", type=datetime.fromisoformat, required=True, - help="Final forecast initiation/cycle, in format YYYY-MM-DD HH:MM:SS, quoted", + help="Final forecast initiation/cycle, in format YYYYMMDDTHHMMZ", ) parser.add_argument( "--cyclefreq", From d6344b25c8fc62f0bf7d8efbeb35477b0ec7bbfa Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 19 Aug 2026 14:04:16 +0100 Subject: [PATCH 54/56] update test --- tests/operators/test_read.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/operators/test_read.py b/tests/operators/test_read.py index ad1b793045..c59a395809 100644 --- a/tests/operators/test_read.py +++ b/tests/operators/test_read.py @@ -1166,11 +1166,11 @@ def test_fix_no_time_coords_callback(cube): assert cube.coord("time").units == "hours since 0001-01-01 00:00:00" -def test_normalise_ML_varname(transect_source_cube): +def test_normalise_longname(transect_source_cube): """Check that pressure varname is changed.""" cube = transect_source_cube.copy() cube.rename = "air_temperature" - read._normalise_ML_varname(cube) + read._normalise_longname(cube) assert cube.long_name == "temperature_at_pressure_levels" From 0f1bf8f6d28e0b4594329bd260ac6172a0c28ff4 Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 19 Aug 2026 14:46:45 +0100 Subject: [PATCH 55/56] fix object types in tests --- utils/proc_reanalysis/test_process_reanalysis.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/utils/proc_reanalysis/test_process_reanalysis.py b/utils/proc_reanalysis/test_process_reanalysis.py index 3e892eda24..cec5b4be54 100644 --- a/utils/proc_reanalysis/test_process_reanalysis.py +++ b/utils/proc_reanalysis/test_process_reanalysis.py @@ -86,7 +86,7 @@ def test_forecast_period_created(tmp_path): proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1, 0)], - forecastlength=4, + forecastlength=timedelta(hours=4), outpath=str(tmp_path), ) @@ -111,7 +111,7 @@ def test_seconds_converted_to_hours(tmp_path): proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], - forecastlength=2, + forecastlength=timedelta(hours=2), outpath=str(tmp_path), ) @@ -131,7 +131,7 @@ def test_minutes_converted_to_hours(tmp_path): proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], - forecastlength=2, + forecastlength=timedelta(hours=2), outpath=str(tmp_path), ) @@ -150,7 +150,7 @@ def test_unknown_time_units_raise(tmp_path): proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], - forecastlength=1, + forecastlength=timedelta(hours=1), outpath=str(tmp_path), ) @@ -164,7 +164,7 @@ def test_forecast_reference_time_created(tmp_path): proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [init_time], - forecastlength=4, + forecastlength=timedelta(hours=4), outpath=str(tmp_path), ) @@ -187,7 +187,7 @@ def test_forecast_attributes_removed(tmp_path): proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], - forecastlength=4, + forecastlength=timedelta(hours=4), outpath=str(tmp_path), ) @@ -207,7 +207,7 @@ def test_cube_skipped_if_insufficient_data(tmp_path): proc_reanalysis.create_forecasts( iris.cube.CubeList([cube]), [datetime(2024, 1, 1)], - forecastlength=10, + forecastlength=timedelta(hours=10), outpath=str(tmp_path), ) @@ -220,7 +220,7 @@ def test_multiple_cubes_processed(tmp_path): proc_reanalysis.create_forecasts( iris.cube.CubeList([cube1, cube2]), [datetime(2024, 1, 1)], - forecastlength=4, + forecastlength=timedelta(hours=4), outpath=str(tmp_path), ) From 2985aa86515a34252f528f9fc59f4cd866c62440 Mon Sep 17 00:00:00 2001 From: James Warner Date: Wed, 19 Aug 2026 15:02:43 +0100 Subject: [PATCH 56/56] update type hints --- utils/proc_reanalysis/process_reanalysis.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/utils/proc_reanalysis/process_reanalysis.py b/utils/proc_reanalysis/process_reanalysis.py index 40e9cf1c76..19d8566b99 100755 --- a/utils/proc_reanalysis/process_reanalysis.py +++ b/utils/proc_reanalysis/process_reanalysis.py @@ -54,7 +54,7 @@ def identify_number_of_cycles_required( def create_forecasts( reanalysis: iris.cube.CubeList, forecast_initialisations: list, - forecastlength: datetime.timedelta, + forecastlength: timedelta, outpath: str, ) -> None: """Create forecast files from reanalysis data. @@ -72,7 +72,7 @@ def create_forecasts( will be extracted. forecast_initialisations: list Forecast initialisation times to process. - forecastlength: datetime.timedelta + forecastlength: timedelta Forecast length in hours. outpath: str Directory to which the generated forecast files will be saved.