17.6: Groupby and Resample

Xarray's groupby functionality is similar to that of pandas and it allows one to aggregate over one or multiple dimensions based on a coordinate. In atmospheric sciences, this feature is invaluable, as it allows for the segmentation of complex, high-dimensional datasets into meaningful groups based on specified criteria, such as time periods or spatial regions, and then apply operations to each group independently.

For instance, suppose we have a dataset of monthly global temperature readings, spanning several years. We may wish to compute the monthly average (i.e., the average of all Januaries, the average of all Februaries, ...). In Xarray, we could accomplish this using the `.groupby()` method along with the coordinate we want to group by, which in this case could be 'time.month'. Once grouped, we would then apply an aggregation method `.mean()`.

Here's a short example demonstrating this:

import xarray as xr
import numpy as np
import pandas as pd

# Set the seed for reproducibility
np.random.seed(0)

# Create a range of monthly periods across multiple years (2005-2007, for example)
time = pd.date_range('2005-01-01', '2007-12-31', freq='M')

# Create a DataArray with random temperature data for each month
temperature = xr.DataArray(np.random.rand(len(time), 2, 2),
                           dims=('time', 'lat', 'lon'),
                           coords={'time': time,
                                   'lat': [10, 20],
                                   'lon': [50, 60]})

# Create a DataArray with random precipitation data for each month
precipitation = xr.DataArray(np.random.rand(len(time), 2, 2),
                             dims=('time', 'lat', 'lon'),
                             coords={'time': time,
                                     'lat': [10, 20],
                                     'lon': [50, 60]})

# Combine into a Dataset
ds = xr.Dataset({'temperature': temperature, 'precipitation': precipitation})

# Group by 'time.month' and calculate the mean for each month across all years
monthly_mean = ds.groupby('time.month').mean('time')

print(monthly_mean)

Here is the output:

<xarray.Dataset>
Dimensions:        (lat: 2, lon: 2, month: 12)
Coordinates:
  * lat            (lat) int64 10 20
  * lon            (lon) int64 50 60
  * month          (month) int64 1 2 3 4 5 6 7 8 9 10 11 12
Data variables:
    temperature    (month, lat, lon) float64 0.4836 0.3663 ... 0.3119 0.4013
    precipitation  (month, lat, lon) float64 0.6187 0.5886 ... 0.5526 0.425

Note that what .groupby() has done is take the months from the years and average them, so that the coordinate month = 1, 2, ... is the average of Januarys, average of Februaries, ..., from three years of the dataset.  It has also retained the `lat` and `lon` coordinates.

By changing the function after the .groupby() method, we can get other statistics of the group.  For example, we could get the maximum of each month with the command: monthly_max = ds.groupby('time.month').max('time'). 

A related process is resampling, a powerful tool that enables one to restructure a dataset from one time resolution to another, aggregating the data according to a given rule. It can be particularly useful when looking for long-term trends or smoothing out short-term fluctuations. This process involves defining a new time frequency and then specifying an aggregation method, such as taking the mean, sum, or maximum over each period in the resampled dataset.

If we use the same example Dataset from above, here's how we can take the dataset with monthly data and resample it to get annual averages:

ds.resample(time='A').mean()

It returns:

<xarray.Dataset>
Dimensions:        (lat: 2, lon: 2, time: 3)
Coordinates:
  * lat            (lat) int64 10 20
  * lon            (lon) int64 50 60
  * time           (time) datetime64[ns] 2005-12-31 2006-12-31 2007-12-31
Data variables:
    temperature    (time, lat, lon) float64 0.5198 0.6375 ... 0.6292 0.5636
    precipitation  (time, lat, lon) float64 0.6056 0.5719 ... 0.5868 0.5283

In the above code, `ds.resample(time='A')` changes the dataset's temporal resolution from monthly ('M') to annual ('A'). Following resampling, `.mean()` computes the average of the data points for each year, which could reveal year-over-year changes and trends that are less discernible in monthly variations.

Note that Xarray resample command is .resample(time='A').  This is different from Pandas, which would just require .resample('A') because Xarray needs to know which coordinate to resample over.