import xarray as xr
import rioxarray as rxr
import numpy as np
from rasterio.enums import Resampling

data = [[[1, 2, 3, 4], [5, np.nan, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]],
        [[17, 18, 19, 20], [21, 22, 23, 24], [25, np.nan, 27, 28], [29, 30, 31, np.nan]]]
x = np.arange(100, 116, 4)
y = np.arange(36, 20, -4)
time = [1, 2]

foo = xr.DataArray(data, coords=[time, y, x], dims=['time', 'y', 'x'])
foo = foo.rio.write_crs('EPSG:4326')
foo = foo.rio.write_nodata(np.nan)

foo_resampled = foo.rio.reproject(
    foo.rio.crs,
    resolution=8,
    resampling=Resampling.average
)

print(f'foo: \n {foo}')
print(f'foo_resampled: \n {foo_resampled}')

 <xarray.DataArray (time: 2, y: 4, x: 4)> Size: 256B
array([[[ 1.,  2.,  3.,  4.],
        [ 5., nan,  7.,  8.],
        [ 9., 10., 11., 12.],
        [13., 14., 15., 16.]],

       [[17., 18., 19., 20.],
        [21., 22., 23., 24.],
        [25., nan, 27., 28.],
        [29., 30., 31., nan]]])
Coordinates:
  * time         (time) int64 16B 1 2
  * y            (y) int64 32B 36 32 28 24
  * x            (x) int64 32B 100 104 108 112
    spatial_ref  int64 8B 0
Attributes:
    _FillValue:  nan

 <xarray.DataArray (time: 2, y: 2, x: 2)> Size: 64B
array([[[ nan,  5.5],
        [11.5, 13.5]],

       [[19.5, 21.5],
        [ nan,  nan]]])
Coordinates:
  * x            (x) float64 16B 102.0 110.0
  * y            (y) float64 16B 34.0 26.0
  * time         (time) int64 16B 1 2
    spatial_ref  int64 8B 0
Attributes:
    _FillValue:  nan

<xarray.DataArray (time: 2, y: 2, x: 2)> Size: 64B
array([[[ 2.66666667,  5.5       ],
        [11.5       , 13.5       ]],

       [[19.5       , 21.5       ],
        [28.        , 28.66666667]]])
Coordinates:
  * x            (x) float64 16B 102.0 110.0
  * y            (y) float64 16B 34.0 26.0
  * time         (time) int64 16B 1 2
    spatial_ref  int64 8B 0
Attributes:
    _FillValue:  nan