Q:

Xarray get hourly means for each month in a dataset (without manual looping?)

In xarray it is easy to get the monthly means:
ds.groupby('time.month').mean(dim='time')
and the hourly means:
ds.groupby('time.hour').mean(dim='time')
but I can't find an efficient way to get the hourly means for each month... you can create your own index combining month and hour but that is extremely slow (about 12 times longer than the grouping by hours alone due to dask...)
I want something like this:
<xarray.Dataset>
Dimensions:    (hour: 24, latitude: 721, longitude: 1440, month: 12)
Coordinates:
  * longitude  (longitude) float32 0.0 0.25 0.5 0.75 ... 359.25 359.5    359.75
  * latitude   (latitude) float32 90.0 89.75 89.5 89.25 ... -89.5 -89.75 -90.0
  * month      (month) int64 1 2 3 4 5 6 7 8 9 10 11 12
  * hour       (hour) int64 0 1 2 3 4 5 6 7 8 9 ... 15 16 17 18 19 20 21 22 23

Data variables:
    value       (hour, month, latitude, longitude) float32 dask.array<shape=(24, 12, 721, 1440), chunksize=(1, 1, 721, 1440)>

is this possible?
EDIT: Perhaps setting values to the first day of each month like 2014-01-01 01:00, 2014-02-01 02:00 etc. would be possible?

A:

You can assign a stacked MultiIndex as a new coordinate, then group by this coordinate:
da.coords['monthhour'] = (
    ('time', ),
    pd.MultiIndex.from_arrays([da.time.dt.month, da.time.dt.hour])

da.groupby('monthhour').mean(dim='time')

the resulting coordinate can then be unstacked (da.unstack('monthour')). I haven't found a way to preserve the dim names through a groupby(MultiIndex) operation, but it does most of what you're looking for.

