[Bug]: Fixed order and number of dimensions assumed by spatial averager API

[tomvothecoder]

What happened?

Related to

• [Bug]: Coordinate ValueError using .group_average #748
• Fix incorrect dimension used for temporal weights generation #749 (comment)

The spatial averaging API assumes a fixed order and number of dimensions while handling bounds. For example, for lat_bnds, it assumes (lat, bnds).

This assumption can possibly presents issues if the order of dims are different or there are more dimensions. We assume fixed order of dimensions using positional indexing in spatial.py here:

xcdat/xcdat/spatial.py

Line 424 in c198620

pm_cells = np.where(domain_bounds[:, 1] - domain_bounds[:, 0] < 0)[0]

xcdat/xcdat/spatial.py

Line 496 in c198620

return np.abs(domain_bounds[:, 1] - domain_bounds[:, 0])

xcdat/xcdat/spatial.py

Lines 587 to 601 in c198620

	if r_bounds[1] >= r_bounds[0]:
	# Case 1 (simple case): not wrapping around prime meridian (or
	# latitude axis).
	# Adjustments for above / right of region.
	d_bounds[d_bounds[:, 0] > r_bounds[1], 0] = r_bounds[1]
	d_bounds[d_bounds[:, 1] > r_bounds[1], 1] = r_bounds[1]
	# Adjustments for below / left of region.
	d_bounds[d_bounds[:, 0] < r_bounds[0], 0] = r_bounds[0]
	d_bounds[d_bounds[:, 1] < r_bounds[0], 1] = r_bounds[0]
	else:
	# Case 2: wrapping around prime meridian [for longitude only]
	domain_lowers = d_bounds[:, 0]
	domain_uppers = d_bounds[:, 1]
	region_lower, region_upper = r_bounds

What did you expect to happen? Are there are possible answers you came across?

Dynamically handle bounds and their dimensions successfully

Minimal Complete Verifiable Example (MVCE)

Relevant log output

Anything else we need to know?

No response

Environment

Latest main and xCDAT <=0.8.0

[pochedls]

I'm not sure how straightforward this issue will be to address. An alternate strategy would be to use xESMF for spatial averages (and abandon the existing algorithm).

[tomvothecoder]

I'm not sure how straightforward this issue will be to address. An alternate strategy would be to use xESMF for spatial averages (and abandon the existing algorithm).

I think the solution would be similar to #749, where we replace numpy index-based arithmetic with xarray label-based arithmetic. We get the label (dimension names) dynamically retrieved using a new function called get_bounds_dim(), which we pass to the required operations.

[tomvothecoder]

This issue is probably not too urgent because nobody reported it before me. I don't know how common it is for lat and lon bounds to have additional dimensions beyond something like (lat, bnds) (lon,bnds), but it seems rare.