Compute volume from cell_area if available

[enekomartinmartinez]

Description

This PR uses _try_adding_calculated_cell_area in calculate_volume, so the volume can be computed from cell_area if available in the cube. This opens the possibility of working with volume_statistics in irregular grids with areacello but not volcello.

Closes #2310

Link to documentation:

---

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
• 🛠 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

Note

I have tried to run tests locally so I can develop a new test but I am having some issues. After installing the environment.yml with conda (conda env create -f environment.yml) and ESMValCore in developer mode with pip (pip install -e .), I am having the following error message when trying to run the tests:

________________________________________________________________ ERROR collecting test session _________________________________________________________________
../../.miniconda/envs/esmvaltool/lib/python3.11/site-packages/pluggy/_hooks.py:501: in __call__
    return self._hookexec(self.name, self._hookimpls.copy(), kwargs, firstresult)
../../.miniconda/envs/esmvaltool/lib/python3.11/site-packages/pluggy/_manager.py:119: in _hookexec
    return self._inner_hookexec(hook_name, methods, kwargs, firstresult)
../../.miniconda/envs/esmvaltool/lib/python3.11/site-packages/_pytest/python.py:210: in pytest_collect_directory
    if pkginit.is_file():
../../.miniconda/envs/esmvaltool/lib/python3.11/pathlib.py:1267: in is_file
    return S_ISREG(self.stat().st_mode)
../../.miniconda/envs/esmvaltool/lib/python3.11/pathlib.py:1013: in stat
    return os.stat(self, follow_symlinks=follow_symlinks)
E   PermissionError: [Errno 13] Permission denied: '/home/lost+found/__init__.py'

---

To help with the number pull requests:

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

@schlunma

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 94.11%. Comparing base (3a0b49f) to head (9230b82).

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2318      +/-   ##
==========================================
+ Coverage   94.10%   94.11%   +0.01%     
==========================================
  Files         246      246              
  Lines       13477    13494      +17     
==========================================
+ Hits        12682    12700      +18     
+ Misses        795      794       -1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[schlunma]

Thanks for the PR!

I am not quite sure what's causing the error here. How are you running the tests? For me, a cd /your/path/to/ESMValCore && pytest usually works.

[enekomartinmartinez]

Thanks @schlunma

Yes, I am doing a cd /path/toESMValCore, and then I have tried several options:

• pytest tests/unit/preprocessor/_volume/test_volume.py (this also produces some errors due to raised warnings so I used the default configuration later)

======================================== ERRORS =========================================
__________________________________ ERROR collecting . ___________________________________
../../.miniconda/envs/esmvaltool/lib/python3.11/site-packages/pluggy/_hooks.py:501: in __call__
    return self._hookexec(self.name, self._hookimpls.copy(), kwargs, firstresult)
../../.miniconda/envs/esmvaltool/lib/python3.11/site-packages/pluggy/_manager.py:119: in _hookexec
    return self._inner_hookexec(hook_name, methods, kwargs, firstresult)
../../.miniconda/envs/esmvaltool/lib/python3.11/site-packages/pytest_mypy.py:111: in pytest_collect_file
    return MypyFile.from_parent(parent=parent, fspath=path)
../../.miniconda/envs/esmvaltool/lib/python3.11/site-packages/pytest_mypy.py:123: in from_parent
    return getattr(super(), 'from_parent', cls)(*args, **kwargs)
../../.miniconda/envs/esmvaltool/lib/python3.11/site-packages/_pytest/nodes.py:661: in from_parent
    return super().from_parent(parent=parent, fspath=fspath, path=path, **kw)
../../.miniconda/envs/esmvaltool/lib/python3.11/site-packages/_pytest/nodes.py:275: in from_parent
    return cls._create(parent=parent, **kw)
../../.miniconda/envs/esmvaltool/lib/python3.11/site-packages/_pytest/nodes.py:152: in _create
    return super().__call__(*k, **kw)
../../.miniconda/envs/esmvaltool/lib/python3.11/site-packages/_pytest/nodes.py:616: in __init__
    path = _imply_path(type(self), path, fspath=fspath)
