Use smarter (units-aware) weights

[schlunma]

Description

Since version 3.5.0, iris allows the usage of cell measures (and a lot of other objects) as weights. This has the great benefit that units in statistical calculations are handled properly, e.g., for an area-weighted sum, resulting units are multiplied by m2.

In addition, this PR adds type hints for all functions of the affected modules.

Closes #1613
Closes #708

Link to documentation:

• https://esmvaltool--2139.org.readthedocs.build/projects/ESMValCore/en/2139/recipe/preprocessor.html#area-statistics
• https://esmvaltool--2139.org.readthedocs.build/projects/ESMValCore/en/2139/recipe/preprocessor.html#climate-statistics
• https://esmvaltool--2139.org.readthedocs.build/projects/ESMValCore/en/2139/recipe/preprocessor.html#axis-statistics
• https://esmvaltool--2139.org.readthedocs.build/projects/ESMValCore/en/2139/recipe/preprocessor.html#depth-integration

Backwards-compatibility

The changes of these PR are NOT backwards-compatible. Since this PR fixes calculations that have been scientifically wrong in previous versions, no deprecation cycle is used. The following preprocessors are affected:

Preprocessor	Affected option	Old behavior	New behavior	Affected recipe(s)
area_statistics	operator='sum'	Cube units remained, e.g., kg m-2 s-1	Cube units are multiplied by m2, e.g., kg s-1	recipe_wenzel16nat, recipe_wenzel14jgr
climate_statistics	period='full' and operator='sum'	Cube units remained, e.g., kg m-2 s-1	Cube units are multiplied by time units, e.g., kg m-2	-
axis_statistics	operator='sum'	Cube units remained, e.g., kg m-2 s-1	Cube units are multiplied by coordinate units, e.g., kg m-1 s-1	-
depth_integration	all	Cube units remained, e.g., kg m-2 s-1	Cube units are multiplied by coordinate units, e.g., kg m-1 s-1	recipe_ocean_example

In addition, there has been a bug in climate_statistics, which lead to a wrong calculation of weighted sums if the weights are equal for all time points. These are the relevant lines of code:

ESMValCore/esmvalcore/preprocessor/_time.py

Lines 740 to 746 in 1abaf8e

	if time_weights.min() == time_weights.max():
	# No weighting needed.
	clim_cube = cube.collapsed('time', operator_method)
	else:
	clim_cube = cube.collapsed('time',
	operator_method,
	weights=time_weights)

While this simplification is correct for calculating means and RMSE, the weighted sum actually differs in both cases. Thus, we need to always consider the weights here, regardless if they are equal or not. This PR fixes that.

I will test all affected recipes and open corresponding PR in ESMValTool if any changes are necessary.

EDIT: I tested all affected recipes; they run fine with minimal changes (ESMValGroup/ESMValTool#3300).

---

Before you get started

• ☝ Create an issue to discuss what you are going to do

Checklist

It is the responsibility of the author to make sure the pull request is ready to review. The icons indicate whether the item will be subject to the 🛠 Technical or 🧪 Scientific review.

• 🧪 The new functionality is relevant and scientifically sound
• 🛠 This pull request has a descriptive title and labels
• 🛠 Code is written according to the code quality guidelines
• 🧪 and 🛠 Documentation is available
• 🛠 Unit tests have been added
• 🛠 Changes are backward compatible (see note above)
• 🛠 Any changed dependencies have been added or removed correctly
• 🛠 The list of authors is up to date
• 🛠 All checks below this pull request were successful (I think assert in production code is fine to make other developers aware of something; this is not intended to catch errors by users)

---

To help with the number pull requests:

• 🙏 We kindly ask you to review two other open pull requests in this repository

[codecov]

Codecov Report

Merging #2139 (9dbc2e5) into main (579bb61) will increase coverage by 0.03%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##             main    #2139      +/-   ##
==========================================
+ Coverage   93.10%   93.14%   +0.03%     
==========================================
  Files         237      237              
  Lines       12826    12828       +2     
==========================================
+ Hits        11942    11948       +6     
+ Misses        884      880       -4     

Files Changed	Coverage Δ
esmvalcore/preprocessor/_area.py	94.54% <100.00%> (+0.88%)	⬆️
esmvalcore/preprocessor/_shared.py	100.00% <100.00%> (ø)
esmvalcore/preprocessor/_time.py	97.72% <100.00%> (-0.01%)	⬇️
esmvalcore/preprocessor/_volume.py	91.59% <100.00%> (+1.16%)	⬆️

[valeriupredoi]

this is a very nice PR and introduced functionality, very cool, @schlunma 🍺 Proper testing too 🎖️ (I don't see any reason) but do we expect any changes in results for recipes that use this set of preprocessors? Also, were there any distinct changes in performance you might have noticed? Cheers!

[valeriupredoi]

PS: I meant actual result changes, not if the recipes just run. I guess we could find that out at v2.10 release time - which reminds me, maybe we should curate a list of PRs that may lead to changes in results and have that at hand when testing for a release

[schlunma]

Good point!

• recipe_wenzel16nat: no change (needs Allow NCL unit conversion kg s-1 -> GtC y-1 ESMValTool#3300).
• recipe_wenzel14jgr: no change (needs Allow NCL unit conversion kg s-1 -> GtC y-1 ESMValTool#3300).
• recipe_ocean_example: the units of the plots of diag_depthInt_1 change from K K m, but I think these new units are correct (@ledm please correct me if I am wrong here). Does not need any changes in ESMValTool to run.

Old version:

New version:

[valeriupredoi]

awesome you doing testing, Manu! About those K -> m.K units change - it's a bit fishy (excuse the pun) at first glance since the numbers are the same, only units changed? You wanna do more testing before we merge this? 🍺

[schlunma]

No, nothing fishy at all, that's exactly what was to be expected: In this preprocessor (depth_integration), a weighted sum of potential temperature along the depth coordinate is taken; i.e., potential temperature (units K) is multiplied by the depth layer thickness (units m) and then summed up, which results in the units K m. 🐟 This PR just fixes the units handling; the actual calculation was correct before (hence the numbers didn't change).

This has also been expected by the author of depth_integration; see the following comment I removed now:

ESMValCore/esmvalcore/preprocessor/_volume.py

Lines 256 to 258 in 4c2578b

	# result.units = Unit('m') * result.units # This doesn't work:
	# TODO: Change units on cube to reflect 2D concentration (not 3D)
	# Waiting for news from iris community.

[valeriupredoi]

excellent! I'll take my 🐟 back then 😁