Support analysis/reanalysis comparisons

[James Warner (jwarner8)]

Addresses #2016

Contribution checklist

Aim to have all relevant checks ticked off before merging. See the developer's guide for more detail.

• Documentation has been updated to reflect change.
• New code has tests, and affected old tests have been updated.
• All tests and CI checks pass.
• Ensured the pull request title is descriptive.
• Attributed any Generative AI, such as GitHub Copilot, used in this PR.
• Marked the PR as ready to review.

[github-actions]

Total coverage: 89% (HTML report)

Name                                                              Stmts   Miss Branch BrPart  Cover
---------------------------------------------------------------------------------------------------
src/CSET/__init__.py                                                 93      2     12      0    98%
src/CSET/_common.py                                                 149      0     52      0   100%
src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py             115     27     24      0    79%
src/CSET/cset_workflow/app/finish_website/bin/finish_website.py      70      0      4      0   100%
src/CSET/cset_workflow/app/parbake_recipes/bin/parbake.py            29      0      8      0   100%
src/CSET/cset_workflow/app/send_email/bin/send_email.py              25      0      4      0   100%
src/CSET/cset_workflow/lib/python/jinja_utils.py                     17      0      6      0   100%
src/CSET/extract_workflow.py                                         47      0     16      0   100%
src/CSET/graph.py                                                    43      0     14      0   100%
src/CSET/operators/__init__.py                                       89      0     26      0   100%
src/CSET/operators/_atmospheric_constants.py                          9      0      0      0   100%
src/CSET/operators/_colormaps.py                                    229      4     62      4    97%
src/CSET/operators/_stash_to_lfric.py                                 3      0      0      0   100%
src/CSET/operators/_utils.py                                        183      8     74      6    95%
src/CSET/operators/ageofair.py                                      141      7     64      5    94%
src/CSET/operators/aggregate.py                                      76      1     22      1    98%
src/CSET/operators/aviation.py                                       60      0     18      0   100%
src/CSET/operators/collapse.py                                      212     72     88      7    65%
src/CSET/operators/constraints.py                                   111      7     48      2    93%
src/CSET/operators/convection.py                                     37      4     10      2    87%
src/CSET/operators/ensembles.py                                      27      0     14      0   100%
src/CSET/operators/feature.py                                        41      0     10      0   100%
src/CSET/operators/filters.py                                        66      2     30      0    98%
src/CSET/operators/fluxes.py                                         41      0     10      0   100%
src/CSET/operators/humidity.py                                      139      0     56      0   100%
src/CSET/operators/imageprocessing.py                                56      0     16      0   100%
src/CSET/operators/mesoscale.py                                      17      0      2      0   100%
src/CSET/operators/misc.py                                          158      1     66      4    98%
src/CSET/operators/plot.py                                          929    219    318     61    73%
src/CSET/operators/power_spectrum.py                                 97      3     30      3    95%
src/CSET/operators/precipitation.py                                 136      0     72      0   100%
src/CSET/operators/pressure.py                                       41      0     12      0   100%
src/CSET/operators/read.py                                          409     38    178     14    89%
src/CSET/operators/regrid.py                                        122     18     64      1    83%
src/CSET/operators/scoreswrappers.py                                 47      6     12      3    85%
src/CSET/operators/temperature.py                                   121      0     32      0   100%
src/CSET/operators/transect.py                                       62      0     24      0   100%
src/CSET/operators/wind.py                                           45      3     10      2    91%
src/CSET/operators/write.py                                          15      0      6      0   100%
src/CSET/recipes/__init__.py                                        101      0     28      0   100%
---------------------------------------------------------------------------------------------------
TOTAL                                                              4408    422   1542    115    89%

[James Warner (jwarner8)]

I now have a function in read.py that detects if the cube forecast_period is static (T+0h) and forecast_reference_time is changing, and thus it is likely to be analysis/reanalysis. If so, then it will modify the cube to recreate a varying auxcoord forecast_period and set the forecast_reference_time to be static. Valid time is unchanged throughout. This will allow us to compare models as a function of lead-time with reanalysis.

I.e. before

    Dimension coordinates:
        time                             x             -                -
        latitude                         -             x                -
        longitude                        -             -                x
    Auxiliary coordinates:
        forecast_reference_time          x             -                -
        grid_latitude                    -             x                -
        grid_longitude                   -             -                x
    Scalar coordinates:
        forecast_period             0.0 hours
        height                      1.5 m
        realization                 0
    Attributes:
        Conventions                 'CF-1.7'
        um_version                  '11.9'
DimCoord :  time / (hours since 1970-01-01 00:00:00, standard calendar)
    points: [
        2024-01-19 00:00:00, 2024-01-19 06:00:00, ...,
        2024-01-26 18:00:00, 2024-01-27 00:00:00]
    shape: (33,)
    dtype: float64
    standard_name: 'time'
    var_name: 'time'
DimCoord :  forecast_reference_time / (seconds since 1970-01-01 00:00:00, standard calendar)
    points: [
        2024-01-19 00:00:00, 2024-01-19 06:00:00, ...,
        2024-01-26 18:00:00, 2024-01-27 00:00:00]
    shape: (33,)
    dtype: int64
    standard_name: 'forecast_reference_time'
    var_name: 'forecast_reference_time'
DimCoord :  forecast_period / (hours)
    points: [0.]
    shape: (1,)
    dtype: float64
    standard_name: 'forecast_period'
    var_name: 'forecast_period'

after