../../.miniconda/envs/esmvaltool/lib/python3.11/site-packages/_pytest/nodes.py:108: in _imply_path
    warnings.warn(
E   pytest.PytestRemovedIn8Warning: The (fspath: py.path.local) argument to MypyFile is deprecated. Please use the (path: pathlib.Path) argument instead.
E   See https://docs.pytest.org/en/latest/deprecations.html#fspath-argument-for-node-constructors-replaced-with-pathlib-path

This both options gave me the error from above:

• pytest -c /dev/null tests/unit/preprocessor/_volume/test_volume.py

• pytest -c /dev/null /path/toESMValCore/tests/unit/preprocessor/_volume/test_volume.py

I don't know where the problem can come from. I have another conda environment with pytest for another library and it works fine. So, it could be related to the specific configuration of ESMValCore tests; maybe I did something wrong? I will try to reinstall the environment and go deeper if it doesn't work when I have more time. I was asking in case the error was familiar to you.

[schlunma]

Ah, your first error actually looks familiar: #2314

This has been fixed, but requires setting up a new environemt. Sorry, this is really unlucky timing.

[enekomartinmartinez]

Thanks, I will then create the new environment and try to add new tests

[enekomartinmartinez]

Hi @schlunma

I am still working on this, here are some updates with the last commit:

1. Now calculate_volume is more flexible, previous implementation needed to have always z in the second dimensions z_dim == 1. I have made that more generic so this preprocessor can run in models where z is in another position (e.g. a cube after a climate_statistics with no time dimension).
2. It also works with non-regular data if areacello is available

I have tested running volume_statistics over regular data with no areacello and with irregular data with areacello and seems to run fine.

I still have to include the test and correct the ones that failed that were related to the condition of z_dim ==1.

[enekomartinmartinez]

Fixed tests and added new ones, I made some new error messages for the following cases:

• z-axis bounds not 2 in last dim -> cannot compute thickness
• xyz_dims > xy_dims +1 -> z depends on other non-spatial dimensions (like time) this would provoke aggregating the cube also in that dimension (this is an adaptation of a previous error)

Still missing test when cell_measure is available and when have to guess bounds of z. I will add them the following days.

[enekomartinmartinez]

Hi @schlunma,

This is ready to merge by my side

[schlunma]

Thanks for these change @enekomartinmartinez, looks really good already!

I have a couple of comments directly on the code.

In addition, it would be great if you could make sure that areacello is loaded automatically if no supplementary is specified (this will not raise an error if that is not found). You can do that by adding it to the decorator here. Also, please make sure the documentation is up-to-date (i.e., here and here).

Thanks!!

[enekomartinmartinez]

Thank you @schlunma for so many comments and suggestions, I will be accepting the changes and answering them one by one!

[enekomartinmartinez]

Thanks for these change @enekomartinmartinez, looks really good already!

I have a couple of comments directly on the code.

In addition, it would be great if you could make sure that areacello is loaded automatically if no supplementary is specified (this will not raise an error if that is not found). You can do that by adding it to the decorator here. Also, please make sure the documentation is up-to-date (i.e., here and here).

Thanks!!

I am working in this other suggestions, I still have few questions. Should I register as well arecella in the case volume_statistics is used with atmospheric variables?

Shouldn't we check in calculate_volume that z-coord is in meters, so we ensure always that the returned volume is in cubic meters?

[schlunma]

Should I register as well arecella in the case volume_statistics is used with atmospheric variables?

Good point. The function has clearly been designed for ocean data, but it would of course be nice to also support atmospheric variables. However, to properly support this we would probably need to consider pressure levels in addition to altitudes. Let's do that in a separate PR.

Shouldn't we check in calculate_volume that z-coord is in meters, so we ensure always that the returned volume is in cubic meters?

Also good point! You could also return the multiplied units of the area and z dimension in calculate_volume and then use that as units for the cell measure. Then it would also work for z-coords with different units.

[enekomartinmartinez]

Good point. The function has clearly been designed for ocean data, but it would of course be nice to also support atmospheric variables. However, to properly support this we would probably need to consider pressure levels in addition to altitudes. Let's do that in a separate PR.

I cannot currently work on further development. I don't know if this is a requirement but there may be people already using the preprocessor with atmospheric data as it may not fail and use the thickness of pressure_levels automatically. Should this be considered risky?

Also good point! You could also return the multiplied units of the area and z dimension in calculate_volume and then use that as units for the cell measure. Then it would also work for z-coords with different units.

My question is should we allow using things like area * thickness in sigma space as a correct volume measure? Wouldn't the weights (kg/m) be incorrect and produce an unrealistic product?

[schlunma]

I cannot currently work on further development. I don't know if this is a requirement but there may be people already using the preprocessor with atmospheric data as it may not fail and use the thickness of pressure_levels automatically. Should this be considered risky?

I am not saying you need to do this 😉 If you think adding areacella is useful, I am happy with it. In this case it might be a good idea to add the ['volcello', 'areacello', 'areacella'] combination to this function, so that for ocean variables we get ['volcello', 'areacello'], and for others ['areacella']. It would probably also be a good idea to a add a note to the preprocessor that states that this preprocessor has been designed for oceanic variables, but it might be applicable to atmospheric data as well.

My question is should we allow using things like area * thickness in sigma space as a correct volume measure? Wouldn't the weights (kg/m) be incorrect and produce an unrealistic product?

Ah, I didn't consider sigma space. If the units can be converted to m then it shouldn't be a problem (the iris function will automatically adapt the target units if necessary), but if that's not the case this is problematic. So a test if the units are convertible to m would be nice indeed!

[enekomartinmartinez]

Hi @schlunma

I have added a conversion of depth units to meters in calculate_volume and a test for a cube in sigma space with no volcello given, where the conversion fails. I have documented it in the same function, volume_statistics and in doc/recipe/preprocessor.rst. I have also added a note in volume_statistics for the atmospheric case.

About adding areacella, as you prefer, maybe we can leave it without it. Then, in case anyone wants to include a conversion for pressure levels to meters, include it.

[schlunma]

Thanks @enekomartinmartinez for these changes! I have one remaining comment about not changing the input cube in-place, then I will approve. Thanks for your patience!!

[enekomartinmartinez]

Hi @schlunma

Thanks for the suggestions. I have accepted them and removed the change of units as when doing the copy is unnecessary.

[schlunma]

Great! Thanks for your contribution 🚀