Anybody run into the need to average over Z axis (e.g., cdutil.averager(tvar, axis="z")), and what did you use?

[tomvothecoder]

Question criteria

• I added a descriptive title here.
• I searched the xCDAT GitHub Discussions to find a similar question and didn't find it.
• I searched the xCDAT documentation and/or Xarray documentation.

Describe your question

I'm running into a case where I need to average over the Z axis. The CDAT code being used is cdutil.averager(tvar, axis="z"), which I need to replace.

@xCDAT/core-developers, has anybody else found the need to average over the Z axis? And if so, what package or custom code are you using to do this?

xCDAT's spatial averager, which is based on cdutil.averager(), only supports rectilinear grids ("X" and "Y"). However, we do mention that the get_weights() method can be extended to support other axes:

xcdat/xcdat/spatial.py

Lines 246 to 252 in fbf1db6

	Notes
	-----
	This method was developed for rectilinear grids only. ``get_weights()``
	recognizes and operate on latitude and longitude, but could be extended
	to work with other standard geophysical dimensions (e.g., time, depth,
	and pressure).
	"""

Are there are any possible answers you came across?

For the short-term, I might be able to use xarray.DataArray.weighted directly. However, I still need to generate the weights xr.DataArray beforehand.

datarray.weighted(weights).mean("height")

Minimal Complete Verifiable Example (MVCE)

No response

Relevant log output

No response

Environment

No response

Anything else we need to know?

No response

[pochedls]

@tomvothecoder – I have not done this recently, but I think the same algorithm would be valid for depth / pressure coordinates (weights should generally scale as the difference between the bound values for each cell). I think it would make sense to extend get_weights to be able to handle vertical axes.

[lee1043]

@tomvothecoder -- I think PMP is now running into the need of average over Z axis as a part of its ENSO metrics migration to xCDAT from CDAT. Without recalling this discussion was already posted, I advised my colleague to try spatial.average with axis=["Z"] which of course did not work.

filename = "/bdd/CMIP6/CMIP/IPSL/IPSL-CM6A-LR/historical/r1i1p1f1/Omon/thetao/gn/latest/thetao_Omon_IPSL-CM6A-LR_historical_r1i1p1f1_gn_185001-194912.nc"
ds = xcdat.open_mfdataset(filename, add_bounds=["T", "X", "Y", "Z"])
ds_ave = ds.spatial.average("thetao", axis=["Z"])

Which returns the expected error:

ValueError: Incorrect `axis` argument value. Supported values include: X, Y.

I think it might worth to tackle this. Maybe additional item for the Dev Day?

This might also requires generation of weighting if not provided.

[pochedls]

I think this is done in this PR, but it needs testing.