Which scipy.stats functions can be used with xarray.reduce()?

What is your issue?

I would like to use scipy.stats.skew and scipy.stats.kurt with multidimensional xarrays, and also with results of the "coarsen" function.

I don't see information about which functions can be used in the documentation, and I find that the errors that I get are difficult to understand (for me).

I'm not sure it is a bug, though, so I report it as an issue. Here is a test code:

import xarray as xr
import scipy.stats
import numpy as np
#  define random dataarrays, 1D and 2D -- 
da = xr.DataArray( np.random.rand(100),dims=['x'])
da2 = xr.DataArray( np.random.rand(10,10),dims=['x','y'])
#  Work on the whole one dimensional array, OK -- 
print('skewness 1D:',da.reduce(func=scipy.stats.skew).values)
#
#  with 2D array it does not work -- 
print('skewness 2D:',da2.reduce(func=scipy.stats.skew).values)
# This gives error: ValueError: dimensions () must have the same length as the number of data dimensions, ndim=1 --
#
#  try to compute the skewness in bins of 10 using coarsen: does not work --
daskew=da.coarsen(x=10).reduce(func=scipy.stats.skew)
#  gives error: TypeError: tuple indices must be integers or slices, not tuple -- 
#
daskew=da.coarsen(x=10).reduce(func=scipy.stats.skew,axis=None)
#  gives error:  TypeError: skew() got multiple values for keyword argument 'axis'

Start with the supplied Python examples and read the xarray.reduce and coarsen entry points alongside scipy.stats.skew and kurt. Determine which call signatures are supported and what users should expect from multidimensional and coarsened reductions; done means the documentation clearly states supported functions and explains the reported errors.