You can run this notebook in a live session or view it on Github.

# Calculating Seasonal Averages from Timeseries of Monthly Means

Author: Joe Hamman

The data used for this example can be found in the xarray-data repository. You may need to change the path to `rasm.nc` below.

Suppose we have a netCDF or `xarray.Dataset` of monthly mean data and we want to calculate the seasonal average. To do this properly, we need to calculate the weighted average considering that each month has a different number of days.

    [1]:

    %matplotlib inline
    import numpy as np
    import pandas as pd
    import xarray as xr
    from netCDF4 import num2date
    import matplotlib.pyplot as plt

## Some calendar information so we can support any netCDF calendar.

    [2]:

    dpm = {'noleap': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
           '365_day': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
           'standard': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
           'gregorian': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
           'proleptic_gregorian': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
           'all_leap': [0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
           '366_day': [0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
           '360_day': [0, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30]}

## A few calendar functions to determine the number of days in each month

If you were just using the standard calendar, it would be easy to use the `calendar.month_range` function.

    [3]:

    def leap_year(year, calendar='standard'):
        """Determine if year is a leap year"""
        leap = False
        if ((calendar in ['standard', 'gregorian',
            'proleptic_gregorian', 'julian']) and
            (year % 4 == 0)):
            leap = True
            if ((calendar == 'proleptic_gregorian') and
                (year % 100 == 0) and
                (year % 400 != 0)):
                leap = False
            elif ((calendar in ['standard', 'gregorian']) and
                     (year % 100 == 0) and (year % 400 != 0) and
                     (year < 1583)):
                leap = False
        return leap

    def get_dpm(time, calendar='standard'):
        """
        return a array of days per month corresponding to the months provided in `months`
        """
        month_length = np.zeros(len(time), dtype=np.int)

        cal_days = dpm[calendar]

        for i, (month, year) in enumerate(zip(time.month, time.year)):
            month_length[i] = cal_days[month]
            if leap_year(year, calendar=calendar) and month == 2:
                month_length[i] += 1
        return month_length

## Open the `Dataset`

    [4]:

    ds = xr.tutorial.open_dataset('rasm').load()
    print(ds)

    ---------------------------------------------------------------------------
    FileNotFoundError                         Traceback (most recent call last)
    <ipython-input-4-3a5e1c385b97> in <module>
    ----> 1 ds = xr.tutorial.open_dataset('rasm').load()
          2 print(ds)

    /build/python-xarray-kTVyut/python-xarray-0.15.0/xarray/tutorial.py in open_dataset(name, cache, cache_dir, github_url, branch, **kws)
         76         # May want to add an option to remove it.
         77         if not _os.path.isdir(longdir):
    ---> 78             _os.mkdir(longdir)
         79
         80         url = "/".join((github_url, "raw", branch, fullname))

    FileNotFoundError: [Errno 2] No such file or directory: '/sbuild-nonexistent/.xarray_tutorial_data'

## Now for the heavy lifting:

We first have to come up with the weights, - calculate the month lengths for each monthly data record - calculate weights using `groupby('time.season')`

Finally, we just need to multiply our weights by the `Dataset` and sum allong the time dimension.

    [5]:

    # Make a DataArray with the number of days in each month, size = len(time)
    month_length = xr.DataArray(get_dpm(ds.time.to_index(), calendar='noleap'),
                                coords=[ds.time], name='month_length')

    # Calculate the weights by grouping by 'time.season'.
    # Conversion to float type ('astype(float)') only necessary for Python 2.x
    weights = month_length.groupby('time.season') / month_length.astype(float).groupby('time.season').sum()

    # Test that the sum of the weights for each season is 1.0
    np.testing.assert_allclose(weights.groupby('time.season').sum().values, np.ones(4))

    # Calculate the weighted average
    ds_weighted = (ds * weights).groupby('time.season').sum(dim='time')

    ---------------------------------------------------------------------------
    NameError                                 Traceback (most recent call last)
    <ipython-input-5-7e7fceb25104> in <module>
          1 # Make a DataArray with the number of days in each month, size = len(time)
    ----> 2 month_length = xr.DataArray(get_dpm(ds.time.to_index(), calendar='noleap'),
          3                             coords=[ds.time], name='month_length')
          4
          5 # Calculate the weights by grouping by 'time.season'.

    NameError: name 'ds' is not defined

    [6]:

    print(ds_weighted)

    ---------------------------------------------------------------------------
    NameError                                 Traceback (most recent call last)
    <ipython-input-6-eb2b7f87ce0b> in <module>
    ----> 1 print(ds_weighted)

    NameError: name 'ds_weighted' is not defined

    [7]:

    # only used for comparisons
    ds_unweighted = ds.groupby('time.season').mean('time')
    ds_diff = ds_weighted - ds_unweighted

    ---------------------------------------------------------------------------
    NameError                                 Traceback (most recent call last)
    <ipython-input-7-993ccc0e89d8> in <module>
          1 # only used for comparisons
    ----> 2 ds_unweighted = ds.groupby('time.season').mean('time')
          3 ds_diff = ds_weighted - ds_unweighted

    NameError: name 'ds' is not defined

    [8]:

    # Quick plot to show the results
    notnull = pd.notnull(ds_unweighted['Tair'][0])

    fig, axes = plt.subplots(nrows=4, ncols=3, figsize=(14,12))
    for i, season in enumerate(('DJF', 'MAM', 'JJA', 'SON')):
        ds_weighted['Tair'].sel(season=season).where(notnull).plot.pcolormesh(
            ax=axes[i, 0], vmin=-30, vmax=30, cmap='Spectral_r',
            add_colorbar=True, extend='both')

        ds_unweighted['Tair'].sel(season=season).where(notnull).plot.pcolormesh(
            ax=axes[i, 1], vmin=-30, vmax=30, cmap='Spectral_r',
            add_colorbar=True, extend='both')

        ds_diff['Tair'].sel(season=season).where(notnull).plot.pcolormesh(
            ax=axes[i, 2], vmin=-0.1, vmax=.1, cmap='RdBu_r',
            add_colorbar=True, extend='both')

        axes[i, 0].set_ylabel(season)
        axes[i, 1].set_ylabel('')
        axes[i, 2].set_ylabel('')

    for ax in axes.flat:
        ax.axes.get_xaxis().set_ticklabels([])
        ax.axes.get_yaxis().set_ticklabels([])
        ax.axes.axis('tight')
        ax.set_xlabel('')

    axes[0, 0].set_title('Weighted by DPM')
    axes[0, 1].set_title('Equal Weighting')
    axes[0, 2].set_title('Difference')

    plt.tight_layout()

    fig.suptitle('Seasonal Surface Air Temperature', fontsize=16, y=1.02)

    ---------------------------------------------------------------------------
    NameError                                 Traceback (most recent call last)
    <ipython-input-8-ded3065dafd3> in <module>
          1 # Quick plot to show the results
    ----> 2 notnull = pd.notnull(ds_unweighted['Tair'][0])
          3
          4 fig, axes = plt.subplots(nrows=4, ncols=3, figsize=(14,12))
          5 for i, season in enumerate(('DJF', 'MAM', 'JJA', 'SON')):

    NameError: name 'ds_unweighted' is not defined

    [9]:

    # Wrap it into a simple function
    def season_mean(ds, calendar='standard'):
        # Make a DataArray of season/year groups
        year_season = xr.DataArray(ds.time.to_index().to_period(freq='Q-NOV').to_timestamp(how='E'),
                                   coords=[ds.time], name='year_season')

        # Make a DataArray with the number of days in each month, size = len(time)
        month_length = xr.DataArray(get_dpm(ds.time.to_index(), calendar=calendar),
                                    coords=[ds.time], name='month_length')
        # Calculate the weights by grouping by 'time.season'
        weights = month_length.groupby('time.season') / month_length.groupby('time.season').sum()

        # Test that the sum of the weights for each season is 1.0
        np.testing.assert_allclose(weights.groupby('time.season').sum().values, np.ones(4))

        # Calculate the weighted average
        return (ds * weights).groupby('time.season').sum(dim='time')

    [ ]: