[Bug]: Axis order dependence in spatial averaging

[pochedls]

What happened?

I was looking at this issue and noticed that the results were sensitive to the order of the axis specification in the spatial averager (["X", "Y"] versus ["Y", "X"]). Admittedly the differences are very small (0.001% difference), but not completely round-off and might be larger for datasets with missing data.

What did you expect to happen?

I think order shouldn't matter here. I wonder if the issue is with the xarray averager doing operations iteratively (here).

Minimal Complete Verifiable Example

# note that add_bounds/decode_times are True by default (I dropped it)
ds = xc.open_dataset('/p/css03/esgf_publish/CMIP6/ScenarioMIP/IPSL/IPSL-CM6A-LR/ssp585/r2i1p1f1/Amon/tas/gr/v20191121/tas_Amon_IPSL-CM6A-LR_ssp585_r2i1p1f1_gr_201501-210012.nc')
ts_global_yx = ds.spatial.average("tas", axis=['Y', 'X'], lat_bounds=(-90, 90), lon_bounds=(0, 360))["tas"]
ts_global_xy = ds.spatial.average("tas", axis=['X', 'Y'], lat_bounds=(-90, 90), lon_bounds=(0, 360))["tas"]

print(np.allclose(ts_global_yx, ts_global_xy))
print((ts_global_yx - ts_global_xy).values)
print((ts_global_yx - ts_global_xy).values/ts_global_xy.values*100)

False [0.00418091 0.00296021 0.00378418 ... 0.00350952 0.00372314 0.00375366] [0.0014663 0.00103671 0.00132236 ... 0.00119679 0.00127493 0.00128867]

Relevant log output

No response

Anything else we need to know?

I looked at one other dataset and the differences were round-off level (10^-12: /p/user_pub/work/CMIP6/CMIP/E3SM-Project/E3SM-2-0/historical/r2i1p1f1/Amon/tas/gr/v20220831/). It would be useful to swap in a simple _averager routine that simultaneously weights X/Y for testing. It would also be useful to loop over a bunch of data to see how often this happens (particularly with missing data).

Environment

xcdat 0.5.0 xarray: 2023.1.0

[tomvothecoder]

Notes from past meetings:

• Create a minimal xarray example and post on xarray repo
• Ask if this is how it “should” behave