Support analysis/reanalysis comparisons

[James Warner (jwarner8)]

Addresses #2016

Microsoft CoPilot was used to suggest some of the unit tests.

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: 93% (HTML report)

Name                                                              Stmts   Miss Branch BrPart  Cover
---------------------------------------------------------------------------------------------------
src/CSET/__init__.py                                                105      0     14      0   100%
src/CSET/_common.py                                                 156      0     54      0   100%
src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py             117     28     26      0    78%
src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py          81      8     28     11    83%
src/CSET/cset_workflow/app/finish_website/bin/finish_website.py      79      2      8      2    95%
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                                        103      1     26      1    98%
src/CSET/graph.py                                                    44      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                                    249      3     72      4    98%
src/CSET/operators/_stash_to_lfric.py                                 3      0      0      0   100%
src/CSET/operators/_utils.py                                        190      8     72      6    95%
src/CSET/operators/ageofair.py                                      142      7     64      5    94%
src/CSET/operators/aggregate.py                                      77      1     22      1    98%
src/CSET/operators/aviation.py                                       61      0     18      0   100%
src/CSET/operators/collapse.py                                      155      8     72      3    93%
src/CSET/operators/constraints.py                                   113      7     50      3    93%
src/CSET/operators/convection.py                                     38      4     10      2    88%
src/CSET/operators/ensembles.py                                      27      0     14      0   100%
src/CSET/operators/feature.py                                        44      0     10      0   100%
src/CSET/operators/filters.py                                        67      2     30      0    98%
src/CSET/operators/fluxes.py                                         41      0     10      0   100%
src/CSET/operators/humidity.py                                      135      0     52      0   100%
src/CSET/operators/imageprocessing.py                                57      0     16      0   100%
src/CSET/operators/mesoscale.py                                      18      0      2      0   100%
src/CSET/operators/misc.py                                          172      1     72      3    98%
src/CSET/operators/plot.py                                         1117    155    412     71    83%
src/CSET/operators/power_spectrum.py                                 98      3     30      3    95%
src/CSET/operators/precipitation.py                                 204      2     92      2    99%
src/CSET/operators/pressure.py                                       41      0     12      0   100%
src/CSET/operators/read.py                                          432     21    184     16    94%
src/CSET/operators/regrid.py                                        147      1     70      3    98%
src/CSET/operators/scoreswrappers.py                                208     16     58      6    91%
src/CSET/operators/temperature.py                                   121      0     32      0   100%
src/CSET/operators/transect.py                                       63      0     24      0   100%
src/CSET/operators/wind.py                                           46      3     10      2    91%
src/CSET/operators/write.py                                          15      0      6      0   100%
src/CSET/recipes/__init__.py                                        104      0     28      0   100%
src/CSET/sample_data/__init__.py                                      0      0      0      0   100%
---------------------------------------------------------------------------------------------------
TOTAL                                                              5039    281   1758    144    93%

[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.

[James Warner (jwarner8)]

Going to refactor this branch so we have standalone code in contrib that a user points at reanalysis, and write out intermediate output that is suitable for CSET ingestion. PR will also include some small changes to recognise if model name is 'ERA5' or 'UM Analysis' and plot these using a different line style/marker across plot.py. Broader progress in #2341
Another motivation to treating analysis is another model is able to use cset bake, whereas GUI implementation restricts us to the workflow only.

[James Warner (jwarner8)]

TODO:

• Add README documentation
• Add comments to code and docstrings, type hints, tidy up code.
• Show both UM analysis and ERA5 working for some examples, using data from this, in workflow.
• Create unit tests (James Frost (@jfrost-mo) should this live in separate file tests_.py or a folder ./tests in this util? Realise not following testing structure of rest of CSET, and this is a bit separate, but good to have tests). Even just having a script to run tests on and check it performs as expected would be good.

[James Warner (jwarner8)]

Comments addressed, tests passed, and have re-run script with ERA5/UM analysis to check its still actually producing sensible output. Ready for final review James Frost (@jfrost-mo)

[James Frost (jfrost-mo)]

Looks good. A couple optional helptext queries, but happy for you to resolve and merge once happy.

[James Frost (jfrost-mo)]

Looks good. The examples in the README should help people get the quoting right. The tests do need fixing before this can be merged; this looks like you just need to replace the integer hours in the tests with timedelta(hours=N).