air_temperature / (K)               (time: 33; latitude: 1920; longitude: 2560)
    Dimension coordinates:
        time                             x             -                -
        latitude                         -             x                -
        longitude                        -             -                x
    Auxiliary coordinates:
        forecast_period                  x             -                -
        grid_latitude                    -             x                -
        grid_longitude                   -             -                x
    Scalar coordinates:
        forecast_reference_time     2024-01-19 00:00:00
        height                      1.5 m
        realization                 0
    Attributes:
        Conventions                 'CF-1.7'
        um_version                  '11.9'
DimCoord :  time / (hours since 1970-01-01 00:00:00, standard calendar)
    points: [
        2024-01-19 00:00:00, 2024-01-19 06:00:00, ...,
        2024-01-26 18:00:00, 2024-01-27 00:00:00]
    shape: (33,)
    dtype: float64
    standard_name: 'time'
    var_name: 'time'
AuxCoord :  forecast_reference_time / (hours since 1970-01-01 00:00:00, standard calendar)
    points: [2024-01-19 00:00:00]
    shape: (1,)
    dtype: float64
    standard_name: 'forecast_reference_time'
AuxCoord :  forecast_period / (hours)
    points: [  0.,   6., ..., 186., 192.]
    shape: (33,)
    dtype: float64
    standard_name: 'forecast_period'

[James Warner (jwarner8)]

[James Warner (jwarner8)]

Challenges:

• How to identify what the acceptable forecast_reference_time should be? As we do not have access to cycle time in read.py (?). If we know the other models being compared (their forecast_reference_time and forecast lead-time), we can make sure we pick up the appropriate reanalysis.

I could back out the path/cycle time, but that would not work if we are baking a recipe on the command line, so avoiding this. I am going to investigate fixing the reanalysis later at a point where I have access to the other models, so I can be a bit more clever with setting the forecast_reference_time to be more appropriate - possibly at the point where we need to collapse and find common time points?

[James Warner (jwarner8)]

Have a function now working for aggregation, creating two analyses cubes that match the reference time and forecast period of the forecast (if two cases are being aggregated)

[James Warner (jwarner8)]

Current thinking. The purpose of this PR is to support the use of reanalysis/analysis as another model across CSET, for all functionality where applicable. I started off creating a new section in the GUI where a user would toggle reanalysis on/off, but I've started backing away from that idea. This is because of the extra overhead of duplicating a huge set of recipes to have recipes that include analysis as a comparison. The current solution I'm leaning towards is that the user simply treats reanalysis as another model. This would also allow us to use multiple reanalyses (UM analysis, alongside ERA5, for example). The main two pieces of work I can think of is

1. Adapting plot.py to search for 'analyses' in the model name and do something special with analysis (i.e. force to be black, first line plotted, etc.)
2. A function that examines what forecast cubes we have, and construct a similar cube from the reanalysis. I have made quite a bit of progress with this. The idea currently (in collapse.py at present), is that after we've read the cubes and done the appropriate metadata cleaning, we will be handed a cubelist containing each model. This could be an forecast, alongside some reanalysis. If there is only one forecast (i.e. no aggregation happening), the cube will be structured valid time, lat, lon (with forecast reference time and forecast period as scalar and auxcoords respectively). If there is multiple forecasts, the cube will be structured forecast_period, forecast_reference_time, lat, lon, where valid time is an auxcoord. We want to be able to plot models as a function of lead-time, alongside reanalysis (which doesn't have a forecast_period, or rather forecast_period is 0). The current approach in the function is as follows:

• a. The cubelist is separated into forecast cubes and analysis cubes, by determining if their maximum forecast period is zero (if so, likely to be analysis).
• b. These should both contain one cube.
• c. If we are aggregating over multiple times, then we will get multiple analysis that look the same, and CSET assumes its an ensemble and gives it a realization dimension, so we just extract the index of this as its not an ensemble.
• d. We check the analysis covers the forecast duration.
• e. For each forecast_reference_time (if there are multiple), we work out the valid times, extract this from the reanalysis, and construct a forecast_period and forecast_reference_time, and append this to a fixed_analysis_cubes list.
• Once we've done this, if it is length 1, then there is only one forecast_reference_time and we can add time as an aux coord. If it is length 2, it means it cant merge as there are multiple forecast_reference_times, and so we need to create time as a 2D aux coordinate that maps to forecast_period and forecast_reference_time.
• We then append this fixed analysis cube with the forecast cubes and return.

There are a couple of options here. I did look for a place to put this in read.py, once everything is read in, but not sure how this would work with load_model and generators (as it needs all model cubes in one place).

OptionA: Embed it somewhere in read.py or similar to its accessible to all subsequent recipes/operators, whether its timeseries, power spectra, transects, profiles etc.
OptionB: Make it a separate function in misc.py and get every recipe to add it as a step such as ensure_aggregatable.
OptionC: Pass more arguments to read.py so it can restructure the analysis cubes without a forecast cube to identify forecast_reference_time and forecast_period from, but this would mean updating all recipes with more arguments, along with making cset bake more complicated.
OptionC: Something else??

Any opinions welcome. I am trying to reduce the additional duplication of having analysis as a separate section by including it as a model.

[ukmo-huw-lewis]

Any opinions welcome. I am trying to reduce the additional duplication of having analysis as a separate section by including it as a model.

Have not digested detail of commentary here, but for feedback (as both 'user' - what's simplest, and 'developer' - what's simplest, and what's possible) would strongly support the 'treat' analysis/reanalysis as model concept. It may be that some of the nuance of working with (re)analysis data need to be handled somewhere, but re-using existing structures and recipes etc will greatly reduce duplication as outlined above.

Suggesting there are some parallels here with experience working with point observations (even including the plot order/style comment), where balance of some obs-specific processing needed, but once 'everything looks like a cube', then real power in no longer needing to track variety of output characteristics.

Think the implied flexibility will be highly useful (and expected) for users. Happy to support / feedback where useful to you James on this one.