Add support for rainfall amount diagnostics

[ukmo-huw-lewis]

Addresses #1367 and #1479.

a) Adds colorbar settings for rainfall_amount diagnostics (would capture rainfall and snowfall amount).
b) Introduces a new histogram method for rainfall amount, based on https://gmd.copernicus.org/articles/10/57/2017/gmd-10-57-2017.html (and subsequent iteration by Segolene Berthou in published work, used in K-Scale etc).
c) Moved unit conversion to read.py callback rather than embedded in plotting code.

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.
• Conda lock files have been updated if dependencies have changed.
• Attributed any Generative AI, such as GitHub Copilot, used in this PR.
• Marked the PR as ready to review.

[github-actions]

[ukmo-huw-lewis]

Link to full cset-workflow outputs (based on 2 West Africa UM vs LFRic cases):
https://wwwspice/~huw.lewis/CSET/CSET_RAL3LFRic/WestAfrica_4p4_708_LOCALDATA_1479/

Plots below illustrate alternative precip pdf methods proposed for amount (left) and rate (right). There is nothing particularly "unique" about one needing to be in one form, the other another, but I suggest there may be value for building confidence and insight on the alternative methods to generate both (i.e. exploit availability of both amount and rate diagnostics) to compare over time etc:

Examples below highlight near-but-not-exact equivalence between accumulated rainfall amount in hour (LHS) and the instantaneous rainfall rate on hour (RHS).

Marking as ready for review, propose by @Sylviabohnenstengel given has looked after initial precip pdf treatment.

[Sylviabohnenstengel]

Could you please try to run a timeseries? I am getting the following error message for this branch:
ValueError: Constraint doesn't produce single cube: ConstraintCombination(ConstraintCombination(ConstraintCombination(Constraint(), Constraint(name='microphysical_rainfall_rate'), ), Constraint(cube_func=<function generate_cell_methods_constraint..check_no_aggregation at 0x7f0a9e6b37e0>), ), Constraint(cube_func=<function generate_level_constraint..no_levels at 0x7f0a9e6b3600>), )
< No cubes >

[Sylviabohnenstengel]

I can see that you managed to run the timeseries on your branch. I think the units for the timeseries could also benefit from conversion into mmhr-1. Currently they give kg m-2 s-1 while other plots give kg m-2 hr-1 and we probably want mm hr-1 as unit instead. Given that the conversion to mmhr-1 has moved back into the read function is this a problem with the plotting function or the conversion? Have the data been mutiplied with 3600. to get kg m-2 hr-1 or mm hr-1?

[Sylviabohnenstengel]

Changes look correct and happy to approve when comments in PR are addressed.

[Sylviabohnenstengel]

I can successfully run the branch now

[Sylviabohnenstengel]

Spatial plots
For the spatial plot do we want to keep units of kg m-2 or change that to mm as that would be more intuitive?

Histograms
Regarding the alternative methods for PDFs for amount and rate, respectively, I am happy to go ahead with it for now and see which style is more useful. For the test cases I have run the y_range jumps significantly between frames at the moment. By allowing to tailor the y_range for each frame we can do a better deep dive to see differences between models, however it will be more difficult to see the time evolution, which is why we fixed the y_range in the first place.
I would be useful to pull through a bit more information on the different pdf methods into the website i.e. how to interpret Klingaman et al.. And actually the current text displayed next to the pdf for amounts is probably incorrect for Klingaman et al.

The histogram for microphysical_rainfall_amount is now in units of kg m-2 instead of mm, while the spatial plots are in mm hr-1.
It would be more intuitive to have mm instead of kg m-2 on the axis.

As you can see in the figure above the y_range has a tiny 1e-6 at the top. This could lead to users reading the wrong amount if they only look at the ticks.

Aggregation
We are topping out the y_range for the aggregation for both rainfall_rate and rainfall_amount

[ukmo-huw-lewis]

Review responses

Have pushed updates for:
i) merge of latest trunk into branch,
ii) add unit conversion of amount to translate "kg m-2" to "mm",
iii) updated histogram recipe description text.

Re-running a cset-workflow example and will post link to output here.

On other points of feedback

1) Timeseries plotting

Below illustrates a rainfall rate timeseries (generated using cset bake on command line) that has correct mm/hr units. Colorbar file has been updated to set y-axis max to 1.0 mm/hr. Not sure of difference in setup with example posted by Sylvia further above in this conversation. As aside here, we could/should push small PR to update timeseries axis labelling (currently based on cube.name() as can often pick up different variable name to title, and the use of 'units' for time axis label pulls out the new basing of all time dimensions to hours since...., which is not reflecting the dates shown

2) Spatial plots

Have put in a new layer to unit conversion to translate kg m-2 to mm (just a 1:1 conversion), and added relevant testing. I also updated the logging.WARNING to logging.INFO here (and updated tests), partly as lots of warning messages result for reading in probability_of_visibility cubes. Below shows (left) rainfall amount, now in mm and (right) rainfall rate. _note different cube.name() based colorbar title to the varname-based title. I think this suggestion therefore addressed.

3) Histograms

Note #1377 will offer some tidying of histogram default behaviour (i.e. 1-histogram per case, or select to output hourly histograms if required). Note the 'jumping' of y-range for histograms is common to ALL variables (and so more noticeable when scanning through hour-by-hour). When developing #1396, I left the y-axis of histograms alone. Given the area under curve illustrates the mean rainfall when using 'Klingaman' method, there may be value in setting a ax.set_ylim(), but I just note for now that this has not been standard previously.

Note rainfall amount histograms now also in units "mm".

Useful to also highlight #1303 as in-progress PR to add greater range of methods for histogram plotting. I have built this PR on top of the update for rainfall_rate, but should resolve choices on ymin, ymax and flexibility of methods in considering how to merge this PR. Note there is still likely to be a need to control different histogram choices for precip variables than for others however.

Histogram recipe descriptions have been updated as suggested (bearing in mind further revision proposed via #1303). Note also that the description suggested for main histogram recipe in #1303 discusses noisy outputs.

Wondering if 'jumpiness' might be a separate issue/PR rather than resolving here, and re-review after #1377 if histograms more likely used for more aggregated inputs?

4) Aggregation

colorbar.json file updated to set ymax=1.0 (some trial and error may be needed to find appropriate value, but proposing not to set to "auto" as default to reduce jumpiness. It is now easier to set ymax values in knowledge that rainfall variables will ALWAYS have unit conversion (i.e. in callback) whereas previously spatial plots had unit conversion, but not sure if timeseries did??

Example below shows (left) rainfall amount in mm and (right) rainfall rate time series aggregation by hour of day.

Suggest the aggregation challenge is now resolved?

@Sylviabohnenstengel - could you please check responses above and latest code changes?

[ukmo-huw-lewis]

Output directory from testing with SWUK 300m data with this branch: https://wwwspice/~huw.lewis/CSET/CSET_testing/

[Sylviabohnenstengel]

happy with the changes and ready to merge.