aggregate_downsample has unexpected behaviour for multidimensional variables #13225

Tested against the development version.

Description

When using aggregate_resample to downsample a TimeSeries, it reduces any multidimensional columns to 1-d.

Expected behavior

I expected aggregate_resample to downsample the time axis, but not any other axes of the data.

Actual behavior

The aggregate_resample function also downsampled the column of multidimensional data on the other axes, which in my case was frequency.

Steps to Reproduce

1. Create a TimeSeries where one column is a multidimensional array.
2. Aggregate downsample that TimeSeries to reduce the time resolution.
3. The aggregate operation is applied across the whole multidimensional array, reducing it to 1-d.

import numpy as np
from astropy import units as u
from astropy.timeseries import TimeSeries, aggregate_downsample

t = TimeSeries(time_start="2022-01-01", time_delta=u.s*1, n_samples=30)
t['data'] = np.ones([30, 3])
print(t)

d = aggregate_downsample(t, time_bin_size=3*u.s)
print(d)

I think that the default behaviour for aggregate_downsample should be to call np.nanmean with axis=1.

System Details

Linux-5.13.0-40-generic-x86_64-with-glibc2.29
Python 3.8.10 (default, Mar 15 2022, 12:22:08)
[GCC 9.4.0]
Numpy 1.20.2
pyerfa 2.0.0.1
astropy 5.0.1
Scipy 1.6.3
Matplotlib 3.4.2